add interflow

README.md

@@ -1,7 +1,7 @@
 # Catchment Arid/Semi-arid Model (CASAM) for streamflow simulation
 CASAM (formerly known as LASAM) is a catchment scale hydrolgic model originally designed for arid or semi arid areas, in which the partitioning of precipitation into infiltration and runoff is important. CASAM is composed of a vadose zone model, which is Layered Green & Ampt with redistribution (LGAR) and represents a layered soil matrix, and a nonlinear reservoir which is filled via a simple formulation for preferential flow. LGAR is a model which partitions precipitation into infiltration and runoff, and is designed for use in arid or semi-arid climates. LGAR closely mimics precipitation partitioning results simulated by the Richards/Richardson equation (RRE), without the inherent reliability and stability challenges the RRE poses. Therefore, this model is useful when accurate, stable precipitation partitioning simulations are desired in arid or semi-arid areas. LGAR in Python (no longer supported) is available [here](https://github.com/NOAA-OWP/LGAR-Py).
 
-CASAM is theoretically a skillful catchment scale hydrolgic model in the event that precipitation partitioning into infiltration and runoff is the most important process for streamflow generation in a given catchment. Nonetheless, CASAM also includes a nonlinear reservoir so that some water stored in the catchment can contribute directly to streamflow (currently, we assume soil water does not directly contribute to streamflow).
+CASAM is theoretically a skillful catchment scale hydrolgic model in the event that precipitation partitioning into infiltration and runoff is the most important process for streamflow generation in a given catchment. Nonetheless, CASAM also includes a nonlinear reservoir and interflow so that some water stored in the catchment can contribute directly to streamflow.
 
 ![Alt text](CASAM_diagram.png)
 

data/README.md

@@ -9,13 +9,11 @@ It is recommended to calibrate:
 - `Ks` for the top two layers  
 - `n` for the top two layers  
 - `field_capacity_psi`  
-- `a`  
-- `b`  
+- `a_con_res`
+- `b_con_res`
 - `frac_to_CR`  
 - `spf_factor`
 
-Note that it is not necessarily recommended to calibrate maximum ponded head or theta_r.
-
 | Parameter name | Units |  Physical limits  | Range tested for stability | Applies to individual soil layers or entire model domain | Description |
 | --- | --- | --- | --- | --- | --------------- |
 | theta_r | - | 0< theta_r <1, <br> theta_r < theta_e | 0.01 < theta_r < 0.15 | Soil layer | theta_r is the residual water content, or the minimum volumetric water content that a soil layer can naturally attain. Note that theta_r must be less than than theta_e. This is set per soil layer, in the .dat file in the data directory.|
@@ -25,10 +23,12 @@ Note that it is not necessarily recommended to calibrate maximum ponded head or
 | Ks | cm/h | Ks > 0 | 0.0001 < K_s < 100 | Soil layer | Ks is the saturated hydraulic conductivity of a soil. Note that in nature, expected values of Ks are distributed logarithmically, so calibrating on the log of Ks rather than Ks directly is likely a better choice for most calibration algorithms. This is set per soil layer, in the .dat file in the data directory.|
 | ponded_head_max | cm | ponded_head_max >= 0 | 0 <= ponded_head_max <= 5 | Entire model domain | This is the maximum amount of ponded water that is allowed to accumulate on the soil surface. While stability tests have only included a maximum value of 5 cm, any value greater than or equal to 0 should be acceptable. A common choice will be 0. This parameter can be set in the config file.  |
 | field_capacity_psi | cm | 0 < field_capacity_psi, <br>field_capacity_psi < wilting_point_psi | 10 < field_capacity_psi < 500 | Entire model domain | This is the wilting point of the model domain, expressed as a capillary head. Together with wilting_point_psi, the field capacity is used to determine the intensity of the reduction of PET to become AET. The numbers 10.3 cm and 516.6 cm correspond to pressures of 1/100 atm and 1/2 atm of water. Note that the model generally uses absolute values of capillary head; in this case, these limits are absolute values of negative numbers and physically represent unsaturated soil. While field capacity will vary per soil type, we use a single value for the entire model domain, following the method for PET->AET correction used by HYDRUS. This parameter can be set in the config file. |
-| a | cm^(1-b) h^-1 | 0 < a | 1E-8 < a < 1E-1 | Applies to nonlinear reservoir | As of this writing (27 June 2025), CASAM fundamentally has two different types of water storage: water stored in the vadose zone that does not contribute to streamflow, and water stored in a nonlinear reservoir that does contribute to streamflow. The nonlinear reservoir has an input of simple bypass through the vadose zone (controlled by frac_to_CR and spf_factor). The nonlinear reservoir releases water to the stream at a rate of a*S^b, where S is the water stored in the reserovir in cm, and a and b are nonlinear reservoir parameters. Note that the units of a depend on the value of b. Defaults to 0. Also note that because a, b, and frac_to_CR all default to 0, all of these must be specified in the config file in order to enable the nonlinear reservoir. |
-| b | - | - | 0.01 < b < 5 | Applies to nonlinear reservoir | As of this writing (27 June 2025), CASAM fundamentally has two different types of water storage: water stored in the vadose zone that does not contribute to streamflow, and water stored in a nonlinear reservoir that does contribute to streamflow. The nonlinear reservoir has an input of simple bypass through the vadose zone (controlled by frac_to_CR and spf_factor). The nonlinear reservoir releases water to the stream at a rate of a*S^b, where S is the water stored in the reserovir in cm, and a and b are nonlinear reservoir parameters. Defaults to 0.  |
-| frac_to_CR | - | 0 <= frac_to_CR <= 1 | 1E-4 <= frac_to_CR <= 1 - 1E-4 | Applies to nonlinear reservoir | storage that contributes directly to streamflow | Simple bypass of water at the soil surface to the nonlinear conceptual reservoir will occur when the most superficial surface wetting front achieves the theta_e value of its layer times spf_factor. When this occurs, the amount of water sent to the nonlinear reservoir is equal to the precipitation plus any ponded water times frac_to_CR. This is a rather simple representation of preferential flow that intends to simulate the episodic nature of streamflow events in arid or semi arid environments. Note that either all or none of a, b, and frac_to_CR must be specified. If none are specified then the model will not simulate a nonlinear reservoir. Further, as of this writing (27 June 2025), fluxes from the lower boundary of the vadose zone do not contribute to the nonlinear reservoir and are said to be lost to deep GW. These fluxes should be rare in arid or semi arid environments and represent a wetting front partially crossing the model lower boundary, so the fraction that did technically contributes to fluxes through LGAR's lower boundary. Defaults to 0. Note that if a second reservoir is desired, a_slow, b_slow, and frac_slow can all be specified with the same ranges as a, b, and frac_to_CR, respectively, with the expection that frac_slow can not equal exactly 0.|
-| spf_factor | - | 0.0 <= spf_factor <= 1 | 0.1 <= spf_factor <= 1 | Applies to nonlinear reservoir | storage that contributes directly to streamflow | Simple bypass of surface water to the nonlinear reservoir will occur when the most superficial wetting front achieves the theta_e value of its layer times spf_factor. When this occurs, the amount of water sent to the nonlinear reservoir is equal to the precipitation plus any ponded water times frac_to_CR. This is a rather simple representation of preferential flow that intends to simulate the episodic nature of streamflow events in arid or semi arid environments. Defaults to 0.98. Specified in the config file. |
+| a_con_res | cm^(1-b_con_res) h^-1 | 0 < a_con_res | 1E-8 < a_con_res < 1E-1 | Applies to nonlinear reservoir | The nonlinear reservoir is one route by which catchment water storage contributes to streamflow. The nonlinear reservoir has an input of simple bypass through the vadose zone, controlled by frac_to_CR and spf_factor, and optionally free drainage if `free_drainage_to_CR=true` in the config file. The nonlinear reservoir releases water to the stream at a rate of a_con_res*S^b_con_res, where S is the water stored in the reservoir in cm, and a_con_res and b_con_res are nonlinear reservoir parameters. Note that the units of a_con_res depend on the value of b_con_res. Defaults to 0. Also note that because a_con_res, b_con_res, and frac_to_CR all default to 0, all of these must be specified in the config file in order to enable the nonlinear reservoir. The legacy name `a` is still accepted. |
+| b_con_res | - | - | 0.01 < b_con_res < 5 | Applies to nonlinear reservoir | The nonlinear reservoir is one route by which catchment water storage contributes to streamflow. The nonlinear reservoir has an input of simple bypass through the vadose zone, controlled by frac_to_CR and spf_factor, and optionally free drainage if `free_drainage_to_CR=true` in the config file. The nonlinear reservoir releases water to the stream at a rate of a_con_res*S^b_con_res, where S is the water stored in the reservoir in cm, and a_con_res and b_con_res are nonlinear reservoir parameters. Defaults to 0. The legacy name `b` is still accepted. |
+| frac_to_CR | - | 0 <= frac_to_CR <= 1 | 1E-4 <= frac_to_CR <= 1 - 1E-4 | Applies to nonlinear reservoir | Simple bypass of water at the soil surface to the nonlinear conceptual reservoir will occur when the most superficial surface wetting front achieves the theta_e value of its layer times spf_factor. When this occurs, the amount of water sent to the nonlinear reservoir is equal to the precipitation plus any ponded water times frac_to_CR. This is a simple representation of preferential flow that intends to simulate the episodic nature of streamflow events in arid or semi arid environments. Note that either all or none of a_con_res, b_con_res, and frac_to_CR must be specified. If none are specified then the model will not simulate a nonlinear reservoir. Free drainage contributes to the nonlinear reservoir only when `free_drainage_enabled=true` and `free_drainage_to_CR=true`; otherwise lower boundary free drainage is accounted as percolation/recharge leaving the LGAR domain. Defaults to 0. Note that if a second reservoir is desired, a_con_res_slow, b_con_res_slow, and frac_slow can all be specified with the same ranges as a_con_res, b_con_res, and frac_to_CR, respectively, with the exception that frac_slow can not equal exactly 0. The legacy names `a_slow` and `b_slow` are still accepted.|
+| spf_factor | - | 0.0 <= spf_factor <= 1 | 0.1 <= spf_factor <= 1 | Applies to nonlinear reservoir | Simple bypass of surface water to the nonlinear reservoir will occur when the most superficial wetting front achieves the theta_e value of its layer times spf_factor. When this occurs, the amount of water sent to the nonlinear reservoir is equal to the precipitation plus any ponded water times frac_to_CR. This is a simple representation of preferential flow that intends to simulate the episodic nature of streamflow events in arid or semi arid environments. Defaults to 0.98. Specified in the config file. |
+| interflow_psi_threshold | cm | interflow_psi_threshold >= 0 | 10000 >= interflow_psi_threshold >= 1 | Entire model domain | Subsurface interflow from LGAR wetting fronts is enabled only when both interflow_psi_threshold and interflow_factor are specified. A wetting front can contribute interflow when its capillary head is less than or equal to this threshold. LGAR uses positive absolute capillary head values, so smaller psi values represent wetter conditions. If log_mode is on, this parameter is interpreted as a log10 value. Defaults to disabled. The legacy name `lateral_flow_psi_threshold` is still accepted. While the tested limits are large, a sensible calibration range would probably not let this go much above the field capacity. Suggested calibration upper limit is 1000 cm.|
+| interflow_factor | - | interflow_factor >= 0 | 1.E-4 <= interflow_factor <= 1.E4 | Entire model domain | Subsurface interflow from LGAR wetting fronts is enabled only when both interflow_psi_threshold and interflow_factor are specified. Interflow for an eligible wetting front is proportional to K(theta) times interflow_factor, scaled by the fraction of the LGAR domain represented by that front. Interflow is removed from LGAR storage and routed through the GIUH with surface runoff, so water stored in the LGAR domain can conditionally contribute to streamflow when this process is enabled. If log_mode is on, this parameter is interpreted as a log10 value. Defaults to disabled. The legacy name `lateral_flow_factor` is still accepted. While the tested physical limits are large, a sensible calibration range is probably 1.E-3 - 1.E2.|
 
 Parameters that are specified per soil layer can either be scalar values with double precision (in the event the model is run with 1 layer) or vectors of doubles (in the event that the model is run with more than 1 layer), whereas parameters that are specified for the entire model domain are scalars with double precision. 
 
@@ -52,4 +52,3 @@ Below is a table of parameters for soils from the HYDRUS soils catalog, which ca
 | Sandy Clay          | 0.1   | 0.38 | 0.027 | 1.23 | 0.12        |
 | Silty Clay          | 0.07  | 0.36 | 0.005 | 1.09 | 0.02        |
 | Clay                | 0.068 | 0.38 | 0.008 | 1.09 | 0.2         |
-

include/all.hxx

@@ -139,12 +139,15 @@ struct lgar_bmi_parameters
   double mbal_tol;                       // if a substep's mass balance error is larger than this number, the model will abort. By default it is set to a large value (10 cm).
   double ponded_depth_cm;                // amount of water on the surface unavailable for surface runoff
   double ponded_depth_max_cm;            // maximum amount of water on the surface unavailable for surface runoff
-  double a = 0.0;                        // parameter for nonlinear reservoir
-  double b = 0.0;                        // parameter for nonlinear reservoir
+  double a_con_res = 0.0;                // parameter for nonlinear reservoir
+  double b_con_res = 0.0;                // parameter for nonlinear reservoir
   double frac_to_CR = 0.0;               // parameter for nonlinear reservoir 
-  double a_slow = 0.0;                        // parameter for nonlinear reservoir
-  double b_slow = 0.0;                        // parameter for nonlinear reservoir
+  double a_con_res_slow = 0.0;           // parameter for nonlinear reservoir
+  double b_con_res_slow = 0.0;           // parameter for nonlinear reservoir
   double frac_slow = 0.0;               // parameter for nonlinear reservoir 
+  bool   interflow_enabled = false;        // if true, wetting fronts can contribute interflow to the GIUH queue
+  double interflow_psi_threshold_cm = 0.0; // interflow is zero when a wetting front's psi is above this threshold [cm]
+  double interflow_factor = 0.0;           // multiplier applied to K(theta) to calculate wetting front interflow
   double spf_factor = 0.98;              // parameter that controls the theta value above which contributions to the nonlinear reservoir will be made
   double precip_previous_timestep_cm;    // amount of rainfall (previous time step)
 
@@ -183,6 +186,7 @@ struct lgar_mass_balance_variables
   double volAET_timestep_cm;         // volume of AET at each timestep
   double volPET_timestep_cm;         // volume of PET at each timestep
   double volrech_timestep_cm;        // volume of water leaving soil to the ground water (ground water recharge)
+  double volinterflow_timestep_cm;// volume of wetting front interflow at each timestep
   double volrunoff_giuh_timestep_cm; // volume of giuh runoff at each timestep
   double volQ_timestep_cm;           // total outgoing water (surface runoff + water from conceptual reservoirs, both of which go through GIUH)
   double volQ_CR_timestep_cm;        // outgoing water just from conceptual reservoirs
@@ -211,6 +215,7 @@ struct lgar_mass_balance_variables
   double accumulated_free_drainage = 0.0;
 
   double volrech_cm;          // volume of water leaving soil through the bottom of the domain (ground water recharge)
+  double volinterflow_cm;  // volume of wetting front interflow routed through GIUH
   double volrunoff_giuh_cm;   // volume of giuh runoff
   double volQ_cm;             // total outgoing water
   double volQ_CR_cm;          // water outgoing just from conceptual reservoirs
@@ -250,12 +255,14 @@ struct lgar_calib_parameters
 
   double field_capacity_psi;    // field capacity in capillary head [cm]
   double ponded_depth_max;      // maximum ponded depth of surface water [cm]
-  double a;                      // parameter for nonlinear reservoir
-  double b;                      // parameter for nonlinear reservoir
+  double a_con_res;              // parameter for nonlinear reservoir
+  double b_con_res;              // parameter for nonlinear reservoir
   double frac_to_CR;             // parameter for nonlinear reservoir
-  double a_slow;                      // parameter for nonlinear reservoir
-  double b_slow;                      // parameter for nonlinear reservoir
+  double a_con_res_slow;         // parameter for nonlinear reservoir
+  double b_con_res_slow;         // parameter for nonlinear reservoir
   double frac_slow;             // parameter for nonlinear reservoir
+  double interflow_psi_threshold_cm; // interflow is zero when a wetting front's psi is above this threshold [cm]
+  double interflow_factor;    // multiplier applied to K(theta) to calculate wetting front interflow
   double spf_factor;             // parameter for nonlinear reservoir
 
 };
