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

@@ -959,6 +959,50 @@ W_{\text{WHC}} - W_{\text{soil}} + F^W_{\text{irrigation}} & \text{flooding}
 
 -->
 
+## Mass Balance Checks
+
+SIPNET includes internal mass balance checks to ensure conservation of carbon and nitrogen throughout the simulation. These checks verify that changes in pool sizes match the difference between system inputs and outputs.
+
+### Carbon Balance
+
+The carbon balance check verifies that:
+
+\begin{equation}
+\Delta C_\text{pools} = C_\text{inputs} - C_\text{outputs}
+\end{equation}
+
+where:
+
+- $\Delta C_\text{pools}$ is the change in total carbon across all pools (wood, leaf, fine root, coarse root, soil, and litter if enabled)
+- $C_\text{inputs}$ includes photosynthesis (GPP) and carbon additions from agronomic events
+- $C_\text{outputs}$ includes autotrophic respiration ($R_a$: vegetation, fine root, and coarse root), heterotrophic respiration ($R_h$: soil and litter if enabled), and carbon removals from agronomic events
+
+### Nitrogen Balance
+
+When the nitrogen cycle is enabled (`nitrogenCycle=1`), a nitrogen balance check verifies that:
+
+\begin{equation}
+\Delta N_\text{pools} = N_\text{inputs} - N_\text{outputs}
+\end{equation}
+
+where:
+
+- $\Delta N_\text{pools}$ is the change in total nitrogen across all pools (plant biomass, soil organic, litter, and mineral nitrogen)
+- $N_\text{inputs}$ includes nitrogen fixation (when implemented) and nitrogen additions from fertilization and organic matter events
+- $N_\text{outputs}$ includes nitrogen volatilization ($F^N_\text{vol}$), leaching ($F^N_\text{leach}$), and nitrogen removals from harvest events
+
+### Implementation Details
+
+The balance checks are performed at each timestep after all pool updates are complete. The checks account for:
+
+- **Wood Carbon Storage**: To properly track carbon balance with the five-day NPP averaging, the wood pool is split into two components: `plantWoodC` and `plantWoodCStorageDelta`. The storage delta tracks the lag between NPP input and allocation output, ensuring accurate carbon accounting. The total wood carbon is the sum of these two pools.
+
+- **Climate Length Adjustment**: All fluxes are adjusted for the timestep length (`climate.length`) to ensure proper comparison between pool deltas and flux totals.
+
+- **Numerical Precision**: Small discrepancies ($< 10^{-8}$) are treated as zero to avoid false positives from floating-point rounding errors.
+
+If a balance check fails (i.e., $|\Delta C_\text{pools} - (C_\text{inputs} - C_\text{outputs})| > 10^{-8}$), SIPNET logs an internal error with the timestep information and the magnitude of the imbalance.
+
 ## References
 
 Braswell, Bobby H., William J. Sacks, Ernst Linder, and David S. Schimel. 2005. Estimating Diurnal to Annual Ecosystem Parameters by Synthesis of a Carbon Flux Model with Eddy Covariance Net Ecosystem Exchange Observations. Global Change Biology 11 (2): 335–55. https://doi.org/10.1111/j.1365-2486.2005.00897.x.

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 use plantWoodC by itself; it's the
+    // reason plantWoodCStorageDelta 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");