From b6463c9c9c711dc05fe6d7d961800bd4cfe8af14 Mon Sep 17 00:00:00 2001 From: mswilburn Date: Thu, 26 Mar 2026 19:29:40 -0500 Subject: [PATCH 01/13] Adding C saturation flag --- src/common/context.c | 7 +++++++ src/common/context.h | 3 ++- src/sipnet/cli.c | 4 ++++ src/sipnet/cli.h | 2 +- src/sipnet/restart.c | 5 +++++ 5 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/common/context.c b/src/common/context.c index 4477271e6..6ff5c9d3e 100644 --- a/src/common/context.c +++ b/src/common/context.c @@ -42,6 +42,7 @@ void initContext(void) { CREATE_INT_CONTEXT(waterHResp, "WATER_HRESP", ARG_ON, FLAG_YES); CREATE_INT_CONTEXT(nitrogenCycle, "NITROGEN_CYCLE", ARG_OFF, FLAG_YES); CREATE_INT_CONTEXT(anaerobic, "ANAEROBIC", ARG_OFF, FLAG_YES); + CREATE_INT_CONTEXT(carbonSaturation,"CARBON_SATURATION",ARG_OFF, FLAG_YES); // Flags, I/O CREATE_INT_CONTEXT(doMainOutput, "DO_MAIN_OUTPUT", ARG_ON, FLAG_YES); @@ -206,6 +207,12 @@ void validateContext(void) { hasError = 1; } + if (ctx.carbonSaturation && !ctx.litterPool) { + logError("carbon-saturation requires litter-pool to be " + "turned on\n"); + hasError = 1; + } + if (hasError) { exit(EXIT_CODE_BAD_PARAMETER_VALUE); } diff --git a/src/common/context.h b/src/common/context.h index 71e986ec8..d7a740f55 100644 --- a/src/common/context.h +++ b/src/common/context.h @@ -35,7 +35,7 @@ struct context_metadata { UT_hash_handle hh; // makes this structure hashable }; -#define NUM_CONTEXT_MODEL_FLAGS 10 +#define NUM_CONTEXT_MODEL_FLAGS 11 // See docs/developer-guide/cli-options.md for details on how to add a new // Context entry struct Context { @@ -51,6 +51,7 @@ struct Context { int waterHResp; int nitrogenCycle; int anaerobic; + int carbonSaturation; // IF ADDING A NEW MODEL FLAG, update NUM_CONTEXT_MODEL_FLAGS above and // relevant code in restart.c diff --git a/src/sipnet/cli.c b/src/sipnet/cli.c index c2d40cb89..2f4ef0a9c 100644 --- a/src/sipnet/cli.c +++ b/src/sipnet/cli.c @@ -31,6 +31,7 @@ static struct option long_options[] = { // NOLINT DECLARE_FLAG(water-hresp), DECLARE_FLAG(nitrogen-cycle), DECLARE_FLAG(anaerobic), + DECLARE_FLAG(carbon-saturation), DECLARE_FLAG(do-main-output), DECLARE_FLAG(do-single-outputs), @@ -62,6 +63,7 @@ char *argNameMap[] = { DECLARE_ARG_FOR_MAP(litterPool), DECLARE_ARG_FOR_MAP(snow), DECLARE_ARG_FOR_MAP(soilPhenol), DECLARE_ARG_FOR_MAP(waterHResp), DECLARE_ARG_FOR_MAP(nitrogenCycle), DECLARE_ARG_FOR_MAP(anaerobic), + DECLARE_ARG_FOR_MAP(carbonSaturation), // I/O DECLARE_ARG_FOR_MAP(doMainOutput), DECLARE_ARG_FOR_MAP(doSingleOutputs), @@ -92,6 +94,7 @@ void usage(char *progName) { printf(" --snow Keep track of snowpack, rather than assuming all precipitation is liquid (1)\n"); printf(" --soil-phenol Use soil temperature to determine leaf growth (0)\n"); printf(" --water-hresp Whether soil moisture affects heterotrophic respiration (1)\n"); + printf(" --carbon-saturation Enable maximum storage limit of soil organic carbon (0)\n"); 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 .out (1)\n"); @@ -115,6 +118,7 @@ void usage(char *progName) { printf(" --soil-phenol and --gdd may not both be turned on\n"); printf(" --anaerobic requires --water-hresp\n"); printf(" --nitrogen-cycle requires both --litter-pool and --anaerobic\n"); + printf(" --carbon-saturation requires --litter-pool\n"); // clang-format on } diff --git a/src/sipnet/cli.h b/src/sipnet/cli.h index 93fb83342..b1b1ff0fe 100644 --- a/src/sipnet/cli.h +++ b/src/sipnet/cli.h @@ -14,7 +14,7 @@ // The run-time option names do not match their corresponding fields in Context, // so we need a way to get from one to the other. -#define NUM_FLAG_OPTIONS 15 +#define NUM_FLAG_OPTIONS 16 extern char *argNameMap[2 * NUM_FLAG_OPTIONS]; /*! diff --git a/src/sipnet/restart.c b/src/sipnet/restart.c index a503c35da..9fccf6d94 100644 --- a/src/sipnet/restart.c +++ b/src/sipnet/restart.c @@ -86,6 +86,7 @@ typedef struct RestartContextModelFlags { int waterHResp; int nitrogenCycle; int anaerobic; + int carbonSaturation; } RestartContextModelFlags; static RestartContextModelFlags modelFlags; @@ -160,6 +161,7 @@ void initResetState(RestartState *state, MeanTracker *npp) { state->flagsPF[7] = (StateField){"flags.waterHResp", FT_INT, &modelFlags.waterHResp, 0}; state->flagsPF[8] = (StateField){"flags.nitrogenCycle", FT_INT, &modelFlags.nitrogenCycle, 0}; state->flagsPF[9] = (StateField){"flags.anaerobic", FT_INT, &modelFlags.anaerobic, 0}; + state->flagsPF[8] = (StateField){"flags.carbonSaturation", FT_INT, &modelFlags.carbonSaturation, 0}; state->flagsPF[10] = (StateField){"flags.invalid", FT_INVALID, NULL, FIELD_INVALID}; state->boundaryPF[0] = (StateField){"boundary.year", FT_INT, &boundaryClimate.year, 0}; @@ -773,6 +775,7 @@ static void checkRestartContextCompatibility(void) { mismatch |= (ctx.waterHResp != modelFlags.waterHResp); mismatch |= (ctx.nitrogenCycle != modelFlags.nitrogenCycle); mismatch |= (ctx.anaerobic != modelFlags.anaerobic); + mismatch |= (ctx.carbonSaturation != modelFlags.carbonSaturation); if (mismatch) { logError("Restart context mismatch: model flags must match checkpoint " @@ -894,6 +897,8 @@ void restartWriteCheckpoint(const char *restartOut, MeanTracker *meanNPP) { ++numFlagsSet; modelFlags.anaerobic = ctx.anaerobic; ++numFlagsSet; + modelFlags.carbonSaturation = ctx.carbonSaturation; + ++numFlagsSet; if (numFlagsSet != NUM_CONTEXT_MODEL_FLAGS) { logInternalError("Not all model flags set while writing checkpoint\n"); exit(EXIT_CODE_INTERNAL_ERROR); From 54c071a3cf7574773683a7306b24132c10ff78c9 Mon Sep 17 00:00:00 2001 From: mswilburn Date: Fri, 27 Mar 2026 10:32:32 -0500 Subject: [PATCH 02/13] Fix flag numbering error --- src/sipnet/restart.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sipnet/restart.c b/src/sipnet/restart.c index 9fccf6d94..f331a5272 100644 --- a/src/sipnet/restart.c +++ b/src/sipnet/restart.c @@ -161,8 +161,8 @@ void initResetState(RestartState *state, MeanTracker *npp) { state->flagsPF[7] = (StateField){"flags.waterHResp", FT_INT, &modelFlags.waterHResp, 0}; state->flagsPF[8] = (StateField){"flags.nitrogenCycle", FT_INT, &modelFlags.nitrogenCycle, 0}; state->flagsPF[9] = (StateField){"flags.anaerobic", FT_INT, &modelFlags.anaerobic, 0}; - state->flagsPF[8] = (StateField){"flags.carbonSaturation", FT_INT, &modelFlags.carbonSaturation, 0}; - state->flagsPF[10] = (StateField){"flags.invalid", FT_INVALID, NULL, FIELD_INVALID}; + state->flagsPF[10] = (StateField){"flags.carbonSaturation", FT_INT, &modelFlags.carbonSaturation, 0}; + state->flagsPF[11] = (StateField){"flags.invalid", FT_INVALID, NULL, FIELD_INVALID}; state->boundaryPF[0] = (StateField){"boundary.year", FT_INT, &boundaryClimate.year, 0}; state->boundaryPF[1] = (StateField){"boundary.day", FT_INT, &boundaryClimate.day, 0}; From 8e174fcbc3e4505e7cebc0184c447639e533555c Mon Sep 17 00:00:00 2001 From: mswilburn Date: Fri, 8 May 2026 15:55:40 -0500 Subject: [PATCH 03/13] Soil carbon pool calculation and fluxes --- src/sipnet/cli.c | 3 +-- src/sipnet/sipnet.c | 23 +++++++++++++++++++---- src/sipnet/state.h | 15 +++++++++++++++ tests/smoke/russell_1/sipnet.param | 1 + tests/smoke/russell_2/sipnet.param | 1 + tests/smoke/russell_3/sipnet.param | 1 + tests/smoke/russell_4/sipnet.param | 1 + 7 files changed, 39 insertions(+), 6 deletions(-) diff --git a/src/sipnet/cli.c b/src/sipnet/cli.c index c7772dd50..231c835ef 100644 --- a/src/sipnet/cli.c +++ b/src/sipnet/cli.c @@ -66,8 +66,7 @@ char *argNameMap[] = { DECLARE_ARG_FOR_MAP(litterPool), DECLARE_ARG_FOR_MAP(snow), DECLARE_ARG_FOR_MAP(soilPhenol), DECLARE_ARG_FOR_MAP(waterHResp), DECLARE_ARG_FOR_MAP(nitrogenCycle), DECLARE_ARG_FOR_MAP(anaerobic), - DECLARE_ARG_FOR_MAP(flooding), - DECLARE_ARG_FOR_MAP(carbonSaturation), + DECLARE_ARG_FOR_MAP(flooding), DECLARE_ARG_FOR_MAP(carbonSaturation), // I/O DECLARE_ARG_FOR_MAP(doMainOutput), DECLARE_ARG_FOR_MAP(doSingleOutputs), diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 8edfdbec8..6bb3a6ca8 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -388,6 +388,9 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) { // Water drainage initializeOneModelParam(modelParams, "waterDrainFrac", &(params.waterDrainFrac), ctx.flooding); + // Soil carbon saturation + initializeOneModelParam(modelParams, "soilCSaturation", &(params.soilCSaturation), ctx.carbonSaturation); + // NOLINTEND // clang-format on @@ -1937,10 +1940,22 @@ void updatePoolsForSoil(void) { fluxes.rLitter - fluxes.litterMethane) * climate->length; - // from [2] and [3], litter and root terms respectively - envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss + - fluxes.litterToSoil - fluxes.rSoil - fluxes.soilMethane) * - climate->length; + if (ctx.carbonSaturation) { + double soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; + // + envi.soilC += (soilInputs * (1 - envi.soilC / params.soilCSaturation) - + fluxes.rSoil - fluxes.soilMethane) * + climate->length; + // send unstabilized carbon back to litter pool + fluxes.soilToLitter = envi.soilC / params.soilCSaturation; + envi.litterC += fluxes.soilToLitter * climate->length; + } else { + // from [2] and [3], litter and root terms respectively + envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss + + fluxes.litterToSoil - fluxes.rSoil - fluxes.soilMethane) * + climate->length; + } } else { // Normal pool (single pool, no microbes) // :: from [1] (and others, TBD), eq (A3), where: diff --git a/src/sipnet/state.h b/src/sipnet/state.h index 58225aa5f..0b5cb720b 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -387,6 +387,14 @@ typedef struct Parameters { // Relative methane production rate in the litter pool, in [0, 1), per day double litterMethaneRate; + + // ****** + // Soil carbon saturation + // ****** + + // Maximum threshold for stabilizing carbon in soil organic pool as + // slow-turnover pool units: g C * m^-2 ground area + double soilCSaturation; } Params; #define NUM_PARAMS (sizeof(Params) / sizeof(double)) @@ -602,6 +610,13 @@ typedef struct FluxVars { double soilMethane; // Methane produced from litter double litterMethane; + + // **************************************** + // Soil carbon saturation + + // unstabilized soil C that exceeds soil C saturation threshold gets sent back + // to litter C pool as a fast-turnover pool (g C * m^-2 ground area * day^-1) + double soilToLitter; } Fluxes; // Global var diff --git a/tests/smoke/russell_1/sipnet.param b/tests/smoke/russell_1/sipnet.param index 3b1df98db..1566bda4c 100644 --- a/tests/smoke/russell_1/sipnet.param +++ b/tests/smoke/russell_1/sipnet.param @@ -70,3 +70,4 @@ anaerobicTransExp 2.0 nFixationFracMax 0.5 halfNFixationMax 1.0 waterDrainFrac 1.0 +soilCSaturation 1.0 diff --git a/tests/smoke/russell_2/sipnet.param b/tests/smoke/russell_2/sipnet.param index 056a8b2a5..53734091c 100644 --- a/tests/smoke/russell_2/sipnet.param +++ b/tests/smoke/russell_2/sipnet.param @@ -72,3 +72,4 @@ anaerobicTransExp 2.0 soilMethaneRate 0.01 litterMethaneRate 0.01 waterDrainFrac 1.0 +soilCSaturation 1.0 diff --git a/tests/smoke/russell_3/sipnet.param b/tests/smoke/russell_3/sipnet.param index 6e33e3d9a..b6b80ebff 100644 --- a/tests/smoke/russell_3/sipnet.param +++ b/tests/smoke/russell_3/sipnet.param @@ -67,3 +67,4 @@ kCN 80.0 nFixationFracMax 0.5 halfNFixationMax 1.0 waterDrainFrac 1.0 +soilCSaturation 1.0 diff --git a/tests/smoke/russell_4/sipnet.param b/tests/smoke/russell_4/sipnet.param index ac1a3f723..0de7bfec2 100644 --- a/tests/smoke/russell_4/sipnet.param +++ b/tests/smoke/russell_4/sipnet.param @@ -65,3 +65,4 @@ woodCN 100.0 fineRootCN 40.0 kCN 80.0 waterDrainFrac 1.0 +soilCSaturation 1.0 From baae4f2e450d546b7dd35198161f1bab24f5f2fb Mon Sep 17 00:00:00 2001 From: mswilburn Date: Fri, 8 May 2026 16:05:09 -0500 Subject: [PATCH 04/13] Updating flag count --- src/sipnet/cli.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sipnet/cli.h b/src/sipnet/cli.h index b1b1ff0fe..fab31d3e8 100644 --- a/src/sipnet/cli.h +++ b/src/sipnet/cli.h @@ -14,7 +14,7 @@ // The run-time option names do not match their corresponding fields in Context, // so we need a way to get from one to the other. -#define NUM_FLAG_OPTIONS 16 +#define NUM_FLAG_OPTIONS 17 extern char *argNameMap[2 * NUM_FLAG_OPTIONS]; /*! From 443ccb015978399a2c1ebbd23c795718fb94eb2d Mon Sep 17 00:00:00 2001 From: mswilburn Date: Mon, 11 May 2026 10:58:11 -0500 Subject: [PATCH 05/13] Updating flux state variable calculation --- src/sipnet/sipnet.c | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 6bb3a6ca8..e325336d0 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -1309,6 +1309,10 @@ void calcLitterFluxes() { fluxes.rLitter = litterBreakdown * params.fracLitterRespired; fluxes.litterToSoil = litterBreakdown * (1.0 - params.fracLitterRespired); + if (ctx.carbonSaturation) { + // unstabilized carbon returns to litter pool from soil + fluxes.soilToLitter = envi.soilC / params.soilCSaturation; + } } else { // litterBreakdown = 0; fluxes.rLitter = 0; @@ -1934,23 +1938,22 @@ void updateMainPools() { */ void updatePoolsForSoil(void) { if (ctx.litterPool) { - // :: from [2], litter model description - envi.litterC += - (fluxes.woodLitter + fluxes.leafLitter - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - if (ctx.carbonSaturation) { + envi.litterC += + (fluxes.woodLitter + fluxes.leafLitter + fluxes.soilToLitter - + fluxes.litterToSoil - fluxes.rLitter - fluxes.litterMethane) * + climate->length; double soilInputs = fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - // - envi.soilC += (soilInputs * (1 - envi.soilC / params.soilCSaturation) - - fluxes.rSoil - fluxes.soilMethane) * + envi.soilC += (soilInputs * (1 - fluxes.soilToLitter) - fluxes.rSoil - + fluxes.soilMethane) * climate->length; - // send unstabilized carbon back to litter pool - fluxes.soilToLitter = envi.soilC / params.soilCSaturation; - envi.litterC += fluxes.soilToLitter * climate->length; } else { + // :: from [2], litter model description + envi.litterC += + (fluxes.woodLitter + fluxes.leafLitter - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * + climate->length; // from [2] and [3], litter and root terms respectively envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil - fluxes.rSoil - fluxes.soilMethane) * From cf727e2b448b345b2295c60a4b449ddf354f5748 Mon Sep 17 00:00:00 2001 From: mswilburn Date: Tue, 19 May 2026 14:00:50 -0500 Subject: [PATCH 06/13] Create unit test --- src/sipnet/sipnet.c | 15 +- src/sipnet/state.h | 6 - .../test_modeling/testCarbonSaturation.c | 143 ++++++++++++++++++ 3 files changed, 150 insertions(+), 14 deletions(-) create mode 100644 tests/sipnet/test_modeling/testCarbonSaturation.c diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index e325336d0..3319d1e26 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -1309,10 +1309,6 @@ void calcLitterFluxes() { fluxes.rLitter = litterBreakdown * params.fracLitterRespired; fluxes.litterToSoil = litterBreakdown * (1.0 - params.fracLitterRespired); - if (ctx.carbonSaturation) { - // unstabilized carbon returns to litter pool from soil - fluxes.soilToLitter = envi.soilC / params.soilCSaturation; - } } else { // litterBreakdown = 0; fluxes.rLitter = 0; @@ -1939,13 +1935,16 @@ void updateMainPools() { void updatePoolsForSoil(void) { if (ctx.litterPool) { if (ctx.carbonSaturation) { + double saturationFraction = envi.soilC / params.soilCSaturation; + double soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; + envi.litterC += - (fluxes.woodLitter + fluxes.leafLitter + fluxes.soilToLitter - + (fluxes.woodLitter + fluxes.leafLitter + (soilInputs * saturationFraction) - fluxes.litterToSoil - fluxes.rLitter - fluxes.litterMethane) * climate->length; - double soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - envi.soilC += (soilInputs * (1 - fluxes.soilToLitter) - fluxes.rSoil - + + envi.soilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - fluxes.soilMethane) * climate->length; } else { diff --git a/src/sipnet/state.h b/src/sipnet/state.h index 0b5cb720b..fb4fe1a4f 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -611,12 +611,6 @@ typedef struct FluxVars { // Methane produced from litter double litterMethane; - // **************************************** - // Soil carbon saturation - - // unstabilized soil C that exceeds soil C saturation threshold gets sent back - // to litter C pool as a fast-turnover pool (g C * m^-2 ground area * day^-1) - double soilToLitter; } Fluxes; // Global var diff --git a/tests/sipnet/test_modeling/testCarbonSaturation.c b/tests/sipnet/test_modeling/testCarbonSaturation.c new file mode 100644 index 000000000..619c5922e --- /dev/null +++ b/tests/sipnet/test_modeling/testCarbonSaturation.c @@ -0,0 +1,143 @@ +#include "utils/tUtils.h" +#include "sipnet/sipnet.c" + +void resetState(void) { + ctx.litterPool = 1; + ctx.anaerobic = 1; + ctx.carbonSaturation = 1; + + envi.soilWater = 7.5; + + params.soilMethaneRate = 0.05; + params.litterMethaneRate = 0.1; +} + +void setupTests(void) { + // set up dummy climate + climate = (ClimateNode *)malloc(sizeof(ClimateNode)); + climate->year = 2024; + climate->day = 70; + climate->length = 0.125; + climate->tsoil = 20.0; + + // Environment + envi.soilWater = 7.5; + // soil CN 7.5 + envi.soilC = 15.0; + envi.soilOrgN = 2.0; + // litter CN 5.0 + envi.litterC = 7.5; + envi.litterN = 1.5; + + params.soilWHC = 10.0; + params.soilRespQ10 = 3.0; + params.fAnoxia = 0.6; + params.anaerobicDecompRate = 0.5; + params.anaerobicTransExp = 2.0; + + // Initialize general state + resetState(); + + // Set up the context + initContext(); +} + +int checkFlux(double calcFlux, double expFlux, const char *label) { + // Make sure we didn't forget to update context, in case dependencies changed + validateContext(); + + int status = 0; + if (!compareDoubles(calcFlux, expFlux)) { + status = 1; + logTest("Calculated %s flux is %f, expected %f\n", label, calcFlux, + expFlux); + } + + return status; +} + +int checkPool(double calcPool, double expPool, const char *label) { + // Make sure we didn't forget to update context, in case dependencies changed + validateContext(); + + int status = 0; + if (!compareDoubles(calcPool, expPool)) { + status = 1; + logTest("Calculated %s pool is %f, expected %f\n", label, calcPool, + expPool); + } + + return status; +} + +int checkSoil(double flux, double pool) { + int status = 0; + status |= checkFlux(fluxes.soilMethane, flux, "soil"); + status |= checkPool(envi.soilC, pool, "soil"); + return status; +} + +int checkLitter(double flux, double pool) { + int status = 0; + status |= checkFlux(fluxes.litterMethane, flux, "litter"); + status |= checkPool(envi.litterC, pool, "litter"); + return status; +} + +void callCalcMethaneFlux(void) { + // Make sure we didn't make a dumb mistake with tests points + // That is, make sure we HAVE some methane flux + if (envi.soilWater / params.soilWHC <= params.fAnoxia + 0.1) { + logTest("Soil water is not high enough to trigger methane production!\n"); + logTest("SoilW %8.4f WHC %8.4f f_a %8.4f\n", envi.soilWater, params.soilWHC, + params.fAnoxia); + exit(1); + } + + calcMethaneFlux(); +} + +int testCarbonSaturation(void) { + int status = 0; + double expSoilF, expSoilC; + double expLitterF, expLitterC; + + double tempEffect = calcTempEffect(climate->tsoil); + double moistEfffect = calcMethaneMoistEffect(envi.soilWater, params.soilWHC); + + // Standard + logTest(" Test: standard\n"); + resetState(); + callCalcMethaneFlux(); + updatePoolsForSoil(); + expSoilF = 0.75 * tempEffect * moistEfffect; + expLitterF = 0.75 * tempEffect * moistEfffect; + expSoilC = 15 - expSoilF * climate->length; + expLitterC = 7.5 - expLitterF * climate->length; + status |= checkSoil(expSoilF, expSoilC); + status |= checkLitter(expLitterF, expLitterC); + return status; +} + +int run(void) { + int status = 0; + + setupTests(); + + status |= testCarbonSaturation(); + + return status; +} + +int main(void) { + int status; + + logTest("Starting testCarbonSaturation:run()\n"); + status = run(); + if (status) { + logTest("FAILED testCarbonSaturation with status %d\n", status); + exit(status); + } + + logTest("PASSED testCarbonSaturation\n"); +} From 68b472cec01710093957fa7ce6b7e3cba94f6ae8 Mon Sep 17 00:00:00 2001 From: mswilburn Date: Tue, 19 May 2026 14:01:45 -0500 Subject: [PATCH 07/13] clang-format --- src/sipnet/sipnet.c | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 3319d1e26..b12497cb3 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -1937,13 +1937,13 @@ void updatePoolsForSoil(void) { if (ctx.carbonSaturation) { double saturationFraction = envi.soilC / params.soilCSaturation; double soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; + + envi.litterC += (fluxes.woodLitter + fluxes.leafLitter + + (soilInputs * saturationFraction) - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * + climate->length; - envi.litterC += - (fluxes.woodLitter + fluxes.leafLitter + (soilInputs * saturationFraction) - - fluxes.litterToSoil - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - envi.soilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - fluxes.soilMethane) * climate->length; From 43d72a0e01969aab44a60e877a23048d54ea0b1a Mon Sep 17 00:00:00 2001 From: mswilburn Date: Thu, 21 May 2026 12:04:10 -0500 Subject: [PATCH 08/13] Updating unit testing and config files --- tests/sipnet/test_modeling/Makefile | 2 +- .../test_modeling/testCarbonSaturation.c | 168 +++++++++++------- tests/smoke/niwot/sipnet.config | 1 + tests/smoke/russell_1/sipnet.config | 3 +- tests/smoke/russell_2/sipnet.config | 3 +- tests/smoke/russell_3/sipnet.config | 3 +- 6 files changed, 113 insertions(+), 67 deletions(-) diff --git a/tests/sipnet/test_modeling/Makefile b/tests/sipnet/test_modeling/Makefile index 964d4cb46..0110f7eab 100644 --- a/tests/sipnet/test_modeling/Makefile +++ b/tests/sipnet/test_modeling/Makefile @@ -8,7 +8,7 @@ LDFLAGS=-L$(ROOT_DIR)/libs LDLIBS=-lsipnet -lsipnet_common -lm # List test files in this directory here -TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c +TEST_CFILES=testNitrogenCycle.c testDependencyFunctions.c testBalance.c testMethane.c testSoilMoisture.c testCarbonSaturation.c # The rest is boilerplate, likely copyable as is to a new test directory TEST_OBJ_FILES=$(TEST_CFILES:%.c=%.o) diff --git a/tests/sipnet/test_modeling/testCarbonSaturation.c b/tests/sipnet/test_modeling/testCarbonSaturation.c index 619c5922e..dd2f03624 100644 --- a/tests/sipnet/test_modeling/testCarbonSaturation.c +++ b/tests/sipnet/test_modeling/testCarbonSaturation.c @@ -3,13 +3,19 @@ void resetState(void) { ctx.litterPool = 1; - ctx.anaerobic = 1; ctx.carbonSaturation = 1; - envi.soilWater = 7.5; + envi.litterC = 10; - params.soilMethaneRate = 0.05; - params.litterMethaneRate = 0.1; + params.soilCSaturation = 10; + + fluxes.woodLitter = 0; + fluxes.leafLitter = 0; + fluxes.litterToSoil = 0; + fluxes.rLitter = 0; + fluxes.litterMethane = 0; + fluxes.rSoil = 0; + fluxes.soilMethane = 0; } void setupTests(void) { @@ -20,21 +26,6 @@ void setupTests(void) { climate->length = 0.125; climate->tsoil = 20.0; - // Environment - envi.soilWater = 7.5; - // soil CN 7.5 - envi.soilC = 15.0; - envi.soilOrgN = 2.0; - // litter CN 5.0 - envi.litterC = 7.5; - envi.litterN = 1.5; - - params.soilWHC = 10.0; - params.soilRespQ10 = 3.0; - params.fAnoxia = 0.6; - params.anaerobicDecompRate = 0.5; - params.anaerobicTransExp = 2.0; - // Initialize general state resetState(); @@ -42,20 +33,6 @@ void setupTests(void) { initContext(); } -int checkFlux(double calcFlux, double expFlux, const char *label) { - // Make sure we didn't forget to update context, in case dependencies changed - validateContext(); - - int status = 0; - if (!compareDoubles(calcFlux, expFlux)) { - status = 1; - logTest("Calculated %s flux is %f, expected %f\n", label, calcFlux, - expFlux); - } - - return status; -} - int checkPool(double calcPool, double expPool, const char *label) { // Make sure we didn't forget to update context, in case dependencies changed validateContext(); @@ -70,52 +47,117 @@ int checkPool(double calcPool, double expPool, const char *label) { return status; } -int checkSoil(double flux, double pool) { +int checkSoil(double pool) { int status = 0; - status |= checkFlux(fluxes.soilMethane, flux, "soil"); status |= checkPool(envi.soilC, pool, "soil"); return status; } -int checkLitter(double flux, double pool) { +int checkLitter(double pool) { int status = 0; - status |= checkFlux(fluxes.litterMethane, flux, "litter"); status |= checkPool(envi.litterC, pool, "litter"); return status; } -void callCalcMethaneFlux(void) { - // Make sure we didn't make a dumb mistake with tests points - // That is, make sure we HAVE some methane flux - if (envi.soilWater / params.soilWHC <= params.fAnoxia + 0.1) { - logTest("Soil water is not high enough to trigger methane production!\n"); - logTest("SoilW %8.4f WHC %8.4f f_a %8.4f\n", envi.soilWater, params.soilWHC, - params.fAnoxia); - exit(1); - } - - calcMethaneFlux(); -} - int testCarbonSaturation(void) { int status = 0; - double expSoilF, expSoilC; - double expLitterF, expLitterC; + double expSoilC; + double expLitterC; + double saturationFraction; + double soilInputs; - double tempEffect = calcTempEffect(climate->tsoil); - double moistEfffect = calcMethaneMoistEffect(envi.soilWater, params.soilWHC); + // Low C inputs, low saturationFraction + logTest(" Test: low soilInputs, low satFrac\n"); + resetState(); + + expSoilC = 2.5; // 25% saturation + envi.soilC = 2.5; + expLitterC = 10; + soilInputs = 100; + saturationFraction = expSoilC / params.soilCSaturation; + + expLitterC += (fluxes.woodLitter + fluxes.leafLitter + + (soilInputs * saturationFraction) - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * + climate->length; - // Standard - logTest(" Test: standard\n"); + expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - + fluxes.soilMethane) * + climate->length; + + updatePoolsForSoil(); + status |= checkSoil(expSoilC); + status |= checkLitter(expLitterC); + + // High C inputs, low saturationFraction + logTest(" Test: high soilInputs, low satFrac\n"); + resetState(); + + expSoilC = 2.5; // 25% saturation + envi.soilC = 2.5; + expLitterC = 10; + soilInputs = 200; + saturationFraction = expSoilC / params.soilCSaturation; + + expLitterC += (fluxes.woodLitter + fluxes.leafLitter + + (soilInputs * saturationFraction) - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * + climate->length; + + expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - + fluxes.soilMethane) * + climate->length; + + updatePoolsForSoil(); + status |= checkSoil(expSoilC); + status |= checkLitter(expLitterC); + + // Low C inputs, high saturationFraction + logTest(" Test: low soilInputs, high satFrac\n"); + resetState(); + + expSoilC = 7.5; // 75% saturation + envi.soilC = 7.5; + expLitterC = 10; + soilInputs = 100; + saturationFraction = expSoilC / params.soilCSaturation; + + expLitterC += (fluxes.woodLitter + fluxes.leafLitter + + (soilInputs * saturationFraction) - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * + climate->length; + + expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - + fluxes.soilMethane) * + climate->length; + + updatePoolsForSoil(); + status |= checkSoil(expSoilC); + status |= checkLitter(expLitterC); + + // High C inputs, high saturationFraction + logTest(" Test: high soilInputs, high satFrac\n"); resetState(); - callCalcMethaneFlux(); + + expSoilC = 7.5; // 75% saturation + envi.soilC = 7.5; + expLitterC = 10; + soilInputs = 200; + saturationFraction = expSoilC / params.soilCSaturation; + + expLitterC += (fluxes.woodLitter + fluxes.leafLitter + + (soilInputs * saturationFraction) - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * + climate->length; + + expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - + fluxes.soilMethane) * + climate->length; + updatePoolsForSoil(); - expSoilF = 0.75 * tempEffect * moistEfffect; - expLitterF = 0.75 * tempEffect * moistEfffect; - expSoilC = 15 - expSoilF * climate->length; - expLitterC = 7.5 - expLitterF * climate->length; - status |= checkSoil(expSoilF, expSoilC); - status |= checkLitter(expLitterF, expLitterC); + status |= checkSoil(expSoilC); + status |= checkLitter(expLitterC); + return status; } diff --git a/tests/smoke/niwot/sipnet.config b/tests/smoke/niwot/sipnet.config index abd352f8c..590efc17c 100644 --- a/tests/smoke/niwot/sipnet.config +++ b/tests/smoke/niwot/sipnet.config @@ -1,4 +1,5 @@ ANAEROBIC DEFAULT 0 + CARBON_SATURATION DEFAULT 0 CLIM_FILE CALCULATED sipnet.clim DO_MAIN_OUTPUT INPUT_FILE 1 DO_SINGLE_OUTPUT INPUT_FILE 0 diff --git a/tests/smoke/russell_1/sipnet.config b/tests/smoke/russell_1/sipnet.config index 881e02bc5..c0db7fe1a 100644 --- a/tests/smoke/russell_1/sipnet.config +++ b/tests/smoke/russell_1/sipnet.config @@ -1,6 +1,7 @@ -Final config for SIPNET run at 2026-04-24 03:35:15 UTC +Final config for SIPNET run at 2026-05-21 16:24:01 UTC Name Source Value ANAEROBIC DEFAULT 0 + CARBON_SATURATION DEFAULT 0 CLIM_FILE CALCULATED sipnet.clim DO_MAIN_OUTPUT DEFAULT 1 DO_SINGLE_OUTPUT DEFAULT 0 diff --git a/tests/smoke/russell_2/sipnet.config b/tests/smoke/russell_2/sipnet.config index cc39144f0..6b2201dd1 100644 --- a/tests/smoke/russell_2/sipnet.config +++ b/tests/smoke/russell_2/sipnet.config @@ -1,6 +1,7 @@ -Final config for SIPNET run at 2026-04-24 03:35:15 UTC +Final config for SIPNET run at 2026-05-21 16:24:01 UTC Name Source Value ANAEROBIC INPUT_FILE 1 + CARBON_SATURATION DEFAULT 0 CLIM_FILE CALCULATED sipnet.clim DO_MAIN_OUTPUT INPUT_FILE 1 DO_SINGLE_OUTPUT INPUT_FILE 0 diff --git a/tests/smoke/russell_3/sipnet.config b/tests/smoke/russell_3/sipnet.config index 5e178df9a..27d465f25 100644 --- a/tests/smoke/russell_3/sipnet.config +++ b/tests/smoke/russell_3/sipnet.config @@ -1,6 +1,7 @@ -Final config for SIPNET run at 2026-04-24 03:35:15 UTC +Final config for SIPNET run at 2026-05-21 16:24:01 UTC Name Source Value ANAEROBIC DEFAULT 0 + CARBON_SATURATION DEFAULT 0 CLIM_FILE CALCULATED sipnet.clim DO_MAIN_OUTPUT DEFAULT 1 DO_SINGLE_OUTPUT DEFAULT 0 From 243c2755749f347c48d803bd277817635261440f Mon Sep 17 00:00:00 2001 From: mswilburn Date: Fri, 22 May 2026 07:33:06 -0500 Subject: [PATCH 09/13] Update to unit test --- .../test_modeling/testCarbonSaturation.c | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/tests/sipnet/test_modeling/testCarbonSaturation.c b/tests/sipnet/test_modeling/testCarbonSaturation.c index dd2f03624..a87bd8c30 100644 --- a/tests/sipnet/test_modeling/testCarbonSaturation.c +++ b/tests/sipnet/test_modeling/testCarbonSaturation.c @@ -16,6 +16,7 @@ void resetState(void) { fluxes.litterMethane = 0; fluxes.rSoil = 0; fluxes.soilMethane = 0; + fluxes.fineRootLoss = 0; } void setupTests(void) { @@ -73,7 +74,10 @@ int testCarbonSaturation(void) { expSoilC = 2.5; // 25% saturation envi.soilC = 2.5; expLitterC = 10; - soilInputs = 100; + fluxes.coarseRootLoss = 100; + + soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; saturationFraction = expSoilC / params.soilCSaturation; expLitterC += (fluxes.woodLitter + fluxes.leafLitter + @@ -96,7 +100,10 @@ int testCarbonSaturation(void) { expSoilC = 2.5; // 25% saturation envi.soilC = 2.5; expLitterC = 10; - soilInputs = 200; + fluxes.coarseRootLoss = 200; + + soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; saturationFraction = expSoilC / params.soilCSaturation; expLitterC += (fluxes.woodLitter + fluxes.leafLitter + @@ -119,7 +126,10 @@ int testCarbonSaturation(void) { expSoilC = 7.5; // 75% saturation envi.soilC = 7.5; expLitterC = 10; - soilInputs = 100; + fluxes.coarseRootLoss = 100; + + soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; saturationFraction = expSoilC / params.soilCSaturation; expLitterC += (fluxes.woodLitter + fluxes.leafLitter + @@ -142,7 +152,10 @@ int testCarbonSaturation(void) { expSoilC = 7.5; // 75% saturation envi.soilC = 7.5; expLitterC = 10; - soilInputs = 200; + fluxes.coarseRootLoss = 200; + + soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; saturationFraction = expSoilC / params.soilCSaturation; expLitterC += (fluxes.woodLitter + fluxes.leafLitter + From 7ff37328855c35e05fbf85c4c60a785cf62c1b8b Mon Sep 17 00:00:00 2001 From: mswilburn Date: Thu, 28 May 2026 11:37:13 -0500 Subject: [PATCH 10/13] Updating documentation and changelog --- docs/CHANGELOG.md | 1 + docs/model-structure.md | 16 ++++++++++++++++ docs/parameters.md | 1 + 3 files changed, 18 insertions(+) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 877e668cb..11e3971ad 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -28,6 +28,7 @@ sections to include in release notes: - `sipnet-view` tool for visualizing SIPNET output files (#317) - `leafon` and `leafoff` events for tracking phenological transitions (#326) +- `carbonSaturation` flag and calculation of soil and litter carbon pools that observes soil carbon saturating behavior (#301) ### Fixed diff --git a/docs/model-structure.md b/docs/model-structure.md index 80f506c47..b07c6da18 100644 --- a/docs/model-structure.md +++ b/docs/model-structure.md @@ -309,6 +309,13 @@ F^C_\text{fert,org} Where $K_{\text{plant},i}$ is the turnover rate of plant pool $i$ that controls the rate at which plant biomass is transferred to litter. +If soil carbon saturation is enabled, the litter carbon pool has an additional input of carbon that was not stabilized for long-term storage by the soil. This functionality is described in more detail below in the Soil Carbon section \eqref{eq:soil_carbon_saturation}. + +\begin{equation} +\frac{dC_\text{litter}}{dt} = F^C_\text{litter} + F^C_{\text{soil}} \cdot \frac{C_{\text{soil}}}{f_{\text{C,soil,saturation}}} - F^C_{\text{decomp}} - F^C_{\text{CH}_4\text{,litter}} +\label{eq:soil_carbon_to_litter} +\end{equation} + $F^C_{\text{decomp}}$ represents the rate at which litter carbon is processed by microbial activity. Litter decomposition is modeled as a first-order process proportional to litter carbon content and modified by temperature and moisture: @@ -370,6 +377,15 @@ F^C_{\text{soil}} = F^C_{\text{soil,litter}} + F^C_{\text{soil,roots}} \label{eq:soil_carbon_flux} \end{equation} +If soil carbon saturation is enabled, only a fraction of the total carbon input will be added to the soil. The amount added to the soil carbon pool is a function of the proximity of the pool to the soil carbon saturation limit. As the soil carbon pool aqcuires more carbon and approaches its saturation limit, the amount of carbon stored in the pool decreases. This represents the physical and chemical constraints of the mineralogical surfaces of soil to stabilize carbon for long-term storage. The remaining carbon that is not stabilized in the soil carbon pool due to saturation constraints will be returned to the litter carbon pool to act as fast-turnover carbon \eqref{eq:soil_carbon_to_litter}. + +\begin{equation} +\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} \cdot (1 - \frac{C_{\text{soil}}}{f_{\text{C,soil,saturation}}}) - R_{\text{soil}} - F^C_{\text{CH}_4\text{,soil}} +\label{eq:soil_carbon_saturation} +\end{equation} + +where $f_{\text{C,soil,saturation}}$ is the soil carbon saturation limit entered as an input paramter. + Soil heterotrophic respiration is modeled as a first-order process proportional to soil organic carbon content and modified by environmental and management factors: diff --git a/docs/parameters.md b/docs/parameters.md index 41318140d..91135d20a 100644 --- a/docs/parameters.md +++ b/docs/parameters.md @@ -197,6 +197,7 @@ Run-time parameters can change from one run to the next, or when the model is st | $Q_{10s}$ | soilRespQ10 | Soil respiration Q10 | unitless | scalar determining effect of temp on soil respiration | | $D_{\text{moisture}}$ | soilRespMoistEffect | scalar determining effect of moisture on soil resp. | unitless | | | $f_{\text{till}}$ | tillageEff | Effect of tillage on decomposition that exponentially decays over time | fraction | Documented in model structure; event-level term in `events.in` | +| $f_{\text{C,soil,saturation}}$ | soilCSaturation | Maximum amount of carbon that can be stabilized in the soil determined by soil texture | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ | | ### Nitrogen Cycle Parameters From 3055e6476eb67e455f946fa8383594b250b2356a Mon Sep 17 00:00:00 2001 From: mswilburn Date: Thu, 28 May 2026 12:11:42 -0500 Subject: [PATCH 11/13] Updating more documentation --- docs/user-guide/model-inputs.md | 1 + docs/user-guide/running-sipnet.md | 2 ++ 2 files changed, 3 insertions(+) diff --git a/docs/user-guide/model-inputs.md b/docs/user-guide/model-inputs.md index b1be9896c..ad4eff7a6 100644 --- a/docs/user-guide/model-inputs.md +++ b/docs/user-guide/model-inputs.md @@ -242,6 +242,7 @@ Thus, command-line arguments override settings in the configuration file, and co | `snow` | on | Keep track of snowpack, rather than assuming all precipitation is liquid | | `soil-phenol` | off | Use soil temperature to determine leaf growth | | `water-hresp` | on | Whether soil moisture affects heterotrophic respiration | +| `carbon-saturation`| off | Enable soil carbon saturation behavior to constrain carbon stored in soil | Note the following restrictions on these options: - `soil-phenol` and `gdd` may not both be turned on diff --git a/docs/user-guide/running-sipnet.md b/docs/user-guide/running-sipnet.md index c2f2e5771..4b5e02bfb 100644 --- a/docs/user-guide/running-sipnet.md +++ b/docs/user-guide/running-sipnet.md @@ -48,6 +48,7 @@ These flags enable or disable optional model processes. Prepend `no-` to the fla | `--snow` | ON (1) | Track snowpack separately; if disabled, all precipitation is treated as liquid | | `--soil-phenol` | OFF (0) | Use soil temperature (instead of growing degree days) to determine leaf growth | | `--water-hresp` | ON (1) | Allow soil moisture to affect heterotrophic respiration rates | +| `--carbon-saturation`| OFF (0) | Enable soil carbon saturation behavior to constrain carbon stored in soil | #### Model Flag Restrictions @@ -56,6 +57,7 @@ The following flag constraints are enforced: - `--soil-phenol` and `--gdd` cannot both be enabled - `--anaerobic` requires `--water-hresp` - `--nitrogen-cycle` requires both `--litter-pool` and `--anaerobic` +- `--carbon-saturation` requires `--litter-pool` ### Output Flags From 8acb7c6d79b22e2c3315f2b3872790cdaef76319 Mon Sep 17 00:00:00 2001 From: mswilburn Date: Fri, 29 May 2026 19:00:37 -0500 Subject: [PATCH 12/13] Revised documentation --- docs/model-structure.md | 10 +++++----- docs/parameters.md | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/model-structure.md b/docs/model-structure.md index b07c6da18..aed428f08 100644 --- a/docs/model-structure.md +++ b/docs/model-structure.md @@ -309,10 +309,10 @@ F^C_\text{fert,org} Where $K_{\text{plant},i}$ is the turnover rate of plant pool $i$ that controls the rate at which plant biomass is transferred to litter. -If soil carbon saturation is enabled, the litter carbon pool has an additional input of carbon that was not stabilized for long-term storage by the soil. This functionality is described in more detail below in the Soil Carbon section \eqref{eq:soil_carbon_saturation}. +When soil carbon saturation is enabled, a fraction of soil carbon inputs may be redirected to the litter pool as fast-turnover carbon. This functionality is described in more detail below in the Soil Carbon section \eqref{eq:soil_carbon_saturation}. \begin{equation} -\frac{dC_\text{litter}}{dt} = F^C_\text{litter} + F^C_{\text{soil}} \cdot \frac{C_{\text{soil}}}{f_{\text{C,soil,saturation}}} - F^C_{\text{decomp}} - F^C_{\text{CH}_4\text{,litter}} +\frac{dC_\text{litter}}{dt} = F^C_\text{litter} + F^C_{\text{soil}} \cdot \frac{C_{\text{soil}}}{C_{\text{soil,saturation}}} - F^C_{\text{decomp}} - F^C_{\text{CH}_4\text{,litter}} \label{eq:soil_carbon_to_litter} \end{equation} @@ -377,14 +377,14 @@ F^C_{\text{soil}} = F^C_{\text{soil,litter}} + F^C_{\text{soil,roots}} \label{eq:soil_carbon_flux} \end{equation} -If soil carbon saturation is enabled, only a fraction of the total carbon input will be added to the soil. The amount added to the soil carbon pool is a function of the proximity of the pool to the soil carbon saturation limit. As the soil carbon pool aqcuires more carbon and approaches its saturation limit, the amount of carbon stored in the pool decreases. This represents the physical and chemical constraints of the mineralogical surfaces of soil to stabilize carbon for long-term storage. The remaining carbon that is not stabilized in the soil carbon pool due to saturation constraints will be returned to the litter carbon pool to act as fast-turnover carbon \eqref{eq:soil_carbon_to_litter}. +When soil carbon saturation is enabled, only a saturation-dependent fraction of gross soil C inputs is added to the soil pool. This fraction declines as $C_{\text{soil}}$ approaches the specified soil C saturation limit. The remaining input C is redirected to the litter pool \eqref{eq:soil_carbon_to_litter} as fast-turnover carbon rather than being added to the soil pool. \begin{equation} -\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} \cdot (1 - \frac{C_{\text{soil}}}{f_{\text{C,soil,saturation}}}) - R_{\text{soil}} - F^C_{\text{CH}_4\text{,soil}} +\frac{dC_\text{soil}}{dt} = F^C_{\text{soil}} \cdot (1 - \frac{C_{\text{soil}}}{C_{\text{soil,saturation}}}) - R_{\text{soil}} - F^C_{\text{CH}_4\text{,soil}} \label{eq:soil_carbon_saturation} \end{equation} -where $f_{\text{C,soil,saturation}}$ is the soil carbon saturation limit entered as an input paramter. +where $C_{\text{soil,saturation}}$ is the soil carbon saturation limit entered as an input paramter. This is based on equation (3) from Stewart et al. (2007). Soil heterotrophic respiration is modeled as a first-order process proportional to soil organic carbon content and modified by environmental and management factors: diff --git a/docs/parameters.md b/docs/parameters.md index 91135d20a..15cfe1411 100644 --- a/docs/parameters.md +++ b/docs/parameters.md @@ -197,7 +197,7 @@ Run-time parameters can change from one run to the next, or when the model is st | $Q_{10s}$ | soilRespQ10 | Soil respiration Q10 | unitless | scalar determining effect of temp on soil respiration | | $D_{\text{moisture}}$ | soilRespMoistEffect | scalar determining effect of moisture on soil resp. | unitless | | | $f_{\text{till}}$ | tillageEff | Effect of tillage on decomposition that exponentially decays over time | fraction | Documented in model structure; event-level term in `events.in` | -| $f_{\text{C,soil,saturation}}$ | soilCSaturation | Maximum amount of carbon that can be stabilized in the soil determined by soil texture | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ | | +| $C_{\text{soil,saturation}}$ | soilCSaturation | Maximum amount of carbon that can be stabilized in the soil determined by soil texture | $\text{g C} \cdot \text{m}^{-2} \text{ ground area}$ | | ### Nitrogen Cycle Parameters From 2ba7714c7d561d07f9bb6df922e164af08b5601c Mon Sep 17 00:00:00 2001 From: mswilburn Date: Fri, 29 May 2026 20:02:48 -0500 Subject: [PATCH 13/13] Capping satFrac and limiting org N --- src/sipnet/nitrogen.c | 24 +++-- src/sipnet/sipnet.c | 44 +++++---- .../test_modeling/testCarbonSaturation.c | 90 ++++++++----------- 3 files changed, 76 insertions(+), 82 deletions(-) diff --git a/src/sipnet/nitrogen.c b/src/sipnet/nitrogen.c index 626a807c4..d0af41c92 100644 --- a/src/sipnet/nitrogen.c +++ b/src/sipnet/nitrogen.c @@ -51,22 +51,32 @@ static void calcNPoolFluxes(void) { double litterMin = fluxes.rLitter / litterCN; double soilMin = fluxes.rSoil / soilCN; + // Adding soil carbon saturation functionality so organic N fluxes to soil + // and litter are proportional to respective carbon fluxes dependent on + // soil carbon saturation + double soilNInputs = fluxes.litterToSoil / litterCN + + fluxes.fineRootLoss / params.fineRootCN + + fluxes.coarseRootLoss / params.woodCN; + // saturationFraction capped between zero and one + double saturationFraction = + (ctx.carbonSaturation) + ? (fmin(1.0, fmax(0.0, envi.soilC / params.soilCSaturation))) + : 0.0; + // litter // The litter org N flux is determined by the carbon fluxes from wood and leaf // litter (modified by leaf N resorption), and N loss due to mineralization. // N added via fertilization is handled elsewhere. - fluxes.nOrgLitter = fluxes.leafLitter / params.leafCN - - fluxes.leafOffNResorption + - fluxes.woodLitter / params.woodCN - litterMin - - fluxes.litterToSoil / litterCN; + fluxes.nOrgLitter = + fluxes.leafLitter / params.leafCN - fluxes.leafOffNResorption + + fluxes.woodLitter / params.woodCN - litterMin - + fluxes.litterToSoil / litterCN + (soilNInputs * saturationFraction); // soil // The soil org N flux is determined by the carbon flux from the litter pool, // carbon fluxes from roots, and N loss due to mineralization // (Note: woodCN is used for coarse roots) - fluxes.nOrgSoil = fluxes.litterToSoil / litterCN - soilMin + - fluxes.fineRootLoss / params.fineRootCN + - fluxes.coarseRootLoss / params.woodCN; + fluxes.nOrgSoil = soilNInputs * (1 - saturationFraction) - soilMin; // mineralization fluxes.nMin = litterMin + soilMin; diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 93fd30fa4..ce4b8652c 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -1650,30 +1650,28 @@ void updateMainPools(void) { */ void updatePoolsForSoil(void) { if (ctx.litterPool) { - if (ctx.carbonSaturation) { - double saturationFraction = envi.soilC / params.soilCSaturation; - double soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - - envi.litterC += (fluxes.woodLitter + fluxes.leafLitter + - (soilInputs * saturationFraction) - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - - envi.soilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - - fluxes.soilMethane) * + double soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; + // Adding soil carbon saturation functionality + // soilFraction capped between zero and one + // if envi.soilC > params.soilCSaturation (due to initialization or events) + // saturationFraction will return 1 and no additional soil C inputs will be + // accepted to envi.soilC until losses due to respiration and methane reduce + // pool below saturation level, all inputs will be directed to envi.litterC + // pool while envi.soilC is above saturation level + double saturationFraction = + (ctx.carbonSaturation) + ? (fmin(1.0, fmax(0.0, envi.soilC / params.soilCSaturation))) + : 0.0; + // :: from [2], litter model description + envi.litterC += (fluxes.woodLitter + fluxes.leafLitter + + (soilInputs * saturationFraction) - fluxes.litterToSoil - + fluxes.rLitter - fluxes.litterMethane) * climate->length; - } else { - // :: from [2], litter model description - envi.litterC += - (fluxes.woodLitter + fluxes.leafLitter - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - // from [2] and [3], litter and root terms respectively - envi.soilC += (fluxes.coarseRootLoss + fluxes.fineRootLoss + - fluxes.litterToSoil - fluxes.rSoil - fluxes.soilMethane) * - climate->length; - } + // from [2] and [3], litter and root terms respectively + envi.soilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - + fluxes.soilMethane) * + climate->length; } else { // Normal pool (single pool, no microbes) // :: from [1] (and others, TBD), eq (A3), where: diff --git a/tests/sipnet/test_modeling/testCarbonSaturation.c b/tests/sipnet/test_modeling/testCarbonSaturation.c index a87bd8c30..6d703d1bf 100644 --- a/tests/sipnet/test_modeling/testCarbonSaturation.c +++ b/tests/sipnet/test_modeling/testCarbonSaturation.c @@ -14,7 +14,6 @@ void resetState(void) { fluxes.litterToSoil = 0; fluxes.rLitter = 0; fluxes.litterMethane = 0; - fluxes.rSoil = 0; fluxes.soilMethane = 0; fluxes.fineRootLoss = 0; } @@ -60,12 +59,26 @@ int checkLitter(double pool) { return status; } +void calcCSaturation(double *expSoilC, double *expLitterC) { + double saturationFraction; + double soilInputs = + fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; + + saturationFraction = + (ctx.carbonSaturation) + ? (fmin(1.0, fmax(0.0, *expSoilC / params.soilCSaturation))) + : 0.0; + + *expLitterC += soilInputs * saturationFraction * climate->length; + + *expSoilC += + (soilInputs * (1 - saturationFraction) - fluxes.rSoil) * climate->length; +} + int testCarbonSaturation(void) { int status = 0; double expSoilC; double expLitterC; - double saturationFraction; - double soilInputs; // Low C inputs, low saturationFraction logTest(" Test: low soilInputs, low satFrac\n"); @@ -75,20 +88,9 @@ int testCarbonSaturation(void) { envi.soilC = 2.5; expLitterC = 10; fluxes.coarseRootLoss = 100; + fluxes.rSoil = 0; - soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - saturationFraction = expSoilC / params.soilCSaturation; - - expLitterC += (fluxes.woodLitter + fluxes.leafLitter + - (soilInputs * saturationFraction) - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - - expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - - fluxes.soilMethane) * - climate->length; - + calcCSaturation(&expSoilC, &expLitterC); updatePoolsForSoil(); status |= checkSoil(expSoilC); status |= checkLitter(expLitterC); @@ -101,20 +103,9 @@ int testCarbonSaturation(void) { envi.soilC = 2.5; expLitterC = 10; fluxes.coarseRootLoss = 200; + fluxes.rSoil = 0; - soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - saturationFraction = expSoilC / params.soilCSaturation; - - expLitterC += (fluxes.woodLitter + fluxes.leafLitter + - (soilInputs * saturationFraction) - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - - expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - - fluxes.soilMethane) * - climate->length; - + calcCSaturation(&expSoilC, &expLitterC); updatePoolsForSoil(); status |= checkSoil(expSoilC); status |= checkLitter(expLitterC); @@ -127,20 +118,9 @@ int testCarbonSaturation(void) { envi.soilC = 7.5; expLitterC = 10; fluxes.coarseRootLoss = 100; + fluxes.rSoil = 0; - soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - saturationFraction = expSoilC / params.soilCSaturation; - - expLitterC += (fluxes.woodLitter + fluxes.leafLitter + - (soilInputs * saturationFraction) - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; - - expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - - fluxes.soilMethane) * - climate->length; - + calcCSaturation(&expSoilC, &expLitterC); updatePoolsForSoil(); status |= checkSoil(expSoilC); status |= checkLitter(expLitterC); @@ -153,21 +133,27 @@ int testCarbonSaturation(void) { envi.soilC = 7.5; expLitterC = 10; fluxes.coarseRootLoss = 200; + fluxes.rSoil = 0; - soilInputs = - fluxes.coarseRootLoss + fluxes.fineRootLoss + fluxes.litterToSoil; - saturationFraction = expSoilC / params.soilCSaturation; + calcCSaturation(&expSoilC, &expLitterC); + updatePoolsForSoil(); + status |= checkSoil(expSoilC); + status |= checkLitter(expLitterC); - expLitterC += (fluxes.woodLitter + fluxes.leafLitter + - (soilInputs * saturationFraction) - fluxes.litterToSoil - - fluxes.rLitter - fluxes.litterMethane) * - climate->length; + // Soil C Pool > saturationFraction + logTest(" Test: soilC > saturation\n"); + resetState(); - expSoilC += (soilInputs * (1 - saturationFraction) - fluxes.rSoil - - fluxes.soilMethane) * - climate->length; + expSoilC = 12.5; // 125% saturation + envi.soilC = 12.5; + expLitterC = 10; + fluxes.coarseRootLoss = 200; + fluxes.rSoil = 50; + calcCSaturation(&expSoilC, &expLitterC); updatePoolsForSoil(); + // logTest(" expSoilC is %f\n", expSoilC); + // logTest(" expLitterC is %f\n", expLitterC); status |= checkSoil(expSoilC); status |= checkLitter(expLitterC);