diff --git a/CASAM_diagram.png b/CASAM_diagram.png index dad16ff..9e4cb38 100644 Binary files a/CASAM_diagram.png and b/CASAM_diagram.png differ diff --git a/README.md b/README.md index 95964bf..10fd1b3 100644 --- a/README.md +++ b/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) diff --git a/configs/README.md b/configs/README.md index 5bdae7c..ef6c8a9 100644 --- a/configs/README.md +++ b/configs/README.md @@ -3,6 +3,10 @@ Example configuration files are provided in this directory. To build and run the A detailed description of the parameters for model configuration (i.e., initialize/setup) is provided below. +CASAM streamflow can be produced by surface runoff routed through the GIUH, by optional nonlinear conceptual reservoir outflow, and by optional subsurface interflow from LGAR wetting fronts. When interflow is disabled, water stored in the LGAR vadose zone domain does not directly contribute to streamflow except through configured lower boundary or bypass pathways. When interflow is enabled, eligible wetting fronts can lose water as subsurface interflow, and that interflow volume is routed through the GIUH with surface runoff. + +The interflow parameters `interflow_psi_threshold` and `interflow_factor` are model wide calibratable parameters exposed through BMI with those names; they are not soil layer specific. There is no separate boolean switch for interflow. Both parameters must be specified together to enable interflow, and omitting both disables it. Legacy config and BMI names `lateral_flow_psi_threshold` and `lateral_flow_factor` are still accepted. + | Variable | Datatype | Limits | Units | Role | Process | Description | | -------- | -------- | ------ | ----- | ---- | ------- | ----------- | @@ -27,14 +31,16 @@ A detailed description of the parameters for model configuration (i.e., initiali | adaptive_timestep | Boolean | true, false | - | adaptive timestep flag | impacts timestep | If set to true, LGAR will use an internal adaptive timestep, and the above timestep is used as a minimum timestep (recommended value of 300 seconds). The adaptive timestep will never be larger than the forcing resolution. If set to false, LGAR will use the above specified timestep as a fixed timestep. Testing indicates that setting this value to true substantially decreases runtime while negligibly changing the simulation. We recommend this to be set to true. | | free_drainage_enabled | Boolean | true, false | - | controls lower boundary condition | affects recharge | If free_drainage_enabled is true, then free drainage will be enabled as the lower boundary condition, where fluxes from the vadose zone to groundwater are controlled by the hydraulic conductivity at the bottom of the model domain. If free_drainage_enabled is set to false, then the lower boundary condition will be no flow. If LGAR is used in an area with substantially more PET than precipitation, the choice of no flow vs free drainage should be less impactful, because the majority of water that leaves the vadose zone will do so as AET. Defaults to false.| | free_drainage_to_CR | Boolean | true, false | - | controls lower boundary condition | affects recharge | If this is set to true, then free drainage water will contribute to the conceptual reservoir. Defaults to false.| +| interflow_psi_threshold | double (scalar) | >=0 | cm | model wide calibratable parameter for wetting front interflow | vadose storage that contributes directly to streamflow | 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. The legacy config and BMI name `lateral_flow_psi_threshold` is still accepted. If log_mode is on, this parameter is interpreted as a log10 value. Defaults to disabled. A practical calibration range should reflect the intended wetness activation threshold, for example near field capacity if only relatively wet fronts should contribute. Suggested calibration upper limit is 1000 cm.| +| interflow_factor | double (scalar) | 1.E-4 <= interflow_factor <= 1.E4 | - | model wide calibratable parameter for wetting front interflow | vadose storage that contributes directly to streamflow | Candidate 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. This factor can be interpreted as an effective subsurface conductance multiplier, with effects such as horizontal to vertical conductivity anisotropy, unresolved hillslope geometry, and catchment connectivity folded into one scalar. Interflow is removed from LGAR storage and routed through the GIUH with surface runoff. If log_mode is on, this parameter is interpreted as a log10 value. Defaults to disabled. A suggested calibration range is 1.E-3 <= interflow_factor <= 1.E2, narrower than the full supported tested range. The legacy config and BMI name `lateral_flow_factor` is still accepted.| | mbal_tol | double (scalar) | >0 | cm | mass balance error resulting from a substep that will trigger a model crash | mass balance accounting | If the mass balance error is greater than this number in a single substep, then the model will abort. Global mass balance errors over the course of year long simulations (i.e. the sum of all mass balance errors over the course of a long simulation) tend to be small, where a value greater than 1E-4 cm tends to be rare, occuring less than 1 in 10000 prameter sets. LGAR in theory should both be mass conservative, however the presence of unlikely but possible edge cases can cause mass balance errors. Flux caching can also cause small mass balance errors. While the model will usually converge with a small mbal_tol value, if general convergence is desired across a large number of parameter sets / forcing datasets, then we recommend that this value is not specified and therefore the default value of 1E1cm will be used. | | PET_affects_precip | Boolean | true, false | - | specifies whether PET is subtracted from precip | forcing data manipulation | If enabled, then PET will be subtracted from precipitation. Defaults to false.| -| a | double (scalar) | 1E-8 < a < 1E-1 | cm^(1-b) h^-1 | parameter for nonlinear reservoir | storage that contributes directly to streamflow | 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), and free drainage if desired. 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.| -| b | double (scalar) | 0.01 < b < 5 | - | parameter for nonlinear reservoir | storage that contributes directly to streamflow | 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), and free drainage if desired. 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 | double (scalar) | 0.0 <= frac_to_CR <= 1 | - | parameter for 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. Defaults to 0.| -| spf_factor | double (scalar) | 0.1 <= spf_factor <= 1 | - | parameter for fluxes 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. | -| allow_flux_caching | Boolean | true, false | - | trades a small amount of accuracy for a lot of speed | flux caching | During dry periods, it is often the case that wetting fronts will move very slowly and AET will be significantly less than PET. In these cases, in the context of streamflow simulation, it is not efficient to recompute fluxes and soil moisture dynamics for each time step. If this is set to true, then fluxes and wetting front movement will only be recomputed once every 24 hours, or when the conditions resulting in dry and slow wetting fronts and low AET cease. During the times for which fluxes are not recomputed, instead they are stored in a cache and fluxes for subsequent time steps are set using this cache. Sligtly different strategies are used for fluxes through the lower boundary and AET. Also note that because NextGen models should ideally provide output for each hour, simply setting an adaptive time step to be larger than one hour is not a preferred runtime reduction method here. Note that this can cause small mass balance errors when the lower boundary condition is set to free drainage. Defaults to false. | -| log_mode | Boolean | true, false | - | helps calibration search space exploration | log transform of parameters | When this is set to true, then all inputs for the van Genuchten parameter alpha, saturated hydraulic conductivity, and the nonlinear reservoir parameter a must be input as their log values rather than the normal values. For example, if an saturated hydraulic conductivity of 0.1 cm/h is desired, then the input value must be -1 because 10^-1 = 0.1. The reasoning for this is that these parameters are not distributed normally in nature but rather are distributed log normally, such that simply sampling the parameter space normally during calibration will vastly undersample a big region of the parameter space in which we expect useful parameter sets to be. Defaults to false. | -| a_slow | double (scalar) | 1E-8 < a_slow < 1E-1 | cm^(1-b) h^-1 | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter a, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0.| -| b_slow | double (scalar) | 0.01 < b_slow < 5 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter b, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0.| -| frac_slow | double (scalar) | 0.0 < frac_slow <= 1 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This describes the partitioning of water to the two reservoris, where the the input to the slow reservoir is equal to the total input for the nonlinear reservoirs times frac_slow. Note that either all or none of a_slow, b_slow, and frac_slow must be specified. If none are specified then the model will not simulate a second nonlinear reservoir. Defaults to 0.| \ No newline at end of file +| a_con_res | double (scalar) | 1E-8 < a_con_res < 1E-1 | cm^(1-b_con_res) h^-1 | parameter for nonlinear reservoir | storage that contributes directly to streamflow | The nonlinear reservoir is one route by which catchment water storage contributes to streamflow. Its input can include simple bypass through the vadose zone, controlled by frac_to_CR and spf_factor, and free drainage if desired. 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. The legacy name `a` is still accepted.| +| b_con_res | double (scalar) | 0.01 < b_con_res < 5 | - | parameter for nonlinear reservoir | storage that contributes directly to streamflow | The nonlinear reservoir is one route by which catchment water storage contributes to streamflow. Its input can include simple bypass through the vadose zone, controlled by frac_to_CR and spf_factor, and free drainage if desired. 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 | double (scalar) | 0.0 <= frac_to_CR <= 1 | - | parameter for 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_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. Defaults to 0.| +| spf_factor | double (scalar) | 0.1 <= spf_factor <= 1 | - | parameter for fluxes to nonlinear reservoir | storage that contributes directly to streamflow | Simple preferential flow (SPF) factor: 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. | +| allow_flux_caching | Boolean | true, false | - | trades a small amount of accuracy for a lot of speed | flux caching | During dry periods, it is often the case that wetting fronts will move very slowly and AET will be significantly less than PET. In these cases, in the context of streamflow simulation, it is not efficient to recompute fluxes and soil moisture dynamics for each time step. If this is set to true, then fluxes and wetting front movement will only be recomputed once every 24 hours, or when the conditions resulting in dry and slow wetting fronts and low AET cease. During the times for which fluxes are not recomputed, instead they are stored in a cache and fluxes for subsequent time steps are set using this cache. Sligtly different strategies are used for fluxes through the lower boundary and AET. Flux caching is disabled whenever interflow is enabled and at least one wetting front is eligible to contribute interflow. Also note that because NextGen models should ideally provide output for each hour, simply setting an adaptive time step to be larger than one hour is not a preferred runtime reduction method here. Note that this can cause small mass balance errors when the lower boundary condition is set to free drainage. Defaults to false. | +| log_mode | Boolean | true, false | - | helps calibration search space exploration | log transform of parameters | When this is set to true, then all inputs for the van Genuchten parameter alpha, saturated hydraulic conductivity, the nonlinear reservoir parameter a_con_res, interflow_psi_threshold, and interflow_factor must be input as their log10 values rather than the normal values. For example, if an saturated hydraulic conductivity of 0.1 cm/h is desired, then the input value must be -1 because 10^-1 = 0.1. The reasoning for this is that these parameters are not distributed normally in nature but rather are distributed log normally, such that simply sampling the parameter space normally during calibration will vastly undersample a big region of the parameter space in which we expect useful parameter sets to be. Defaults to false. | +| a_con_res_slow | double (scalar) | 1E-8 < a_con_res_slow < 1E-1 | cm^(1-b_con_res_slow) h^-1 | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter a_con_res, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0. The legacy name `a_slow` is still accepted.| +| b_con_res_slow | double (scalar) | 0.01 < b_con_res_slow < 5 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This is exactly like the parameter b_con_res, except it corresponds to a second nonlinear reservoir, which was added to simulate cases where receding limbs have behaviors that can not easily be captured by one reservoir. Defaults to 0. The legacy name `b_slow` is still accepted.| +| frac_slow | double (scalar) | 0.0 < frac_slow <= 1 | - | parameter for second nonlinear reservoir | storage that contributes directly to streamflow | This describes the partitioning of water to the two reservoris, where the the input to the slow reservoir is equal to the total input for the nonlinear reservoirs times frac_slow. Note that either all or none of a_con_res_slow, b_con_res_slow, and frac_slow must be specified. If none are specified then the model will not simulate a second nonlinear reservoir. Defaults to 0.| diff --git a/data/README.md b/data/README.md index 8084545..7a83311 100644 --- a/data/README.md +++ b/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,
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,
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 | - diff --git a/include/all.hxx b/include/all.hxx index aacbf6e..7976ba6 100755 --- a/include/all.hxx +++ b/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, diff --git a/include/bmi_lgar.hxx b/include/bmi_lgar.hxx index 4136791..f47cd78 100644 --- a/include/bmi_lgar.hxx +++ b/include/bmi_lgar.hxx @@ -84,10 +84,12 @@ public: this->calib_var_names[5] = "van_genuchten_alpha_2"; this->calib_var_names[6] = "hydraulic_conductivity_2"; this->calib_var_names[7] = "field_capacity"; - this->calib_var_names[8] = "a"; - this->calib_var_names[9] = "b"; + this->calib_var_names[8] = "a_con_res"; + this->calib_var_names[9] = "b_con_res"; this->calib_var_names[10] = "frac_to_CR"; - this->calib_var_names[11] = "spf_factor"; + this->calib_var_names[11] = "interflow_psi_threshold"; + this->calib_var_names[12] = "interflow_factor"; + this->calib_var_names[13] = "spf_factor"; }; @@ -153,7 +155,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 = 12; + static const int calib_var_name_count = 14; std::string input_var_names[input_var_name_count]; std::string output_var_names[output_var_name_count]; diff --git a/src/bmi_lgar.cxx b/src/bmi_lgar.cxx index e47930b..4613f59 100644 --- a/src/bmi_lgar.cxx +++ b/src/bmi_lgar.cxx @@ -56,6 +56,16 @@ string verbosity="none"; // small epsillon that is used to determine if the difference between two quantities is 0 while avoiding machine precision errors #define SMALL_EPS 1.E-12 +// checks if there is at least one WF contributing to interflow. If so then flux caching will be disabled. +static bool any_wetting_front_can_interflow(struct wetting_front *head, double interflow_psi_threshold_cm) +{ + for (struct wetting_front *current = head; current != NULL; current = current->next) { + if (current->psi_cm <= interflow_psi_threshold_cm) + return true; + } + return false; +} + /** * @brief Delete dynamic arrays allocated in Initialize() and held by this object @@ -151,6 +161,7 @@ Update() state->lgar_mass_balance.volCRend_cm = state->lgar_mass_balance.volCRstart_cm; state->lgar_mass_balance.volAET_cm = 0.0; state->lgar_mass_balance.volrech_cm = 0.0; + state->lgar_mass_balance.volinterflow_timestep_cm = 0.0; state->lgar_mass_balance.volrunoff_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; state->lgar_mass_balance.volQ_cm += state->lgar_bmi_input_params->precipitation_mm_per_h * mm_to_cm; state->lgar_mass_balance.volQ_CR_cm = 0.0; @@ -205,6 +216,7 @@ Update() double volon_timestep_cm = state->lgar_mass_balance.volon_timestep_cm; double volrunoff_timestep_cm = 0.0; double volrech_timestep_cm = 0.0; + double volinterflow_timestep_cm = 0.0; double volrunoff_giuh_timestep_cm = 0.0; double volQ_timestep_cm = 0.0; double volQ_CR_timestep_cm = 0.0; @@ -230,12 +242,14 @@ Update() 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; - double a = state->lgar_bmi_params.a; - double b = state->lgar_bmi_params.b; + double a_con_res = state->lgar_bmi_params.a_con_res; + double b_con_res = state->lgar_bmi_params.b_con_res; double frac_to_CR = state->lgar_bmi_params.frac_to_CR; - double a_slow = state->lgar_bmi_params.a_slow; - double b_slow = state->lgar_bmi_params.b_slow; + double a_con_res_slow = state->lgar_bmi_params.a_con_res_slow; + double b_con_res_slow = state->lgar_bmi_params.b_con_res_slow; double frac_slow = state->lgar_bmi_params.frac_slow; + double interflow_psi_threshold_cm = state->lgar_bmi_params.interflow_psi_threshold_cm; + double interflow_factor = state->lgar_bmi_params.interflow_factor; double spf_factor = state->lgar_bmi_params.spf_factor; bool use_closed_form_G = state->lgar_bmi_params.use_closed_form_G; bool adaptive_timestep = state->lgar_bmi_params.adaptive_timestep; @@ -306,6 +320,11 @@ Update() } } + if (state->lgar_bmi_params.interflow_enabled + && any_wetting_front_can_interflow(state->head, state->lgar_bmi_params.interflow_psi_threshold_cm)) { + state->lgar_mass_balance.cache_fluxes = FALSE; + } + if (caching_at_start && !state->lgar_mass_balance.cache_fluxes){ switch_caching = TRUE;//if you switch from cached to not, you need to add the "missing" PET back into the mass balance and AET calculation } @@ -403,6 +422,7 @@ Update() double temp_rch = 0.0; //handles case when a fraction of a wetting front technically crosses the lower boundary of the vadose zone double free_drainage_subtimestep_cm = 0.0; double free_drainage_for_CR = 0.0; + double interflow_subtimestep_cm = 0.0; PET_subtimestep_cm_per_h = state->lgar_bmi_input_params->PET_mm_per_h * mm_to_cm; @@ -550,7 +570,8 @@ Update() // move the wetting fronts without adding any water; this is done to close the mass balance // and also to merge / cross if necessary - temp_rch = lgar_move_wetting_fronts(subtimestep_h, &free_drainage_subtimestep_cm, &temp_pd, wf_free_drainage_demand, volend_subtimestep_cm, mass_correction_for_cached_free_drainage_fluxes, + temp_rch = lgar_move_wetting_fronts(subtimestep_h, &free_drainage_subtimestep_cm, &interflow_subtimestep_cm, + interflow_psi_threshold_cm, interflow_factor, &temp_pd, wf_free_drainage_demand, volend_subtimestep_cm, mass_correction_for_cached_free_drainage_fluxes, 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); @@ -639,7 +660,8 @@ Update() double volin_subtimestep_cm_temp = volin_subtimestep_cm; /* passing this for mass balance only, the method modifies it and returns percolated value, so we need to keep its original value stored to copy it back*/ - temp_rch = lgar_move_wetting_fronts(subtimestep_h, &free_drainage_subtimestep_cm, &volin_subtimestep_cm, wf_free_drainage_demand, volend_subtimestep_cm, mass_correction_for_cached_free_drainage_fluxes, + temp_rch = lgar_move_wetting_fronts(subtimestep_h, &free_drainage_subtimestep_cm, &interflow_subtimestep_cm, + interflow_psi_threshold_cm, interflow_factor, &volin_subtimestep_cm, wf_free_drainage_demand, volend_subtimestep_cm, mass_correction_for_cached_free_drainage_fluxes, 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); @@ -700,7 +722,7 @@ Update() volCRstart_subtimestep_cm = state->lgar_mass_balance.CR_fast_storage_cm + state->lgar_mass_balance.CR_slow_storage_cm; double volin_CR_subtimestep_cm = (precip_for_CR_subtimestep_cm_per_h + ponded_flux_for_CR)*subtimestep_h + free_drainage_for_CR; - double volQ_CR_subtimestep_cm = calc_CR_Q(subtimestep_h, a, a_slow, b, b_slow, frac_slow, precip_for_CR_subtimestep_cm_per_h + ponded_flux_for_CR + free_drainage_for_CR/subtimestep_h, &state->lgar_mass_balance.CR_fast_storage_cm, &state->lgar_mass_balance.CR_slow_storage_cm); + double volQ_CR_subtimestep_cm = calc_CR_Q(subtimestep_h, a_con_res, a_con_res_slow, b_con_res, b_con_res_slow, frac_slow, precip_for_CR_subtimestep_cm_per_h + ponded_flux_for_CR + free_drainage_for_CR/subtimestep_h, &state->lgar_mass_balance.CR_fast_storage_cm, &state->lgar_mass_balance.CR_slow_storage_cm); state->lgar_mass_balance.volrunoff_CR_cm += volQ_CR_subtimestep_cm; volQ_CR_timestep_cm += volQ_CR_subtimestep_cm; volCRend_subtimestep_cm = state->lgar_mass_balance.CR_fast_storage_cm + state->lgar_mass_balance.CR_slow_storage_cm; @@ -721,7 +743,7 @@ Update() // mass balance at the subtimestep (local mass balance) double local_mb = volstart_subtimestep_cm + precip_subtimestep_cm + volon_timestep_cm - volrunoff_subtimestep_cm - volQ_CR_subtimestep_cm - volCRend_subtimestep_cm + volCRstart_subtimestep_cm - - AET_subtimestep_cm - volon_subtimestep_cm - volrech_subtimestep_cm - volend_subtimestep_cm; + - AET_subtimestep_cm - volon_subtimestep_cm - volrech_subtimestep_cm - interflow_subtimestep_cm - volend_subtimestep_cm; /*----------------------------------------------------------------------*/ @@ -730,6 +752,7 @@ Update() //separating code such that most non substep vars (so xxx_timestep and not xxx_subtimestep) are updated in just one place. not all, because some must be set before substepping. volin_timestep_cm += volin_subtimestep_cm; volrech_timestep_cm += volrech_subtimestep_cm; + volinterflow_timestep_cm += interflow_subtimestep_cm; volrunoff_timestep_cm += volrunoff_subtimestep_cm; @@ -758,9 +781,10 @@ Update() Runoff = %14.10f \n\ AET = %14.10f \n\ Percolation = %14.10f \n\ + Interflow = %14.10f \n\ Final water = %14.10f \n", local_mb, volstart_subtimestep_cm, precip_subtimestep_cm, volon_subtimestep_cm, volin_subtimestep_cm, volrunoff_subtimestep_cm, AET_subtimestep_cm, volrech_subtimestep_cm, - volend_subtimestep_cm); + interflow_subtimestep_cm, volend_subtimestep_cm); } else { printf("\nLocal mass balance at this timestep... \n\ @@ -772,13 +796,14 @@ Update() Runoff (total) = %14.10f \n\ AET = %14.10f \n\ Percolation = %14.10f \n\ + Interflow = %14.10f \n\ Final water (LGAR) = %14.10f \n\ Water added (con res) = %14.10f \n\ Initial water (con res) = %14.10f \n\ Final water (con res) = %14.10f \n\ Runoff (con res) = %14.10f \n", local_mb, volstart_subtimestep_cm, precip_subtimestep_cm, volon_subtimestep_cm, volin_subtimestep_cm, volrunoff_subtimestep_cm, AET_subtimestep_cm, volrech_subtimestep_cm, - volend_subtimestep_cm, volin_CR_subtimestep_cm, volCRstart_subtimestep_cm, + interflow_subtimestep_cm, volend_subtimestep_cm, volin_CR_subtimestep_cm, volCRstart_subtimestep_cm, volCRend_subtimestep_cm, volQ_CR_subtimestep_cm); } @@ -805,7 +830,7 @@ Update() } // end of subcycling //update giuh at the time step level (was previously updated at the sub time step level) - volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm + volQ_CR_timestep_cm , num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue); + volrunoff_giuh_timestep_cm = giuh_convolution_integral(volrunoff_timestep_cm + volQ_CR_timestep_cm + volinterflow_timestep_cm, num_giuh_ordinates, giuh_ordinates, giuh_runoff_queue); // total mass of water leaving the system, at this time it is the giuh-only, but later will add groundwater component as well. // when groundwater component is added, it should probably happen inside of the subcycling loop. @@ -850,6 +875,7 @@ Update() state->lgar_mass_balance.volCRend_timestep_cm = volCRend_timestep_cm; state->lgar_mass_balance.volAET_timestep_cm = AET_timestep_cm; state->lgar_mass_balance.volrech_timestep_cm = volrech_timestep_cm; + state->lgar_mass_balance.volinterflow_timestep_cm = volinterflow_timestep_cm; state->lgar_mass_balance.volrunoff_timestep_cm = volrunoff_timestep_cm; state->lgar_mass_balance.volQ_timestep_cm = volQ_timestep_cm; state->lgar_mass_balance.volQ_CR_timestep_cm = volQ_CR_timestep_cm; @@ -871,6 +897,7 @@ Update() state->lgar_mass_balance.volCRend_cm = volCRend_timestep_cm; state->lgar_mass_balance.volAET_cm += AET_timestep_cm; state->lgar_mass_balance.volrech_cm += volrech_timestep_cm; + state->lgar_mass_balance.volinterflow_cm += volinterflow_timestep_cm; state->lgar_mass_balance.volrunoff_cm += volrunoff_timestep_cm; state->lgar_mass_balance.volQ_cm += volQ_timestep_cm; state->lgar_mass_balance.volQ_CR_cm += volQ_CR_timestep_cm; @@ -1012,15 +1039,17 @@ update_calibratable_parameters() 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 - <<", a = " << state->lgar_bmi_params.a - <<", b = " << state->lgar_bmi_params.b + <<", a_con_res = " << state->lgar_bmi_params.a_con_res + <<", b_con_res = " << state->lgar_bmi_params.b_con_res <<", frac_to_CR = " << state->lgar_bmi_params.frac_to_CR + <<", interflow_psi_threshold = " << state->lgar_bmi_params.interflow_psi_threshold_cm + <<", interflow_factor = " << state->lgar_bmi_params.interflow_factor <<", spf_factor = " << state->lgar_bmi_params.spf_factor << "\n"; if (state->lgar_bmi_params.frac_slow){ std::cerr - <<", a_slow = " << state->lgar_bmi_params.a_slow - <<", b_slow = " << state->lgar_bmi_params.b_slow + <<", a_con_res_slow = " << state->lgar_bmi_params.a_con_res_slow + <<", b_con_res_slow = " << state->lgar_bmi_params.b_con_res_slow <<", frac__slow = " << state->lgar_bmi_params.frac_slow << "\n"; } @@ -1030,21 +1059,29 @@ update_calibratable_parameters() 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; - state->lgar_bmi_params.a = state->lgar_calib_params.a; - state->lgar_bmi_params.b = state->lgar_calib_params.b; + state->lgar_bmi_params.a_con_res = state->lgar_calib_params.a_con_res; + state->lgar_bmi_params.b_con_res = state->lgar_calib_params.b_con_res; state->lgar_bmi_params.frac_to_CR = state->lgar_calib_params.frac_to_CR; + if (state->lgar_bmi_params.interflow_enabled){ + state->lgar_bmi_params.interflow_psi_threshold_cm = state->lgar_calib_params.interflow_psi_threshold_cm; + state->lgar_bmi_params.interflow_factor = state->lgar_calib_params.interflow_factor; + } state->lgar_bmi_params.spf_factor = state->lgar_calib_params.spf_factor; if (state->lgar_bmi_params.frac_slow){ - state->lgar_bmi_params.a_slow = state->lgar_calib_params.a_slow; - state->lgar_bmi_params.b_slow = state->lgar_calib_params.b_slow; + state->lgar_bmi_params.a_con_res_slow = state->lgar_calib_params.a_con_res_slow; + state->lgar_bmi_params.b_con_res_slow = state->lgar_calib_params.b_con_res_slow; state->lgar_bmi_params.frac_slow = state->lgar_calib_params.frac_slow; } if (state->lgar_bmi_params.log_mode){ - state->lgar_bmi_params.a = pow(10.0, state->lgar_calib_params.a); + state->lgar_bmi_params.a_con_res = pow(10.0, state->lgar_calib_params.a_con_res); if (state->lgar_bmi_params.frac_slow){ - state->lgar_bmi_params.a_slow = pow(10.0, state->lgar_calib_params.a_slow); + state->lgar_bmi_params.a_con_res_slow = pow(10.0, state->lgar_calib_params.a_con_res_slow); + } + if (state->lgar_bmi_params.interflow_enabled){ + state->lgar_bmi_params.interflow_psi_threshold_cm = pow(10.0, state->lgar_calib_params.interflow_psi_threshold_cm); + state->lgar_bmi_params.interflow_factor = pow(10.0, state->lgar_calib_params.interflow_factor); } } @@ -1052,15 +1089,17 @@ update_calibratable_parameters() 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 - <<", a = " << state->lgar_bmi_params.a - <<", b = " << state->lgar_bmi_params.b + <<", a_con_res = " << state->lgar_bmi_params.a_con_res + <<", b_con_res = " << state->lgar_bmi_params.b_con_res <<", frac_to_CR = " << state->lgar_bmi_params.frac_to_CR + <<", interflow_psi_threshold = " << state->lgar_bmi_params.interflow_psi_threshold_cm + <<", interflow_factor = " << state->lgar_bmi_params.interflow_factor <<", spf_factor = " << state->lgar_bmi_params.spf_factor << "\n"; if (state->lgar_bmi_params.frac_slow){ std::cerr - <<", a_slow = " << state->lgar_bmi_params.a_slow - <<", b_slow = " << state->lgar_bmi_params.b_slow + <<", a_con_res_slow = " << state->lgar_bmi_params.a_con_res_slow + <<", b_con_res_slow = " << state->lgar_bmi_params.b_con_res_slow <<", frac__slow = " << state->lgar_bmi_params.frac_slow << "\n"; } @@ -1120,8 +1159,9 @@ GetVarGrid(std::string name) || name.compare("potential_evapotranspiration") == 0 || name.compare("actual_evapotranspiration") == 0) // double return 1; - else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0 || name.compare("a") == 0 || name.compare("b") == 0 || name.compare("frac_to_CR") == 0 || name.compare("spf_factor") == 0 - || name.compare("a_slow") == 0 || name.compare("b_slow") == 0 || name.compare("frac_slow") == 0 || name.compare("soil_storage") == 0 || name.compare("field_capacity") == 0 || name.compare("ponded_depth_max") == 0)// double + else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0 || name.compare("a_con_res") == 0 || name.compare("a") == 0 || name.compare("b_con_res") == 0 || name.compare("b") == 0 || name.compare("frac_to_CR") == 0 || name.compare("spf_factor") == 0 + || name.compare("interflow_psi_threshold") == 0 || name.compare("lateral_flow_psi_threshold") == 0 || name.compare("interflow_factor") == 0 || name.compare("lateral_flow_factor") == 0 + || name.compare("a_con_res_slow") == 0 || name.compare("a_slow") == 0 || name.compare("b_con_res_slow") == 0 || name.compare("b_slow") == 0 || name.compare("frac_slow") == 0 || 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("conceptual_reservoir_to_stream_discharge") == 0) // double @@ -1229,8 +1269,9 @@ GetVarLocation(std::string name) name.compare("potential_evapotranspiration") == 0 || name.compare("potential_evapotranspiration_rate") == 0 || name.compare("actual_evapotranspiration") == 0) // double return "node"; - else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0 || name.compare("a") == 0 || name.compare("b") == 0 || name.compare("frac_to_CR") == 0 || name.compare("spf_factor") == 0 - || name.compare("a_slow") == 0 || name.compare("b_slow") == 0 || name.compare("frac_slow") == 0 || name.compare("soil_storage") == 0) // double + else if (name.compare("surface_runoff") == 0 || name.compare("giuh_runoff") == 0 || name.compare("a_con_res") == 0 || name.compare("a") == 0 || name.compare("b_con_res") == 0 || name.compare("b") == 0 || name.compare("frac_to_CR") == 0 || name.compare("spf_factor") == 0 + || name.compare("interflow_psi_threshold") == 0 || name.compare("lateral_flow_psi_threshold") == 0 || name.compare("interflow_factor") == 0 || name.compare("lateral_flow_factor") == 0 + || name.compare("a_con_res_slow") == 0 || name.compare("a_slow") == 0 || name.compare("b_con_res_slow") == 0 || name.compare("b_slow") == 0 || name.compare("frac_slow") == 0 || name.compare("soil_storage") == 0) // double return "node"; else if (name.compare("total_discharge") == 0 || name.compare("infiltration") == 0 || name.compare("percolation") == 0 || name.compare("conceptual_reservoir_to_stream_discharge") == 0) // double @@ -1404,18 +1445,22 @@ GetValuePtr (std::string name) 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 if (name.compare("a") == 0) - return (void*)&this->state->lgar_calib_params.a; - else if (name.compare("b") == 0) - return (void*)&this->state->lgar_calib_params.b; + else if (name.compare("a_con_res") == 0 || name.compare("a") == 0) + return (void*)&this->state->lgar_calib_params.a_con_res; + else if (name.compare("b_con_res") == 0 || name.compare("b") == 0) + return (void*)&this->state->lgar_calib_params.b_con_res; else if (name.compare("frac_to_CR") == 0) return (void*)&this->state->lgar_calib_params.frac_to_CR; - else if (name.compare("a_slow") == 0) - return (void*)&this->state->lgar_calib_params.a_slow; - else if (name.compare("b_slow") == 0) - return (void*)&this->state->lgar_calib_params.b_slow; + else if (name.compare("a_con_res_slow") == 0 || name.compare("a_slow") == 0) + return (void*)&this->state->lgar_calib_params.a_con_res_slow; + else if (name.compare("b_con_res_slow") == 0 || name.compare("b_slow") == 0) + return (void*)&this->state->lgar_calib_params.b_con_res_slow; else if (name.compare("frac_slow") == 0) return (void*)&this->state->lgar_calib_params.frac_slow; + else if (name.compare("interflow_psi_threshold") == 0 || name.compare("lateral_flow_psi_threshold") == 0) + return (void*)&this->state->lgar_calib_params.interflow_psi_threshold_cm; + else if (name.compare("interflow_factor") == 0 || name.compare("lateral_flow_factor") == 0) + return (void*)&this->state->lgar_calib_params.interflow_factor; else if (name.compare("spf_factor") == 0) return (void*)&this->state->lgar_calib_params.spf_factor; else { diff --git a/src/conceptual_reservoir.cxx b/src/conceptual_reservoir.cxx index f2fd07f..6d37389 100644 --- a/src/conceptual_reservoir.cxx +++ b/src/conceptual_reservoir.cxx @@ -6,8 +6,8 @@ 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, @@ -18,7 +18,7 @@ extern double calc_CR_Q( double input_fast = precip_for_CR_subtimestep_cm_per_h - input_slow; // implicit (1 - frac_slow) // === FAST reservoir outflow === - double Q_fast = subtimestep_h * (a_fast * pow(*CR_fast_storage_cm, b_fast)); + double Q_fast = subtimestep_h * (a_con_res * pow(*CR_fast_storage_cm, b_con_res)); if (*CR_fast_storage_cm < 0.01) Q_fast = 0.0; double delta_fast = subtimestep_h * input_fast - Q_fast; @@ -30,7 +30,7 @@ extern double calc_CR_Q( } // === SLOW reservoir outflow === - double Q_slow = subtimestep_h * (a_slow * pow(*CR_slow_storage_cm, b_slow)); + double Q_slow = subtimestep_h * (a_con_res_slow * pow(*CR_slow_storage_cm, b_con_res_slow)); if (*CR_slow_storage_cm < 0.01) Q_slow = 0.0; double delta_slow = subtimestep_h * input_slow - Q_slow; diff --git a/src/lgar.cxx b/src/lgar.cxx index 0e8a966..5c0fcd0 100755 --- a/src/lgar.cxx +++ b/src/lgar.cxx @@ -104,12 +104,14 @@ extern void lgar_initialize(string config_file, struct model_state *state) // also initialize calibratable parameters state->lgar_calib_params.field_capacity_psi = state->lgar_bmi_params.field_capacity_psi_cm; state->lgar_calib_params.ponded_depth_max = state->lgar_bmi_params.ponded_depth_max_cm; - state->lgar_calib_params.a = state->lgar_bmi_params.a; - state->lgar_calib_params.b = state->lgar_bmi_params.b; + state->lgar_calib_params.a_con_res = state->lgar_bmi_params.a_con_res; + state->lgar_calib_params.b_con_res = state->lgar_bmi_params.b_con_res; state->lgar_calib_params.frac_to_CR = state->lgar_bmi_params.frac_to_CR; - state->lgar_calib_params.a_slow = state->lgar_bmi_params.a_slow; - state->lgar_calib_params.b_slow = state->lgar_bmi_params.b_slow; + state->lgar_calib_params.a_con_res_slow = state->lgar_bmi_params.a_con_res_slow; + state->lgar_calib_params.b_con_res_slow = state->lgar_bmi_params.b_con_res_slow; state->lgar_calib_params.frac_slow = state->lgar_bmi_params.frac_slow; + state->lgar_calib_params.interflow_psi_threshold_cm = state->lgar_bmi_params.interflow_psi_threshold_cm; + state->lgar_calib_params.interflow_factor = state->lgar_bmi_params.interflow_factor; state->lgar_calib_params.spf_factor = state->lgar_bmi_params.spf_factor; struct wetting_front *current = state->head; @@ -181,7 +183,9 @@ extern void lgar_initialize(string config_file, struct model_state *state) state->lgar_mass_balance.volCRend_cm = 0.0; state->lgar_mass_balance.volAET_cm = 0.0; state->lgar_mass_balance.volrech_cm = 0.0; + state->lgar_mass_balance.volinterflow_timestep_cm = 0.0; state->lgar_mass_balance.volrunoff_cm = 0.0; + state->lgar_mass_balance.volinterflow_cm = 0.0; state->lgar_mass_balance.volrunoff_giuh_cm = 0.0; state->lgar_mass_balance.volQ_cm = 0.0; state->lgar_mass_balance.volQ_CR_cm = 0.0; @@ -278,6 +282,9 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) state->lgar_bmi_params.log_mode = false; state->lgar_bmi_params.free_drainage_enabled = false; state->lgar_bmi_params.free_drainage_to_CR = false; + state->lgar_bmi_params.interflow_enabled = false; + state->lgar_bmi_params.interflow_psi_threshold_cm = 0.0; + state->lgar_bmi_params.interflow_factor = 0.0; // setting mass balance tolerance to be large by default; this can be specified in the config file state->lgar_bmi_params.mbal_tol = 1.E1; @@ -289,12 +296,14 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) 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_a_set = false; - bool is_b_set = false; + bool is_a_con_res_set = false; + bool is_b_con_res_set = false; bool is_frac_to_CR_set = false; - bool is_a_slow_set = false; - bool is_b_slow_set = false; + bool is_a_con_res_slow_set = false; + bool is_b_con_res_slow_set = false; bool is_frac_slow_set = false; + bool is_interflow_psi_threshold_set = false; + bool is_interflow_factor_set = false; bool is_soil_params_file_set = false; bool is_max_valid_soil_types_set = false; bool is_giuh_ordinates_set = false; @@ -457,22 +466,22 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) continue; } - else if (param_key == "a") { - state->lgar_bmi_params.a = stod(param_value); - is_a_set = true; + else if (param_key == "a_con_res" || param_key == "a") { + state->lgar_bmi_params.a_con_res = stod(param_value); + is_a_con_res_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"a : "<lgar_bmi_params.a<<"\n"; + std::cerr<<"a_con_res"<<(param_key == "a" ? " (using old name in config)" : "")<<" : "<lgar_bmi_params.a_con_res<<"\n"; std::cerr<<" ***** \n"; } continue; } - else if (param_key == "b") { - state->lgar_bmi_params.b = stod(param_value); - is_b_set = true; + else if (param_key == "b_con_res" || param_key == "b") { + state->lgar_bmi_params.b_con_res = stod(param_value); + is_b_con_res_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"b : "<lgar_bmi_params.b<<"\n"; + std::cerr<<"b_con_res"<<(param_key == "b" ? " (using old name in config)" : "")<<" : "<lgar_bmi_params.b_con_res<<"\n"; std::cerr<<" ***** \n"; } continue; @@ -487,22 +496,22 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } continue; } - else if (param_key == "a_slow") { - state->lgar_bmi_params.a_slow = stod(param_value); - is_a_slow_set = true; + else if (param_key == "a_con_res_slow" || param_key == "a_slow") { + state->lgar_bmi_params.a_con_res_slow = stod(param_value); + is_a_con_res_slow_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"a_slow : "<lgar_bmi_params.a_slow<<"\n"; + std::cerr<<"a_con_res_slow"<<(param_key == "a_slow" ? " (using old name in config)" : "")<<" : "<lgar_bmi_params.a_con_res_slow<<"\n"; std::cerr<<" ***** \n"; } continue; } - else if (param_key == "b_slow") { - state->lgar_bmi_params.b_slow = stod(param_value); - is_b_slow_set = true; + else if (param_key == "b_con_res_slow" || param_key == "b_slow") { + state->lgar_bmi_params.b_con_res_slow = stod(param_value); + is_b_con_res_slow_set = true; if (verbosity.compare("high") == 0) { - std::cerr<<"b_slow : "<lgar_bmi_params.b_slow<<"\n"; + std::cerr<<"b_con_res_slow"<<(param_key == "b_slow" ? " (using old name in config)" : "")<<" : "<lgar_bmi_params.b_con_res_slow<<"\n"; std::cerr<<" ***** \n"; } continue; @@ -517,6 +526,26 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) } continue; } + else if (param_key == "interflow_psi_threshold" || param_key == "lateral_flow_psi_threshold") { + state->lgar_bmi_params.interflow_psi_threshold_cm = stod(param_value); + is_interflow_psi_threshold_set = true; + + if (verbosity.compare("high") == 0) { + std::cerr<<"interflow_psi_threshold"<<(param_key.rfind("lateral_flow", 0) == 0 ? " (using old name in config)" : "")<<" [cm] : "<lgar_bmi_params.interflow_psi_threshold_cm<<"\n"; + std::cerr<<" ***** \n"; + } + continue; + } + else if (param_key == "interflow_factor" || param_key == "lateral_flow_factor") { + state->lgar_bmi_params.interflow_factor = stod(param_value); + is_interflow_factor_set = true; + + if (verbosity.compare("high") == 0) { + std::cerr<<"interflow_factor"<<(param_key == "lateral_flow_factor" ? " (using old name in config)" : "")<<" : "<lgar_bmi_params.interflow_factor<<"\n"; + std::cerr<<" ***** \n"; + } + continue; + } else if (param_key == "frac_to_GW") { state->lgar_bmi_params.frac_to_CR = stod(param_value); is_frac_to_CR_set = true; @@ -613,7 +642,7 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) else if ( (param_value == "true") || (param_value == "1")) { state->lgar_bmi_params.log_mode = true; if (verbosity.compare("high") == 0) { - printf("log_mode enabled. So K_s for each layer, alpha for each layer, and a for the nonlinear reservoir(s) will use the log of their input values. \n"); + printf("log_mode enabled. So K_s for each layer, alpha for each layer, a_con_res for the nonlinear reservoir(s), interflow_psi_threshold, and interflow_factor will use the log of their input values. \n"); } } else { @@ -793,11 +822,37 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) throw runtime_error(errMsg.str()); } + // if log_mode enabled, K_s for each layer, alpha for each layer, a_con_res for the nonlinear reservoir(s), interflow_psi_threshold, and interflow_factor will use the log of their input values. + // so for example if a value of 0.01 cm/h is desired for a K_s value, you should use -2.0 for that K_s, because 10^-2 = 0.01. if (state->lgar_bmi_params.log_mode){ - state->lgar_bmi_params.a = pow(10.0, state->lgar_bmi_params.a); - if (is_a_slow_set){ - state->lgar_bmi_params.a_slow = pow(10.0, state->lgar_bmi_params.a_slow); + state->lgar_bmi_params.a_con_res = pow(10.0, state->lgar_bmi_params.a_con_res); + if (is_a_con_res_slow_set){ + state->lgar_bmi_params.a_con_res_slow = pow(10.0, state->lgar_bmi_params.a_con_res_slow); + } + if (is_interflow_psi_threshold_set && is_interflow_factor_set){ + state->lgar_bmi_params.interflow_psi_threshold_cm = pow(10.0, state->lgar_bmi_params.interflow_psi_threshold_cm); + state->lgar_bmi_params.interflow_factor = pow(10.0, state->lgar_bmi_params.interflow_factor); + } + } + + if (is_interflow_psi_threshold_set != is_interflow_factor_set) { + stringstream errMsg; + errMsg << "The configuration file \'" << config_file <<"\' does not correctly set interflow_psi_threshold and interflow_factor. Either both or neither must be set. In log_mode, both interflow parameters are interpreted as log10 values. Legacy names lateral_flow_psi_threshold and lateral_flow_factor are still accepted. \n"; + throw runtime_error(errMsg.str()); + } + + if (is_interflow_factor_set) { + if (state->lgar_bmi_params.interflow_psi_threshold_cm < 0.0) { + stringstream errMsg; + errMsg << "The configuration file \'" << config_file <<"\' sets interflow_psi_threshold below zero. The threshold must be >= 0 cm. \n"; + throw runtime_error(errMsg.str()); } + if (state->lgar_bmi_params.interflow_factor < 0.0) { + stringstream errMsg; + errMsg << "The configuration file \'" << config_file <<"\' sets interflow_factor below 0. \n"; + throw runtime_error(errMsg.str()); + } + state->lgar_bmi_params.interflow_enabled = true; } if(is_soil_params_file_set) { @@ -886,20 +941,20 @@ extern void InitFromConfigFile(string config_file, struct model_state *state) throw runtime_error(errMsg.str()); } - if (! ( (is_a_set == is_b_set) && (is_frac_to_CR_set == is_b_set)) ){ + if (! ( (is_a_con_res_set == is_b_con_res_set) && (is_frac_to_CR_set == is_b_con_res_set)) ){ //in this case, it must be either the case that all of these have been set (the user wants a nonlinear reservoir), or that none of these are set (the user does not want this). //it can not be the case that only one or two of these three have been set. stringstream errMsg; - errMsg << "The configuration file \'" << config_file <<"\' does not correctly set a, b, and frac_to_CR. Either all or none must be set. a and b must be 0 or greater and frac_to_CR must be between 0 and 1. \n"; + errMsg << "The configuration file \'" << config_file <<"\' does not correctly set a_con_res, b_con_res, and frac_to_CR. Either all or none must be set. a_con_res and b_con_res must be 0 or greater and frac_to_CR must be between 0 and 1. Legacy names a and b are still accepted. \n"; throw runtime_error(errMsg.str()); } - if (! ( (is_a_slow_set == is_b_slow_set) && (is_frac_slow_set == is_b_slow_set)) ){ + if (! ( (is_a_con_res_slow_set == is_b_con_res_slow_set) && (is_frac_slow_set == is_b_con_res_slow_set)) ){ //in this case, it must be either the case that all of these have been set (the user wants a second nonlinear reservoir), or that none of these are set (the user does not want this). //technically you can set the "slow" reservoir and not the other one -- in either case it amounts to 1 nonlinear reservoir. //it can not be the case that only one or two of these three have been set. stringstream errMsg; - errMsg << "The configuration file \'" << config_file <<"\' does not correctly set a_slow, b_slow, and frac_slow. Either all or none must be set. a_slow and b_slow must be 0 or greater and frac_slow must be between 0 and 1 (but greater than 0). \n"; + errMsg << "The configuration file \'" << config_file <<"\' does not correctly set a_con_res_slow, b_con_res_slow, and frac_slow. Either all or none must be set. a_con_res_slow and b_con_res_slow must be 0 or greater and frac_slow must be between 0 and 1 (but greater than 0). Legacy names a_slow and b_slow are still accepted. \n"; throw runtime_error(errMsg.str()); } @@ -1165,6 +1220,7 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run double volprecip = state->lgar_mass_balance.volprecip_cm; double volrunoff = state->lgar_mass_balance.volrunoff_cm; double volrunoff_CR = state->lgar_mass_balance.volrunoff_CR_cm; + double volinterflow = state->lgar_mass_balance.volinterflow_cm; double volAET = state->lgar_mass_balance.volAET_cm; double volPET = state->lgar_mass_balance.volPET_cm; double volon = state->lgar_mass_balance.volon_cm; @@ -1182,7 +1238,7 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run for(int i=0; i <= state->lgar_bmi_params.num_giuh_ordinates; i++) volend_giuh_cm += giuh_runoff_queue_cm[i]; - double global_error_cm = volstart + volprecip - volrunoff - volAET - volon - volrech - volend + volchange_calib_cm - volrunoff_CR - volCRend; + double global_error_cm = volstart + volprecip - volrunoff - volAET - volon - volrech - volinterflow - volend + volchange_calib_cm - volrunoff_CR - volCRend; printf("\n********************************************************* \n"); printf("-------------------- Simulation Summary ----------------- \n"); @@ -1197,6 +1253,7 @@ extern void lgar_global_mass_balance(struct model_state *state, double *giuh_run printf("GIUH runoff = %14.10f cm\n", volrunoff_giuh); printf("GIUH water (in array) = %14.10f cm\n", volend_giuh_cm); printf("Total percolation = %14.10f cm\n", volrech); + printf("Total interflow = %14.10f cm\n", volinterflow); printf("Total AET = %14.10f cm\n", volAET); printf("Total PET = %14.10f cm\n", volPET); if (state->lgar_bmi_params.frac_to_CR){ @@ -1253,6 +1310,366 @@ extern int wetting_front_free_drainage(struct wetting_front* head) { return wf_that_supplies_free_drainage_demand; } +// ############################################################################################ +/* + Compute the vertical support depth represented by one wetting front within its soil layer. + + The support depth is the portion of the layer assigned to current for interflow scaling. + If another wetting front exists immediately above current in the same layer, the support + depth starts at that upper front's depth. Otherwise it starts at the top of the layer. For a + to_bottom front, the support depth extends to the bottom of the layer; for a normal front, it + extends to current->depth_cm. + + The returned depth is later divided by total column depth to scale candidate interflow by + the fraction of the LGAR domain represented by this wetting front. +*/ +// ############################################################################################ +double lgar_interflow_support_depth_cm(double layer_top_cm, double layer_bottom_cm, + const struct wetting_front *previous, + const struct wetting_front *current) +{ + if (current == NULL) + return 0.0; + + double support_top_cm = layer_top_cm; + if (previous != NULL && previous->layer_num == current->layer_num) + support_top_cm = fmax(layer_top_cm, previous->depth_cm); + + double support_bottom_cm = current->to_bottom ? layer_bottom_cm : current->depth_cm; + support_bottom_cm = fmin(support_bottom_cm, layer_bottom_cm); + + return fmax(0.0, support_bottom_cm - support_top_cm); +} + + +/* + Compute candidate interflow for each wetting front for the current subtimestep. + + A wetting front is eligible when its capillary head is less than or equal to + interflow_psi_threshold_cm. The candidate interflow is K(theta) times + interflow_factor, scaled by the fraction of the LGAR column represented by + that wetting front. The function only fills interflow_flux_cm_by_front; it does + not remove water from the wetting front state. + + For non-deepest to_bottom fronts, the interflow amount is assigned to the + first non-to_bottom wetting front below it. If no such front exists, the amount + is skipped here; the deepest to_bottom front can still be assigned interflow + directly and is handled later by the to_bottom stack solver. + + Also note that in the interflow code, "stack" refers to a to_bottom wetting + front and any consecutive to_bottom wetting fronts directly above it. +*/ +static void lgar_calc_interflow_fluxes_by_front(double timestep_h, int num_layers, double interflow_psi_threshold_cm, + double interflow_factor, double *cum_layer_thickness_cm, + struct wetting_front* head, std::vector& interflow_flux_cm_by_front) +{ + if (interflow_factor <= 0.0 || head == NULL) + return; + + double column_depth_cm = cum_layer_thickness_cm[num_layers]; + if (column_depth_cm <= 0.0) + return; + + for (struct wetting_front *previous = NULL, *current = head; current != NULL; previous = current, current = current->next) { + double interflow_flux_cm_per_h = 0.0; + + if (current->psi_cm <= interflow_psi_threshold_cm) { + double layer_top_cm = cum_layer_thickness_cm[current->layer_num - 1]; + double layer_bottom_cm = cum_layer_thickness_cm[current->layer_num]; + double support_depth_cm = lgar_interflow_support_depth_cm(layer_top_cm, layer_bottom_cm, previous, current); + + support_depth_cm = fmax(0.0, fmin(support_depth_cm, column_depth_cm)); + double vadose_fraction = support_depth_cm / column_depth_cm; + + interflow_flux_cm_per_h = fmax(0.0, current->K_cm_per_h) * interflow_factor * vadose_fraction; + } + + if (interflow_flux_cm_per_h <= 0.0) + continue; + + struct wetting_front *target = current; + + if (current->to_bottom && current->layer_num < num_layers) { + target = current->next; + while (target != NULL && target->to_bottom) + target = target->next; + } + + if (target == NULL) + continue; + + if (target->front_num >= 0 && target->front_num < (int)interflow_flux_cm_by_front.size()) { + interflow_flux_cm_by_front[target->front_num] += interflow_flux_cm_per_h * timestep_h; + if (verbosity.compare("high") == 0) { + printf("Interflow assigned from WF %d to WF %d: rate = %.10e cm/h, amount = %.10e cm\n", + current->front_num, target->front_num, interflow_flux_cm_per_h, interflow_flux_cm_per_h * timestep_h); + } + } + } +} + + +/* + Compute the mass expression used by lgar_theta_mass_balance at a specified + capillary head. + + The returned value is the amount of water represented by the current wetting + front mass balance term, in cm over the model column. This lets interflow + cap removals at the driest theta the model can represent before the iterative + mass balance solve is called. +*/ +static double lgar_prior_mass_for_psi(double psi_cm, int layer_num, double *delta_theta, + double *delta_thickness, int *soil_type, + struct soil_properties_ *soil_properties) +{ + double prior_mass = 0.0; + + for (int k = 1; k <= layer_num; k++) { + int soil_num = soil_type[k]; + double theta = calc_theta_from_h(psi_cm, soil_properties[soil_num].vg_alpha_per_cm, + soil_properties[soil_num].vg_m, soil_properties[soil_num].vg_n, + soil_properties[soil_num].theta_e, soil_properties[soil_num].theta_r); + prior_mass += delta_thickness[k] * (theta - delta_theta[k]); + } + + return prior_mass; +} + + +/* + Return the driest capillary head that interflow is allowed to create. + + Interflow is activated at interflow_psi_threshold_cm, and a single + interflow removal is capped so the updated wetting front does not dry more + than 1 cm past that threshold in terms of capillary head. The tiny buffer + keeps numerical solves inside that limit. +*/ +static double lgar_interflow_psi_cap_cm(double interflow_psi_threshold_cm) +{ + const double interflow_psi_cap_buffer_cm = 1.0e-3; + return fmax(0.0, fmin(interflow_psi_threshold_cm + 1.0 - interflow_psi_cap_buffer_cm, + PSI_UPPER_LIM)); +} + + +/* + Apply interflow to a wetting front mass balance term. + + The requested interflow is capped by the removable mass above + minimum_prior_mass_cm. That minimum is computed by the caller from the same + mass expression that will be used to update theta, so interflow cannot ask + the following mass balance solve to dry the wetting front beyond the chosen + interflow psi cap. The applied amount is subtracted from prior_mass, added + to the cumulative interflow for the subtimestep, and returned. +*/ +static double lgar_apply_interflow_flux_to_prior_mass(double requested_interflow_flux_cm, double *prior_mass, + double minimum_prior_mass_cm, + double *interflow_subtimestep_cm) +{ + if (requested_interflow_flux_cm <= 0.0 || prior_mass == NULL) + return 0.0; + + double applied_interflow_flux_cm = fmin(requested_interflow_flux_cm, + fmax(*prior_mass - fmax(minimum_prior_mass_cm, 0.0), 0.0)); + *prior_mass -= applied_interflow_flux_cm; + + if (interflow_subtimestep_cm != NULL) + *interflow_subtimestep_cm += applied_interflow_flux_cm; + + return applied_interflow_flux_cm; +} + + +/* + Compute the one-layer mass expression used by the to_bottom stack solver. + + The portion below the nearest same-layer non-to_bottom front is controlled by + stack_theta. Any portion above that front is included only as a fixed offset + and is not changed by the stack solve. This helper does not reconstruct an + arbitrary multi-front profile above the to_bottom region. +*/ +double lgar_to_bottom_stack_layer_mass_cm(int layer_num, double layer_top_cm, double layer_bottom_cm, + double stack_theta, const struct wetting_front *previous_front) +{ + if (previous_front != NULL && !previous_front->to_bottom && previous_front->layer_num == layer_num) { + double upper_depth_cm = fmax(layer_top_cm, fmin(previous_front->depth_cm, layer_bottom_cm)); + return (upper_depth_cm - layer_top_cm) * previous_front->theta + (layer_bottom_cm - upper_depth_cm) * stack_theta; + } + + return (layer_bottom_cm - layer_top_cm) * stack_theta; +} + + +/* + Compute the current water mass represented by a consecutive stack of to_bottom wetting fronts. + + This helper is used when interflow is applied to the deepest to_bottom wetting front stack. Each + front in the stack represents the water content of its soil layer down to that layer's bottom. If + there is a non-to_bottom wetting front immediately above a stack front in the same layer, only the + portion below that upper front is counted. This prevents the stack mass from including water that is + represented by a separate wetting front above the to_bottom stack. + + The returned value is an equivalent water depth in cm over the model column. +*/ +static double lgar_to_bottom_stack_mass_from_profile(int stack_start_front_num, int stack_end_front_num, + double *cum_layer_thickness_cm, struct wetting_front* head) +{ + double stack_mass_cm = 0.0; + + for (int front_num = stack_start_front_num; front_num <= stack_end_front_num; front_num++) { + struct wetting_front *front = listFindFront(front_num, head, NULL); + if (front == NULL) + return 0.0; + + int layer_num = front->layer_num; + double layer_top_cm = cum_layer_thickness_cm[layer_num - 1]; + double layer_bottom_cm = cum_layer_thickness_cm[layer_num]; + struct wetting_front *previous_front = front_num > 1 ? listFindFront(front_num - 1, head, NULL) : NULL; + stack_mass_cm += lgar_to_bottom_stack_layer_mass_cm(layer_num, layer_top_cm, layer_bottom_cm, + front->theta, previous_front); + } + + return stack_mass_cm; +} + +/* + Compute the water mass that a consecutive to_bottom wetting front stack would + have if all fronts in the stack shared the specified capillary head. + + This is used by the stack interflow solver to search for the new common psi + that gives the target stack mass after interflow has been removed. +*/ +static double lgar_to_bottom_stack_mass_for_psi(double psi_cm, int stack_start_front_num, int stack_end_front_num, + double *cum_layer_thickness_cm, int *soil_type, + struct wetting_front* head, struct soil_properties_ *soil_properties) +{ + double stack_mass_cm = 0.0; + + for (int front_num = stack_start_front_num; front_num <= stack_end_front_num; front_num++) { + struct wetting_front *front = listFindFront(front_num, head, NULL); + if (front == NULL) + continue; + + int layer_num = front->layer_num; + int soil_num = soil_type[layer_num]; + double layer_top_cm = cum_layer_thickness_cm[layer_num - 1]; + double layer_bottom_cm = cum_layer_thickness_cm[layer_num]; + struct wetting_front *previous_front = front_num > 1 ? listFindFront(front_num - 1, head, NULL) : NULL; + double theta = calc_theta_from_h(psi_cm, soil_properties[soil_num].vg_alpha_per_cm, + soil_properties[soil_num].vg_m, soil_properties[soil_num].vg_n, + soil_properties[soil_num].theta_e, soil_properties[soil_num].theta_r); + stack_mass_cm += lgar_to_bottom_stack_layer_mass_cm(layer_num, layer_top_cm, layer_bottom_cm, + theta, previous_front); + } + + return stack_mass_cm; +} + +/* + Return whether a wetting front was updated by the deepest to_bottom stack + interflow solve during this subtimestep. + + Upper wetting front mass balance calculations use this flag to decide whether + the lower reference front should come from the previous state or from the + already-updated current stack state. +*/ +static bool lgar_front_changed_by_interflow_stack(const struct wetting_front *front, + const std::vector& interflow_stack_changed_by_front) +{ + return front != NULL && front->front_num >= 0 + && front->front_num < (int)interflow_stack_changed_by_front.size() + && interflow_stack_changed_by_front[front->front_num] != 0; +} + +/* + Apply interflow to the deepest to_bottom wetting front stack. + + A deepest to_bottom front can share one capillary head with a consecutive stack + of to_bottom fronts above it. In that case, removing interflow from only one + front would not conserve the stack mass consistently. This function finds the + consecutive to_bottom stack, caps the requested interflow by the removable + stack mass, solves for the new common psi that gives the reduced stack mass + without drying past the interflow psi cap, updates theta, psi, and K for + every front in the stack, marks those fronts as changed, and returns the + applied interflow. +*/ +static double lgar_apply_interflow_flux_to_deepest_to_bottom_stack(double requested_interflow_flux_cm, double interflow_psi_cap_cm, + int stack_end_front_num, + int num_layers, double *cum_layer_thickness_cm, + int *soil_type, double *frozen_factor, + struct wetting_front** head, + struct wetting_front* state_previous, + struct soil_properties_ *soil_properties, + double *interflow_subtimestep_cm, + std::vector *interflow_stack_changed_by_front) +{ + if (requested_interflow_flux_cm <= 0.0 || head == NULL || *head == NULL) + return 0.0; + + struct wetting_front *stack_end = listFindFront(stack_end_front_num, *head, NULL); + if (stack_end == NULL || !stack_end->to_bottom || stack_end->layer_num != num_layers) + return 0.0; + + int stack_start_front_num = stack_end_front_num; + while (stack_start_front_num > 1) { + struct wetting_front *previous_front = listFindFront(stack_start_front_num - 1, *head, NULL); + if (previous_front == NULL || !previous_front->to_bottom) + break; + stack_start_front_num--; + } + + double prior_stack_mass_cm = lgar_to_bottom_stack_mass_from_profile(stack_start_front_num, stack_end_front_num, + cum_layer_thickness_cm, state_previous); + + double minimum_stack_mass_cm = lgar_to_bottom_stack_mass_for_psi(interflow_psi_cap_cm, stack_start_front_num, stack_end_front_num, + cum_layer_thickness_cm, soil_type, *head, soil_properties); + double applied_interflow_flux_cm = fmin(requested_interflow_flux_cm, fmax(prior_stack_mass_cm - minimum_stack_mass_cm, 0.0)); + if (applied_interflow_flux_cm <= 0.0) + return 0.0; + + double target_stack_mass_cm = prior_stack_mass_cm - applied_interflow_flux_cm; + double psi_low_cm = 0.0; + double psi_high_cm = interflow_psi_cap_cm; + + // Bisection search for the common psi that produces the target post-interflow stack mass. + for (int iter = 0; iter < 120; iter++) { + double psi_mid_cm = 0.5 * (psi_low_cm + psi_high_cm); + double stack_mass_cm = lgar_to_bottom_stack_mass_for_psi(psi_mid_cm, stack_start_front_num, stack_end_front_num, + cum_layer_thickness_cm, soil_type, *head, soil_properties); + if (stack_mass_cm > target_stack_mass_cm) + psi_low_cm = psi_mid_cm; + else + psi_high_cm = psi_mid_cm; + } + + double psi_new_cm = 0.5 * (psi_low_cm + psi_high_cm); + for (int front_num = stack_start_front_num; front_num <= stack_end_front_num; front_num++) { + struct wetting_front *front = listFindFront(front_num, *head, NULL); + if (front == NULL) + continue; + + int layer_num = front->layer_num; + int soil_num = soil_type[layer_num]; + front->theta = calc_theta_from_h(psi_new_cm, soil_properties[soil_num].vg_alpha_per_cm, + soil_properties[soil_num].vg_m, soil_properties[soil_num].vg_n, + soil_properties[soil_num].theta_e, soil_properties[soil_num].theta_r); + front->psi_cm = psi_new_cm; + + double Se = calc_Se_from_theta(front->theta, soil_properties[soil_num].theta_e, soil_properties[soil_num].theta_r); + front->K_cm_per_h = calc_K_from_Se(Se, frozen_factor[layer_num] * soil_properties[soil_num].Ksat_cm_per_h, + soil_properties[soil_num].vg_m); + + if (interflow_stack_changed_by_front != NULL && front->front_num >= 0 + && front->front_num < (int)interflow_stack_changed_by_front->size()) + (*interflow_stack_changed_by_front)[front->front_num] = 1; + } + + if (interflow_subtimestep_cm != NULL) + *interflow_subtimestep_cm += applied_interflow_flux_cm; + + return applied_interflow_flux_cm; +} + // ####################################################################################################### /* the function moves wetting fronts, merge wetting fronts and does the mass balance correction when needed @@ -1266,7 +1683,8 @@ extern int wetting_front_free_drainage(struct wetting_front* head) { Note: '_old' denotes the wetting_front or variables at the previous timestep (or state) */ // ####################################################################################################### -extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_subtimestep_cm, double *volin_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 *volin_cm, int wf_free_drainage_demand, double old_mass, double mass_correction_for_cached_free_drainage_fluxes, int num_layers, double *AET_demand_cm, double *cum_layer_thickness_cm, int *soil_type, double *frozen_factor, struct wetting_front** head, struct wetting_front* state_previous, struct soil_properties_ *soil_properties) @@ -1292,6 +1710,12 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ int layer_num, soil_num; int number_of_wetting_fronts = listLength(*head); + std::vector interflow_flux_cm_by_front(number_of_wetting_fronts + 1, 0.0); + std::vector interflow_stack_changed_by_front(number_of_wetting_fronts + 1, 0); + lgar_calc_interflow_fluxes_by_front(timestep_h, num_layers, interflow_psi_threshold_cm, + interflow_factor, cum_layer_thickness_cm, + state_previous, interflow_flux_cm_by_front); + double interflow_psi_cap_cm = lgar_interflow_psi_cap_cm(interflow_psi_threshold_cm); current = *head; @@ -1460,6 +1884,15 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ if (wf_free_drainage_demand == wf) prior_mass += precip_mass_to_add - (free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes + actual_ET_demand); + double minimum_prior_mass_cm = lgar_prior_mass_for_psi(interflow_psi_cap_cm, layer_num, delta_thetas, + delta_thickness, soil_type, soil_properties); + double applied_interflow_flux_cm = lgar_apply_interflow_flux_to_prior_mass(interflow_flux_cm_by_front[wf], &prior_mass, + minimum_prior_mass_cm, + interflow_subtimestep_cm); + if (applied_interflow_flux_cm > 0.0 && verbosity.compare("high") == 0) { + printf("Applied interflow to WF %d mass balance: %.10e cm\n", wf, applied_interflow_flux_cm); + } + // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, precip_mass_to_add, AET_demand_cm, @@ -1480,6 +1913,26 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ } + // case to apply interflow to the deepest to_bottom wetting front when other wetting fronts + // exist above it. This deepest front has no next wetting front, so it is not handled by the + // within-layer mass balance cases below, but it may still contribute interflow. Because + // interface to_bottom fronts above it share its psi, solve the whole connected to_bottom stack + // so the storage change equals the interflow flux counted in the mass balance. + /*************************************************************************************/ + if (wf == number_of_wetting_fronts && current->to_bottom && current->layer_num == num_layers && number_of_wetting_fronts > num_layers) { + double applied_interflow_flux_cm = lgar_apply_interflow_flux_to_deepest_to_bottom_stack(interflow_flux_cm_by_front[wf], + interflow_psi_cap_cm, wf, + num_layers, cum_layer_thickness_cm, + soil_type, frozen_factor, head, + state_previous, soil_properties, + interflow_subtimestep_cm, + &interflow_stack_changed_by_front); + if (applied_interflow_flux_cm > 0.0 && verbosity.compare("high") == 0) { + printf("Applied interflow to deepest to_bottom WF %d stack mass balance: %.10e cm\n", wf, applied_interflow_flux_cm); + } + } + + // case to check if the 'current' wetting front is within the layer and not at the layer's interface // layer_num == layer_num_below means there is another wetting front below the current wetting front // and they both belong to the same layer (in simple words, wetting fronts not at the interface) @@ -1498,11 +1951,24 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ // double free_drainage_demand = 0; // prior mass = mass contained in the current old wetting front - double prior_mass = current_old->depth_cm * (current_old->theta - next_old->theta); + double next_reference_theta = lgar_front_changed_by_interflow_stack(next, interflow_stack_changed_by_front) ? next->theta : next_old->theta; + double prior_mass = current_old->depth_cm * (current_old->theta - next_reference_theta); if (wf_free_drainage_demand == wf) prior_mass += precip_mass_to_add - (free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes + actual_ET_demand); + double depth_after_movement_cm = current->depth_cm + current->dzdt_cm_per_h * timestep_h; + if (depth_after_movement_cm > column_depth) + depth_after_movement_cm = column_depth + TRUNCATION_DEPTH; + double minimum_theta = calc_theta_from_h(interflow_psi_cap_cm, vg_a, vg_m, vg_n, theta_e, theta_r); + double minimum_prior_mass_cm = depth_after_movement_cm * (minimum_theta - next->theta); + double applied_interflow_flux_cm = lgar_apply_interflow_flux_to_prior_mass(interflow_flux_cm_by_front[wf], &prior_mass, + minimum_prior_mass_cm, + interflow_subtimestep_cm); + if (applied_interflow_flux_cm > 0.0 && verbosity.compare("high") == 0) { + printf("Applied interflow to WF %d mass balance: %.10e cm\n", wf, applied_interflow_flux_cm); + } + current->depth_cm += current->dzdt_cm_per_h * timestep_h; /* condition to bound the wetting front depth, if depth of a wf, at this timestep, @@ -1528,7 +1994,14 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ current = listDeleteFront(current->front_num, head, soil_type, soil_properties); current = next; double mass_after_theta_went_below_theta_r = lgar_calc_mass_bal(cum_layer_thickness_cm, *head); - *AET_demand_cm = *AET_demand_cm - fabs(mass_before_theta_went_below_theta_r - mass_after_theta_went_below_theta_r); + // Reduce free drainage first so deleting this tiny front does not force negative AET. + double removal_correction_cm = fabs(mass_before_theta_went_below_theta_r - mass_after_theta_went_below_theta_r); + double free_drainage_reduction_cm = fmin(removal_correction_cm, fmax(free_drainage_demand, 0.0)); + free_drainage_demand -= free_drainage_reduction_cm; + if (free_drainage_subtimestep_cm != NULL) + *free_drainage_subtimestep_cm -= free_drainage_reduction_cm; + removal_correction_cm -= free_drainage_reduction_cm; + *AET_demand_cm = *AET_demand_cm - removal_correction_cm; actual_ET_demand = *AET_demand_cm; if (verbosity.compare("high") == 0) { printf("Deleting WF that will go below theta_r (after)...\n"); @@ -1572,14 +2045,16 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ double psi_cm_old = current_old->psi_cm; - double psi_cm_below_old = current_old->next->psi_cm; + bool next_changed_by_interflow_stack = lgar_front_changed_by_interflow_stack(next, interflow_stack_changed_by_front); + double psi_cm_below_old = next_changed_by_interflow_stack ? next->psi_cm : current_old->next->psi_cm; double psi_cm = current->psi_cm; double psi_cm_below = next->psi_cm; // mass = delta(depth) * delta(theta) // = difference in current and next wetting front thetas times depth of the current wetting front - double prior_mass = (current_old->depth_cm - cum_layer_thickness_cm[layer_num-1]) * (current_old->theta - next_old->theta); + double next_reference_theta = next_changed_by_interflow_stack ? next->theta : next_old->theta; + double prior_mass = (current_old->depth_cm - cum_layer_thickness_cm[layer_num-1]) * (current_old->theta - next_reference_theta); double new_mass = (current->depth_cm - cum_layer_thickness_cm[layer_num-1]) * (current->theta - next->theta); // compute mass in the layers above the current wetting front @@ -1620,6 +2095,15 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ if (wf_free_drainage_demand == wf) prior_mass += precip_mass_to_add - (free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes + actual_ET_demand); + + double minimum_prior_mass_cm = lgar_prior_mass_for_psi(interflow_psi_cap_cm, layer_num, delta_thetas, + delta_thickness, soil_type, soil_properties); + double applied_interflow_flux_cm = lgar_apply_interflow_flux_to_prior_mass(interflow_flux_cm_by_front[wf], &prior_mass, + minimum_prior_mass_cm, + interflow_subtimestep_cm); + if (applied_interflow_flux_cm > 0.0 && verbosity.compare("high") == 0) { + printf("Applied interflow to WF %d mass balance: %.10e cm\n", wf, applied_interflow_flux_cm); + } // theta mass balance computes new theta that conserves the mass; new theta is assigned to the current wetting front double theta_new = lgar_theta_mass_balance(layer_num, soil_num, psi_cm, new_mass, prior_mass, precip_mass_to_add, AET_demand_cm, delta_thetas, delta_thickness, soil_type, soil_properties); @@ -1649,7 +2133,8 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ int soil_num_k1 = soil_type[wf_free_drainage->layer_num]; double theta_e_k1 = soil_properties[soil_num_k1].theta_e; - double mass_timestep = (old_mass + precip_mass_to_add) - (actual_ET_demand + free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes); + double interflow_for_mass_balance_cm = interflow_subtimestep_cm == NULL ? 0.0 : *interflow_subtimestep_cm; + double mass_timestep = (old_mass + precip_mass_to_add) - (actual_ET_demand + free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes + interflow_for_mass_balance_cm); assert (old_mass > 0.0); @@ -1745,7 +2230,8 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ //in layered soils, this can cause a mass balance error. It is fairly rare and only seems to impact cases where the model domain is entirely saturated, which shouldn't happen when LGAR is applied in the correct environment / with sufficient layer thicknesses. if (break_flag) { current_mass = lgar_calc_mass_bal(cum_layer_thickness_cm, *head); - mass_timestep = (old_mass + precip_mass_to_add) - (actual_ET_demand + free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes); + interflow_for_mass_balance_cm = interflow_subtimestep_cm == NULL ? 0.0 : *interflow_subtimestep_cm; + mass_timestep = (old_mass + precip_mass_to_add) - (actual_ET_demand + free_drainage_demand + mass_correction_for_cached_free_drainage_fluxes + interflow_for_mass_balance_cm); mass_balance_error = mass_timestep - current_mass; bottom_boundary_flux_cm += mass_balance_error; } @@ -1861,7 +2347,8 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ current->depth_cm = cum_layer_thickness_cm[1]; } } - return(bottom_boundary_flux_cm); + return(bottom_boundary_flux_cm); //the small amount of water that left the LGARTO domain through the lower boundary because it was truncated (WF got slightly too deep), where exceeding the lower boundary by only a small amount is enforced, + //or the small amount of water needed to close mass balance in the very rare event break_flag was set to true when updating a completely saturated WF. Note that this is different than free drainage, which is explicitly accounted for elsewhere. } @@ -1872,6 +2359,38 @@ extern double lgar_move_wetting_fronts(double timestep_h, double *free_drainage_ */ // ############################################################################################ +static bool lgar_wetting_fronts_can_merge(struct wetting_front *current) +{ + const double theta_tolerance = 1.0e-12; + + if (current == NULL || current->next == NULL || current->next->next == NULL) + return false; + + struct wetting_front *next = current->next; + struct wetting_front *next_to_next = next->next; + + if (current->depth_cm <= next->depth_cm) + return false; + if (current->layer_num != next->layer_num) + return false; + if (next->to_bottom) + return false; + if (current->theta <= next->theta + theta_tolerance) + return false; + if (next->theta <= next_to_next->theta + theta_tolerance) + return false; + + double denominator = current->theta - next_to_next->theta; + if (denominator <= theta_tolerance) + return false; + + double current_mass_this_layer = current->depth_cm * (current->theta - next->theta) + + next->depth_cm * (next->theta - next_to_next->theta); + double merged_depth_cm = current_mass_this_layer / denominator; + + return isfinite(merged_depth_cm) && merged_depth_cm > 0.0; +} + extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, struct wetting_front** head, struct soil_properties_ *soil_properties) { @@ -1908,13 +2427,13 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str // 'current->depth_cm > next->depth_cm' ensures that merging is needed // 'current->layer_num == next->layer_num' ensures wetting fronts are in the same layer // '!next->to_bottom' ensures that the next wetting front is not the deepest wetting front in the layer - if ( (current->depth_cm > next->depth_cm) && (current->layer_num == next->layer_num) && !next->to_bottom) { + // The theta checks protect the merge formula from dry-over-wet or equal-theta states, which are + // handled by later correction passes. + if (lgar_wetting_fronts_can_merge(current)) { double current_mass_this_layer = current->depth_cm * (current->theta - next->theta) + next->depth_cm*(next->theta - next_to_next->theta); current->depth_cm = current_mass_this_layer / (current->theta - next_to_next->theta); - assert (current->depth_cm > 0.0); - layer_num = current->layer_num; soil_num = soil_type[layer_num]; theta_e = soil_properties[soil_num].theta_e; @@ -1940,6 +2459,8 @@ extern void lgar_merge_wetting_fronts(int *soil_type, double *frozen_factor, str printf ("Deleting wetting front (after) ... \n"); listPrint(*head); } + + break; } current = current->next; @@ -2048,6 +2569,8 @@ extern void lgar_wetting_fronts_cross_layer_boundary(int num_layers, if (depth_new>cum_layer_thickness_cm[num_layers] && next->layer_num==num_layers){ theta_correction_necessary = true; } + + break; } @@ -2289,6 +2812,7 @@ extern void lgar_fix_dry_over_wet_wetting_fronts(double *mass_change, double* cu double mass_after = lgar_calc_mass_bal(cum_layer_thickness_cm, *head); *mass_change += (mass_after - prior_mass); + break; } current = current->next; @@ -2550,12 +3074,13 @@ extern void lgar_create_surficial_front(int num_layers, double *ponded_depth_cm, if (current->next!=NULL){// sometimes a new WF immediately has to merge with another WF bool had_to_merge = false; - while ( (current->depth_cm > current->next->depth_cm) && (current->layer_num == current->next->layer_num) && !(current->next->to_bottom)){ + while (lgar_wetting_fronts_can_merge(current)){ // Technically this should be replaced with the function that iteratively checks for merging, layer crossing, lower boundary crossing, and dry over wet, // but because all we are doing is adding a new WF onto a linked list that will not need correction because correction was just done on it, it could be that all we need to do here is merge // because the resulting depths should all be in the top layer lgar_merge_wetting_fronts(soil_type, frozen_factor, head, soil_properties); had_to_merge = true; + current = *head; } if (had_to_merge){ //pretty sure this is not necessary but keeping it in lgar_wetting_fronts_cross_layer_boundary(num_layers, cum_layer_thickness_cm, soil_type, frozen_factor, @@ -3078,6 +3603,13 @@ extern double lgar_theta_mass_balance(int layer_num, int soil_num, double psi_cm } if ( (psi_cm_loc > PSI_UPPER_LIM) && (iter > first_speedup_thresh) ){ //unrealistic pressures, but there are some cases where convergence is possible even at large psi values, and there is a case where AET, free drainage, or WF movement can bring psi above PSI_UPPER_LIM, so we do want to allow a few iterations + // Return the dry-limit mass, not the slightly drier mass from the last trial step. + psi_cm_loc = PSI_UPPER_LIM; + theta = calc_theta_from_h(psi_cm_loc, soil_properties[soil_num].vg_alpha_per_cm, soil_properties[soil_num].vg_m, + soil_properties[soil_num].vg_n,soil_properties[soil_num].theta_e, + soil_properties[soil_num].theta_r); + new_mass = lgar_prior_mass_for_psi(psi_cm_loc, layer_num, delta_theta, delta_thickness, soil_type, soil_properties); + delta_mass = fabs(new_mass - prior_mass); break; } @@ -3134,7 +3666,7 @@ extern int lgarto_correction_type_surf(int num_layers, double* cum_layer_thickne if (next!=NULL){ // if ( (current->is_WF_GW==0) && (next->is_WF_GW==0) && (current->theta>next->theta) && (current->depth_cm > next->depth_cm) && (current->layer_num == next->layer_num) && (!next->to_bottom) ){ - if ( (current->theta>next->theta) && (current->depth_cm > next->depth_cm) && (current->layer_num == next->layer_num) && (!next->to_bottom) ){ + if (lgar_wetting_fronts_can_merge(current)){ correction_type_surf = 1; //this is surface-surface WF merging break; } diff --git a/tests/main_unit_test_bmi.cxx b/tests/main_unit_test_bmi.cxx index 9e3974d..f93b711 100644 --- a/tests/main_unit_test_bmi.cxx +++ b/tests/main_unit_test_bmi.cxx @@ -589,9 +589,14 @@ int main(int argc, char *argv[]) double Ksat_1; double Ksat_2; double field_capacity; - double a; - double b; + double a_con_res; + double b_con_res; + double a_con_res_legacy; double frac_to_CR; + double interflow_psi_threshold; + double interflow_factor; + double interflow_psi_threshold_legacy; + double interflow_factor_legacy; double spf_factor; // double smcmax_set[] = {0.3513, 0.3773, 0.3617}; @@ -607,9 +612,11 @@ int main(int argc, char *argv[]) double hydraulic_conductivity_1_set = 0.446; double hydraulic_conductivity_2_set = 0.0743; double field_capacity_set = 103.3; - double a_set = 0.003; - double b_set = 2.5; + double a_con_res_set = 0.003; + double b_con_res_set = 2.5; double frac_to_CR_set = 0.05; + double interflow_psi_threshold_set = 1500.0; + double interflow_factor_set = 1.0; double spf_factor_set = 0.92; // Get the initial values set through the config file @@ -625,9 +632,11 @@ int main(int argc, char *argv[]) model_calib.GetValue("van_genuchten_alpha_2", &vg_alpha_2); model_calib.GetValue("hydraulic_conductivity_2", &Ksat_2); model_calib.GetValue("field_capacity", &field_capacity); - model_calib.GetValue("a", &a); - model_calib.GetValue("b", &b); + model_calib.GetValue("a_con_res", &a_con_res); + model_calib.GetValue("b_con_res", &b_con_res); model_calib.GetValue("frac_to_CR", &frac_to_CR); + model_calib.GetValue("interflow_psi_threshold", &interflow_psi_threshold); + model_calib.GetValue("interflow_factor", &interflow_factor); model_calib.GetValue("spf_factor", &spf_factor); // for (int i=0; i < num_layers; i++) @@ -643,8 +652,8 @@ int main(int argc, char *argv[]) // printf("a: %lf \n", vg_alpha_2); // printf("a: %lf \n", Ksat_2); // printf("field_capacity: %lf \n", field_capacity); - // printf("a: %lf \n", a); - // printf("a: %lf \n", b); + // printf("a_con_res: %lf \n", a_con_res); + // printf("b_con_res: %lf \n", b_con_res); // printf("a: %lf \n", frac_to_CR); // printf("a: %lf \n", spf_factor); @@ -665,9 +674,11 @@ int main(int argc, char *argv[]) model_calib.SetValue("van_genuchten_alpha_2", &van_genuchten_alpha_2_set); model_calib.SetValue("hydraulic_conductivity_2", &hydraulic_conductivity_2_set); model_calib.SetValue("field_capacity", &field_capacity_set); - model_calib.SetValue("a", &a_set); - model_calib.SetValue("b", &b_set); + model_calib.SetValue("a_con_res", &a_con_res_set); + model_calib.SetValue("b_con_res", &b_con_res_set); model_calib.SetValue("frac_to_CR", &frac_to_CR_set); + model_calib.SetValue("interflow_psi_threshold", &interflow_psi_threshold_set); + model_calib.SetValue("interflow_factor", &interflow_factor_set); model_calib.SetValue("spf_factor", &spf_factor_set); // // get the new/updated values @@ -686,9 +697,14 @@ int main(int argc, char *argv[]) model_calib.GetValue("van_genuchten_alpha_2", &vg_alpha_2); model_calib.GetValue("hydraulic_conductivity_2", &Ksat_2); model_calib.GetValue("field_capacity", &field_capacity); - model_calib.GetValue("a", &a); - model_calib.GetValue("b", &b); + model_calib.GetValue("a_con_res", &a_con_res); + model_calib.GetValue("b_con_res", &b_con_res); + model_calib.GetValue("a", &a_con_res_legacy); model_calib.GetValue("frac_to_CR", &frac_to_CR); + model_calib.GetValue("interflow_psi_threshold", &interflow_psi_threshold); + model_calib.GetValue("interflow_factor", &interflow_factor); + model_calib.GetValue("lateral_flow_psi_threshold", &interflow_psi_threshold_legacy); + model_calib.GetValue("lateral_flow_factor", &interflow_factor_legacy); model_calib.GetValue("spf_factor", &spf_factor); @@ -760,30 +776,65 @@ int main(int argc, char *argv[]) throw std::runtime_error(errMsg.str()); } - if (fabs(a - a_set) > 1.E-5) { + if (fabs(a_con_res - a_con_res_set) > 1.E-5) { std::stringstream errMsg; - errMsg << "Mismatch between a calibrated values set and get "<< a_set<<" "<< a + errMsg << "Mismatch between a_con_res calibrated values set and get "<< a_con_res_set<<" "<< a_con_res << " which is unexpected. \n"; throw std::runtime_error(errMsg.str()); } - if (fabs(b - b_set) > 1.E-5) { + if (fabs(b_con_res - b_con_res_set) > 1.E-5) { std::stringstream errMsg; - errMsg << "Mismatch between a calibrated values set and get "<< a_set<<" "<< a + errMsg << "Mismatch between b_con_res calibrated values set and get "<< b_con_res_set<<" "<< b_con_res + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + + if (fabs(a_con_res_legacy - a_con_res_set) > 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between a legacy alias and a_con_res calibrated values "<< a_con_res_set<<" "<< a_con_res_legacy << " which is unexpected. \n"; throw std::runtime_error(errMsg.str()); } if (fabs(frac_to_CR - frac_to_CR_set) > 1.E-5) { std::stringstream errMsg; - errMsg << "Mismatch between a calibrated values set and get "<< a_set<<" "<< a + errMsg << "Mismatch between frac_to_CR calibrated values set and get "<< frac_to_CR_set<<" "<< frac_to_CR + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + + if (fabs(interflow_psi_threshold - interflow_psi_threshold_set) > 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between interflow_psi_threshold calibrated values set and get "<< interflow_psi_threshold_set<<" "<< interflow_psi_threshold + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + + if (fabs(interflow_factor - interflow_factor_set) > 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between interflow_factor calibrated values set and get "<< interflow_factor_set<<" "<< interflow_factor + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + + if (fabs(interflow_psi_threshold_legacy - interflow_psi_threshold_set) > 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between lateral_flow_psi_threshold legacy alias and interflow_psi_threshold calibrated values "<< interflow_psi_threshold_set<<" "<< interflow_psi_threshold_legacy + << " which is unexpected. \n"; + throw std::runtime_error(errMsg.str()); + } + + if (fabs(interflow_factor_legacy - interflow_factor_set) > 1.E-5) { + std::stringstream errMsg; + errMsg << "Mismatch between lateral_flow_factor legacy alias and interflow_factor calibrated values "<< interflow_factor_set<<" "<< interflow_factor_legacy << " which is unexpected. \n"; throw std::runtime_error(errMsg.str()); } if (fabs(spf_factor - spf_factor_set) > 1.E-5) { std::stringstream errMsg; - errMsg << "Mismatch between a calibrated values set and get "<< a_set<<" "<< a + errMsg << "Mismatch between spf_factor calibrated values set and get "<< spf_factor_set<<" "<< spf_factor << " which is unexpected. \n"; throw std::runtime_error(errMsg.str()); } @@ -833,9 +884,11 @@ int main(int argc, char *argv[]) printf("vg_alpha_1 = %lf \n", vg_alpha_1); printf("vg_alpha_2 = %lf \n", vg_alpha_2); printf("field_capacity = %lf \n", field_capacity); - printf("a = %lf \n", a); - printf("b = %lf \n", b); + printf("a_con_res = %lf \n", a_con_res); + printf("b_con_res = %lf \n", b_con_res); printf("frac_to_CR = %lf \n", frac_to_CR); + printf("interflow_psi_threshold = %lf \n", interflow_psi_threshold); + printf("interflow_factor = %lf \n", interflow_factor); printf("spf_factor = %lf \n", spf_factor); std::cout<<"| *************************************** \n"; std::cout<<"| LASAM Calibration test passed? YES \n";