diff --git a/NEWS.md b/NEWS.md
index 5791bf5b4..7cbd5bbc9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,9 +1,13 @@
# NEWS
# SOILWAT2 v8.1.0-devel
-* This version produces nearly identical simulation output as previously.
- Small deviations arise from replacing all remaining variables of type float
- with type double (see commit 62237ae on 2024-July-30).
+* This version produces similar but not identical simulation output
+ as previously because of the following changes:
+ * Small deviations arise from replacing all remaining variables of
+ type float with type double (see commit 62237ae on 2024-July-30).
+ * Saturated percolation is now limited which leads to different outcomes
+ during periods of high infiltration (e.g., snowmelt) and during conditions
+ of low hydraulic conductivity (e.g., frozen soils, sapric organic matter).
* The two models of SOILWAT2 can now be compiled with the following flags:
* `make CPPFLAGS=-DSWTXT` (or as previously `make all`) for txt-based
@@ -11,6 +15,35 @@
* Tests now require `c++17` and utilize `googletest` `v1.15.2` (issue #427).
+* SOILWAT2 can now represent the influence of soil organic matter on the
+ soil water retention curve and the saturated hydraulic conductivity
+ parameter (#397; @dschlaep). The implemented approach first determines
+ organic matter properties for the soil layers assuming fibric peat
+ characteristics at the soil surface and characteristics of sapric peat
+ at a user-specified depth. Then, bulk soil parameters of the soil water
+ retention curve are estimated as linear combinations of properties for
+ the mineral soil component and of properties for the organic matter soil
+ component using the proportion of organic matter in the bulk soil as weights.
+ The bulk soil saturated hydraulic conductivity parameter accounts for
+ flow pathways through organic matter above a threshold and assumes
+ conductivities through mineral and organic components in series outside of
+ those pathways.
+
+* Saturated percolation is now limited. The upper bound is a function based on
+ the saturated hydraulic conductivity parameter
+ (which includes effects of organic matter), frozen soils, and a
+ user-specified `"permeability"` factor.
+
+## Changes to inputs
+* New input via `"siteparam.in"` to specify the depth at which characteristics
+ of organic matter have completely switched from fibric to sapric peat
+ (default is 50 cm).
+* New input via `"soils.in"` to provide the proportion of soil organic matter
+ to bulk soil by weight for each soil layer.
+* New input via `"swrc_params*.in"` to provide parameter values of the
+ soil water retention curve representing fibric and sapric peat.
+ Note: Some parameter values for the `"FXW"` SWRC are missing.
+
# SOILWAT2 v8.0.0
* Simulation output remains the same as the previous version.
diff --git a/include/SW_Flow_lib.h b/include/SW_Flow_lib.h
index f8b311600..734ac8c31 100644
--- a/include/SW_Flow_lib.h
+++ b/include/SW_Flow_lib.h
@@ -95,6 +95,7 @@ void infiltrate_water_high(
unsigned int nlyrs,
const double swcfc[],
double swcsat[],
+ double ksat[],
const double impermeability[],
double *standingWater,
double lyrFrozen[]
diff --git a/include/SW_Site.h b/include/SW_Site.h
index f97183486..21e01e8f5 100644
--- a/include/SW_Site.h
+++ b/include/SW_Site.h
@@ -205,7 +205,7 @@ unsigned int encode_str2ptf(char *ptf_name);
void SWRC_PTF_estimate_parameters(
unsigned int ptf_type,
- double *swrcp,
+ double *swrcpMineralSoil,
double sand,
double clay,
double gravel,
@@ -214,7 +214,7 @@ void SWRC_PTF_estimate_parameters(
);
void SWRC_PTF_Cosby1984_for_Campbell1974(
- double *swrcp, double sand, double clay
+ double *swrcpMineralSoil, double sand, double clay
);
@@ -240,6 +240,7 @@ double SW_swcBulk_saturated(
unsigned int ptf_type,
double sand,
double clay,
+ double fom,
LOG_INFO *LogInfo
);
@@ -252,23 +253,36 @@ double SW_swcBulk_minimum(
double ui_sm_min,
double sand,
double clay,
+ double fom,
double swcBulk_sat,
double SWCMinVal,
LOG_INFO *LogInfo
);
+
void PTF_Saxton2006(
- double *theta_sat, double sand, double clay, LOG_INFO *LogInfo
+ double *theta_sat, double sand, double clay, double fom, LOG_INFO *LogInfo
);
void PTF_RawlsBrakensiek1985(
double *theta_min,
double sand,
double clay,
+ double fom,
double porosity,
LOG_INFO *LogInfo
);
+void SWRC_bulkSoilParameters(
+ unsigned int swrc_type,
+ double *swrcp,
+ const double *swrcpMS,
+ double swrcpOM[][SWRC_PARAM_NMAX],
+ double fom,
+ double depthSapric,
+ double depthT,
+ double depthB
+);
double calculate_soilBulkDensity(double matricDensity, double fractionGravel);
@@ -322,6 +336,7 @@ void set_soillayers(
const double *pclay,
const double *imperm,
const double *soiltemp,
+ const double *pom,
int nRegions,
double *regionLowerBounds,
LOG_INFO *LogInfo
diff --git a/include/SW_datastructs.h b/include/SW_datastructs.h
index 19fa22ebe..0f9260d7c 100644
--- a/include/SW_datastructs.h
+++ b/include/SW_datastructs.h
@@ -299,8 +299,9 @@ typedef struct {
char site_swrc_name[64], site_ptf_name[64];
- Bool site_has_swrcp; /**< Are `swrcp` already (TRUE) or not yet estimated
- (FALSE)? */
+ /** Are `swrcp` of the mineral soil already (TRUE) or not yet estimated
+ (FALSE)? */
+ Bool site_has_swrcpMineralSoil;
/* transpiration regions shallow, moderately shallow, */
/* deep and very deep. units are in layer numbers. */
@@ -309,10 +310,13 @@ typedef struct {
SWCWetVal, /* value for a "wet" day, */
SWCMinVal; /* lower bound on swc. */
- /* bulk = relating to the whole soil, i.e., matric + rock/gravel/coarse
- * fragments */
- /* matric = relating to the < 2 mm fraction of the soil, i.e., sand, clay,
- * and silt */
+ /* Soil components
+ * bulk = relating to the whole soil
+ i.e., matric + coarse fragment (gravel)
+ * matric component = relating to the < 2 mm fraction of the soil
+ * mineral component = sand, clay, silt
+ * organic component = organic matter
+ */
double
/* Inputs */
@@ -335,6 +339,8 @@ typedef struct {
(0=permeable, 1=impermeable) */
avgLyrTempInit[MAX_LAYERS], /* initial soil temperature for each soil
layer */
+ /** Organic matter content as weight fraction of bulk soil [g g-1] */
+ fractionWeight_om[MAX_LAYERS],
/* Derived soil characteristics */
soilMatric_density[MAX_LAYERS], /* matric soil density of the < 2 mm
@@ -360,7 +366,11 @@ typedef struct {
/* Saxton et al. 2006 */
swcBulk_saturated[MAX_LAYERS]; /* saturated bulk SWC [cm] */
- // currently, not used;
+
+ /** Saturated hydraulic conductivity of the bulk soil */
+ double ksat[MAX_LAYERS];
+
+ // currently, not used;
// Saxton2006_K_sat_matric, /* saturated matric conductivity [cm / day] */
// Saxton2006_K_sat_bulk, /* saturated bulk conductivity [cm / day] */
// Saxton2006_fK_gravel, /* gravel-correction factor for conductivity [1] */
@@ -368,6 +378,9 @@ typedef struct {
double depths[MAX_LAYERS]; // soil layer depths of SoilWat soil
+ /** Depth [cm] at which soil properties reach values of sapric peat */
+ double depthSapric;
+
/* Soil water retention curve (SWRC) */
unsigned int swrc_type[MAX_LAYERS], /**< Type of SWRC (see #swrc2str) */
ptf_type[MAX_LAYERS]; /**< Type of PTF (see #ptf2str) */
@@ -377,9 +390,19 @@ typedef struct {
`swrcp` but we need to loop over soil layers for every
vegetation type in `my_transp_rng`
*/
+ /** SWRC parameters of the bulk soil
+ (weighted average of mineral and organic SWRC).
+
+ Note: parameter interpretation varies with selected SWRC,
+ see `SWRC_check_parameters()`
+ */
double swrcp[MAX_LAYERS][SWRC_PARAM_NMAX];
- /**< Parameters of SWRC: parameter interpretation varies with selected SWRC,
- * see `SWRC_check_parameters()` */
+ /** SWRC parameters of the mineral soil component */
+ double swrcpMineralSoil[MAX_LAYERS][SWRC_PARAM_NMAX];
+ /** SWRC parameters of the organic soil component
+ for (1) fibric and (2) sapric peat. */
+ double swrcpOM[2][SWRC_PARAM_NMAX];
+
LyrIndex my_transp_rgn[NVEGTYPES][MAX_LAYERS]; /* which transp zones from
Site am I in? */
diff --git a/src/SW_Flow.c b/src/SW_Flow.c
index ea87f435f..a12b1eb63 100644
--- a/src/SW_Flow.c
+++ b/src/SW_Flow.c
@@ -486,6 +486,7 @@ void SW_Water_Flow(SW_RUN *sw, LOG_INFO *LogInfo) {
n_layers,
sw->Site.swcBulk_fieldcap,
sw->Site.swcBulk_saturated,
+ sw->Site.ksat,
sw->Site.impermeability,
&UpNeigh_standingWater,
sw->SoilWat.lyrFrozen
@@ -517,6 +518,7 @@ void SW_Water_Flow(SW_RUN *sw, LOG_INFO *LogInfo) {
n_layers,
sw->Site.swcBulk_fieldcap,
sw->Site.swcBulk_saturated,
+ sw->Site.ksat,
sw->Site.impermeability,
standingWaterToday,
sw->SoilWat.lyrFrozen
diff --git a/src/SW_Flow_lib.c b/src/SW_Flow_lib.c
index 8dd254bea..ce050dae9 100644
--- a/src/SW_Flow_lib.c
+++ b/src/SW_Flow_lib.c
@@ -376,6 +376,7 @@ void litter_intercepted_water(
(cm H2O).
@param swcsat Soilwater content in each layer at saturation
(m3H2O).
+@param ksat Saturated hydraulic conductivity [cm day-1]
@param impermeability Impermeability measures for each layer.
@param *standingWater Remaining water on the surface
(m3H2O).
@@ -404,6 +405,7 @@ void infiltrate_water_high(
unsigned int nlyrs,
const double swcfc[],
double swcsat[],
+ double ksat[],
const double impermeability[],
double *standingWater,
double lyrFrozen[]
@@ -411,7 +413,6 @@ void infiltrate_water_high(
unsigned int i;
unsigned int j;
- double d[MAX_LAYERS] = {0};
double push;
double ksat_rel;
@@ -427,19 +428,21 @@ void infiltrate_water_high(
ksat_rel = 1.;
}
+ ksat_rel *= (1. - impermeability[i]);
+
/* calculate potential saturated percolation */
- d[i] =
- fmax(0., ksat_rel * (1. - impermeability[i]) * (swc[i] - swcfc[i]));
- drain[i] = d[i];
+ drain[i] = ksat_rel * fmin(ksat[i], swc[i] - swcfc[i]);
+ drain[i] = fmax(0., drain[i]);
+
if (i < nlyrs - 1) {
/* percolate up to next-to-last layer */
- swc[i + 1] += d[i];
- swc[i] -= d[i];
+ swc[i + 1] += drain[i];
+ swc[i] -= drain[i];
} else {
/* percolate last layer */
- (*drainout) = d[i];
- swc[i] -= (*drainout);
+ (*drainout) = drain[i];
+ swc[i] -= drain[i];
}
}
diff --git a/src/SW_Site.c b/src/SW_Site.c
index 6c7507987..462001d08 100644
--- a/src/SW_Site.c
+++ b/src/SW_Site.c
@@ -139,6 +139,16 @@ const char *const swrc2str[N_SWRCs] = {
"Campbell1974", "vanGenuchten1980", "FXW"
};
+
+/** Index to saturated hydraulic conductivity parameter for each SWRC
+
+@note Code maintenance:
+ - Order must exactly match "indices of `swrc2str`" (see also \ref swrc2str)
+ - See details in section \ref swrc_ptf
+*/
+const unsigned int swrcp2ksat[N_SWRCs] = {3, 4, 4};
+
+
/** Character representation of implemented Pedotransfer Functions (PTF)
@note Code maintenance:
@@ -214,6 +224,12 @@ static Bool SW_check_soil_properties(
fval = SW_Site->impermeability[layerno];
errtype = Str_Dup("impermeability", LogInfo);
+ } else if (LT(SW_Site->fractionWeight_om[layerno], 0.) ||
+ GT(SW_Site->fractionWeight_om[layerno], 1.)) {
+ res = swFALSE;
+ fval = SW_Site->fractionWeight_om[layerno];
+ errtype = Str_Dup("organic matter content", LogInfo);
+
} else if (LT(SW_Site->evap_coeff[layerno], 0.) ||
GT(SW_Site->evap_coeff[layerno], 1.)) {
res = swFALSE;
@@ -283,7 +299,7 @@ by user input via #SW_SITE.SWCMinVal as
then calculated as maximum of
`lower_limit_of_theta_min()` and `PTF_RawlsBrakensiek1985()` with
pedotransfer function by Rawls & Brakensiek 1985 \cite rawls1985WmitE
- (independently of the selected SWRC).
+ (independently of the selected SWRC) if there is no organic matter
@param[in] ui_sm_min User input of requested minimum soil moisture,
see SWCMinVal
@@ -292,6 +308,7 @@ by user input via #SW_SITE.SWCMinVal as
@param[in] width Soil layer width [cm]
@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
+@param[in] fom Proportion by weight of organic matter to bulk soil [g g-1]
@param[in] swcBulk_sat Saturated water content of the bulk soil [cm]
@param[in] swrc_type Identification number of selected SWRC
@param[in] *swrcp Vector of SWRC parameters
@@ -306,6 +323,7 @@ static double ui_theta_min(
double width,
double sand,
double clay,
+ double fom,
double swcBulk_sat,
unsigned int swrc_type,
double *swrcp,
@@ -332,9 +350,13 @@ static double ui_theta_min(
&tmp_vwcmin,
sand,
clay,
+ fom,
swcBulk_sat / ((1. - gravel) * width),
LogInfo
);
+ if (LogInfo->stopRun) {
+ return SW_MISSING; // Exit function prematurely due to error
+ }
/* if `PTF_RawlsBrakensiek1985()` was successful, then take max */
if (!missing(tmp_vwcmin)) {
@@ -358,6 +380,17 @@ static double ui_theta_min(
return vwc_min;
}
+/** Extract saturated hydraulic conductivity from SWRC parameters
+
+@param[in] swrc_type Identification number of selected SWRC
+@param[in] *swrcp Vector of SWRC parameters
+
+@return Saturated hydraulic conductivity [cm day-1]
+*/
+static double SWRC_get_ksat(unsigned int swrc_type, double *swrcp) {
+ return swrcp[swrcp2ksat[swrc_type]];
+}
+
/* =================================================== */
/* Global Function Definitions */
/* --------------------------------------------------- */
@@ -418,7 +451,8 @@ unsigned int encode_str2ptf(char *ptf_name) {
See #ptf2str for implemented PTFs.
@param[in] ptf_type Identification number of selected PTF
-@param[out] *swrcp Vector of SWRC parameters to be estimated
+@param[out] *swrcpMineralSoil Vector of SWRC parameters to be estimated
+ for the mineral soil component
@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
@param[in] gravel Coarse fragments (> 2 mm; e.g., gravel)
@@ -430,7 +464,7 @@ See #ptf2str for implemented PTFs.
*/
void SWRC_PTF_estimate_parameters(
unsigned int ptf_type,
- double *swrcp,
+ double *swrcpMineralSoil,
double sand,
double clay,
double gravel,
@@ -439,10 +473,10 @@ void SWRC_PTF_estimate_parameters(
) {
/* Initialize swrcp[] to 0 */
- memset(swrcp, 0, SWRC_PARAM_NMAX * sizeof(swrcp[0]));
+ memset(swrcpMineralSoil, 0, SWRC_PARAM_NMAX * sizeof(swrcpMineralSoil[0]));
if (ptf_type == sw_Cosby1984AndOthers || ptf_type == sw_Cosby1984) {
- SWRC_PTF_Cosby1984_for_Campbell1974(swrcp, sand, clay);
+ SWRC_PTF_Cosby1984_for_Campbell1974(swrcpMineralSoil, sand, clay);
} else {
LogError(LogInfo, LOGERROR, "PTF is not implemented in SOILWAT2.");
@@ -452,13 +486,11 @@ void SWRC_PTF_estimate_parameters(
/**********************************/
/* TODO: remove once PTFs are implemented that utilize gravel */
/* avoiding `error: unused parameter 'gravel' [-Werror=unused-parameter]` */
- if (gravel < 0.) {
- };
+ (void) gravel;
/* TODO: remove once PTFs are implemented that utilize bdensity */
/* avoiding `error: unused parameter 'gravel' [-Werror=unused-parameter]` */
- if (bdensity < 0.) {
- };
+ (void) bdensity;
/**********************************/
}
@@ -478,28 +510,36 @@ but they are not used here.
See `SWRC_SWCtoSWP_Campbell1974()` and `SWRC_SWPtoSWC_Campbell1974()`
for implementation of Campbell's 1974 SWRC (\cite Campbell1974).
-@param[out] *swrcp Vector of SWRC parameters to be estimated
+@param[out] *swrcpMineralSoil Vector of SWRC parameters to be estimated
+ for the mineral soil component
@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
*/
void SWRC_PTF_Cosby1984_for_Campbell1974(
- double *swrcp, double sand, double clay
+ double *swrcpMineralSoil, double sand, double clay
) {
/* Table 4 */
/* Original equations formulated with percent sand (%) and clay (%),
here re-formulated for fraction of sand and clay */
- /* swrcp[0] = psi_saturated: originally formulated as function of silt
- here re-formulated as function of clay */
- swrcp[0] = powe(10.0, -1.58 * sand - 0.63 * clay + 2.17);
- /* swrcp[1] = theta_saturated: originally with units [100 * cm / cm]
- here re-formulated with units [cm / cm] */
- swrcp[1] = -0.142 * sand - 0.037 * clay + 0.505;
- /* swrcp[2] = b */
- swrcp[2] = -0.3 * sand + 15.7 * clay + 3.10;
- /* swrcp[3] = K_saturated: originally with units [inches / day]
- here re-formulated with units [cm / day] */
- swrcp[3] = 2.54 * 24. * powe(10.0, 1.26 * sand - 0.64 * clay - 0.60);
+ /* swrcpMineralSoil[0] = psi_saturated:
+ originally formulated as function of silt
+ here re-formulated as function of clay */
+ swrcpMineralSoil[0] = powe(10.0, -1.58 * sand - 0.63 * clay + 2.17);
+
+ /* swrcpMineralSoil[1] = theta_saturated:
+ originally with units [100 * cm / cm]
+ here re-formulated with units [cm / cm] */
+ swrcpMineralSoil[1] = -0.142 * sand - 0.037 * clay + 0.505;
+
+ /* swrcpMineralSoil[2] = b */
+ swrcpMineralSoil[2] = -0.3 * sand + 15.7 * clay + 3.10;
+
+ /* swrcpMineralSoil[3] = K_saturated:
+ originally with units [inches / day]
+ here re-formulated with units [cm / day] */
+ swrcpMineralSoil[3] =
+ 2.54 * 24. * powe(10.0, 1.26 * sand - 0.64 * clay - 0.60);
}
/**
@@ -512,13 +552,14 @@ Saturated volumetric water content is usually estimated as one of the
SWRC parameters; this is what this function returns.
For historical reasons, if `swrc_name` is "Campbell1974", then a
-`ptf_name` of "Cosby1984AndOthers" will reproduce `SOILWAT2` legacy mode
-(`SOILWAT2` prior to v7.0.0) and return saturated soil water content
-estimated by Saxton et al. 2006 (\cite Saxton2006) PTF instead; `ptf_name` of
+`ptf_name` of "Cosby1984AndOthers" and zero organic matter will reproduce
+`SOILWAT2` legacy mode (`SOILWAT2` prior to v7.0.0) and
+return saturated soil water content estimated by
+Saxton et al. 2006 (\cite Saxton2006) PTF instead; `ptf_name` of
"Cosby1984" will return saturated soil water content estimated by Cosby et
al. 1984 (\cite Cosby1984) PTF.
-The arguments `ptf_type`, `sand`, and `clay` are utilized only if
+The arguments `ptf_type`, `sand`, `clay`, and `fom` are utilized only if
`ptf_name` is "Cosby1984AndOthers" (see #ptf2str).
@param[in] swrc_type Identification number of selected SWRC
@@ -529,6 +570,7 @@ The arguments `ptf_type`, `sand`, and `clay` are utilized only if
@param[in] ptf_type Identification number of selected PTF
@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
+@param[in] fom Proportion by weight of organic matter to bulk soil [g g-1]
@param[out] LogInfo Holds information on warnings and errors
@return Estimated saturated water content of the bulk soil [cm]
@@ -541,6 +583,7 @@ double SW_swcBulk_saturated(
unsigned int ptf_type,
double sand,
double clay,
+ double fom,
LOG_INFO *LogInfo
) {
double theta_sat = SW_MISSING;
@@ -549,7 +592,7 @@ double SW_swcBulk_saturated(
case sw_Campbell1974:
if (ptf_type == sw_Cosby1984AndOthers) {
// Cosby1984AndOthers (backwards compatible)
- PTF_Saxton2006(&theta_sat, sand, clay, LogInfo);
+ PTF_Saxton2006(&theta_sat, sand, clay, fom, LogInfo);
} else {
theta_sat = swrcp[1];
}
@@ -590,7 +633,7 @@ use a realistic lower limit that does not crash the simulation.
The lower limit is determined via `ui_theta_min()` using user input
and is strictly larger than the theoretical SWRC value.
-The arguments `sand`, `clay`, and `swcBulk_sat`
+The arguments `sand`, `clay`, `fom`, and `swcBulk_sat`
are utilized only in legacy mode, i.e., `ptf` equal to
"Cosby1984AndOthers".
@@ -604,6 +647,7 @@ are utilized only in legacy mode, i.e., `ptf` equal to
see #SW_SITE.SWCMinVal
@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
+@param[in] fom Proportion by weight of organic matter to bulk soil [g g-1]
@param[in] swcBulk_sat Saturated water content of the bulk soil [cm]
@param[in] SWCMinVal Lower bound on swc.
@param[out] LogInfo Holds information on warnings and errors
@@ -619,6 +663,7 @@ double SW_swcBulk_minimum(
double ui_sm_min,
double sand,
double clay,
+ double fom,
double swcBulk_sat,
double SWCMinVal,
LOG_INFO *LogInfo
@@ -658,13 +703,12 @@ double SW_swcBulk_minimum(
width,
sand,
clay,
+ fom,
swcBulk_sat,
swrc_type,
swrcp,
- // `(Bool) ptf_type == sw_Cosby1984AndOthers` doesn't work for unit
- // test:
- // error: "no known conversion from 'bool' to 'Bool'"
- ptf_type == sw_Cosby1984AndOthers ? swTRUE : swFALSE,
+ // legacy mode? i.e., consider PTF_RawlsBrakensiek1985()
+ (ptf_type == sw_Cosby1984AndOthers) ? swTRUE : swFALSE,
SWCMinVal,
LogInfo
);
@@ -1017,31 +1061,43 @@ Bool SWRC_check_parameters_for_FXW(double *swrcp, LOG_INFO *LogInfo) {
@brief Saxton et al. 2006 PTFs \cite Saxton2006
to estimate saturated soil water content
-@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
-@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
@param[out] *theta_sat Estimated saturated volumetric water content
of the matric soil [cm/cm]
+@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
+@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
+@param[in] fom Proportion by weight of organic matter to bulk soil [g g-1]
@param[out] LogInfo Holds information on warnings and errors
*/
void PTF_Saxton2006(
- double *theta_sat, double sand, double clay, LOG_INFO *LogInfo
+ double *theta_sat, double sand, double clay, double fom, LOG_INFO *LogInfo
) {
- double OM = 0.;
double theta_33;
double theta_33t;
double theta_S33;
double theta_S33t;
+ double om = 100. * fom; /* equations use percent and not proportion OM */
+
+ if (fom > 0.08) {
+ LogError(
+ LogInfo,
+ LOGERROR,
+ "PTF_Saxton2006(): invalid value of "
+ "organic matter content = %f (must be within 0-0.08)\n",
+ fom
+ );
+ return; // Exit function prematurely due to error
+ }
/* Eq. 2: 33 kPa moisture */
- theta_33t = +0.299 - 0.251 * sand + 0.195 * clay + 0.011 * OM +
- 0.006 * sand * OM - 0.027 * clay * OM + 0.452 * sand * clay;
+ theta_33t = +0.299 - 0.251 * sand + 0.195 * clay + 0.011 * fom +
+ 0.006 * sand * om - 0.027 * clay * om + 0.452 * sand * clay;
theta_33 =
theta_33t + (1.283 * squared(theta_33t) - 0.374 * theta_33t - 0.015);
/* Eq. 3: SAT-33 kPa moisture */
- theta_S33t = +0.078 + 0.278 * sand + 0.034 * clay + 0.022 * OM -
- 0.018 * sand * OM - 0.027 * clay * OM - 0.584 * sand * clay;
+ theta_S33t = +0.078 + 0.278 * sand + 0.034 * clay + 0.022 * om -
+ 0.018 * sand * om - 0.027 * clay * om - 0.584 * sand * clay;
theta_S33 = theta_S33t + (0.636 * theta_S33t - 0.107);
@@ -1067,8 +1123,8 @@ void PTF_Saxton2006(
double R_w, alpha, theta_1500, theta_1500t;
/* Eq. 1: 1500 kPa moisture */
- theta_1500t = +0.031 - 0.024 * sand + 0.487 * clay + 0.006 * OM +
- 0.005 * sand * OM - 0.013 * clay * OM + 0.068 * sand * clay;
+ theta_1500t = +0.031 - 0.024 * sand + 0.487 * clay + 0.006 * fom +
+ 0.005 * sand * om - 0.013 * clay * om + 0.068 * sand * clay;
theta_1500 = theta_1500t + (0.14 * theta_1500t - 0.02);
@@ -1122,6 +1178,7 @@ void PTF_Saxton2006(
@param[in] sand Sand content of the matric soil (< 2 mm fraction) [g/g]
@param[in] clay Clay content of the matric soil (< 2 mm fraction) [g/g]
+@param[in] fom Proportion by weight of organic matter to bulk soil [g g-1]
@param[in] porosity Pore space of the matric soil (< 2 mm fraction) [cm3/cm3]
@param[out] *theta_min Estimated residual volumetric water content
of the matric soil [cm/cm]
@@ -1131,11 +1188,12 @@ void PTF_RawlsBrakensiek1985(
double *theta_min,
double sand,
double clay,
+ double fom,
double porosity,
LOG_INFO *LogInfo
) {
if (GE(clay, 0.05) && LE(clay, 0.6) && GE(sand, 0.05) && LE(sand, 0.7) &&
- GE(porosity, 0.1) && LT(porosity, 1.)) {
+ GE(porosity, 0.1) && LT(porosity, 1.) && ZRO(fom)) {
sand *= 100.;
clay *= 100.;
/* Note: the equation requires sand and clay in units of [100 * g / g];
@@ -1152,14 +1210,161 @@ void PTF_RawlsBrakensiek1985(
);
} else {
- LogError(
- LogInfo,
- LOGWARN,
- "`PTF_RawlsBrakensiek1985()`: sand, clay, or porosity out of valid "
- "range."
- );
-
*theta_min = SW_MISSING;
+
+ if (!ZRO(fom)) {
+ LogError(
+ LogInfo,
+ LOGERROR,
+ "PTF_RawlsBrakensiek1985(): fom = %f (must be 0)",
+ fom
+ );
+ } else {
+ LogError(
+ LogInfo,
+ LOGWARN,
+ "PTF_RawlsBrakensiek1985(): "
+ "sand = %f (must be in 0.05-0.7), "
+ "clay = %f (must be in 0.05-0.6), and/or "
+ "porosity = %f (must be in 0.1-0.1)"
+ "out of valid range",
+ sand,
+ clay,
+ porosity
+ );
+ }
+ }
+}
+
+/** Interpolate organic matter parameter between fibric and sapric
+characteristics
+
+Parameters of the organic component are interpolate between
+fibric peat characteristics assumed at surface and sapric peat at depth.
+
+@param[in] pomFibric Parameter of fibric peat (surface conditions)
+@param[in] pomSapric Parameter of sapric peat (at depth `depthSapric`)
+@param[in] depthSapric Soil depth [`cm`] at which soil properties reach
+ values of sapric peat
+@param[in] depthT Depth at top of soil layer [`cm`]
+@param[in] depthB Depth at bottom of soil layer [`cm`]
+
+@return Parameter value of the organic component
+*/
+static double interpolateFibricSapric(
+ double pomFibric,
+ double pomSapric,
+ double depthSapric,
+ double depthT,
+ double depthB
+) {
+ double res;
+ double b;
+
+ if (ZRO(depthSapric)) {
+ res = pomSapric;
+
+ } else {
+
+ b = (pomSapric - pomFibric) / depthSapric;
+
+ if (pomFibric > pomSapric) {
+ res = 0.5 * (fmax(pomSapric, pomFibric + b * depthT) +
+ fmax(pomSapric, pomFibric + b * depthB));
+
+ } else {
+ res = 0.5 * (fmin(pomSapric, pomFibric + b * depthT) +
+ fmin(pomSapric, pomFibric + b * depthB));
+ }
+ }
+
+ return res;
+}
+
+/** Calculate bulk soil SWRC parameters
+
+Parameters of the bulk soil are calculated as the weighted average of
+parameters of the organic and mineral components.
+
+Parameters of the organic component are interpolate between
+fibric peat characteristics assumed at surface and sapric peat at depth.
+
+Bulk soil saturated hydraulic conductivity accounts for connected pathways
+that only consist of organic matter above a threshold value of organic matter.
+
+@param[in] swrc_type Identification number of selected SWRC
+@param[out] *swrcp Vector of SWRC parameters of the bulk soil
+@param[in] *swrcpMS Vector of SWRC parameters of the mineral component
+@param[in] swrcpOM Array of length two of vectors of SWRC parameters
+ of fibric peat (surface conditions) and
+ of sapric peat (at depth `depthSapric`)
+@param[in] fom Proportion by weight of organic matter to bulk soil
+@param[in] depthSapric Soil depth [`cm`] at which soil properties reach
+ values of sapric peat
+@param[in] depthT Depth at top of soil layer [`cm`]
+@param[in] depthB Depth at bottom of soil layer [`cm`]
+*/
+void SWRC_bulkSoilParameters(
+ unsigned int swrc_type,
+ double *swrcp,
+ const double *swrcpMS,
+ double swrcpOM[][SWRC_PARAM_NMAX],
+ double fom,
+ double depthSapric,
+ double depthT,
+ double depthB
+) {
+ unsigned int k;
+ unsigned int iksat = swrcp2ksat[swrc_type];
+ double pOM;
+ double fperc;
+ double unconnectedKsat;
+
+ static const double percBeta = 0.139; /* percolation exponent */
+ static const double fthreshold = 0.5; /* percolation threshold */
+
+ if (fom > 0.) {
+ /* Has organic matter:
+ interpolate between organic and mineral components */
+
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ /* Interpolate organic parameter from surface to depth conditions */
+ pOM = interpolateFibricSapric(
+ swrcpOM[0][k], swrcpOM[1][k], depthSapric, depthT, depthB
+ );
+
+ if (k == iksat) {
+ /* ksat: account for effects of connected flow pathways */
+
+ if (fom >= fthreshold) {
+ /* (1 - fperc) is the sum of the fraction of mineral soil
+ and the fraction of non-percolating organic component */
+ fperc = fom * pow(1. - fthreshold, -percBeta) *
+ pow(fom - fthreshold, percBeta);
+ } else {
+ fperc = 0.;
+ }
+
+ /* non-connected fraction assumes conductivities
+ through mineral and organic components in series */
+ unconnectedKsat =
+ (fom < 1.) ?
+ 1. / ((1. - fom) / swrcpMS[k] + (fom - fperc) / pOM) :
+ 0.;
+
+ swrcp[k] = fperc * pOM + (1. - fperc) * unconnectedKsat;
+
+ } else {
+ /* other swrcp: weighted average of organic and mineral param */
+ swrcp[k] = fom * pOM + (1. - fom) * swrcpMS[k];
+ }
+ }
+
+ } else {
+ /* No organic matter: bulk soil corresponds to mineral components */
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ swrcp[k] = swrcpMS[k];
+ }
}
}
@@ -1348,13 +1553,14 @@ void SW_SIT_read(
while (GetALine(f, inbuf, MAX_FILENAMESIZE)) {
- strLine = (Bool) (lineno == 35 || lineno == 37 || lineno == 38);
+ strLine = (Bool) (lineno == 35 || lineno == 38 || lineno == 39);
- if (!strLine && lineno <= 38) {
+ if (!strLine && lineno <= 39) {
/* Check to see if the line number contains a double or integer
* value */
- doDoubleConv = (Bool) ((lineno >= 0 && lineno <= 2) ||
- (lineno >= 5 && lineno <= 31));
+ doDoubleConv =
+ (Bool) ((lineno >= 0 && lineno <= 2) ||
+ (lineno >= 5 && lineno <= 31) || lineno == 37);
if (doDoubleConv) {
doubleRes = sw_strtod(inbuf, MyFileName, LogInfo);
@@ -1514,7 +1720,12 @@ void SW_SIT_read(
case 36:
SW_Site->type_soilDensityInput = intRes;
break;
+
case 37:
+ SW_Site->depthSapric = doubleRes;
+ break;
+
+ case 38:
resSNP = snprintf(
SW_Site->site_swrc_name,
sizeof SW_Site->site_swrc_name,
@@ -1534,7 +1745,7 @@ void SW_SIT_read(
goto closeFile;
}
break;
- case 38:
+ case 39:
resSNP = snprintf(
SW_Site->site_ptf_name,
sizeof SW_Site->site_ptf_name,
@@ -1550,12 +1761,12 @@ void SW_SIT_read(
}
SW_Site->site_ptf_type = encode_str2ptf(SW_Site->site_ptf_name);
break;
- case 39:
- SW_Site->site_has_swrcp = itob(intRes);
+ case 40:
+ SW_Site->site_has_swrcpMineralSoil = itob(intRes);
break;
default:
- if (lineno > 39 + MAX_TRANSP_REGIONS) {
+ if (lineno > 40 + MAX_TRANSP_REGIONS) {
break; /* skip extra lines */
}
@@ -1620,6 +1831,18 @@ void SW_SIT_read(
goto closeFile;
}
+ if (LT(SW_Site->depthSapric, 0.)) {
+ LogError(
+ LogInfo,
+ LOGERROR,
+ "%s : depth at which organic matter has characteristics of "
+ "sapric peat = %f (value ranges between 0 and +inf)\n",
+ MyFileName,
+ SW_Site->depthSapric
+ );
+ goto closeFile;
+ }
+
if (too_many_regions) {
LogError(
LogInfo,
@@ -1676,8 +1899,9 @@ void SW_LYR_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
double imperm;
double soiltemp;
double f_gravel;
+ double fom;
char inbuf[MAX_FILENAMESIZE];
- char inDoubleStrs[12][20] = {{'\0'}};
+ char inDoubleStrs[13][20] = {{'\0'}};
double *inDoubleVals[] = {
&dmax,
&soildensity,
@@ -1690,9 +1914,10 @@ void SW_LYR_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
&psand,
&pclay,
&imperm,
- &soiltemp
+ &soiltemp,
+ &fom
};
- const int numDoubleInStrings = 12;
+ const int numDoubleInStrings = 13;
/* note that Files.read() must be called prior to this. */
char *MyFileName = InFiles[eLayers];
@@ -1707,7 +1932,7 @@ void SW_LYR_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
x = sscanf(
inbuf,
- "%19s %19s %19s %19s %19s %19s %19s %19s %19s %19s %19s %19s",
+ "%19s %19s %19s %19s %19s %19s %19s %19s %19s %19s %19s %19s %19s",
inDoubleStrs[0],
inDoubleStrs[1],
inDoubleStrs[2],
@@ -1719,12 +1944,13 @@ void SW_LYR_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
inDoubleStrs[8],
inDoubleStrs[9],
inDoubleStrs[10],
- inDoubleStrs[11]
+ inDoubleStrs[11],
+ inDoubleStrs[12]
);
- /* Check that we have 12 values per layer */
+ /* Check that we have 13 values per layer */
/* Adjust number if new variables are added */
- if (x != 12) {
+ if (x != 13) {
LogError(
LogInfo,
LOGERROR,
@@ -1759,6 +1985,7 @@ void SW_LYR_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
SW_Site->fractionWeightMatric_clay[lyrno] = pclay;
SW_Site->impermeability[lyrno] = imperm;
SW_Site->avgLyrTempInit[lyrno] = soiltemp;
+ SW_Site->fractionWeight_om[lyrno] = fom;
if (lyrno >= MAX_LAYERS) {
LogError(
@@ -1805,6 +2032,7 @@ them from an input file as _read_layers() does)
soil matrix [g3 / g3]
@param[in] imperm Array of size \p nlyrs with impermeability coefficients [0, 1]
@param[in] soiltemp Array of size \p nlyrs with initial soil temperature [C]
+@param[in] pom Array of size \p nlyrs for organic matter [g3 / g3]
@param nRegions The number of transpiration regions to create. Must be between
1 and \ref MAX_TRANSP_REGIONS.
@param[in] regionLowerBounds Array of size \p nRegions containing the lower
@@ -1839,6 +2067,7 @@ void set_soillayers(
const double *pclay,
const double *imperm,
const double *soiltemp,
+ const double *pom,
int nRegions,
double *regionLowerBounds,
LOG_INFO *LogInfo
@@ -1888,6 +2117,7 @@ void set_soillayers(
SW_Site->fractionWeightMatric_clay[lyrno] = pclay[i];
SW_Site->impermeability[lyrno] = imperm[i];
SW_Site->avgLyrTempInit[lyrno] = soiltemp[i];
+ SW_Site->fractionWeight_om[lyrno] = pom[i];
}
@@ -1989,28 +2219,27 @@ void derive_soilRegions(
}
/** Obtain soil water retention curve parameters from disk
- *
- * @param[in,out] SW_Site Struct of type SW_SITE describing the simulated site
- * @param[in] InFiles Array of program in/output files
- * @param[out] LogInfo Holds information on warnings and errors
- *
- */
-void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
- /* Don't read values from disk if they will be estimated via a PTF */
- if (!SW_Site->site_has_swrcp) {
- return;
- }
+The first set (row) of parameters represent fibric peat;
+the second set of parameters represent sapric peat;
+the remaining rows represent parameters of the mineral soil in each soil layer.
+
+@param[in,out] SW_Site Struct of type SW_SITE describing the simulated site
+@param[in] InFiles Array of program in/output files
+@param[out] LogInfo Holds information on warnings and errors
+*/
+void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
FILE *f;
LyrIndex lyrno = 0;
LyrIndex k;
+ /* Indicator if inputs are for organic or mineral soil components */
+ Bool isMineral = swFALSE;
int x;
int index;
double tmp_swrcp[SWRC_PARAM_NMAX];
char inbuf[MAX_FILENAMESIZE];
- char swrcpDoubleStrs[6][20] = {{'\0'}};
- const int numDoubleInStrings = 6;
+ char swrcpDoubleStrs[SWRC_PARAM_NMAX][20] = {{'\0'}};
/* note that Files.read() must be called prior to this. */
char *MyFileName = InFiles[eSWRCp];
@@ -2021,6 +2250,11 @@ void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
}
while (GetALine(f, inbuf, MAX_FILENAMESIZE)) {
+ /* PTF used for mineral swrcp: read only organic parameters from disk */
+ if (isMineral && !SW_Site->site_has_swrcpMineralSoil) {
+ goto closeFile;
+ }
+
x = sscanf(
inbuf,
"%19s %19s %19s %19s %19s %19s",
@@ -2032,8 +2266,7 @@ void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
swrcpDoubleStrs[5]
);
- /* Note: `SW_SIT_init_run()` will call function to check for valid
- * values */
+ /* Note: `SW_SIT_init_run()` will check for valid values */
/* Check that we have n = `SWRC_PARAM_NMAX` values per layer */
if (x != SWRC_PARAM_NMAX) {
@@ -2049,7 +2282,7 @@ void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
}
/* Convert strings to doubles */
- for (index = 0; index < numDoubleInStrings; index++) {
+ for (index = 0; index < SWRC_PARAM_NMAX; index++) {
tmp_swrcp[index] =
sw_strtod(swrcpDoubleStrs[index], MyFileName, LogInfo);
if (LogInfo->stopRun) {
@@ -2058,7 +2291,7 @@ void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
}
/* Check that we are within `SW_Site.n_layers` */
- if (lyrno >= SW_Site->n_layers) {
+ if (isMineral && lyrno >= SW_Site->n_layers) {
LogError(
LogInfo,
LOGERROR,
@@ -2073,10 +2306,21 @@ void SW_SWRC_read(SW_SITE *SW_Site, char *InFiles[], LOG_INFO *LogInfo) {
/* Copy values into structure */
for (k = 0; k < SWRC_PARAM_NMAX; k++) {
- SW_Site->swrcp[lyrno][k] = tmp_swrcp[k];
+ if (isMineral) {
+ SW_Site->swrcpMineralSoil[lyrno][k] = tmp_swrcp[k];
+ } else {
+ SW_Site->swrcpOM[lyrno][k] = tmp_swrcp[k];
+ }
}
lyrno++;
+
+ if (!isMineral && lyrno > 1) {
+ /* Fibric and sapric peat are completed.
+ Now: reset and restart for swrcp of the mineral component */
+ isMineral = swTRUE;
+ lyrno = 0;
+ }
}
closeFile: { CloseFile(&f, LogInfo); }
@@ -2144,7 +2388,7 @@ void SW_SIT_init_run(
/* Check compatibility between selected SWRC and PTF */
- if (!SW_Site->site_has_swrcp) {
+ if (!SW_Site->site_has_swrcpMineralSoil) {
if (!check_SWRC_vs_PTF(
SW_Site->site_swrc_name, SW_Site->site_ptf_name
)) {
@@ -2159,6 +2403,31 @@ void SW_SIT_init_run(
}
}
+ /* Check parameters of organic SWRC */
+ if (!SWRC_check_parameters(
+ SW_Site->site_swrc_type, SW_Site->swrcpOM[0], LogInfo
+ )) {
+ LogError(
+ LogInfo,
+ LOGERROR,
+ "Checks of parameters for SWRC '%s' in fibric peat failed.",
+ swrc2str[SW_Site->site_swrc_type]
+ );
+ return; // Exit function prematurely due to error
+ }
+
+ if (!SWRC_check_parameters(
+ SW_Site->site_swrc_type, SW_Site->swrcpOM[1], LogInfo
+ )) {
+ LogError(
+ LogInfo,
+ LOGERROR,
+ "Checks of parameters for SWRC '%s' in sapric peat failed.",
+ swrc2str[SW_Site->site_swrc_type]
+ );
+ return; // Exit function prematurely due to error
+ }
+
/* Loop over soil layers check variables and calculate parameters */
ForEachSoilLayer(s, SW_Site->n_layers) {
@@ -2229,15 +2498,16 @@ void SW_SIT_init_run(
}
- if (!SW_Site->site_has_swrcp) {
+ if (!SW_Site->site_has_swrcpMineralSoil) {
/* Use pedotransfer function PTF
to estimate parameters of soil water retention curve (SWRC) for
- layer. If `has_swrcp`, then parameters already obtained from disk
- by `SW_SWRC_read()`
+ the mineral component of the layer.
+ If `site_has_swrcpMineralSoil`, then parameters have already
+ been obtained from disk by `SW_SWRC_read()`
*/
SWRC_PTF_estimate_parameters(
SW_Site->ptf_type[s],
- SW_Site->swrcp[s],
+ SW_Site->swrcpMineralSoil[s],
SW_Site->fractionWeightMatric_sand[s],
SW_Site->fractionWeightMatric_clay[s],
SW_Site->fractionVolBulk_gravel[s],
@@ -2249,20 +2519,52 @@ void SW_SIT_init_run(
}
}
- /* Check parameters of selected SWRC */
+ /* Check parameters of mineral soil SWRC */
+ if (!SWRC_check_parameters(
+ SW_Site->swrc_type[s], SW_Site->swrcpMineralSoil[s], LogInfo
+ )) {
+ LogError(
+ LogInfo,
+ LOGERROR,
+ "Checks of parameters for SWRC '%s' "
+ "in the mineral component of layer %d failed.",
+ swrc2str[SW_Site->swrc_type[s]],
+ s + 1
+ );
+ return; // Exit function prematurely due to error
+ }
+
+ /* Calculate bulk soil SWRCp from organic and mineral soil components */
+ SWRC_bulkSoilParameters(
+ SW_Site->swrc_type[s],
+ SW_Site->swrcp[s],
+ SW_Site->swrcpMineralSoil[s],
+ SW_Site->swrcpOM,
+ SW_Site->fractionWeight_om[s],
+ SW_Site->depthSapric,
+ (s > 0) ? SW_Site->depths[s - 1] : 0,
+ SW_Site->depths[s]
+ );
+
+ /* Check parameters of bulk soil SWRC */
if (!SWRC_check_parameters(
SW_Site->swrc_type[s], SW_Site->swrcp[s], LogInfo
)) {
LogError(
LogInfo,
LOGERROR,
- "Checks of parameters for SWRC '%s' in layer %d failed.",
+ "Checks of parameters for SWRC '%s' "
+ "in the bulk soil layer %d failed.",
swrc2str[SW_Site->swrc_type[s]],
s + 1
);
return; // Exit function prematurely due to error
}
+ /* Extract ksat from swrcp */
+ SW_Site->ksat[s] =
+ SWRC_get_ksat(SW_Site->swrc_type[s], SW_Site->swrcp[s]);
+
/* Calculate SWC at field capacity and at wilting point */
SW_Site->swcBulk_fieldcap[s] =
SW_SWRC_SWPtoSWC(0.333, SW_Site, s, LogInfo);
@@ -2305,6 +2607,7 @@ void SW_SIT_init_run(
SW_Site->ptf_type[s],
SW_Site->fractionWeightMatric_sand[s],
SW_Site->fractionWeightMatric_clay[s],
+ SW_Site->fractionWeight_om[s],
LogInfo
);
if (LogInfo->stopRun) {
@@ -2320,6 +2623,7 @@ void SW_SIT_init_run(
SW_Site->SWCMinVal,
SW_Site->fractionWeightMatric_sand[s],
SW_Site->fractionWeightMatric_clay[s],
+ SW_Site->fractionWeight_om[s],
SW_Site->swcBulk_saturated[s],
SW_Site->SWCMinVal,
LogInfo
@@ -2348,16 +2652,13 @@ void SW_SIT_init_run(
}
- /* test validity of values */
+ /* test validity of calculated values */
if (LT(SW_Site->swcBulk_init[s], SW_Site->swcBulk_min[s])) {
LogError(
LogInfo,
LOGERROR,
- "%s : Layer %d\n"
- " calculated `swcBulk_init` (%.4f cm) <= `swcBulk_min` (%.4f "
- "cm).\n"
- " Recheck parameters and try again.\n",
- "Site",
+ "Soil layer %d: swcBulk_init (%f cm) <= "
+ "swcBulk_min (%f cm) but should be larger",
s + 1,
SW_Site->swcBulk_init[s],
SW_Site->swcBulk_min[s]
@@ -2369,11 +2670,8 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGERROR,
- "%s : Layer %d\n"
- " calculated `swcBulk_wiltpt` (%.4f cm) <= `swcBulk_min` "
- "(%.4f cm).\n"
- " Recheck parameters and try again.\n",
- "Site",
+ "Soil layer %d: swcBulk_wiltpt (%f cm) <= "
+ "swcBulk_min (%f cm) but should be larger",
s + 1,
SW_Site->swcBulk_wiltpt[s],
SW_Site->swcBulk_min[s]
@@ -2385,11 +2683,10 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGWARN,
- "%s : Layer %d\n"
- " calculated `swcBulk_halfwiltpt` (%.4f cm / %.2f MPa)\n"
- " <= `swcBulk_min` (%.4f cm / %.2f MPa).\n"
- " `swcBulk_halfwiltpt` was set to `swcBulk_min`.\n",
- "Site",
+ "Soil layer %d: calculated "
+ "swcBulk_halfwiltpt (%f cm / %f MPa) <= "
+ "swcBulk_min (%f cm / %f MPa); therefore, "
+ "swcBulk_halfwiltpt was set to the value of swcBulk_min",
s + 1,
SW_Site->swcBulk_halfwiltpt[s],
-0.1 * SW_SWRC_SWCtoSWP(
@@ -2411,11 +2708,8 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGERROR,
- "%s : Layer %d\n"
- " calculated `swcBulk_wet` (%.4f cm) <= `swcBulk_min` (%.4f "
- "cm).\n"
- " Recheck parameters and try again.\n",
- "Site",
+ "Soil layer %d: calculated swcBulk_wet (%f cm) <= "
+ "swcBulk_min (%f cm)",
s + 1,
SW_Site->swcBulk_wet[s],
SW_Site->swcBulk_min[s]
@@ -2482,13 +2776,12 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGWARN,
- "%s : Layer %d - vegtype %d\n"
- " calculated `swcBulk_atSWPcrit` (%.4f cm / %.4f MPa)\n"
- " <= `swcBulk_min` (%.4f cm / %.4f MPa).\n"
- " `SWcrit` adjusted to %.4f MPa "
- "(and swcBulk_atSWPcrit in every layer will be "
- "re-calculated).\n",
- "Site",
+ "Soil layer %d - vegtype %d: "
+ "calculated swcBulk_atSWPcrit (%f cm / %f MPa) "
+ "<= `swcBulk_min` (%f cm / %f MPa); thus, "
+ "SWcrit was adjusted to %f MPa "
+ "(and swcBulk_atSWPcrit in every soil layer is "
+ "re-calculated)",
s + 1,
k + 1,
SW_Site->swcBulk_atSWPcrit[k][s],
@@ -2528,8 +2821,8 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGERROR,
- ": Top soil layer must be included\n"
- " in %s tranpiration regions.\n",
+ "Top soil layer must be included "
+ "in %s tranpiration region",
key2veg[k]
);
return; // Exit function prematurely due to error
@@ -2537,10 +2830,9 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGERROR,
- ": Transpiration region %d \n"
- " is deeper than the deepest layer with a\n"
- " %s transpiration coefficient > 0 (%d).\n"
- " Please fix the discrepancy and try again.\n",
+ "Transpiration region %d "
+ "is deeper than the deepest layer with a "
+ "%s transpiration coefficient > 0 (%d)",
r + 1,
key2veg[k],
s
@@ -2572,11 +2864,10 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGERROR,
- "%s : Layer %d - vegtype %d\n"
- " calculated `swcBulk_atSWPcrit` (%.4f cm)\n"
- " <= `swcBulk_min` (%.4f cm).\n"
- " even with adjusted `SWcrit` (%.4f MPa).\n",
- "Site",
+ "Soil layer %d - vegtype %d: "
+ "calculated swcBulk_atSWPcrit (%f cm) "
+ "<= swcBulk_min (%f cm) "
+ "despite adjusted SWcrit (%f MPa)",
s + 1,
k + 1,
SW_Site->swcBulk_atSWPcrit[k][s],
@@ -2602,10 +2893,9 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGWARN,
- "%s : Evaporation coefficients were normalized:\n"
- "\tSum of coefficients was %.4f, but must be 1.0. "
+ "Soils: Evaporation coefficients were normalized: "
+ "sum of coefficients was %f, but must be 1.0. "
"New coefficients are:",
- "Site",
evsum
);
@@ -2614,7 +2904,7 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGWARN,
- " Layer %2d : %.4f",
+ " Layer %2d : evco = %.4f",
s + 1,
SW_Site->evap_coeff[s]
);
@@ -2627,10 +2917,9 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGWARN,
- "%s : Transpiration coefficients were normalized for %s:\n"
- "\tSum of coefficients was %.4f, but must be 1.0. "
+ "Soils: Transpiration coefficients were normalized for %s: "
+ "sum of coefficients was %f, but must be 1.0. "
"New coefficients are:",
- "Site",
key2veg[k],
trsum_veg[k]
);
@@ -2641,7 +2930,7 @@ void SW_SIT_init_run(
LogError(
LogInfo,
LOGWARN,
- " Layer %2d : %.4f",
+ " Layer %2d : trco = %.4f",
s + 1,
SW_Site->transp_coeff[k][s]
);
diff --git a/src/SW_SoilWater.c b/src/SW_SoilWater.c
index 9f7ebcade..a9015339a 100644
--- a/src/SW_SoilWater.c
+++ b/src/SW_SoilWater.c
@@ -1264,8 +1264,9 @@ void SW_SWC_end_day(SW_SOILWAT *SW_SoilWat, LyrIndex n_layers) {
/* =================================================== */
LyrIndex i;
- ForEachSoilLayer(i, n_layers) SW_SoilWat->swcBulk[Yesterday][i] =
- SW_SoilWat->swcBulk[Today][i];
+ ForEachSoilLayer(i, n_layers) {
+ SW_SoilWat->swcBulk[Yesterday][i] = SW_SoilWat->swcBulk[Today][i];
+ }
SW_SoilWat->snowpack[Yesterday] = SW_SoilWat->snowpack[Today];
}
diff --git a/tests/example/Input/siteparam.in b/tests/example/Input/siteparam.in
index 9932cee09..641275be2 100644
--- a/tests/example/Input/siteparam.in
+++ b/tests/example/Input/siteparam.in
@@ -82,6 +82,11 @@ RCP85
0
+# --- Organic matter ---
+# Depth [cm] at which characteristics of organic matter have completely switched from fibric to sapric peat
+50
+
+
#--- Soil water retention curve (SWRC) ------
#
# Implemented options (`swrc_name`/`ptf_name`, see `swrc2str[]`/`ptf2str[]`):
@@ -99,7 +104,7 @@ Campbell1974 # Specify soil water retention curve
Cosby1984AndOthers # Specify pedotransfer function
# (if not implemented, then provide SWRC parameters via "swrc_params.in")
-0 # Has SWRC parameters (see `has_swrcp`)?
+0 # Has SWRC parameters for the mineral soil component (see `has_swrcp`)?
# 0: Estimate with specified pedotransfer function
# 1: Use values from "swrc_params.in"
diff --git a/tests/example/Input/soils.in b/tests/example/Input/soils.in
index dec62a939..c950de2cb 100644
--- a/tests/example/Input/soils.in
+++ b/tests/example/Input/soils.in
@@ -15,18 +15,19 @@
# - imperm = (frac, 0-1) proportion of 'impermeability' to water
# percolation(/infiltration/drainage)
# - soiltemp = (C) initial temperature
+# - %om = (frac, 0-1) proportion of organic matter by weight
# USER: the evco and trco columns must sum to 1.0 or they will be normalized.
# The minimum number of soil layers is 1 and the maximum number is `MAX_LAYERS`
# which is defined in `SW_Defines.h`.
-# depth matricd gravel_content evco trco_grass trco_shrub trco_tree trco_forb %sand %clay imperm soiltemp
- 5.0 1.430 0.100 0.813 0.0496 0.134 0.0496 0.134 0.510 0.150 0.000 -1.000
- 10.0 1.410 0.100 0.153 0.0495 0.094 0.0495 0.094 0.440 0.260 0.000 -1.000
- 20.0 1.390 0.100 0.034 0.1006 0.176 0.1006 0.176 0.350 0.410 0.000 -1.000
- 30.0 1.390 0.100 0.000 0.1006 0.175 0.1006 0.175 0.320 0.450 0.000 -1.000
- 40.0 1.380 0.200 0.000 0.1006 0.110 0.1006 0.110 0.310 0.470 0.000 0.000
- 60.0 1.150 0.200 0.000 0.1997 0.109 0.1997 0.109 0.320 0.470 0.000 0.000
- 80.0 1.310 0.200 0.000 0.1997 0.101 0.1997 0.101 0.570 0.280 0.000 1.000
-100.0 1.310 0.200 0.000 0.1997 0.101 0.1997 0.101 0.570 0.280 0.000 1.000
+# depth matricd gravel_content evco trco_grass trco_shrub trco_tree trco_forb %sand %clay imperm soiltemp %om
+ 5.0 1.430 0.100 0.813 0.0496 0.134 0.0496 0.134 0.510 0.150 0.000 -1.000 0.0
+ 10.0 1.410 0.100 0.153 0.0495 0.094 0.0495 0.094 0.440 0.260 0.000 -1.000 0.0
+ 20.0 1.390 0.100 0.034 0.1006 0.176 0.1006 0.176 0.350 0.410 0.000 -1.000 0.0
+ 30.0 1.390 0.100 0.000 0.1006 0.175 0.1006 0.175 0.320 0.450 0.000 -1.000 0.0
+ 40.0 1.380 0.200 0.000 0.1006 0.110 0.1006 0.110 0.310 0.470 0.000 0.000 0.0
+ 60.0 1.150 0.200 0.000 0.1997 0.109 0.1997 0.109 0.320 0.470 0.000 0.000 0.0
+ 80.0 1.310 0.200 0.000 0.1997 0.101 0.1997 0.101 0.570 0.280 0.000 1.000 0.0
+100.0 1.310 0.200 0.000 0.1997 0.101 0.1997 0.101 0.570 0.280 0.000 1.000 0.0
diff --git a/tests/example/Input/swrc_params.in b/tests/example/Input/swrc_params.in
index 744a38f98..e52041419 100644
--- a/tests/example/Input/swrc_params.in
+++ b/tests/example/Input/swrc_params.in
@@ -1,7 +1,7 @@
-#------ Input for Soil Water Retention Curves (by soil layer) ------
+#------ Input for Soil Water Retention Curves ------
-# A table with up to `MAX_LAYERS` rows (soil layers) and 6 columns:
-# - the soil layers must match `soils.in`
+
+# Tables with 6 columns:
# - the interpretation of columns (SWRC parameters) depends on the
# selected SWRC (see `siteparam.in`)
# - unused columns are ignored (if selected SWRC uses fewer than 6 parameters)
@@ -28,6 +28,17 @@
# * param6 = L, tortuosity/connectivity parameter [-]
+# Table with two sets of six SWRC parameters
+# * first set (row): characteristics of fibric peat
+# * second set (row): characteristics of sapric peat
+# * source: Letts et al. 2000, doi:10.1080/07055900.2000.9649643
+# param1 param2 param3 param4 param5 param6
+ 1.03 0.93 2.7 2419.2 0.0000 0.0000
+ 1.01 0.83 12.0 0.864 0.0000 0.0000
+
+
+# Table with six SWRC parameters for up to `MAX_LAYERS` rows (soil layers)
+# * the soil layers must match `soils.in`
# param1 param2 param3 param4 param5 param6
18.6080 0.42703 5.3020 53.90697 0.0000 0.0000
20.4644 0.43290 7.0500 37.41493 0.0000 0.0000
diff --git a/tests/example/Input/swrc_params_FXW.in b/tests/example/Input/swrc_params_FXW.in
index 7f1bc8df3..139628b9f 100644
--- a/tests/example/Input/swrc_params_FXW.in
+++ b/tests/example/Input/swrc_params_FXW.in
@@ -1,7 +1,6 @@
#------ Input for Soil Water Retention Curves (by soil layer) ------
-# A table with up to `MAX_LAYERS` rows (soil layers) and 6 columns:
-# - the soil layers must match `soils.in`
+# Tables with 6 columns:
# - the interpretation of columns (SWRC parameters) depends on the
# selected SWRC (see `siteparam.in`)
# - unused columns are ignored (if selected SWRC uses fewer than 6 parameters)
@@ -15,6 +14,19 @@
# * param6 = L, tortuosity/connectivity parameter [-]
+# Table with two sets of six SWRC parameters
+# * first set (row): characteristics of fibric peat
+# * second set (row): characteristics of sapric peat
+# * source: Letts et al. 2000, doi:10.1080/07055900.2000.9649643
+# * Note: values for the van Genuchten SWRC assumed for
+# comparable parameters; however, param4 (m) and param6 (L) are missing
+# param1 param2 param3 param4 param5 param6
+ 0.93 0.08 1.9 NaN 2419.2 NaN
+ 0.83 0.003 1.6 NaN 0.864 NaN
+
+
+# Table with six SWRC parameters for up to `MAX_LAYERS` rows (soil layers)
+# * the soil layers must match `soils.in`
# param1 param2 param3 param4 param5 param6
0.437461 0.050757 1.247689 0.308681 22.985379 2.697338
0.452401 0.103033 1.146533 0.195394 89.365566 2.843288
diff --git a/tests/example/Input/swrc_params_vanGenuchten1980.in b/tests/example/Input/swrc_params_vanGenuchten1980.in
index d31396bcc..5a3da5472 100644
--- a/tests/example/Input/swrc_params_vanGenuchten1980.in
+++ b/tests/example/Input/swrc_params_vanGenuchten1980.in
@@ -1,7 +1,6 @@
#------ Input for Soil Water Retention Curves (by soil layer) ------
-# A table with up to `MAX_LAYERS` rows (soil layers) and 6 columns:
-# - the soil layers must match `soils.in`
+# Tables with 6 columns:
# - the interpretation of columns (SWRC parameters) depends on the
# selected SWRC (see `siteparam.in`)
# - unused columns are ignored (if selected SWRC uses fewer than 6 parameters)
@@ -13,6 +12,18 @@
# * param4 = n, measure of the pore-size distribution [-]
# * param5 = saturated hydraulic conductivity [cm/day]
+
+# Table with two sets of six SWRC parameters
+# * first set (row): characteristics of fibric peat
+# * second set (row): characteristics of sapric peat
+# * source: Letts et al. 2000, doi:10.1080/07055900.2000.9649643
+# param1 param2 param3 param4 param5 param6
+ 0.04 0.93 0.08 1.9 2419.2 0.0000
+ 0.22 0.83 0.003 1.6 0.864 0.0000
+
+
+# Table with six SWRC parameters for up to `MAX_LAYERS` rows (soil layers)
+# * the soil layers must match `soils.in`
# param1 param2 param3 param4 param5 param6
0.07564425 0.3925437 0.010035788 1.412233 19.871040 0
0.10061329 0.4011315 0.009425738 1.352274 9.090754 0
diff --git a/tests/gtests/sw_testhelpers.cc b/tests/gtests/sw_testhelpers.cc
index bbda90ef4..b28bb8467 100644
--- a/tests/gtests/sw_testhelpers.cc
+++ b/tests/gtests/sw_testhelpers.cc
@@ -92,6 +92,8 @@ void create_test_soillayers(
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
double soiltemp[MAX_LAYERS] = {-1, -1, -1, -1, 0, 0, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2};
+ double om[MAX_LAYERS] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
int const nRegions = 3;
double regionLowerBounds[3] = {20., 50., 100.};
@@ -112,6 +114,7 @@ void create_test_soillayers(
pclay,
imperm,
soiltemp,
+ om,
nRegions,
regionLowerBounds,
LogInfo
@@ -134,7 +137,7 @@ void setup_SW_Site_for_tests(SW_SITE *SW_Site) {
SW_Site->slow_drain_coeff = 0.02;
- SW_Site->site_has_swrcp = swFALSE;
+ SW_Site->site_has_swrcpMineralSoil = swFALSE;
(void) snprintf(
SW_Site->site_swrc_name,
sizeof SW_Site->site_swrc_name,
@@ -150,6 +153,18 @@ void setup_SW_Site_for_tests(SW_SITE *SW_Site) {
"Cosby1984AndOthers"
);
SW_Site->site_ptf_type = encode_str2ptf(SW_Site->site_ptf_name);
+
+ SW_Site->swrcpOM[0][0] = 1.03;
+ SW_Site->swrcpOM[1][0] = 1.01;
+
+ SW_Site->swrcpOM[0][1] = 0.93;
+ SW_Site->swrcpOM[1][1] = 0.83;
+
+ SW_Site->swrcpOM[0][2] = 2.7;
+ SW_Site->swrcpOM[1][2] = 12.0;
+
+ SW_Site->swrcpOM[0][3] = 2419.2;
+ SW_Site->swrcpOM[1][3] = 0.864;
}
/* Set up global variables for testing and read in values from SOILWAT2 example
diff --git a/tests/gtests/test_SW_Flow_Lib.cc b/tests/gtests/test_SW_Flow_Lib.cc
index 9b5abc189..c881ddaa2 100644
--- a/tests/gtests/test_SW_Flow_Lib.cc
+++ b/tests/gtests/test_SW_Flow_Lib.cc
@@ -165,17 +165,19 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
// declare inputs
double pptleft = 5.0;
- double standingWater;
+ double standingWater = 0.;
double drainout;
+ const double swc0init = 0.8;
// ***** Tests when nlyrs = 1 ***** //
/// provide inputs
int nlyrs = 1;
- double swc[1] = {0.8};
+ double swc[1] = {swc0init};
double swcfc[1] = {1.1};
double swcsat[1] = {1.6};
double impermeability[1] = {0.};
double drain[1] = {0.};
+ double ksat[1] = {1e6}; // very large number
infiltrate_water_high(
swc,
@@ -185,6 +187,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc,
swcsat,
+ ksat,
impermeability,
&standingWater,
lyrFrozen
@@ -201,7 +204,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
EXPECT_DOUBLE_EQ(drainout, drain[0]);
/* Test when pptleft and standingWater are 0 (No drainage) */
- pptleft = 0.0, standingWater = 0.0, drain[0] = 0., swc[0] = 0.8,
+ pptleft = 0.0, standingWater = 0.0, drain[0] = 0., swc[0] = swc0init,
swcfc[0] = 1.1, swcsat[0] = 1.6;
infiltrate_water_high(
@@ -212,6 +215,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc,
swcsat,
+ ksat,
impermeability,
&standingWater,
lyrFrozen
@@ -223,8 +227,9 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
/* Test when impermeability is greater than 0 and large precipitation */
pptleft = 20.0;
+ standingWater = 0.;
impermeability[0] = 1.;
- swc[0] = 0.8, drain[0] = 0.;
+ swc[0] = swc0init, drain[0] = 0.;
infiltrate_water_high(
swc,
@@ -234,6 +239,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc,
swcsat,
+ ksat,
impermeability,
&standingWater,
lyrFrozen
@@ -242,17 +248,49 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
// When impermeability is 1, drainage should be 0
EXPECT_DOUBLE_EQ(0., drain[0]);
- /* When impermeability is 1,standingWater should be equivalent to
- * pptLeft + swc[0] - swcsat[0]) */
- EXPECT_DOUBLE_EQ(standingWater, (pptleft + 0.8) - swcsat[0]);
+ /* When impermeability is 1, standingWater should be equivalent to
+ * pptLeft + swc0init - swcsat[0]) */
+ EXPECT_DOUBLE_EQ(standingWater, (pptleft + swc0init) - swcsat[0]);
+
+
+ /* Test when ksat is 0 and large precipitation */
+ ksat[0] = 0.;
+ standingWater = 0.;
+ pptleft = 20.0;
+ impermeability[0] = 0.;
+ swc[0] = swc0init, drain[0] = 0.;
+
+ infiltrate_water_high(
+ swc,
+ drain,
+ &drainout,
+ pptleft,
+ nlyrs,
+ swcfc,
+ swcsat,
+ ksat,
+ impermeability,
+ &standingWater,
+ lyrFrozen
+ );
+
+ // When ksat is 0, drainage should be 0
+ EXPECT_DOUBLE_EQ(0., drain[0]);
+
+ /* When ksat is 0, standingWater should be equivalent to
+ * pptLeft + swc0init - swcsat[0]) */
+ EXPECT_DOUBLE_EQ(standingWater, (pptleft + swc0init) - swcsat[0]);
+
// ***** Test when nlyrs = MAX_LAYERS (SW_Defines.h) ***** //
/// generate inputs using a for loop
unsigned int i;
nlyrs = MAX_LAYERS, pptleft = 5.0;
+ standingWater = 0.;
double *swc2 = new double[nlyrs];
double *swcfc2 = new double[nlyrs];
double *swcsat2 = new double[nlyrs];
+ double *ksat2 = new double[nlyrs];
double *impermeability2 = new double[nlyrs];
double *drain2 = new double[nlyrs];
@@ -264,6 +302,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
swcfc2[i] = RandNorm(1, .5, &infiltrate_rng);
// swcsat will always be greater than swcfc in each layer
swcsat2[i] = swcfc2[i] + .1;
+ ksat2[i] = 1e6; // very large number
impermeability2[i] = 0.0;
}
@@ -275,6 +314,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc2,
swcsat2,
+ ksat2,
impermeability2,
&standingWater,
lyrFrozen
@@ -316,6 +356,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc3,
swcsat3,
+ ksat2,
impermeability2,
&standingWater,
lyrFrozen
@@ -333,6 +374,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
double *swcfc4 = new double[nlyrs];
double *swcsat4 = new double[nlyrs];
pptleft = 20.0;
+ standingWater = 0.;
for (i = 0; i < MAX_LAYERS; i++) {
swc4[i] = RandNorm(1., 0.5, &infiltrate_rng);
@@ -344,7 +386,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
}
// Need to hard code this value because swc4 is altered by function
- swc4[0] = 0.8;
+ swc4[0] = swc0init;
infiltrate_water_high(
swc4,
@@ -354,14 +396,15 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc4,
swcsat4,
+ ksat2,
impermeability4,
&standingWater,
lyrFrozen
);
- /* When impermeability is 1,standingWater should be equivalent to pptLeft +
- * swc[0] - swcsat[0]) */
- EXPECT_DOUBLE_EQ(standingWater, (pptleft + 0.8) - swcsat4[0]);
+ /* When impermeability is 1, standingWater should be equivalent to
+ pptLeft + swc0init - swcsat[0]) */
+ EXPECT_DOUBLE_EQ(standingWater, (pptleft + swc0init) - swcsat4[0]);
for (i = 0; i < MAX_LAYERS; i++) {
// When impermeability is 1, drainage should be 0
@@ -375,6 +418,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
double *swcfc5 = new double[nlyrs];
double *swcsat5 = new double[nlyrs];
pptleft = 5.0;
+ standingWater = 0.;
for (i = 0; i < MAX_LAYERS; i++) {
swc5[i] = RandNorm(1.2, 0.5, &infiltrate_rng);
@@ -394,6 +438,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
nlyrs,
swcfc5,
swcsat5,
+ ksat2,
impermeability5,
&standingWater,
lyrFrozen
@@ -410,27 +455,10 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) {
// deallocate pointers
- double *array_list[] = {
- impermeability2,
- drain2,
- swc2,
- swcfc2,
- swcsat2,
- drain3,
- swc3,
- swcfc3,
- swcsat3,
- impermeability4,
- drain4,
- swc4,
- swcfc4,
- swcsat4,
- impermeability5,
- drain5,
- swc5,
- swcfc5,
- swcsat5
- };
+ double *array_list[] = {drain2, swc2, swcfc2, swcsat2, impermeability2,
+ drain3, swc3, swcfc3, swcsat3, impermeability4,
+ drain4, swc4, swcfc4, swcsat4, impermeability5,
+ drain5, swc5, swcfc5, swcsat5, ksat2};
for (i = 0; i < sw_length(array_list); i++) {
delete[] array_list[i];
diff --git a/tests/gtests/test_SW_Site.cc b/tests/gtests/test_SW_Site.cc
index c42b86a3b..159d2c502 100644
--- a/tests/gtests/test_SW_Site.cc
+++ b/tests/gtests/test_SW_Site.cc
@@ -371,6 +371,87 @@ TEST(SiteTest, SiteSWRCpChecks) {
swrcp[5] = tmp;
}
+// Test 'SWRC_bulkSoilParameters'
+TEST(SiteTest, SWRCBulkSoilParameters) {
+ double swrcp[SWRC_PARAM_NMAX];
+ double swrcpMin[SWRC_PARAM_NMAX];
+ double swrcpOrg[2][SWRC_PARAM_NMAX];
+ double fom;
+ const double depthSapric = 50.;
+ double depthT = 0.;
+ double depthB = 10.;
+
+ unsigned int k;
+ const unsigned int swrc_type = 1;
+
+ // Initialize swrcps
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ swrcpMin[k] = 1.;
+ swrcpOrg[0][k] = 10.;
+ swrcpOrg[1][k] = 20.;
+ }
+
+ // Expect swrcp = mineral if organic matter is 0
+ fom = 0.;
+ SWRC_bulkSoilParameters(
+ swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB
+ );
+
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ EXPECT_DOUBLE_EQ(swrcp[k], swrcpMin[k]);
+ }
+
+ // Expect swrcp = fibric if organic matter is 1 and layer at surface
+ fom = 1.;
+ depthT = 0.;
+ depthB = 0.;
+ SWRC_bulkSoilParameters(
+ swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB
+ );
+
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ EXPECT_DOUBLE_EQ(swrcp[k], swrcpOrg[0][k]);
+ }
+
+ // Expect fibric < swrcp < sapric if organic matter is 1 and layer medium
+ fom = 1.;
+ depthT = depthSapric / 4.;
+ depthB = depthT + depthSapric / 4.;
+ SWRC_bulkSoilParameters(
+ swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB
+ );
+
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ EXPECT_GT(swrcp[k], swrcpOrg[0][k]);
+ EXPECT_LT(swrcp[k], swrcpOrg[1][k]);
+ }
+
+ // Expect swrcp = sapric if organic matter is 1 and layer is at depth
+ fom = 1.;
+ depthT = depthSapric;
+ depthB = depthT + 10.;
+ SWRC_bulkSoilParameters(
+ swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB
+ );
+
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ EXPECT_DOUBLE_EQ(swrcp[k], swrcpOrg[1][k]);
+ }
+
+ // Expect min < swrcp < fibric if organic matter is 0-1 and layer at surface
+ fom = 0.5;
+ depthT = 0.;
+ depthB = 0.;
+ SWRC_bulkSoilParameters(
+ swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB
+ );
+
+ for (k = 0; k < SWRC_PARAM_NMAX; k++) {
+ EXPECT_GT(swrcp[k], swrcpMin[k]);
+ EXPECT_LT(swrcp[k], swrcpOrg[0][k]);
+ }
+}
+
// Test 'PTF_RawlsBrakensiek1985'
TEST(SiteTest, SitePTFRawlsBrakensiek1985) {
LOG_INFO LogInfo;
@@ -381,6 +462,7 @@ TEST(SiteTest, SitePTFRawlsBrakensiek1985) {
double theta_min;
double clay = 0.1;
double sand = 0.6;
+ const double fom = 0.;
double porosity = 0.4;
int k1;
int k2;
@@ -388,27 +470,27 @@ TEST(SiteTest, SitePTFRawlsBrakensiek1985) {
//--- EXPECT SW_MISSING if soil texture is out of range
// within range: sand [0.05, 0.7], clay [0.05, 0.6], porosity [0.1, 1[
- PTF_RawlsBrakensiek1985(&theta_min, 0., clay, porosity, &LogInfo);
+ PTF_RawlsBrakensiek1985(&theta_min, 0., clay, fom, porosity, &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, SW_MISSING);
- PTF_RawlsBrakensiek1985(&theta_min, 0.75, clay, porosity, &LogInfo);
+ PTF_RawlsBrakensiek1985(&theta_min, 0.75, clay, fom, porosity, &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, SW_MISSING);
- PTF_RawlsBrakensiek1985(&theta_min, sand, 0., porosity, &LogInfo);
+ PTF_RawlsBrakensiek1985(&theta_min, sand, 0., fom, porosity, &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, SW_MISSING);
- PTF_RawlsBrakensiek1985(&theta_min, sand, 0.65, porosity, &LogInfo);
+ PTF_RawlsBrakensiek1985(&theta_min, sand, 0.65, fom, porosity, &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, SW_MISSING);
- PTF_RawlsBrakensiek1985(&theta_min, sand, clay, 0., &LogInfo);
+ PTF_RawlsBrakensiek1985(&theta_min, sand, clay, fom, 0., &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, SW_MISSING);
- PTF_RawlsBrakensiek1985(&theta_min, sand, clay, 1., &LogInfo);
+ PTF_RawlsBrakensiek1985(&theta_min, sand, clay, fom, 1., &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, SW_MISSING);
@@ -424,7 +506,7 @@ TEST(SiteTest, SitePTFRawlsBrakensiek1985) {
porosity = 0.1 + (double) k3 / 5. * (0.99 - 0.1);
PTF_RawlsBrakensiek1985(
- &theta_min, sand, clay, porosity, &LogInfo
+ &theta_min, sand, clay, fom, porosity, &LogInfo
);
// exit test program if unexpected error
sw_fail_on_error(&LogInfo);
@@ -435,8 +517,8 @@ TEST(SiteTest, SitePTFRawlsBrakensiek1985) {
}
}
- // Expect theta_min = 0 if sand = 0.4, clay = 0.5, and porosity = 0.1
- PTF_RawlsBrakensiek1985(&theta_min, 0.4, 0.5, 0.1, &LogInfo);
+ // Expect theta_min = 0 if sand = 0.4, clay = 0.5, fom = 0., porosity = 0.1
+ PTF_RawlsBrakensiek1985(&theta_min, 0.4, 0.5, 0.0, 0.1, &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
EXPECT_DOUBLE_EQ(theta_min, 0);
}
diff --git a/tests/gtests/test_WaterBalance.cc b/tests/gtests/test_WaterBalance.cc
index 8e4395599..8ec7edf3f 100644
--- a/tests/gtests/test_WaterBalance.cc
+++ b/tests/gtests/test_WaterBalance.cc
@@ -192,6 +192,7 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithHighGravelVolume) {
// Re-calculate soils
SW_SIT_init_run(&SW_Run.VegProd, &SW_Run.Site, &LogInfo);
+ SW_SWC_init_run(&SW_Run.SoilWat, &SW_Run.Site, &SW_Run.Weather.temp_snow);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
// Run the simulation
@@ -272,6 +273,46 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithVegetationFromClimate1) {
}
}
+TEST_F(WaterBalanceFixtureTest, WaterBalanceWithOrganicMatter) {
+ unsigned int i;
+
+ // Set PTF (Cosby1984AndOthers handles OM only up to 8%)
+ (void) snprintf(
+ SW_Run.Site.site_ptf_name,
+ sizeof SW_Run.Site.site_ptf_name,
+ "%s",
+ "Cosby1984"
+ );
+ SW_Run.Site.site_ptf_type = encode_str2ptf(SW_Run.Site.site_ptf_name);
+ SW_Run.Site.site_has_swrcpMineralSoil = swTRUE;
+
+ // Set organic matter > 0
+ SW_Run.Site.fractionWeight_om[0] = 1.;
+ for (i = 1; i < SW_Run.Site.n_layers; i++) {
+ SW_Run.Site.fractionWeight_om[i] = 0.5;
+ }
+
+ // Update soils
+ SW_SIT_init_run(&SW_Run.VegProd, &SW_Run.Site, &LogInfo);
+ SW_SWC_init_run(&SW_Run.SoilWat, &SW_Run.Site, &SW_Run.Weather.temp_snow);
+ sw_fail_on_error(&LogInfo); // exit test program if unexpected error
+
+ // Two simulation years are sufficient
+ SW_Run.Model.startyr = 1980;
+ SW_Run.Model.endyr = 1981;
+
+ // Run the simulation
+ SW_CTL_main(&SW_Run, &SW_Domain.OutDom, &LogInfo);
+ sw_fail_on_error(&LogInfo); // exit test program if unexpected error
+
+ // Collect and output from daily checks
+ for (i = 0; i < N_WBCHECKS; i++) {
+ EXPECT_EQ(0, SW_Run.SoilWat.wbError[i])
+ << "Water balance error in test " << i << ": "
+ << SW_Run.SoilWat.wbErrorNames[i];
+ }
+}
+
TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCvanGenuchten1980) {
int i;
@@ -292,7 +333,7 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCvanGenuchten1980) {
"Rosetta3"
);
SW_Run.Site.site_ptf_type = encode_str2ptf(SW_Run.Site.site_ptf_name);
- SW_Run.Site.site_has_swrcp = swTRUE;
+ SW_Run.Site.site_has_swrcpMineralSoil = swTRUE;
free(SW_Domain.PathInfo.InFiles[eSWRCp]);
SW_Domain.PathInfo.InFiles[eSWRCp] =
@@ -305,6 +346,7 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCvanGenuchten1980) {
// Update soils
SW_SIT_init_run(&SW_Run.VegProd, &SW_Run.Site, &LogInfo);
+ SW_SWC_init_run(&SW_Run.SoilWat, &SW_Run.Site, &SW_Run.Weather.temp_snow);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
// Run the simulation
@@ -320,7 +362,7 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCvanGenuchten1980) {
}
TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCFXW) {
- int i;
+ unsigned int i;
// Set SWRC and PTF (and SWRC parameter input filename)
(void) snprintf(
@@ -339,7 +381,7 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCFXW) {
"neuroFX2021"
);
SW_Run.Site.site_ptf_type = encode_str2ptf(SW_Run.Site.site_ptf_name);
- SW_Run.Site.site_has_swrcp = swTRUE;
+ SW_Run.Site.site_has_swrcpMineralSoil = swTRUE;
free(SW_Domain.PathInfo.InFiles[eSWRCp]);
SW_Domain.PathInfo.InFiles[eSWRCp] =
@@ -350,8 +392,16 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCFXW) {
SW_SWRC_read(&SW_Run.Site, SW_Domain.PathInfo.InFiles, &LogInfo);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
+ // FXW doesn't yet handle organic matter:
+ // not all values for organic SWRC parameters have been determined
+ // (see "tests/example/Input/swrc_params_FXW.in")
+ for (i = 0; i < SW_Run.Site.n_layers; i++) {
+ SW_Run.Site.fractionWeight_om[i] = 0.;
+ }
+
// Update soils
SW_SIT_init_run(&SW_Run.VegProd, &SW_Run.Site, &LogInfo);
+ SW_SWC_init_run(&SW_Run.SoilWat, &SW_Run.Site, &SW_Run.Weather.temp_snow);
sw_fail_on_error(&LogInfo); // exit test program if unexpected error
// Run the simulation