Skip to content

Commit

Permalink
tested LASAM over approximately 40k parameter sets. Resolved mass bal…
Browse files Browse the repository at this point in the history 
…ance errors and infinite loop errors that became apparent when LASAM was run with many parameter sets. Also updated field capacity formulation to use a head rather than a relative water content, and updated caluclation for capillary drive G to use an adaptive integral, both as in LGARTO. Cleaned up some commenting and redundancy after edits were made.
  • Loading branch information
Peter La Follette authored and Peter La Follette committed Apr 2, 2024
1 parent 82d63c3 commit 2b0bb13
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 114 deletions.
3 changes: 0 additions & 3 deletions src/aet.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,6 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin
struct wetting_front *current;

double theta_wp;

// double relative_moisture_at_which_PET_equals_AET = 0.75;

double Se,theta_e,theta_r;
double vg_a, vg_m, vg_n;
Expand All @@ -46,7 +44,6 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin
vg_n = soil_properties[soil_num].vg_n;

// compute theta field capacity
// double theta_fc = (theta_e - theta_r) * relative_moisture_at_which_PET_equals_AET + theta_r;
double head_at_which_PET_equals_AET_cm = 340.9;//*10/33; //340.9 is 0.33 atm, expressed in water depth, which is a good field capacity for most soils.
//Coarser soils like sand will have a field capacity of 0.1 atm or so.
double theta_fc = calc_theta_from_h(head_at_which_PET_equals_AET_cm, vg_a,vg_m, vg_n, theta_e, theta_r);
Expand Down
30 changes: 9 additions & 21 deletions src/lgar.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1154,7 +1154,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
the layer depth */
current->depth_cm = fmin(current->depth_cm, column_depth);

double theta_old = current->theta;
double theta_old = current->theta; //might not be necessary
if (current->dzdt_cm_per_h == 0.0 && current->to_bottom == FALSE) // a new front was just created, so don't update it.
current->theta = current->theta;
else
Expand Down Expand Up @@ -1337,11 +1337,10 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf

// loop to adjust the depth for mass balance
int iter = 0;
bool is_WF_deeper_than_domain = FALSE;
while (fabs(mass_balance_error - tolerance) > 1.E-10) {
iter++;
if (iter>10000000) {
printf("mass balance closure not possible within 10000000 iterations. case 1. Timeout \n");
printf("mass balance closure not possible within 10000000 iterations. Timeout \n");
// printf("depth_new: %lf \n", depth_new);
printf("states: \n");
listPrint(*head);
Expand All @@ -1366,12 +1365,6 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
current_mass = lgar_calc_mass_bal(cum_layer_thickness_cm, *head);
mass_balance_error = fabs(current_mass - mass_timestep);

if (depth_new > cum_layer_thickness_cm[num_layers]){
is_WF_deeper_than_domain = TRUE;
// *volin_cm = *volin_cm - mass_balance_error;
// break;
}

}

}
Expand Down Expand Up @@ -1401,9 +1394,10 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
all wetting fronts, and then a new loop does lower boundary crossing for all wetting fronts. this is a thorough
way to deal with these scenarios. */

//As of 20 March 2024, we approach wetting front merging / layer boundary crossing / lower boundary crossing in a new way.
//As of 20 March 2024, we have added a new method by which wetting fronts merge or cross a layer or the model lower boundary, although it is commented out for now.
//we now have a while loop that detects if these scenarios are necessary and then only does them if so, instead of relying on
//a fixed order of merging and crossing (which was previously merge->cross layers->merge->cross lower bdy).
//a fixed order of merging and crossing (which was previously merge->cross layers->merge->cross lower bdy), it can automatically be detected which mergeing / crossing events are necessary.
//however, thius code slowed down some simulations by orders of magnitude. Since the old version of relying on a fixed, predetermined order of merging and crossing seems to work, I'd say don't fix it.

