SIP180 Adding N demand, fixation, and uptake

docs/CHANGELOG.md

@@ -45,7 +45,11 @@ 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)
+<<<<<<< HEAD
+- Adds Nitrogen demand, fixation flux, and uptake flux (#265)
+=======
 - Methane production option, controlled by the `anaerobic` cli option (#269)
+>>>>>>> master
 
 ### Fixed
 

docs/model-structure.md

@@ -593,20 +593,33 @@ drainage.
 
 ### Plant Nitrogen Demand  $F^{N}_{\text{demand}}$
 
-Plant N demand is the amount of N required to support plant growth. This is calculated as the sum of changes in plant N
-pools:
+Plant N demand is the amount of N required to support plant growth. This is calculated as the sum of carbon creation fluxes divided by their respective C:N ratios:
 
 \begin{equation}
-F^N_\text{demand}=\frac{dN_\text{plant}}{dt} = \sum_{i} \frac{dN_{\text{plant,}i}}{dt}
-\label{eq:plant_n_demand}
+F^N_{\text{demand,}leafOn} = \frac{F^C_{\text{creation,}leafOn}}{CN_{\text{leaf}}} -
+\frac{F^C_{\text{creation,}leafOn}}{CN_{\text{wood}}}
+\label{eq:leaf_on_n_demand}
+\end{equation}
+
+\begin{equation}
+F^N_{\text{demand,}creation} = \sum_{i} \frac{F^C_{\text{creation,}i}}{CN_{\text{i}}} 
+\label{eq:creation_n_demand}
 \end{equation}
 
 \begin{equation*}
-\small i \in \{\text{leaf, wood, fine root, coarse root}\}
+\small i \in \{\text{wood, leaf, fine root, coarse root}\}
 \end{equation*}
 
-Each term in the sum is calculated according to \eqref{eq:plant_n}. Total plant N demand $F^N_\text{demand}$ is then
-partitioned between fixation and soil N uptake using \eqref{eq:n_fix_demand} and \eqref{eq:n_uptake_demand}.
+\begin{equation}
+F^N_{\text{demand,}total} =
+F^N_{\text{demand,}leafOn} +
+F^N_{\text{demand,}creation}
+\label{eq:plant_n_demand}
+\end{equation}
+
+Each term in the sum is calculated according to \eqref{eq:plant_n}. Total plant N demand $F^N_{\text{demand,}total}$ is then partitioned between fixation and soil N uptake using \eqref{eq:n_fix_demand} and \eqref{eq:n_uptake_demand}.
+
+**TODO:** possibly include more context about leaf on events
 
 ### Nitrogen Fixation and Uptake $F^N_\text{fix}, F^N_\text{uptake}$
 

src/sipnet/balance.c

@@ -56,9 +56,7 @@ void updateBalanceTrackerPostUpdate(void) {
 
   // NITROGEN
   if (ctx.nitrogenCycle) {
-    balanceTracker.inputsN =
-        // TODO: fluxes.fixation +
-        fluxes.eventInputN;
+    balanceTracker.inputsN = fluxes.nFixation + fluxes.eventInputN;
     balanceTracker.outputsN =
         fluxes.nLeaching + fluxes.nVolatilization + fluxes.eventOutputN;
 
@@ -112,14 +110,14 @@ void checkBalance(void) {
         "Carbon balance check failed (delta=%8.4f, Y: %d D: %d T: %4.2f)\n",
         balanceTracker.deltaC, climate->year, climate->day, climate->time);
   }
-  // RE-ENABLE WHEN N IS FIXED
-  // if (fabs(balanceTracker.deltaN) > 0.0) {
-  //  err = 1;
-  //  logInternalError("Nitrogen balance check failed (delta=%8.4f)\n",
-  //                   balanceTracker.deltaN);
-  //}
+  if (fabs(balanceTracker.deltaN) > 0.0) {
+    // err = 1;
+    logInternalError(
+        "Nitrogen balance check failed (delta=%8.4f, Y: %d D: %d T: %4.2f)\n",
+        balanceTracker.deltaN, climate->year, climate->day, climate->time);
+  }
   if (err) {
     logInternalError("Exiting\n");
-    //  exit(EXIT_CODE_INTERNAL_ERROR);
+    exit(EXIT_CODE_INTERNAL_ERROR);
   }
 }

src/sipnet/sipnet.c

@@ -380,6 +380,8 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) {
   initializeOneModelParam(modelParams, "woodCN", &(params.woodCN), ctx.nitrogenCycle);
   initializeOneModelParam(modelParams, "fineRootCN", &(params.fineRootCN), ctx.nitrogenCycle);
   initializeOneModelParam(modelParams, "kCN", &(params.kCN), ctx.nitrogenCycle);
+  initializeOneModelParam(modelParams, "nFixationFracMax", &(params.nFixationFracMax), ctx.nitrogenCycle);
+  initializeOneModelParam(modelParams, "halfNFixationMax", &(params.halfNFixationMax), ctx.nitrogenCycle);
 
   // New moisture dependency params
   initializeOneModelParam(modelParams, "fAnoxia", &(params.fAnoxia), ctx.anaerobic || ctx.nitrogenCycle);
@@ -433,7 +435,8 @@ void outputHeader(FILE *out) {
   fprintf(out, "npp      nee   cumNEE      gpp rAboveground    rSoil    "
                "rRoot       ra       rh     rtot evapotranspiration ");
   fprintf(out, "fluxestranspiration     minN  soilOrgN    litterN       n2o "
-               "nLeachFlux      ch4  nppStorage  bcdeltaC  bcdeltaN\n");
+               "nLeaching  nFixation  nUptake      ch4  nppStorage  bcdeltaC  "
+               "bcdeltaN\n");
 }
 /*!
  * Print current state values to output file
@@ -457,9 +460,10 @@ void outputState(FILE *out, int year, int day, double time) {
       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.4f %9.4f %10.4f %9.6f %10.4f %8.4f",
+  fprintf(out, "%19.4f %8.4f %9.4f %10.4f %9.6f %9.4f %10.4f %8.4f %8.4f",
           fluxes.transpiration, envi.minN, envi.soilOrgN, envi.litterN,
-          fluxes.nVolatilization, fluxes.nLeaching, trackers.methane);
+          fluxes.nVolatilization, trackers.nLeaching, trackers.nFixation,
+          trackers.nUptake, trackers.methane);
   fprintf(out, "%12.4f %9.5f %9.5f\n", envi.plantWoodCStorageDelta,
           balanceTracker.deltaC, balanceTracker.deltaN);
 }
@@ -1351,6 +1355,47 @@ void calcNLeachingFlux(void) {
   fluxes.nLeaching = envi.minN * phi * params.nLeachingFrac;
 }
 
+/*!
+ * Calculate plant N demand
+ */
+double calcPlantNDemand() {
+  // leafOnCreation is a transfer from the wood pool so
+  // demand should only count the difference in the N
+  // between those two pools, hence substracting
+  // params.woodCN from params.leafCN for leafOnCreation
+  double leafOnDemand = fluxes.leafOnCreation / params.leafCN -
+                        fluxes.leafOnCreation / params.woodCN;
+  // calculate demand from all other creation terms
+  double creationDemand = fluxes.woodCreation / params.woodCN +
+                          fluxes.leafCreation / params.leafCN +
+                          fluxes.fineRootCreation / params.fineRootCN +
+                          fluxes.coarseRootCreation / params.woodCN;
+  // total demand is leafOnDemand plus creationDemand
+  double totalDemand = leafOnDemand + creationDemand;
+  return totalDemand;
+}
+
+/*!
+ * Calculate plant N fixation and uptake fluxes
+ */
+void calcNFixationAndUptakeFluxes() {
+  double nFixationInhibition;
+  double nFixationFrac;
+  // Calculate inhibition of N fixation by soil mineral N
+  // using down-regulation function with increasing soil min N
+  // dimensionless between 0 and 1
+  nFixationInhibition =
+      params.halfNFixationMax / (params.halfNFixationMax + envi.minN);
+  // Calculate fraction of plant N demand met by fixation
+  // dimensionless
+  nFixationFrac = params.nFixationFracMax * nFixationInhibition;
+  // Calculate N fixation flux
+  double nDemand = calcPlantNDemand();
+  fluxes.nFixation = nFixationFrac * nDemand;
+  // Calculate N uptake flux
+  fluxes.nUptake = (1 - nFixationFrac) * nDemand;
+}
+
 /**
  * Calculate nitrogen fluxes for soil and litter pools
  */
@@ -1514,6 +1559,7 @@ void calculateFluxes(void) {
   if (ctx.nitrogenCycle) {
     calcNVolatilizationFlux();
     calcNLeachingFlux();
+    calcNFixationAndUptakeFluxes();
     calcNPoolFluxes();
   }
 }
@@ -1550,6 +1596,9 @@ void initTrackers(void) {
   trackers.rAboveground = 0.0;
 
   trackers.yearlyLitter = 0.0;
+  trackers.nLeaching = 0.0;
+  trackers.nFixation = 0.0;
+  trackers.nUptake = 0.0;
 }
 
 // If var < minVal, then set var = 0
@@ -1662,6 +1711,11 @@ void updateTrackers(double oldSoilWater) {
       (oldSoilWater + envi.soilWater) / (2.0 * params.soilWHC);
 
   trackers.yearlyLitter += fluxes.leafLitter;
+
+  // N cycle trackers
+  trackers.nLeaching = fluxes.nLeaching * climate->length;
+  trackers.nFixation = fluxes.nFixation * climate->length;
+  trackers.nUptake = fluxes.nUptake * climate->length;
 }
 
 void updateMeanTrackers(void) {
@@ -1782,9 +1836,8 @@ void updateNitrogenPools(void) {
 
   // Soil mineral N (note we have one mineral pool for soil+litter)
   // Mineral N additions from fertilization are handled with the events
-  //
-  // TODO: add plant uptake flux once implemented
-  envi.minN += (fluxes.nMin - fluxes.nVolatilization - fluxes.nLeaching) *
+  envi.minN += (fluxes.nMin - fluxes.nVolatilization - fluxes.nLeaching -
+                fluxes.nUptake) *
                climate->length;
 
   // Soil organic N

src/sipnet/state.h

@@ -349,6 +349,14 @@ typedef struct Parameters {
   // C:N ratio at which D_CN is 1/2 for soil and litter
   double kCN;
 
+  // Maximum fraction of plant N demand that can be met by fixation
+  // under low soil N, dimensionless between 0 and 1
+  double nFixationFracMax;
+
+  // Amount of mineral N at which fixation is reduced by half
+  // g N * m^-2 ground area
+  double halfNFixationMax;
+
   // ******
   // New moisture dependency functions
   // ******
@@ -532,6 +540,14 @@ typedef struct FluxVars {
   // g N * m^-2 ground area * day^-1
   double nMin;
 
+  // Plant N gained from fixation
+  // g N * m^-2 ground area * day^-1
+  double nFixation;
+
+  // Mineral N drawn from soil pool to meet plant N demand
+  // g N * m^-2 ground area * day^-1
+  double nUptake;
+
   // ****************************************
   // Fluxes for event handling
   //  - event fluxes tracked as part of modeling from [4]
@@ -641,6 +657,15 @@ typedef struct TrackerVars {  // variables to track various things
   // g N * m^-2 ground area, Mineral N lost to volatilization
   double n2o;
 
+  // g N * m^-2 N leached from soil mineral N pool
+  double nLeaching;
+
+  // g N * m^-2 N fixed by plants
+  double nFixation;
+
+  // g N * m^-2 N taken up by plants from soil mineral N pool
+  double nUptake;
+
 } Trackers;
 
 // Global var

tests/sipnet/test_modeling/balance.param

@@ -64,6 +64,8 @@ leafCN 20.0
 woodCN 100.0
 fineRootCN 40.0
 kCN 80.0
+nFixationFracMax 0.5
+halfNFixationMax 1.0
 fAnoxia 0.6
 anaerobicDecompRate 0.01
 anaerobicTransExp 2.0

tests/sipnet/test_modeling/testNitrogenCycle.c

@@ -32,6 +32,11 @@ void resetState() {
   // Leaching
   params.nLeachingFrac = 0;
   fluxes.nLeaching = 0;
+  // Fixation/Uptake
+  params.nFixationFracMax = 0;
+  params.halfNFixationMax = 0;
+  fluxes.nFixation = 0;
+  fluxes.nUptake = 0;
 }
 
 void setupTests(void) {
@@ -217,6 +222,81 @@ int testNLeaching(void) {
   return status;
 }
 
+/////
+// N Fixation and Uptake
+void initNFixationState(double initN, double nFixFracMax, double nFixHalved) {
+  resetState();
+  envi.minN = initN;
+  params.nFixationFracMax = nFixFracMax;
+  params.halfNFixationMax = nFixHalved;
+}
+
+int testNFixation(void) {
+  int status = 0;
+  double minN;
+  double nFixFracMax;
+  double nFixHalved;
+  double expNFixation;
+  double expNUptake;
+  double nDemand;
+  logTest("Running testNFixation\n");
+
+  // Half demand met by fixation, half met by uptake
+  minN = 2;
+  nFixFracMax = 1;
+  nFixHalved = 2;
+  nDemand = 10;
+  initNFixationState(minN, nFixFracMax, nFixHalved);
+  double nFixInhib = nFixHalved / (nFixHalved + minN);
+  double nFixFrac = nFixFracMax * nFixInhib;
+  expNFixation = nFixFrac * nDemand;
+  expNUptake = (1 - nFixFrac) * nDemand;
+  calcNFixationAndUptakeFluxes();
+  status |= checkFlux(fluxes.nFixation, expNFixation, "N fixation");
+  status |= checkFlux(fluxes.nUptake, expNUptake, "N uptake");
+
+  // Zero demand met by fixation, all demand met by uptake
+  minN = 1;
+  nFixFracMax = 0.5;
+  nFixHalved = 0;
+  nDemand = 10;
+  initNFixationState(minN, nFixFracMax, nFixHalved);
+  nFixInhib = nFixHalved / (nFixHalved + minN);
+  nFixFrac = nFixFracMax * nFixInhib;
+  expNFixation = nFixFrac * nDemand;
+  expNUptake = (1 - nFixFrac) * nDemand;
+  calcNFixationAndUptakeFluxes();
+  status |= checkFlux(fluxes.nFixation, expNFixation, "N fixation");
+  status |= checkFlux(fluxes.nUptake, expNUptake, "N uptake");
+
+  // Zero demand by uptake due to zero minN pool, this test will
+  // need to change when limitation is implemented
+  minN = 0;
+  nFixFracMax = 0.5;
+  nFixHalved = 2;
+  nDemand = 10;
+  initNFixationState(minN, nFixFracMax, nFixHalved);
+  nFixInhib = nFixHalved / (nFixHalved + minN);
+  nFixFrac = nFixFracMax * nFixInhib;
+  expNFixation = nFixFrac * nDemand;
+  expNUptake = (1 - nFixFrac) * nDemand;
+  calcNFixationAndUptakeFluxes();
+  status |= checkFlux(fluxes.nFixation, expNFixation, "N fixation");
+  status |= checkFlux(fluxes.nUptake, expNUptake, "N uptake");
+
+  // Check minN for the last
+  updateNitrogenPools();
+  double expMinN = minN - (expNUptake * climate->length);
+  int minStatus = 0;
+  if (!compareDoubles(envi.minN, expMinN)) {
+    minStatus = 1;
+    logTest("minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN);
+  }
+  status |= minStatus;
+
+  return status;
+}
+
 /////
 // Organic N pools
 void initOrganicNState(double initLitterN, double initSoilN) {