Address PR #248 review feedback: document wood carbon storage pool

.gitignore

@@ -9,6 +9,8 @@ transpose
 .mkdocs.stamp
 site/*
 site_preview/*
+_codeql_build_dir/
+_codeql_detected_source_root
 
 # documentation
 docs/api

CMakeLists.txt

@@ -37,6 +37,7 @@ add_library(sipnetlib
         src/sipnet/runmean.c
         src/sipnet/sipnet.c
         src/sipnet/state.c
+        src/sipnet/balance.c
 )
 
 add_library(tests
@@ -52,4 +53,6 @@ add_library(tests
         tests/sipnet/test_events_types/testEventTillage.c
         tests/sipnet/test_bugfixes/testEventFileOrderChecks.c
         tests/sipnet/test_modeling/testNitrogenCycle.c
+        tests/sipnet/test_modeling/testDependencyFunctions.c
+        tests/sipnet/test_modeling/testBalance.c
 )

Makefile

@@ -11,7 +11,7 @@ COMMON_CFILES:=context.c logging.c modelParams.c util.c
 COMMON_CFILES:=$(addprefix src/common/, $(COMMON_CFILES))
 COMMON_OFILES=$(COMMON_CFILES:.c=.o)
 
-SIPNET_CFILES:=sipnet.c cli.c events.c frontend.c outputItems.c runmean.c state.c
+SIPNET_CFILES:=sipnet.c cli.c events.c frontend.c outputItems.c runmean.c state.c balance.c
 SIPNET_CFILES:=$(addprefix src/sipnet/, $(SIPNET_CFILES))
 SIPNET_OFILES=$(SIPNET_CFILES:.c=.o)
 SIPNET_LIBS=-lsipnet_common

docs/CHANGELOG.md

@@ -40,6 +40,9 @@ sections to include in release notes:
 - Organic N handling for fertilization and soil dynamics (#199)
 - Dynamic C:N for soil and litter (#214)
 - Updates for soil and litter respiration (tillage and C:N effects) (#214)
+- Carbon and nitrogen mass balance checks (#248)
+- Event tracking for leaf-on and leaf-off phenological transitions (#248)
+- Split wood carbon pool into `plantWoodC` and `plantWoodCStorageDelta` to track NPP storage lag (#248)
 
 ### Fixed
 

docs/model-structure.md

@@ -132,25 +132,32 @@ Plant death is implemented as a harvest event with the fraction of biomass trans
 
 As stated above, SIPNET uses a five-day averaged NPP when allocating gained carbon to plant growth. To implement this, 
 the current timestep's net primary production (adjusted GPP - autotrophic respiration) is added to the wood carbon pool
-where it acts as an _implicit_ storage pool, and all allocations from the averaged NPP are deducted from that pool.
-We can represent this storage of carbon conceptually as:
+where it acts as a storage pool, and all allocations from the averaged NPP are deducted from that pool.
+
+Starting in SIPNET v2.1, to support accurate mass balance tracking, this storage is explicitly tracked as a separate pool called $C_{\text{wood,storage}}$ (`plantWoodCStorageDelta` in code). We can represent this storage of carbon as:
 
 \begin{equation}
-NPP_\text{storage} = (GPP - R_a) - \overline{\text{NPP}}_\text{alloc}
+C_{\text{wood,storage}} = (GPP - R_a) - \overline{\text{NPP}}_\text{alloc}
 \end{equation}
 
-where $\overline{NPP}_\text{alloc}$ is the sum of the carbon allocated to the biomass pools as growth. Note that we do not explicitly track this storage term.
+where $\overline{NPP}_\text{alloc}$ is the sum of the carbon allocated to the biomass pools as growth. This storage term represents the lag between NPP input and allocation output due to the five-day averaging.
+
+The total wood carbon is the sum of the structural wood carbon and the storage pool:
+
+\begin{equation}
+C_{\text{wood,total}} = C_{\text{wood}} + C_{\text{wood,storage}}
+\end{equation}
 
 Thus, changes to wood carbon over time are determined by:
 
 \begin{equation}
-\frac{dC_\text{wood}}{dt} = NPP_\text{storage} + \alpha_\text{wood} \cdot \overline{\text{NPP}} - F^C_\text{litter,wood}
+\frac{dC_\text{wood}}{dt} = C_{\text{wood,storage}} + \alpha_\text{wood} \cdot \overline{\text{NPP}} - F^C_\text{litter,wood}
 \label{eq:Braswell_A1}
 \end{equation}
 
 where $\alpha_\text{wood}\cdot\overline{\text{NPP}}$ represents the amount of carbon allocated to growth and $(F^C_\text{litter,wood})$ is the wood litter production.
 
-This is equation (A1) from Braswell, et al. (2005), manipulated to use NPP_\text{storage}.
+This is equation (A1) from Braswell, et al. (2005), modified to explicitly track the storage component.
 
 ### Leaf Carbon
 

docs/parameters.md

@@ -92,6 +92,7 @@ Run-time parameters can change from one run to the next, or when the model is st
 |     | Symbol                     | Parameter Name  | Definition                                                               | Units                                                | notes                                                                              |
 | --- | -------------------------- | --------------- | ------------------------------------------------------------------------ | ---------------------------------------------------- | ---------------------------------------------------------------------------------- |
 | 1   | $C_{\text{wood},0}$        | plantWoodInit   | Initial wood carbon                                                      | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ | above-ground + roots                                                               |
+|     | $C_{\text{wood,storage}}$  |                 | Wood carbon storage pool (delta)                                         | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ | Tracks lag between NPP input and allocation output due to 5-day averaging; initialized to 0; can be negative. In code: `plantWoodCStorageDelta`. Total wood C = $C_{\text{wood}} + C_{\text{wood,storage}}$ |
 | 2   | $LAI_0$                    | laiInit         | Initial leaf area                                                        | m^2 leaves \* m^-2 ground area                       | multiply by SLW to get initial plant leaf C: $C_{\text{leaf},0} = LAI_0 \cdot SLW$ |
 | 3   | $C_{\text{litter},0}$      | litterInit      | Initial litter carbon                                                    | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ |                                                                                    |
 | 4   | $C_{\text{soil},0}$        | soilInit        | Initial soil carbon                                                      | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ |                                                                                    |

src/sipnet/balance.c

@@ -0,0 +1,124 @@
+#include "balance.h"
+
+#include <math.h>
+
+#include "common/context.h"
+#include "common/exitCodes.h"
+#include "common/logging.h"
+#include "state.h"
+
+// Definition of global balance tracker struct
+BalanceTracker balanceTracker;
+
+void getMassTotals(double *carbon, double *nitrogen) {
+  *carbon = (envi.plantWoodC + envi.plantWoodCStorageDelta) + envi.plantLeafC +
+            envi.fineRootC + envi.coarseRootC + envi.soilC;
+  if (ctx.litterPool) {
+    *carbon += envi.litterC;
+  }
+
+  if (ctx.nitrogenCycle) {
+    // Note: this is the one place where we ust plantWoodC by itself; it's the
+    // reason plantNSCWoodCDelta was created, so that we can ignore it here.
+    *nitrogen =
+        envi.plantWoodC / params.woodCN + envi.plantLeafC / params.leafCN +
+        envi.fineRootC / params.fineRootCN + envi.coarseRootC / params.woodCN +
+        envi.soilOrgN + envi.litterN + envi.minN;
+  } else {
+    *nitrogen = 0.0;
+  }
+}
+
+void updateBalanceTrackerPreUpdate(void) {
+  // Set the pre-update pool totals
+  getMassTotals(&balanceTracker.preTotalC, &balanceTracker.preTotalN);
+}
+
+void updateBalanceTrackerPostUpdate(void) {
+  // Set the post-update pool totals
+  getMassTotals(&balanceTracker.postTotalC, &balanceTracker.postTotalN);
+
+  // Calculate the system inputs and outputs
+  // CARBON
+  balanceTracker.inputsC = fluxes.photosynthesis +  // GPP
+                           fluxes.eventInputC;  // agro event additions
+  balanceTracker.outputsC = fluxes.rVeg + fluxes.rFineRoot +
+                            fluxes.rCoarseRoot +  // R_a
+                            fluxes.rSoil +  // R_h
+                            fluxes.eventOutputC;  // agro event removals
+  if (ctx.litterPool) {
+    balanceTracker.outputsC += fluxes.rLitter;
+  }
+  // Account for climate length
+  balanceTracker.inputsC *= climate->length;
+  balanceTracker.outputsC *= climate->length;
+
+  // NITROGEN
+  if (ctx.nitrogenCycle) {
+    balanceTracker.inputsN =
+        // TODO: fluxes.fixation +
+        fluxes.eventInputN;
+    balanceTracker.outputsN =
+        fluxes.nLeaching + fluxes.nVolatilization + fluxes.eventOutputN;
+
+    // Account for climate length
+    balanceTracker.inputsN *= climate->length;
+    balanceTracker.outputsN *= climate->length;
+  }
+}
+
+void initBalanceTracker(void) {
+  // Initialize all to zero
+  balanceTracker.preTotalC = 0.0;
+  balanceTracker.postTotalC = 0.0;
+  balanceTracker.inputsC = 0.0;
+  balanceTracker.outputsC = 0.0;
+  balanceTracker.preTotalN = 0.0;
+  balanceTracker.postTotalN = 0.0;
+  balanceTracker.inputsN = 0.0;
+  balanceTracker.outputsN = 0.0;
+  balanceTracker.deltaC = 0.0;
+  balanceTracker.deltaN = 0.0;
+}
+
+void checkBalance(void) {
+  // CARBON
+  // Pool delta
+  double poolCDelta = balanceTracker.postTotalC - balanceTracker.preTotalC;
+  // System delta
+  double systemCDelta = balanceTracker.inputsC - balanceTracker.outputsC;
+  balanceTracker.deltaC = poolCDelta - systemCDelta;
+
+  // NITROGEN
+  // Pool delta
+  double poolNDelta = balanceTracker.postTotalN - balanceTracker.preTotalN;
+  // System delta
+  double systemNDelta = balanceTracker.inputsN - balanceTracker.outputsN;
+  balanceTracker.deltaN = poolNDelta - systemNDelta;
+
+  // To avoid weird negative-zero issues...
+  if (fabs(balanceTracker.deltaC) < 1e-8) {
+    balanceTracker.deltaC = 0.0;
+  }
+  if (fabs(balanceTracker.deltaN) < 1e-8) {
+    balanceTracker.deltaN = 0.0;
+  }
+
+  int err = 0;
+  if (fabs(balanceTracker.deltaC) > 0.0) {
+    err = 1;
+    logInternalError(
+        "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 (err) {
+    logInternalError("Exiting\n");
+    //  exit(EXIT_CODE_INTERNAL_ERROR);
+  }
+}

src/sipnet/balance.h

@@ -0,0 +1,41 @@
+#ifndef BALANCE_H
+#define BALANCE_H
+
+typedef struct BalanceTrackerStruct {
+  // Mass balance checks:
+  //   X_t - X_(t-1) = inputs - outputs + tolerance
+  // or, we check that
+  //   (pool delta) + (system delta) < tolerance
+  // where
+  //   pool delta = X_t - X_(t-1)
+  //   system delta = outputs - inputs
+
+  // Carbon balance
+  double preTotalC;
+  double postTotalC;
+  double inputsC;
+  double outputsC;
+
+  // Nitrogen balance
+  double preTotalN;
+  double postTotalN;
+  double inputsN;
+  double outputsN;
+
+  // Checks
+  double deltaC;
+  double deltaN;
+} BalanceTracker;
+
+// Global var
+extern BalanceTracker balanceTracker;
+
+void updateBalanceTrackerPreUpdate(void);
+
+void updateBalanceTrackerPostUpdate(void);
+
+void initBalanceTracker(void);
+
+void checkBalance(void);
+
+#endif  // BALANCE_H

src/sipnet/cli.c

@@ -89,7 +89,7 @@ void usage(char *progName) {
   printf("\n");
   printf("Output flags: (prepend flag with 'no-' to force off, eg '--no-print-header')\n");
   printf("  --do-main-output     Print time series of all output variables to <file-name>.out (1)\n");
-  printf("  --do-single-outputs  Print outputs one variable per file (e.g. <file-name>.NEE)\n");
+  printf("  --do-single-outputs  Print selection* of outputs one variable per file (e.g. <file-name>.NEE)\n");
   printf("  --dump-config        Print final config to <file-name>.config (0)\n");
   printf("  --print-header       Whether to print header row in output files (1)\n");
   printf("  --quiet              Suppress info and warning message (0)\n");
@@ -98,6 +98,8 @@ void usage(char *progName) {
   printf("  -h, --help           Print this message and exit\n");
   printf("  -v, --version        Print version information and exit\n");
   printf("\n");
+  printf("*do-single-outputs option outputs are: NEE, NEE_cum, GPP, GPP_cum\n");
+  printf("\n");
   printf("Configuration options are read from <input_file>. Other options specified on the command\n");
   printf("line override settings from that file.\n");
   printf("\n");