/* Note: we check for dry over wet case before we call merge_wetting_fronts to avoid negative wetting fronts
due to unknown corner/rare cases */
Expand Down Expand Up @@ -1509,7 +1503,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf
// make sure all psi values are updated
// there is a rare error where, after a wetting front crosses a layer boundary to a soil layer with a much more sensitive soil water retention curve, and the wetting front is near but not at saturation,
// the barely unsatruated wetting front in the new layer will update its psi value as 0. This causes unequal psi values among adjacent wetting fronts in different layers and then a mass balance error.
// For example, consider the following: (1.0/(pow(1.0+pow(0.001039*0.00296620706633,2.938410),0.659680))*(0.618625-0.052623)+0.052623). This expression, using specific values in calc_theta_from_h, will yield a theta of theta_e (0.618625), while having a psi value that is slightly greater than 0, which is 0.00296620706633.
// For example, consider the following: (1.0/(pow(1.0+pow(0.001039*0.00296620706633,2.938410),0.659680))*(0.618625-0.052623)+0.052623). This expression, using specific values in calc_theta_from_h, will yield a theta of theta_e (0.618625), while having a psi value that is slightly greater than 0, which is 0.00296620706633 in this case.
// While that is ok for this soil layer in particular, adjacent wetting fronts above this one with a less sensitive soil water retention curve will yield a non-theta_e value for the psi value that is slightly above 0.
// A solution is to either not run the following code, or to not run it when the wetting front is very close to saturation with a very sensitive soil water retention curve. Adding code that runs the following only for psi>1.

Expand Down Expand Up @@ -2076,7 +2070,7 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste

// if ( (current_mass)/max_storage > 0.99 ){
// printf("warning: vadose zone is 99 percent full or greater. If you are using the model in an environment with more precipitation than PET, LGAR is not an appropriate model because its lower boundary condition is no flow. \n ");
// } //turning off for now
// } //turning off for now; should really print like once per model run or so, not every time step

double ponded_depth_temp = *ponded_depth_cm;

Expand Down Expand Up @@ -2115,12 +2109,6 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste
// order is important here; assign zero to ponded depth once we compute volume in and runoff
*volin_this_timestep = fmin(*ponded_depth_cm, fp_cm); //
runoff = *ponded_depth_cm < fp_cm ? 0.0 : (*ponded_depth_cm - *volin_this_timestep);
// if (head->psi_cm<1e-8 && head->depth_cm==(cum_layer_thickness_cm[-1])){// ok this works for 1 layer, but generally for multilayer you want to check if the storage is saturated
//replaced this code with code that reduces f_p after if's calculated if there is not enough room in the soil (currently line 1915)
// runoff = *ponded_depth_cm;
// *volin_this_timestep = 0;
// listPrint(head);
// }
*ponded_depth_cm = 0.0;

}
Expand Down Expand Up @@ -2577,7 +2565,7 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm
iter++;

// if (iter>10000000) {
// printf("mass balance closure not possible within 10000000 iterations. case 2. Timeout \n");
// printf("mass balance closure not possible within 10000000 iterations. Timeout \n");
// printf("psi_cm_loc: %lf \n", psi_cm_loc);
// listPrint(*head);
// abort();
Expand Down Expand Up @@ -2693,7 +2681,7 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm

}


