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.
 
+If soil carbon saturation is enabled, the litter carbon pool has an additional input of carbon that was not stabilized for long-term storage by the soil. 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}}}{f_{\text{C,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}
 
+If soil carbon saturation is enabled, only a fraction of the total carbon input will be added to the soil. The amount added to the soil carbon pool is a function of the proximity of the pool to the soil carbon saturation limit. As the soil carbon pool aqcuires more carbon and approaches its saturation limit, the amount of carbon stored in the pool decreases. This represents the physical and chemical constraints of the mineralogical surfaces of soil to stabilize carbon for long-term storage. The remaining carbon that is not stabilized in the soil carbon pool due to saturation constraints will be returned to the litter carbon pool to act as fast-turnover carbon \eqref{eq:soil_carbon_to_litter}.
+
+\begin{equation}
+\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} \cdot (1 - \frac{C_{\text{soil}}}{f_{\text{C,soil,saturation}}}) - R_{\text{soil}} - F^C_{\text{CH}_4\text{,soil}}
+\label{eq:soil_carbon_saturation}
+\end{equation}
+
+where $f_{\text{C,soil,saturation}}$ is the soil carbon saturation limit entered as an input paramter.
+
 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` |
+| $f_{\text{C,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,6 +242,7 @@ 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

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/restart.c

@@ -87,6 +87,7 @@ typedef struct RestartContextModelFlags {
   int nitrogenCycle;
   int anaerobic;
   int flooding;
+  int carbonSaturation;
 } RestartContextModelFlags;
 static RestartContextModelFlags modelFlags;
 
@@ -162,7 +163,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};
@@ -777,6 +779,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 "
@@ -900,6 +903,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,16 +1650,30 @@ void updateMainPools(void) {
  */
 void updatePoolsForSoil(void) {
   if (ctx.litterPool) {
-    // :: from [2], litter model description
-    envi.litterC +=
-        (fluxes.woodLitter + fluxes.leafLitter - fluxes.litterToSoil -
-         fluxes.rLitter - fluxes.litterMethane) *
-        climate->length;
+    if (ctx.carbonSaturation) {
+      double saturationFraction = envi.soilC / params.soilCSaturation;
+      double soilInputs =
+          fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil;
+
+      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) *
-                  climate->length;
+      envi.soilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil -
+                     fluxes.soilMethane) *
+                    climate->length;
+    } else {
+      // :: from [2], litter model description
+      envi.litterC +=
+          (fluxes.woodLitter + fluxes.leafLitter - 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) *
+                    climate->length;
+    }
   } else {
     // Normal pool (single pool, no microbes)
     // :: from [1] (and others, TBD), eq (A3), where:

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)