From 82d63c3fecf269c31db954153e91bce2f1a2fbd5 Mon Sep 17 00:00:00 2001 From: Peter La Follette Date: Mon, 25 Mar 2024 10:41:21 -0700 Subject: [PATCH] tested LASAM over approximately 40k parameter sets. Resolved mass balance 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. Rebasing to incorporate recent changes. --- include/all.hxx | 11 +- src/aet.cxx | 8 +- src/bmi_lgar.cxx | 24 ++- src/lgar.cxx | 415 ++++++++++++++++++++++++++++++++++++++------- src/soil_funcs.cxx | 179 +++++++++++++++---- 5 files changed, 540 insertions(+), 97 deletions(-) diff --git a/include/all.hxx b/include/all.hxx index d998540..ac71b82 100755 --- a/include/all.hxx +++ b/include/all.hxx @@ -272,12 +272,12 @@ extern int lgar_read_vG_param_file(char const* vG_param_file_name, int num_soil_ struct soil_properties_ *soil_properties); // creates a surficial front (new top most wetting front) -extern void lgar_create_surficial_front(double *ponded_depth_cm, double *volin, double dry_depth, +extern void lgar_create_surficial_front(int num_layers, double *ponded_depth_cm, double *volin, double dry_depth, double theta1, int *soil_type, double *cum_layer_thickness_cm, double *frozen_factor, struct wetting_front **head, struct soil_properties_ *soil_properties); // computes the infiltration capacity, fp, of the soil -extern double lgar_insert_water(bool use_closed_form_G, int nint, double timestep_h, double *ponded_depth, +extern double lgar_insert_water(bool use_closed_form_G, int nint, double timestep_h, double AET_demand_cm, double *ponded_depth, double *volin_this_timestep, double precip_timestep_cm, int wf_free_drainge_demand, int num_layers, double ponded_depth_max_cm, int *soil_type, double *cum_layer_thickness_cm, double *frozen_factor, struct wetting_front* head, struct soil_properties_ *soil_properties); @@ -288,6 +288,9 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *ponded_depth_cm, double *cum_layer_thickness_cm, int *soil_type_by_layer, double *frozen_factor, struct wetting_front** head, struct wetting_front* state_previous, struct soil_properties_ *soil_properties); +//this allows for detection of the type of correction necessary (merge, layer cross, or lower bdy cross) +extern int lgar_correction_type(int num_layers, double* cum_layer_thickness_cm, struct wetting_front* head); + // the subroutine merges the wetting fronts; called from lgar_move_wetting_fronts extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, struct wetting_front** head, struct soil_properties_ *soil_properties); @@ -315,8 +318,8 @@ extern int wetting_front_free_drainage(struct wetting_front* head); // computes updated theta (soil moisture content) after moving down a wetting front; called for each wetting front to ensure mass is conserved extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm, double new_mass, - double prior_mass, double *delta_theta, double *layer_thickness_cm, - int *soil_type, struct soil_properties_ *soil_properties); + double prior_mass, double *AET_demand_cm, double *delta_theta, double *layer_thickness_cm, + int *soil_type, struct soil_properties_ *soil_properties, struct wetting_front** head); /********************************************************************/ // Bmi functions diff --git a/src/aet.cxx b/src/aet.cxx index 4c81bd1..1df0688 100644 --- a/src/aet.cxx +++ b/src/aet.cxx @@ -29,7 +29,7 @@ extern double calc_aet(double PET_timestep_cm, double time_step_h, double wiltin double theta_wp; - double relative_moisture_at_which_PET_equals_AET = 0.75; + // double relative_moisture_at_which_PET_equals_AET = 0.75; double Se,theta_e,theta_r; double vg_a, vg_m, vg_n; @@ -46,10 +46,14 @@ 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 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); double wp_head_theta = calc_theta_from_h(wilting_point_psi_cm, vg_a,vg_m, vg_n, theta_e, theta_r); + theta_wp = (theta_fc - wp_head_theta)*1/2 + wp_head_theta; // theta_50 in python Se = calc_Se_from_theta(theta_wp,theta_e,theta_r); diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index 5a215f7..7d8002b 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -238,20 +238,37 @@ Update() double temp_pd = 0.0; // necessary to assign zero precip due to the creation of new wetting front; AET will still be taken out of the layers // move the wetting fronts without adding any water; this is done to close the mass balance + // and also to merge / cross if necessary lgar_move_wetting_fronts(subtimestep_h, &temp_pd, wf_free_drainage_demand, volend_subtimestep_cm, num_layers, &AET_subtimestep_cm, state->lgar_bmi_params.cum_layer_thickness_cm, state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.frozen_factor, &state->head, state->state_previous, state->soil_properties); + + if (temp_pd != 0.0){ //if temp_pd != 0.0, that means that some water left the model through the lower model bdy + volrech_subtimestep_cm = temp_pd; + volrech_timestep_cm += volrech_subtimestep_cm; + temp_pd = 0.0; + } // depth of the surficial front to be created dry_depth = lgar_calc_dry_depth(use_closed_form_G, nint, subtimestep_h, &delta_theta, state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.cum_layer_thickness_cm, state->lgar_bmi_params.frozen_factor, state->head, state->soil_properties); + + if (verbosity.compare("high") == 0) { + printf("State before moving creating new WF...\n"); + listPrint(state->head); + } - lgar_create_surficial_front(&ponded_depth_subtimestep_cm, &volin_subtimestep_cm, dry_depth, state->head->theta, + lgar_create_surficial_front(num_layers, &ponded_depth_subtimestep_cm, &volin_subtimestep_cm, dry_depth, state->head->theta, state->lgar_bmi_params.layer_soil_type, state->lgar_bmi_params.cum_layer_thickness_cm, state->lgar_bmi_params.frozen_factor, &state->head, state->soil_properties); + if (verbosity.compare("high") == 0) { + printf("State after moving creating new WF...\n"); + listPrint(state->head); + } + state->state_previous = NULL; state->state_previous = listCopy(state->head); @@ -269,7 +286,7 @@ Update() if (ponded_depth_subtimestep_cm > 0 && !create_surficial_front) { - volrunoff_subtimestep_cm = lgar_insert_water(use_closed_form_G, nint, subtimestep_h, &ponded_depth_subtimestep_cm, + volrunoff_subtimestep_cm = lgar_insert_water(use_closed_form_G, nint, subtimestep_h, AET_subtimestep_cm, &ponded_depth_subtimestep_cm, &volin_subtimestep_cm, precip_subtimestep_cm_per_h, wf_free_drainage_demand, num_layers, ponded_depth_max_cm, state->lgar_bmi_params.layer_soil_type, @@ -319,7 +336,6 @@ Update() volin_subtimestep_cm = volin_subtimestep_cm_temp; } - /*----------------------------------------------------------------------*/ // calculate derivative (dz/dt) for all wetting fronts lgar_dzdt_calc(use_closed_form_G, nint, ponded_depth_subtimestep_cm, state->lgar_bmi_params.layer_soil_type, @@ -360,7 +376,7 @@ Update() listPrint(state->head); } - bool unexpected_local_error = fabs(local_mb) > 1.0E-6 ? true : false; + bool unexpected_local_error = fabs(local_mb) > 1.0E-4 ? true : false; if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0 || unexpected_local_error) { printf("\nLocal mass balance at this timestep... \n\ diff --git a/src/lgar.cxx b/src/lgar.cxx index 9099db6..ceea916 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -787,7 +787,7 @@ extern void frozen_factor_hydraulic_conductivity(struct lgar_bmi_parameters lgar } - assert (layer_temp > 100.0); // just a check to ensure the while loop executes at least once + // assert (layer_temp > 100.0); // just a check to ensure the while loop executes at least once assert (count > 0); // just a check to ensure the while loop executes at least once layer_temp /= count; // layer-averaged temperature @@ -1107,8 +1107,10 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf prior_mass += precip_mass_to_add - (free_drainage_demand + actual_ET_demand); // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front - double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, - delta_thetas, delta_thickness, soil_type, soil_properties); + + double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm, + delta_thetas, delta_thickness, soil_type, soil_properties, head); + actual_ET_demand = *AET_demand_cm; current->theta = fmin(theta_new, theta_e); @@ -1152,11 +1154,23 @@ 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; 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 current->theta = fmin(theta_e, prior_mass/current->depth_cm + next->theta); + if (current->theta < theta_r){ + //the idea here is that in some cases, the reduction in theta via WF movement or AET will be intense enough such that theta goes below theta_r. + //it requires a fairly unusual soil, which I encountered during random parameter sampling. + + double mass_before_theta_went_below_theta_r = lgar_calc_mass_bal(cum_layer_thickness_cm, *head); + listDeleteFront(current->front_num, head); + double mass_after_theta_went_below_theta_r = lgar_calc_mass_bal(cum_layer_thickness_cm, *head); + *AET_demand_cm = *AET_demand_cm - fabs(mass_before_theta_went_below_theta_r - mass_after_theta_went_below_theta_r); + actual_ET_demand = *AET_demand_cm; + + } } else { @@ -1229,27 +1243,65 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf delta_thickness[layer_num] = current->depth_cm - cum_layer_thickness_cm[layer_num-1]; double free_drainage_demand = 0; + if (wf_free_drainage_demand == wf) prior_mass += precip_mass_to_add - (free_drainage_demand + actual_ET_demand); - // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front - double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, - delta_thetas, delta_thickness, soil_type, soil_properties); + // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front + + double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, AET_demand_cm, + delta_thetas, delta_thickness, soil_type, soil_properties, head); + actual_ET_demand = *AET_demand_cm; current->theta = fmin(theta_new, theta_e); - } + } double Se = calc_Se_from_theta(current->theta,theta_e,theta_r); current->psi_cm = calc_h_from_Se(Se, vg_a, vg_m, vg_n); + } + // if f_p (predicted infiltration) causes theta > theta_e, mass correction is needed. // depth of the wetting front is increased to close the mass balance when theta > theta_e. // l == 1 is the last iteration (top most wetting front), so do a check on the mass balance) // this part should be moved out of here to a subroutine; add a call to that subroutine - if (wf == 1) { - + if (wf == 1) { + //first do a test layer bdy cross and then recalc wf_free_drainage_demand + lgar_wetting_fronts_cross_layer_boundary(num_layers, cum_layer_thickness_cm, soil_type, frozen_factor, + head, soil_properties); //replacing with more complete code that does merging and crossing as necessary + + // int correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // printf("correction_type: %d \n", correction_type); + + // while (correction_type!=0){ + // if (correction_type==1){ + // // *************************** MERGE ************************************ + // lgar_merge_wetting_fronts(soil_type, frozen_factor, head, soil_properties); + // correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // } + + // if (correction_type==2){ + // // ************************ CROSS LAYER ********************************* + // lgar_wetting_fronts_cross_layer_boundary(num_layers, cum_layer_thickness_cm, soil_type, frozen_factor, + // head, soil_properties); + // correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // } + + // if (correction_type==3){ + // // ************************ CROSS BOUNDARY ******************************** + // //lower bound + // bottom_boundary_flux_cm += lgar_wetting_front_cross_domain_boundary(cum_layer_thickness_cm[num_layers], soil_type, + // frozen_factor, head, soil_properties); + // *volin_cm = bottom_boundary_flux_cm; + // correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // } + // } + + + wf_free_drainage_demand = wetting_front_free_drainage(*head); + // if ((wf == wf_free_drainage_demand) && (current->theta>=theta_e) ) { int soil_num_k1 = soil_type[wf_free_drainage_demand]; double theta_e_k1 = soil_properties[soil_num_k1].theta_e; @@ -1266,8 +1318,15 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf double mass_balance_error = fabs(current_mass - mass_timestep); // mass error double factor = 1.0; + if (wf_free_drainage->next!=NULL){ + if (fabs(wf_free_drainage->theta - wf_free_drainage->next->theta)<0.01){// if two adjacent theta values are quite close, an initial factor of 1.0 might not make the mass balance close within 10000 iterations + factor = factor / fabs(wf_free_drainage->theta - wf_free_drainage->next->theta); + } + } + + // double factor = fmax(1,current->psi_cm/100); speed optimization should look at optimal factor values bool switched = false; - double tolerance = 1e-12; + double tolerance = 1e-10; // check if the difference is less than the tolerance if (mass_balance_error <= tolerance) { @@ -1277,7 +1336,17 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf double depth_new = wf_free_drainage->depth_cm; // loop to adjust the depth for mass balance - while (fabs(mass_balance_error - tolerance) > 1.E-12) { + 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("depth_new: %lf \n", depth_new); + printf("states: \n"); + listPrint(*head); + abort(); + } if (current_mass < mass_timestep) { depth_new += 0.01 * factor; @@ -1297,6 +1366,12 @@ 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; + } + } } @@ -1308,6 +1383,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf // end of the for loop /*******************************************************************/ + if (verbosity.compare("high") == 0) { printf("State after moving but before merging wetting fronts...\n"); listPrint(*head); @@ -1325,9 +1401,16 @@ 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. + //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). + /* 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 */ + + + double mass_change = 0.0; // *************************** MERGE ************************************ @@ -1336,7 +1419,8 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf if (is_dry_over_wet_wf) lgar_fix_dry_over_wet_wetting_fronts(&mass_change, cum_layer_thickness_cm, soil_type, head, soil_properties); - + + lgar_merge_wetting_fronts(soil_type, frozen_factor, head, soil_properties); @@ -1362,6 +1446,45 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf *volin_cm = bottom_boundary_flux_cm; + + // double mass_change = 0.0; + + // //20 March 2024: replacing merge/cross/merge/cross lower bdy fixed order code with function that applies only necessary corrections + // //this is code that could / should be implemented at some point, however it looks like running merge / cross / moerge / lower bdy cross in this fixed order does actually work across ~40k parameter sets: no need to change it for now + + // // check if dry over wet wetting fronts exist before calling merge + // bool is_dry_over_wet_wf = lgar_check_dry_over_wet_wetting_fronts(*head); + // if (is_dry_over_wet_wf) + // lgar_fix_dry_over_wet_wetting_fronts(&mass_change, cum_layer_thickness_cm, soil_type, head, soil_properties); + + + // int correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + + // while (correction_type!=0){ + // if (correction_type==1){ + // // *************************** MERGE ************************************ + // lgar_merge_wetting_fronts(soil_type, frozen_factor, head, soil_properties); + // correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // } + + // if (correction_type==2){ + // // ************************ CROSS LAYER ********************************* + // lgar_wetting_fronts_cross_layer_boundary(num_layers, cum_layer_thickness_cm, soil_type, frozen_factor, + // head, soil_properties); + // correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // } + + // if (correction_type==3){ + // // ************************ CROSS BOUNDARY ******************************** + // //lower bound + // bottom_boundary_flux_cm += lgar_wetting_front_cross_domain_boundary(cum_layer_thickness_cm[num_layers], soil_type, + // frozen_factor, head, soil_properties); + // *volin_cm = bottom_boundary_flux_cm; + // correction_type = lgar_correction_type(num_layers, cum_layer_thickness_cm, *head); + // } + // } + + // reset current to head to fix any mass balance issues and dry-over-wet wetting fronts conditions is_dry_over_wet_wf = lgar_check_dry_over_wet_wetting_fronts(*head); @@ -1382,30 +1505,38 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf *AET_demand_cm = *AET_demand_cm - mass_change; } - /***********************************************/ // 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. + // 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. + current = *head; for (int wf=1; wf != listLength(*head); wf++) { - int soil_num_k = soil_type[current->layer_num]; + if (current->psi_cm>1.0){ + int soil_num_k = soil_type[current->layer_num]; - double theta_e_k = soil_properties[soil_num_k].theta_e; - double theta_r_k = soil_properties[soil_num_k].theta_r; - double vg_a_k = soil_properties[soil_num_k].vg_alpha_per_cm; - double vg_m_k = soil_properties[soil_num_k].vg_m; - double vg_n_k = soil_properties[soil_num_k].vg_n; + double theta_e_k = soil_properties[soil_num_k].theta_e; + double theta_r_k = soil_properties[soil_num_k].theta_r; + double vg_a_k = soil_properties[soil_num_k].vg_alpha_per_cm; + double vg_m_k = soil_properties[soil_num_k].vg_m; + double vg_n_k = soil_properties[soil_num_k].vg_n; - double Ksat_cm_per_h_k = frozen_factor[current->layer_num] * soil_properties[soil_num_k].Ksat_cm_per_h; + double Ksat_cm_per_h_k = frozen_factor[current->layer_num] * soil_properties[soil_num_k].Ksat_cm_per_h; - double Se = calc_Se_from_theta(current->theta,theta_e_k,theta_r_k); - current->psi_cm = calc_h_from_Se(Se, vg_a_k, vg_m_k, vg_n_k); - current->K_cm_per_h = calc_K_from_Se(Se, Ksat_cm_per_h_k, vg_m_k); + double Se = calc_Se_from_theta(current->theta,theta_e_k,theta_r_k); + current->psi_cm = calc_h_from_Se(Se, vg_a_k, vg_m_k, vg_n_k); + current->K_cm_per_h = calc_K_from_Se(Se, Ksat_cm_per_h_k, vg_m_k); + } current = current->next; } + if (verbosity.compare("high") == 0) printf("Moving/merging wetting fronts done... \n"); @@ -1437,6 +1568,11 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str struct wetting_front *next_to_next; current = *head; + if (verbosity.compare("high") == 0) { + printf("State before merging wetting fronts...\n"); + listPrint(*head); + } + if (verbosity.compare("high") == 0) { printf("Merging wetting fronts... \n"); } @@ -1527,6 +1663,7 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, if (verbosity.compare("high") == 0) { printf("Layer boundary crossing... \n"); } + for (int wf=1; wf != listLength(*head); wf++) { if (verbosity.compare("high") == 0) { @@ -1556,19 +1693,21 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, double current_theta = fmin(theta_e, current->theta); double overshot_depth = current->depth_cm - next->depth_cm; - double next_theta_e = soil_properties[soil_num+1].theta_e; - double next_theta_r = soil_properties[soil_num+1].theta_r; - double next_vg_a = soil_properties[soil_num+1].vg_alpha_per_cm; - double next_vg_m = soil_properties[soil_num+1].vg_m; - double next_vg_n = soil_properties[soil_num+1].vg_n; - //double next_Ksat_cm_per_h = soil_properties[soil_num+1].Ksat_cm_per_h * frozen_factor[current->layer_num]; + int soil_num_next = soil_type[layer_num+1]; + + double next_theta_e = soil_properties[soil_num_next].theta_e; + double next_theta_r = soil_properties[soil_num_next].theta_r; + double next_vg_a = soil_properties[soil_num_next].vg_alpha_per_cm; + double next_vg_m = soil_properties[soil_num_next].vg_m; + double next_vg_n = soil_properties[soil_num_next].vg_n; + //double next_Ksat_cm_per_h = soil_properties[soil_num_next].Ksat_cm_per_h * frozen_factor[current->layer_num]; double Se = calc_Se_from_theta(current->theta,theta_e, theta_r); current->psi_cm = calc_h_from_Se(Se, vg_a, vg_m, vg_n); current->K_cm_per_h = calc_K_from_Se(Se, Ksat_cm_per_h, vg_m); - // current psi with van Genuchten properties of the next layer to get new theta + // current psi with van Gunechten properties of the next layer to get new theta double theta_new = calc_theta_from_h(current->psi_cm, next_vg_a, next_vg_m, next_vg_n, next_theta_e, next_theta_r); double mbal_correction = overshot_depth * (current_theta - next->theta); @@ -1587,11 +1726,15 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, current->dzdt_cm_per_h = 0; current->to_bottom = TRUE; next->to_bottom = FALSE; + + if (isinf(next->depth_cm)){//in some rare cases, due to (very) different shapes of soil water rentention curves in adjacent soil layers near saturation, the WF that crossed a layer bdy can yield a theta value identical to the one below it, resulting in division by 0 and an infinite WF depth. + next->depth_cm = current->depth_cm + cum_layer_thickness_cm[layer_num+1]/(1e3); + } } if (verbosity.compare("high") == 0) { - printf("State after wetting fronts cross layer boundary...\n"); + printf("States after wetting fronts cross layer boundary...\n"); listPrint(*head); } current = current->next; @@ -1666,22 +1809,23 @@ extern double lgar_wetting_front_cross_domain_boundary(double domain_depth_cm, i break; } - if (verbosity.compare("high") == 0) { - printf("State after lowest wetting front contributes to flux through the bottom boundary...\n"); - listPrint(*head); - } - current = current->next; - - if (verbosity.compare("high") == 0){ - printf("Bottom boundary flux = %lf \n",bottom_flux_cm); } + if (bottom_flux_cm_temp!=0){//PTL: perhaps due to potential problems with precision error, this should just be replaces with a flag that gets toggled when bottom_boundary_flux_cm is changed from 0 to nonzero break; } } + if (verbosity.compare("high") == 0) { + printf("State after lowest wetting front contributes to flux through the bottom boundary...\n"); + listPrint(*head); + } + + if (verbosity.compare("high") == 0){ + printf("Bottom boundary flux = %lf \n",bottom_flux_cm); } + return bottom_flux_cm; - + } @@ -1730,9 +1874,18 @@ extern void lgar_fix_dry_over_wet_wetting_fronts(double *mass_change, double* cu // now this is the wetting front that was below the dry wetting front current->psi_cm = calc_h_from_Se(Se_k, vg_a_k, vg_m_k, vg_n_k); - struct wetting_front *current_local = *head; + struct wetting_front *current_local = *head; //shouldn't be head, should be the highest to_bottom WF above the one that got deleted before the next to_bottom==FALSE one + current_local = listFindFront(current->front_num - 1, *head, NULL); + if (listFindFront(current_local->front_num - 1, *head, NULL)!=NULL){ + while (listFindFront(current_local->front_num - 1, *head, NULL)->to_bottom==TRUE){ + current_local = listFindFront(current_local->front_num - 1, *head, NULL); + if (listFindFront(current_local->front_num - 1, *head, NULL)==NULL){ + break; + } + } + } - // update psi and theta for all wetting fronts above the current wetting front + // update psi and theta for all wetting fronts above the current wetting front //is now for all WFs above the current WF until the next to_bottom == FALSE one while (current_local->layer_num < layer_num_k) { int soil_num_k1 = soil_type[current_local->layer_num]; theta_e_k = soil_properties[soil_num_k1].theta_e; @@ -1740,9 +1893,16 @@ extern void lgar_fix_dry_over_wet_wetting_fronts(double *mass_change, double* cu vg_a_k = soil_properties[soil_num_k1].vg_alpha_per_cm; vg_m_k = soil_properties[soil_num_k1].vg_m; vg_n_k = soil_properties[soil_num_k1].vg_n; - double Se_l = calc_Se_from_theta(current->theta,theta_e_k,theta_r_k); - current_local->psi_cm = calc_h_from_Se(Se_l, vg_a_k, vg_m_k, vg_n_k); - current_local->theta = calc_theta_from_h(current->psi_cm, vg_a_k, vg_m_k, vg_n_k,theta_e_k,theta_r_k); + + current_local->psi_cm = current->psi_cm; + + if (current_local->next!=NULL){ + if (isnan(current_local->psi_cm)){//in some highly unrealistic soils tested (for example, with a volumetric water content of 0.67 with a capillary head of 1e6 cm), this fxn would yield nan. This code corrects that. + //code has been further corrected shuch that this should never happen (fixed an error where Se_l, an old variable, was negative when it couldn't be). Nonetheless this should also redundantly correct the states if it does somehow happen. + current_local->psi_cm = current_local->next->psi_cm; + } + } + current_local->theta = calc_theta_from_h(current->psi_cm, vg_a_k, vg_m_k, vg_n_k,theta_e_k,theta_r_k); current_local = current_local->next; } } @@ -1757,7 +1917,7 @@ extern void lgar_fix_dry_over_wet_wetting_fronts(double *mass_change, double* cu back to AET after deleting the drier front */ } - + current = current->next; if (current == NULL) @@ -1809,7 +1969,7 @@ extern bool lgar_check_dry_over_wet_wetting_fronts(struct wetting_front* head) in the current timestep, that is precipitation in the current and previous timesteps was greater than zero */ // ############################################################################################ -extern double lgar_insert_water(bool use_closed_form_G, int nint, double timestep_h, double *ponded_depth_cm, +extern double lgar_insert_water(bool use_closed_form_G, int nint, double timestep_h, double AET_demand_cm, double *ponded_depth_cm, double *volin_this_timestep, double precip_timestep_cm, int wf_free_drainage_demand, int num_layers, double ponded_depth_max_cm, int *soil_type, double *cum_layer_thickness_cm, double *frozen_factor, @@ -1846,6 +2006,8 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste if (number_of_wetting_fronts == num_layers) { Geff = 0.0; // i.e., case of no capillary suction, dz/dt is also zero for all wetting fronts + soil_num = soil_type[layer_num_fp]; + Ksat_cm_per_h = soil_properties[soil_num].Ksat_cm_per_h * frozen_factor[current->layer_num]; //23 feb 2024 } else { @@ -1864,8 +2026,8 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste double bc_psib_cm = soil_properties[soil_num].bc_psib_cm; Ksat_cm_per_h = soil_properties[soil_num].Ksat_cm_per_h * frozen_factor[current->layer_num]; - //Se = calc_Se_from_theta(theta,theta_e,theta_r); - //psi_cm = calc_h_from_Se(Se, vg_a, vg_m, vg_n); + // Se = calc_Se_from_theta(theta,theta_e,theta_r); + // psi_cm = calc_h_from_Se(Se, vg_a, vg_m, vg_n); Geff = calc_Geff(use_closed_form_G, theta_below, theta_e, theta_e, theta_r, vg_a, vg_n, vg_m, h_min_cm, Ksat_cm_per_h, nint, lambda, bc_psib_cm); @@ -1897,7 +2059,24 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste // if free drainge has to be included, which currently we don't, then the following will be set to hydraulic conductivity // of the deeepest layer if ((layer_num_fp == num_layers) && (current_free_drainage->theta == theta_e1) && (num_layers == number_of_wetting_fronts)) - f_p = 0.0; + f_p = fmin(f_p, AET_demand_cm/timestep_h); //the idea here is that, if the soil is completely saturated, a little bit of water can still enter if AET is sufficiently large + + //this code checks if there is enough storage available for infiltrating water. That is, f_p can only be as big as there is room for water. + double current_mass = lgar_calc_mass_bal(cum_layer_thickness_cm, head); + double max_storage = 0.0; + for (int k = 1; k < num_layers+1; k++) { + int layer_num = k; + soil_num = soil_type[layer_num]; + max_storage += soil_properties[soil_num].theta_e * (cum_layer_thickness_cm[k]-cum_layer_thickness_cm[k-1]); + }// would be more efficient to make max_storage a vairable that is not computed every time lgar_insert_water is called + + if (f_p > (max_storage - current_mass)/timestep_h){ + f_p = (max_storage - current_mass)/timestep_h; + } + + // 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 double ponded_depth_temp = *ponded_depth_cm; @@ -1936,6 +2115,12 @@ 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; } @@ -1947,7 +2132,7 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste /* This subroutine is called iff there is no surfacial front, it creates a new front and inserts ponded depth, and will return some amount if can't fit all water */ // ###################################################################################### -extern void lgar_create_surficial_front(double *ponded_depth_cm, double *volin, double dry_depth, +extern void lgar_create_surficial_front(int num_layers, double *ponded_depth_cm, double *volin, double dry_depth, double theta1, int *soil_type, double *cum_layer_thickness_cm, double *frozen_factor, struct wetting_front** head, struct soil_properties_ *soil_properties) { @@ -1990,7 +2175,8 @@ extern void lgar_create_surficial_front(double *ponded_depth_cm, double *volin, if (dry_depth < cum_layer_thickness_cm[1]) listInsertFirst(dry_depth, theta_e, front_num, layer_num, to_bottom, head); else - listInsertFirst(dry_depth, theta_e, front_num, layer_num, 1, head); + // listInsertFirst(dry_depth, theta_e, front_num, layer_num, 1, head); + listInsertFirst(dry_depth, theta_e, front_num, layer_num, to_bottom, head); //the idea here is that a new WF should never have to_bottom as 1 -- if it needs to merge with the one below it, it will //hp_cm = *ponded_depth_cm; } @@ -2007,6 +2193,16 @@ extern void lgar_create_surficial_front(double *ponded_depth_cm, double *volin, current->dzdt_cm_per_h = 0.0; //for now assign 0 to dzdt as it will be computed/updated in lgar_dzdt_calc function + if (current->next!=NULL){// sometimes a new WF immediately has to merge with another WF or cross a layer bdy + if (current->depth_cm>current->next->depth_cm){ + //do a merge cross merge + // I'm thinking it might not be necessary -- this should happen later. Leaving in for now however, + lgar_merge_wetting_fronts(soil_type, frozen_factor, head, soil_properties); + lgar_wetting_fronts_cross_layer_boundary(num_layers, cum_layer_thickness_cm, soil_type, frozen_factor, + head, soil_properties); + } + } + return; } @@ -2240,8 +2436,13 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so layer_num = current->layer_num; // what layer the front is in K_cm_per_h = current->K_cm_per_h; // K(theta) - if (K_cm_per_h <=0) { - printf("K is zero: %d %lf \n", layer_num, K_cm_per_h); + if (K_cm_per_h < 0) { + printf("K is negative: %d %lf \n", layer_num, K_cm_per_h); + printf("current->front_num: %d \n", current->front_num); + listPrint(head); + printf("Is your n value very close to 1? Very small n values can cause K to become 0. \n"); + //The parameter n must physically attain a value greater than 1. However, when n is small, and apparently less than 1.02, sometimes n can make K evaluate to 0, for larger values of psi. + //So, checking for K_cm_per_h <= 0 has been replaced by checking if K_cm_per_h is negative. K_cm_per_h should never be negative (although perhaps machine precision could make this occur, although we haven't seen it yet), but mathematically can be 0 in some rare cases. abort(); } @@ -2323,6 +2524,9 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so dzdt = 0.0; } + // if (dzdt>1e4){//In some cases, dzdt can mathematically be very large, for example if theta is very close to theta_b, and so the theta-theta_b term in the denominator of the dzdt calculation makes dzdt large. Because this loses physical meaning, we restrict dzdt. + // dzdt = 1e4; + // } current->dzdt_cm_per_h = dzdt; current = current->next; // point to the next link @@ -2337,15 +2541,16 @@ extern void lgar_dzdt_calc(bool use_closed_form_G, int nint, double h_p, int *so is within a tolerance. */ // ############################################################################################ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm, double new_mass, - double prior_mass, double *delta_theta, double *delta_thickness, - int *soil_type, struct soil_properties_ *soil_properties) + double prior_mass, double *AET_demand_cm, double *delta_theta, double *delta_thickness, + int *soil_type, struct soil_properties_ *soil_properties, struct wetting_front** head) { double psi_cm_loc = psi_cm; // location psi double delta_mass = fabs(new_mass - prior_mass); // mass different between the new and prior + // double original_delta_mass = delta_mass; double tolerance = 1e-12; - double factor = 1.0; + double factor = fmax(1,psi_cm/100);//was 1.0 previously. This code is far faster and seems to avoid loops with >10000 iterations. Low n values can cause this bool switched = false; // flag that determines capillary head to be incremented or decremented double theta = 0; // this will be updated and returned @@ -2364,8 +2569,25 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm // the loop increments/decrements the capillary head until mass difference between // the new and prior is within the tolerance + int iter = 0; + bool iter_aug_flag = FALSE; + bool went_below_theta_r_flag = FALSE; + while (delta_mass > tolerance) { - + iter++; + + // if (iter>10000000) { + // printf("mass balance closure not possible within 10000000 iterations. case 2. Timeout \n"); + // printf("psi_cm_loc: %lf \n", psi_cm_loc); + // listPrint(*head); + // abort(); + // } + + if (iter>1000 && iter_aug_flag==FALSE){ + factor = factor*100; + iter_aug_flag = TRUE; + } + if (new_mass > prior_mass) { psi_cm_loc += 0.1 * factor; switched = false; @@ -2395,6 +2617,15 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm theta = calc_theta_from_h(psi_cm_loc, soil_properties[soil_num].vg_alpha_per_cm, soil_properties[soil_num].vg_m, soil_properties[soil_num].vg_n,soil_properties[soil_num].theta_e, soil_properties[soil_num].theta_r); + + if (theta tolerance)" will attempt to increase psi until mass balance closure occurs, but if the + //AET value is large enough to make theta less than theta_r (or if dzdt was large enough to make the WF move such that the updated theta would be less than theta_r), + //the loop would never be able to achieve mass balance closure, and it will break when the one of above conditions is met. + //So, the next few lines beginning with "if (delta_mass > tolerance)" addresses the rare case where water would not be able to be extracted without theta becoming less than theta_r anyway + } + + if (delta_mass > tolerance){ + *AET_demand_cm = *AET_demand_cm - fabs(delta_mass - tolerance); + } + + if (thetanext; + struct wetting_front *next_to_next = NULL; + if (next!=NULL){ + next_to_next = current->next->next; + } + // need to consider case where next==NULL + + for (int wf = 1; wf != (listLength(head)); wf++) { + struct wetting_front *previous = listFindFront(wf - 1, head, NULL); + + if (next!=NULL){ + + if ( (current->depth_cm > next->depth_cm) && (current->layer_num == next->layer_num) && (!next->to_bottom) ){ + correction_type = 1; //this is surface-surface WF merging + } + + if ( (current->depth_cm > cum_layer_thickness_cm[current->layer_num]) && (next->depth_cm == cum_layer_thickness_cm[current->layer_num]) && (current->layer_num!=num_layers) ){ + correction_type = 2; //this is surface WF layer bdy crossing + } + } + + if ( (next_to_next == NULL) && (current->depth_cm > cum_layer_thickness_cm[current->layer_num]) ){ + correction_type = 3; //this is a surface WF crossing the model lower bdy + } + + if (correction_type!=0){ + break; + } + + current = next; + next = current->next; + if (next!=NULL){ + next_to_next = current->next->next; + } + } + + return correction_type; +} + + #endif diff --git a/src/soil_funcs.cxx b/src/soil_funcs.cxx index 628f54b..b252282 100755 --- a/src/soil_funcs.cxx +++ b/src/soil_funcs.cxx @@ -26,15 +26,107 @@ /* used for redistribution following the equation published by Ogden and Saghafian 1995. */ /* to compile: "gcc one_block.c -lm" */ /* */ -/* author: Fred Ogden, June, 2021, */ +/* author: Fred Ogden, June, 2021, edited by Peter La Follette in 2023 */ +//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) + 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) { + 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 @@ -44,7 +136,6 @@ extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, do 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) [-] @@ -63,28 +154,54 @@ extern double calc_Geff(bool use_closed_form_G, double theta1, double theta2, do 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; + dh = -1*(h_f-h_i)/(double)nint; + dh = dh*0.01; + 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++) { + h2 = h_f + dh; + + while(h2 0.02 ){//if K1 disagrees with K2 by more than this fraction, then dh is made smaller + dh = dh*0.5; + } + else if ( (K1-K2)/K2 < 0.03 ){//but if K1 and K2 are within a certain fraction of each other, then dh is made larger + dh = dh*10.0; + } + else{ + dh = dh*2.0; //but if K1 and K2 are within a certain fraction of each other, then dh is made larger + } + + if (h2