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

…le, in the same style as wilting point. I also updated the calibratable params list in bmi_lgar.hxx to include field capacaity, as well as max ponded head. I also edited that list to not include vg_m and instead use vg_n, because the parameters file provides n rather than m.

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
@@ -188,8 +189,8 @@ 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 *vg_alpha;    // Van Genuchten alpha
+  double *vg_m;        // Van Genuchten m
   double *Ksat;        // Hydraulic conductivity [cm/hr]
 };
 
@@ -332,7 +333,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];

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);
     }

src/lgar.cxx

@@ -159,6 +159,7 @@ extern void lgar_initialize(string config_file, struct model_state *state)
   @param frozen_factor          : frozen factor causing the hydraulic conductivity to decrease due to frozen soil
                                   (when coupled to soil freeze thaw model)
   @param wilting_point_psi_cm   : wilting point (the amount of water not available for plants or not accessible by plants)
+  @param field_capacity_psi_cm  : field capacity, represented with a capillary head (head above which drainage is much faster)
   @param ponded_depth_cm        : amount of water on the surface not available for surface drainage (initialized to zero)
   @param ponded_depth_max cm    : maximum amount of water on the surface not available for surface drainage (default is zero)
   @param nint                   : number of trapezoids used in integrating the Geff function (set to 120)
@@ -210,18 +211,19 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
   state->lgar_bmi_params.sft_coupled       = false;
   state->lgar_bmi_params.use_closed_form_G = false;
 
-  bool is_layer_thickness_set      = false;
-  bool is_initial_psi_set          = false;
-  bool is_timestep_set             = false;
-  bool is_endtime_set              = false;
-  bool is_forcing_resolution_set   = false;
-  bool is_layer_soil_type_set      = false;
-  bool is_wilting_point_psi_cm_set = false;
-  bool is_soil_params_file_set     = false;
-  bool is_max_soil_types_set       = false;
-  bool is_giuh_ordinates_set       = false;
-  bool is_soil_z_set               = false;
-  bool is_ponded_depth_max_cm_set  = false;
+  bool is_layer_thickness_set       = false;
+  bool is_initial_psi_set           = false;
+  bool is_timestep_set              = false;
+  bool is_endtime_set               = false;
+  bool is_forcing_resolution_set    = false;
+  bool is_layer_soil_type_set       = false;
+  bool is_wilting_point_psi_cm_set  = false;
+  bool is_field_capacity_psi_cm_set = false;
+  bool is_soil_params_file_set      = false;
+  bool is_max_soil_types_set        = false;
+  bool is_giuh_ordinates_set        = false;
+  bool is_soil_z_set                = false;
+  bool is_ponded_depth_max_cm_set   = false;
 
   string soil_params_file;
 
@@ -368,6 +370,17 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
       continue;
     }
+    else if (param_key == "field_capacity_psi") {
+      state->lgar_bmi_params.field_capacity_psi_cm = stod(param_value);
+      is_field_capacity_psi_cm_set = true;
+
+      if (verbosity.compare("high") == 0) {
+	std::cerr<<"Field capacity Psi [cm] : "<<state->lgar_bmi_params.field_capacity_psi_cm<<"\n";
+	std::cerr<<"          *****         \n";
+      }
+
+      continue;
+    }
     else if (param_key == "use_closed_form_G") { 
       if (param_value == "false") {
         state->lgar_bmi_params.use_closed_form_G = false;
@@ -522,6 +535,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
     state->soil_properties = new soil_properties_[state->lgar_bmi_params.num_soil_types+1];
     int num_soil_types = state->lgar_bmi_params.num_soil_types;
     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;
     int max_num_soil_in_file = lgar_read_vG_param_file(soil_params_file.c_str(), num_soil_types,
 						       wilting_point_psi_cm, state->soil_properties);
 
@@ -572,7 +586,13 @@ extern void InitFromConfigFile(string config_file, struct model_state *state)
 
   if(!is_wilting_point_psi_cm_set) {
     stringstream errMsg;
-    errMsg << "The configuration file \'" << config_file <<"\' does not set wilting_point_psi. \n";
+    errMsg << "The configuration file \'" << config_file <<"\' does not set wilting_point_psi. \n Recommended value of 15495.0[cm], corresponding to 15 atm. \n";
+    throw runtime_error(errMsg.str());
+  }
+
+  if(!is_field_capacity_psi_cm_set) {
+    stringstream errMsg;
+    errMsg << "The configuration file \'" << config_file <<"\' does not set field_capacity_psi. \n Recommended value of 340.9[cm] for most soils, corresponding to 1/3 atm, or 103.3[cm] for sands, corresponding to 1/10 atm. \n";
     throw runtime_error(errMsg.str());
   }
 

tests/configs/config_lasam_synth_0.txt

@@ -9,5 +9,6 @@ forcing_resolution=3600[sec]
 layer_soil_type=13,14,15
 max_soil_types=15
 wilting_point_psi=15495.0[cm]
+field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
 ponded_depth_max=1.0[cm]

tests/configs/config_lasam_synth_1.txt

@@ -9,4 +9,5 @@ forcing_resolution=300[sec]
 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.0,0.0,0.0,0.0,0.0

tests/configs/config_lasam_synth_2.txt

@@ -9,4 +9,5 @@ forcing_resolution=300[sec]
 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.0,0.0,0.0,0.0,0.0

tests/configs/unittest.txt

@@ -9,5 +9,6 @@ forcing_resolution=3600[sec]
 layer_soil_type=13,14,15
 max_soil_types=15
 wilting_point_psi=15495.0[cm]
+field_capacity_psi=340.9[cm]
 giuh_ordinates=0.06,0.51,0.28,0.12,0.03
 calib_params=true