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.

Rebasing to incorporate recent changes.
  • Loading branch information
Peter La Follette authored and Peter La Follette committed Apr 2, 2024
1 parent 7d0d559 commit 82d63c3
Show file tree
Hide file tree
Showing 5 changed files with 540 additions and 97 deletions.
11 changes: 7 additions & 4 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 All @@ -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);
Expand Down Expand Up @@ -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
Expand Down
8 changes: 6 additions & 2 deletions src/aet.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);
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 82d63c3

Please sign in to comment.