Feature soil organic content

NEWS.md

@@ -1,16 +1,49 @@
 # NEWS
 
 # SOILWAT2 v8.1.0-devel
-* This version produces nearly identical simulation output as previously.
-  Small deviations arise from replacing all remaining variables of type float
-  with type double (see commit 62237ae on 2024-July-30).
+* This version produces similar but not identical simulation output
+  as previously because of the following changes:
+    * Small deviations arise from replacing all remaining variables of
+      type float with type double (see commit 62237ae on 2024-July-30).
+    * Saturated percolation is now limited which leads to different outcomes
+      during periods of high infiltration (e.g., snowmelt) and during conditions
+      of low hydraulic conductivity (e.g., frozen soils, sapric organic matter).
 
 * The two models of SOILWAT2 can now be compiled with the following flags:
     * `make CPPFLAGS=-DSWTXT` (or as previously `make all`) for txt-based
     * `make CPPFLAGS=-DSWNC` for nc-based SOILWAT2.
 
 * Tests now require `c++17` and utilize `googletest` `v1.15.2` (issue #427).
 
+* SOILWAT2 can now represent the influence of soil organic matter on the
+  soil water retention curve and the saturated hydraulic conductivity
+  parameter (#397; @dschlaep). The implemented approach first determines
+  organic matter properties for the soil layers assuming fibric peat
+  characteristics at the soil surface and characteristics of sapric peat
+  at a user-specified depth. Then, bulk soil parameters of the soil water
+  retention curve are estimated as linear combinations of properties for
+  the mineral soil component and of properties for the organic matter soil
+  component using the proportion of organic matter in the bulk soil as weights.
+  The bulk soil saturated hydraulic conductivity parameter accounts for
+  flow pathways through organic matter above a threshold and assumes
+  conductivities through mineral and organic components in series outside of
+  those pathways.
+
+* Saturated percolation is now limited. The upper bound is a function based on
+  the saturated hydraulic conductivity parameter
+  (which includes effects of organic matter), frozen soils, and a
+  user-specified `"permeability"` factor.
+
+## Changes to inputs
+* New input via `"siteparam.in"` to specify the depth at which characteristics
+  of organic matter have completely switched from fibric to sapric peat
+  (default is 50 cm).
+* New input via `"soils.in"` to provide the proportion of soil organic matter
+  to bulk soil by weight for each soil layer.
+* New input via `"swrc_params*.in"` to provide parameter values of the
+  soil water retention curve representing fibric and sapric peat.
+  Note: Some parameter values for the `"FXW"` SWRC are missing.
+
 
 # SOILWAT2 v8.0.0
 * Simulation output remains the same as the previous version.

include/SW_Flow_lib.h

@@ -95,6 +95,7 @@ void infiltrate_water_high(
     unsigned int nlyrs,
     const double swcfc[],
     double swcsat[],
+    double ksat[],
     const double impermeability[],
     double *standingWater,
     double lyrFrozen[]

include/SW_Site.h

@@ -205,7 +205,7 @@ unsigned int encode_str2ptf(char *ptf_name);
 
 void SWRC_PTF_estimate_parameters(
     unsigned int ptf_type,
-    double *swrcp,
+    double *swrcpMineralSoil,
     double sand,
     double clay,
     double gravel,
@@ -214,7 +214,7 @@ void SWRC_PTF_estimate_parameters(
 );
 
 void SWRC_PTF_Cosby1984_for_Campbell1974(
-    double *swrcp, double sand, double clay
+    double *swrcpMineralSoil, double sand, double clay
 );
 
 
@@ -240,6 +240,7 @@ double SW_swcBulk_saturated(
     unsigned int ptf_type,
     double sand,
     double clay,
+    double fom,
     LOG_INFO *LogInfo
 );
 
@@ -252,23 +253,36 @@ double SW_swcBulk_minimum(
     double ui_sm_min,
     double sand,
     double clay,
+    double fom,
     double swcBulk_sat,
     double SWCMinVal,
     LOG_INFO *LogInfo
 );
 
+
 void PTF_Saxton2006(
-    double *theta_sat, double sand, double clay, LOG_INFO *LogInfo
+    double *theta_sat, double sand, double clay, double fom, LOG_INFO *LogInfo
 );
 
 void PTF_RawlsBrakensiek1985(
     double *theta_min,
     double sand,
     double clay,
+    double fom,
     double porosity,
     LOG_INFO *LogInfo
 );
 
+void SWRC_bulkSoilParameters(
+    unsigned int swrc_type,
+    double *swrcp,
+    const double *swrcpMS,
+    double swrcpOM[][SWRC_PARAM_NMAX],
+    double fom,
+    double depthSapric,
+    double depthT,
+    double depthB
+);
 
 double calculate_soilBulkDensity(double matricDensity, double fractionGravel);
 
@@ -322,6 +336,7 @@ void set_soillayers(
     const double *pclay,
     const double *imperm,
     const double *soiltemp,
+    const double *pom,
     int nRegions,
     double *regionLowerBounds,
     LOG_INFO *LogInfo

include/SW_datastructs.h

@@ -299,8 +299,9 @@ typedef struct {
 
     char site_swrc_name[64], site_ptf_name[64];
 
-    Bool site_has_swrcp; /**< Are `swrcp` already (TRUE) or not yet estimated
-                            (FALSE)? */
+    /** Are `swrcp` of the mineral soil already (TRUE) or not yet estimated
+        (FALSE)? */
+    Bool site_has_swrcpMineralSoil;
 
     /* transpiration regions  shallow, moderately shallow,  */
     /* deep and very deep. units are in layer numbers. */
@@ -309,10 +310,13 @@ typedef struct {
         SWCWetVal,     /* value for a "wet" day,       */
         SWCMinVal;     /* lower bound on swc.          */
 
-    /* bulk = relating to the whole soil, i.e., matric + rock/gravel/coarse
-     * fragments */
-    /* matric = relating to the < 2 mm fraction of the soil, i.e., sand, clay,
-     * and silt */
+    /* Soil components
+            * bulk = relating to the whole soil
+              i.e., matric + coarse fragment (gravel)
+            * matric component = relating to the < 2 mm fraction of the soil
+            * mineral component = sand, clay, silt
+            * organic component = organic matter
+     */
 
     double
         /* Inputs */
@@ -335,6 +339,8 @@ typedef struct {
                                        (0=permeable, 1=impermeable)    */
         avgLyrTempInit[MAX_LAYERS], /* initial soil temperature for each soil
                                        layer */
+        /** Organic matter content as weight fraction of bulk soil [g g-1] */
+        fractionWeight_om[MAX_LAYERS],
 
         /* Derived soil characteristics */
         soilMatric_density[MAX_LAYERS], /* matric soil density of the < 2 mm
@@ -360,14 +366,21 @@ typedef struct {
 
         /* Saxton et al. 2006 */
         swcBulk_saturated[MAX_LAYERS]; /* saturated bulk SWC [cm] */
-                                       // currently, not used;
+
+    /** Saturated hydraulic conductivity of the bulk soil */
+    double ksat[MAX_LAYERS];
+
+    // currently, not used;
     // Saxton2006_K_sat_matric, /* saturated matric conductivity [cm / day] */
     // Saxton2006_K_sat_bulk, /* saturated bulk conductivity [cm / day] */
     // Saxton2006_fK_gravel, /* gravel-correction factor for conductivity [1] */
     // Saxton2006_lambda; /* Slope of logarithmic tension-moisture curve */
 
     double depths[MAX_LAYERS]; // soil layer depths of SoilWat soil
 
+    /** Depth [cm] at which soil properties reach values of sapric peat */
+    double depthSapric;
+
     /* Soil water retention curve (SWRC) */
     unsigned int swrc_type[MAX_LAYERS], /**< Type of SWRC (see #swrc2str) */
         ptf_type[MAX_LAYERS];           /**< Type of PTF (see #ptf2str) */
@@ -377,9 +390,19 @@ typedef struct {
                      `swrcp` but we need to loop over soil layers for every
                      vegetation type in `my_transp_rng`
     */
+    /** SWRC parameters of the bulk soil
+        (weighted average of mineral and organic SWRC).
+
+        Note: parameter interpretation varies with selected SWRC,
+        see `SWRC_check_parameters()`
+    */
     double swrcp[MAX_LAYERS][SWRC_PARAM_NMAX];
-    /**< Parameters of SWRC: parameter interpretation varies with selected SWRC,
-     * see `SWRC_check_parameters()` */
+    /** SWRC parameters of the mineral soil component */
+    double swrcpMineralSoil[MAX_LAYERS][SWRC_PARAM_NMAX];
+    /** SWRC parameters of the organic soil component
+        for (1) fibric and (2) sapric peat. */
+    double swrcpOM[2][SWRC_PARAM_NMAX];
+
     LyrIndex my_transp_rgn[NVEGTYPES][MAX_LAYERS]; /* which transp zones from
                                                       Site am I in? */
 

src/SW_Flow.c

@@ -486,6 +486,7 @@ void SW_Water_Flow(SW_RUN *sw, LOG_INFO *LogInfo) {
             n_layers,
             sw->Site.swcBulk_fieldcap,
             sw->Site.swcBulk_saturated,
+            sw->Site.ksat,
             sw->Site.impermeability,
             &UpNeigh_standingWater,
             sw->SoilWat.lyrFrozen
@@ -517,6 +518,7 @@ void SW_Water_Flow(SW_RUN *sw, LOG_INFO *LogInfo) {
         n_layers,
         sw->Site.swcBulk_fieldcap,
         sw->Site.swcBulk_saturated,
+        sw->Site.ksat,
         sw->Site.impermeability,
         standingWaterToday,
         sw->SoilWat.lyrFrozen

src/SW_Flow_lib.c

@@ -376,6 +376,7 @@ void litter_intercepted_water(
 (cm H<SUB>2</SUB>O).
 @param swcsat Soilwater content in each layer at saturation
     (m<SUP>3</SUP>H<SUB>2</SUB>O).
+@param ksat Saturated hydraulic conductivity [cm day-1]
 @param impermeability Impermeability measures for each layer.
 @param *standingWater Remaining water on the surface
     (m<SUP>3</SUP>H<SUB>2</SUB>O).
@@ -404,14 +405,14 @@ void infiltrate_water_high(
     unsigned int nlyrs,
     const double swcfc[],
     double swcsat[],
+    double ksat[],
     const double impermeability[],
     double *standingWater,
     double lyrFrozen[]
 ) {
 
     unsigned int i;
     unsigned int j;
-    double d[MAX_LAYERS] = {0};
     double push;
     double ksat_rel;
 
@@ -427,19 +428,21 @@ void infiltrate_water_high(
             ksat_rel = 1.;
         }
 
+        ksat_rel *= (1. - impermeability[i]);
+
         /* calculate potential saturated percolation */
-        d[i] =
-            fmax(0., ksat_rel * (1. - impermeability[i]) * (swc[i] - swcfc[i]));
-        drain[i] = d[i];
+        drain[i] = ksat_rel * fmin(ksat[i], swc[i] - swcfc[i]);
+        drain[i] = fmax(0., drain[i]);
+
 
         if (i < nlyrs - 1) {
             /* percolate up to next-to-last layer */
-            swc[i + 1] += d[i];
-            swc[i] -= d[i];
+            swc[i + 1] += drain[i];
+            swc[i] -= drain[i];
         } else {
             /* percolate last layer */
-            (*drainout) = d[i];
-            swc[i] -= (*drainout);
+            (*drainout) = drain[i];
+            swc[i] -= drain[i];
         }
     }