updating LASAM to take field capacity as a parameter in the config file

configs/README.md

@@ -16,10 +16,11 @@ A detailed description of the parameters for model configuration (i.e., initiali
 | endtime | double (scalar)| >0 | sec, min, hr, d | simulation duration | - | time at which model simulation ends |
 | layer_soil_type | int (1D array) | - | - | state variable | - | layer soil type (read from the database file soil_params_file) |
 | max_soil_types | int | >1 | - | - | - | maximum number of soil types read from the file soil_params_file (default is set to 15) |
-| wilting_point_psi | double (scalar) | - | cm | state variable | - | wilting point (the amount of water not available for plants) used in computing AET |
+| wilting_point_psi | double (scalar) | - | cm | state variable | - | wilting point (the amount of water not available for plants) used in computing AET. Suggested value is 15495.0 cm, corresponding to 15 atm. |
+| field_capacity_psi | double (scalar) | - | cm | state variable | - | capillary head corresponding to volumetric water content at which gravity drainage becomes slower, used in computing AET. Suggested value is 340.9 cm for most soils, corresponding to 1/3 atm, and 103.3 cm for sands, corresponding to 1/10 atm. |
 | use_closed_form_G | bool | true or false | - | - | - | determines whether the numeric integral or closed form for G is used; a value of true will use the closed form. This defaults to false. |
 | giuh_ordinates | double (1D array)| - | - | state parameter | - | GIUH ordinates (for giuh based surface runoff) |
 | verbosity | string | high, low, none | - | debugging | - | controls IO (screen outputs and writing to disk) |
 | sft_coupled | Boolean | true, false | - | model coupling | impacts hydraulic conductivity | couples LASAM to SFT. Coupling to SFT reduces hydraulic conducitivity, and hence infiltration, when soil is frozen|
 | soil_z | double (1D array) | - | cm | spatial resolution | - | vertical resolution of the soil column (computational domain of the SFT model) |
-| calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_m`, and `vg_alpha` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |
+| calib_params | Boolean | true, false | - | calibratable params flag | impacts soil properties | If set to true, soil `smcmax`, `smcmin`, `vg_n`, `vg_alpha`, `hydraulic_conductivity`, `field_capacity_psi`, and `ponded_depth_max` are calibrated. defualt is false. vg = van Genuchten, SMC= soil moisture content |

configs/config_lasam_Bushland.txt

@@ -11,4 +11,5 @@ use_closed_form_G=false
 layer_soil_type=16,17,18
 max_soil_types=25
 wilting_point_psi=15495.0[cm]
+field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03

configs/config_lasam_Phillipsburg.txt

@@ -11,4 +11,5 @@ use_closed_form_G=false
 layer_soil_type=13,14,15
 max_soil_types=25
 wilting_point_psi=15495.0[cm]
+field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03

configs/config_lasam_sft_ngen.txt

@@ -11,6 +11,7 @@ use_closed_form_G=false
 layer_soil_type=13,14,15
 max_soil_types=25
 wilting_point_psi=15495.0[cm]
+field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
 sft_coupled=true
 soil_z=10,20,30,40,50,60,70,80,90,100.0,110.,120,130.,140.,150.,160.,170.,180.,190.,200.0[cm]

include/all.hxx

@@ -3,7 +3,7 @@
 #define _ALL_HXX
 
 /*
-  authors : Ahmad Jan and Fred Ogden
+  authors : Ahmad Jan and Fred Ogden and Peter La Follette
   year    : 2022
   email   : [email protected]
   - This header file constains functions' definitions used in the lgar.cxx, bmi_lasam.cxx and in other files.
@@ -128,6 +128,7 @@ struct lgar_bmi_parameters
 					    depth from the surface in meters */
   double *frozen_factor;                 // frozen factor added to the hydraulic conductivity due to coupling to soil freeze-thaw
   double  wilting_point_psi_cm;          // wilting point (the amount of water not available for plants or not accessible by plants)
+  double  field_capacity_psi_cm;          // field capacity represented as a capillary head. Note that both wilting point and field capacity are specified for the whole model domain with single values
   bool   use_closed_form_G = false;      /* true if closed form of capillary drive calculation is desired, false if numeric integral
 					    for capillary drive calculation is desired */
   double ponded_depth_cm;                // amount of water on the surface unavailable for surface runoff
@@ -186,11 +187,14 @@ struct lgar_mass_balance_variables
 // the structure holds pointer bmi output variables
 struct lgar_calib_parameters
 {
-  double *theta_e;     // theta_e = smcmax
-  double *theta_r;     // theta_r = smcmin
-  double *vg_alpha;    // Van Genuchton alpha
-  double *vg_m;        // Van Genuchton m
-  double *Ksat;        // Hydraulic conductivity [cm/hr]
+  double *theta_e;               // theta_e = smcmax [-]
+  double *theta_r;               // theta_r = smcmin [-]
+  double *vg_n;                  // Van Genuchten n [-]
+  double *vg_alpha;              // Van Genuchten alpha [1/cm]
+  double *Ksat;                  // Hydraulic conductivity [cm/hr]
+  double field_capacity_psi;    // field capacity in capillary head [cm]
+  double ponded_depth_max;      // maximum ponded depth of surface water [cm]
+
 };
 
 // nested structure of structures; main structure for the use in bmi
@@ -332,7 +336,7 @@ extern void InitializeWettingFronts(int num_layers, double initial_psi_cm, int *
 /*Other function prototypes for doing hydrology calculations, etc.  */
 /********************************************************************/
 
-extern double calc_aet(double PET_timestep_cm, double timestep_h, double wilting_point_psi_cm, int *soil_type,
+extern double calc_aet(double PET_timestep_cm, double timestep_h, double wilting_point_psi_cm, double field_capacity_psi_cm, int *soil_type,
 		       double AET_thresh_Theta, double AET_expon, struct wetting_front* head, struct soil_properties_ *soil_props);
 
 /********************************************************************/

include/bmi_lgar.hxx

@@ -67,9 +67,11 @@ public:
     // calibratable parameters
     this->calib_var_names[0] = "smcmax";
     this->calib_var_names[1] = "smcmin";
-    this->calib_var_names[2] = "van_genuchten_m";
+    this->calib_var_names[2] = "van_genuchten_n";
     this->calib_var_names[3] = "van_genuchten_alpha";
     this->calib_var_names[4] = "hydraulic_conductivity";
+    this->calib_var_names[5] = "field_capacity";
+    this->calib_var_names[6] = "ponded_depth_max";
   };
 
   void Initialize(std::string config_file);
@@ -132,7 +134,7 @@ private:
   struct model_state* state;
   static const int input_var_name_count  = 3;
   static const int output_var_name_count = 15;
-  static const int calib_var_name_count  = 5;
+  static const int calib_var_name_count  = 7;
 
   std::string input_var_names[input_var_name_count];
   std::string output_var_names[output_var_name_count];

realizations/realization_config_lasam.json

@@ -53,6 +53,15 @@
 					"precipitation_rate" : "P",
 					"potential_evapotranspiration_rate" : "PET"
 				    },
+					"model_params": {
+						"smcmax": [0.50, 0.49, 0.48],
+						"smcmin": [0.09, 0.08, 0.07],
+						"van_genuchten_n": [2.0, 2.1, 2.2],
+						"van_genuchten_alpha": [0.009, 0.008, 0.007],
+						"hydraulic_conductivity": [10.0, 11.0, 12.0],
+						"field_capacity": 333.3,
+						"ponded_depth_max": 1.1
+					},
 				    "output_variables" : [
 					"precipitation",
 					"potential_evapotranspiration",

realizations/realization_config_lasam_sft.json

@@ -88,6 +88,15 @@
 					"precipitation_rate" : "P",
 					"potential_evapotranspiration_rate" : "PET"
 				    },
+					"model_params": {
+						"smcmax": [0.50, 0.49, 0.48],
+						"smcmin": [0.09, 0.08, 0.07],
+						"van_genuchten_n": [2.0, 2.1, 2.2],
+						"van_genuchten_alpha": [0.009, 0.008, 0.007],
+						"hydraulic_conductivity": [10.0, 11.0, 12.0],
+						"field_capacity": 333.3,
+						"ponded_depth_max": 1.1
+					},
 				    "output_variables" : [
 					"precipitation",
 					"potential_evapotranspiration",

realizations/realization_config_lasam_smp.json

@@ -74,6 +74,15 @@
 					"precipitation_rate" : "P",
 					"potential_evapotranspiration_rate" : "PET"
 				    },
+					"model_params": {
+						"smcmax": [0.50, 0.49, 0.48],
+						"smcmin": [0.09, 0.08, 0.07],
+						"van_genuchten_n": [2.0, 2.1, 2.2],
+						"van_genuchten_alpha": [0.009, 0.008, 0.007],
+						"hydraulic_conductivity": [10.0, 11.0, 12.0],
+						"field_capacity": 333.3,
+						"ponded_depth_max": 1.1
+					},
 				    "output_variables" : [
 					"precipitation",
 					"potential_evapotranspiration",

src/aet.cxx

@@ -4,7 +4,7 @@
 #include "../include/all.hxx"
 
 //################################################################################
-/* authors : Fred Ogden and Ahmad Jan
+/* authors : Fred Ogden and Ahmad Jan and Peter La Follette
    year    : 2022
    the code computes actual evapotranspiration given PET.
    It uses an S-shaped function used in HYDRUS-1D (Simunek & Sejna, 2018).
@@ -14,7 +14,7 @@
 //################################################################################
 
 
-extern double calc_aet(double PET_timestep_cm, double time_step_h, double wilting_point_psi_cm,
+extern double calc_aet(double PET_timestep_cm, double time_step_h, double wilting_point_psi_cm, double field_capacity_psi_cm,
 		       int *soil_type, double AET_thresh_Theta, double AET_expon,
 		       struct wetting_front* head, struct soil_properties_ *soil_properties)
 {
@@ -44,8 +44,8 @@ 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 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 head_at_which_PET_equals_AET_cm = field_capacity_psi_cm; //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, which would be 103.3 cm.
   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);

src/bmi_lgar.cxx

@@ -115,6 +115,7 @@ Update()
   double subtimestep_h = state->lgar_bmi_params.timestep_h;
   int nint = state->lgar_bmi_params.nint;
   double wilting_point_psi_cm = state->lgar_bmi_params.wilting_point_psi_cm;
+  double field_capacity_psi_cm = state->lgar_bmi_params.field_capacity_psi_cm;
   bool use_closed_form_G = state->lgar_bmi_params.use_closed_form_G; 
 
   // constant value used in the AET function
@@ -192,7 +193,7 @@ Update()
 
     // Calculate AET from PET if PET is non-zero
     if (PET_subtimestep_cm_per_h > 0.0) {
-      AET_subtimestep_cm = calc_aet(PET_subtimestep_cm_per_h, subtimestep_h, wilting_point_psi_cm,
+      AET_subtimestep_cm = calc_aet(PET_subtimestep_cm_per_h, subtimestep_h, wilting_point_psi_cm, field_capacity_psi_cm,
                                     state->lgar_bmi_params.layer_soil_type, AET_thresh_Theta, AET_expon,
                                     state->head, state->soil_properties);
     }
@@ -509,27 +510,27 @@ update_calibratable_parameters()
 
   double volstart_before = lgar_calc_mass_bal(state->lgar_bmi_params.cum_layer_thickness_cm, state->head);
 
-  for (int i=0; i<state->lgar_bmi_params.num_wetting_fronts; i++) {
+  for (int i=0; i<state->lgar_bmi_params.num_wetting_fronts; i++) {//first we update the parameters that depend on soil layer, for each layer
     layer_num  = current->layer_num;
     soil = state->lgar_bmi_params.layer_soil_type[layer_num];
 
     assert (current != NULL);
 
     if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) {
-      std::cerr<<"----------- Calibratable parameters (initial values) ----------- \n";
+      std::cerr<<"----------- Calibratable parameters depending on soil layer (initial values) ----------- \n";
       std::cerr<<"| soil_type = "<< soil <<", layer = "<<layer_num
 	       <<", smcmax = "   << state->soil_properties[soil].theta_e
 	       <<", smcmin = "   << state->soil_properties[soil].theta_r
-	       <<", vg_m = "     << state->soil_properties[soil].vg_m
+	       <<", vg_n = "     << state->soil_properties[soil].vg_n
 	       <<", vg_alpha = " << state->soil_properties[soil].vg_alpha_per_cm
 	       <<", Ksat = "     << state->soil_properties[soil].Ksat_cm_per_h
-	       <<", theta = "   << current->theta <<"\n";
+	       <<", theta = "    << current->theta <<"\n";
     }
 
     state->soil_properties[soil].theta_e = state->lgar_calib_params.theta_e[layer_num-1];
     state->soil_properties[soil].theta_r = state->lgar_calib_params.theta_r[layer_num-1];
-    state->soil_properties[soil].vg_m    = state->lgar_calib_params.vg_m[layer_num-1];
-    state->soil_properties[soil].vg_n    = 1.0/(1.0 - state->soil_properties[soil].vg_m);
+    state->soil_properties[soil].vg_n    = state->lgar_calib_params.vg_n[layer_num-1];
+    state->soil_properties[soil].vg_m    = 1.0 - 1.0/state->soil_properties[soil].vg_n;
     state->soil_properties[soil].vg_alpha_per_cm = state->lgar_calib_params.vg_alpha[layer_num-1];
     state->soil_properties[soil].Ksat_cm_per_h   = state->lgar_calib_params.Ksat[layer_num-1];
 
@@ -538,18 +539,34 @@ update_calibratable_parameters()
 				       state->soil_properties[soil].theta_e, state->soil_properties[soil].theta_r);
 
     if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) {
-      std::cerr<<"----------- Calibratable parameters (updated values) ----------- \n";
+      std::cerr<<"----------- Calibratable parameters depending on soil layer (updated values) ----------- \n";
       std::cerr<<"| soil_type = "<< soil <<", layer = "<<layer_num
 	       <<", smcmax = "   << state->soil_properties[soil].theta_e
 	       <<", smcmin = "   << state->soil_properties[soil].theta_r
-	       <<", vg_m = "     << state->soil_properties[soil].vg_m
+	       <<", vg_n = "     << state->soil_properties[soil].vg_n
 	       <<", vg_alpha = " << state->soil_properties[soil].vg_alpha_per_cm
 	       <<", Ksat = "     << state->soil_properties[soil].Ksat_cm_per_h
-	       <<", theta = "   << current->theta <<"\n";
+	       <<", theta = "    << current->theta <<"\n";
     }
 
     current = current->next;
   }
+
+  //next we update the parameters that apply to the whole model domain and do not depend on soil layer
+  if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) {
+    std::cerr<<"----------- Calibratable parameters independent of soil layer (initial values) ----------- \n";
+    std::cerr<<"field_capacity_psi = "   << state->lgar_bmi_params.field_capacity_psi_cm
+      <<", ponded_depth_max = "     << state->lgar_bmi_params.ponded_depth_max_cm <<"\n";
+  }
+
+  state->lgar_bmi_params.field_capacity_psi_cm = state->lgar_calib_params.field_capacity_psi;
+  state->lgar_bmi_params.ponded_depth_max_cm   = state->lgar_calib_params.ponded_depth_max;
+
+  if (verbosity.compare("high") == 0 || verbosity.compare("low") == 0) {
+    std::cerr<<"----------- Calibratable parameters independent of soil layer (updated values) ----------- \n";
+    std::cerr<<"field_capacity_psi = "   << state->lgar_bmi_params.field_capacity_psi_cm
+      <<", ponded_depth_max = "     << state->lgar_bmi_params.ponded_depth_max_cm <<"\n";
+  }
 
   if (verbosity.compare("high") == 0)
     listPrint(state->head);
@@ -581,15 +598,15 @@ GetVarGrid(std::string name)
 	   || name.compare("actual_evapotranspiration") == 0) // double
     return 1;
   else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0
-	   || name.compare("soil_storage") == 0) // double
+	   || name.compare("soil_storage") == 0 || name.compare("field_capacity") == 0 || name.compare("ponded_depth_max") == 0)// double
     return 1;
   else if (name.compare("total_discharge") == 0 || name.compare("infiltration") == 0
 	   || name.compare("percolation") == 0  || name.compare("groundwater_to_stream_recharge") == 0) // double
     return 1;
   else if (name.compare("mass_balance") == 0)
     return 1;
   else if (name.compare("soil_depth_layers") == 0  || name.compare("smcmax") == 0 || name.compare("smcmin") == 0
-	   || name.compare("van_genuchten_m") == 0 || name.compare("van_genuchten_alpha") == 0
+	   || name.compare("van_genuchten_m") == 0 || name.compare("van_genuchten_alpha") == 0 || name.compare("van_genuchten_n") == 0 
 	   || name.compare("hydraulic_conductivity") == 0) // array of doubles (fixed length)
     return 2;
   else if (name.compare("soil_moisture_wetting_fronts") == 0 || name.compare("soil_depth_wetting_fronts") == 0) // array of doubles (dynamic length)
@@ -806,12 +823,16 @@ GetValuePtr (std::string name)
     return (void*)this->state->lgar_calib_params.theta_e;
   else if (name.compare("smcmin") == 0)
     return (void*)this->state->lgar_calib_params.theta_r;
-  else if (name.compare("van_genuchten_m") == 0)
-    return (void*)this->state->lgar_calib_params.vg_m;
+  else if (name.compare("van_genuchten_n") == 0)
+    return (void*)this->state->lgar_calib_params.vg_n;
   else if (name.compare("van_genuchten_alpha") == 0)
     return (void*)this->state->lgar_calib_params.vg_alpha;
   else if (name.compare("hydraulic_conductivity") == 0)
     return (void*)this->state->lgar_calib_params.Ksat;
+  else if (name.compare("ponded_depth_max") == 0)
+    return (void*)&this->state->lgar_calib_params.ponded_depth_max;
+  else if (name.compare("field_capacity") == 0)
+    return (void*)&this->state->lgar_calib_params.field_capacity_psi;
   else {
     std::stringstream errMsg;
     errMsg << "variable "<< name << " does not exist";