Skip to content

Commit

Permalink
Ptl mbal testing (#21)
Browse files Browse the repository at this point in the history 
* 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.

* 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. Cleaned up some commenting and redundancy after edits were made.

* updated unit test

* updating with suggestions from PR review

* Updated lgar.cxx in response to a few comments

* Removed a few comments that were not necessary

* when theta updates, it is now enforced to be greater than theta_r and less than theta_e

* increased clarity of comment detailing how lgar_theta_mass_balance will never yield theta<theta_r

---------

Co-authored-by: Peter La Follette <[email protected]>
  • Loading branch information
@peterlafollette
peterlafollette and Peter La Follette authored Apr 4, 2024
1 parent 7d0d559 commit 4c21c46
Show file tree
Hide file tree
Showing 6 changed files with 260 additions and 104 deletions.
6 changes: 3 additions & 3 deletions include/all.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);

/********************************************************************/
Expand Down
7 changes: 4 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,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);
Expand Down
24 changes: 20 additions & 4 deletions src/bmi_lgar.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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,
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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\
Expand Down
Loading

0 comments on commit 4c21c46

Please sign in to comment.