SIP90 Add N loss to leaching

.gitignore

@@ -27,6 +27,7 @@ tests/sipnet/*/*
 .DS_Store
 .Rproj.user
 .Rhistory
+src/sipnet/.vscode/*
 
 # Temporary and backup files
 *.bak

Makefile

@@ -81,7 +81,7 @@ clean:
 	rm -f .doxygen.stamp .mkdocs.stamp
 
 # UNIT TESTS
-SIPNET_TEST_DIRS:=$(shell find tests/sipnet -type d -mindepth 1 -maxdepth 1)
+SIPNET_TEST_DIRS:=$(shell find tests/sipnet -mindepth 1 -maxdepth 1 -type d)
 SIPNET_TEST_DIRS_RUN:= $(addsuffix .run, $(SIPNET_TEST_DIRS))
 SIPNET_TEST_DIRS_CLEAN:= $(addsuffix .clean, $(SIPNET_TEST_DIRS))
 

src/sipnet/sipnet.c

@@ -387,8 +387,9 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) {
   initializeOneModelParam(modelParams, "microbePulseEff", &(params.microbePulseEff), ctx.microbes);
 
   // Nitrogen cycle params from [5] LeBauer et al. (unpublished)
-  initializeOneModelParam(modelParams, "mineralNInit", &(params.minNInit), ctx.nitrogenCycle);
-  initializeOneModelParam(modelParams, "nVolatilizationFrac", &(params.nVolatilizationFrac), ctx.nitrogenCycle);
+    initializeOneModelParam(modelParams, "mineralNInit", &(params.minNInit), ctx.nitrogenCycle);
+    initializeOneModelParam(modelParams, "nVolatilizationFrac", &(params.nVolatilizationFrac), ctx.nitrogenCycle);
+    initializeOneModelParam(modelParams, "nLeachingFrac", &(params.nLeachingFrac), ctx.nitrogenCycle);
 
   // NOLINTEND
   // clang-format on
@@ -409,8 +410,9 @@ void outputHeader(FILE *out) {
   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\n");
+  fprintf(out,
+          "npp nee cumNEE gpp rAboveground rSoil rRoot ra rh rtot "
+          "evapotranspiration fluxestranspiration minN n2oFlux nLeachFlux\n");
 }
 
 /*!
@@ -431,11 +433,12 @@ void outputState(FILE *out, int year, int day, double time) {
           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\n",
+          "%8.4f %8.3f %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, fluxes.nVolatilization);
+          fluxes.transpiration, envi.minN, fluxes.nVolatilization,
+          fluxes.nLeaching);
 }
 
 // de-allocate space used for climate linked list
@@ -1223,6 +1226,21 @@ void calcNVolatilizationFlux() {
       params.nVolatilizationFrac * envi.minN * d_temp * d_water;
 }
 
+/*!
+ * Calculate mineral N leaching flux
+ */
+void calcNLeachingFlux() {
+  double phi;
+  // phi is (drainage / soilWHC) between 0 and 1
+  if ((fluxes.drainage / params.soilWHC) < 1) {
+    phi = fluxes.drainage / params.soilWHC;
+  } else {
+    phi = 1;
+  }
+  // flux = nMin * phi * leaching fraction, g N * m^-2 * day^-1
+  fluxes.nLeaching = envi.minN * phi * params.nLeachingFrac;
+}
+
 /*!
  * Calculate flux terms for sipnet as part of main model flow
  *
@@ -1289,6 +1307,9 @@ void calculateFluxes(void) {
   // Nitrogen cycle
   if (ctx.nitrogenCycle) {
     calcNVolatilizationFlux();
+    // Leaching depends on drainage flux so makes sure calcNLeachingFlux
+    // occurs after calcSoilWaterFluxes
+    calcNLeachingFlux();
   }
 }
 
@@ -1577,7 +1598,7 @@ void updatePoolsForSoil(void) {
       climate->length;
 
   // Nitrogen Cycle
-  envi.minN -= fluxes.nVolatilization * climate->length;
+  envi.minN -= (fluxes.nVolatilization + fluxes.nLeaching) * climate->length;
 }
 
 // !!! main runner function !!!

src/sipnet/state.h

@@ -358,12 +358,15 @@ typedef struct Parameters {
   // ****************************************
   // Nitrogen Cycle
 
-  // Initial soil mineral nitrogen pool amount, g C * m^-2 ground area
+  // Initial soil mineral nitrogen pool amount, g N * m^-2 ground area
   double minNInit;
 
   // Fraction of mineral N available to be volatilized per day, d^-1
   double nVolatilizationFrac;
 
+  // Fraction of mineral N lost to leaching per day
+  double nLeachingFrac;
+
 } Params;
 
 #define NUM_PARAMS (sizeof(Params) / sizeof(double))
@@ -530,6 +533,8 @@ typedef struct FluxVars {
 
   // Mineral N lost to volatilization
   double nVolatilization;
+  // Mineral N lost to leaching
+  double nLeaching;
 
   // ****************************************
   // Fluxes for event handling

tests/sipnet/test_modeling/testNitrogenCycle.c

@@ -18,10 +18,11 @@ void setupTests(void) {
   ctx.nitrogenCycle = 1;
 }
 
-void initState(double initN, double nVol) {
+void initState(double initN, double nVol, double nLeachFrac) {
   // per test
   envi.minN = initN;
   params.nVolatilizationFrac = nVol;
+  params.nLeachingFrac = nLeachFrac;
 
   // static
   envi.soilWater = 5.0;
@@ -32,6 +33,7 @@ void initState(double initN, double nVol) {
   params.soilWHC = 10.0;
 
   fluxes.nVolatilization = 0;
+  fluxes.nLeaching = 0;
 
   // Values from russell_2 smoke test
   params.soilRespMoistEffect = 1.0;
@@ -50,6 +52,63 @@ int checkVolatilizationFlux(double expNVolFlux) {
   return status;
 }
 
+int checkNLeachingFlux(double expNLeachingFlux) {
+  int status = 0;
+  if (!compareDoubles(fluxes.nLeaching, expNLeachingFlux)) {
+    status = 1;
+    logTest("N leaching flux is %f, expected %f\n", fluxes.nLeaching,
+            expNLeachingFlux);
+  }
+
+  return status;
+}
+
+int testNLeaching(void) {
+  int status = 0;
+  double minN;
+  double nLeachFrac;
+  double expNLeaching;
+  double phi;
+
+  minN = 1;
+  nLeachFrac = 0.5;
+  fluxes.drainage = 5;
+  initState(minN, 0, nLeachFrac);
+  if ((fluxes.drainage / params.soilWHC) < 1) {
+    phi = fluxes.drainage / params.soilWHC;
+  } else {
+    phi = 1;
+  }
+  expNLeaching = minN * phi * nLeachFrac;
+  calcNLeachingFlux();
+  status |= checkNLeachingFlux(expNLeaching);
+
+  minN = 1;
+  nLeachFrac = 0.5;
+  fluxes.drainage = 20;
+  initState(minN, 0, nLeachFrac);
+  if ((fluxes.drainage / params.soilWHC) < 1) {
+    phi = fluxes.drainage / params.soilWHC;
+  } else {
+    phi = 1;
+  }
+  expNLeaching = minN * phi * nLeachFrac;
+  calcNLeachingFlux();
+  status |= checkNLeachingFlux(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);
+  }
+  status |= minStatus;
+
+  return status;
+}
+
 int testNVolatilization(void) {
   int status = 0;
   double minN;
@@ -58,7 +117,7 @@ int testNVolatilization(void) {
 
   minN = 2;
   nVolFrac = 0.1;
-  initState(minN, nVolFrac);
+  initState(minN, nVolFrac, 0);
   double tEffect = calcTempEffect(climate->tsoil);
   double mEffect = calcMoistEffect(envi.soilWater, params.soilWHC);
   expNVolFlux = nVolFrac * minN * tEffect * mEffect;
@@ -68,7 +127,7 @@ int testNVolatilization(void) {
   // easy proportionality test - doubling params should double output
   minN *= 2;
   expNVolFlux *= 2;
-  initState(minN, nVolFrac);
+  initState(minN, nVolFrac, 0);
   calcNVolatilizationFlux();
   status |= checkVolatilizationFlux(expNVolFlux);
 
@@ -92,7 +151,7 @@ int testFertilization(void) {
   double expMinN, expEventMinNFlux, expNVolFlux;
 
   // init minN 2, nVol 0.1
-  initState(initN, nVolFrac);
+  initState(initN, nVolFrac, 0);
 
   // fert event: 15 5 10
   double fertMinN = 10;
@@ -140,6 +199,8 @@ int run(void) {
 
   status |= testNVolatilization();
 
+  status |= testNLeaching();
+
   return status;
 }