@@ -352,7 +359,8 @@ extern double lgar_insert_water(bool use_closed_form_G, int nint, double timeste
 				double *frozen_factor, struct wetting_front* head, struct soil_properties_ *soil_properties);
 
 // the subroutine moves wetting fronts, merges wetting fronts, and does the mass balance correction if needed
-extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_subtimestep_cm, double *ponded_depth_cm, int wf_free_drainage_demand,
+extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_subtimestep_cm, double *interflow_subtimestep_cm, double interflow_psi_threshold_cm,
+				     double interflow_factor, double *ponded_depth_cm, int wf_free_drainage_demand,
 				     double old_mass, double mass_correction_for_cached_free_drainage_fluxes, int number_of_layers, double *actual_ET_demand,
 				     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);
@@ -450,8 +458,8 @@ extern bool is_epsilon_less_than(double a, double eps);
 //function for contribtion to streamflow from conceptual reservoir
 extern double calc_CR_Q(
     double subtimestep_h,
-    double a_fast, double a_slow,
-    double b_fast, double b_slow,
+    double a_con_res, double a_con_res_slow,
+    double b_con_res, double b_con_res_slow,
     double frac_slow,  // fraction (0 - 1) of recharge going to slow reservoir
     double precip_for_CR_subtimestep_cm_per_h,
     double *CR_fast_storage_cm,