SIP261 Implement methane production

docs/CHANGELOG.md

@@ -45,6 +45,7 @@ sections to include in release notes:
 - Split wood carbon pool into `plantWoodC` and `plantWoodCStorageDelta` to track NPP storage lag (#248)
 - New moisture dependency function option for soil respiration, controlled by the `anaerobic` cli option (#259)
 - New moisture dependency functions for N volatilization, methane (#259)
+- Methane production option, controlled by the `anaerobic` cli option (#269)
 
 ### Fixed
 

docs/model-structure.md

@@ -259,15 +259,17 @@ This is from Zobitz, et al. (2008), Sec. 5.4.2.
 
 ### Litter Carbon
 
-The change in the litter carbon pool over time is defined by the input of new litter and losses due to decomposition:
+The change in the litter carbon pool over time is defined by the input of new litter and losses due to decomposition 
+and methane production:
 
 \begin{equation}
-\frac{dC_\text{litter}}{dt} = F^C_\text{litter} - F^C_{\text{decomp}}
+\frac{dC_\text{litter}}{dt} = F^C_\text{litter} - F^C_{\text{decomp}} - F^C_\mathit{CH_4,litter}
 \end{equation}
 
-Where $F^C_\text{litter}$ is the carbon flux from aboveground plant biomass \eqref{eq:litter_flux} and
-$F^C_{\text{decomp}}$ is the total litter decomposition flux \eqref{eq:decomp_rate}. Note that belowground turnover
-is routed directly to the soil carbon pool (see Soil Carbon).
+Where $F^C_\text{litter}$ is the carbon flux from aboveground plant biomass \eqref{eq:litter_flux},
+$F^C_{\text{decomp}}$ is the total litter decomposition flux \eqref{eq:decomp_rate}, and $F^C_\mathit{CH_4,litter}$ is
+the methane flux from the litter \eqref{eq:ch4}. Note that belowground turnover is routed directly to the soil carbon
+pool (see Soil Carbon).
 
 $F^C_\text{litter}$ is the sum of litter produced through aboveground senescence, transfer of biomass during harvest,
 and organic matter amendments:
@@ -295,7 +297,7 @@ is modeled as a first-order process proportional to litter carbon content and mo
 
 \begin{equation}
 F^C_{\text{decomp}} =
-K_\text{litter} \cdot C_\text{litter} \cdot D_{\text{temp}} \cdot D_{\text{water}R_H} \cdot
+K_\text{litter} \cdot C_\text{litter} \cdot D_{\text{temp}} \cdot D_{\text{water},R_H} \cdot
 D_{CN} \cdot
 D_{tillage}
 \label{eq:decomp_rate}
@@ -329,16 +331,19 @@ F^C_{\text{soil,litter}} =
 
 The change in the soil organic carbon (SOC) pool over time is determined by: (i) inputs of carbon transferred from
 litter during decomposition, (ii) direct inputs from belowground plant turnover, and (iii) losses of carbon due to
-heterotrophic respiration:
+heterotrophic respiration and methane production:
 
 \begin{equation}
-\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} - R_{\text{soil}}
+\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} - R_{\text{soil}} - F^C_\mathit{CH_4,soil}
 \label{eq:Braswell_A3}
 \end{equation}
 
-This is equation (A3) from Braswell, et al. (2005).
+where $F^C_{\text{soil}}$ is the total carbon input to the soil, $R_{\text{soil}}$ is the soil heterotrophic 
+respiration, and $F^C_\mathit{CH_4,soil}$ is the methane flux from the soil. 
 
-Total carbon input to the soil, $F^C_{\text{soil}}$, includes both
+This is equation (A3) from Braswell, et al. (2005), with the addition of the methane flux.
+
+Total carbon input to the soil includes both
 (i) carbon transferred from the litter pool during decomposition \eqref{eq:soil_carbon} and
 (ii) inputs from root turnover:
 
@@ -347,7 +352,7 @@ F^C_{\text{soil}} = F^C_{\text{soil,litter}} + F^C_{\text{soil,roots}}
 \label{eq:soil_carbon_flux}
 \end{equation}
 
-Soil heterotrophic respiration, $R_{\text{soil}}$, is modeled as a first-order process proportional
+Soil heterotrophic respiration is modeled as a first-order process proportional
 to soil organic carbon content and modified by environmental and management factors:
 
 \begin{equation}
@@ -383,18 +388,18 @@ heterotrophic respiration is then defined as a fixed fraction of this decomposit
 (\eqref{eq:decomp_carbon}--\eqref{eq:r_litter}), with the remainder transferred to the soil carbon pool
 (\eqref{eq:soil_carbon}).
 
-### $\frak{Methane \ Production \ (C \rightarrow CH_4)}$
+### Methane Production $(C \rightarrow CH_4)$
 
 \begin{equation}
-F^C_\mathit{CH_4} = \left(\sum_{j} K_{CH_4,j} \cdot C_\text{j}\right) \cdot D_\mathrm{water, CH_4} \cdot D_\text{temp}
+F^C_\mathit{CH_4,j} = K_{CH_4,j} \cdot C_\text{j} \cdot D_\mathrm{water, CH_4} \cdot D_\text{temp}
 \label{eq:ch4}
 \end{equation}
 
 \begin{equation*}
 \small j \in \{\text{soil, litter}\}
 \end{equation*}
 
-The calculation of methane flux $(F^C_{CH_4})$ is analogous to that of $R_H$. It uses the same carbon pools as substrate
+The calculation of methane flux $(F^C_{CH_4})$ for soil and litter is analogous to that of $R_H$. It uses the same carbon pools as substrate
 and temperature dependence but has specific rate parameters $(K_{\mathit{CH_4,}j})$, a moisture dependence function 
 based on oxygen availability \eqref{eq:water_ch4}, and no direct dependence on tillage.
 

docs/parameters.md

@@ -223,7 +223,9 @@ pool stoichiometry.
 | $f_{\text{fastflow}}$     | fastFlowFrac        | fraction of water entering soil that goes directly to drainage                                                         | unitless                                       |                                                                                                                                                                                         |
 | $f_a$                     | fAnoxia             | Soil wetness fraction at which oxygen diffusion begins to limit aerobic respiration                                    | unitless                                       | Used in moisture partitioning between aerobic and anaerobic pathways (\eqref{eq:water_rh_2} in model structure)                                                                         |
 | $\eta$                    | anaerobicDecompRate | Relative anaerobic decomposition rate                                                                                  | unitless                                       | Active when anaerobic pathway is enabled; expected range $(0,1]$                                                                                                                        |
-| $p$                       | anaerobicTransExp   | Methane anoxia sensitivity exponent in $D_{\text{water},CH_4}=A^p$                                                     | unitless                                       |                                                                                                                                                                                         |
+| $p$                       | anaerobicTransExp   | Methane anoxia sensitivity exponent in $D_{\text{water},CH_4}=A^p$                                                     | unitless                                       | Expected range $\ge 1$                                                                                                                                                                  |
+| $K_{\text{CH_4,soil}}$    | soilMethaneRate     | First-order methane production rate constant for the soil C pool                                                       | $\text{day}^{-1}$                              | Expected range $\ge 0$                                                                                                                                                                  |
+| $K_{\text{CH_4,litter}}$  | litterMethaneRate   | First-order methane production rate constant for the litter C pool                                                     | $\text{day}^{-1}$                              | Expected range $\ge 0$                                                                                                                                                                  |
 | $M_{\text{snow}}$         | snowMelt            | rate at which snow melts                                                                                               | cm water equivalent per degree Celsius per day |                                                                                                                                                                                         |
 | $r_{d}$                   | rdConst             | scalar determining amount of aerodynamic resistance                                                                    | unitless                                       |                                                                                                                                                                                         |
 | $r_{s,1}$                 | rSoilConst1         |                                                                                                                        | unitless                                       | soil resistance = e^(rSoilConst1 - rSoilConst2 \* W1) , where W1 = (water/soilWHC)                                                                                                      |

src/sipnet/balance.c

@@ -45,9 +45,10 @@ void updateBalanceTrackerPostUpdate(void) {
   balanceTracker.outputsC = fluxes.rVeg + fluxes.rFineRoot +
                             fluxes.rCoarseRoot +  // R_a
                             fluxes.rSoil +  // R_h
+                            fluxes.soilMethane +  // methane
                             fluxes.eventOutputC;  // agro event removals
   if (ctx.litterPool) {
-    balanceTracker.outputsC += fluxes.rLitter;
+    balanceTracker.outputsC += fluxes.rLitter + fluxes.litterMethane;
   }
   // Account for climate length
   balanceTracker.inputsC *= climate->length;

src/sipnet/sipnet.c

@@ -384,8 +384,11 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) {
   // New moisture dependency params
   initializeOneModelParam(modelParams, "fAnoxia", &(params.fAnoxia), ctx.anaerobic || ctx.nitrogenCycle);
   initializeOneModelParam(modelParams, "anaerobicDecompRate", &(params.anaerobicDecompRate), ctx.anaerobic);
-  // TODO: this should depend on methane, when that gets in
-  initializeOneModelParam(modelParams, "anaerobicTransExp", &(params.anaerobicTransExp), 0);
+
+  // Methane
+  initializeOneModelParam(modelParams, "anaerobicTransExp", &(params.anaerobicTransExp), ctx.anaerobic);
+  initializeOneModelParam(modelParams, "soilMethaneRate", &(params.soilMethaneRate), ctx.anaerobic);
+  initializeOneModelParam(modelParams, "litterMethaneRate", &(params.litterMethaneRate), ctx.anaerobic);
 
   // NOLINTEND
   // clang-format on
@@ -456,7 +459,7 @@ void outputState(FILE *out, int year, int day, double time) {
       trackers.rh, trackers.rtot, trackers.evapotranspiration);
   fprintf(out, "%19.4f %8.4f %9.4f %10.4f %9.6f %10.4f %8.4f",
           fluxes.transpiration, envi.minN, envi.soilOrgN, envi.litterN,
-          fluxes.nVolatilization, fluxes.nLeaching, 0.0);
+          fluxes.nVolatilization, fluxes.nLeaching, trackers.methane);
   fprintf(out, "%12.4f %9.5f %9.5f\n", envi.plantWoodCStorageDelta,
           balanceTracker.deltaC, balanceTracker.deltaN);
 }
@@ -1054,7 +1057,11 @@ void vegResp2(double *folResp, double *woodResp, double *growthResp,
 /////////////////
 
 // ensure that all the allocation to wood + leaves + fine roots < 1,
-// and calculate coarse root allocation
+// calculate coarse root allocation as:
+//   coarse = 1 - leaf - wood - fine, making sum(params) = 1
+// require:
+//   leaf, wood, fine root < 1 individually
+//   coarse root >=0 which enforces leaf + wood + fine root <= 1
 void ensureAllocation(void) {
   // :: from [3], root model description
   params.coarseRootAllocation = 1 - params.leafAllocation -
@@ -1122,9 +1129,8 @@ double calcRespMoistEffect(double water, double whc) {
       // Unimodal moisture response: suppressed under dry conditions, maximal
       // at intermediate moisture, reduced under saturated/anoxic conditions
 
-      double f_a = params.fAnoxia;
       // Aerobic water availability (dry limitation)
-      double D_aer = fmin(fmax(f_whc / f_a, 0), 1);
+      double D_aer = fmin(fmax(f_whc / params.fAnoxia, 0), 1);
       // Anaerobic index (oxygen limitation proxy)
       double A = calcAnaerobicIndex(water, whc);
       // Uni-modal moisture response
@@ -1318,7 +1324,7 @@ void calcRootAndWoodFluxes(void) {
 /*!
  * Calculate mineral N volatilization flux
  */
-void calcNVolatilizationFlux() {
+void calcNVolatilizationFlux(void) {
   // flux = k_vol * nMin * Dtemp * Dwater
   // Note k_vol is in units of day^-1, so we do not need to divide
   // by climate length to make this a flux
@@ -1333,7 +1339,7 @@ void calcNVolatilizationFlux() {
 /*!
  * Calculate mineral N leaching flux
  */
-void calcNLeachingFlux() {
+void calcNLeachingFlux(void) {
   double phi;
   // phi is (drainage / soilWHC) between 0 and 1
   if ((fluxes.drainage / params.soilWHC) < 1) {
@@ -1348,7 +1354,7 @@ void calcNLeachingFlux() {
 /**
  * Calculate nitrogen fluxes for soil and litter pools
  */
-void calcNPoolFluxes() {
+void calcNPoolFluxes(void) {
   // C:N ratios for litter and soil, needed in most of the succeeding calcs
   double litterCN = calcLitterCN();
   double soilCN = calcSoilCN();
@@ -1379,6 +1385,60 @@ void calcNPoolFluxes() {
   fluxes.nMin = litterMin + soilMin;
 }
 
+/**
+ * Calculate methane flux
+ */
+void calcMethaneFlux(void) {
+  // Like soil respiration, but with own moisture dep and no tillage or CN
+  double tempEffect = calcTempEffect(climate->tsoil);
+  double moistEffect = calcMethaneMoistEffect(envi.soilWater, params.soilWHC);
+
+  fluxes.soilMethane =
+      params.soilMethaneRate * envi.soilC * tempEffect * moistEffect;
+  if (ctx.litterPool) {
+    fluxes.litterMethane =
+        params.litterMethaneRate * envi.litterC * tempEffect * moistEffect;
+  } else {
+    fluxes.litterMethane = 0.0;
+  }
+}
+
+void resetFluxes(void) {
+  fluxes.photosynthesis = 0.0;
+  fluxes.leafLitter = 0.0;
+  fluxes.woodLitter = 0.0;
+  fluxes.rVeg = 0.0;
+  fluxes.rSoil = 0.0;
+  fluxes.rain = 0.0;
+  fluxes.transpiration = 0.0;
+  fluxes.drainage = 0.0;
+  fluxes.litterToSoil = 0.0;
+  fluxes.rLitter = 0.0;
+  fluxes.snowFall = 0.0;
+  fluxes.snowMelt = 0.0;
+  fluxes.sublimation = 0.0;
+  fluxes.immedEvap = 0.0;
+  fluxes.fastFlow = 0.0;
+  fluxes.evaporation = 0.0;
+  fluxes.fineRootLoss = 0.0;
+  fluxes.coarseRootLoss = 0.0;
+  fluxes.fineRootCreation = 0.0;
+  fluxes.coarseRootCreation = 0.0;
+  fluxes.rCoarseRoot = 0.0;
+  fluxes.rFineRoot = 0.0;
+  fluxes.leafCreation = 0.0;
+  fluxes.leafOnCreation = 0.0;
+  fluxes.woodCreation = 0.0;
+  fluxes.nVolatilization = 0.0;
+  fluxes.nLeaching = 0.0;
+  fluxes.nOrgSoil = 0.0;
+  fluxes.nOrgLitter = 0.0;
+  fluxes.nMin = 0.0;
+  fluxes.soilMethane = 0.0;
+  fluxes.litterMethane = 0.0;
+  // event fluxes are handled in events.c:resetEventFluxes()
+}
+
 /*!
  * Calculate flux terms for sipnet as part of main model flow
  *
@@ -1401,6 +1461,9 @@ void calculateFluxes(void) {
   // net rain, equal to (rain - immedEvap) (cm/day)
   double netRain;
 
+  // Let's make sure to get all fluxes to zero before starting this
+  resetFluxes();
+
   // Psn, moisture and water fluxes
   lai = envi.plantLeafC / params.leafCSpWt;  // current lai
 
@@ -1438,6 +1501,11 @@ void calculateFluxes(void) {
   // Soil respiration
   calcSoilRespiration(climate->tsoil, envi.soilWater, params.soilWHC);
 
+  // Methane
+  if (ctx.anaerobic) {
+    calcMethaneFlux();
+  }
+
   // Nitrogen cycle
   //
   // Many of the nitrogen fluxes depend on carbon flux calculations, so make
@@ -1580,6 +1648,9 @@ void updateTrackers(double oldSoilWater) {
   trackers.woodCreation = fluxes.woodCreation * climate->length;
   trackers.n2o = fluxes.nVolatilization * climate->length;
 
+  trackers.methane =
+      (fluxes.soilMethane + fluxes.litterMethane) * climate->length;
+
   // evapotranspiration includes water lost to evaporation from canopy
   // irrigation (fluxes.eventEvap)
   trackers.evapotranspiration =
@@ -1673,13 +1744,14 @@ void updateMainPools() {
 void updatePoolsForSoil(void) {
   if (ctx.litterPool) {
     // :: from [2], litter model description
-    envi.litterC += (fluxes.woodLitter + fluxes.leafLitter -
-                     fluxes.litterToSoil - fluxes.rLitter) *
-                    climate->length;
+    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.litterToSoil - fluxes.rSoil - fluxes.soilMethane) *
                   climate->length;
   } else {
     // Normal pool (single pool, no microbes)
@@ -1688,9 +1760,10 @@ void updatePoolsForSoil(void) {
     //     L_l = fluxes.leafLitter
     //     R_h = fluxes.rSoil
     // :: from [3], root terms
-    envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss +
-                   fluxes.woodLitter + fluxes.leafLitter - fluxes.rSoil) *
-                  climate->length;
+    envi.soilC +=
+        (fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.woodLitter +
+         fluxes.leafLitter - fluxes.rSoil - fluxes.soilMethane) *
+        climate->length;
   }
 
   // :: from [3], root model description, except that we deal with

src/sipnet/state.h

@@ -363,6 +363,16 @@ typedef struct Parameters {
   // Methane anoxia sensitivity, controls sharpness of anaerobic transition
   // dimensionless, greater than 1
   double anaerobicTransExp;
+
+  // ******
+  // Methane production
+  // ******
+
+  // Relative methane production rate in the soil pool, in [0, 1), per day
+  double soilMethaneRate;
+
+  // Relative methane production rate in the litter pool, in [0, 1), per day
+  double litterMethaneRate;
 } Params;
 
 #define NUM_PARAMS (sizeof(Params) / sizeof(double))
@@ -554,6 +564,14 @@ typedef struct FluxVars {
   double eventInputN;
   // Total system nitrogen output, for mass balance checks
   double eventOutputN;
+
+  // ****************************************
+  // Methane production
+
+  // Methane produced from soil
+  double soilMethane;
+  // Methane produced from litter
+  double litterMethane;
 } Fluxes;
 
 // Global var
@@ -617,6 +635,9 @@ typedef struct TrackerVars {  // variables to track various things
   // g C * m^-2 wood creation
   double woodCreation;
 
+  // g C * m^-2 methane production
+  double methane;
+
   // g N * m^-2 ground area, Mineral N lost to volatilization
   double n2o;
 

tests/sipnet/test_events_types/testEventFertilization.c

@@ -70,6 +70,7 @@ int run(void) {
 
   //// TWO HARVEST EVENTS
   updateIntContext("litterPool", 1, CTX_TEST);
+  updateIntContext("anaerobic", 1, CTX_TEST);
   updateIntContext("nitrogenCycle", 1, CTX_TEST);
   logTest("Litter pool is %s\n", ctx.litterPool ? "on" : "off");
   logTest("Nitrogen cycle is %s\n", ctx.nitrogenCycle ? "on" : "off");

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
+TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c
 
 # The rest is boilerplate, likely copyable as is to a new test directory
 TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o)
@@ -40,6 +40,6 @@ $(RUN_EXECUTABLES):
 	./$(basename $@)
 
 clean:
-	rm -f $(TEST_OBJ_FILES) $(TEST_EXECUTABLES) events.in events.out
+	rm -f $(TEST_OBJ_FILES) $(TEST_EXECUTABLES) events.in events.out balance.out
 
 .PHONY: all tests clean run $(RUN_EXECUTABLES)