SIP180 Add N gain from fixation

src/sipnet/sipnet.c

@@ -1,3 +1,3 @@
 /* SiPnET: Simplified PnET
 
    Author: Bill Sacks  [email protected]
@@ -390,13 +390,14 @@
     initializeOneModelParam(modelParams, "mineralNInit", &(params.minNInit), ctx.nitrogenCycle);
     initializeOneModelParam(modelParams, "nVolatilizationFrac", &(params.nVolatilizationFrac), ctx.nitrogenCycle);
     initializeOneModelParam(modelParams, "nLeachingFrac", &(params.nLeachingFrac), ctx.nitrogenCycle);
+    initializeOneModelParam(modelParams, "nFixationFrac", &(params.nFixationFrac), ctx.nitrogenCycle);
 
   // NOLINTEND
-  // clang-format on
+    // clang-format on
 
-  readModelParams(modelParams, paramF);
+    readModelParams(modelParams, paramF);
 
-  fclose(paramF);
+    fclose(paramF);
 }
 
 /*!
@@ -410,9 +411,9 @@
   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 n2oFlux nLeachFlux\n");
+  fprintf(out, "npp nee cumNEE gpp rAboveground rSoil rRoot ra rh rtot "
+               "evapotranspiration fluxestranspiration minN n2oFlux nLeachFlux "
+               "nFixationFlux\n");
 }
 
 /*!
@@ -433,12 +434,12 @@
           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.6f %8.4f\n",
+          "%8.4f %8.3f %8.6f %8.4f %8.6f\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, fluxes.nVolatilization,
-          fluxes.nLeaching);
+          fluxes.nLeaching, fluxes.nFixation);
 }
 
 // de-allocate space used for climate linked list
@@ -1241,6 +1242,16 @@
   fluxes.nLeaching = envi.minN * phi * params.nLeachingFrac;
 }
 
+/*!
+ * Calculate mineral N fixation flux
+ */
+void calcNFixationFlux() {
+  // flux = k_fix * NPP, g N * m^-2 * day^-1
+  // k_fix represents C cost to fix N in range 0.07-0.17 gN/gC
+
+  fluxes.nFixation = params.nFixationFrac * getMeanTrackerMean(meanNPP);
+}
+
 /*!
  * Calculate flux terms for sipnet as part of main model flow
  *
@@ -1310,6 +1321,7 @@
     // Leaching depends on drainage flux so makes sure calcNLeachingFlux
     // occurs after calcSoilWaterFluxes
     calcNLeachingFlux();
+    calcNFixationFlux();
   }
 }
 
@@ -1598,7 +1610,8 @@
       climate->length;
 
   // Nitrogen Cycle
-  envi.minN -= (fluxes.nVolatilization + fluxes.nLeaching) * climate->length;
+  envi.minN += (fluxes.nFixation - fluxes.nVolatilization - fluxes.nLeaching) *
+               climate->length;
 }
 
 // !!! main runner function !!!

src/sipnet/state.h

@@ -367,6 +367,9 @@ typedef struct Parameters {
   // Fraction of mineral N lost to leaching per day
   double nLeachingFrac;
 
+  // Fraction of NPP available to be fixed by vegetation
+  double nFixationFrac;
+
 } Params;
 
 #define NUM_PARAMS (sizeof(Params) / sizeof(double))
@@ -535,6 +538,8 @@ typedef struct FluxVars {
   double nVolatilization;
   // Mineral N lost to leaching
   double nLeaching;
+  // Mineral N gained through plant fixation
+  double nFixation;
 
   // ****************************************
   // Fluxes for event handling

tests/sipnet/test_modeling/testNitrogenCycle.c

@@ -17,7 +17,7 @@ void setupTests(void) {
   ctx.litterPool = 1;
   ctx.nitrogenCycle = 1;
 }
-
+/*
 void initState(double initN, double nVol, double nLeachFrac) {
   // per test
   envi.minN = initN;
@@ -40,7 +40,83 @@ void initState(double initN, double nVol, double nLeachFrac) {
   params.baseSoilResp = 0.06;
   params.soilRespQ10 = 2.9;
 }
+*/
+void initNLeachingState(double initN, double nLeachFrac) {
+  // per test
+  envi.minN = initN;
+  params.nLeachingFrac = nLeachFrac;
+
+  // static
+  envi.soilWater = 5.0;
+  envi.soil = 1.5;
+  envi.litter = 1;
+
+  params.minNInit = 0.0;
+  params.soilWHC = 10.0;
+
+  fluxes.nVolatilization = 0;
+  fluxes.nLeaching = 0;
+
+  // Values from russell_2 smoke test
+  params.soilRespMoistEffect = 1.0;
+  params.baseSoilResp = 0.06;
+  params.soilRespQ10 = 2.9;
+}
+
+void initNVolatilizationState(double initN, double nVol) {
+  // per test
+  envi.minN = initN;
+  params.nVolatilizationFrac = nVol;
+
+  // static
+  envi.soilWater = 5.0;
+  envi.soil = 1.5;
+  envi.litter = 1;
+
+  params.minNInit = 0.0;
+  params.soilWHC = 10.0;
+
+  fluxes.nVolatilization = 0;
+
+  // Values from russell_2 smoke test
+  params.soilRespMoistEffect = 1.0;
+  params.baseSoilResp = 0.06;
+  params.soilRespQ10 = 2.9;
+}
+
+void initNFixationState(double initN, double nFixFrac) {
+  // per test
+  envi.minN = initN;
+  params.nFixationFrac = nFixFrac;
+
+  // static
+  envi.soilWater = 5.0;
+  envi.soil = 1.5;
+  envi.litter = 1;
+
+  params.minNInit = 0.0;
+  params.soilWHC = 10.0;
+
+  fluxes.nVolatilization = 0;
+  fluxes.nLeaching = 0;
+  fluxes.nFixation = 0;
+
+  // Values from russell_2 smoke test
+  params.soilRespMoistEffect = 1.0;
+  params.baseSoilResp = 0.06;
+  params.soilRespQ10 = 2.9;
+}
+
+int checkFlux(double calcFlux, double expFlux) {
+  int status = 0;
+  if (!compareDoubles(calcFlux, expFlux)) {
+    status = 1;
+    logTest("Calculated N flux is %f, expected %f\n", calcFlux, expFlux);
+  }
 
+  return status;
+}
+/*
 int checkVolatilizationFlux(double expNVolFlux) {
   int status = 0;
   if (!compareDoubles(fluxes.nVolatilization, expNVolFlux)) {
@@ -62,6 +138,41 @@ int checkNLeachingFlux(double expNLeachingFlux) {
 
   return status;
 }
+*/
+
+int testNFixation(void) {
+  int status = 0;
+  double minN;
+  double nFixFrac;
+  double expNFixation;
+
+  minN = 1;
+  nFixFrac = 0.25;
+  initNFixationState(minN, nFixFrac);
+  expNFixation = nFixFrac;
+  calcNFixationFlux();
+  status |= checkFlux(fluxes.nFixation, expNFixation);
+
+  minN = 1;
+  nFixFrac = 0.5;
+  initNFixationState(minN, nFixFrac);
+  expNFixation = nFixFrac;
+  calcNFixationFlux();
+  status |= checkFlux(fluxes.nFixation, expNFixation);
+
+  // Check minN for the last
+  updatePoolsForSoil();
+  double expMinN = minN + (expNFixation * climate->length);
+  int minStatus = 0;
+  if (!compareDoubles(envi.minN, expMinN)) {
+    minStatus = 1;
+    logTest("Fixation minN pool is %8.3f, expected %8.3f\n", envi.minN,
+            expMinN);
+  }
+  status |= minStatus;
+
+  return status;
+}
 
 int testNLeaching(void) {
   int status = 0;
@@ -73,36 +184,39 @@ int testNLeaching(void) {
   minN = 1;
   nLeachFrac = 0.5;
   fluxes.drainage = 5;
-  initState(minN, 0, nLeachFrac);
+  initNLeachingState(minN, nLeachFrac);
   if ((fluxes.drainage / params.soilWHC) < 1) {
     phi = fluxes.drainage / params.soilWHC;
   } else {
     phi = 1;
   }
   expNLeaching = minN * phi * nLeachFrac;
   calcNLeachingFlux();
-  status |= checkNLeachingFlux(expNLeaching);
+  // status |= checkNLeachingFlux(expNLeaching);
+  status |= checkFlux(fluxes.nLeaching, expNLeaching);
 
   minN = 1;
   nLeachFrac = 0.5;
   fluxes.drainage = 20;
-  initState(minN, 0, nLeachFrac);
+  initNLeachingState(minN, nLeachFrac);
   if ((fluxes.drainage / params.soilWHC) < 1) {
     phi = fluxes.drainage / params.soilWHC;
   } else {
     phi = 1;
   }
   expNLeaching = minN * phi * nLeachFrac;
   calcNLeachingFlux();
-  status |= checkNLeachingFlux(expNLeaching);
+  // status |= checkNLeachingFlux(expNLeaching);
+  status |= checkFlux(fluxes.nLeaching, expNLeaching);
 
   // Check minN for the last
   updatePoolsForSoil();
   double expMinN = minN - (expNLeaching * 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);
+    logTest("Leaching minN pool is %8.3f, expected %8.3f\n", envi.minN,
+            expMinN);
   }
   status |= minStatus;
 
@@ -117,27 +231,30 @@ int testNVolatilization(void) {
 
   minN = 2;
   nVolFrac = 0.1;
-  initState(minN, nVolFrac, 0);
+  initNVolatilizationState(minN, nVolFrac);
   double tEffect = calcTempEffect(climate->tsoil);
   double mEffect = calcMoistEffect(envi.soilWater, params.soilWHC);
   expNVolFlux = nVolFrac * minN * tEffect * mEffect;
   calcNVolatilizationFlux();
-  status |= checkVolatilizationFlux(expNVolFlux);
+  // status |= checkVolatilizationFlux(expNVolFlux);
+  status |= checkFlux(fluxes.nVolatilization, expNVolFlux);
 
   // easy proportionality test - doubling params should double output
   minN *= 2;
   expNVolFlux *= 2;
-  initState(minN, nVolFrac, 0);
+  initNVolatilizationState(minN, nVolFrac);
   calcNVolatilizationFlux();
-  status |= checkVolatilizationFlux(expNVolFlux);
+  // status |= checkVolatilizationFlux(expNVolFlux);
+  status |= checkFlux(fluxes.nVolatilization, expNVolFlux);
 
   // Check minN for the last
   updatePoolsForSoil();
   double expMinN = minN - expNVolFlux * 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);
+    logTest("Volatilization minN pool is %8.3f, expected %8.3f\n", envi.minN,
+            expMinN);
   }
   status |= minStatus;
 
@@ -151,7 +268,7 @@ int testFertilization(void) {
   double expMinN, expEventMinNFlux, expNVolFlux;
 
   // init minN 2, nVol 0.1
-  initState(initN, nVolFrac, 0);
+  initNVolatilizationState(initN, nVolFrac);
 
   // fert event: 15 5 10
   double fertMinN = 10;
@@ -201,6 +318,8 @@ int run(void) {
 
   status |= testNLeaching();
 
+  status |= testNFixation();
+
   return status;
 }
 

tests/smoke/russell_2/sipnet.param

@@ -80,3 +80,4 @@ microbePulseEff 0.45
 mineralNInit 1.0
 nVolatilizationFrac 0.1
 nLeachingFrac 0.25
+nFixationFrac 0.11