//while this function may or may not be used in LASAM, it will be used in LGARTO
extern int lgar_correction_type(int num_layers, double* cum_layer_thickness_cm, struct wetting_front* head){
int correction_type = 0;
struct wetting_front *current = head;
Expand Down
91 changes: 1 addition & 90 deletions src/soil_funcs.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,100 +26,11 @@
/* used for redistribution following the equation published by Ogden and Saghafian 1995. */
/* to compile: "gcc one_block.c -lm" */
/* */
/* author: Fred Ogden, June, 2021, edited by Peter La Follette in 2023 */
/* author: Fred Ogden, June, 2021, edited by Peter La Follette in 2023 for adaptive integral */
//updated to calculate G using an adaptive integral method, rather than always using a fixed number of trapezoids.
//This really helps when the limits of integration in terms of psi are quite far apart, and also helps to save runtime.
/***********************************************************************************************/

// extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, double theta_e, double theta_r,
// double vg_alpha, double vg_n, double vg_m, double h_min, double Ksat, int nint,
// double lambda, double bc_psib_cm)
// {
// double Geff; // this is the result to be returned.

// if (use_closed_form_G==false) {
// // local variables
// // note: units of h in cm. units of K in cm/s
// double h_i,h_f,Se_i,Se_f; // variables to store initial and final values
// double Se;
// double h2; // the head at the right-hand side of the trapezoid being integrated [m]
// double dh; // the delta h over which integration is performed [m]
// double Se1,Se2; // the scaled moisture content on left- and right-hand side of trapezoid
// double K1,K2; // the K(h) values on the left and right of the region dh integrated [m]


// Se_i = calc_Se_from_theta(theta1,theta_e,theta_r); // scaled initial water content (0-1) [-]
// Se_f = calc_Se_from_theta(theta2,theta_e,theta_r); // scaled final water content (0-1) [-]

// h_i = calc_h_from_Se(Se_i,vg_alpha,vg_m,vg_n); // capillary head associated with Se_i [cm]
// h_f = calc_h_from_Se(Se_f,vg_alpha,vg_m,vg_n); // capillary head associated with Se_f [cm]

// if(h_i < h_min) {/* if the lower limit of integration is less than h_min FIXME */
// //return h_min; // commenting out as this is not used in the Python version
// }

// if (verbosity.compare("high") == 0) {
// // debug statements to see if calc_Se_from_h function is working properly
// Se = calc_Se_from_h(h_i,vg_alpha,vg_m,vg_n);
// printf("Se_i = %8.6lf, Se_inverse = %8.6lf\n", Se_i, Se);

// Se = calc_Se_from_h(h_f,vg_alpha,vg_m,vg_n);
// printf("Se_f = %8.6lf, Se_inverse = %8.6lf\n", Se_f, Se);
// }

// // nint = number of "dh" intervals to integrate over using trapezoidal rule
// dh = (h_f-h_i)/(double)nint;
// Geff = 0.0;

// // integrate k(h) dh from h_i to h_f, using trapezoidal rule, with subscript
// // 1 denoting the left-hand side of the trapezoid, and 2 denoting the right-hand side

// Se1 = Se_i; // could just use Se_i in next statement. Done 4 completeness.
// K1 = calc_K_from_Se(Se1, Ksat, vg_m);
// h2 = h_i + dh;

// for(int i=1; i<=nint; i++) {

// Se2 = calc_Se_from_h(h2, vg_alpha, vg_m, vg_n);
// K2 = calc_K_from_Se(Se2, Ksat, vg_m);
// Geff += (K1+K2)*dh/2.0; // trapezoidal rule
// // reset for next time through loop
// K1 = K2;
// h2 += dh;
// }

// Geff = fabs(Geff/Ksat); // by convention Geff is a positive quantity

// if (verbosity.compare("high") == 0){
// printf ("Capillary suction (G) = %8.6lf \n", Geff);
// }

// return Geff;

// }
// else {

// double Se_f = calc_Se_from_theta(theta1,theta_e,theta_r); // the scaled moisture content of the wetting front
// double Se_i = calc_Se_from_theta(theta2,theta_e,theta_r); // the scaled moisture content below the wetting front
// double H_c = bc_psib_cm*((2+3*lambda) / (1+3*lambda)); // Green ampt capillary drive parameter, which can be used in the approximation of G with the Brooks-Corey model (See Ogden and Saghafian, 1997)

// Geff = H_c*(pow(Se_i,(3+1/lambda))-pow(Se_f,(3+1/lambda)))/(1-pow(Se_f,(3+1/lambda)));

// if (isinf(Geff)) {
// Geff = H_c;
// }
// if (isnan(Geff)) {
// Geff = H_c;
// }

// return Geff;

// }

// }



extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, double theta_e, double theta_r,
double vg_alpha, double vg_n, double vg_m, double h_min, double Ksat, int nint, double lambda, double bc_psib_cm)

Expand Down

0 comments on commit 2b0bb13

Please sign in to comment.