SIP128 Add fert event handling (C only) (#131)

* Interim commit

* Updates for fert event, harvest fix, and tests

* Updates for fert event, harvest fix, and tests

* Updates testEventOutputFile after harvest changes

* updates test output for consistency

* Adds test for fertilization event

* Update model structure docs to explain what happens when LITTER_POOL is off (#135)

* updated documentation to describe what happens when LITTER_POOL=0

* remove redundant explanation

---------

Co-authored-by: David LeBauer <[email protected]>

Author: Alomir

docs/model-structure.md

@@ -54,6 +54,17 @@ component out of
 -->
 ## Carbon Dynamics
 
+
+### Litter Pool
+
+SIPNET can be run with or without a separate litter pool (LITTER_POOL=1 or 0). 
+Equations in this document assume LITTER_POOL=1 unless otherwise noted.
+
+When LITTER_POOL=0:
+- Carbon fluxes that would go to the litter pool ($F^C_\text{litter}$) are routed directly to the soil carbon pool
+- All decomposition occurs in the soil pool
+- This affects carbon routing from harvest events, organic matter additions, plant senescence, and other processes involving $F^C_\text{litter}$
+
 ### Maximum Photosynthetic Rate
 
 $$
@@ -168,7 +179,8 @@ $$
 R_\text{wood} = K_\text{wood} \cdot C_\text{wood} \cdot D_{\text{temp,Q10}_v} \tag{Braswell A19}\label{eq:A19}
 $$
 
-Wood maintenance respiration  $(R_m)$ depends on the wood carbon content  $(C_\text{wood})$, a scaling constant  $(k_\text{wood})$, and the temperature sensitivity scaling function $D_{\text{temp,Q10}_v}$.
+Wood maintenance respiration $(R_m)$ depends on the wood carbon content  $(C_\text{wood})$, 
+a scaling constant  $(k_\text{wood})$, and the temperature sensitivity scaling function $D_{\text{temp,Q10}_v}$.
 
 
 ### Litter Carbon
@@ -221,7 +233,7 @@ The rate of decomposition is a function of the litter carbon content and the dec
 ### Soil Carbon
 
 $$
-\frac{dC_\text{soil}}{dt} = F^C_\text{soil} - R_{H_\text{soil}} \tag{Braswell A3}\label{eq:A3}
+\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} - R_{H_\text{soil}} \tag{Braswell A3}\label{eq:A3}
 $$
 
 The change in the SOC pool over time $\frac{dC_\text{soil}}{dt}$ is determined by the addition of litter carbon and the loss of carbon to heterotrophic respiration. This model assumes no loss of SOC to leaching or erosion.
@@ -579,25 +591,30 @@ For the relationship between $N_2O$ flux and soil moisture, Wang et al (2023) su
 
 ## $\frak{Agronomic \ Management \ Events}$
 
-All management events are specified in the `events.in`. Each event is a separate record that includes the date of the event, the type of event, and associated parameters.
+All management events are specified in the `events.in`. Each event is a separate record that includes the 
+date of the event, the type of event, and associated parameters.
 
-### $\frak{Fertilizer \ and \ Organic \ Matter \ Additions}$ 
+### $\frak{Fertilizer}$ and Organic Matter Additions 
 
-Additions of Mineral N, Organic N, and Organic C are represented by the fluxes $F^N_{\text{fert,min}}$, $F^N_{\text{fert,org}},$ and $F^C_{\text{fert,org}}$ that are specified in the `events.in` configuration file.
+Additions of Mineral N, Organic N, and Organic C are added directly to their respective pools via the 
+fluxes $F^N_{\text{fert,min}}$, $F^N_{\text{fert,org}},$ and $F^C_{\text{fert,org}}$ that are specified 
+in the `events.in` configuration file.
 
 Event parameters specified in the `events.in` file:
-
 - Organic N added  $(F^N_{\text{fert,org}})$
 - Organic C added  $(F^C_{\text{fert,org}})$
 - Mineral N added  $(F^N_{\text{fert,min}})$
 
-These are added to the litter C and N and mineral N pools, respectively.
-
-Mineral N includes fertilizer supplied as NO3, NH4, and Urea-N. Urea-N is assumed to hydrolyze to ammonium and bicarbonate rapidly and is treated as a mineral N pool. This is a common assumption because of the high rate of this conversion, and is consistent the DayCent formulation (Parton et al TK-ref, other models and refs?). Only relatively recently did DayCent explicitly model Urea-N to NH4 in order to represent the impact of urease inhibitors (Gurung et al 2021) that slow down the rate.
+Mineral N includes fertilizer supplied as NO3, NH4, and Urea-N. Urea-N is assumed to hydrolyze to ammonium 
+and bicarbonate rapidly and is treated as a mineral N pool. This is a common model assumption because of 
+the fast conversion of Urea to ammonium, and is consistent with the DayCent formulation (Parton et al, 2001). 
+Only relatively recently did DayCent explicitly model Urea-N to NH4 in order to represent the impact of 
+urease inhibitors (Gurung et al 2021) that slow down the rate.
 
 ### $\frak{Tillage}$
 
-To represent tillage, we define two new adjustment factors that modify the decomposition rates of litter $K_{\text{litter}}$ and soil organic matter $K_{\text{som}}$:
+To represent tillage, we define two new adjustment factors that modify the decomposition rates
+ of litter $K_{\text{litter}}$ and soil organic matter $K_{\text{som}}$:
 
 Event parameters from the `events.in` file:
 

src/common/context.c

@@ -137,7 +137,7 @@ void updateCharContext(const char *name, const char *value,
 int hasSourcePrecedence(struct context_metadata *s,
                         context_source_t newSource) {
   // If newSource is greater (or equal) to old source, then it's good
-  return (s->source < newSource);
+  return (s->source <= newSource);
 }
 
 // Get a printable version of source enum for dumpConfig()

src/sipnet/sipnet.c

@@ -1742,6 +1742,7 @@ void processEvents(void) {
         const double fracTA = harvParams->fractionTransferredAbove;
         const double fracRB = harvParams->fractionRemovedBelow;
         const double fracTB = harvParams->fractionTransferredBelow;
+        char carbonPool[20] = "";
 
         // Litter increase
         const double litterAdd = fracTA * (envi.plantLeafC + envi.plantWoodC) +
@@ -1755,15 +1756,22 @@ void processEvents(void) {
         const double coarseDelta = -envi.coarseRootC * (fracRB + fracTB);
 
         // Pool updates:
-        envi.litter += litterAdd;
+        if (ctx.litterPool) {
+          envi.litter += litterAdd;
+          strcpy(carbonPool, "env.litter");
+        } else {
+          // If the litter pool is not turned on, add it to the soil pool
+          envi.soil += litterAdd;
+          strcpy(carbonPool, "env.soil");
+        }
         envi.plantLeafC += leafDelta;
         envi.plantWoodC += woodDelta;
         envi.fineRootC += fineDelta;
         envi.coarseRootC += coarseDelta;
 
         // FUTURE: move/remove biomass in N pools
 
-        writeEventOut(eventOutFile, event, 5, "env.litter", litterAdd,
+        writeEventOut(eventOutFile, event, 5, carbonPool, litterAdd,
                       "envi.plantLeafC", leafDelta, "envi.plantWoodC",
                       woodDelta, "envi.fineRootC", fineDelta,
                       "envi.coarseRootC", coarseDelta);
@@ -1772,10 +1780,32 @@ void processEvents(void) {
         // TBD
         printf("Tillage events not yet implemented\n");
         break;
-      case FERTILIZATION:
-        // TBD
-        printf("Fertilization events not yet implemented\n");
-        break;
+      case FERTILIZATION: {
+        const FertilizationParams *fertParams = event->eventParams;
+        // const double orgN = fertParams->orgN;
+        const double orgC = fertParams->orgC;
+        // const double minN = fertParams->minN;
+        char carbonPool[20] = "";
+
+        // Pool updates:
+        if (ctx.litterPool) {
+          envi.litter += orgC;
+          strcpy(carbonPool, "env.litter");
+          // orgN
+          // minN
+        } else {
+          // If the litter pool is not turned on, add it to the soil pool
+          envi.soil += orgC;
+          strcpy(carbonPool, "env.soil");
+          // orgN
+          // minN
+        }
+
+        // FUTURE: allocate to N pools
+
+        // This will (likely) be 3 params eventually
+        writeEventOut(eventOutFile, event, 1, carbonPool, orgC);
+      } break;
       default:
         printf("Unknown event type (%d) in processEvents()\n", event->type);
         exit(EXIT_CODE_UNKNOWN_EVENT_TYPE_OR_PARAM);

tests/sipnet/test_bugfixes/testEventFileOrderChecks.c

@@ -28,7 +28,7 @@ int run(void) {
   test_assert(jmp_rval == 0);
   status |= !exit_result;
   if (!exit_result) {
-    printf("FAIL with infra_events_year_boundary.in\n");
+    logTest("FAIL with infra_events_year_boundary.in\n");
   }
 
   // Second test - this should exit
@@ -42,7 +42,7 @@ int run(void) {
   test_assert(jmp_rval == 1);
   status |= !exit_result;
   if (!exit_result) {
-    printf("FAIL with infra_events_bad_order.in\n");
+    logTest("FAIL with infra_events_bad_order.in\n");
   }
 
   // Allow a real exit, not that this is really needed
@@ -54,13 +54,13 @@ int run(void) {
 int main(void) {
   int status;
 
-  printf("Starting testEventFileOrderChecks:run()\n");
+  logTest("Starting testEventFileOrderChecks:run()\n");
   status = run();
   if (status) {
     really_exit = 1;
-    printf("FAILED testEventFileOrderChecks with status %d\n", status);
+    logTest("FAILED testEventFileOrderChecks with status %d\n", status);
     exit(status);
   }
 
-  printf("PASSED testEventFileOrderChecks\n");
+  logTest("PASSED testEventFileOrderChecks\n");
 }

tests/sipnet/test_events_infrastructure/events_output_no_header.out

@@ -1,6 +1,6 @@
 2023   65  plant    envi.plantLeafC=10.00,envi.plantWoodC=5.00,envi.fineRootC=4.00,envi.coarseRootC=3.00
 2023   70  irrig    envi.soilWater=5.00,fluxes.immedEvap=0.00
-2023  200  harv     env.litter=12.40,envi.plantLeafC=-4.80,envi.plantWoodC=-3.20,envi.fineRootC=-4.80,envi.coarseRootC=-4.80
+2023  200  harv     env.soil=12.40,envi.plantLeafC=-4.80,envi.plantWoodC=-3.20,envi.fineRootC=-4.80,envi.coarseRootC=-4.80
 2024   65  plant    envi.plantLeafC=10.00,envi.plantWoodC=5.00,envi.fineRootC=4.00,envi.coarseRootC=3.00
 2024   70  irrig    envi.soilWater=2.50,fluxes.immedEvap=2.50
-2024  200  harv     env.litter=12.14,envi.plantLeafC=-10.32,envi.plantWoodC=-5.88,envi.fineRootC=-2.88,envi.coarseRootC=-2.48
+2024  200  harv     env.soil=12.14,envi.plantLeafC=-10.32,envi.plantWoodC=-5.88,envi.fineRootC=-2.88,envi.coarseRootC=-2.48

tests/sipnet/test_events_infrastructure/testEventOutputFile.c

@@ -84,12 +84,14 @@ int run(void) {
   initContext();
   updateIntContext("events", 1, CTX_TEST);
 
+  updateIntContext("litterPool", 0, CTX_TEST);  // default, but to be sure
   status = runTest("events_output_no_header", 0);
   if (status) {
     logTest("runTest(no_header) failed\n");
   }
   testStatus |= status;
 
+  updateIntContext("litterPool", 1, CTX_TEST);
   status = runTest("events_output_header", 1);
   if (status) {
     logTest("runTest(header) failed\n");

tests/sipnet/test_events_types/Makefile

@@ -8,7 +8,7 @@ LDFLAGS=-L$(ROOT_DIR)/libs
 LDLIBS=-lsipnet -lsipnet_common -lm
 
 # List test files in this directory here
-TEST_CFILES=testEventIrrigation.c testEventPlanting.c testEventHarvest.c
+TEST_CFILES=testEventIrrigation.c testEventPlanting.c testEventHarvest.c testEventFertilization.c
 
 # The rest is boilerplate, likely copyable as is to a new test directory
 TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o)

tests/sipnet/test_events_types/events_one_fert.in

@@ -0,0 +1 @@
+2024 70 fert 15  5 10 # fertilized with 15 g/m2 organic N, 5 g/m2 organic C, and 10 g/m2 mineral N on day 40

tests/sipnet/test_events_types/events_two_fert.in

@@ -0,0 +1,2 @@
+2024 70 fert 15  5 10 # fertilized with 15 g/m2 organic N, 5 g/m2 organic C, and 10 g/m2 mineral N on day 40
+2024 70 fert  5  2  3 # fertilized with 15 g/m2 organic N, 5 g/m2 organic C, and 10 g/m2 mineral N on day 40

tests/sipnet/test_events_types/testEventFertilization.c

@@ -0,0 +1,84 @@
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "utils/tUtils.h"
+
+#include "typesUtils.h"
+
+int checkOutput(double litterC) {
+  int status = 0;
+  double curLitterC = 0;
+  if (ctx.litterPool) {
+    logTest("Checking litter pool\n");
+    curLitterC = envi.litter;
+  } else {
+    logTest("Checking soil pool\n");
+    curLitterC = envi.soil;
+    litterC += 0.5;  // We bumped init soil C to distinguish
+  }
+  if (!compareDoubles(litterC, curLitterC)) {
+    logTest("Litter/soil C is %f, expected %f\n", curLitterC, litterC);
+    status = 1;
+  }
+  return status;
+}
+
+void initEnv(void) {
+  envi.soil = 1.5;
+  envi.litter = 1;
+  // Others to be added for N
+}
+
+int run(void) {
+  int status = 0;
+  double expLitterC;
+
+  // We will need to switch back and forth between litter pool and soil manually
+  prepTypesTest();
+
+  // init values
+  initEnv();
+
+  //// ONE PLANTING EVENT
+  updateIntContext("litterPool", 0, CTX_TEST);
+  logTest("Litter pool is %s\n", ctx.litterPool ? "on" : "off");
+  initEvents("events_one_fert.in", 0);
+  setupEvents();
+  processEvents();
+
+  // First fert: (15-5-10)
+  expLitterC = 1 + 5;
+  // litterN + 15
+  // minN + 10
+  status |= checkOutput(expLitterC);
+
+  //// TWO HARVEST EVENTS
+  updateIntContext("litterPool", 1, CTX_TEST);
+  logTest("Litter pool is %s\n", ctx.litterPool ? "on" : "off");
+  initEnv();
+  initEvents("events_two_fert.in", 1);
+  setupEvents();
+  processEvents();
+  // First event same as above (15-5-10)
+  expLitterC = 1 + 5;
+  // litterN
+  // minN
+  // Second fert (5-2-3)
+  expLitterC += 2;
+  // litterN += 5
+  // minN += 3
+  status |= checkOutput(expLitterC);
+
+  return status;
+}
+
+int main(void) {
+  logTest("Starting run()\n");
+  int status = run();
+  if (status) {
+    logTest("FAILED testEventFertilization with status %d\n", status);
+    exit(status);
+  }
+
+  logTest("PASSED testEventFertilization\n");
+}