SIP287 Add C Saturation

docs/CHANGELOG.md

@@ -30,6 +30,7 @@ sections to include in release notes:
 - `leafon` and `leafoff` events for tracking phenological transitions (#326)
 - `leafon` limited by available carbon and nitrogen; N storage pool; N resorption on `leafoff` (#337)
 - New required parameter `leafOnReallocFrac` to control how much of wood and coarse root carbon is reallocated to leaves on `leafon` (#337)
+- `carbonSaturation` flag and calculation of soil and litter carbon pools that observes soil carbon saturating behavior (#301)
 
 ### Fixed
 

docs/model-structure.md

@@ -309,6 +309,13 @@ F^C_\text{fert,org}
 Where $K_{\text{plant},i}$ is the turnover rate of plant pool $i$ that controls the rate at which plant biomass is
 transferred to litter.
 
+When soil carbon saturation is enabled, a fraction of soil carbon inputs may be redirected to the litter pool as fast-turnover carbon. This functionality is described in more detail below in the Soil Carbon section \eqref{eq:soil_carbon_saturation}.
+
+\begin{equation}
+\frac{dC_\text{litter}}{dt} = F^C_\text{litter} + F^C_{\text{soil}} \cdot \frac{C_{\text{soil}}}{C_{\text{soil,saturation}}} - F^C_{\text{decomp}} - F^C_{\text{CH}_4\text{,litter}}
+\label{eq:soil_carbon_to_litter}
+\end{equation}
+
 $F^C_{\text{decomp}}$ represents the rate at which litter carbon is processed by microbial activity. Litter
 decomposition
 is modeled as a first-order process proportional to litter carbon content and modified by temperature and moisture:
@@ -370,6 +377,15 @@ F^C_{\text{soil}} = F^C_{\text{soil,litter}} + F^C_{\text{soil,roots}}
 \label{eq:soil_carbon_flux}
 \end{equation}
 
+When soil carbon saturation is enabled, only a saturation-dependent fraction of gross soil C inputs is added to the soil pool. This fraction declines as $C_{\text{soil}}$ approaches the specified soil C saturation limit. The remaining input C is redirected to the litter pool \eqref{eq:soil_carbon_to_litter} as fast-turnover carbon rather than being added to the soil pool.
+
+\begin{equation}
+\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} \cdot (1 - \frac{C_{\text{soil}}}{C_{\text{soil,saturation}}}) - R_{\text{soil}} - F^C_{\text{CH}_4\text{,soil}}
+\label{eq:soil_carbon_saturation}
+\end{equation}
+
+where $C_{\text{soil,saturation}}$ is the soil carbon saturation limit entered as an input paramter. This is based on equation (3) from Stewart et al. (2007).
+
 Soil heterotrophic respiration is modeled as a first-order process proportional
 to soil organic carbon content and modified by environmental and management factors:
 

docs/parameters.md

@@ -197,6 +197,7 @@ Run-time parameters can change from one run to the next, or when the model is st
 | $Q_{10s}$             | soilRespQ10         | Soil respiration Q10                                                                | unitless                                                                            | scalar determining effect of temp on soil respiration          |
 | $D_{\text{moisture}}$ | soilRespMoistEffect | scalar determining effect of moisture on soil resp.                                 | unitless                                                                            |                                                                |
 | $f_{\text{till}}$     | tillageEff          | Effect of tillage on decomposition that exponentially decays over time              | fraction                                                                            | Documented in model structure; event-level term in `events.in` |
+| $C_{\text{soil,saturation}}$ | soilCSaturation | Maximum amount of carbon that can be stabilized in the soil determined by soil texture        | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$                                                                                   |                                                                |
 
 ### Nitrogen Cycle Parameters
 

docs/user-guide/model-inputs.md

@@ -242,11 +242,13 @@ Thus, command-line arguments override settings in the configuration file, and co
 | `snow`           | on      | Keep track of snowpack, rather than assuming all precipitation is liquid                |
 | `soil-phenol`    | off     | Use soil temperature to determine leaf growth                                           |
 | `water-hresp`    | on      | Whether soil moisture affects heterotrophic respiration                                 |
+| `carbon-saturation`| off   | Enable soil carbon saturation behavior to constrain carbon stored in soil               |
 
 Note the following restrictions on these options:
  - `soil-phenol` and `gdd` may not both be turned on
  - `anaerobic` requires `water-hresp`
  - 'nitrogen-cycle' requires both 'litter-pool' and 'anaerobic'
+ - `carbon-saturation` requires `litter-pool`
 
 ### Command Line Arguments
 

docs/user-guide/running-sipnet.md

@@ -48,6 +48,7 @@ These flags enable or disable optional model processes. Prepend `no-` to the fla
 | `--snow`           | ON (1)  | Track snowpack separately; if disabled, all precipitation is treated as liquid |
 | `--soil-phenol`    | OFF (0) | Use soil temperature (instead of growing degree days) to determine leaf growth |
 | `--water-hresp`    | ON (1)  | Allow soil moisture to affect heterotrophic respiration rates                  |
+| `--carbon-saturation`| OFF (0) | Enable soil carbon saturation behavior to constrain carbon stored in soil    |
 
 #### Model Flag Restrictions
 
@@ -56,6 +57,7 @@ The following flag constraints are enforced:
 - `--soil-phenol` and `--gdd` cannot both be enabled
 - `--anaerobic` requires `--water-hresp`
 - `--nitrogen-cycle` requires both `--litter-pool` and `--anaerobic`
+- `--carbon-saturation` requires `--litter-pool`
 
 ### Output Flags
 

src/common/context.c

@@ -43,6 +43,7 @@ void initContext(void) {
   CREATE_INT_CONTEXT(nitrogenCycle,   "NITROGEN_CYCLE",   ARG_OFF, FLAG_YES);
   CREATE_INT_CONTEXT(anaerobic,       "ANAEROBIC",        ARG_OFF, FLAG_YES);
   CREATE_INT_CONTEXT(flooding,        "FLOODING",         ARG_OFF, FLAG_YES);
+  CREATE_INT_CONTEXT(carbonSaturation,"CARBON_SATURATION",ARG_OFF, FLAG_YES);
 
   // Flags, I/O
   CREATE_INT_CONTEXT(doMainOutput,    "DO_MAIN_OUTPUT",   ARG_ON,  FLAG_YES);
@@ -209,6 +210,12 @@ void validateContext(void) {
     hasError = 1;
   }
 
+  if (ctx.carbonSaturation && !ctx.litterPool) {
+    logError("carbon-saturation requires litter-pool to be "
+             "turned on\n");
+    hasError = 1;
+  }
+
   if (hasError) {
     exit(EXIT_CODE_BAD_PARAMETER_VALUE);
   }

src/common/context.h

@@ -36,7 +36,7 @@ struct context_metadata {
   UT_hash_handle hh;  // makes this structure hashable
 };
 
-#define NUM_CONTEXT_MODEL_FLAGS 11
+#define NUM_CONTEXT_MODEL_FLAGS 12
 // See docs/developer-guide/cli-options.md for details on how to add a new
 // Context entry
 struct Context {
@@ -53,7 +53,7 @@ struct Context {
   int nitrogenCycle;
   int anaerobic;
   int flooding;
-
+  int carbonSaturation;
   // IF ADDING A NEW MODEL FLAG, update NUM_CONTEXT_MODEL_FLAGS above and
   // relevant code in restart.c
 

src/sipnet/cli.c

@@ -32,6 +32,7 @@ static struct option long_options[] = {  // NOLINT
     DECLARE_FLAG(nitrogen-cycle),
     DECLARE_FLAG(anaerobic),
     DECLARE_FLAG(flooding),
+    DECLARE_FLAG(carbon-saturation),
 
     DECLARE_FLAG(do-main-output),
     DECLARE_FLAG(do-single-outputs),
@@ -65,7 +66,7 @@ char *argNameMap[] = {
     DECLARE_ARG_FOR_MAP(litterPool), DECLARE_ARG_FOR_MAP(snow),
     DECLARE_ARG_FOR_MAP(soilPhenol), DECLARE_ARG_FOR_MAP(waterHResp),
     DECLARE_ARG_FOR_MAP(nitrogenCycle), DECLARE_ARG_FOR_MAP(anaerobic),
-    DECLARE_ARG_FOR_MAP(flooding),
+    DECLARE_ARG_FOR_MAP(flooding), DECLARE_ARG_FOR_MAP(carbonSaturation),
 
     // I/O
     DECLARE_ARG_FOR_MAP(doMainOutput), DECLARE_ARG_FOR_MAP(doSingleOutputs),
@@ -98,6 +99,7 @@ void usage(char *progName) {
   printf("  --snow               Keep track of snowpack, rather than assuming all precipitation is liquid (1)\n");
   printf("  --soil-phenol        Use soil temperature to determine leaf growth (0)\n");
   printf("  --water-hresp        Whether soil moisture affects heterotrophic respiration (1)\n");
+  printf("  --carbon-saturation  Enable maximum storage limit of soil organic carbon (0)\n");
   printf("\n");
   printf("Output flags: (prepend flag with 'no-' to force off, eg '--no-print-header')\n");
   printf("  --do-main-output     Print time series of all output variables to <file-prefix>.out (1)\n");
@@ -121,6 +123,7 @@ void usage(char *progName) {
   printf(" --soil-phenol and --gdd may not both be turned on\n");
   printf(" --anaerobic requires --water-hresp\n");
   printf(" --nitrogen-cycle requires both --litter-pool and --anaerobic\n");
+  printf(" --carbon-saturation requires --litter-pool\n");
   // clang-format on
 }
 

src/sipnet/cli.h

@@ -14,7 +14,7 @@
 
 // The run-time option names do not match their corresponding fields in Context,
 // so we need a way to get from one to the other.
-#define NUM_FLAG_OPTIONS 16
+#define NUM_FLAG_OPTIONS 17
 extern char *argNameMap[2 * NUM_FLAG_OPTIONS];
 
 /*!

src/sipnet/nitrogen.c

@@ -51,22 +51,32 @@ static void calcNPoolFluxes(void) {
   double litterMin = fluxes.rLitter / litterCN;
   double soilMin = fluxes.rSoil / soilCN;
 
+  // Adding soil carbon saturation functionality so organic N fluxes to soil
+  // and litter are proportional to respective carbon fluxes dependent on
+  // soil carbon saturation
+  double soilNInputs = fluxes.litterToSoil / litterCN +
+                       fluxes.fineRootLoss / params.fineRootCN +
+                       fluxes.coarseRootLoss / params.woodCN;
+  // saturationFraction capped between zero and one
+  double saturationFraction =
+      (ctx.carbonSaturation)
+          ? (fmin(1.0, fmax(0.0, envi.soilC / params.soilCSaturation)))
+          : 0.0;
+
   // litter
   // The litter org N flux is determined by the carbon fluxes from wood and leaf
   // litter (modified by leaf N resorption), and N loss due to mineralization.
   // N added via fertilization is handled elsewhere.
-  fluxes.nOrgLitter = fluxes.leafLitter / params.leafCN -
-                      fluxes.leafOffNResorption +
-                      fluxes.woodLitter / params.woodCN - litterMin -
-                      fluxes.litterToSoil / litterCN;
+  fluxes.nOrgLitter =
+      fluxes.leafLitter / params.leafCN - fluxes.leafOffNResorption +
+      fluxes.woodLitter / params.woodCN - litterMin -
+      fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction);
 
   // soil
   // The soil org N flux is determined by the carbon flux from the litter pool,
   // carbon fluxes from roots, and N loss due to mineralization
   // (Note: woodCN is used for coarse roots)
-  fluxes.nOrgSoil = fluxes.litterToSoil / litterCN - soilMin +
-                    fluxes.fineRootLoss / params.fineRootCN +
-                    fluxes.coarseRootLoss / params.woodCN;
+  fluxes.nOrgSoil = soilNInputs * (1 - saturationFraction) - soilMin;
 
   // mineralization
   fluxes.nMin = litterMin + soilMin;

src/sipnet/restart.c

@@ -85,6 +85,7 @@ typedef struct RestartContextModelFlags {
   int nitrogenCycle;
   int anaerobic;
   int flooding;
+  int carbonSaturation;
 } RestartContextModelFlags;
 static RestartContextModelFlags modelFlags;
 
@@ -160,7 +161,8 @@ void initResetState(RestartState *state, MeanTracker *npp) {
   state->flagsPF[8] = (StateField){"flags.nitrogenCycle", FT_INT, &modelFlags.nitrogenCycle, 0};
   state->flagsPF[9] = (StateField){"flags.anaerobic",     FT_INT, &modelFlags.anaerobic,     0};
   state->flagsPF[10] = (StateField){"flags.flooding",     FT_INT, &modelFlags.flooding,      0};
-  state->flagsPF[11] = (StateField){"flags.invalid",      FT_INVALID, NULL, FIELD_INVALID};
+  state->flagsPF[11] = (StateField){"flags.carbonSaturation", FT_INT, &modelFlags.carbonSaturation, 0};
+  state->flagsPF[12] = (StateField){"flags.invalid",      FT_INVALID, NULL, FIELD_INVALID};
 
   state->boundaryPF[0] = (StateField){"boundary.year",    FT_INT,     &boundaryClimate.year,   0};
   state->boundaryPF[1] = (StateField){"boundary.day",     FT_INT,     &boundaryClimate.day,    0};
@@ -766,6 +768,7 @@ static void checkRestartContextCompatibility(void) {
   mismatch |= (ctx.nitrogenCycle != modelFlags.nitrogenCycle);
   mismatch |= (ctx.anaerobic != modelFlags.anaerobic);
   mismatch |= (ctx.flooding != modelFlags.flooding);
+  mismatch |= (ctx.carbonSaturation != modelFlags.carbonSaturation);
 
   if (mismatch) {
     logError("Restart context mismatch: model flags must match checkpoint "
@@ -889,6 +892,8 @@ void restartWriteCheckpoint(const char *restartOut, MeanTracker *meanNPP) {
   ++numFlagsSet;
   modelFlags.flooding = ctx.flooding;
   ++numFlagsSet;
+  modelFlags.carbonSaturation = ctx.carbonSaturation;
+  ++numFlagsSet;
   if (numFlagsSet != NUM_CONTEXT_MODEL_FLAGS) {
     logInternalError("Not all model flags set while writing checkpoint\n");
     exit(EXIT_CODE_INTERNAL_ERROR);

src/sipnet/sipnet.c

@@ -392,6 +392,9 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) {
   // Water drainage
   initializeOneModelParam(modelParams, "waterDrainFrac", &(params.waterDrainFrac), ctx.flooding);
 
+  // Soil carbon saturation
+  initializeOneModelParam(modelParams, "soilCSaturation", &(params.soilCSaturation), ctx.carbonSaturation);
+
   // NOLINTEND
   // clang-format on
 
@@ -1647,15 +1650,27 @@ void updateMainPools(void) {
  */
 void updatePoolsForSoil(void) {
   if (ctx.litterPool) {
+    double soilInputs =
+        fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil;
+    // Adding soil carbon saturation functionality
+    // soilFraction capped between zero and one
+    // if envi.soilC > params.soilCSaturation (due to initialization or events)
+    // saturationFraction will return 1 and no additional soil C inputs will be
+    // accepted to envi.soilC until losses due to respiration and methane reduce
+    // pool below saturation level, all inputs will be directed to envi.litterC
+    // pool while envi.soilC is above saturation level
+    double saturationFraction =
+        (ctx.carbonSaturation)
+            ? (fmin(1.0, fmax(0.0, envi.soilC / params.soilCSaturation)))
+            : 0.0;
     // :: from [2], litter model description
-    envi.litterC +=
-        (fluxes.woodLitter + fluxes.leafLitter - fluxes.litterToSoil -
-         fluxes.rLitter - fluxes.litterMethane) *
-        climate->length;
-
+    envi.litterC += (fluxes.woodLitter + fluxes.leafLitter +
+                     (soilInputs * saturationFraction) - fluxes.litterToSoil -
+                     fluxes.rLitter - fluxes.litterMethane) *
+                    climate->length;
     // from [2] and [3], litter and root terms respectively
-    envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss +
-                   fluxes.litterToSoil - fluxes.rSoil - fluxes.soilMethane) *
+    envi.soilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil -
+                   fluxes.soilMethane) *
                   climate->length;
   } else {
     // Normal pool (single pool, no microbes)

src/sipnet/state.h

@@ -397,6 +397,14 @@ typedef struct Parameters {
 
   // Relative methane production rate in the litter pool, in [0, 1), per day
   double litterMethaneRate;
+
+  // ******
+  // Soil carbon saturation
+  // ******
+
+  // Maximum threshold for stabilizing carbon in soil organic pool as
+  // slow-turnover pool units: g C * m^-2 ground area
+  double soilCSaturation;
 } Params;
 
 #define NUM_PARAMS (sizeof(Params) / sizeof(double))
@@ -628,6 +636,7 @@ typedef struct FluxVars {
   double soilMethane;
   // Methane produced from litter
   double litterMethane;
+
 } Fluxes;
 
 // Global var

tests/sipnet/test_modeling/Makefile

@@ -8,7 +8,7 @@ LDFLAGS=-L$(ROOT_DIR)/libs
 LDLIBS=-lsipnet -lsipnet_common -lm
 
 # List test files in this directory here
-TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c
+TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c testCarbonSaturation.c
 
 # The rest is boilerplate, likely copyable as is to a new test directory
 TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o)