Skip to content

SIP310 Implement H2O mass-balance checking #343

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ANAMASGARD wants to merge 2 commits into
base: master
Choose a base branch
from
Open

SIP310 Implement H2O mass-balance checking #343

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8f4987b
Add water mass-balance checking alongside existing C/N checks.
ANAMASGARD May 27, 2026
566c8bd
Fix smoke baselines to match clean rebuild output.
ANAMASGARD May 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
57 changes: 57 additions & 0 deletions src/sipnet/balance.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,11 @@
// Definition of global balance tracker struct
BalanceTracker balanceTracker;

// No leaf-water pool is tracked in envi; ctx.leafWater only caps the immedEvap
// flux rate (see calcPrecip()). Excess intercepted water goes to soil via
// netRain.
static double getWaterTotal(void) { return envi.soilWater + envi.snow; }

void getMassTotals(double *carbon, double *nitrogen) {
*carbon = (envi.plantWoodC + envi.plantWoodCStorageDelta) + envi.plantLeafC +
envi.fineRootC + envi.coarseRootC + envi.soilC;
Expand All @@ -32,15 +37,18 @@ void getMassTotals(double *carbon, double *nitrogen) {
void updateBalanceTrackerPreUpdate(void) {
// Set the pre-update pool totals
getMassTotals(&balanceTracker.preTotalC, &balanceTracker.preTotalN);
balanceTracker.preTotalWater = getWaterTotal();
}

void updateBalanceTrackerPostUpdate(void) {
// Set the post-update pool totals
getMassTotals(&balanceTracker.postTotalC, &balanceTracker.postTotalN);
balanceTracker.postTotalWater = getWaterTotal();
}

void updateBalanceTrackerPostClamp(void) {
getMassTotals(&balanceTracker.finalC, &balanceTracker.finalN);
balanceTracker.finalWater = getWaterTotal();

// Difference between post-update and post-clamp is the amount we "gained" by
// setting negative values to zero
Expand All @@ -62,6 +70,16 @@ void updateBalanceTrackerPostClamp(void) {
balanceTracker.clampedN = 0;
}

balanceTracker.clampedWater =
balanceTracker.finalWater - balanceTracker.postTotalWater;
if (balanceTracker.clampedWater < -EPS) {
// This shouldn't happen, by construction
logInternalError("Non-negative clamping has cause water loss\n");
}
if (balanceTracker.clampedWater < EPS) {
balanceTracker.clampedWater = 0;
}

// Calculate the system inputs and outputs
// CARBON
balanceTracker.inputsC = fluxes.photosynthesis + // GPP
Expand Down Expand Up @@ -89,11 +107,26 @@ void updateBalanceTrackerPostClamp(void) {
balanceTracker.outputsN *= climate->length;
}

// WATER
// eventEvap: irrigation water immediately evaporated; counted as both input
// (water entered the system) and output (left as vapour). Net to soil is
// eventSoilWater only. snowMelt is excluded (internal snow -> soil transfer).
balanceTracker.inputsWater =
fluxes.rain + fluxes.snowFall + fluxes.eventSoilWater + fluxes.eventEvap;
balanceTracker.outputsWater =
fluxes.immedEvap + fluxes.transpiration + fluxes.evaporation +
fluxes.sublimation + fluxes.drainage + fluxes.fastFlow + fluxes.eventEvap;

// Account for climate length
balanceTracker.inputsWater *= climate->length;
balanceTracker.outputsWater *= climate->length;

// Account for gains from clamping
balanceTracker.inputsC += balanceTracker.clampedC;
if (ctx.nitrogenCycle) {
balanceTracker.inputsN += balanceTracker.clampedN;
}
balanceTracker.inputsWater += balanceTracker.clampedWater;
}

void initBalanceTracker(void) {
Expand All @@ -107,16 +140,23 @@ void initBalanceTracker(void) {
balanceTracker.postTotalN = 0.0;
balanceTracker.clampedN = 0.0;
balanceTracker.finalN = 0.0;
balanceTracker.preTotalWater = 0.0;
balanceTracker.postTotalWater = 0.0;
balanceTracker.clampedWater = 0.0;
balanceTracker.finalWater = 0.0;

// System I/O
balanceTracker.inputsC = 0.0;
balanceTracker.outputsC = 0.0;
balanceTracker.inputsN = 0.0;
balanceTracker.outputsN = 0.0;
balanceTracker.inputsWater = 0.0;
balanceTracker.outputsWater = 0.0;

// Checks
balanceTracker.deltaC = 0.0;
balanceTracker.deltaN = 0.0;
balanceTracker.deltaWater = 0.0;
}

void checkBalance(void) {
Expand All @@ -134,13 +174,25 @@ void checkBalance(void) {
double systemNDelta = balanceTracker.inputsN - balanceTracker.outputsN;
balanceTracker.deltaN = poolNDelta - systemNDelta;

// WATER
// Pool delta
double poolWaterDelta =
balanceTracker.finalWater - balanceTracker.preTotalWater;
// System delta
double systemWaterDelta =
balanceTracker.inputsWater - balanceTracker.outputsWater;
balanceTracker.deltaWater = poolWaterDelta - systemWaterDelta;

// To avoid weird negative-zero issues...
if (fabs(balanceTracker.deltaC) < EPS) {
balanceTracker.deltaC = 0.0;
}
if (fabs(balanceTracker.deltaN) < EPS) {
balanceTracker.deltaN = 0.0;
}
if (fabs(balanceTracker.deltaWater) < EPS) {
balanceTracker.deltaWater = 0.0;
}

int err = 0;
if (fabs(balanceTracker.deltaC) > 0.0) {
Expand All @@ -162,6 +214,11 @@ void checkBalance(void) {
balanceTracker.inputsN, balanceTracker.outputsN,
balanceTracker.clampedN, balanceTracker.deltaN);
}
if (fabs(balanceTracker.deltaWater) > 0.0) {
logWarning(
"Water balance check failed (delta=%8.4f, Y: %d D: %d T: %4.2f)\n",
balanceTracker.deltaWater, climate->year, climate->day, climate->time);
}
if (err) {
logInternalError("Exiting\n");
exit(EXIT_CODE_INTERNAL_ERROR);
Expand Down
9 changes: 9 additions & 0 deletions src/sipnet/balance.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,15 +28,24 @@ typedef struct BalanceTrackerStruct {
double inputsN;
double outputsN;

// Water balance
double preTotalWater;
double postTotalWater;
double finalWater;
double inputsWater;
double outputsWater;

// Stock clamping adjustments: mass added by clamping negative stocks to zero.
// Positive values indicate that stocks were clamped (mass was "created" to
// prevent negative pools).
double clampedC;
double clampedN;
double clampedWater;

// Checks
double deltaC;
double deltaN;
double deltaWater;
} BalanceTracker;

// Global var
Expand Down
6 changes: 4 additions & 2 deletions src/sipnet/sipnet.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -437,7 +437,7 @@ void outputHeader(FILE *out) {
fprintf(out,
"fluxestranspiration minN soilOrgN litterN "
"plantStorageN n2o nLeaching nFixation nUptake ch4 "
"nppStorage\n");
"nppStorage bcdeltaC bcdeltaN bcdeltaW\n");
}

/*!
Expand Down Expand Up @@ -467,7 +467,9 @@ void outputState(FILE *out, int year, int day, double time) {
fprintf(out, "%9.6f %9.4f %10.4f %8.4f %8.4f", trackers.n2o,
trackers.nLeaching, trackers.nFixation, trackers.nUptake,
trackers.methane);
fprintf(out, "%12.4f\n", envi.plantWoodCStorageDelta);
fprintf(out, "%12.4f %9.5f %9.5f %9.5f\n", envi.plantWoodCStorageDelta,
balanceTracker.deltaC, balanceTracker.deltaN,
balanceTracker.deltaWater);
}

// de-allocate space used for climate linked list
Expand Down
30 changes: 29 additions & 1 deletion tests/sipnet/test_modeling/testBalance.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
// For help with debugging failures
#define DEBUG_C 0
#define DEBUG_N 0
#define DEBUG_W 0
#define DO_OUTPUT 0
#define OUTPUT_FILE "balance.out"

Expand Down Expand Up @@ -74,6 +75,19 @@ int checkNitrogen(void) {
return status;
}

int checkWater(void) {
validateContext();

int status = 0;
if (fabs(balanceTracker.deltaWater) > TEST_EPS) {
status = 1;
logTest("water balance check delta %8.5f (Y %d D %d T %4.2f)\n",
balanceTracker.deltaWater, climate->year, climate->day,
climate->time);
}
return status;
}

int testBalanceSimple(void) {
logTest("Starting testBalanceSimple()\n");
int status = 0;
Expand All @@ -83,6 +97,7 @@ int testBalanceSimple(void) {

status |= checkCarbon();
status |= checkNitrogen();
status |= checkWater();

return status;
}
Expand All @@ -101,6 +116,7 @@ int testBalanceLeaf(void) {

status |= checkCarbon();
status |= checkNitrogen();
status |= checkWater();

climate = climate->nextClim;
}
Expand All @@ -121,6 +137,7 @@ int testBalanceLeafEvents(void) {

status |= checkCarbon();
status |= checkNitrogen();
status |= checkWater();

climate = climate->nextClim;
}
Expand All @@ -142,6 +159,7 @@ int testBalanceNoLitterPool(void) {
// No nitrogen with no litter pool, but this is a good check that we
// aren't mishandling the case
status |= checkNitrogen();
status |= checkWater();

climate = climate->nextClim;
}
Expand Down Expand Up @@ -191,7 +209,7 @@ int main(void) {
}

void dumpState(void) {
#if DEBUG_C | DEBUG_N
#if DEBUG_C | DEBUG_N | DEBUG_W
double len = climate->length;
#endif
#if DEBUG_C
Expand Down Expand Up @@ -237,4 +255,14 @@ void dumpState(void) {
fluxes.fineRootLoss * len / params.fineRootCN,
fluxes.coarseRootLoss * len / params.woodCN);
#endif
#if DEBUG_W
logTest("preW %10.6f postW %10.6f inW %10.6f outW %10.6f pDelta %10.6f "
"sDelta %10.6f delta %10.6f\n",
balanceTracker.preTotalWater, balanceTracker.postTotalWater,
balanceTracker.inputsWater, balanceTracker.outputsWater,
balanceTracker.postTotalWater - balanceTracker.preTotalWater,
balanceTracker.inputsWater - balanceTracker.outputsWater,
balanceTracker.deltaWater);
logTest(" W: soil %10.6f snow %10.6f\n", envi.soilWater, envi.snow);
#endif
}
Loading
Loading