From d2f23f35e4d18c19f473d6e0c9769fbb29895732 Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Wed, 25 Sep 2024 16:54:19 -0400 Subject: [PATCH 1/9] Fix WaterBalance tests - SWC-related variables need to be (re-)initialized after soils are updated --- tests/gtests/test_WaterBalance.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/gtests/test_WaterBalance.cc b/tests/gtests/test_WaterBalance.cc index 8e4395599..bd8a8b85d 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 @@ -305,6 +306,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 @@ -352,6 +354,7 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithSWRCFXW) { // 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 From e3d8a38c99c7b5c205627dcfd093c44c12e7b9ea Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Thu, 26 Sep 2024 16:22:38 -0400 Subject: [PATCH 2/9] New tests for organic SWRC parameters - tests for SWRCBulkSoilParameters() - water balance tests with soil organic matter > 0 --- tests/gtests/test_SW_Site.cc | 81 +++++++++++++++++++++++++++++++ tests/gtests/test_WaterBalance.cc | 30 ++++++++++++ 2 files changed, 111 insertions(+) diff --git a/tests/gtests/test_SW_Site.cc b/tests/gtests/test_SW_Site.cc index c42b86a3b..8b708d72a 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; + double depthSapric = 50.; + double depthT = 0.; + double depthB = 10.; + + unsigned int k; + + // 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( + 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( + 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( + 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( + 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( + 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; diff --git a/tests/gtests/test_WaterBalance.cc b/tests/gtests/test_WaterBalance.cc index bd8a8b85d..db61d8e4e 100644 --- a/tests/gtests/test_WaterBalance.cc +++ b/tests/gtests/test_WaterBalance.cc @@ -273,6 +273,36 @@ TEST_F(WaterBalanceFixtureTest, WaterBalanceWithVegetationFromClimate1) { } } +TEST_F(WaterBalanceFixtureTest, WaterBalanceWithOrganicMatter) { + unsigned int i; + + // 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; From 7d548d5ec6e30222e40b45930d1790709c61977f Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Thu, 26 Sep 2024 15:35:04 -0400 Subject: [PATCH 3/9] FXW doesn't yet handle organic matter: undetermined organic SWRC paramaters - for now, set om to 0 for water balance test --- tests/gtests/test_WaterBalance.cc | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tests/gtests/test_WaterBalance.cc b/tests/gtests/test_WaterBalance.cc index db61d8e4e..7976ed606 100644 --- a/tests/gtests/test_WaterBalance.cc +++ b/tests/gtests/test_WaterBalance.cc @@ -352,7 +352,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( @@ -382,6 +382,13 @@ 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); From 8f895992f208793d33d92746c3eca9382febec9f Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Wed, 25 Sep 2024 16:52:58 -0400 Subject: [PATCH 4/9] Soil organic matter influence on soil water retention Close #397 The implemented approach has two steps: 1) Determine 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 (with a default of 50 cm). 2) Estimate bulk soil parameters of the soil water retention curve 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. New inputs * siteparam.in: depthSapric * soils.in: fractionWeight_om * swrc_params*.in: swrcpOM[fibric], swrcpOM[sapric] New functions * interpolateFibricSapric() to linearly interpolate organic matter parameter between fibric and sapric characteristics * SWRC_bulkSoilParameters() to calculate parameters of the bulk soil are as the weighted average of parameters of the organic and mineral components Functions with new arguments/behavior * SW_SIT_read() now reads in value for depthSapric from "siteparam.in" * SW_LYR_read() now reads in values for fractionWeight_om from "soils.in" * SW_SWRC_read() now reads in values for swrcpOM[fibric] and swrcpOM[sapric] from "swrc_params*.in" * SW_SIT_init_run() now runs checks for swrcpOM[fibric] and swrcpOM[sapric] and calculates bulk soil SWRCp from organic and mineral soil components via a call to SWRC_bulkSoilParameters() * SW_swcBulk_saturated() gained argument `"fom"` which is passed on to PTF_Saxton2006() if legacy `"Cosby1984AndOthers"` is selected as PTF * SW_swcBulk_minimum() and ui_theta_min() gained argument `"fom"` which is passed on to PTF_RawlsBrakensiek1985() if legacy `"Cosby1984AndOthers"` is selected as PTF * PTF_Saxton2006() gained argument `"fom"` and throws an error if fom > 0.08 (these published equations incorporate influence of organic matter up to 8% which was previously fixed at 0%) * PTF_RawlsBrakensiek1985() gained argument `"fom"` and throws an error if fom > 0 * SW_check_soil_properties() now checks that user-input `fractionWeight_om` is within 0-1 * set_soillayers() gained argument `"pom"` to set a value for `fractionWeight_om` --- NEWS.md | 20 + include/SW_Site.h | 19 +- include/SW_datastructs.h | 35 +- src/SW_Site.c | 502 +++++++++++++----- src/SW_SoilWater.c | 5 +- tests/example/Input/siteparam.in | 7 +- tests/example/Input/soils.in | 19 +- tests/example/Input/swrc_params.in | 17 +- tests/example/Input/swrc_params_FXW.in | 16 +- .../Input/swrc_params_vanGenuchten1980.in | 15 +- tests/gtests/sw_testhelpers.cc | 17 +- tests/gtests/test_SW_Site.cc | 20 +- tests/gtests/test_WaterBalance.cc | 14 +- 13 files changed, 526 insertions(+), 180 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5791bf5b4..a80a99c26 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,26 @@ * Tests now require `c++17` and utilize `googletest` `v1.15.2` (issue #427). +* SOILWAT2 can now represent the influence of soil organic matter on + soil water retention (#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. + +## 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_Site.h b/include/SW_Site.h index f97183486..2cb308e70 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,34 @@ 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( + double *swrcp, + const double *swrcpMS, + const double swrcOM[2][SWRC_PARAM_NMAX], + double fom, + double depthSapric, + double depthT, + double depthB +); double calculate_soilBulkDensity(double matricDensity, double fractionGravel); @@ -322,6 +334,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..2c661680a 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 @@ -368,6 +374,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 +386,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_Site.c b/src/SW_Site.c index 6c7507987..554391f1a 100644 --- a/src/SW_Site.c +++ b/src/SW_Site.c @@ -214,6 +214,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 +289,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 +298,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 +313,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 +340,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)) { @@ -418,7 +430,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 +443,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 +452,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 +465,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 +489,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 +531,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 +549,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 +562,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 +571,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 +612,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 +626,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 +642,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 +682,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 +1040,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 +1102,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 +1157,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 +1167,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 +1189,126 @@ 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. + +@param[out] *swrcp Vector of SWRC parameters of the bulk soil +@param[in] *swrcpMS Vector of SWRC parameters of the mineral component +@param[in] swrcOM 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( + double *swrcp, + const double *swrcpMS, + const double swrcOM[2][SWRC_PARAM_NMAX], + double fom, + double depthSapric, + double depthT, + double depthB +) { + int k; + double pOM; + + if (GT(fom, 0.)) { + /* Has organic matter: + interpolate between organic and mineral components */ + for (k = 0; k < SWRC_PARAM_NMAX; k++) { + pOM = interpolateFibricSapric( + swrcOM[0][k], swrcOM[1][k], depthSapric, depthT, depthB + ); + /* Note: this approach does not account for effects of + connected flow pathways if fom > threshold for + saturated hydraulic conductivity */ + swrcp[k] = fom * pOM + (1 - fom) * swrcpMS[k]; + } + + } else { + /* No organic matter: use mineral components */ + for (k = 0; k < SWRC_PARAM_NMAX; k++) { + swrcp[k] = swrcpMS[k]; + } } } @@ -1348,13 +1497,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 +1664,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 +1689,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 +1705,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 +1775,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 +1843,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 +1858,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 +1876,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 +1888,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 +1929,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 +1976,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 +2011,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 +2061,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 +2163,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 +2194,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 +2210,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 +2226,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 +2235,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 +2250,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 +2332,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 +2347,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 +2442,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,14 +2463,41 @@ 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->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 ); @@ -2305,6 +2546,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 +2562,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 +2591,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 +2609,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 +2622,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 +2647,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 +2715,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 +2760,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 +2769,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 +2803,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 +2832,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 +2843,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 +2856,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 +2869,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_Site.cc b/tests/gtests/test_SW_Site.cc index 8b708d72a..db033dab9 100644 --- a/tests/gtests/test_SW_Site.cc +++ b/tests/gtests/test_SW_Site.cc @@ -449,7 +449,6 @@ TEST(SiteTest, SWRCBulkSoilParameters) { EXPECT_GT(swrcp[k], swrcpMin[k]); EXPECT_LT(swrcp[k], swrcpOrg[0][k]); } - } // Test 'PTF_RawlsBrakensiek1985' @@ -462,6 +461,7 @@ TEST(SiteTest, SitePTFRawlsBrakensiek1985) { double theta_min; double clay = 0.1; double sand = 0.6; + double fom = 0.; double porosity = 0.4; int k1; int k2; @@ -469,27 +469,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); @@ -505,7 +505,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); @@ -516,8 +516,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 7976ed606..8ec7edf3f 100644 --- a/tests/gtests/test_WaterBalance.cc +++ b/tests/gtests/test_WaterBalance.cc @@ -276,6 +276,16 @@ 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++) { @@ -323,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] = @@ -371,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] = From 8029a9d1a70e53405d509ff22fb9608613c38792 Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Fri, 11 Oct 2024 10:10:54 -0400 Subject: [PATCH 5/9] Limit saturated percolation ** This changes simulation output ** * Different outcomes during periods of high infiltration (e.g., snowmelt) and during conditions of low hydraulic conductivity (e.g., frozen soils, sapric organic matter). * Saturated percolation is now limited -- previously, it was unbounded. The upper bound is a function based on saturated hydraulic conductivity (which includes effects of organic matter), frozen soils, and a user-specified `"permeability"` factor. * Bulk soil saturated hydraulic conductivity parameter now accounts for flow pathways through organic matter above a threshold and assumes conductivities through mineral and organic components in series outside of those pathways. --- NEWS.md | 37 ++++++++++------ include/SW_Flow_lib.h | 1 + include/SW_Site.h | 2 + include/SW_datastructs.h | 6 ++- src/SW_Flow.c | 2 + src/SW_Flow_lib.c | 19 ++++---- src/SW_Site.c | 74 +++++++++++++++++++++++++++++--- tests/gtests/test_SW_Flow_Lib.cc | 63 ++++++++++++++++++++++----- tests/gtests/test_SW_Site.cc | 11 ++--- 9 files changed, 172 insertions(+), 43 deletions(-) diff --git a/NEWS.md b/NEWS.md index a80a99c26..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,15 +15,24 @@ * Tests now require `c++17` and utilize `googletest` `v1.15.2` (issue #427). -* SOILWAT2 can now represent the influence of soil organic matter on - soil water retention (#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. +* 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 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 2cb308e70..7a9f84bd2 100644 --- a/include/SW_Site.h +++ b/include/SW_Site.h @@ -259,6 +259,7 @@ double SW_swcBulk_minimum( LOG_INFO *LogInfo ); + void PTF_Saxton2006( double *theta_sat, double sand, double clay, double fom, LOG_INFO *LogInfo ); @@ -273,6 +274,7 @@ void PTF_RawlsBrakensiek1985( ); void SWRC_bulkSoilParameters( + unsigned int swrc_type, double *swrcp, const double *swrcpMS, const double swrcOM[2][SWRC_PARAM_NMAX], diff --git a/include/SW_datastructs.h b/include/SW_datastructs.h index 2c661680a..0f9260d7c 100644 --- a/include/SW_datastructs.h +++ b/include/SW_datastructs.h @@ -366,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] */ 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 554391f1a..0b745a345 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: @@ -370,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 */ /* --------------------------------------------------- */ @@ -1268,6 +1289,9 @@ 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[out] *swrcp Vector of SWRC parameters of the bulk soil @param[in] *swrcpMS Vector of SWRC parameters of the mineral component @param[in] swrcOM Array of length two of vectors of SWRC parameters @@ -1280,6 +1304,7 @@ fibric peat characteristics assumed at surface and sapric peat at depth. @param[in] depthB Depth at bottom of soil layer [`cm`] */ void SWRC_bulkSoilParameters( + unsigned int swrc_type, double *swrcp, const double *swrcpMS, const double swrcOM[2][SWRC_PARAM_NMAX], @@ -1288,24 +1313,54 @@ void SWRC_bulkSoilParameters( double depthT, double depthB ) { - int k; + 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 (GT(fom, 0.)) { + 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( swrcOM[0][k], swrcOM[1][k], depthSapric, depthT, depthB ); - /* Note: this approach does not account for effects of - connected flow pathways if fom > threshold for - saturated hydraulic conductivity */ - swrcp[k] = fom * pOM + (1 - fom) * swrcpMS[k]; + + 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: use mineral components */ + /* No organic matter: bulk soil corresponds to mineral components */ for (k = 0; k < SWRC_PARAM_NMAX; k++) { swrcp[k] = swrcpMS[k]; } @@ -2480,6 +2535,7 @@ void SW_SIT_init_run( /* 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, @@ -2504,6 +2560,10 @@ void SW_SIT_init_run( 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); diff --git a/tests/gtests/test_SW_Flow_Lib.cc b/tests/gtests/test_SW_Flow_Lib.cc index 9b5abc189..d471f4983 100644 --- a/tests/gtests/test_SW_Flow_Lib.cc +++ b/tests/gtests/test_SW_Flow_Lib.cc @@ -167,15 +167,17 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { double pptleft = 5.0; double standingWater; double drainout; + double swc0init = 0.8; // ***** Tests when nlyrs = 1 ***** // /// provide inputs int nlyrs = 1; - double swc[1] = {0.8}; + double swc[1]; 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,9 +248,39 @@ 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 @@ -253,6 +289,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { 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 +301,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 +313,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { nlyrs, swcfc2, swcsat2, + ksat2, impermeability2, &standingWater, lyrFrozen @@ -316,6 +355,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { nlyrs, swcfc3, swcsat3, + ksat2, impermeability2, &standingWater, lyrFrozen @@ -344,7 +384,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 +394,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 @@ -394,6 +435,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { nlyrs, swcfc5, swcsat5, + ksat2, impermeability5, &standingWater, lyrFrozen @@ -411,6 +453,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { // deallocate pointers double *array_list[] = { + ksat2, impermeability2, drain2, swc2, diff --git a/tests/gtests/test_SW_Site.cc b/tests/gtests/test_SW_Site.cc index db033dab9..dee63575e 100644 --- a/tests/gtests/test_SW_Site.cc +++ b/tests/gtests/test_SW_Site.cc @@ -382,6 +382,7 @@ TEST(SiteTest, SWRCBulkSoilParameters) { double depthB = 10.; unsigned int k; + unsigned int swrc_type = 1; // Initialize swrcps for (k = 0; k < SWRC_PARAM_NMAX; k++) { @@ -393,7 +394,7 @@ TEST(SiteTest, SWRCBulkSoilParameters) { // Expect swrcp = mineral if organic matter is 0 fom = 0.; SWRC_bulkSoilParameters( - swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB + swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB ); for (k = 0; k < SWRC_PARAM_NMAX; k++) { @@ -405,7 +406,7 @@ TEST(SiteTest, SWRCBulkSoilParameters) { depthT = 0.; depthB = 0.; SWRC_bulkSoilParameters( - swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB + swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB ); for (k = 0; k < SWRC_PARAM_NMAX; k++) { @@ -417,7 +418,7 @@ TEST(SiteTest, SWRCBulkSoilParameters) { depthT = depthSapric / 4.; depthB = depthT + depthSapric / 4.; SWRC_bulkSoilParameters( - swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB + swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB ); for (k = 0; k < SWRC_PARAM_NMAX; k++) { @@ -430,7 +431,7 @@ TEST(SiteTest, SWRCBulkSoilParameters) { depthT = depthSapric; depthB = depthT + 10.; SWRC_bulkSoilParameters( - swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB + swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB ); for (k = 0; k < SWRC_PARAM_NMAX; k++) { @@ -442,7 +443,7 @@ TEST(SiteTest, SWRCBulkSoilParameters) { depthT = 0.; depthB = 0.; SWRC_bulkSoilParameters( - swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB + swrc_type, swrcp, swrcpMin, swrcpOrg, fom, depthSapric, depthT, depthB ); for (k = 0; k < SWRC_PARAM_NMAX; k++) { From a2034abf6c860f440fb1af89c5592353854b06ae Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Mon, 21 Oct 2024 17:10:27 -0400 Subject: [PATCH 6/9] Fix soc commits - SWRC_bulkSoilParameters(): added documentation of argument "swrc_type" to fix documentation - SWRC_bulkSoilParameters(): corrected spelling of argument "swrcOM" to "swrcpOM" - tests "SWFlowSaturatedPercolation", "SWRCBulkSoilParameters" and "SitePTFRawlsBrakensiek1985": clang-tidy: use const variables where appropriate - test "SWFlowSaturatedPercolation": clang-format --- include/SW_Site.h | 2 +- src/SW_Site.c | 7 ++++--- tests/gtests/test_SW_Flow_Lib.cc | 28 +++++----------------------- tests/gtests/test_SW_Site.cc | 6 +++--- 4 files changed, 13 insertions(+), 30 deletions(-) diff --git a/include/SW_Site.h b/include/SW_Site.h index 7a9f84bd2..b0fe521d1 100644 --- a/include/SW_Site.h +++ b/include/SW_Site.h @@ -277,7 +277,7 @@ void SWRC_bulkSoilParameters( unsigned int swrc_type, double *swrcp, const double *swrcpMS, - const double swrcOM[2][SWRC_PARAM_NMAX], + const double swrcpOM[2][SWRC_PARAM_NMAX], double fom, double depthSapric, double depthT, diff --git a/src/SW_Site.c b/src/SW_Site.c index 0b745a345..567472a3a 100644 --- a/src/SW_Site.c +++ b/src/SW_Site.c @@ -1292,9 +1292,10 @@ 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] swrcOM Array of length two of vectors of SWRC parameters +@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 @@ -1307,7 +1308,7 @@ void SWRC_bulkSoilParameters( unsigned int swrc_type, double *swrcp, const double *swrcpMS, - const double swrcOM[2][SWRC_PARAM_NMAX], + const double swrcpOM[2][SWRC_PARAM_NMAX], double fom, double depthSapric, double depthT, @@ -1329,7 +1330,7 @@ void SWRC_bulkSoilParameters( for (k = 0; k < SWRC_PARAM_NMAX; k++) { /* Interpolate organic parameter from surface to depth conditions */ pOM = interpolateFibricSapric( - swrcOM[0][k], swrcOM[1][k], depthSapric, depthT, depthB + swrcpOM[0][k], swrcpOM[1][k], depthSapric, depthT, depthB ); if (k == iksat) { diff --git a/tests/gtests/test_SW_Flow_Lib.cc b/tests/gtests/test_SW_Flow_Lib.cc index d471f4983..1ba17039e 100644 --- a/tests/gtests/test_SW_Flow_Lib.cc +++ b/tests/gtests/test_SW_Flow_Lib.cc @@ -167,7 +167,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { double pptleft = 5.0; double standingWater; double drainout; - double swc0init = 0.8; + const double swc0init = 0.8; // ***** Tests when nlyrs = 1 ***** // /// provide inputs @@ -452,28 +452,10 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { // deallocate pointers - double *array_list[] = { - ksat2, - 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 dee63575e..159d2c502 100644 --- a/tests/gtests/test_SW_Site.cc +++ b/tests/gtests/test_SW_Site.cc @@ -377,12 +377,12 @@ TEST(SiteTest, SWRCBulkSoilParameters) { double swrcpMin[SWRC_PARAM_NMAX]; double swrcpOrg[2][SWRC_PARAM_NMAX]; double fom; - double depthSapric = 50.; + const double depthSapric = 50.; double depthT = 0.; double depthB = 10.; unsigned int k; - unsigned int swrc_type = 1; + const unsigned int swrc_type = 1; // Initialize swrcps for (k = 0; k < SWRC_PARAM_NMAX; k++) { @@ -462,7 +462,7 @@ TEST(SiteTest, SitePTFRawlsBrakensiek1985) { double theta_min; double clay = 0.1; double sand = 0.6; - double fom = 0.; + const double fom = 0.; double porosity = 0.4; int k1; int k2; From 3916d7d99519f29bf23bed1cfab4413f85edf9ed Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Wed, 23 Oct 2024 07:26:28 -0400 Subject: [PATCH 7/9] Fix SWFlowTest.SWFlowSaturatedPercolation - the variable standingWater was previously used uninitialized which lead to NaN when compiled with g++ --- tests/gtests/test_SW_Flow_Lib.cc | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/tests/gtests/test_SW_Flow_Lib.cc b/tests/gtests/test_SW_Flow_Lib.cc index 1ba17039e..ddb114804 100644 --- a/tests/gtests/test_SW_Flow_Lib.cc +++ b/tests/gtests/test_SW_Flow_Lib.cc @@ -165,7 +165,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { // declare inputs double pptleft = 5.0; - double standingWater; + double standingWater = 0.; double drainout; const double swc0init = 0.8; @@ -286,6 +286,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { /// 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]; @@ -373,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); @@ -400,7 +402,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { lyrFrozen ); - /* When impermeability is 1,standingWater should be equivalent to + /* When impermeability is 1, standingWater should be equivalent to pptLeft + swc0init - swcsat[0]) */ EXPECT_DOUBLE_EQ(standingWater, (pptleft + swc0init) - swcsat4[0]); @@ -416,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); From 846f4074908906547cb8139c98682d1ece540314 Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Wed, 23 Oct 2024 11:05:24 -0400 Subject: [PATCH 8/9] Fix swrcpOM 2-dim array function argument Previously, SWRC_bulkSoilParameters() was called with "SW_Site->swrcpOM" as argument for "const double swrcpOM[2][SWRC_PARAM_NMAX]" to pass fibric and sapric organic matter SWRC parameters. Problems (here using gcc/g++ v14.2.0): - gcc warning (`CC=gcc make clean all`): "'SWRC_bulkSoilParameters' accessing 96 bytes in a region of size 48 [-Wstringop-overflow=]" - gcc error (`CC=gcc make clean bin_debug_severe`): "invalid use of pointers to arrays with different qualifiers in ISO C before C23 [-Werror=pedantic]" - g++ error (`CXX=g++ make clean test`): "cannot convert 'const double (*)[]' to 'const double (*)[6]'" Option 1: coerce call with "(const double (*)[SWRC_PARAM_NMAX]) SW_Site->swrcpOM" as argument for "const double swrcpOM[][SWRC_PARAM_NMAX]" - gcc warning (`CC=gcc make clean all`) resolved - gcc error (`CC=gcc make clean bin_debug_severe`) resolved - g++ error (`CXX=g++ make clean test`) resolved Option 2: remove "const", i.e., call with "SW_Site->swrcpOM" as argument for "double swrcpOM[][SWRC_PARAM_NMAX]" - gcc warning (`CC=gcc make clean all`) resolved - gcc error (`CC=gcc make clean bin_debug_severe`) resolved - g++ error (`CXX=g++ make clean test`) resolved This commit implements option 2. Thanks to @N1ckP3rsl3y --- include/SW_Site.h | 2 +- src/SW_Site.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/include/SW_Site.h b/include/SW_Site.h index b0fe521d1..21e01e8f5 100644 --- a/include/SW_Site.h +++ b/include/SW_Site.h @@ -277,7 +277,7 @@ void SWRC_bulkSoilParameters( unsigned int swrc_type, double *swrcp, const double *swrcpMS, - const double swrcpOM[2][SWRC_PARAM_NMAX], + double swrcpOM[][SWRC_PARAM_NMAX], double fom, double depthSapric, double depthT, diff --git a/src/SW_Site.c b/src/SW_Site.c index 567472a3a..462001d08 100644 --- a/src/SW_Site.c +++ b/src/SW_Site.c @@ -1308,7 +1308,7 @@ void SWRC_bulkSoilParameters( unsigned int swrc_type, double *swrcp, const double *swrcpMS, - const double swrcpOM[2][SWRC_PARAM_NMAX], + double swrcpOM[][SWRC_PARAM_NMAX], double fom, double depthSapric, double depthT, From 861950ef074f6254079a1b1f582376efb876326d Mon Sep 17 00:00:00 2001 From: Daniel Schlaepfer Date: Wed, 23 Oct 2024 13:49:07 -0400 Subject: [PATCH 9/9] Fix SWFlowTest.SWFlowSaturatedPercolation (part2) - commit 8029a9d1a70e53405d509ff22fb9608613c38792 introduced for the test SWFlowTest.SWFlowSaturatedPercolation a new variable "swc0init' to initialize soil water content - by mistake, the first use of "swc" was, however, not initialized -> this commit now initializes "swc[]" as intended --- tests/gtests/test_SW_Flow_Lib.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/gtests/test_SW_Flow_Lib.cc b/tests/gtests/test_SW_Flow_Lib.cc index ddb114804..c881ddaa2 100644 --- a/tests/gtests/test_SW_Flow_Lib.cc +++ b/tests/gtests/test_SW_Flow_Lib.cc @@ -172,7 +172,7 @@ TEST(SWFlowTest, SWFlowSaturatedPercolation) { // ***** Tests when nlyrs = 1 ***** // /// provide inputs int nlyrs = 1; - double swc[1]; + double swc[1] = {swc0init}; double swcfc[1] = {1.1}; double swcsat[1] = {1.6}; double impermeability[1] = {0.};