diff --git a/.LSAN_suppr.txt b/.LSAN_suppr.txt index 6384b856e..92813ffca 100644 --- a/.LSAN_suppr.txt +++ b/.LSAN_suppr.txt @@ -5,3 +5,6 @@ # https://github.com/DrylandEcology/SOILWAT2/issues/205 leak:SW_VPD_construct + +# on Apple arm64: https://github.com/google/sanitizers/issues/1501 +leak:realizeClassWithoutSwift diff --git a/NEWS.md b/NEWS.md index 9bcfcaa9d..06be52f77 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,21 @@ # NEWS # SOILWAT2 v7.0.0-9000 +* This version produces nearly identical simulation output + as the previous release under default values for the new inputs. + Small deviations arise due to a fix in the handling of soil moisture values + when between field capacity and saturation. + +* Multiple soil water release curves (`SWRC`) are now implemented and can be + selected with new input `swrc_name`. Implemented `SWRCs` currently include + `"Campbell1974"`, `"vanGenuchten1980"`, and `"FXW"`. New input `has_swrcp` + determines if parameters for a `SWRC` are estimated at run-time via an + implemented pedotransfer function (`PTF`) based on new input `ptf_name` or + if they are provided as inputs via new input file `"swrc_params.in"`. + `rSOILWAT2` implements additional pedotransfer functions. See documentation + entry of `"swrc_ptf"` for additional details and for guidance on how to + implement additional `SWRCs` and `PTFs` (issue #315; @dschlaep). + * Soil density inputs can now represent either matric or bulk density (issue #280; @dschlaep). * Automatic conversion between matric and bulk density as needed @@ -56,6 +71,15 @@ ## Changes to inputs +* New inputs via `"siteparam.in"` select a soil water release curve `swrc_name` + and determine parameter source `has_swrcp`, i.e., + either estimated via selected pedotransfer function `ptf_name` or + values obtained from new input file `"swrc_params.in"`. + Default values `"Campbell1974"`, `"Cosby1984AndOthers"`, and 0 + (i.e., use `PTF` to estimate paramaters) reproduce previous behavior. +* New input file `"swrc_params.in"` to provide parameters of the selected + soil water release curve (if not estimated via a pedotransfer function) + for each soil layer. * SOILWAT2 gains `type_soilDensityInput` as new user input (`siteparam.in`) with default value 0 (i.e., matric soil density) that reproduces previous behavior. diff --git a/SW_Control.c b/SW_Control.c index 960012265..710afa6f0 100644 --- a/SW_Control.c +++ b/SW_Control.c @@ -275,9 +275,19 @@ void SW_CTL_read_inputs_from_disk(void) { if (debug) swprintf(" > 'veg'"); #endif - SW_SIT_read(); // inputs also soil layer data + SW_SIT_read(); #ifdef SWDEBUG - if (debug) swprintf(" > 'site' + 'soils'"); + if (debug) swprintf(" > 'site'"); + #endif + + SW_LYR_read(); + #ifdef RSWDEBUG + if (debug) swprintf(" > 'soils'"); + #endif + + SW_SWRC_read(); + #ifdef SWDEBUG + if (debug) swprintf(" > 'swrc parameters'"); #endif SW_VES_read(); diff --git a/SW_Defines.h b/SW_Defines.h index 586ef1c57..2380eba6d 100644 --- a/SW_Defines.h +++ b/SW_Defines.h @@ -51,13 +51,21 @@ extern "C" { #define SLOW_DRAIN_DEPTH 15. /* numerator over depth in slow drain equation */ /* some basic constants */ -#define MAX_LAYERS 25 -#define MAX_TRANSP_REGIONS 4 -#define MAX_ST_RGR 100 +#define MAX_LAYERS 25 /**< Maximum number of soil layers */ +#define MAX_TRANSP_REGIONS 4 /**< Maximum number of transpiration regions */ +#define MAX_ST_RGR 100 /**< Maximum number of soil temperature nodes */ #define MAX_NYEAR 2500 /**< An integer representing the max calendar year that is supported. The number just needs to be reasonable, it is an artifical limit. */ -#define SW_MISSING 999. /* value to use as MISSING */ +#define SW_MISSING 999. /**< Value to use as MISSING */ + + +// Euler's constant +#ifdef M_E + #define swE M_E +#else + #define swE 2.71828182845904523536028747135266249 +#endif /* M_PI and M_PI_2 from if implementation conforms to POSIX extension @@ -80,7 +88,6 @@ extern "C" { #define deg_to_rad 0.0174532925199433 /**< Convert arc-degrees to radians, i.e., x * deg_to_rad with deg_to_rad = pi / 180 */ #define rad_to_deg 57.29577951308232 /**< Convert radians to arc-degrees, i.e., x * rad_to_deg with rad_to_deg = 180 / pi */ -#define BARCONV 1024. #define SEC_PER_DAY 86400. // the # of seconds in a day... (24 hrs * 60 mins/hr * 60 sec/min = 86400 seconds) @@ -104,7 +111,7 @@ extern "C" { #define SW_MAX 1 /* indices to vegetation types */ -#define NVEGTYPES 4 +#define NVEGTYPES 4 /**< Number of vegetation types implemented */ #define SW_TREES 0 #define SW_SHRUB 1 #define SW_FORBS 2 diff --git a/SW_Files.c b/SW_Files.c index 6b42a83ab..69a615c98 100644 --- a/SW_Files.c +++ b/SW_Files.c @@ -148,53 +148,53 @@ void SW_F_read(const char *s) { #endif switch (lineno) { - case 5: + case 6: strcpy(weather_prefix, inbuf); break; - case 13: + case 14: strcpy(output_prefix, inbuf); break; - case 15: + case 16: InFiles[eOutputDaily] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputDaily]); break; - case 16: + case 17: InFiles[eOutputWeekly] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputWeekly]); //printf("filename: %s \n",InFiles[eOutputWeekly]); break; - case 17: + case 18: InFiles[eOutputMonthly] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputMonthly]); //printf("filename: %s \n",InFiles[eOutputMonthly]); break; - case 18: + case 19: InFiles[eOutputYearly] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputYearly]); break; - case 19: + case 20: InFiles[eOutputDaily_soil] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputDaily_soil]); //printf("filename: %s \n",InFiles[eOutputDaily]); break; - case 20: + case 21: InFiles[eOutputWeekly_soil] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputWeekly_soil]); //printf("filename: %s \n",InFiles[eOutputWeekly]); break; - case 21: + case 22: InFiles[eOutputMonthly_soil] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputMonthly_soil]); //printf("filename: %s \n",InFiles[eOutputMonthly]); break; - case 22: + case 23: InFiles[eOutputYearly_soil] = Str_Dup(inbuf); ++fileno; SW_CSV_F_INIT(InFiles[eOutputYearly_soil]); diff --git a/SW_Files.h b/SW_Files.h index 6c3e10c99..c0b97d50a 100644 --- a/SW_Files.h +++ b/SW_Files.h @@ -20,16 +20,30 @@ extern "C" { #endif -#define SW_NFILES 22 +#define SW_NFILES 23 /* The number of enum elements between eNoFile and * eEndFile (not inclusive) must match SW_NFILES. * also, these elements must match the order of - * input from files.in. + * input from `files.in`. */ typedef enum { - eNoFile = -1, eFirst = 0, eModel, eLog, eSite, eLayers, eWeather, - eMarkovProb, eMarkovCov, eSky, eVegProd, eVegEstab, eCarbon, eSoilwat, + eNoFile = -1, + /* List of all input files */ + eFirst = 0, + /* Description of a model run */ + eModel, eLog, + /* Description of simulated site */ + eSite, eLayers, eSWRCp, + /* Weather and climate forcing */ + eWeather, eMarkovProb, eMarkovCov, eSky, + /* Description of vegetation */ + eVegProd, eVegEstab, + /* Description of CO2 effects */ + eCarbon, + /* (optional) soil moisture measurements */ + eSoilwat, + /* Simulation outputs */ eOutput, eOutputDaily, eOutputWeekly, eOutputMonthly, eOutputYearly, eOutputDaily_soil, eOutputWeekly_soil, eOutputMonthly_soil, eOutputYearly_soil, eEndFile diff --git a/SW_Flow_lib.c b/SW_Flow_lib.c index 9c6215a28..480594907 100644 --- a/SW_Flow_lib.c +++ b/SW_Flow_lib.c @@ -372,7 +372,7 @@ void transp_weighted_avg(double *swp_avg, unsigned int n_tr_rgns, unsigned int n for (i = 0; i < n_layers; i++) { if (tr_regions[i] == r) { - swp += tr_coeff[i] * SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i); + swp += tr_coeff[i] * SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]); sumco += tr_coeff[i]; } } @@ -484,10 +484,10 @@ void pot_soil_evap(double *bserate, unsigned int nelyrs, double ecoeff[], double } x = width[i] * ecoeff[i]; sumwidth += x; - avswp += x * SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i); + avswp += x * SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]); } - // Note: avswp = 0 if swc = 0 because that is the return value of SW_SWCbulk2SWPmatric + // Note: avswp = 0 if swc = 0 because that is the return value of SW_SWRC_SWCtoSWP avswp /= (ZRO(sumwidth)) ? 1 : sumwidth; /* 8/27/92 (SLC) if totagb > Es_param_limit, assume soil surface is @@ -542,7 +542,7 @@ void pot_soil_evap_bs(double *bserate, unsigned int nelyrs, double ecoeff[], dou for (i = 0; i < nelyrs; i++) { x = width[i] * ecoeff[i]; sumwidth += x; - avswp += x * SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i); + avswp += x * SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]); } avswp /= sumwidth; @@ -742,12 +742,18 @@ void remove_from_soil(double swc[], double qty[], double *aet, unsigned int nlyr **********************************************************************/ unsigned int i; - double swpfrac[MAX_LAYERS], sumswp = 0.0, swc_avail, q; + double swpfrac[MAX_LAYERS], sumswp = 0.0, swc_avail, q, tmpswp; SW_SOILWAT *st = &SW_Soilwat; for (i = 0; i < nlyrs; i++) { - swpfrac[i] = coeff[i] / SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, swc[i], i); + tmpswp = SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]); + if (GT(tmpswp, 0.)) { + swpfrac[i] = coeff[i] / tmpswp; + } else { + // saturated conditions -> divide by SWP [-bar] at field capacity + swpfrac[i] = coeff[i] / 0.333; + } sumswp += swpfrac[i]; } @@ -969,7 +975,7 @@ void hydraulic_redistribution( swc[i] - fmin(lyr[i]->swcBulk_wiltpt, lyr[i]->swcBulk_atSWPcrit[vegk]) ); - swp[i] = SW_SWCbulk2SWPmatric(lyr[i]->fractionVolBulk_gravel, swc[i], i); + swp[i] = SW_SWRC_SWCtoSWP(swc[i], SW_Site.lyr[i]); /* Ryel et al. 2002: eq. 7 relative soil-root conductance */ relCondroot[i] = fmin( diff --git a/SW_Output.h b/SW_Output.h index e83539c10..8a7e7bae6 100644 --- a/SW_Output.h +++ b/SW_Output.h @@ -234,6 +234,7 @@ extern IntUS ncol_OUT[SW_OUTNKEYS]; extern char const *key2str[]; extern char const *pd2longstr[]; +extern char const *styp2str[]; diff --git a/SW_Output_get_functions.c b/SW_Output_get_functions.c index d31fdc482..4fe1baccd 100644 --- a/SW_Output_get_functions.c +++ b/SW_Output_get_functions.c @@ -979,8 +979,8 @@ void get_swpMatric_text(OutPeriod pd) ForEachSoilLayer(i) { /* swpMatric at this point is identical to swcBulk */ - val = SW_SWCbulk2SWPmatric(SW_Site.lyr[i]->fractionVolBulk_gravel, - vo->swpMatric[i], i); + val = SW_SWRC_SWCtoSWP(vo->swpMatric[i], SW_Site.lyr[i]); + snprintf(str, OUTSTRLEN, "%c%.*f", _Sep, OUT_DIGITS, val); strcat(sw_outstr, str); @@ -1006,8 +1006,7 @@ void get_swpMatric_mem(OutPeriod pd) ForEachSoilLayer(i) { /* swpMatric at this point is identical to swcBulk */ - p[iOUT(i, pd)] = SW_SWCbulk2SWPmatric( - SW_Site.lyr[i]->fractionVolBulk_gravel, vo->swpMatric[i], i); + p[iOUT(i, pd)] = SW_SWRC_SWCtoSWP(vo->swpMatric[i], SW_Site.lyr[i]); } } @@ -1031,8 +1030,7 @@ void get_swpMatric_agg(OutPeriod pd) ForEachSoilLayer(i) { /* swpMatric at this point is identical to swcBulk */ - val = SW_SWCbulk2SWPmatric( - SW_Site.lyr[i]->fractionVolBulk_gravel, vo->swpMatric[i], i); + val = SW_SWRC_SWCtoSWP(vo->swpMatric[i], SW_Site.lyr[i]); do_running_agg(p, psd, iOUT(i, pd), Globals->currIter, val); } diff --git a/SW_Site.c b/SW_Site.c index b0b59fe13..821840da9 100644 --- a/SW_Site.c +++ b/SW_Site.c @@ -91,6 +91,32 @@ RealD _SWCMinVal; /* lower bound on swc. */ +/** Character representation of implemented Soil Water Retention Curves (SWRC) + + @note Code maintenance: + - Values must exactly match those provided in `siteparam.in`. + - Order must exactly match "indices of `swrc2str`" + - See details in section \ref swrc_ptf +*/ +char const *swrc2str[N_SWRCs] = { + "Campbell1974", + "vanGenuchten1980", + "FXW" +}; + +/** Character representation of implemented Pedotransfer Functions (PTF) + + @note Code maintenance: + - Values must exactly match those provided in `siteparam.in`. + - Order must exactly match "indices of `ptf2str`" + - See details in section \ref swrc_ptf + - `rSOILWAT2` may implemented additional PTFs +*/ +char const *ptf2str[N_PTFs] = { + "Cosby1984AndOthers", + "Cosby1984" +}; + /* =================================================== */ /* Local Variables */ @@ -102,186 +128,887 @@ static char *MyFileName; /* Local Function Definitions */ /* --------------------------------------------------- */ -static void _read_layers(void) { - /* =================================================== */ - /* 5-Feb-2002 (cwb) removed dmin requirement in input file */ +/** Check validity of soil properties - SW_SITE *v = &SW_Site; - FILE *f; - LyrIndex lyrno; - int x, k; - RealF dmin = 0.0, dmax, evco, trco_veg[NVEGTYPES], psand, pclay, soildensity, imperm, - soiltemp, f_gravel; + @param[in] *lyr Soil information. - /* note that Files.read() must be called prior to this. */ - MyFileName = SW_F_name(eLayers); + @return A logical value indicating if soil properties passed the checks. +*/ +static Bool SW_check_soil_properties(SW_LAYER_INFO *lyr) { + int k; + RealD fval = 0; + const char *errtype = "\0"; + Bool res = swTRUE; - f = OpenFile(MyFileName, "r"); - while (GetALine(f, inbuf)) { - lyrno = _newlayer(); + if (LE(lyr->width, 0.)) { + res = swFALSE; + fval = lyr->width; + errtype = Str_Dup("layer width"); - x = sscanf( - inbuf, - "%f %f %f %f %f %f %f %f %f %f %f %f", - &dmax, - &soildensity, - &f_gravel, - &evco, - &trco_veg[SW_GRASS], &trco_veg[SW_SHRUB], &trco_veg[SW_TREES], &trco_veg[SW_FORBS], - &psand, - &pclay, - &imperm, - &soiltemp + } else if ( + LT(lyr->soilDensityInput, 0.) || + GT(lyr->soilDensityInput, 2.65) + ) { + res = swFALSE; + fval = lyr->soilDensityInput; + errtype = Str_Dup("soil density"); + + } else if ( + LT(lyr->fractionVolBulk_gravel, 0.) || + GE(lyr->fractionVolBulk_gravel, 1.) + ) { + res = swFALSE; + fval = lyr->fractionVolBulk_gravel; + errtype = Str_Dup("gravel content"); + + } else if ( + LE(lyr->fractionWeightMatric_sand, 0.) || + GE(lyr->fractionWeightMatric_sand, 1.) + ) { + res = swFALSE; + fval = lyr->fractionWeightMatric_sand; + errtype = Str_Dup("sand proportion"); + + } else if ( + LE(lyr->fractionWeightMatric_clay, 0.) || + GE(lyr->fractionWeightMatric_clay, 1.) + ) { + res = swFALSE; + fval = lyr->fractionWeightMatric_clay; + errtype = Str_Dup("clay proportion"); + + } else if ( + GE(lyr->fractionWeightMatric_sand + lyr->fractionWeightMatric_clay, 1.) + ) { + res = swFALSE; + fval = lyr->fractionWeightMatric_sand + lyr->fractionWeightMatric_clay; + errtype = Str_Dup("sand+clay proportion"); + + } else if ( + LT(lyr->impermeability, 0.) || + GT(lyr->impermeability, 1.) + ) { + res = swFALSE; + fval = lyr->impermeability; + errtype = Str_Dup("impermeability"); + + } else if ( + LT(lyr->evap_coeff, 0.) || + GT(lyr->evap_coeff, 1.) + ) { + res = swFALSE; + fval = lyr->evap_coeff; + errtype = Str_Dup("bare-soil evaporation coefficient"); + + } else { + ForEachVegType(k) { + if (LT(lyr->transp_coeff[k], 0.) || GT(lyr->transp_coeff[k], 1.)) { + res = swFALSE; + fval = lyr->transp_coeff[k]; + errtype = Str_Dup("transpiration coefficient"); + break; + } + } + } + + if (!res) { + LogError( + logfp, + LOGFATAL, + "'%s' has an invalid value (%5.4f) in layer %d.\n", + errtype, fval, lyr->id + 1 ); + } - /* Check that we have 12 values per layer */ - /* Adjust number if new variables are added */ - if (x != 12) { - CloseFile(&f); + return res; +} + + + + +/** A realistic lower limit for minimum `theta` + +@note Currently, -30 MPa + (based on limited test runs across western US including hot deserts) + lower than "air-dry" = hygroscopic point (-10. MPa; Porporato et al. 2001) + not as extreme as "oven-dry" (-1000. MPa; Fredlund et al. 1994) +*/ +static double lower_limit_of_theta_min( + unsigned int swrc_type, + double *swrcp, + double gravel, + double width +) { + double res = SWRC_SWPtoSWC(300., swrc_type, swrcp, gravel, width, LOGFATAL); + + // convert bulk [cm] to matric [cm / cm] + return res / ((1. - gravel) * width); +} + +/** + @brief User-input based minimum volumetric water content + + This function returns minimum soil water content `theta_min` determined + by user input via #_SWCMinVal as + * a fixed VWC value, + * a fixed SWP value, or + * a realistic lower limit from `lower_limit_of_theta_min()`, + unless legacy mode (`ptf` equal to "Cosby1984AndOthers") is used + (reproducing behavior prior to v7.0.0), 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). + + @param[in] ui_sm_min User input of requested minimum soil moisture, + see #_SWCMinVal + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @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] 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 + @param[in] legacy_mode If true then legacy behavior (see details) + + @return Minimum volumetric water content of the matric soil [cm / cm] +*/ +static double ui_theta_min( + double ui_sm_min, + double gravel, + double width, + double sand, + double clay, + double swcBulk_sat, + unsigned int swrc_type, + double *swrcp, + Bool legacy_mode +) { + double vwc_min = SW_MISSING, tmp_vwcmin; + + if (LT(ui_sm_min, 0.0)) { + /* user input: request to estimate minimum theta */ + + /* realistic lower limit for theta_min */ + vwc_min = lower_limit_of_theta_min(swrc_type, swrcp, gravel, width); + + if (legacy_mode) { + /* residual theta estimated with Rawls & Brakensiek (1985) PTF */ + PTF_RawlsBrakensiek1985( + &tmp_vwcmin, + sand, + clay, + swcBulk_sat / ((1. - gravel) * width) + ); + + /* if `PTF_RawlsBrakensiek1985()` was successful, then take max */ + if (!missing(tmp_vwcmin)) { + vwc_min = fmax(vwc_min, tmp_vwcmin); + } + } + + } else if (GE(_SWCMinVal, 1.0)) { + /* user input: fixed (matric) SWP value; unit(_SWCMinVal) == -bar */ + vwc_min = SWRC_SWPtoSWC( + _SWCMinVal, + swrc_type, + swrcp, + gravel, + width, + LOGFATAL + ) / ((1. - gravel) * width); + + } else { + /* user input: fixed matric VWC; unit(_SWCMinVal) == cm/cm */ + vwc_min = _SWCMinVal / width; + } + + return vwc_min; +} + + + +/* =================================================== */ +/* Global Function Definitions */ +/* --------------------------------------------------- */ + + +/** + @brief Translate a SWRC name into a SWRC type number + + See #swrc2str and `siteparam.in`. + Throws an error if `SWRC` is not implemented. + + @param[in] *swrc_name Name of a SWRC + + @return Internal identification number of selected SWRC +*/ +unsigned int encode_str2swrc(char *swrc_name) { + unsigned int k; + + for ( + k = 0; + k < N_SWRCs && Str_CompareI(swrc_name, (char *)swrc2str[k]); + k++ + ); + + if (k == N_SWRCs) { + LogError( + logfp, + LOGFATAL, + "SWRC '%s' is not implemented.", + swrc_name + ); + } + + return k; +} + +/** + @brief Translate a PTF name into a PTF type number + + See #ptf2str and `siteparam.in`. + + @param[in] *ptf_name Name of a PTF + + @return Internal identification number of selected PTF; + #SW_MISSING if not implemented. +*/ +unsigned int encode_str2ptf(char *ptf_name) { + unsigned int k; + + for ( + k = 0; + k < N_PTFs && Str_CompareI(ptf_name, (char *)ptf2str[k]); + k++ + ); + + if (k == N_PTFs) { + k = (unsigned int) SW_MISSING; + } + + return k; +} + + + +/** + @brief Estimate parameters of selected soil water retention curve (SWRC) + using selected pedotransfer function (PTF) + + 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[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) + of the whole soil [m3/m3] + @param[in] bdensity Density of the whole soil + (matric soil plus coarse fragments) [g/cm3]; + accepts #SW_MISSING if not used by selected PTF + +*/ +void SWRC_PTF_estimate_parameters( + unsigned int ptf_type, + double *swrcp, + double sand, + double clay, + double gravel, + double bdensity +) { + + /* Initialize swrcp[] to 0 */ + memset(swrcp, 0., SWRC_PARAM_NMAX * sizeof(swrcp[0])); + + if (ptf_type == sw_Cosby1984AndOthers) { + SWRC_PTF_Cosby1984_for_Campbell1974(swrcp, sand, clay); + + } else if (ptf_type == sw_Cosby1984) { + SWRC_PTF_Cosby1984_for_Campbell1974(swrcp, sand, clay); + + } else { + LogError( + logfp, + LOGFATAL, + "PTF is not implemented in SOILWAT2.", + ptf_type + ); + } + + /**********************************/ + /* TODO: remove once PTFs are implemented that utilize gravel */ + /* avoiding `error: unused parameter 'gravel' [-Werror=unused-parameter]` */ + if (gravel < 0.) {}; + + /* TODO: remove once PTFs are implemented that utilize bdensity */ + /* avoiding `error: unused parameter 'gravel' [-Werror=unused-parameter]` */ + if (bdensity < 0.) {}; + /**********************************/ +} + + +/** + @brief Estimate Campbell's 1974 SWRC parameters + using Cosby et al. 1984 multivariate PTF + + Estimation of four SWRC parameter values `swrcp` + based on sand, clay, and (silt). + + Parameters are explained in `SWRC_check_parameters_for_Campbell1974()`. + + Multivariate PTFs are from Cosby et al. 1984 (\cite Cosby1984) Table 4; + Cosby et al. 1984 provided also univariate PTFs in Table 5 + 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[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 +) { + /* 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 - 6.4 * clay - 0.60); +} + + +/** + @brief Saturated soil water content + + See #ptf2str for implemented PTFs. + See #swrc2str for implemented SWRCs. + + 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 "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 + `ptf_name` is "Cosby1984AndOthers" (see #ptf2str). + + @param[in] swrc_type Identification number of selected SWRC + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @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] + + @return Estimated saturated water content of the bulk soil [cm] +*/ +double SW_swcBulk_saturated( + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + unsigned int ptf_type, + double sand, + double clay +) { + double theta_sat = SW_MISSING; + + switch (swrc_type) { + case sw_Campbell1974: + if (ptf_type == sw_Cosby1984AndOthers) { + // Cosby1984AndOthers (backwards compatible) + PTF_Saxton2006(&theta_sat, sand, clay); + } else { + theta_sat = swrcp[1]; + } + break; + + case sw_vanGenuchten1980: + theta_sat = swrcp[1]; + break; + + case sw_FXW: + theta_sat = swrcp[0]; + break; + + default: LogError( logfp, LOGFATAL, - "%s : Incomplete record %d.\n", - MyFileName, lyrno + 1 + "`SW_swcBulk_saturated()`: SWRC (type %d) is not implemented.", + swrc_type ); - } + break; + } - v->lyr[lyrno]->width = dmax - dmin; + // Convert from matric [cm/cm] to bulk [cm] + return theta_sat * width * (1. - gravel); +} - /* checks for valid values now carried out by `SW_SIT_init_run()` */ - dmin = dmax; - v->lyr[lyrno]->fractionVolBulk_gravel = f_gravel; - v->lyr[lyrno]->soilDensityInput = soildensity; - v->lyr[lyrno]->evap_coeff = evco; +/** + @brief Lower limit of simulated soil water content + + See #ptf2str for implemented PTFs. + See #swrc2str for implemented SWRCs. + + SWRCs provide a theoretical minimum/residual volumetric water content; + however, water potential at that minimum value is infinite. Instead, we + 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` + are utilized only in legacy mode, i.e., `ptf` equal to "Cosby1984AndOthers". + + @param[in] swrc_type Identification number of selected SWRC + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @param[in] ptf_type Identification number of selected PTF + @param[in] ui_sm_min User input of requested minimum soil moisture, + see #_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] swcBulk_sat Saturated water content of the bulk soil [cm] + + @return Minimum water content of the bulk soil [cm] +*/ +double SW_swcBulk_minimum( + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + unsigned int ptf_type, + double ui_sm_min, + double sand, + double clay, + double swcBulk_sat +) { + double theta_min_sim, theta_min_theoretical = SW_MISSING; + + /* `theta_min` based on theoretical SWRC */ + switch (swrc_type) { + case sw_Campbell1974: // phi = infinity at theta_min + theta_min_theoretical = 0.; + break; - ForEachVegType(k) - { - v->lyr[lyrno]->transp_coeff[k] = trco_veg[k]; - } + case sw_vanGenuchten1980: // phi = infinity at theta_min + theta_min_theoretical = swrcp[0]; + break; - v->lyr[lyrno]->fractionWeightMatric_sand = psand; - v->lyr[lyrno]->fractionWeightMatric_clay = pclay; - v->lyr[lyrno]->impermeability = imperm; - v->lyr[lyrno]->avgLyrTemp = soiltemp; + case sw_FXW: // phi = 6.3 x 10^6 cm at theta_min + theta_min_theoretical = 0.; + break; - if (lyrno >= MAX_LAYERS) { - CloseFile(&f); + default: LogError( logfp, LOGFATAL, - "%s : Too many layers specified (%d).\n" - "Maximum number of layers is %d\n", - MyFileName, lyrno + 1, MAX_LAYERS + "`SW_swcBulk_minimum()`: SWRC (type %d) is not implemented.", + swrc_type ); + break; + } + + /* `theta_min` based on user input `ui_sm_min` */ + theta_min_sim = ui_theta_min( + ui_sm_min, + gravel, + width, + sand, + clay, + 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 + ); + + /* `theta_min_sim` must be strictly larger than `theta_min_theoretical` */ + theta_min_sim = fmax(theta_min_sim, theta_min_theoretical + D_DELTA); + + /* Convert from matric [cm/cm] to bulk [cm] */ + return theta_min_sim * width * (1. - gravel); +} + + + +/** + @brief Check whether selected PTF and SWRC are compatible + + See #ptf2str for implemented PTFs. + See #swrc2str for implemented SWRCs. + + @param[in] *swrc_name Name of selected SWRC + @param[in] *ptf_name Name of selected PTF + + @return A logical value indicating if SWRC and PTF are compatible. +*/ +Bool check_SWRC_vs_PTF(char *swrc_name, char *ptf_name) { + Bool res = swFALSE; + + if (!missing((double) encode_str2ptf(ptf_name))) { + if ( + Str_CompareI(swrc_name, (char *) "Campbell1974") == 0 && + ( + Str_CompareI(ptf_name, (char *) "Cosby1984AndOthers") == 0 || + Str_CompareI(ptf_name, (char *) "Cosby1984") == 0 + ) + ) { + res = swTRUE; } } - CloseFile(&f); + return res; } +/** + @brief Check Soil Water Retention Curve (SWRC) parameters + + See #swrc2str for implemented SWRCs. + + @param[in] swrc_type Identification number of selected SWRC + @param[in] *swrcp Vector of SWRC parameters + + @return A logical value indicating if parameters passed the checks. +*/ +Bool SWRC_check_parameters(unsigned int swrc_type, double *swrcp) { + Bool res = swFALSE; + + switch (swrc_type) { + case sw_Campbell1974: + res = SWRC_check_parameters_for_Campbell1974(swrcp); + break; + + case sw_vanGenuchten1980: + res = SWRC_check_parameters_for_vanGenuchten1980(swrcp); + break; + + case sw_FXW: + res = SWRC_check_parameters_for_FXW(swrcp); + break; + + default: + LogError( + logfp, + LOGFATAL, + "Selected SWRC (type %d) is not implemented.", + swrc_type + ); + break; + } + + return res; +} + -/* =================================================== */ -/* Global Function Definitions */ -/* --------------------------------------------------- */ /** - \brief Calculate soil moisture characteristics for each layer. + @brief Check Campbell's 1974 SWRC parameters - Bulk refers to the whole soil, i.e., including the rock/gravel component, - whereas matric refers to the < 2 mm fraction. + See `SWRC_SWCtoSWP_Campbell1974()` and `SWRC_SWPtoSWC_Campbell1974()` + for implementation of Campbell's 1974 SWRC (\cite Campbell1974). - Saturated moisture content of the matric component (thetasMatric), saturation matric - potential (psisMatric), and the slope of the retention curve (bMatric) for each - layer are calculated using equations found in Cosby et al. (1984). \cite Cosby1984 - The saturated moisture content of the whole soil (bulk) for each layer (swcBulk_saturated) - is calculated using equations found in Saxton and Rawls (2006; Equations 2, 3 & 5). - \cite Saxton2006 + Campbell1974 has four parameters (three are used for the SWRC): + - `swrcp[0]` (`psisMatric`): air-entry suction [cm] + - `swrcp[1]` (previously named `thetasMatric`): + saturated volumetric water content for the matric component [cm/cm] + - `swrcp[2]` (`bMatric`): slope of the linear log-log retention curve [-] + - `swrcp[3]` (`K_sat`): saturated hydraulic conductivity `[cm / day]` - Return from the function is void. Calculated values stored in SW_Site object. + @param[in] *swrcp Vector of SWRC parameters - SOILWAT2 calculates internally based on soil bulk density of the whole soil, - i.e., including rock/gravel component. However, inputs are expected to - represent soil (matric) density of the < 2 mm fraction. + @return A logical value indicating if parameters passed the checks. +*/ +Bool SWRC_check_parameters_for_Campbell1974(double *swrcp) { + Bool res = swTRUE; - sand + clay + silt must equal one. Fraction silt is calculated: 1 - (sand + clay). + if (LE(swrcp[0], 0.0)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_Campbell1974(): invalid value of " + "psi(saturated, matric, [cm]) = %f (must > 0)\n", + swrcp[1] + ); + } + + if (LE(swrcp[1], 0.0) || GT(swrcp[1], 1.0)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_Campbell1974(): invalid value of " + "theta(saturated, matric, [cm/cm]) = %f (must be within 0-1)\n", + swrcp[1] + ); + } + + if (ZRO(swrcp[2])) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_Campbell1974(): invalid value of " + "beta = %f (must be != 0)\n", + swrcp[2] + ); + } + + if (LE(swrcp[3], 0.0)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_Campbell1974(): invalid value of " + "K_sat = %f (must be > 0)\n", + swrcp[3] + ); + } + + return res; +} + +/** + @brief Check van Genuchten 1980 SWRC parameters + + See `SWRC_SWCtoSWP_vanGenuchten1980()` and `SWRC_SWPtoSWC_vanGenuchten1980()` + for implementation of van Genuchten's 1980 SWRC (\cite vanGenuchten1980). - \param fractionGravel The fraction of gravel in a layer by volume. - \param sand The fraction of sand in a layer by weight. - \param clay The fraction of clay in a layer by weight. - \param n Soil layer index. + "vanGenuchten1980" has five parameters (four are used for the SWRC): + - `swrcp[0]` (`theta_r`): residual volumetric water content + of the matric component [cm/cm] + - `swrcp[1]` (`theta_s`): saturated volumetric water content + of the matric component [cm/cm] + - `swrcp[2]` (`alpha`): related to the inverse of air entry suction [cm-1] + - `swrcp[3]` (`n`): measure of the pore-size distribution [-] + - `swrcp[4]` (`K_sat`): saturated hydraulic conductivity `[cm / day]` - \sideeffect - - thetasMatric Saturated water content for the matric component (m^3/m^3). - - psisMatric Saturation matric potential (MPa). - - bMatric Slope of the linear log-log retention curve (unitless). - - swcBulk_saturated The saturated water content for the whole soil (bulk) (cm/layer). + @param[in] *swrcp Vector of SWRC parameters + @return A logical value indicating if parameters passed the checks. */ +Bool SWRC_check_parameters_for_vanGenuchten1980(double *swrcp) { + Bool res = swTRUE; -void water_eqn(RealD fractionGravel, RealD sand, RealD clay, LyrIndex n) { + if (LE(swrcp[0], 0.0) || GT(swrcp[0], 1.)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_vanGenuchten1980(): invalid value of " + "theta(residual, matric, [cm/cm]) = %f (must be within 0-1)\n", + swrcp[0] + ); + } - /* Cosby, B. J., G. M. Hornberger, R. B. Clapp, and T. R. Ginn. 1984. - A statistical exploration of the relationships of soil moisture - characteristics to the physical properties of soils. - Water Resources Research 20:682–690. - https://doi.org/10.1029/WR020i006p00682. */ + if (LE(swrcp[1], 0.0) || GT(swrcp[1], 1.)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_vanGenuchten1980(): invalid value of " + "theta(saturated, matric, [cm/cm]) = %f (must be within 0-1)\n", + swrcp[1] + ); + } - /* Table 4 */ - SW_Site.lyr[n]->thetasMatric = -14.2 * sand - 3.7 * clay + 50.5; - SW_Site.lyr[n]->psisMatric = powe(10.0, -1.58 * sand - 0.63 * clay + 2.17); - SW_Site.lyr[n]->bMatric = -0.3 * sand + 15.7 * clay + 3.10; + if (LE(swrcp[1], swrcp[0])) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_vanGenuchten1980(): invalid values for " + "theta(residual, matric, [cm/cm]) = %f and " + "theta(saturated, matric, [cm/cm]) = %f " + "(must be theta_r < theta_s)\n", + swrcp[0], swrcp[1] + ); + } + if (LE(swrcp[2], 0.0)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_vanGenuchten1980(): invalid value of " + "alpha([1 / cm]) = %f (must > 0)\n", + swrcp[2] + ); + } - if ( - LE(SW_Site.lyr[n]->thetasMatric, 0.0) || - GT(SW_Site.lyr[n]->thetasMatric, 100.0) - ) { + if (LE(swrcp[3], 1.0)) { + res = swFALSE; LogError( logfp, - LOGFATAL, - "water_eqn(): invalid value of " - "theta(saturated, matric, [%]; Cosby et al. 1984) = %f " - "(must within 0-100%)\n", - SW_Site.lyr[n]->thetasMatric + LOGWARN, + "SWRC_check_parameters_for_vanGenuchten1980(): invalid value of " + "n = %f (must be > 1)\n", + swrcp[3] ); } - if (ZRO(SW_Site.lyr[n]->bMatric)) { + if (LE(swrcp[4], 0.0)) { + res = swFALSE; LogError( logfp, - LOGFATAL, - "water_eqn(): invalid value of beta = %f (must be != 0)\n", - SW_Site.lyr[n]->bMatric + LOGWARN, + "SWRC_check_parameters_for_vanGenuchten1980(): invalid value of " + "K_sat = %f (must be > 0)\n", + swrcp[4] ); } - SW_Site.lyr[n]->binverseMatric = 1.0 / SW_Site.lyr[n]->bMatric; + return res; +} - /* Saxton, K. E. and W. J. Rawls. 2006. Soil water characteristic estimates - by texture and organic matter for hydrologic solutions. - Soil Science Society of America Journal 70:1569-1578. - */ +/** + @brief Check FXW SWRC parameters - RealD - OM = 0., - theta_S, theta_33, theta_33t, theta_S33, theta_S33t, theta_1500, theta_1500t, - R_w, alpha; + See `SWRC_SWCtoSWP_FXW()` and `SWRC_SWPtoSWC_FXW()` + for implementation of FXW's SWRC + (\cite fredlund1994CGJa, \cite wang2018wrr). - /* 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; + "FXWJD" is a synonym for the FXW SWRC (\cite wang2018wrr). - theta_1500 = theta_1500t + (0.14 * theta_1500t - 0.02); + FXW has six parameters (four are used for the SWRC): + - `swrcp[0]` (`theta_s`): saturated volumetric water content + of the matric component [cm/cm] + - `swrcp[1]` (`alpha`): shape parameter [cm-1] + - `swrcp[2]` (`n`): shape parameter [-] + - `swrcp[3]` (`m`): shape parameter [-] + - `swrcp[4]` (`K_sat`): saturated hydraulic conductivity [cm / day] + - `swrcp[5]` (`L`): tortuosity/connectivity parameter [-] + + FXW SWRC does not include a non-zero residual water content parameter; + instead, it assumes a pressure head #FXW_h0 of 6.3 x 10^6 cm + (equivalent to c. -618 MPa) at zero water content + and a "residual" pressure head #FXW_hr of 1500 cm + (\cite fredlund1994CGJa, \cite wang2018wrr, \cite rudiyanto2021G). + + Acceptable parameter values are partially based on + Table 1 in \cite wang2022WRRa. + + @param[in] *swrcp Vector of SWRC parameters + + @return A logical value indicating if parameters passed the checks. +*/ +Bool SWRC_check_parameters_for_FXW(double *swrcp) { + Bool res = swTRUE; + + if (LE(swrcp[0], 0.0) || GT(swrcp[0], 1.)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_FXW(): invalid value of " + "theta(saturated, matric, [cm/cm]) = %f (must be within 0-1)\n", + swrcp[0] + ); + } + + if (LE(swrcp[1], 0.0)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_FXW(): invalid value of " + "alpha([1 / cm]) = %f (must be > 0)\n", + swrcp[1] + ); + } + + if (LE(swrcp[2], 1.0) || GT(swrcp[2], 10.)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_FXW(): invalid value of " + "n = %f (must be within 1-10)\n", + swrcp[2] + ); + } + + if (LE(swrcp[3], 0.0) || GT(swrcp[3], 1.5)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_FXW(): invalid value of " + "m = %f (must be within 0-1.5)\n", + swrcp[3] + ); + } + + if (LE(swrcp[4], 0.0)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_FXW(): invalid value of " + "K_sat = %f (must be > 0)\n", + swrcp[4] + ); + } + + if (LE(swrcp[5], 0.)) { + res = swFALSE; + LogError( + logfp, + LOGWARN, + "SWRC_check_parameters_for_FXW(): invalid value of " + "L = %f (must be > 0)\n", + swrcp[5] + ); + } + + return res; +} + + +/** + @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] +*/ +void PTF_Saxton2006( + double *theta_sat, + double sand, + double clay +) { + double + OM = 0., + theta_33, theta_33t, theta_S33, theta_S33t; /* Eq. 2: 33 kPa moisture */ theta_33t = @@ -293,7 +1020,8 @@ void water_eqn(RealD fractionGravel, RealD sand, RealD clay, LyrIndex n) { - 0.027 * clay * OM \ + 0.452 * sand * clay; - theta_33 = theta_33t + (1.283 * squared(theta_33t) - 0.374 * theta_33t - 0.015); + theta_33 = + theta_33t + (1.283 * squared(theta_33t) - 0.374 * theta_33t - 0.015); /* Eq. 3: SAT-33 kPa moisture */ theta_S33t = @@ -309,24 +1037,39 @@ void water_eqn(RealD fractionGravel, RealD sand, RealD clay, LyrIndex n) { /* Eq. 5: saturated moisture */ - theta_S = theta_33 + theta_S33 - 0.097 * sand + 0.043; + *theta_sat = theta_33 + theta_S33 - 0.097 * sand + 0.043; if ( - LE(theta_S, 0.) || - GT(theta_S, 1.) + LE(*theta_sat, 0.) || + GT(*theta_sat, 1.) ) { LogError( logfp, LOGFATAL, - "water_eqn(): invalid value of " - "theta(saturated, [cm / cm]; Saxton et al. 2006) = %f " - "(must be within 0-1)\n", - theta_S + "PTF_Saxton2006(): invalid value of " + "theta(saturated, [cm / cm]) = %f (must be within 0-1)\n", + *theta_sat ); } - SW_Site.lyr[n]->swcBulk_saturated = - SW_Site.lyr[n]->width * (1. - fractionGravel) * theta_S; + + +// currently, unused and defunct code: +// (e.g., SW_Site.lyr[n] assignments don't work here anymore!): +#ifdef UNUSED_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_1500 = theta_1500t + (0.14 * theta_1500t - 0.02); /* Eq. 18: slope of logarithmic tension-moisture curve */ @@ -361,12 +1104,71 @@ void water_eqn(RealD fractionGravel, RealD sand, RealD clay, LyrIndex n) { SW_Site.lyr[n]->Saxton2006_K_sat_bulk *= SW_Site.lyr[n]->Saxton2006_fK_gravel; - } else { - SW_Site.lyr[n]->Saxton2006_fK_gravel = 1.; + } else { + SW_Site.lyr[n]->Saxton2006_fK_gravel = 1.; + } +#endif +} + + + +/** + @brief Rawls and Brakensiek 1985 PTFs \cite rawls1985WmitE + to estimate residual soil water content for the Brooks-Corey SWRC + \cite brooks1964a + + @note This function was previously named "SW_VWCBulkRes". + @note This function is only well-defined for `clay` values in 5-60%, + `sand` values in 5-70%, and `porosity` must in 10-<100%. + + @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] 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] +*/ +void PTF_RawlsBrakensiek1985( + double *theta_min, + double sand, + double clay, + double porosity +) { + if ( + GE(clay, 0.05) && LE(clay, 0.6) && + GE(sand, 0.05) && LE(sand, 0.7) && + GE(porosity, 0.1) && LT(porosity, 1.) + ) { + sand *= 100.; + clay *= 100.; + /* Note: the equation requires sand and clay in units of [100 * g / g]; + porosity, however, must be in units of [cm3 / cm3] + */ + *theta_min = fmax( + 0., + - 0.0182482 \ + + 0.00087269 * sand \ + + 0.00513488 * clay \ + + 0.02939286 * porosity \ + - 0.00015395 * squared(clay) \ + - 0.0010827 * sand * porosity \ + - 0.00018233 * squared(clay) * squared(porosity) \ + + 0.00030703 * squared(clay) * porosity \ + - 0.0023584 * squared(porosity) * clay + ); + + } else { + LogError( + logfp, + LOGWARN, + "`PTF_RawlsBrakensiek1985()`: sand, clay, or porosity out of valid range." + ); + + *theta_min = SW_MISSING; } - } + + /** @brief Estimate soil density of the whole soil (bulk). @@ -505,6 +1307,7 @@ LyrIndex _newlayer(void) { : (SW_LAYER_INFO **) Mem_ReAlloc(v->lyr, sizeof(SW_LAYER_INFO *) * (v->n_layers)); /* else realloc() */ v->lyr[v->n_layers - 1] = (SW_LAYER_INFO *) Mem_Calloc(1, sizeof(SW_LAYER_INFO), "_newlayer()"); + v->lyr[v->n_layers - 1]->id = v->n_layers - 1; return v->n_layers - 1; } @@ -709,9 +1512,20 @@ void SW_SIT_read(void) { case 41: v->type_soilDensityInput = atoi(inbuf); break; + case 42: + strcpy(v->site_swrc_name, inbuf); + v->site_swrc_type = encode_str2swrc(v->site_swrc_name); + break; + case 43: + strcpy(v->site_ptf_name, inbuf); + v->site_ptf_type = encode_str2ptf(v->site_ptf_name); + break; + case 44: + v->site_has_swrcp = itob(atoi(inbuf)); + break; default: - if (lineno > 41 + MAX_TRANSP_REGIONS) + if (lineno > 44 + MAX_TRANSP_REGIONS) break; /* skip extra lines */ if (MAX_TRANSP_REGIONS < v->n_transp_rgn) { @@ -735,31 +1549,136 @@ void SW_SIT_read(void) { CloseFile(&f); if (LT(v->percentRunoff, 0.) || GT(v->percentRunoff, 1.)) { - LogError(logfp, LOGFATAL, "%s : proportion of ponded surface water removed as daily" - "runoff = %f (value ranges between 0 and 1)\n", MyFileName, v->percentRunoff); + LogError( + logfp, + LOGFATAL, + "%s : proportion of ponded surface water removed as daily" + "runoff = %f (value ranges between 0 and 1)\n", + MyFileName, v->percentRunoff + ); } if (LT(v->percentRunon, 0.)) { - LogError(logfp, LOGFATAL, "%s : proportion of water that arrives at surface added " - "as daily runon = %f (value ranges between 0 and +inf)\n", MyFileName, v->percentRunon); + LogError( + logfp, + LOGFATAL, + "%s : proportion of water that arrives at surface added " + "as daily runon = %f (value ranges between 0 and +inf)\n", + MyFileName, v->percentRunon + ); } if (too_many_regions) { - LogError(logfp, LOGFATAL, "%s : Number of transpiration regions" - " exceeds maximum allowed (%d > %d)\n", MyFileName, v->n_transp_rgn, MAX_TRANSP_REGIONS); + LogError( + logfp, + LOGFATAL, + "%s : Number of transpiration regions" + " exceeds maximum allowed (%d > %d)\n", + MyFileName, v->n_transp_rgn, MAX_TRANSP_REGIONS + ); } /* check for any discontinuities (reversals) in the transpiration regions */ for (r = 1; r < v->n_transp_rgn; r++) { if (_TranspRgnBounds[r - 1] >= _TranspRgnBounds[r]) { - LogError(logfp, LOGFATAL, "%s : Discontinuity/reversal in transpiration regions.\n", SW_F_name(eSite)); + LogError( + logfp, + LOGFATAL, + "%s : Discontinuity/reversal in transpiration regions.\n", + SW_F_name(eSite) + ); + } + } +} + + + +/** Reads soil layers and soil properties from input file + + @note Previously, the function was static and named `_read_layers()`. +*/ +void SW_LYR_read(void) { + /* =================================================== */ + /* 5-Feb-2002 (cwb) removed dmin requirement in input file */ + + SW_SITE *v = &SW_Site; + FILE *f; + LyrIndex lyrno; + int x, k; + RealF dmin = 0.0, dmax, evco, trco_veg[NVEGTYPES], psand, pclay, soildensity, imperm, + soiltemp, f_gravel; + + /* note that Files.read() must be called prior to this. */ + MyFileName = SW_F_name(eLayers); + + f = OpenFile(MyFileName, "r"); + + while (GetALine(f, inbuf)) { + lyrno = _newlayer(); + + x = sscanf( + inbuf, + "%f %f %f %f %f %f %f %f %f %f %f %f", + &dmax, + &soildensity, + &f_gravel, + &evco, + &trco_veg[SW_GRASS], &trco_veg[SW_SHRUB], &trco_veg[SW_TREES], &trco_veg[SW_FORBS], + &psand, + &pclay, + &imperm, + &soiltemp + ); + /* Check that we have 12 values per layer */ + /* Adjust number if new variables are added */ + if (x != 12) { + CloseFile(&f); + LogError( + logfp, + LOGFATAL, + "%s : Incomplete record %d.\n", + MyFileName, lyrno + 1 + ); + } + + v->lyr[lyrno]->width = dmax - dmin; + + /* checks for valid values now carried out by `SW_SIT_init_run()` */ + + dmin = dmax; + v->lyr[lyrno]->fractionVolBulk_gravel = f_gravel; + v->lyr[lyrno]->soilDensityInput = soildensity; + v->lyr[lyrno]->evap_coeff = evco; + + ForEachVegType(k) + { + v->lyr[lyrno]->transp_coeff[k] = trco_veg[k]; + } + + v->lyr[lyrno]->fractionWeightMatric_sand = psand; + v->lyr[lyrno]->fractionWeightMatric_clay = pclay; + v->lyr[lyrno]->impermeability = imperm; + v->lyr[lyrno]->avgLyrTemp = soiltemp; + + if (lyrno >= MAX_LAYERS) { + CloseFile(&f); + LogError( + logfp, + LOGFATAL, + "%s : Too many layers specified (%d).\n" + "Maximum number of layers is %d\n", + MyFileName, lyrno + 1, MAX_LAYERS + ); } } - _read_layers(); + CloseFile(&f); } + + + /** @brief Creates soil layers based on function arguments (instead of reading them from an input file as _read_layers() does) @@ -940,6 +1859,90 @@ void derive_soilRegions(int nRegions, RealD *regionLowerBounds){ } } + + +/** Obtain soil water retention curve parameters from disk +*/ +void SW_SWRC_read(void) { + + /* Don't read values from disk if they will be estimated via a PTF */ + if (!SW_Site.site_has_swrcp) return; + + FILE *f; + LyrIndex lyrno = 0, k; + int x; + RealF tmp_swrcp[SWRC_PARAM_NMAX]; + + /* note that Files.read() must be called prior to this. */ + MyFileName = SW_F_name(eSWRCp); + + f = OpenFile(MyFileName, "r"); + + while (GetALine(f, inbuf)) { + x = sscanf( + inbuf, + "%f %f %f %f %f %f", + &tmp_swrcp[0], + &tmp_swrcp[1], + &tmp_swrcp[2], + &tmp_swrcp[3], + &tmp_swrcp[4], + &tmp_swrcp[5] + ); + + /* Note: `SW_SIT_init_run()` will call function to check for valid values */ + + /* Check that we have n = `SWRC_PARAM_NMAX` values per layer */ + if (x != SWRC_PARAM_NMAX) { + CloseFile(&f); + LogError( + logfp, + LOGFATAL, + "%s : Bad number of SWRC parameters %d -- must be %d.\n", + MyFileName, x, SWRC_PARAM_NMAX + ); + } + + /* Check that we are within `SW_Site.n_layers` */ + if (lyrno >= SW_Site.n_layers) { + CloseFile(&f); + LogError( + logfp, + LOGFATAL, + "%s : Number of layers with SWRC parameters (%d) " + "must match number of soil layers (%d)\n", + MyFileName, lyrno + 1, SW_Site.n_layers + ); + } + + /* Copy values into structure */ + for (k = 0; k < SWRC_PARAM_NMAX; k++) { + SW_Site.lyr[lyrno]->swrcp[k] = tmp_swrcp[k]; + } + + lyrno++; + } + + CloseFile(&f); +} + + +/** + @brief Derive and check soil properties from inputs + + Bulk refers to the whole soil, + i.e., including the rock/gravel component (coarse fragments), + whereas matric refers to the < 2 mm fraction. + + Internally, SOILWAT2 calculates based on bulk soil, i.e., the whole soil. + However, sand and clay inputs are expected to represent the soil matric, + i.e., the < 2 mm fraction. + + sand + clay + silt must equal one. + Fraction of silt is calculated: 1 - (sand + clay). + + @sideeffect Values stored in global variable `SW_Site`. +*/ void SW_SIT_init_run(void) { /* =================================================== */ /* potentially this routine can be called whether the @@ -953,12 +1956,8 @@ void SW_SIT_init_run(void) { SW_LAYER_INFO *lyr; LyrIndex s, r, curregion; int k, flagswpcrit = 0; - Bool fail = swFALSE; RealD - fval = 0, - evsum = 0., trsum_veg[NVEGTYPES] = {0.}, - swcmin_help1, swcmin_help2, tmp; - const char *errtype = "\0"; + evsum = 0., trsum_veg[NVEGTYPES] = {0.}, tmp; #ifdef SWDEBUG int debug = 0; @@ -975,93 +1974,40 @@ void SW_SIT_init_run(void) { add_deepdrain_layer(); + /* Check compatibility between selected SWRC and PTF */ + if (!sp->site_has_swrcp) { + if (!check_SWRC_vs_PTF(sp->site_swrc_name, sp->site_ptf_name)) { + LogError( + logfp, + LOGFATAL, + "Selected PTF '%s' is incompatible with selected SWRC '%s'\n", + sp->site_ptf_name, + sp->site_swrc_name + ); + } + } + + /* Loop over soil layers check variables and calculate parameters */ ForEachSoilLayer(s) { lyr = sp->lyr[s]; - /* Check validity of soil variables: - previously, checked by code in `_read_layers()`, - erroneously skipped by `set_soillayers()`, - and checked by code in rSOILWAT2's `onSet_SW_LYR()` + /* Copy site-level SWRC/PTF information to each layer: + We currently allow specifying one SWRC/PTF for a site for all layers; + remove in the future if we allow SWRC/PTF to vary by soil layer */ - if (LE(lyr->width, 0.)) { - fail = swTRUE; - fval = lyr->width; - errtype = Str_Dup("layer width"); - - } else if ( - LT(lyr->soilDensityInput, 0.) || - GT(lyr->soilDensityInput, 2.65) - ) { - fail = swTRUE; - fval = lyr->soilDensityInput; - errtype = Str_Dup("soil density"); - - } else if ( - LT(lyr->fractionVolBulk_gravel, 0.) || - GE(lyr->fractionVolBulk_gravel, 1.) - ) { - fail = swTRUE; - fval = lyr->fractionVolBulk_gravel; - errtype = Str_Dup("gravel content"); - - } else if ( - LE(lyr->fractionWeightMatric_sand, 0.) || - GE(lyr->fractionWeightMatric_sand, 1.) - ) { - fail = swTRUE; - fval = lyr->fractionWeightMatric_sand; - errtype = Str_Dup("sand proportion"); - - } else if ( - LE(lyr->fractionWeightMatric_clay, 0.) || - GE(lyr->fractionWeightMatric_clay, 1.) - ) { - fail = swTRUE; - fval = lyr->fractionWeightMatric_clay; - errtype = Str_Dup("clay proportion"); - - } else if ( - GE(lyr->fractionWeightMatric_sand + lyr->fractionWeightMatric_clay, 1.) - ) { - fail = swTRUE; - fval = lyr->fractionWeightMatric_sand + lyr->fractionWeightMatric_clay; - errtype = Str_Dup("sand+clay proportion"); - - } else if ( - LT(lyr->impermeability, 0.) || - GT(lyr->impermeability, 1.) - ) { - fail = swTRUE; - fval = lyr->impermeability; - errtype = Str_Dup("impermeability"); + lyr->swrc_type = sp->site_swrc_type; + lyr->ptf_type = sp->site_ptf_type; - } else if ( - LT(lyr->evap_coeff, 0.) || - GT(lyr->evap_coeff, 1.) - ) { - fail = swTRUE; - fval = lyr->evap_coeff; - errtype = Str_Dup("bare-soil evaporation coefficient"); - - } else { - ForEachVegType(k) { - if (LT(lyr->transp_coeff[k], 0.) || GT(lyr->transp_coeff[k], 1.)) { - fail = swTRUE; - fval = lyr->transp_coeff[k]; - errtype = Str_Dup("transpiration coefficient"); - break; - } - } - } - if (fail) { + /* Check soil properties for valid values */ + if (!SW_check_soil_properties(lyr)) { LogError( logfp, LOGFATAL, - "Invalid %s (%5.4f) in layer %d.\n", - errtype, fval, s + 1 + "Invalid soil properties in layer %d.\n", + lyr->id + 1 ); } @@ -1100,26 +2046,37 @@ void SW_SIT_init_run(void) { - /* Calculate pedotransfer function parameters */ - water_eqn( - lyr->fractionVolBulk_gravel, - lyr->fractionWeightMatric_sand, - lyr->fractionWeightMatric_clay, - s - ); - /* Calculate SWC at field capacity and at wilting point */ - lyr->swcBulk_fieldcap = lyr->width * SW_SWPmatric2VWCBulk( - lyr->fractionVolBulk_gravel, - 0.333, - s - ); + if (!sp->site_has_swrcp) { + /* 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()` + */ + SWRC_PTF_estimate_parameters( + lyr->ptf_type, + lyr->swrcp, + lyr->fractionWeightMatric_sand, + lyr->fractionWeightMatric_clay, + lyr->fractionVolBulk_gravel, + lyr->soilBulk_density + ); + } - lyr->swcBulk_wiltpt = lyr->width * SW_SWPmatric2VWCBulk( - lyr->fractionVolBulk_gravel, - 15, - s - ); + /* Check parameters of selected SWRC */ + if (!SWRC_check_parameters(lyr->swrc_type, lyr->swrcp)) { + LogError( + logfp, + LOGFATAL, + "Checks of parameters for SWRC '%s' in layer %d failed.", + swrc2str[lyr->swrc_type], + lyr->id + 1 + ); + } + + /* Calculate SWC at field capacity and at wilting point */ + lyr->swcBulk_fieldcap = SW_SWRC_SWPtoSWC(0.333, lyr); + lyr->swcBulk_wiltpt = SW_SWRC_SWPtoSWC(15., lyr); /* Calculate lower SWC limit of bare-soil evaporation as `max(0.5 * wiltpt, SWC@hygroscopic)` @@ -1132,69 +2089,43 @@ void SW_SIT_init_run(void) { */ lyr->swcBulk_halfwiltpt = fmax( 0.5 * lyr->swcBulk_wiltpt, - lyr->width * SW_SWPmatric2VWCBulk(lyr->fractionVolBulk_gravel, 100., s) + SW_SWRC_SWPtoSWC(100., lyr) ); - /* Compute swc wet and dry limits and init value */ - if (LT(_SWCMinVal, 0.0)) { - /* input: estimate mininum SWC */ - - /* residual SWC of Rawls & Brakensiek (1985) */ - swcmin_help1 = SW_VWCBulkRes( - lyr->fractionVolBulk_gravel, - lyr->fractionWeightMatric_sand, - lyr->fractionWeightMatric_clay, - lyr->swcBulk_saturated / ((1. - lyr->fractionVolBulk_gravel) * lyr->width) - ); - - /* Lower limit for swc_min - Notes: - - used in case the equation for residual SWC doesn't work or - produces unrealistic small values - - currently, -30 MPa - (based on limited test runs across western US including hot deserts) - lower than "air-dry" = hygroscopic point (-10. MPa; Porporato et al. 2001) - not as extreme as "oven-dry" (-1000. MPa; Fredlund et al. 1994) - */ - swcmin_help2 = SW_SWPmatric2VWCBulk(lyr->fractionVolBulk_gravel, 300., s); - - // if `SW_VWCBulkRes()` returns SW_MISSING then use `swcmin_help2` - if (missing(swcmin_help1)) { - lyr->swcBulk_min = swcmin_help2; - - } else { - lyr->swcBulk_min = fmax(swcmin_help1, swcmin_help2); - } - - } else if (GE(_SWCMinVal, 1.0)) { - /* input: fixed SWP value as minimum SWC; unit(_SWCMinVal) == -bar */ - lyr->swcBulk_min = SW_SWPmatric2VWCBulk( - lyr->fractionVolBulk_gravel, - _SWCMinVal, - s - ); - - } else { - /* input: fixed VWC value as minimum SWC; unit(_SWCMinVal) == cm/cm */ - lyr->swcBulk_min = _SWCMinVal; - } + /* Extract or estimate additional properties */ + lyr->swcBulk_saturated = SW_swcBulk_saturated( + lyr->swrc_type, + lyr->swrcp, + lyr->fractionVolBulk_gravel, + lyr->width, + lyr->ptf_type, + lyr->fractionWeightMatric_sand, + lyr->fractionWeightMatric_clay + ); - /* Convert VWC to SWC */ - lyr->swcBulk_min *= lyr->width; + lyr->swcBulk_min = SW_swcBulk_minimum( + lyr->swrc_type, + lyr->swrcp, + lyr->fractionVolBulk_gravel, + lyr->width, + lyr->ptf_type, + _SWCMinVal, + lyr->fractionWeightMatric_sand, + lyr->fractionWeightMatric_clay, + lyr->swcBulk_saturated + ); /* Calculate wet limit of SWC for what inputs defined as wet */ - lyr->swcBulk_wet = lyr->width * (GE(_SWCWetVal, 1.0) ? - SW_SWPmatric2VWCBulk(lyr->fractionVolBulk_gravel, _SWCWetVal, s) : - _SWCWetVal - ); + lyr->swcBulk_wet = GE(_SWCWetVal, 1.0) ? + SW_SWRC_SWPtoSWC(_SWCWetVal, lyr) : + _SWCWetVal * lyr->width; /* Calculate initial SWC based on inputs */ - lyr->swcBulk_init = lyr->width * (GE(_SWCInitVal, 1.0) ? - SW_SWPmatric2VWCBulk(lyr->fractionVolBulk_gravel, _SWCInitVal, s) : - _SWCInitVal - ); + lyr->swcBulk_init = GE(_SWCInitVal, 1.0) ? + SW_SWRC_SWPtoSWC(_SWCInitVal, lyr) : + _SWCInitVal * lyr->width; /* test validity of values */ @@ -1204,7 +2135,7 @@ void SW_SIT_init_run(void) { "%s : Layer %d\n" " calculated `swcBulk_init` (%.4f cm) <= `swcBulk_min` (%.4f cm).\n" " Recheck parameters and try again.\n", - MyFileName, s + 1, lyr->swcBulk_init, lyr->swcBulk_min + MyFileName, lyr->id + 1, lyr->swcBulk_init, lyr->swcBulk_min ); } @@ -1214,7 +2145,7 @@ void SW_SIT_init_run(void) { "%s : Layer %d\n" " calculated `swcBulk_wiltpt` (%.4f cm) <= `swcBulk_min` (%.4f cm).\n" " Recheck parameters and try again.\n", - MyFileName, s + 1, lyr->swcBulk_wiltpt, lyr->swcBulk_min + MyFileName, lyr->id + 1, lyr->swcBulk_wiltpt, lyr->swcBulk_min ); } @@ -1225,11 +2156,11 @@ void SW_SIT_init_run(void) { " calculated `swcBulk_halfwiltpt` (%.4f cm / %.2f MPa)\n" " <= `swcBulk_min` (%.4f cm / %.2f MPa).\n" " `swcBulk_halfwiltpt` was set to `swcBulk_min`.\n", - MyFileName, s + 1, + MyFileName, lyr->id + 1, lyr->swcBulk_halfwiltpt, - -0.1 * SW_SWCbulk2SWPmatric(lyr->fractionVolBulk_gravel, lyr->swcBulk_halfwiltpt, s), + -0.1 * SW_SWRC_SWCtoSWP(lyr->swcBulk_halfwiltpt, lyr), lyr->swcBulk_min, - -0.1 * SW_SWCbulk2SWPmatric(lyr->fractionVolBulk_gravel, lyr->swcBulk_min, s) + -0.1 * SW_SWRC_SWCtoSWP(lyr->swcBulk_min, lyr) ); lyr->swcBulk_halfwiltpt = lyr->swcBulk_min; @@ -1241,7 +2172,7 @@ void SW_SIT_init_run(void) { "%s : Layer %d\n" " calculated `swcBulk_wet` (%.4f cm) <= `swcBulk_min` (%.4f cm).\n" " Recheck parameters and try again.\n", - MyFileName, s + 1, lyr->swcBulk_wet, lyr->swcBulk_min + MyFileName, lyr->id + 1, lyr->swcBulk_wet, lyr->swcBulk_min ); } @@ -1252,18 +2183,14 @@ void SW_SIT_init_run(void) { "L[%d] swcmin=%f = swpmin=%f\n", s, lyr->swcBulk_min, - SW_SWCbulk2SWPmatric(lyr->fractionVolBulk_gravel, lyr->swcBulk_min, s) + SW_SWRC_SWCtoSWP(lyr->swcBulk_min, lyr) ); swprintf( "L[%d] SWC(HalfWiltpt)=%f = swp(hw)=%f\n", s, lyr->swcBulk_halfwiltpt, - SW_SWCbulk2SWPmatric( - lyr->fractionVolBulk_gravel, - lyr->swcBulk_halfwiltpt, - s - ) + SW_SWRC_SWCtoSWP(lyr->swcBulk_halfwiltpt, lyr) ); } #endif @@ -1276,10 +2203,9 @@ void SW_SIT_init_run(void) { { trsum_veg[k] += lyr->transp_coeff[k]; /* calculate soil water content at SWPcrit for each vegetation type */ - lyr->swcBulk_atSWPcrit[k] = lyr->width * SW_SWPmatric2VWCBulk( - lyr->fractionVolBulk_gravel, + lyr->swcBulk_atSWPcrit[k] = SW_SWRC_SWPtoSWC( SW_VegProd.veg[k].SWPcrit, - s + lyr ); if (LT(lyr->swcBulk_atSWPcrit[k], lyr->swcBulk_min)) { @@ -1288,11 +2214,7 @@ void SW_SIT_init_run(void) { // lower SWcrit [-bar] to SWP-equivalent of swBulk_min tmp = fmin( SW_VegProd.veg[k].SWPcrit, - SW_SWCbulk2SWPmatric( - lyr->fractionVolBulk_gravel, - lyr->swcBulk_min, - s - ) + SW_SWRC_SWCtoSWP(lyr->swcBulk_min, lyr) ); LogError( @@ -1302,11 +2224,11 @@ void SW_SIT_init_run(void) { " <= `swcBulk_min` (%.4f cm / %.4f MPa).\n" " `SWcrit` adjusted to %.4f MPa " "(and swcBulk_atSWPcrit in every layer will be re-calculated).\n", - MyFileName, s + 1, k + 1, + MyFileName, lyr->id + 1, k + 1, lyr->swcBulk_atSWPcrit[k], -0.1 * SW_VegProd.veg[k].SWPcrit, lyr->swcBulk_min, - -0.1 * SW_SWCbulk2SWPmatric(lyr->fractionVolBulk_gravel, lyr->swcBulk_min, s), + -0.1 * SW_SWRC_SWCtoSWP(lyr->swcBulk_min, lyr), -0.1 * tmp ); @@ -1359,10 +2281,9 @@ void SW_SIT_init_run(void) { ForEachVegType(k) { /* calculate soil water content at adjusted SWPcrit */ - lyr->swcBulk_atSWPcrit[k] = lyr->width * SW_SWPmatric2VWCBulk( - lyr->fractionVolBulk_gravel, + lyr->swcBulk_atSWPcrit[k] = SW_SWRC_SWPtoSWC( SW_VegProd.veg[k].SWPcrit, - s + lyr ); if (LT(lyr->swcBulk_atSWPcrit[k], lyr->swcBulk_min)) { @@ -1372,7 +2293,7 @@ void SW_SIT_init_run(void) { " calculated `swcBulk_atSWPcrit` (%.4f cm)\n" " <= `swcBulk_min` (%.4f cm).\n" " even with adjusted `SWcrit` (%.4f MPa).\n", - MyFileName, s + 1, k + 1, + MyFileName, lyr->id + 1, k + 1, lyr->swcBulk_atSWPcrit[k], lyr->swcBulk_min, -0.1 * SW_VegProd.veg[k].SWPcrit @@ -1392,16 +2313,26 @@ void SW_SIT_init_run(void) { /* normalize the evap and transp coefficients separately * to avoid obfuscation in the above loop */ if (!EQ_w_tol(evsum, 1.0, 1e-4)) { // inputs are not more precise than at most 3-4 digits - LogError(logfp, LOGWARN, - "%s : Evaporation coefficients were normalized:\n" \ - "\tSum of coefficients was %.4f, but must be 1.0. " \ - "New coefficients are:", MyFileName, evsum); + LogError( + logfp, + LOGWARN, + "%s : Evaporation coefficients were normalized:\n" + "\tSum of coefficients was %.4f, but must be 1.0. " + "New coefficients are:", + MyFileName, + evsum + ); ForEachEvapLayer(s) { SW_Site.lyr[s]->evap_coeff /= evsum; - LogError(logfp, LOGNOTE, " Layer %2d : %.4f", - s + 1, SW_Site.lyr[s]->evap_coeff); + LogError( + logfp, + LOGNOTE, + " Layer %2d : %.4f", + SW_Site.lyr[s]->id + 1, + SW_Site.lyr[s]->evap_coeff + ); } LogError(logfp, LOGQUIET, ""); @@ -1409,19 +2340,29 @@ void SW_SIT_init_run(void) { ForEachVegType(k) { - if (!EQ_w_tol(trsum_veg[k], 1.0, 1e-4)) { // inputs are not more precise than at most 3-4 digits - LogError(logfp, LOGWARN, - "%s : Transpiration coefficients were normalized for %s:\n" \ - "\tSum of coefficients was %.4f, but must be 1.0. " \ - "New coefficients are:", MyFileName, key2veg[k], trsum_veg[k]); + // inputs are not more precise than at most 3-4 digits + if (!EQ_w_tol(trsum_veg[k], 1.0, 1e-4)) { + LogError( + logfp, + LOGWARN, + "%s : Transpiration coefficients were normalized for %s:\n" + "\tSum of coefficients was %.4f, but must be 1.0. " + "New coefficients are:", + MyFileName, key2veg[k], trsum_veg[k] + ); ForEachSoilLayer(s) { if (GT(SW_Site.lyr[s]->transp_coeff[k], 0.)) { SW_Site.lyr[s]->transp_coeff[k] /= trsum_veg[k]; - LogError(logfp, LOGNOTE, " Layer %2d : %.4f", - s + 1, SW_Site.lyr[s]->transp_coeff[k]); + LogError( + logfp, + LOGNOTE, + " Layer %2d : %.4f", + SW_Site.lyr[s]->id + 1, + SW_Site.lyr[s]->transp_coeff[k] + ); } } @@ -1437,17 +2378,23 @@ void SW_SIT_init_run(void) { if (too_many_RGR) { // because we will use loops such `for (i = 0; i <= nRgr + 1; i++)` - LogError(logfp, LOGWARN, - "\nSOIL_TEMP FUNCTION ERROR: the number of regressions is > the "\ - "maximum number of regressions. resetting max depth, deltaX, nRgr "\ - "values to 180, 15, & 11 respectively\n"); + LogError( + logfp, + LOGWARN, + "\nSOIL_TEMP FUNCTION ERROR: the number of regressions is > the " + "maximum number of regressions. resetting max depth, deltaX, nRgr " + "values to 180, 15, & 11 respectively\n" + ); } else { // because we don't deal with partial layers - LogError(logfp, LOGWARN, - "\nSOIL_TEMP FUNCTION ERROR: max depth is not evenly divisible by "\ - "deltaX (ie the remainder != 0). resetting max depth, deltaX, nRgr "\ - "values to 180, 15, & 11 respectively\n"); + LogError( + logfp, + LOGWARN, + "\nSOIL_TEMP FUNCTION ERROR: max depth is not evenly divisible by " + "deltaX (ie the remainder != 0). resetting max depth, deltaX, nRgr " + "values to 180, 15, & 11 respectively\n" + ); } // resets it to the default values @@ -1607,22 +2554,72 @@ void _echo_inputs(void) { ForEachSoilLayer(i) { - LogError(logfp, LOGNOTE, " %3d %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f\n", i + 1, - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_fieldcap, i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_wiltpt, i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_atSWPcrit[SW_FORBS], i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_atSWPcrit[SW_TREES], i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_atSWPcrit[SW_SHRUB], i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_atSWPcrit[SW_GRASS], i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_wet, i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_min, i), - SW_SWCbulk2SWPmatric(s->lyr[i]->fractionVolBulk_gravel, s->lyr[i]->swcBulk_init, i)); + LogError( + logfp, + LOGNOTE, + " %3d %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f %15.4f\n", + i + 1, + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_fieldcap, s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_wiltpt, s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_atSWPcrit[SW_FORBS], s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_atSWPcrit[SW_TREES], s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_atSWPcrit[SW_SHRUB], s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_atSWPcrit[SW_SHRUB], s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_atSWPcrit[SW_GRASS], s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_min, s->lyr[i]), + SW_SWRC_SWCtoSWP(s->lyr[i]->swcBulk_init, s->lyr[i]) + ); + } + LogError( + logfp, + LOGNOTE, + "\nSoil Water Retention Curve:\n---------------------------\n" + ); + LogError( + logfp, + LOGNOTE, + " SWRC type: %d (%s)\n", + s->site_swrc_type, + s->site_swrc_name + ); + LogError( + logfp, + LOGNOTE, + " PTF type: %d (%s)\n", + s->site_ptf_type, + s->site_ptf_name + ); + + LogError( + logfp, + LOGNOTE, + " Lyr Param1 Param2 Param3 Param4 Param5 Param6\n" + ); + ForEachSoilLayer(i) + { + LogError( + logfp, + LOGNOTE, + " %3d%11.4f%11.4f%11.4f%11.4f%11.4f%11.4f\n", + i + 1, + s->lyr[i]->swrcp[0], + s->lyr[i]->swrcp[1], + s->lyr[i]->swrcp[2], + s->lyr[i]->swrcp[3], + s->lyr[i]->swrcp[4], + s->lyr[i]->swrcp[5] + ); } - LogError(logfp, LOGNOTE, "\n------------ End of Site Parameters ------------------\n"); - //fflush(logfp); + + LogError( + logfp, + LOGNOTE, + "\n------------ End of Site Parameters ------------------\n" + ); + //fflush(logfp); } #ifdef DEBUG_MEM diff --git a/SW_Site.h b/SW_Site.h index 6008c0f04..8f1f5de39 100644 --- a/SW_Site.h +++ b/SW_Site.h @@ -56,12 +56,118 @@ extern "C" { #define SW_MATRIC 0 #define SW_BULK 1 + +/** + @defgroup swrc_ptf Soil Water Retention Curves + + __Soil Water Retention Curves (SWRCs) -- Pedotransfer functions (PTFs)__ + + __Overview__ + + Historically (before v7.0.0), `SOILWAT2` utilized a hard-coded SWRC + by Campbell 1974 (\cite Campbell1974) and estimated SWRC parameters at + run-time from soil texture using PTFs by Cosby et al. 1984 (\cite Cosby1984). + This behavior can be reproduced with "Campbell1974" and "Cosby1984AndOthers" + (see input file `siteparam.in`). + + Now, users of `SOILWAT2` can choose from a range of implemented SWRCs + (see input file `siteparam.in`); + SWRC parameters can be estimated at run-time from soil properties by + selecting a matching PTF (see input file `siteparam.in`) or, + alternatively (`has_swrcp`), provide adequate SWRC parameter values + (see input file `swrc_params.in`). + Please note that `rSOILWAT2` may provide additional PTF functionality. + + + __Approach__ + + -# User selections of SWRC and PTF are read in from input file `siteparam.in` + by `SW_SIT_read()` and, if `has_swrcp`, SWRC parameters are read from + input file `swrc_params.in` by `SW_SWRC_read()`. + + -# `SW_SIT_init_run()` + - if not `has_swrcp` + - calls `check_SWRC_vs_PTF()` to check that selected SWRC and PTF are + compatible + - calls `SWRC_PTF_estimate_parameters()` to estimate + SWRC parameter values from soil properties based on selected PTF + - calls `SWRC_check_parameters()` to check that SWRC parameter values + are resonable for the selected SWRC. + + -# `SW_SWRC_SWCtoSWP()` and `SW_SWRC_SWPtoSWC()` are used during simulation + runs to convert between soil water content and soil water potential. + + -# These high-level "wrapper" functions hide details of any specific SWRC/PTF + implementations and are used by SOILWAT2 simulation code. + Thus, most of SOILWAT2 is "unaware" about the selected SWRC/PTF + and how to interpret SWRC parameters. Instead, these "wrapper" functions + know how to call the appropriate SWRC and/or PTF specific functions + which implement SWRC and/or PTF specific details. + + + __Steps to implement a new SWRC "XXX" and corresponding PTF "YYY"__ + + -# Update #N_SWRCs and #N_PTFs + + -# Add new names to #swrc2str and #ptf2str and + add corresponding macros of indices + + -# Update input files `siteparam.in` and `swrc_params.in` + + -# Implement new XXX/YYY-specific functions + - `SWRC_check_parameters_for_XXX()` to validate parameter values, + e.g., `SWRC_check_parameters_for_Campbell1974()` + - `SWRC_PTF_YYY_for_XXX()` to estimate parameter values (if implemented), + e.g., `SWRC_PTF_Cosby1984_for_Campbell1974()` + - `SWRC_SWCtoSWP_XXX()` to translate moisture content to water potential, + e.g., `SWRC_SWCtoSWP_Campbell1974()` + - `SWRC_SWPtoSWC_XXX()` to translate water potential to moisture content, + e.g., `SWRC_SWPtoSWC_Campbell1974()` + + -# Update "wrapper" functions that select and call XXX/YYY-specific functions + and/or parameters + - `check_SWRC_vs_PTF()` + - `SWRC_PTF_estimate_parameters()` (if PTF is implemented) + - `SWRC_check_parameters()` + - `SWRC_SWCtoSWP()` + - `SWRC_SWPtoSWC()` + - `SW_swcBulk_minimum()` + - `SW_swcBulk_saturated()` + + -# Expand existing unit tests and add new unit tests + to utilize new XXX/YYY functions +*/ + +#define SWRC_PARAM_NMAX 6 /**< Maximal number of SWRC parameters implemented */ +#define N_SWRCs 3 /**< Number of SWRCs implemented by SOILWAT2 */ +#define N_PTFs 2 /**< Number of PTFs implemented by SOILWAT2 */ + +// Indices of #swrc2str (for code readability) +#define sw_Campbell1974 0 +#define sw_vanGenuchten1980 1 +#define sw_FXW 2 + +// Indices of #swrc2str (for code readability) +#define sw_Cosby1984AndOthers 0 +#define sw_Cosby1984 1 + + +#define FXW_h0 6.3e6 /**< Pressure head at zero water content [cm] of FWX SWRC */ +#define FXW_hr 1500. /**< Pressure head at residual water content [cm] of FXW SWRC */ + + +/* =================================================== */ +/* TYPEDEFS */ +/* --------------------------------------------------- */ + typedef unsigned int LyrIndex; typedef struct { /* 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 */ + LyrIndex id; /**< Number of soil layer: 1 = most shallow, 2 = second shallowest, etc. up to ::MAX_LAYERS */ + RealD /* Inputs */ width, /* width of the soil layer (cm) */ @@ -86,17 +192,19 @@ typedef struct { swcBulk_atSWPcrit[NVEGTYPES], /* SWC corresponding to critical SWP for transpiration */ /* Saxton et al. 2006 */ - swcBulk_saturated, /* saturated bulk SWC [cm] */ - 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] */ - Saxton2006_lambda, /* Slope of logarithmic tension-moisture curve */ - - /* Cosby et al. (1984): SOILWAT2's soil water retention curve */ - thetasMatric, /* saturated matric SWC [cm] */ - psisMatric, /* saturated matric SWP [cm] */ - bMatric, /* slope of the logarithmic retention curve */ - binverseMatric; /* inverse of bMatric */ + swcBulk_saturated; /* saturated bulk SWC [cm] */ + // 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] */ + //Saxton2006_lambda; /* Slope of logarithmic tension-moisture curve */ + + + /* Soil water retention curve (SWRC) */ + unsigned int + swrc_type, /**< Type of SWRC (see #swrc2str) */ + ptf_type; /**< Type of PTF (see #ptf2str) */ + RealD swrcp[SWRC_PARAM_NMAX]; /**< Parameters of SWRC: parameter interpretation varies with selected SWRC, see `SWRC_check_parameters()` */ LyrIndex my_transp_rgn[NVEGTYPES]; /* which transp zones from Site am I in? */ } SW_LAYER_INFO; @@ -153,6 +261,17 @@ typedef struct { SW_LAYER_INFO **lyr; /* one struct per soil layer pointed to by */ /* a dynamically allocated block of pointers */ + /* Soil water retention curve (SWRC), see `SW_LAYER_INFO` */ + unsigned int + site_swrc_type, + site_ptf_type; + + char + site_swrc_name[64], + site_ptf_name[64]; + + Bool site_has_swrcp; /**< Are `swrcp` already (TRUE) or not yet estimated (FALSE)? */ + } SW_SITE; @@ -164,12 +283,71 @@ extern SW_SITE SW_Site; extern LyrIndex _TranspRgnBounds[MAX_TRANSP_REGIONS]; extern RealD _SWCInitVal, _SWCWetVal, _SWCMinVal; +extern char const *swrc2str[]; +extern char const *ptf2str[]; /* =================================================== */ /* Global Function Declarations */ /* --------------------------------------------------- */ -void water_eqn(RealD fractionGravel, RealD sand, RealD clay, LyrIndex n); + +unsigned int encode_str2swrc(char *swrc_name); +unsigned int encode_str2ptf(char *ptf_name); + +void SWRC_PTF_estimate_parameters( + unsigned int ptf_type, + double *swrcp, + double sand, + double clay, + double gravel, + double bdensity +); +void SWRC_PTF_Cosby1984_for_Campbell1974( + double *swrcp, + double sand, + double clay +); + + +Bool check_SWRC_vs_PTF(char *swrc_name, char *ptf_name); +Bool SWRC_check_parameters(unsigned int swrc_type, double *swrcp); +Bool SWRC_check_parameters_for_Campbell1974(double *swrcp); +Bool SWRC_check_parameters_for_vanGenuchten1980(double *swrcp); +Bool SWRC_check_parameters_for_FXW(double *swrcp); + +double SW_swcBulk_saturated( + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + unsigned int ptf_type, + double sand, + double clay +); +double SW_swcBulk_minimum( + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + unsigned int ptf_type, + double ui_sm_min, + double sand, + double clay, + double swcBulk_sat +); +void PTF_Saxton2006( + double *theta_sat, + double sand, + double clay +); +void PTF_RawlsBrakensiek1985( + double *theta_min, + double sand, + double clay, + double porosity +); + + RealD calculate_soilBulkDensity(RealD matricDensity, RealD fractionGravel); RealD calculate_soilMatricDensity(RealD bulkDensity, RealD fractionGravel); LyrIndex nlayers_bsevap(void); @@ -183,6 +361,8 @@ void SW_SIT_init_run(void); void _echo_inputs(void); /* these used to be in Layers */ +void SW_LYR_read(void); +void SW_SWRC_read(void); void SW_SIT_clear_layers(void); LyrIndex _newlayer(void); void add_deepdrain_layer(void); diff --git a/SW_SoilWater.c b/SW_SoilWater.c index 0467818a9..125d093b5 100644 --- a/SW_SoilWater.c +++ b/SW_SoilWater.c @@ -113,6 +113,237 @@ static void _reset_swc(void) { } +/** + @brief Convert soil water tension to volumetric water content using + FXW Soil Water Retention Curve + \cite fredlund1994CGJa, \cite wang2018wrr + + @note + This is the internal, bare-bones version of FXW, + use `SWRC_SWPtoSWC_FXW()` instead. + + @param[in] phi Pressure head (water tension) [cm of H20] + @param[in] *swrcp Vector of SWRC parameters + + @return Volumetric soil water content of the matric soil [cm / cm] +*/ +static double FXW_phi_to_theta( + double phi, + double *swrcp +) { + double tmp, res, S_e, C_f; + + if (GE(phi, FXW_h0)) { + res = 0.; + + } else { + // Relative saturation function: Rudiyanto et al. 2021: eq. 3 + tmp = pow(fabs(swrcp[1] * phi), swrcp[2]); // `pow()` because phi >= 0 + S_e = powe(log(swE + tmp), -swrcp[3]); + + // Correction factor: Rudiyanto et al. 2021: eq. 4 (with hr = 1500 cm) + // log(1. + 6.3e6 / 1500.) = 8.34307787116938; + C_f = 1. - log(1. + phi / FXW_hr) / 8.34307787116938; + + // Calculate matric theta [cm / cm] which is within [0, theta_sat] + // Rudiyanto et al. 2021: eq. 1 + res = swrcp[0] * S_e * C_f; + } + + // matric theta [cm / cm] + return res; +} + + +/** + \brief Interpolate, Truncate, Project (ITP) method (\cite oliveira2021ATMS) to solve FXW at a given water content + + We estimate `phi` by finding a root of + `f(phi) = theta - FXW(phi) = 0` + with + `a < b && f(a) < 0 && f(b) > 0` + + FXW SWRC (\cite fredlund1994CGJa, \cite wang2018wrr) is defined as + `theta = theta_sat * S_e(phi) * C_f(phi)` + + where + - `phi` is water tension (pressure head) [cm of H20], + - `theta` is volumetric water content [cm / cm], + - `theta_sat` is saturated water content [cm / cm], + - `S_e` is the relative saturation function [-], and + - `C_f` is a correction factor [-]. + + We choose starting values to solve `f(phi)` + - `a = 0` for which `f(a) = theta - theta_sat < 0 (for theta < theta_sat)` + - `b = FXW_h0` for which `f(b) = theta - 0 > 0 (for theta > 0)` + + [Additional info on the ITP method can be found here](https://en.wikipedia.org/wiki/ITP_method) + + + \param[in] theta Volumetric soil water content of the matric soil [cm / cm] + \param[in] *swrcp Vector of FXW SWRC parameters + + \return Water tension [bar] corresponding to theta +*/ + +static double itp_FXW_for_phi(double theta, double *swrcp) { + double + tol2 = 2e-9, // tolerance convergence + a = 0., // lower bound of bracket: a < b && f(a) < 0 + b = FXW_h0, // upper bound of bracket: b > a && f(b) > 0 + diff_ba = FXW_h0, diff_hf, + x_f, x_t, x_itp, + y_a, y_b, y_itp, + k1, k2, + x_half, + r, + delta, + sigma, + phi; + int + j = 0, n0, n_max; + + #ifdef SWDEBUG + short debug = 0; + #endif + + + // Evaluate at starting bracket (and checks) + y_a = theta - FXW_phi_to_theta(a, swrcp); + y_b = theta - FXW_phi_to_theta(b, swrcp); + + + // Set hyper-parameters + k1 = 3.174603e-08; // 0 < k1 = 0.2 / (b - a) < inf + /* + k1 = 2. converges in about 31-33 iterations + k1 = 0.2 converges in 28-30 iterations + k1 = 2.e-2 converges in 25-27 iterations + k1 = 2.e-3 converges in 22-24 iterations + k1 = 2.e-4 converges in about 20 iterations but fails for some + k1 = 2.e-6 converges in about 14 iterations or near n_max or fails + k1 = 0.2 / (b - a) = 3.174603e-08 (suggested value) converges + in about 8 iterations (for very dry soils) or near n_max = 53 or fails + */ + k1 = 2.e-3; + k2 = 2.; // 1 <= k2 < (3 + sqrt(5))/2 + n0 = 1; // 0 < n0 < inf + + // Upper limit of iterations before convergence + n_max = (int) ceil(log2(diff_ba / tol2)) + n0; + + #ifdef SWDEBUG + if (debug) { + swprintf( + "\nitp_FXW_for_phi(theta=%f): init: " + "f(a)=%f, f(b)=%f, n_max=%d, k1=%f\n", + theta, y_a, y_b, n_max, k1 + ); + } + #endif + + + // Iteration + do { + x_half = (a + b) / 2.; + + // Calculate parameters + r = (tol2 * exp2(n_max - j) - diff_ba) / 2.; + delta = k1 * pow(diff_ba, k2); + + // Interpolation + x_f = (y_b * a - y_a * b) / (y_b - y_a); + + // Truncation + diff_hf = x_half - x_f; + if (ZRO(diff_hf)) { + sigma = 0.; // `copysign()` returns -1 or +1 but never 0 + } else { + sigma = copysign(1., diff_hf); + } + + x_t = (delta <= fabs(diff_hf)) ? (x_f + sigma * delta) : x_half; + + // Projection + x_itp = (fabs(x_t - x_half) <= r) ? x_t : (x_half - sigma * r); + + // Evaluate function + y_itp = theta - FXW_phi_to_theta(x_itp, swrcp); + + #ifdef SWDEBUG + if (debug) { + swprintf( + "j=%d:" + "\ta=%f, b=%f, y_a=%f, y_b=%f,\n" + "\tr=%f, delta=%f, x_half=%f, x_f=%f, x_t=%f (s=%.0f),\n" + "\tx_itp=%f, y_itp=%f\n", + j, + a, b, y_a, y_b, + r, delta, x_half, x_f, x_t, sigma, + x_itp, y_itp + ); + } + #endif + + // Update interval brackets + if (GT(y_itp, 0.)) { + b = x_itp; + y_b = y_itp; + diff_ba = b - a; + } else if (LT(y_itp, 0.)) { + a = x_itp; + y_a = y_itp; + diff_ba = b - a; + } else { + a = x_itp; + b = x_itp; + diff_ba = 0.; + } + + j++; + } while (diff_ba > tol2 && j <= n_max); + + phi = (a + b) / 2.; + + if (j > n_max + 1 || fabs(diff_ba) > tol2 || LT(phi, 0.) || GT(phi, FXW_h0)) { + // Error if not converged or if `phi` physically impossible + LogError( + logfp, + LOGERROR, + "itp_FXW_for_phi(theta = %f): convergence failed at phi = %.9f [cm H2O] " + "after %d iterations with b - a = %.9f - %.9f = %.9f !<= tol2 = %.9f.", + theta, phi, j, b, a, b - a, tol2 + ); + + phi = SW_MISSING; + + } else { + // convert [cm of H2O at 4 C] to [bar] + phi /= 1019.716; + } + + + #ifdef SWDEBUG + if (debug) { + if (!missing(phi)) { + swprintf( + "itp_FXW_for_phi() = phi = %f [bar]: converged in j=%d/%d steps\n", + phi, j, n_max + ); + } else { + swprintf( + "itp_FXW_for_phi(): failed to converged in j=%d/%d steps\n", + j, n_max + ); + + } + } + #endif + + return phi; +} + + /* =================================================== */ /* Global Function Definitions */ /* --------------------------------------------------- */ @@ -1011,171 +1242,582 @@ RealD SW_SnowDepth(RealD SWE, RealD snowdensity) { } } + +/** + @brief Convert soil water content to soil water potential using + specified soil water retention curve (SWRC) + + SOILWAT2 convenience wrapper for `SWRC_SWCtoSWP()`. + + See #swrc2str() for implemented SWRCs. + + @param[in] swcBulk Soil water content in the layer [cm] + @param[in] *lyr Soil information including + SWRC type, SWRC parameters, + coarse fragments (e.g., gravel), and soil layer width. + + @return Soil water potential [-bar] +*/ +RealD SW_SWRC_SWCtoSWP(RealD swcBulk, SW_LAYER_INFO *lyr) { + return SWRC_SWCtoSWP( + swcBulk, + lyr->swrc_type, + lyr->swrcp, + lyr->fractionVolBulk_gravel, + lyr->width, + LOGFATAL + ); +} + /** - @brief Calculates the soil water potential from soil water content of the - n-th soil layer. + @brief Convert soil water content to soil water potential using + specified soil water retention curve (SWRC) - The equation and its coefficients are based on a - paper by Cosby,Hornberger,Clapp,Ginn, in WATER RESOURCES RESEARCH - June 1984. Moisture retention data was fit to the power function. + See #swrc2str() for implemented SWRCs. The code assumes the following conditions: - * checked by `SW_SIT_init_run()` - * width > 0 - * fractionGravel, sand, clay, and sand + clay in [0, 1] - * checked by function `water_eqn()` - * thetasMatric > 0 - * bMatric != 0 - - @param fractionGravel Fraction of soil containing gravel. - @param swcBulk Soilwater content of the current layer (cm/layer) - @param n Layer number to index the **lyr pointer - - @return soil water potential -**/ - -RealD SW_SWCbulk2SWPmatric(RealD fractionGravel, RealD swcBulk, LyrIndex n) { -/********************************************************************** - -HISTORY: - DATE: April 2, 1992 - 9/1/92 (SLC) if swc comes in as zero, set swpotentl to - upperbnd. (Previously, we flagged this - as an error, and set swpotentl to zero). - - 27-Aug-03 (cwb) removed the upperbnd business. Except for - missing values, swc < 0 is impossible, so it's an error, - and the previous limit of swp to 80 seems unreasonable. - return 0.0 if input value is MISSING - - These are the values for each layer obtained via lyr[n]: - width - width of current soil layer - psisMatric - "saturation" matric potential - thetasMatric - saturated moisture content. - bMatric - see equation below. - swc_lim - limit for matric potential - - LOCAL VARIABLES: - theta1 - volumetric soil water content - - DEFINED CONSTANTS: - barconv - conversion factor from bars to cm water. (i.e. - 1 bar = 1024cm water) - - COMMENT: - See the routine "watreqn" for a description of how the variables - psisMatric, bMatric, binverseMatric, thetasMatric are initialized - **********************************************************************/ - - SW_LAYER_INFO *lyr = SW_Site.lyr[n]; - RealD theta1, theta2, swp = .0; - - if (missing(swcBulk) || ZRO(swcBulk)) - return 0.0; - - if (GT(swcBulk, 0.0)) { - // we have soil moisture - - // calculate matric VWC [cm / cm %] from bulk VWC - theta1 = (swcBulk / lyr->width) * 100. / (1. - fractionGravel); - - // calculate (VWC / VWC(saturated)) ^ b - theta2 = powe(theta1 / lyr->thetasMatric, lyr->bMatric); - - if (isnan(theta2) || ZRO(theta2)) { - LogError(logfp, LOGFATAL, "SW_SWCbulk2SWPmatric(): Year = %d, DOY=%d, Layer = %d:\n" - "\tinvalid value of (theta / theta(saturated)) ^ b = %f (must be != 0)\n", - SW_Model.year, SW_Model.doy, n, theta2); + - checked by `SW_SIT_init_run()` + - width > 0 + - fractionGravel, sand, clay, and sand + clay in [0, 1] + - SWRC parameters checked by `SWRC_check_parameters()`. + + @param[in] swcBulk Soil water content in the layer [cm] + @param[in] swrc_type Identification number of selected SWRC + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @param[in] errmode An error code passed to `LogError()`. + SOILWAT2 uses `LOGFATAL` and fails but + other applications may want to warn only (`LOGWARN`) and return. + + @return Soil water potential [-bar] +*/ +double SWRC_SWCtoSWP( + double swcBulk, + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + const int errmode +) { + double res = SW_MISSING; + + if (LT(swcBulk, 0.) || EQ(gravel, 1.) || LE(width, 0.)) { + LogError( + logfp, + errmode, + "SWRC_SWCtoSWP(): invalid SWC = %.4f (must be >= 0)\n", + swcBulk + ); + + return res; + } + + switch (swrc_type) { + case sw_Campbell1974: + res = SWRC_SWCtoSWP_Campbell1974( + swcBulk, swrcp, gravel, width, + errmode + ); + break; + + case sw_vanGenuchten1980: + res = SWRC_SWCtoSWP_vanGenuchten1980( + swcBulk, swrcp, gravel, width, + errmode + ); + break; + + case sw_FXW: + res = SWRC_SWCtoSWP_FXW( + swcBulk, swrcp, gravel, width, + errmode + ); + break; + + default: + LogError( + logfp, + errmode, + "SWRC (type %d) is not implemented.", + swrc_type + ); + break; + } + + return res; +} + + +/** + @brief Convert soil water content to soil water potential using + Campbell's 1974 \cite Campbell1974 Soil Water Retention Curve + + Parameters are explained in `SWRC_check_parameters_for_Campbell1974()`. + + @note + This function was previously named `SW_SWCbulk2SWPmatric()`. + + @note + `SWRC_SWPtoSWC_Campbell1974()` and `SWRC_SWCtoSWP_Campbell1974()` + are the inverse of each other + for `(phi, theta)` between `(swrcp[0], `theta` at `swrcp[0]`)` + and `(infinity, 0)`. + + @note + The function has a discontinuity at saturated water content for which + the matric potential at the "air-entry suction" point (see `swrcp[0]`) + is returned whereas 0 bar is returned for larger values. + + @param[in] swcBulk Soil water content in the layer [cm] + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @param[in] errmode An error code passed to `LogError()`. + SOILWAT2 uses `LOGFATAL` and fails but + other applications may want to warn only (`LOGWARN`) and return. + + @return Soil water potential [-bar] +*/ +double SWRC_SWCtoSWP_Campbell1974( + double swcBulk, + double *swrcp, + double gravel, + double width, + const int errmode +) { + // assume that we have soil moisture + double theta, tmp, res; + + // convert bulk SWC [cm] to theta = matric VWC [cm / cm] + theta = swcBulk / (width * (1. - gravel)); + + if (GT(theta, swrcp[1])) { + // `theta` should not become larger than `theta_sat`; + // however, "Cosby1984AndOthers" does not use `swrcp[1]` for `theta_sat` + // which can lead to inconsistencies; thus, + // we return with 0 instead of, correctly, with errmode and SW_MISSING + res = 0.; + + } else { + // calculate (theta / theta_s) ^ b + tmp = powe(theta / swrcp[1], swrcp[2]); + + if (LE(tmp, 0.)) { + LogError( + logfp, + errmode, + "SWRC_SWCtoSWP_Campbell1974(): invalid value of\n" + "\t(theta / theta(saturated)) ^ b = (%f / %f) ^ %f =\n" + "\t= %f (must be > 0)\n", + theta, swrcp[1], swrcp[2], tmp + ); + + return SW_MISSING; + } + + res = swrcp[0] / tmp; + } + + // convert [cm of H20; SOILWAT2 legacy value] to [bar] + return res / 1024.; +} + + +/** + @brief Convert soil water content to soil water potential using + van Genuchten 1980 \cite vanGenuchten1980 Soil Water Retention Curve + + Parameters are explained in `SWRC_check_parameters_for_vanGenuchten1980()`. + + @note + `SWRC_SWPtoSWC_vanGenuchten1980()` and `SWRC_SWCtoSWP_vanGenuchten1980()` + are the inverse of each other + for `(phi, theta)` between `(0, theta_sat)` and `(infinity, theta_min)`. + + @param[in] swcBulk Soil water content in the layer [cm] + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @param[in] errmode An error code passed to `LogError()`. + SOILWAT2 uses `LOGFATAL` and fails but + other applications may want to warn only (`LOGWARN`) and return. + + @return Soil water potential [-bar] +*/ +double SWRC_SWCtoSWP_vanGenuchten1980( + double swcBulk, + double *swrcp, + double gravel, + double width, + const int errmode +) { + double res, tmp, theta; + + // convert bulk SWC [cm] to theta = matric VWC [cm / cm] + theta = swcBulk / (width * (1. - gravel)); + + // calculate if theta in ]theta_min, theta_sat] + if (GT(theta, swrcp[0])) { + if (LT(theta, swrcp[1])) { + // calculate inverse of normalized theta + tmp = (swrcp[1] - swrcp[0]) / (theta - swrcp[0]); + + // calculate tension [cm of H20] + tmp = powe(tmp, 1. / (1. - 1. / swrcp[3])); // tmp values are >= 1 + res = pow(-1. + tmp, 1. / swrcp[3]) / swrcp[2]; // `pow()` because x >= 0 + + // convert [cm of H2O at 4 C; value from `soilDB::KSSL_VG_model()`] to [bar] + res /= 1019.716; + + } else if (EQ(theta, swrcp[1])) { + // theta is theta_sat + res = 0; + } else { - swp = lyr->psisMatric / theta2 / BARCONV; + // theta is > theta_sat + LogError( + logfp, + errmode, + "SWRC_SWCtoSWP_vanGenuchten1980(): invalid value of\n" + "\ttheta = %f (must be <= theta_sat = %f)\n", + theta, swrcp[1] + ); + + res = SW_MISSING; } } else { - LogError(logfp, LOGFATAL, "Invalid SWC value (%.4f) in SW_SWC_swc2potential.\n" - " Year = %d, DOY=%d, Layer = %d\n", swcBulk, SW_Model.year, SW_Model.doy, n); + // theta is <= theta_min + LogError( + logfp, + errmode, + "SWRC_SWCtoSWP_vanGenuchten1980(): invalid value of\n" + "\ttheta = %f (must be > theta_min = %f)\n", + theta, swrcp[0] + ); + + res = SW_MISSING; } - return swp; + return res; } + + /** -@brief Convert soil water potential to bulk volumetric water content. + @brief Convert soil water content to soil water potential using + FXW Soil Water Retention Curve + \cite fredlund1994CGJa, \cite wang2018wrr -@param fractionGravel Fraction of soil containing gravel, percentage. -@param swpMatric lyr->psisMatric calculated in water equation function -@param n Layer of soil. + Parameters are explained in `SWRC_check_parameters_for_FXW()`. -@return Volumentric water content (cm H2O/cm SOIL). -**/ + @note + `SWRC_SWPtoSWC_FXW()` and `SWRC_SWCtoSWP_FXW()` + are the inverse of each other + for `(phi, theta)` between `(0, theta_sat)` and `(6.3 x 10^6 cm, 0)`. + @param[in] swcBulk Soil water content in the layer [cm] + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @param[in] errmode An error code passed to `LogError()`. + SOILWAT2 uses `LOGFATAL` and fails but + other applications may want to warn only (`LOGWARN`) and return. + + @return Soil water potential [-bar] +*/ +double SWRC_SWCtoSWP_FXW( + double swcBulk, + double *swrcp, + double gravel, + double width, + const int errmode +) { + double res, theta; + + // convert bulk SWC [cm] to theta = matric VWC [cm / cm] + theta = swcBulk / (width * (1. - gravel)); + + // calculate if theta in [0, theta_sat] + if (GE(theta, 0.)) { + if (LT(theta, swrcp[0])) { + // calculate tension = phi [bar] + res = itp_FXW_for_phi(theta, swrcp); + + } else if (EQ(theta, swrcp[0])) { + // theta is theta_sat + res = 0.; + + } else { + // theta is > theta_sat + LogError( + logfp, + errmode, + "SWRC_SWCtoSWP_FXW(): invalid value of\n" + "\ttheta = %f (must be <= theta_sat = %f)\n", + theta, swrcp[0] + ); + + res = SW_MISSING; + } + + } else { + // theta is < 0 + LogError( + logfp, + errmode, + "SWRC_SWCtoSWP_FXW(): invalid value of\n" + "\ttheta = %f (must be >= 0)\n", + theta + ); + + res = SW_MISSING; + } + + return res; +} -RealD SW_SWPmatric2VWCBulk(RealD fractionGravel, RealD swpMatric, LyrIndex n) { /** - History: - 27-Aug-03 (cwb) moved from the Site module. -**/ - - SW_LAYER_INFO *lyr = SW_Site.lyr[n]; - RealD t, p; - swpMatric *= BARCONV; - p = powe(lyr->psisMatric / swpMatric, lyr->binverseMatric); // lyr->psisMatric calculated in water equation function | todo: check to make sure these are calculated before - t = lyr->thetasMatric * p * 0.01 * (1 - fractionGravel); - return (t); + @brief Convert soil water potential to soil water content using + specified soil water retention curve (SWRC) + + SOILWAT2 convenience wrapper for `SWRC_SWPtoSWC()`. + + See #swrc2str() for implemented SWRCs. + + @param[in] swpMatric Soil water potential [-bar] + @param[in] *lyr Soil information including + SWRC type, SWRC parameters, + coarse fragments (e.g., gravel), and soil layer width. + + @return Soil water content in the layer [cm] +*/ +RealD SW_SWRC_SWPtoSWC(RealD swpMatric, SW_LAYER_INFO *lyr) { + return SWRC_SWPtoSWC( + swpMatric, + lyr->swrc_type, + lyr->swrcp, + lyr->fractionVolBulk_gravel, + lyr->width, + LOGFATAL + ); } + /** -@brief Calculates 'Brooks-Corey' residual volumetric soil water. + @brief Convert soil water potential to soil water content using + specified soil water retention curve (SWRC) -Equations based on: Rawls WJ, Brakensiek DL (1985) Prediction of soil water properties - for hydrological modeling, based on @cite ASCE1985 + See #swrc2str() for implemented SWRCs. -@param fractionGravel Fraction of soil consisting of gravel, percentage. -@param sand Fraction of soil consisting of sand, percentage. -@param clay Fraction of soil consisting of clay, percentage. -@param porosity Fraction of Soil porosity as the saturated VWC, percentage. + The code assumes the following conditions: + - checked by `SW_SIT_init_run()` + - width > 0 + - fractionGravel, sand, clay, and sand + clay in [0, 1] + - SWRC parameters checked by `SWRC_check_parameters()`. + + @param[in] swpMatric Soil water potential [-bar] + @param[in] swrc_type Identification number of selected SWRC + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + @param[in] errmode An error code passed to `LogError()`. + SOILWAT2 uses `LOGFATAL` and fails but + other applications may want to warn only (`LOGWARN`) and return. + + @return Soil water content in the layer [cm] +*/ +double SWRC_SWPtoSWC( + double swpMatric, + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + const int errmode +) { + double res = SW_MISSING; + + if (LT(swpMatric, 0.)) { + LogError( + logfp, + errmode, + "SWRC_SWPtoSWC(): invalid SWP = %.4f (must be >= 0)\n", + swpMatric + ); + + return res; + } -@returns Residual volumetric soil water (cm/cm) -**/ + switch (swrc_type) { + case sw_Campbell1974: + res = SWRC_SWPtoSWC_Campbell1974(swpMatric, swrcp, gravel, width); + break; -RealD SW_VWCBulkRes(RealD fractionGravel, RealD sand, RealD clay, RealD porosity) { -/*--------------------- -History: - 02/03/2012 (drs) calculates 'Brooks-Corey' residual volumetric soil water based on Rawls WJ, Brakensiek DL (1985) Prediction of soil water properties for hydrological modeling. In Watershed management in the Eighties (eds Jones EB, Ward TJ), pp. 293-299. American Society of Civil Engineers, New York. - however, equation is only valid if (0.05 < clay < 0.6) & (0.05 < sand < 0.7) + case sw_vanGenuchten1980: + res = SWRC_SWPtoSWC_vanGenuchten1980(swpMatric, swrcp, gravel, width); + break; ----------------------*/ + case sw_FXW: + res = SWRC_SWPtoSWC_FXW(swpMatric, swrcp, gravel, width); + break; - if (clay < .05 || clay > .6 || sand < .05 || sand > .7) { - LogError( - logfp, - LOGWARN, - "Sand and/or clay values out of valid range, simulation outputs may differ." - ); - return SW_MISSING; + default: + LogError( + logfp, + errmode, + "SWRC (type %d) is not implemented.", + swrc_type + ); + break; + } - } else { - RealD res; - sand *= 100.; - clay *= 100.; - - res = (1. - fractionGravel) * ( - - 0.0182482 \ - + 0.00087269 * sand \ - + 0.00513488 * clay \ - + 0.02939286 * porosity \ - - 0.00015395 * squared(clay) \ - - 0.0010827 * sand * porosity \ - - 0.00018233 * squared(clay) * squared(porosity) \ - + 0.00030703 * squared(clay) * porosity \ - - 0.0023584 * squared(porosity) * clay - ); - - return (fmax(res, 0.)); - } + return res; +} + + +/** + @brief Convert soil water potential to soil water content using + Campbell's 1974 \cite Campbell1974 Soil Water Retention Curve + + Parameters are explained in `SWRC_check_parameters_for_Campbell1974()`. + + @note + This function was previously named `SW_SWPmatric2VWCBulk()`. + + @note + The function returns saturated water content if `swpMatric` is at or below + the "air-entry suction" point (see `swrcp[0]`). + + @note + `SWRC_SWPtoSWC_Campbell1974()` and `SWRC_SWCtoSWP_Campbell1974()` + are the inverse of each other + for `(phi, theta)` between `(swrcp[0], `theta` at `swrcp[0]`)` + and `(infinity, 0)`. + + @param[in] swpMatric Soil water potential [-bar] + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + + @return Soil water content in the layer [cm] +*/ +double SWRC_SWPtoSWC_Campbell1974( + double swpMatric, + double *swrcp, + double gravel, + double width +) { + double phi, res; + + // convert SWP [-bar] to phi [cm of H20; SOILWAT2 legacy value] + phi = swpMatric * 1024.; + + if (LT(phi, swrcp[0])) { + res = swrcp[1]; // theta_sat + + } else { + // calculate matric theta [cm / cm] + // within [0, theta_sat] because `phi` > swrcp[0] + res = swrcp[1] * powe(swrcp[0] / phi, 1. / swrcp[2]); + } + + // convert matric theta [cm / cm] to bulk SWC [cm] + return (1. - gravel) * width * res; } + +/** + @brief Convert soil water potential to soil water content using + van Genuchten 1980 \cite vanGenuchten1980 Soil Water Retention Curve + + Parameters are explained in `SWRC_check_parameters_for_vanGenuchten1980()`. + + @note + `SWRC_SWPtoSWC_vanGenuchten1980()` and `SWRC_SWCtoSWP_vanGenuchten1980()` + are the inverse of each other + for `(phi, theta)` between `(0, theta_sat)` and `(infinity, theta_min)`. + + @param[in] swpMatric Soil water potential [-bar] + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + + @return Soil water content in the layer [cm] +*/ +double SWRC_SWPtoSWC_vanGenuchten1980( + double swpMatric, + double *swrcp, + double gravel, + double width +) { + double phi, tmp, res; + + // convert SWP [-bar] to phi [cm of H2O at 4 C; + // value from `soilDB::KSSL_VG_model()`] + phi = swpMatric * 1019.716; + + tmp = powe(swrcp[2] * phi, swrcp[3]); + tmp = powe(1. + tmp, 1. - 1. / swrcp[3]); + + // calculate matric theta [cm / cm] which is within [theta_min, theta_sat] + res = swrcp[0] + (swrcp[1] - swrcp[0]) / tmp; + + // convert matric theta [cm / cm] to bulk SWC [cm] + return (1. - gravel) * width * res; +} + + +/** + @brief Convert soil water potential to soil water content using + FXW Soil Water Retention Curve + \cite fredlund1994CGJa, \cite wang2018wrr + + Parameters are explained in `SWRC_check_parameters_for_FXW()`. + + @note + `SWRC_SWPtoSWC_FXW()` and `SWRC_SWCtoSWP_FXW()` + are the inverse of each other + for `(phi, theta)` between `(0, theta_sat)` and `(6.3 x 10^6 cm, 0)`. + + @param[in] swpMatric Soil water potential [-bar] + @param[in] *swrcp Vector of SWRC parameters + @param[in] gravel Coarse fragments (> 2 mm; e.g., gravel) + of the whole soil [m3/m3] + @param[in] width Soil layer width [cm] + + @return Soil water content in the layer [cm] +*/ +double SWRC_SWPtoSWC_FXW( + double swpMatric, + double *swrcp, + double gravel, + double width +) { + double phi, res; + + // convert SWP [-bar] to phi [cm of H2O at 4 C; + // value from `soilDB::KSSL_VG_model()`] + phi = swpMatric * 1019.716; + + res = FXW_phi_to_theta(phi, swrcp); + + // convert matric theta [cm / cm] to bulk SWC [cm] + return (1. - gravel) * width * res; +} + + + + /** @brief This routine sets the known memory refs in this module so they can be checked for leaks, etc. diff --git a/SW_SoilWater.h b/SW_SoilWater.h index ce7a676b6..dc40a94e3 100644 --- a/SW_SoilWater.h +++ b/SW_SoilWater.h @@ -161,11 +161,67 @@ void SW_SWC_adjust_snow(RealD temp_min, RealD temp_max, RealD ppt, RealD *rain, RealD SW_SWC_snowloss(RealD pet, RealD *snowpack); RealD SW_SnowDepth(RealD SWE, RealD snowdensity); void SW_SWC_end_day(void); -RealD SW_SWCbulk2SWPmatric(RealD fractionGravel, RealD swcBulk, LyrIndex n); -RealD SW_SWPmatric2VWCBulk(RealD fractionGravel, RealD swpMatric, LyrIndex n); -RealD SW_VWCBulkRes(RealD fractionGravel, RealD sand, RealD clay, RealD porosity); void get_dSWAbulk(int i); +double SW_SWRC_SWCtoSWP(double swcBulk, SW_LAYER_INFO *lyr); +double SWRC_SWCtoSWP( + double swcBulk, + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + const int errmode +); +double SWRC_SWCtoSWP_Campbell1974( + double swcBulk, + double *swrcp, + double gravel, + double width, + const int errmode +); +double SWRC_SWCtoSWP_vanGenuchten1980( + double swcBulk, + double *swrcp, + double gravel, + double width, + const int errmode +); +double SWRC_SWCtoSWP_FXW( + double swcBulk, + double *swrcp, + double gravel, + double width, + const int errmode +); + +RealD SW_SWRC_SWPtoSWC(RealD swpMatric, SW_LAYER_INFO *lyr); +double SWRC_SWPtoSWC( + double swpMatric, + unsigned int swrc_type, + double *swrcp, + double gravel, + double width, + const int errmode +); +double SWRC_SWPtoSWC_Campbell1974( + double swpMatric, + double *swrcp, + double gravel, + double width +); +double SWRC_SWPtoSWC_vanGenuchten1980( + double swpMatric, + double *swrcp, + double gravel, + double width +); +double SWRC_SWPtoSWC_FXW( + double swpMatric, + double *swrcp, + double gravel, + double width +); + #ifdef SWDEBUG void SW_WaterBalance_Checks(void); #endif diff --git a/SW_VegEstab.c b/SW_VegEstab.c index 4d85f53c2..ed04dc31b 100644 --- a/SW_VegEstab.c +++ b/SW_VegEstab.c @@ -425,14 +425,15 @@ void _spp_init(unsigned int sppnum) { /* The thetas and psis etc should be initialized by now */ /* because init_layers() must be called prior to this routine */ /* (see watereqn() ) */ - v->min_swc_germ = SW_SWPmatric2VWCBulk(lyr[0]->fractionVolBulk_gravel, v->bars[SW_GERM_BARS], 0) * lyr[0]->width; + v->min_swc_germ = SW_SWRC_SWPtoSWC(v->bars[SW_GERM_BARS], lyr[0]); /* due to possible differences in layer textures and widths, we need * to average the estab swc across the given layers to peoperly * compare the actual swc average in the checkit() routine */ v->min_swc_estab = 0.; - for (i = 0; i < v->estab_lyrs; i++) - v->min_swc_estab += SW_SWPmatric2VWCBulk(lyr[i]->fractionVolBulk_gravel, v->bars[SW_ESTAB_BARS], i) * lyr[i]->width; + for (i = 0; i < v->estab_lyrs; i++) { + v->min_swc_estab += SW_SWRC_SWPtoSWC(v->bars[SW_ESTAB_BARS], lyr[i]); + } v->min_swc_estab /= v->estab_lyrs; _sanity_check(sppnum); diff --git a/doc/SOILWAT2.bib b/doc/SOILWAT2.bib index a918d1b06..508883098 100644 --- a/doc/SOILWAT2.bib +++ b/doc/SOILWAT2.bib @@ -1,20 +1,10 @@ -@Book{ASCE1985, - author = {Jones, E.B., Ward, T.J.}, - title = {Watershed management in the eighties}, - publisher = {American Society of civil Engineers}, - year = 1985, - address = {New York}, - isbn = {0872624498}, - pages = {293-299}, -} - @Article{ASCE2000, author = "Walter, I.A., Allen, A.G., et. al", year = 2000, title = "ASCEs Standardized Reference Evapotranspiration Equation", journal = "Watershed Management and Operations Management 2000", pages = {59}, - } +} @report{ASCE2005, @@ -70,6 +60,7 @@ @Article{Cosby1984 journal = "Water Resources Research", volume = 20, pages = {682-690}, + doi = {10.1029/WR020i006p00682} } @Article{Eitzinger2000, @@ -198,7 +189,6 @@ @Article{Parton1984 @phdthesis{Gerrits2010, title = {The role of interception in the hydrological cycle}, ISBN = {9789065622488}, - url = {http://resolver.tudelft.nl/uuid:7dd2523b-2169-4e7e-992c-365d2294d02e}, school = {Technische Universiteit Delft}, author = {Gerrits, A. M. J.}, year = {2010} @@ -333,6 +323,123 @@ @article{huang2018JAMC doi = {10.1175/JAMC-D-17-0334.1} } +@article{Campbell1974, + title = {A Simple Method for Determining Unsaturated Conductivity from Moisture Retention Data}, + author = {Campbell, Gaylon S.}, + date = {1974}, + journaltitle = {Soil Science}, + volume = {117}, + number = {6}, + pages = {311--314}, + doi = {10.1097/00010694-197406000-00001} +} + +@article{vanGenuchten1980, + title = {A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils}, + author = {van Genuchten, M. Th.}, + date = {1980}, + journaltitle = {Soil Science Society of America Journal}, + volume = {44}, + number = {5}, + pages = {892--898}, + doi = {10.2136/sssaj1980.03615995004400050002x} +} + +@article{zhang2017JoH, + title = {Weighted Recalibration of the Rosetta Pedotransfer Model with Improved Estimates of Hydraulic Parameter Distributions and Summary Statistics (Rosetta3)}, + author = {Zhang, Yonggen and Schaap, Marcel G.}, + year = {2017}, + journal = {Journal of Hydrology}, + volume = {547}, + pages = {39--53}, + doi = {10.1016/j.jhydrol.2017.01.004} +} + +@book{brooks1964a, + title = {Hydraulic Properties of Porous Media. Hydrology Papers: 3}, + author = {Brooks, R. H. and Corey, A. T.}, + year = {1964}, + publisher = {Colorado State University}, + address = {Fort Collins, Colorado, USA} +} + +@incollection{rawls1985WmitE, + title = {Prediction of Soil Water Properties for Hydrological Modeling}, + booktitle = {Watershed Management in the Eighties}, + author = {Rawls, W J and Brakensiek, D L}, + editor = {Jones, E B and Ward, T. J.}, + year = {1985}, + pages = {293--299}, + publisher = {American Society of Civil Engineers}, + address = {New York, USA} +} + + +@article{fredlund1994CGJa, + title = {Equations for the soil-water characteristic curve}, + author = {Fredlund, D. G. and Xing, A. Q.}, + year = {1994}, + journal = {Canadian Geotechnical Journal}, + volume = {31}, + pages = {512–-532}, + doi = {10.1139/t94-061} +} + +@article{rudiyanto2020JoH, + title = {Simple Functions for Describing Soil Water Retention and the Unsaturated Hydraulic Conductivity from Saturation to Complete Dryness}, + author = {{Rudiyanto} and Minasny, Budiman and Shah, Ramisah M. and Setiawan, Budi I. and {van Genuchten}, Martinus Th.}, + year = {2020}, + journal = {Journal of Hydrology}, + volume = {588}, + pages = {125041}, + doi = {10.1016/j.jhydrol.2020.125041} +} + +@article{rudiyanto2021G, + title = {Pedotransfer Functions for Estimating Soil Hydraulic Properties from Saturation to Dryness}, + author = {{Rudiyanto} and Minasny, Budiman and Chaney, Nathaniel W. and Maggi, Federico and Goh Eng Giap, Sunny and Shah, Ramisah M. and Fiantis, Dian and Setiawan, Budi I.}, + year = {2021}, + journal = {Geoderma}, + volume = {403}, + pages = {115194}, + issn = {0016-7061}, + doi = {10.1016/j.geoderma.2021.115194} +} + + +@article{wang2018wrr, + title = {Alternative Model for Predicting Soil Hydraulic Conductivity Over the Complete Moisture Range}, + author = {Wang, Yunquan and Jin, Menggui and Deng, Zijuan}, + year = {2018}, + journal = {Water Resources Research}, + volume = {54}, + number = {9}, + pages = {6860--6876}, + doi = {10.1029/2018WR023037} +} + +@article{wang2022WRRa, + title = {Development of a New Pedotransfer Function Addressing Limitations in Soil Hydraulic Models and Observations}, + author = {Wang, Yunquan and Zhou, Jieliang and Ma, Rui and Zhu, Gaofeng and Zhang, Yongyong}, + year = {2022}, + journal = {Water Resources Research}, + volume = {58}, + number = {6}, + pages = {e2021WR031406}, + doi = {10.1029/2021WR031406} +} + +@article{oliveira2021ATMS, + title = {An Enhancement of the Bisection Method Average Performance Preserving Minmax Optimality}, + author = {Oliveira, I. F. D. and Takahashi, R. H. C.}, + year = {2021}, + journal = {ACM Transactions on Mathematical Software}, + volume = {47}, + number = {1}, + pages = {1--24}, + doi = {10.1145/3423597} +} + @article{marsaglia1964SR, title = {A Convenient Method for Generating Normal Variables}, author = {Marsaglia, G. and Bray, T. A.}, diff --git a/doc/additional_pages/SOILWAT2_Inputs.md b/doc/additional_pages/SOILWAT2_Inputs.md index ab7d23509..926d1b3c2 100644 --- a/doc/additional_pages/SOILWAT2_Inputs.md +++ b/doc/additional_pages/SOILWAT2_Inputs.md @@ -41,6 +41,7 @@ SOILWAT2 needs the following input files for a simulation run: \refitem yearsin years.in \refitem siteparamin siteparam.in \refitem soilsin soils.in +\refitem swrcpin swrc_params.in \refitem weathsetupin weathsetup.in \refitem mkvprobin mkv_prob.in \refitem mkvcovarin mkv_covar.in @@ -80,6 +81,12 @@ Go back to the \ref explain_inputs "list of input files" Go back to the \ref explain_inputs "list of input files"
+\section swrcpin swrc_params.in +\verbinclude testing/Input/swrc_params.in + +Go back to the \ref explain_inputs "list of input files" +
+ \section weathsetupin weathsetup.in \verbinclude testing/Input/weathsetup.in diff --git a/test/test_SW_Flow_Lib.cc b/test/test_SW_Flow_Lib.cc index f47bfd4d9..284576820 100644 --- a/test/test_SW_Flow_Lib.cc +++ b/test/test_SW_Flow_Lib.cc @@ -451,7 +451,6 @@ namespace shift = 45, shape = 0.1, inflec = 0.25, range = 0.5; double swc[25], width[25], lyrEvapCo[25]; - double swc0[25] = { 0. }; // for test if swc = 0 // Loop over tests with varying number of soil layers for (k = 0; k < 2; k++) @@ -510,15 +509,6 @@ namespace "pot_soil_evap != 0 if fbse = 0 for " << nelyrs << " soil layers"; - // Begin TEST if (swc = 0) - pot_soil_evap(&bserate, nelyrs, lyrEvapCo, totagb, fbse, petday, shift, - shape, inflec, range, width, swc0, Es_param_limit); - - // expect baresoil evaporation rate = 0 if swc = 0 - EXPECT_DOUBLE_EQ(bserate, 0.) << - "pot_soil_evap != 0 if swc = 0 for " << nelyrs << " soil layers"; - - // Begin TEST if (totagb < Es_param_limit) pot_soil_evap(&bserate, nelyrs, lyrEvapCo, totagb, fbse, petday, shift, shape, inflec, range, width, swc, Es_param_limit); diff --git a/test/test_SW_Site.cc b/test/test_SW_Site.cc index a4d34bc62..cfbbca391 100644 --- a/test/test_SW_Site.cc +++ b/test/test_SW_Site.cc @@ -40,67 +40,404 @@ namespace { - // Test the water equation function 'water_eqn' - TEST(SWSiteTest, WaterEquation) { - - //declare inputs - RealD fractionGravel = 0.1, sand = .33, clay =.33; - LyrIndex n = 1; - - water_eqn(fractionGravel, sand, clay, n); - - // Test swcBulk_saturated - EXPECT_GT(SW_Site.lyr[n]->swcBulk_saturated, 0.); // The swcBulk_saturated should be greater than 0 - EXPECT_LT(SW_Site.lyr[n]->swcBulk_saturated, SW_Site.lyr[n]->width); // The swcBulk_saturated can't be greater than the width of the layer - - // Test thetasMatric - EXPECT_GT(SW_Site.lyr[n]->thetasMatric, 36.3); /* Value should always be greater - than 36.3 based upon complete consideration of potential range of sand and clay values */ - EXPECT_LT(SW_Site.lyr[n]->thetasMatric, 46.8); /* Value should always be less - than 46.8 based upon complete consideration of potential range of sand and clay values */ - EXPECT_DOUBLE_EQ(SW_Site.lyr[n]->thetasMatric, 44.593); /* If sand is .33 and - clay is .33, thetasMatric should be 44.593 */ - - // Test psisMatric - EXPECT_GT(SW_Site.lyr[n]->psisMatric, 3.890451); /* Value should always be greater - than 3.890451 based upon complete consideration of potential range of sand and clay values */ - EXPECT_LT(SW_Site.lyr[n]->psisMatric, 34.67369); /* Value should always be less - than 34.67369 based upon complete consideration of potential range of sand and clay values */ - EXPECT_DOUBLE_EQ(SW_Site.lyr[n]->psisMatric, 27.586715750763947); /* If sand is - .33 and clay is .33, psisMatric should be 27.5867 */ - - // Test bMatric - EXPECT_GT(SW_Site.lyr[n]->bMatric, 2.8); /* Value should always be greater than - 2.8 based upon complete consideration of potential range of sand and clay values */ - EXPECT_LT(SW_Site.lyr[n]->bMatric, 18.8); /* Value should always be less - than 18.8 based upon complete consideration of potential range of sand and clay values */ - EXPECT_DOUBLE_EQ(SW_Site.lyr[n]->bMatric, 8.182); /* If sand is .33 and clay is .33, - thetasMatric should be 8.182 */ + // List SWRC: PTFs + const char *ns_ptfca2C1974[] = { + "Campbell1974", + "Cosby1984AndOthers", "Cosby1984" + }; + const char *ns_ptfa2vG1980[] = { + "vanGenuchten1980", + // all PTFs + "Rosetta3" + }; + const char *ns_ptfc2vG1980[] = { + "vanGenuchten1980" + // PTFs implemented in SOILWAT2 + }; + const char *ns_ptfa2FXW[] = { + "FXW" + // all PTFs + "neuroFX2021" + }; + const char *ns_ptfc2FXW[] = { + "FXW" + // PTFs implemented in SOILWAT2 + }; + + // Test pedotransfer functions + TEST(SiteTest, PTFs) { + // inputs + RealD + swrcp[SWRC_PARAM_NMAX], + sand = 0.33, + clay = 0.33, + gravel = 0.1, + bdensity = 1.4; + unsigned int swrc_type, k; + + + //--- Matching PTF-SWRC pairs + // (k starts at 1 because 0 holds the SWRC) + + swrc_type = encode_str2swrc((char *) ns_ptfca2C1974[0]); + for (k = 1; k < length(ns_ptfca2C1974); k++) { + SWRC_PTF_estimate_parameters( + encode_str2ptf((char *) ns_ptfca2C1974[k]), + swrcp, + sand, + clay, + gravel, + bdensity + ); + EXPECT_TRUE((bool) SWRC_check_parameters(swrc_type, swrcp)); + } + + swrc_type = encode_str2swrc((char *) ns_ptfc2vG1980[0]); + for (k = 1; k < length(ns_ptfc2vG1980); k++) { + SWRC_PTF_estimate_parameters( + encode_str2ptf((char *) ns_ptfc2vG1980[k]), + swrcp, + sand, + clay, + gravel, + bdensity + ); + EXPECT_TRUE((bool) SWRC_check_parameters(swrc_type, swrcp)); + } + + swrc_type = encode_str2swrc((char *) ns_ptfc2FXW[0]); + for (k = 1; k < length(ns_ptfc2FXW); k++) { + SWRC_PTF_estimate_parameters( + encode_str2ptf((char *) ns_ptfc2FXW[k]), + swrcp, + sand, + clay, + gravel, + bdensity + ); + EXPECT_TRUE((bool) SWRC_check_parameters(swrc_type, swrcp)); + } + } + + + // Test fatal failures of PTF estimation + TEST(SiteDeathTest, PTFs) { + + RealD + swrcp[SWRC_PARAM_NMAX], + sand = 0.33, + clay = 0.33, + gravel = 0.1, + bdensity = 1.4; + unsigned int ptf_type; + + + //--- Test unimplemented PTF + ptf_type = N_PTFs + 1; + + EXPECT_DEATH_IF_SUPPORTED( + SWRC_PTF_estimate_parameters( + ptf_type, + swrcp, + sand, + clay, + gravel, + bdensity + ), + "@ generic.c LogError" + ); + } + + + // Test PTF-SWRC pairings + TEST(SiteTest, PTF2SWRC) { + unsigned int k; // `length()` returns "unsigned long" + + for (k = 1; k < length(ns_ptfca2C1974); k++) { + EXPECT_TRUE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfca2C1974[0], + (char *) ns_ptfca2C1974[k] + ) + ); + + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfa2vG1980[0], + (char *) ns_ptfca2C1974[k] + ) + ); + + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfa2FXW[0], + (char *) ns_ptfca2C1974[k] + ) + ); + } + + for (k = 1; k < length(ns_ptfa2vG1980); k++) { + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfa2vG1980[0], + (char *) ns_ptfa2vG1980[k] + ) + ); + + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfca2C1974[0], + (char *) ns_ptfa2vG1980[k] + ) + ); + + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfa2FXW[0], + (char *) ns_ptfa2vG1980[k] + ) + ); + } + + + for (k = 1; k < length(ns_ptfa2FXW); k++) { + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfa2FXW[0], + (char *) ns_ptfa2FXW[k] + ) + ); + + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfca2C1974[0], + (char *) ns_ptfa2FXW[k] + ) + ); + + EXPECT_FALSE( + (bool) check_SWRC_vs_PTF( + (char *) ns_ptfa2vG1980[0], + (char *) ns_ptfa2FXW[k] + ) + ); + } - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); } - // Test that water equation function 'water_eqn' fails - TEST(SWSiteDeathTest, WaterEquation) { - //declare inputs - RealD fractionGravel = 0.1; - LyrIndex n = 1; + // Test fatal failures of SWRC parameter checks + TEST(SiteDeathTest, SWRCpChecks) { + + // inputs + RealD swrcp[SWRC_PARAM_NMAX]; + unsigned int swrc_type; + - // Test that error will be logged when b_matric is 0 - RealD sand = 10. + 1./3.; // So that bmatric will equal 0, even though this is a very irrealistic value - RealD clay = 0; + //--- Test unimplemented SWRC + swrc_type = N_SWRCs + 1; EXPECT_DEATH_IF_SUPPORTED( - water_eqn(fractionGravel, sand, clay, n), + SWRC_check_parameters(swrc_type, swrcp), "@ generic.c LogError" ); } + // Test nonfatal failures of SWRC parameter checks + TEST(SiteTest, SWRCpChecks) { + + // inputs + RealD + swrcp[SWRC_PARAM_NMAX], + tmp; + unsigned int swrc_type; + + + //--- SWRC: Campbell1974 + swrc_type = encode_str2swrc((char *) "Campbell1974"); + memset(swrcp, 0., SWRC_PARAM_NMAX * sizeof(swrcp[0])); + swrcp[0] = 24.2159; + swrcp[1] = 0.4436; + swrcp[2] = 10.3860; + swrcp[3] = 14.14351; + EXPECT_TRUE((bool) SWRC_check_parameters(swrc_type, swrcp)); + + // Param1 = psi_sat (> 0) + tmp = swrcp[0]; + swrcp[0] = -1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[0] = tmp; + + // Param2 = theta_sat (0-1) + tmp = swrcp[1]; + swrcp[1] = -1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[1] = 1.5; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[1] = tmp; + + // Param3 = beta (!= 0) + tmp = swrcp[2]; + swrcp[2] = 0.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[2] = tmp; + + + + //--- Fail SWRC: vanGenuchten1980 + swrc_type = encode_str2swrc((char *) "vanGenuchten1980"); + memset(swrcp, 0., SWRC_PARAM_NMAX * sizeof(swrcp[0])); + swrcp[0] = 0.1246; + swrcp[1] = 0.4445; + swrcp[2] = 0.0112; + swrcp[3] = 1.2673; + swrcp[4] = 7.7851; + EXPECT_TRUE((bool) SWRC_check_parameters(swrc_type, swrcp)); + + + // Param1 = theta_res (0-1) + tmp = swrcp[0]; + swrcp[0] = -1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[0] = 1.5; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[0] = tmp; + + // Param2 = theta_sat (0-1 & > theta_res) + tmp = swrcp[1]; + swrcp[1] = -1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[1] = 1.5; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[1] = 0.5 * swrcp[0]; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[1] = tmp; + + // Param3 = alpha (> 0) + tmp = swrcp[2]; + swrcp[2] = 0.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[2] = tmp; + + // Param4 = n (> 1) + tmp = swrcp[3]; + swrcp[3] = 1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[3] = tmp; + + + + + + //--- Fail SWRC: FXW + swrc_type = encode_str2swrc((char *) "FXW"); + memset(swrcp, 0., SWRC_PARAM_NMAX * sizeof(swrcp[0])); + swrcp[0] = 0.437461; + swrcp[1] = 0.050757; + swrcp[2] = 1.247689; + swrcp[3] = 0.308681; + swrcp[4] = 22.985379; + swrcp[5] = 2.697338; + EXPECT_TRUE((bool) SWRC_check_parameters(swrc_type, swrcp)); + + + // Param1 = theta_sat (0-1) + tmp = swrcp[0]; + swrcp[0] = -1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[0] = 1.5; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[0] = tmp; + + // Param2 = alpha (> 0) + tmp = swrcp[1]; + swrcp[1] = 0.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[1] = tmp; + + // Param3 = n (> 1) + tmp = swrcp[2]; + swrcp[2] = 1.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[2] = tmp; + + // Param4 = m (> 0) + tmp = swrcp[3]; + swrcp[3] = 0.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[3] = tmp; + + // Param5 = Ksat (> 0) + tmp = swrcp[4]; + swrcp[4] = 0.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[4] = tmp; + + // Param6 = L (> 0) + tmp = swrcp[5]; + swrcp[5] = 0.; + EXPECT_FALSE((bool) SWRC_check_parameters(swrc_type, swrcp)); + swrcp[5] = tmp; + } + + + // Test 'PTF_RawlsBrakensiek1985' + TEST(SiteTest, PTFRawlsBrakensiek1985) { + //declare mock INPUTS + double + theta_min, + clay = 0.1, + sand = 0.6, + porosity = 0.4; + int k1, k2, k3; + + //--- 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); + EXPECT_DOUBLE_EQ(theta_min, SW_MISSING); + + PTF_RawlsBrakensiek1985(&theta_min, 0.75, clay, porosity); + EXPECT_DOUBLE_EQ(theta_min, SW_MISSING); + + PTF_RawlsBrakensiek1985(&theta_min, sand, 0., porosity); + EXPECT_DOUBLE_EQ(theta_min, SW_MISSING); + + PTF_RawlsBrakensiek1985(&theta_min, sand, 0.65, porosity); + EXPECT_DOUBLE_EQ(theta_min, SW_MISSING); + + PTF_RawlsBrakensiek1985(&theta_min, sand, clay, 0.); + EXPECT_DOUBLE_EQ(theta_min, SW_MISSING); + + PTF_RawlsBrakensiek1985(&theta_min, sand, clay, 1.); + EXPECT_DOUBLE_EQ(theta_min, SW_MISSING); + + + // Check that `theta_min` is reasonable over ranges of soil properties + for (k1 = 0; k1 <= 5; k1++) { + sand = 0.05 + (double) k1 / 5. * (0.7 - 0.05); + + for (k2 = 0; k2 <= 5; k2++) { + clay = 0.05 + (double) k2 / 5. * (0.6 - 0.05); + + for (k3 = 0; k3 <= 5; k3++) { + porosity = 0.1 + (double) k3 / 5. * (0.99 - 0.1); + + PTF_RawlsBrakensiek1985(&theta_min, sand, clay, porosity); + EXPECT_GE(theta_min, 0.); + EXPECT_LT(theta_min, porosity); + } + } + } + + // 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); + EXPECT_DOUBLE_EQ(theta_min, 0); + } + + // Test that `SW_SIT_init_run` fails on bad soil inputs - TEST(SWSiteDeathTest, SoilParameters) { + TEST(SiteDeathTest, SoilParameters) { LyrIndex n1 = 0, n2 = 1, k = 2; RealD help; diff --git a/test/test_SW_SoilWater.cc b/test/test_SW_SoilWater.cc index 2a6c50275..17f0629ec 100644 --- a/test/test_SW_SoilWater.cc +++ b/test/test_SW_SoilWater.cc @@ -1,8 +1,11 @@ #include "gtest/gtest.h" #include +#include #include #include -#include +#include +#include +#include #include #include "../generic.h" #include "../filefuncs.h" @@ -20,54 +23,8 @@ namespace{ - // Test the 'SW_SoilWater' function 'SW_VWCBulkRes' - TEST(SWSoilWaterTest, VWCBulkRes){ - //declare mock INPUTS - RealD fractionGravel = .1; - RealD clay = .7; - RealD sand = .2; - RealD porosity = 1; - - RealD res = SW_VWCBulkRes(fractionGravel, sand, clay, porosity); - // when clay > .6, we expect res == SW_MISSING since this isn't within reasonable - // range - EXPECT_DOUBLE_EQ(res, SW_MISSING); - - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); - - clay = .5; - sand = .04; - - res = SW_VWCBulkRes(fractionGravel, sand, clay, porosity); - // when sand < .05, we expect res == SW_MISSING since this isn't within reasonable - // range - EXPECT_DOUBLE_EQ(res, SW_MISSING); - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); - - sand = .4; - porosity = .4; - - res = SW_VWCBulkRes(fractionGravel, sand, clay, porosity); - // when sand == .4, clay == .5, porosity == .4 and fractionGravel ==.1, - // we expect res == .088373829599999967 - EXPECT_DOUBLE_EQ(res, .088373829599999967); - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); - - porosity = .1; - - res = SW_VWCBulkRes(fractionGravel, sand, clay, porosity); - // when sand == .4, clay == .5, porosity == .1 and fractionGravel ==.1, - // we expect res == 0 - EXPECT_DOUBLE_EQ(res, 0); - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); - } - // Test the 'SW_SoilWater' function 'SW_SWC_adjust_snow' - TEST(SWSoilWaterTest, SWCdjustSnow){ + TEST(SWSoilWaterTest, SWCadjustSnow){ // setup mock variables SW_Site.TminAccu2 = 0; SW_Model.doy = 1; @@ -112,153 +69,306 @@ namespace{ EXPECT_EQ(snowmelt, 0); } - // Test the 'SW_SoilWater' function 'SW_SWCbulk2SWPmatric' - TEST(SWSoilWaterTest, SWCbulk2SWPmatric){ - // Note: function `SW_SWCbulk2SWPmatric` accesses `SW_Site.lyr[n]` - - RealD tol = 1e-2; // pedotransfer functions are not very exact - RealD fractionGravel = 0.2; - RealD swcBulk; - RealD res; - RealD help; - LyrIndex n = 1; - - // when swc is 0, we expect res == 0 - res = SW_SWCbulk2SWPmatric(fractionGravel, 0., n); - EXPECT_EQ(res, 0.0); - - // when swc is SW_MISSING, we expect res == 0 - res = SW_SWCbulk2SWPmatric(fractionGravel, SW_MISSING, n); - EXPECT_EQ(res, 0.0); - - // if swc > field capacity, then we expect res < 0.33 bar - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, - SW_Site.lyr[n]->swcBulk_fieldcap + 0.1, n); - EXPECT_LT(res, 0.33 + tol); - - // if swc = field capacity, then we expect res == 0.33 bar - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, - SW_Site.lyr[n]->swcBulk_fieldcap, n); - EXPECT_NEAR(res, 0.33, tol); - - // if field capacity > swc > wilting point, then - // we expect 15 bar > res > 0.33 bar - swcBulk = (SW_Site.lyr[n]->swcBulk_fieldcap + - SW_Site.lyr[n]->swcBulk_wiltpt) / 2; - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, - swcBulk, n); - EXPECT_GT(res, 0.33 - tol); - EXPECT_LT(res, 15 + tol); - - // if swc = wilting point, then we expect res == 15 bar - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, - SW_Site.lyr[n]->swcBulk_wiltpt, n); - EXPECT_NEAR(res, 15., tol); - - // if swc < wilting point, then we expect res > 15 bar - swcBulk = (SW_Site.lyr[n]->swcBulk_wiltpt) / 2; - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, - swcBulk, n); - EXPECT_GT(res, 15. - tol); - - - // ------ meddling with internal value - // if fractionGravel == 1: no soil volume available to hold any soil water - // this would also lead to theta1 == 0 and division by zero - // this situation does normally not occur because it is - // checked during input by function `_read_layers` - // Note: this situation is tested by the death test - // `SWSWCbulk2SWPmatricDeathTest`: we cannot test it here because the - // Address Sanitizer would complain with `UndefinedBehaviorSanitizer` - // see [issue #231](https://github.com/DrylandEcology/SOILWAT2/issues/231) - // res = SW_SWCbulk2SWPmatric(1., SW_Site.lyr[n]->swcBulk_fieldcap, n); - // EXPECT_DOUBLE_EQ(res, 0.); // SWP "ought to be" infinity [bar] - - // if theta(sat, matric; Cosby et al. 1984) == 0: would be division by zero - // this situation does normally not occur because it is - // checked during input by function `water_eqn` - help = SW_Site.lyr[n]->thetasMatric; - SW_Site.lyr[n]->thetasMatric = 0.; - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, 0., n); - EXPECT_DOUBLE_EQ(res, 0.); // SWP "ought to be" infinity [bar] - SW_Site.lyr[n]->thetasMatric = help; - - // if lyr->width == 0: would be division by zero - // this situation does normally not occur because it is - // checked during input by function `_read_layers` - help = SW_Site.lyr[n]->bMatric; - SW_Site.lyr[n]->width = 0.; - res = SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, 0., n); - EXPECT_DOUBLE_EQ(res, 0.); // swc < width - SW_Site.lyr[n]->width = help; - - - // No need to reset to previous global states because we didn't change any - // global states - // Reset_SOILWAT2_after_UnitTest(); + + // Test the 'SW_SoilWater' functions 'SWRC_SWCtoSWP' and `SWRC_SWPtoSWC` + TEST(SWSoilWaterTest, TranslateBetweenSWCandSWP) { + // set up mock variables + unsigned int swrc_type, ptf_type, k; + const int em = LOGFATAL; + RealD + phi, + swcBulk, swc_sat, swc_fc, swc_wp, + swp, + swrcp[SWRC_PARAM_NMAX], + sand = 0.33, + clay = 0.33, + gravel = 0.2, + bdensity = 1.4, + width = 10., + // SWP values in [0, Inf[ but FXW maxes out at 6178.19079 bar + swpsb[12] = { + 0., 0.001, 0.01, 0.026, 0.027, 0.33, 15., 30., 100., 300., 1000., 6178. + }, + // SWP values in [fc, Inf[ but FXW maxes out at 6178.19079 bar + swpsi[7] = {0.33, 15., 30., 100., 300., 1000., 6178.}; + + std::ostringstream msg; + + + + // Loop over SWRCs + for (swrc_type = 0; swrc_type < N_SWRCs; swrc_type++) { + memset(swrcp, 0., SWRC_PARAM_NMAX * sizeof(swrcp[0])); + + // Find a suitable PTF to generate `SWRCp` + for ( + ptf_type = 0; + ptf_type < N_PTFs && !check_SWRC_vs_PTF( + (char *) swrc2str[swrc_type], + (char *) ptf2str[ptf_type] + ); + ptf_type++ + ) {} + + + // Obtain SWRCp + if (ptf_type < N_PTFs) { + // PTF implemented in C: estimate parameters + SWRC_PTF_estimate_parameters( + ptf_type, + swrcp, + sand, + clay, + gravel, + bdensity + ); + + } else { + // PTF not implemented in C: provide hard coded values + if ( + Str_CompareI( + (char *) swrc2str[swrc_type], + (char *) "vanGenuchten1980" + ) == 0 + ) { + swrcp[0] = 0.11214750; + swrcp[1] = 0.4213539; + swrcp[2] = 0.007735474; + swrcp[3] = 1.344678; + swrcp[4] = 7.78506; + + } else if ( + Str_CompareI( + (char *) swrc2str[swrc_type], + (char *) "FXW" + ) == 0 + ) { + swrcp[0] = 0.437461; + swrcp[1] = 0.050757; + swrcp[2] = 1.247689; + swrcp[3] = 0.308681; + swrcp[4] = 22.985379; + swrcp[5] = 2.697338; + + } else { + FAIL() << "No SWRC parameters available for " << swrc2str[swrc_type]; + } + } + + + //------ Tests SWC -> SWP + msg.str(""); + msg << "SWRC/PTF = " << swrc_type << "/" << ptf_type; + + // preferably we would print names instead of type codes, but + // this leads to "global-buffer-overflow" + // 0 bytes to the right of global variable 'ptf2str' + //msg << "SWRC/PTF = " << swrc2str[swrc_type] << "/" << ptf2str[ptf_type]; + + swc_sat = SWRC_SWPtoSWC(0., swrc_type, swrcp, gravel, width, em); + swc_fc = SWRC_SWPtoSWC(1. / 3., swrc_type, swrcp, gravel, width, em); + swc_wp = SWRC_SWPtoSWC(15., swrc_type, swrcp, gravel, width, em); + + + // if swc = saturation, then we expect phi in [0, fc] + // for instance, Campbell1974 goes to (theta_sat, swrcp[0]) instead of 0 + swp = SWRC_SWCtoSWP(swc_sat, swrc_type, swrcp, gravel, width, em); + EXPECT_GE(swp, 0.) << msg.str(); + EXPECT_LT(swp, 1. / 3.) << msg.str(); + + + // if swc > field capacity, then we expect phi < 0.33 bar + swcBulk = (swc_sat + swc_fc) / 2.; + EXPECT_LT( + SWRC_SWCtoSWP(swcBulk, swrc_type, swrcp, gravel, width, em), + 1. / 3. + ) << msg.str(); + + // if swc = field capacity, then we expect phi == 0.33 bar + EXPECT_NEAR( + SWRC_SWCtoSWP(swc_fc, swrc_type, swrcp, gravel, width, em), + 1. / 3., + tol9 + ) << msg.str(); + + // if field capacity > swc > wilting point, then + // we expect 15 bar > phi > 0.33 bar + swcBulk = (swc_wp + swc_fc) / 2.; + phi = SWRC_SWCtoSWP(swcBulk, swrc_type, swrcp, gravel, width, em); + EXPECT_GT(phi, 1. / 3.) << msg.str(); + EXPECT_LT(phi, 15.) << msg.str(); + + // if swc = wilting point, then we expect phi == 15 bar + EXPECT_NEAR( + SWRC_SWCtoSWP(swc_wp, swrc_type, swrcp, gravel, width, em), + 15., + tol9 + ) << msg.str(); + + // if swc < wilting point, then we expect phi > 15 bar + swcBulk = SWRC_SWPtoSWC(2. * 15., swrc_type, swrcp, gravel, width, em); + EXPECT_GT( + SWRC_SWCtoSWP(swcBulk, swrc_type, swrcp, gravel, width, em), + 15. + ) << msg.str(); + + + + //------ Tests SWP -> SWC + // when fractionGravel is 1, we expect theta == 0 + EXPECT_EQ( + SWRC_SWPtoSWC(15., swrc_type, swrcp, 1., width, em), + 0. + ) << msg.str(); + + // when width is 0, we expect theta == 0 + EXPECT_EQ( + SWRC_SWPtoSWC(15., swrc_type, swrcp, gravel, 0., em), + 0. + ) << msg.str(); + + // check bounds of swc + for (k = 0; k < length(swpsb); k++) { + swcBulk = SWRC_SWPtoSWC(swpsb[k], swrc_type, swrcp, gravel, width, em); + EXPECT_GE(swcBulk, 0.) << + msg.str() << " at SWP = " << swpsb[k] << " bar"; + EXPECT_LE(swcBulk, width * (1. - gravel)) << + msg.str() << " at SWP = " << swpsb[k] << " bar"; + } + + + //------ Tests that both SWP <-> SWC are inverse of each other + // for phi at 0 (saturation) and phi in [fc, infinity] + // but not necessarily if phi in ]0, fc[; + // for instance, Campbell1974 is not inverse in ]0, swrcp[0][ + for (k = 0; k < length(swpsi); k++) { + swcBulk = SWRC_SWPtoSWC(swpsi[k], swrc_type, swrcp, gravel, width, em); + swp = SWRC_SWCtoSWP(swcBulk, swrc_type, swrcp, gravel, width, em); + + EXPECT_NEAR(swp, swpsi[k], tol9) << + msg.str() << " at SWP = " << swpsi[k] << " bar"; + EXPECT_NEAR( + SWRC_SWPtoSWC(swp, swrc_type, swrcp, gravel, width, em), + swcBulk, + tol9 + ) << msg.str() << " at SWC = " << swcBulk << " cm"; + } + } } - TEST(SWSoilWaterDeathTest, SWCbulk2SWPmatricDeathTest) { - LyrIndex n = 1; - RealD help; - // we expect fatal errors and write to log under two situations: + // Death Tests of 'SW_SoilWater' function 'SWRC_SWCtoSWP' + TEST(SoilWaterDeathTest, SWCtoSWP) { + // set up mock variables + RealD + swrcp[SWRC_PARAM_NMAX], + gravel = 0.1, + width = 10.; + + unsigned int swrc_type; + - // if swc < 0: water content can physically not be negative + //--- we expect (non-)fatal errors in a few situations + // (fatality depends on the error mode) + + //--- 1) Unimplemented SWRC + swrc_type = N_SWRCs + 1; EXPECT_DEATH_IF_SUPPORTED( - SW_SWCbulk2SWPmatric(SW_Site.lyr[n]->fractionVolBulk_gravel, -1., n), + SWRC_SWCtoSWP(1., swrc_type, swrcp, gravel, width, LOGFATAL), "@ generic.c LogError" ); + EXPECT_DOUBLE_EQ( + SWRC_SWCtoSWP(1., swrc_type, swrcp, gravel, width, LOGWARN), + SW_MISSING + ); + + + // --- 2) swc < 0: water content cannot be negative + for (swrc_type = 0; swrc_type < N_SWRCs; swrc_type++) { + EXPECT_DEATH_IF_SUPPORTED( + SWRC_SWCtoSWP(-1., swrc_type, swrcp, gravel, width, LOGFATAL), + "@ generic.c LogError" + ); + EXPECT_DOUBLE_EQ( + SWRC_SWCtoSWP(-1., swrc_type, swrcp, gravel, width, LOGWARN), + SW_MISSING + ); + + EXPECT_DEATH_IF_SUPPORTED( + SWRC_SWCtoSWP(1., swrc_type, swrcp, 1., width, LOGFATAL), + "@ generic.c LogError" + ); + EXPECT_DOUBLE_EQ( + SWRC_SWCtoSWP(1., swrc_type, swrcp, 1., width, LOGWARN), + SW_MISSING + ); + + EXPECT_DEATH_IF_SUPPORTED( + SWRC_SWCtoSWP(1., swrc_type, swrcp, gravel, 0., LOGFATAL), + "@ generic.c LogError" + ); + EXPECT_DOUBLE_EQ( + SWRC_SWCtoSWP(1., swrc_type, swrcp, gravel, 0., LOGWARN), + SW_MISSING + ); + } + + // --- *) if (theta - theta_res) < 0 (specific to vanGenuchten1980) + // note: this case is normally prevented due to SWC checks + swrc_type = encode_str2swrc((char *) "vanGenuchten1980"); + memset(swrcp, 0., SWRC_PARAM_NMAX * sizeof(swrcp[0])); + swrcp[0] = 0.1246; + swrcp[1] = 0.4445; + swrcp[2] = 0.0112; + swrcp[3] = 1.2673; + swrcp[4] = 7.78506; - // if theta1 == 0 (i.e., gravel == 1) && lyr->bMatric == 0: - // would be division by NaN - // note: this case is in normally prevented due to checks of inputs by - // function `water_eqn` for `bMatric` and function `_read_layers` for - // `gravelFraction` - help = SW_Site.lyr[n]->bMatric; - SW_Site.lyr[n]->bMatric = 0.; EXPECT_DEATH_IF_SUPPORTED( - SW_SWCbulk2SWPmatric(1., SW_Site.lyr[n]->swcBulk_fieldcap, n), + SWRC_SWCtoSWP(0.99 * swrcp[0], swrc_type, swrcp, gravel, width, LOGFATAL), "@ generic.c LogError" ); - SW_Site.lyr[n]->bMatric = help; - - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); + EXPECT_DOUBLE_EQ( + SWRC_SWCtoSWP(0.99 * swrcp[0], swrc_type, swrcp, gravel, width, LOGWARN), + SW_MISSING + ); } - // Test the 'SW_SoilWater' function 'SW_SWPmatric2VWCBulk' - TEST(SWSoilWaterTest, SWPmatric2VWCBulk){ + + // Death Tests of 'SW_SoilWater' function 'SWRC_SWPtoSWC' + TEST(SoilWaterDeathTest, SWPtoSWC) { // set up mock variables - RealD fractionGravel = .1, swpMatric = 15.0, p = 0.11656662532982573, - psisMatric = 18.608013, binverseMatric = 0.188608, thetaMatric = 41.37; - RealD tExpect, t, actualExpectDiff; - int i; - LyrIndex n = 0; - SW_Site.lyr[n]->thetasMatric = thetaMatric; - SW_Site.lyr[n]->psisMatric = psisMatric; - SW_Site.lyr[n]->bMatric = binverseMatric; - - // set gravel fractions on the interval [.0, .8], step .05 - for (i = 0; i <= 16; i++){ - fractionGravel = i / 20.; - tExpect = p * (1 - fractionGravel); - - t = SW_SWPmatric2VWCBulk(fractionGravel, swpMatric, n); - actualExpectDiff = fabs(t - tExpect); - - // when fractionGravel is between [.0, .8], we expect t = p * (1 - fractionGravel) - EXPECT_LT(actualExpectDiff, tol6); + RealD + swrcp[SWRC_PARAM_NMAX], + gravel = 0.1, + width = 10.; - } - Reset_SOILWAT2_after_UnitTest(); + unsigned int swrc_type; - fractionGravel = 1; - t = SW_SWPmatric2VWCBulk(fractionGravel, swpMatric, n); - // when fractionGravel is 1, we expect t == 0 - EXPECT_EQ(t, 0); - // Reset to previous global states - Reset_SOILWAT2_after_UnitTest(); + //--- we expect (non-)fatal errors in two situations + // (fatality depends on the error mode) + + //--- 1) Unimplemented SWRC + swrc_type = N_SWRCs + 1; + EXPECT_DEATH_IF_SUPPORTED( + SWRC_SWPtoSWC(15., swrc_type, swrcp, gravel, width, LOGFATAL), + "@ generic.c LogError" + ); + EXPECT_DOUBLE_EQ( + SWRC_SWPtoSWC(15., swrc_type, swrcp, gravel, width, LOGWARN), + SW_MISSING + ); + + // --- 2) swp < 0: water content cannot be negative (any SWRC) + for (swrc_type = 0; swrc_type < N_SWRCs; swrc_type++) { + EXPECT_DEATH_IF_SUPPORTED( + SWRC_SWPtoSWC(-1., swrc_type, swrcp, gravel, width, LOGFATAL), + "@ generic.c LogError" + ); + EXPECT_DOUBLE_EQ( + SWRC_SWPtoSWC(-1., swrc_type, swrcp, gravel, width, LOGWARN), + SW_MISSING + ); + } } } diff --git a/test/test_WaterBalance.cc b/test/test_WaterBalance.cc index c1e8f6657..4472c6917 100644 --- a/test/test_WaterBalance.cc +++ b/test/test_WaterBalance.cc @@ -222,4 +222,72 @@ namespace { Reset_SOILWAT2_after_UnitTest(); } + TEST(WaterBalanceTest, WithSWRCvanGenuchten1980) { + int i; + + // Set SWRC and PTF (and SWRC parameter input filename) + strcpy(SW_Site.site_swrc_name, (char *) "vanGenuchten1980"); + SW_Site.site_swrc_type = encode_str2swrc(SW_Site.site_swrc_name); + strcpy(SW_Site.site_ptf_name, (char *) "Rosetta3"); + SW_Site.site_ptf_type = encode_str2ptf(SW_Site.site_ptf_name); + SW_Site.site_has_swrcp = swTRUE; + + Mem_Free(InFiles[eSWRCp]); + InFiles[eSWRCp] = Str_Dup("Input/swrc_params_vanGenuchten1980.in"); + + // Read SWRC parameter input file (which is not read by default) + SW_SWRC_read(); + + // Update soils + SW_SIT_init_run(); + + // Run the simulation + SW_CTL_main(); + + // Collect and output from daily checks + for (i = 0; i < N_WBCHECKS; i++) { + EXPECT_EQ(0, SW_Soilwat.wbError[i]) << + "Water balance error in test " << + i << ": " << (char*)SW_Soilwat.wbErrorNames[i]; + } + + // Reset to previous global state + Reset_SOILWAT2_after_UnitTest(); + } + + + + TEST(WaterBalanceTest, WithSWRCFXW) { + int i; + + // Set SWRC and PTF (and SWRC parameter input filename) + strcpy(SW_Site.site_swrc_name, (char *) "FXW"); + SW_Site.site_swrc_type = encode_str2swrc(SW_Site.site_swrc_name); + strcpy(SW_Site.site_ptf_name, (char *) "neuroFX2021"); + SW_Site.site_ptf_type = encode_str2ptf(SW_Site.site_ptf_name); + SW_Site.site_has_swrcp = swTRUE; + + Mem_Free(InFiles[eSWRCp]); + InFiles[eSWRCp] = Str_Dup("Input/swrc_params_FXW.in"); + + // Read SWRC parameter input file (which is not read by default) + SW_SWRC_read(); + + // Update soils + SW_SIT_init_run(); + + // Run the simulation + SW_CTL_main(); + + // Collect and output from daily checks + for (i = 0; i < N_WBCHECKS; i++) { + EXPECT_EQ(0, SW_Soilwat.wbError[i]) << + "Water balance error in test " << + i << ": " << (char*)SW_Soilwat.wbErrorNames[i]; + } + + // Reset to previous global state + Reset_SOILWAT2_after_UnitTest(); + } + } // namespace diff --git a/testing/Input/siteparam.in b/testing/Input/siteparam.in index e019ee2bf..bce9a0557 100644 --- a/testing/Input/siteparam.in +++ b/testing/Input/siteparam.in @@ -9,9 +9,11 @@ #---- Soil water content initialization, minimum, and wet condition --1.0 # swc_min : cm/cm if 0 - <1.0, -bars if >= 1.0.; if < 0. then estimate residual water content for each layer -15.0 # swc_init: cm/cm if < 1.0, -bars if >= 1.0. -15.0 # swc_wet : cm/cm if < 1.0, -bars if >= 1.0. +-1.0 # swc_min : [cm/cm] if >= 0. and < 1.0 + # [-bars] if >= 1.0. + # estimate (from realistic limit and Rawls et al. 1985) if < 0. +15.0 # swc_init: [cm/cm] if < 1.0, [-bars] if >= 1.0. +15.0 # swc_wet : [cm/cm] if < 1.0, [-bars] if >= 1.0. #---- Diffuse recharge and runoff/runon 0 # reset (1/0): do/don't reset soil water content for each year @@ -83,10 +85,34 @@ NAN # aspect = surface azimuth angle (degrees): S=0, E=-90, N=±180, W=90; # Name of CO2 scenario: see input file `carbon.in` RCP85 + # --- Soil characterization --- # Are inputs of density representing bulk soil (type 1) or the matric component (type 0)? 0 + +#--- Soil water retention curve (SWRC) ------ +# +# Implemented options (`swrc_name`/`ptf_name`, see `swrc2str[]`/`ptf2str[]`): +# - ptf_name = : SWRC parameters must be provided via "swrc_params.in" +# - swrc_name = "Campbell1974" (Campbell 1974) +# * ptf_name = "Cosby1984AndOthers" (Cosby et al. 1984 but `swc_sat` by Saxton et al. 2006) +# * ptf_name = "Cosby1984" (Cosby et al. 1984) +# - swrc_name = "vanGenuchten1980" (van Genuchten 1980) +# - swrc_name = "FXW" (Fredlund and Xing 1994, Wang et al. 2018) +# +# Note: option "Campbell1974"/"Cosby1984AndOthers" was hard-coded < v7.0.0 +# Note: `rSOILWAT2` may implement additional PTFs + +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: Estimate with specified pedotransfer function + # 1: Use values from "swrc_params.in" + + #---- Transpiration regions # ndx : 1=shallow, 2=medium, 3=deep, 4=very deep # layer: deepest soil layer number of the region. diff --git a/testing/Input/swrc_params.in b/testing/Input/swrc_params.in new file mode 100644 index 000000000..f05fb9e02 --- /dev/null +++ b/testing/Input/swrc_params.in @@ -0,0 +1,39 @@ +#------ 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` +# - 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) + +# swrc = "Campbell1974" (default values below, from "Cosby1984") +# * param1 = air-entry suction [cm] +# * param2 = saturated volumetric water content for the matric component [cm/cm] +# * param3 = b, slope of the linear log-log retention curve [-] +# * param4 = saturated hydraulic conductivity [cm/day] + +# swrc = "vanGenuchten1980" +# * param1 = residual volumetric water content for the matric component [cm/cm] +# * param2 = saturated volumetric water content for the matric component [cm/cm] +# * param3 = alpha, related to the inverse of air entry suction [cm-1] +# * param4 = n, measure of the pore-size distribution [-] +# * param5 = saturated hydraulic conductivity [cm/day] + +# swrc = "FXW" +# * param1 = saturated volumetric water content of the matric component [cm/cm] +# * param2 = alpha, shape parameter [cm-1] +# * param3 = n, shape parameter [-] +# * param4 = m, shape parameter [-] +# * param5 = saturated hydraulic conductivity [cm / day] +# * param6 = L, tortuosity/connectivity parameter [-] + + +# param1 param2 param3 param4 param5 param6 +18.6080 0.42703 5.3020 24.03047 0.0000 0.0000 +20.4644 0.43290 7.0500 14.94351 0.0000 0.0000 +22.8402 0.44013 9.4320 13.71177 0.0000 0.0000 +24.0381 0.44291 10.0690 13.43171 0.0000 0.0000 +24.2159 0.44359 10.3860 14.14351 0.0000 0.0000 +23.3507 0.44217 10.3830 40.63764 0.0000 0.0000 +12.3880 0.41370 7.3250 37.22899 0.0000 0.0000 +12.3880 0.41370 7.3250 37.22899 0.0000 0.0000 diff --git a/testing/Input/swrc_params_FXW.in b/testing/Input/swrc_params_FXW.in new file mode 100644 index 000000000..7f1bc8df3 --- /dev/null +++ b/testing/Input/swrc_params_FXW.in @@ -0,0 +1,26 @@ +#------ 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` +# - 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) + +# swrc = "FXW" (values below, from "neuroFX2021") +# * param1 = saturated volumetric water content of the matric component [cm/cm] +# * param2 = alpha, shape parameter [cm-1] +# * param3 = n, shape parameter [-] +# * param4 = m, shape parameter [-] +# * param5 = saturated hydraulic conductivity [cm / day] +# * param6 = L, tortuosity/connectivity parameter [-] + + +# 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 +0.471163 0.149055 1.143810 0.124494 332.262496 2.988864 +0.475940 0.153117 1.141559 0.112295 420.418728 3.012669 +0.480690 0.157887 1.142653 0.105748 534.172981 3.049937 +0.538088 0.174184 1.124589 0.098441 978.516197 3.287010 +0.453070 0.169900 1.308269 0.182713 672.009929 3.218662 +0.453070 0.169900 1.308269 0.182713 672.009929 3.218662 diff --git a/testing/Input/swrc_params_vanGenuchten1980.in b/testing/Input/swrc_params_vanGenuchten1980.in new file mode 100644 index 000000000..d31396bcc --- /dev/null +++ b/testing/Input/swrc_params_vanGenuchten1980.in @@ -0,0 +1,24 @@ +#------ 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` +# - 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) + +# swrc = "vanGenuchten1980" (values below, from "Rosetta3") +# * param1 = residual volumetric water content for the matric component [cm/cm] +# * param2 = saturated volumetric water content for the matric component [cm/cm] +# * param3 = alpha, related to the inverse of air entry suction [cm-1] +# * param4 = n, measure of the pore-size distribution [-] +# * param5 = saturated hydraulic conductivity [cm/day] + +# 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 +0.12060752 0.4278807 0.010424896 1.287923 7.862807 0 +0.12336711 0.4393192 0.010807529 1.274654 9.100139 0 +0.12461498 0.4444546 0.011155912 1.267327 9.902011 0 +0.12480807 0.4426857 0.011408809 1.264873 9.784818 0 +0.10129327 0.3878319 0.014241212 1.311722 10.985124 0 +0.10129327 0.3878319 0.014241212 1.311722 10.985124 0 diff --git a/testing/files.in b/testing/files.in index 53174483c..fe3aa484e 100755 --- a/testing/files.in +++ b/testing/files.in @@ -11,6 +11,7 @@ Output/logfile.log # Output file to which warnings, errors, and other important #--- Description of simulated site Input/siteparam.in # Input file for site location, initialization, and miscellaneous model parameters Input/soils.in # Input for soil information: soil layer, soil texture, etc. +Input/swrc_params.in # Input for soil water retention curve (used if pdf_type = 0, i.e., pedotransfer functions are not used) #--- Inputs of weather forcing and description of climate conditions Input/weathsetup.in # Input file for weather-related parameters and weather generator settings