SIP182b Add litter and soil org n fluxes

docs/CHANGELOG.md

@@ -34,6 +34,10 @@ sections to include in release notes:
 - Soil mineral pool (#170)
 - Nitrogen effects of fertilization (#173)
 - `logAppend` logging function (#173)
+- Nitrogen volatilization (#175)
+- Nitrogen leaching (#178)
+- Soil and litter organic N pools (#195, #198)
+- Organic N handling for fertilization and soil dynamics (#198)
 
 ### Fixed
 

docs/model-structure.md

@@ -236,7 +236,7 @@ $$
 R_{H} = f_{R_H} \cdot 
   \left(\sum_j  K_j \cdot C_j 
     \right) \cdot 
-    D_{\text{temp}} \cdot D_{\text{water,}R_H} \cdot D_{CN} \mathfrak{\cdot D_{\text{tillage}}} 
+    D_{\text{temp}} \cdot D_{\text{water,}R_H} \cdot D_{CN} \cdot D_{\text{tillage}}
     \tag{7}\label{eq:rh}
 $$
 
@@ -261,7 +261,7 @@ The calculation of methane flux  $(F^C_{CH_4})$ is analagous to to that of $R_H$
 
 The carbon and nitrogen cycle are tightly coupled by the C:N ratios of plant and organic matter pools. The C:N ratio of plant biomass pools is fixed, while the C:N ratio of soil organic matter and litter pools is dynamic.
 
-### $\frak{Fixed \ Plant \ C:N \ Ratios}$
+### Fixed Plant C:N Ratios
 
 Plant biomass pools have a fixed CN ratio and are thus stoichiometrically coupled to carbon:
 
@@ -275,15 +275,15 @@ Where $i$ is the leaf, wood, fine root, or coarse root pool. This relationship a
 
 Soil organic matter and litter pools have dynamic CN that is determined below.
 
-### $\frak{Dynamic \ Soil  \ Organic \ Matter \ and \ Litter \ C:N \ Ratios}$
+### Dynamic Soil Organic Matter and Litter C:N Ratios
 
-The change in the soil C:N ratio over time of soil and litter pools depends on the rate of change of carbon and nitrogen in the pool, normalized by the total nitrogen in the pool. This makes sense as it captures how changes in carbon and nitrogen affect their ratio.
+In SIPNET, the C:N ratio of soil and litter pools is calculated directly from the carbon and nitrogen pools.
 
 $$
-\frac{dCN_{\text{j}}}{dt} = \frac{1}{N_{\text{j}}} \left( \frac{dC_{\text{j}}}{dt} - CN_{\text{j}} \cdot \frac{dN_{\text{j}}}{dt} \right) \tag{10}\label{eq:cn}
+CN_j = \frac{C_j}{N_j}, \qquad j \in \{\text{soil, litter}\}.
 $$
 
-$$\small j \in \{\text{soil, litter}\}$$
+This is used to calculate C:N-dependency $D_{CN}$ used in Eq. \eqref{eq:cn_dep}.
 
 ### $\frak{C:N \ Dependency \ Function \ (D_{CN})}$
 
@@ -308,7 +308,7 @@ $$
 $$\small i \in \{\text{leaf, wood, fine root, coarse root}\}$$
 
 
-### $\frak{Litter \ Nitrogen}$
+### Litter Nitrogen $N_\text{litter}$
 
 The change in litter nitrogen over time, $N_\text{litter}$ is determined by inputs including leaf and wood litter, nitrogen in organic matter amendments, and losses to mineralization:
 
@@ -320,21 +320,26 @@ $$
   F^N_\text{litter,min} \tag{13}\label{eq:litter_dndt}
 $$
 
-$$\small i \in \{\text{leaf, wood, fine root, coarse root}\}$$
+$$\small i \in \{\text{leaf, wood}\}$$
 
 The flux of nitrogen from living biomass to the litter pool is proportional to the carbon content of the biomass, based on the C:N ratio of the biomass pool \eqref{eq:cn_stoich}. Similarly, nitrogen from organic matter amendments is calculated from the carbon content and the C:N ratio of the inputs.
 
-### $\frak{Soil \ Organic \ Nitrogen}$
+### Soil Organic Nitrogen $N_\text{org,soil}$
+
+The change in soil nitrogen over time, $N_\text{org,soil}$ is determined by inputs including root loss, litter decomposition, and losses to mineralization:
 
 $$
-  \frac{dN_\text{org,soil}}{dt} = 
+  \frac{dN_\text{org,soil}}{dt} =
+  \sum_{i} F^N_i +
    F^N_\text{litter} -
    F^N_\text{soil,min} \tag{14}\label{eq:org_soil_dndt}
 $$
 
+$$\small i \in \{\text{fine root, coarse root}\}$$
+
 The change in nitrogen pools in this model is proportional to the ratio of carbon to nitrogen in the pool. Equations for the evolution of soil and litter CN are below.
 
-### $\frak{Soil \ Mineral \ Nitrogen \ F^N_\text{min}}$
+### Soil Mineral Nitrogen $N_\text{min}$
 
 Change in the mineral nitrogen pool over time is determined by inputs from mineralization and fertilization, and losses to volatilization, leaching, and plant uptake:
 
@@ -351,7 +356,7 @@ $$
 
 Mineralization and fertilization add to the mineral nitrogen pool. Losses include volatilization, leaching, and plant uptake, described below. Fixed N enters the plant pool directly (Eq. \eqref{eq:n_fix_demand}).
 
-### $\frak{N \ Mineralization \ (F^N_\text{min})}$
+### Nitrogen Mineralization $F^N_\text{min}$
 
 
 Total nitrogen mineralization is proportional to the total heterotrophic respiration from soil and litter pools, divided by the C:N ratio of the pool. The effects of temperature, moisture, tillage, and C:N ratio on mineralization rate are captured in the calculation of $R_\text{H}$.
@@ -376,7 +381,7 @@ $$
 F^N_\mathrm{vol} = K_\text{vol} \cdot N_\text{min} \cdot D_{\text{temp}} \cdot D_{\text{water}R_H} \tag{17}\label{eq:n_vol}
 $$
 
-### $\frak{Nitrogen \ Leaching \ F^N_\text{leach}}$
+### Nitrogen Leaching $F^N_\text{leach}$
 
 $$
 F^N_\text{leach} = N_\text{min} \cdot F^W_{drainage} \cdot f_{N leach} \tag{18}\label{eq:n_leach}

docs/user-guide/model-inputs.md

@@ -2,7 +2,7 @@
 
 These are the input files needed to run SIPNET:
 
-1. `sipnet.param`: Model parameter file.
+1. `sipnet.param`: Model parameter file. See [Parameters](../parameters.md) for full description of this file.
 2. `<sitename>.clim`: Climate file, provides weather data for each time step of the simulation period.
 3. [optional] `events.in`: Agronomic events.
 4. [optional] Run time options (command line arguments or config file `sipnet.in`)

src/common/context.c

@@ -199,6 +199,11 @@ void validateContext(void) {
     hasError = 1;
   }
 
+  if (ctx.nitrogenCycle && !ctx.litterPool) {
+    logError("nitrogen-cycle require litter-pool to be turned on\n");
+    hasError = 1;
+  }
+
   if (hasError) {
     exit(EXIT_CODE_BAD_PARAMETER_VALUE);
   }

src/sipnet/events.c

@@ -36,6 +36,7 @@ EventNode *createEventNode(int year, int day, int eventType,
   newEvent->year = year;
   newEvent->day = day;
   newEvent->type = eventType;
+  static int nitrogenWarned = 0;
 
   switch (eventType) {
     case HARVEST: {
@@ -97,6 +98,12 @@ EventNode *createEventNode(int year, int day, int eventType,
       fParams->minN = minN;
       // params->nh4_no3_frac = nh4_nos_frac;
       newEvent->eventParams = fParams;
+
+      if (!ctx.nitrogenCycle && (orgN > 0.0 || minN > 0.0) && !nitrogenWarned) {
+        logWarning("Fertilization nitrogen quantities are being ignored since "
+                   "nitrogen cycle modeling is off\n");
+        nitrogenWarned = 1;
+      }
     } break;
     case PLANTING: {
       double leafC, woodC, fineRootC, coarseRootC;
@@ -440,9 +447,16 @@ void processEvents(void) {
         const double orgC = fertParams->orgC;
         const double minN = fertParams->minN;
 
-        fluxes.eventOrgN += orgN / climLen;
         fluxes.eventLitterC += orgC / climLen;
-        fluxes.eventMinN += minN / climLen;
+        if (ctx.nitrogenCycle) {
+          fluxes.eventOrgN += orgN / climLen;
+          fluxes.eventMinN += minN / climLen;
+        } else {
+          // As the warning says in readEventData(), we ignore N when the
+          // nitrogen cycle model is off
+          fluxes.eventOrgN += 0;
+          fluxes.eventMinN += 0;
+        }
 
         writeEventOut(gEvent, 3, "fluxes.eventOrgN", orgN / climLen,
                       "fluxes.eventLitterC", orgC / climLen, "fluxes.eventMinN",
@@ -458,26 +472,34 @@ void processEvents(void) {
 }
 
 void updatePoolsForEvents(void) {
+  // CARBON
   // Harvest and planting events
   envi.plantWoodC += fluxes.eventWoodC * climate->length;
   envi.plantLeafC += fluxes.eventLeafC * climate->length;
 
-  // Irrigation events
-  envi.soilWater += fluxes.eventSoilWater * climate->length;
-
   // Harvest and fertilization events
   if (ctx.litterPool) {
-    envi.litterOrgN += fluxes.eventOrgN * climate->length;
     envi.litter += fluxes.eventLitterC * climate->length;
   } else {
-    envi.soilOrgN += fluxes.eventOrgN * climate->length;
     envi.soil += fluxes.eventLitterC * climate->length;
   }
-  envi.minN += fluxes.eventMinN * climate->length;
 
   // Harvest and planting events
   envi.coarseRootC += fluxes.eventCoarseRootC * climate->length;
   envi.fineRootC += fluxes.eventFineRootC * climate->length;
+
+  // WATER
+  // Irrigation events
+  envi.soilWater += fluxes.eventSoilWater * climate->length;
+
+  // NITROGEN
+  // Fertilization events
+  // (Harvest and planting events TBD)
+  // Note: litter_pool is required for nitrogen_cycle
+  if (ctx.nitrogenCycle) {
+    envi.minN += fluxes.eventMinN * climate->length;
+    envi.litterOrgN += fluxes.eventOrgN * climate->length;
+  }
 }
 
 void freeEventList(void) {

src/sipnet/sipnet.c

@@ -392,6 +392,9 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) {
   initializeOneModelParam(modelParams, "litterOrgNInit", &(params.litterOrgNInit), ctx.nitrogenCycle);
   initializeOneModelParam(modelParams, "nVolatilizationFrac", &(params.nVolatilizationFrac), ctx.nitrogenCycle);
   initializeOneModelParam(modelParams, "nLeachingFrac", &(params.nLeachingFrac), ctx.nitrogenCycle);
+  initializeOneModelParam(modelParams, "leafCNRatio", &(params.leafCNRatio), ctx.nitrogenCycle);
+  initializeOneModelParam(modelParams, "woodCNRatio", &(params.woodCNRatio), ctx.nitrogenCycle);
+  initializeOneModelParam(modelParams, "rootCNRatio", &(params.rootCNRatio), ctx.nitrogenCycle);
 
   // NOLINTEND
   // clang-format on
@@ -409,12 +412,13 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) {
 void outputHeader(FILE *out) {
   fprintf(out, "Notes: (PlantWoodC, PlantLeafC, Soil and Litter in g C/m^2; "
                "Water and Snow in cm; SoilWetness is fraction of WHC;\n");
-  fprintf(out, "year day time plantWoodC plantLeafC woodCreation ");
-  fprintf(out, "soil microbeC coarseRootC fineRootC ");
-  fprintf(out, "litter soilWater soilWetnessFrac snow ");
-  fprintf(out, "npp nee cumNEE gpp rAboveground rSoil rRoot ra rh rtot "
-               "evapotranspiration fluxestranspiration minN soilOrgN "
-               "litterOrgN n2oFlux nLeachFlux\n");
+  fprintf(out, "year day  time plantWoodC plantLeafC woodCreation     ");
+  fprintf(out, "soil microbeC coarseRootC fineRootC   ");
+  fprintf(out, "litter soilWater soilWetnessFrac     snow      ");
+  fprintf(out, "npp      nee   cumNEE      gpp rAboveground    rSoil    "
+               "rRoot       ra       rh     rtot evapotranspiration ");
+  fprintf(out, "fluxestranspiration     minN  soilOrgN litterOrgN   n2oFlux "
+               "nLeachFlux\n");
 }
 
 /*!
@@ -426,21 +430,22 @@ void outputHeader(FILE *out) {
  */
 void outputState(FILE *out, int year, int day, double time) {
 
-  fprintf(out, "%4d %3d %5.2f %8.2f %8.2f %8.2f ", year, day, time,
+  fprintf(out, "%4d %3d %5.2f %10.2f %10.2f %12.2f ", year, day, time,
           envi.plantWoodC, envi.plantLeafC, trackers.woodCreation);
   fprintf(out, "%8.2f ", envi.soil);
   fprintf(out, "%8.2f ", envi.microbeC);
-  fprintf(out, "%8.2f %8.2f", envi.coarseRootC, envi.fineRootC);
-  fprintf(out, " %8.2f %8.2f %8.3f %8.2f ", envi.litter, envi.soilWater,
+  fprintf(out, "%11.2f %9.2f", envi.coarseRootC, envi.fineRootC);
+  fprintf(out, " %8.2f %9.2f %15.3f %8.2f ", envi.litter, envi.soilWater,
           trackers.soilWetnessFrac, envi.snow);
-  fprintf(out,
-          "%8.2f %8.2f %8.2f %8.2f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.8f "
-          "%8.4f %8.3f %8.4f %8.4f %8.6f %8.4f\n",
-          trackers.npp, trackers.nee, trackers.totNee, trackers.gpp,
-          trackers.rAboveground, trackers.rSoil, trackers.rRoot, trackers.ra,
-          trackers.rh, trackers.rtot, trackers.evapotranspiration,
-          fluxes.transpiration, envi.minN, envi.soilOrgN, envi.litterOrgN,
-          fluxes.nVolatilization, fluxes.nLeaching);
+  fprintf(
+      out,
+      "%8.2f %8.2f %8.2f %8.2f %12.3f %8.3f %8.3f %8.3f %8.3f %8.3f %18.8f ",
+      trackers.npp, trackers.nee, trackers.totNee, trackers.gpp,
+      trackers.rAboveground, trackers.rSoil, trackers.rRoot, trackers.ra,
+      trackers.rh, trackers.rtot, trackers.evapotranspiration);
+  fprintf(out, "%19.4f %8.3f %9.4f %10.4f %9.6f %10.4f\n", fluxes.transpiration,
+          envi.minN, envi.soilOrgN, envi.litterOrgN, fluxes.nVolatilization,
+          fluxes.nLeaching);
 }
 
 // de-allocate space used for climate linked list
@@ -1243,6 +1248,33 @@ void calcNLeachingFlux() {
   fluxes.nLeaching = envi.minN * phi * params.nLeachingFrac;
 }
 
+/**
+ * Calculate organic nitrogen fluxes
+ */
+void calcOrgNFluxes() {
+  double litterCN, soilCN;
+  // for both litter and soil, mineralization is calculated as heterotrophic
+  // respiration divided by the C:N ratio of that pool.
+
+  // litter
+  // The litter org N flux is determined by the carbon fluxes from wood and leaf
+  // litter, and N loss due to mineralization. N added via fertilization
+  // is handled elsewhere.
+  litterCN = envi.litter / envi.litterOrgN;
+  fluxes.nOrgLitter = fluxes.leafLitter / params.leafCNRatio +
+                      fluxes.woodLitter / params.woodCNRatio -
+                      fluxes.rLitter / litterCN;
+
+  // 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: woodCNRatio is used for coarse roots)
+  soilCN = envi.soil / envi.soilOrgN;
+  fluxes.nOrgSoil = (fluxes.litterToSoil - fluxes.rSoil) / soilCN +
+                    fluxes.fineRootLoss / params.rootCNRatio +
+                    fluxes.coarseRootLoss / params.woodCNRatio;
+}
+
 /*!
  * Calculate flux terms for sipnet as part of main model flow
  *
@@ -1307,11 +1339,14 @@ void calculateFluxes(void) {
   }
 
   // Nitrogen cycle
+  //
+  // Many of the nitrogen fluxes depend on carbon flux calculations, so make
+  // sure this stays at the bottom of this function (or after the carbon calcs,
+  // at least).
   if (ctx.nitrogenCycle) {
     calcNVolatilizationFlux();
-    // Leaching depends on drainage flux so makes sure calcNLeachingFlux
-    // occurs after calcSoilWaterFluxes
     calcNLeachingFlux();
+    calcOrgNFluxes();
   }
 }
 
@@ -1603,11 +1638,22 @@ void updatePoolsForSoil(void) {
   envi.fineRootC +=
       (fluxes.fineRootCreation - fluxes.fineRootLoss - fluxes.rFineRoot) *
       climate->length;
+}
 
+void updateNitrogenPools(void) {
   // Nitrogen Cycle
+  // :: from [5], nitrogen cycle model
+  // TBD: add equation numbers once published
+
+  // Soil mineral N (note we have one mineral pool for soil+litter)
+  // Mineral N additions from fertilization are handled with the events
   envi.minN -= (fluxes.nVolatilization + fluxes.nLeaching) * climate->length;
-  // envi.soilOrgN += ...  TBD
-  // envi.litterOrgN += ...  TBD
+
+  // Soil organic N
+  envi.soilOrgN += fluxes.nOrgSoil * climate->length;
+
+  // Litter organic N
+  envi.litterOrgN += fluxes.nOrgLitter * climate->length;
 }
 
 // !!! main runner function !!!
@@ -1638,6 +1684,9 @@ void updateState(void) {
   // Update soil carbon pools
   updatePoolsForSoil();
 
+  // Update nitrogen cycle pools
+  updateNitrogenPools();
+
   // Update pools for fluxes from events
   updatePoolsForEvents();