diff --git a/include/all.hxx b/include/all.hxx index d998540..67b7ebc 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); @@ -315,7 +315,7 @@ 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, + double prior_mass, double *AET_demand_cm, double *delta_theta, double *layer_thickness_cm, int *soil_type, struct soil_properties_ *soil_properties); /********************************************************************/ diff --git a/src/aet.cxx b/src/aet.cxx index 4c81bd1..60cb8ab 100644 --- a/src/aet.cxx +++ b/src/aet.cxx @@ -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; @@ -46,10 +44,13 @@ 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); 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..bfcf3b9 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,10 +1107,12 @@ 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, + + 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); + actual_ET_demand = *AET_demand_cm; - current->theta = fmin(theta_new, theta_e); + current->theta = fmax(theta_r, 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); @@ -1152,10 +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; //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 - current->theta = fmin(theta_e, prior_mass/current->depth_cm + next->theta); + else { + if ((prior_mass/current->depth_cm + next->theta)depth_cm*(current->theta - (prior_mass/current->depth_cm + next->theta)); + 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 {//This is the case where theta>theta_r, which will be almost all of the time + current->theta = fmax(theta_r, fmin(theta_e, prior_mass/current->depth_cm + next->theta)); + } + } } else { @@ -1229,27 +1244,37 @@ 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, + // 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); + actual_ET_demand = *AET_demand_cm; - current->theta = fmin(theta_new, theta_e); - } + current->theta = fmax(theta_r, 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 + + 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 +1291,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 +1309,14 @@ 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; + while (fabs(mass_balance_error - tolerance) > 1.E-10) { + iter++; + if (iter>1e4) { + *AET_demand_cm = *AET_demand_cm + mass_balance_error; + actual_ET_demand = *AET_demand_cm; + break; + } if (current_mass < mass_timestep) { depth_new += 0.01 * factor; @@ -1308,6 +1347,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); @@ -1336,7 +1376,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); @@ -1361,7 +1402,7 @@ extern void lgar_move_wetting_fronts(double timestep_h, double *volin_cm, int wf frozen_factor, head, soil_properties); *volin_cm = bottom_boundary_flux_cm; - + // 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 +1423,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 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. + 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"); @@ -1438,6 +1487,8 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str current = *head; if (verbosity.compare("high") == 0) { + printf("State before merging wetting fronts...\n"); + listPrint(*head); printf("Merging wetting fronts... \n"); } @@ -1527,6 +1578,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,27 +1608,31 @@ 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); double mbal_Z_correction = mbal_correction / (theta_new - next_to_next->theta); // this is the new wetting front depth + if (isinf(mbal_Z_correction)){//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. + mbal_Z_correction = cum_layer_thickness_cm[layer_num+1]/(1e3); + } double depth_new = cum_layer_thickness_cm[layer_num] + mbal_Z_correction; // this is the new wetting front absolute depth - current->depth_cm = cum_layer_thickness_cm[layer_num]; next->theta = theta_new; @@ -1591,7 +1647,7 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, } 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 +1722,21 @@ 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); + printf("Bottom boundary flux = %lf \n",bottom_flux_cm); + } + return bottom_flux_cm; - + } @@ -1730,7 +1785,7 @@ 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; //In LGARTO, should be the highest to_bottom WF above the one that got deleted before the next to_bottom==FALSE one, but *head is fine for LGAR // update psi and theta for all wetting fronts above the current wetting front while (current_local->layer_num < layer_num_k) { @@ -1740,9 +1795,10 @@ 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; + + 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 +1813,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 +1865,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 +1902,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 +1922,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 +1955,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; should really print like once per model run or so, not every time step double ponded_depth_temp = *ponded_depth_cm; @@ -1947,7 +2022,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 +2065,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 +2083,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 +2326,12 @@ 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 (layer_num, wf_num, K): %d %d %lf \n", layer_num, current->front_num, K_cm_per_h); + 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(); } @@ -2337,15 +2427,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, + double prior_mass, double *AET_demand_cm, double *delta_theta, double *delta_thickness, int *soil_type, struct soil_properties_ *soil_properties) { 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 +2455,17 @@ 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; + while (delta_mass > tolerance) { - + iter++; + + 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 +2495,7 @@ 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); + mass_layers += delta_thickness[layer_num] * (theta - delta_theta[layer_num]); for (int k=1; k tolerance){ + *AET_demand_cm = *AET_demand_cm - fabs(delta_mass - tolerance); } return theta; diff --git a/src/soil_funcs.cxx b/src/soil_funcs.cxx index 628f54b..bfedb11 100755 --- a/src/soil_funcs.cxx +++ b/src/soil_funcs.cxx @@ -26,15 +26,18 @@ /* 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 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 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){ // 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 +47,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 +65,52 @@ 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 = (h_i-h_f)/(double)nint; + dh = dh*0.01; //factor used to make dh small to begin with; dh begins small and is adaptively changed + 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, K1 and K2 differ by more than 2 percent. The factor of 2 percent seemed to offer the optimal intersection of accuracy and speed. + if ( (K1-K2)/K2 > 0.02 ){//if K1 disagrees with K2 by more than this fraction, then dh is made smaller + dh = dh*0.5; + } + else {//but if K1 and K2 are within a certain fraction of each other, then dh is made larger + dh = dh*10.0; + } + + if (h2 depth_wf_b = {1.873813, 44.00,175.0, 200.0}; // in cm - std::vector depth_wf_b = {4.3762450616157, 44.00,175.0, 200.0}; // in cm - std::vector theta_wf_b = {0.2153965837523, 0.172703948143618, 0.252113867764474, 0.179593529195751}; + std::vector depth_wf_b = {4.55355, 44.00,175.0, 200.0}; // in cm + std::vector theta_wf_b = {0.213716, 0.172703948143618, 0.252113867764474, 0.179593529195751}; int m_to_cm = 100; int m_to_mm = 1000; @@ -518,7 +518,7 @@ int main(int argc, char *argv[]) // check total infiltration, AET, and PET. double infiltration_check_mm = 1.896; // in mm - double AET_check_mm = 0.0288829525; // in mm + double AET_check_mm = 0.029801; // in mm double PET_check_mm = 0.104; // in mm double infiltration_computed = 0.0; double PET_computed = 0.0; @@ -562,6 +562,9 @@ int main(int argc, char *argv[]) std::stringstream errMsg; errMsg << "Error between benchmark and simulated AET is "<< fabs(AET_check_mm - AET_computed * m_to_mm) << " which is unexpected. \n"; + printf("benchmark AET: %lf \n", AET_check_mm); + printf("computed AET: %lf \n", AET_computed*m_to_mm); + throw std::runtime_error(errMsg.str()); }