From ab563eb18a822380239286e218a3102a9e22941b Mon Sep 17 00:00:00 2001 From: mswilburn Date: Mon, 24 Nov 2025 21:07:33 -0600 Subject: [PATCH 1/2] Adding N gain from fixation --- src/sipnet/sipnet.c | 20 ++- src/sipnet/state.h | 5 + .../sipnet/test_modeling/testNitrogenCycle.c | 141 ++++++++++++++++-- tests/smoke/russell_2/sipnet.param | 1 + 4 files changed, 151 insertions(+), 16 deletions(-) diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 17cb6fdbb..3358a5079 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -390,6 +390,7 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) { initializeOneModelParam(modelParams, "mineralNInit", &(params.minNInit), ctx.nitrogenCycle); initializeOneModelParam(modelParams, "nVolatilizationFrac", &(params.nVolatilizationFrac), ctx.nitrogenCycle); initializeOneModelParam(modelParams, "nLeachingFrac", &(params.nLeachingFrac), ctx.nitrogenCycle); + initializeOneModelParam(modelParams, "nFixationFrac", &(params.nFixationFrac), ctx.nitrogenCycle); // NOLINTEND // clang-format on @@ -412,7 +413,7 @@ void outputHeader(FILE *out) { fprintf(out, "litter soilWater soilWetnessFrac snow "); fprintf(out, "npp nee cumNEE gpp rAboveground rSoil rRoot ra rh rtot " - "evapotranspiration fluxestranspiration minN n2oFlux nLeachFlux\n"); + "evapotranspiration fluxestranspiration minN n2oFlux nLeachFlux nFixationFlux\n"); } /*! @@ -433,12 +434,12 @@ void outputState(FILE *out, int year, int day, double time) { trackers.soilWetnessFrac, envi.snow); fprintf(out, "%8.2f %8.2f %8.2f %8.2f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.8f " - "%8.4f %8.3f %8.6f %8.4f\n", + "%8.4f %8.3f %8.6f %8.4f %8.6f\n", trackers.npp, trackers.nee, trackers.totNee, trackers.gpp, trackers.rAboveground, trackers.rSoil, trackers.rRoot, trackers.ra, trackers.rh, trackers.rtot, trackers.evapotranspiration, fluxes.transpiration, envi.minN, fluxes.nVolatilization, - fluxes.nLeaching); + fluxes.nLeaching, fluxes.nFixation); } // de-allocate space used for climate linked list @@ -1241,6 +1242,16 @@ void calcNLeachingFlux() { fluxes.nLeaching = envi.minN * phi * params.nLeachingFrac; } +/*! + * Calculate mineral N fixation flux + */ +void calcNFixationFlux() { + // flux = k_fix * NPP, g N * m^-2 * day^-1 + // k_fix represents C cost to fix N in range 0.07-0.17 gN/gC + + fluxes.nFixation = params.nFixationFrac * getMeanTrackerMean(meanNPP); +} + /*! * Calculate flux terms for sipnet as part of main model flow * @@ -1310,6 +1321,7 @@ void calculateFluxes(void) { // Leaching depends on drainage flux so makes sure calcNLeachingFlux // occurs after calcSoilWaterFluxes calcNLeachingFlux(); + calcNFixationFlux(); } } @@ -1598,7 +1610,7 @@ void updatePoolsForSoil(void) { climate->length; // Nitrogen Cycle - envi.minN -= (fluxes.nVolatilization + fluxes.nLeaching) * climate->length; + envi.minN += (fluxes.nFixation - fluxes.nVolatilization - fluxes.nLeaching) * climate->length; } // !!! main runner function !!! diff --git a/src/sipnet/state.h b/src/sipnet/state.h index d107033d7..328e11f1e 100644 --- a/src/sipnet/state.h +++ b/src/sipnet/state.h @@ -367,6 +367,9 @@ typedef struct Parameters { // Fraction of mineral N lost to leaching per day double nLeachingFrac; + // Fraction of NPP available to be fixed by vegetation + double nFixationFrac; + } Params; #define NUM_PARAMS (sizeof(Params) / sizeof(double)) @@ -535,6 +538,8 @@ typedef struct FluxVars { double nVolatilization; // Mineral N lost to leaching double nLeaching; + // Mineral N gained through plant fixation + double nFixation; // **************************************** // Fluxes for event handling diff --git a/tests/sipnet/test_modeling/testNitrogenCycle.c b/tests/sipnet/test_modeling/testNitrogenCycle.c index e3999dc0a..e94371b1f 100644 --- a/tests/sipnet/test_modeling/testNitrogenCycle.c +++ b/tests/sipnet/test_modeling/testNitrogenCycle.c @@ -17,7 +17,7 @@ void setupTests(void) { ctx.litterPool = 1; ctx.nitrogenCycle = 1; } - +/* void initState(double initN, double nVol, double nLeachFrac) { // per test envi.minN = initN; @@ -40,7 +40,84 @@ void initState(double initN, double nVol, double nLeachFrac) { params.baseSoilResp = 0.06; params.soilRespQ10 = 2.9; } +*/ +void initNLeachingState(double initN, double nLeachFrac) { + // per test + envi.minN = initN; + params.nLeachingFrac = nLeachFrac; + + // static + envi.soilWater = 5.0; + envi.soil = 1.5; + envi.litter = 1; + + params.minNInit = 0.0; + params.soilWHC = 10.0; + + fluxes.nVolatilization = 0; + fluxes.nLeaching = 0; + + // Values from russell_2 smoke test + params.soilRespMoistEffect = 1.0; + params.baseSoilResp = 0.06; + params.soilRespQ10 = 2.9; +} +void initNVolatilizationState(double initN, double nVol) { + // per test + envi.minN = initN; + params.nVolatilizationFrac = nVol; + + // static + envi.soilWater = 5.0; + envi.soil = 1.5; + envi.litter = 1; + + params.minNInit = 0.0; + params.soilWHC = 10.0; + + fluxes.nVolatilization = 0; + + // Values from russell_2 smoke test + params.soilRespMoistEffect = 1.0; + params.baseSoilResp = 0.06; + params.soilRespQ10 = 2.9; +} + +void initNFixationState(double initN, double nFixFrac) { + // per test + envi.minN = initN; + params.nFixationFrac = nFixFrac; + + // static + envi.soilWater = 5.0; + envi.soil = 1.5; + envi.litter = 1; + + params.minNInit = 0.0; + params.soilWHC = 10.0; + + fluxes.nVolatilization = 0; + fluxes.nLeaching = 0; + fluxes.nFixation = 0; + + // Values from russell_2 smoke test + params.soilRespMoistEffect = 1.0; + params.baseSoilResp = 0.06; + params.soilRespQ10 = 2.9; +} + +int checkFlux(double calcFlux, double expFlux) { + int status = 0; + if (!compareDoubles(calcFlux, expFlux)) { + status = 1; + logTest("Calculated N flux is %f, expected %f\n", + calcFlux, expFlux); + } + + return status; +} +/* int checkVolatilizationFlux(double expNVolFlux) { int status = 0; if (!compareDoubles(fluxes.nVolatilization, expNVolFlux)) { @@ -62,6 +139,40 @@ int checkNLeachingFlux(double expNLeachingFlux) { return status; } +*/ + +int testNFixation(void) { + int status = 0; + double minN; + double nFixFrac; + double expNFixation; + + minN = 1; + nFixFrac = 0.25; + initNFixationState(minN, nFixFrac); + expNFixation = nFixFrac; + calcNFixationFlux(); + status |= checkFlux(fluxes.nFixation, expNFixation); + + minN = 1; + nFixFrac = 0.5; + initNFixationState(minN, nFixFrac); + expNFixation = nFixFrac; + calcNFixationFlux(); + status |= checkFlux(fluxes.nFixation, expNFixation); + + // Check minN for the last + updatePoolsForSoil(); + double expMinN = minN + (expNFixation * climate->length); + int minStatus = 0; + if (!compareDoubles(envi.minN, expMinN)) { + minStatus = 1; + logTest("Fixation minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); + } + status |= minStatus; + + return status; +} int testNLeaching(void) { int status = 0; @@ -73,7 +184,7 @@ int testNLeaching(void) { minN = 1; nLeachFrac = 0.5; fluxes.drainage = 5; - initState(minN, 0, nLeachFrac); + initNLeachingState(minN, nLeachFrac); if ((fluxes.drainage / params.soilWHC) < 1) { phi = fluxes.drainage / params.soilWHC; } else { @@ -81,12 +192,13 @@ int testNLeaching(void) { } expNLeaching = minN * phi * nLeachFrac; calcNLeachingFlux(); - status |= checkNLeachingFlux(expNLeaching); + //status |= checkNLeachingFlux(expNLeaching); + status |= checkFlux(fluxes.nLeaching, expNLeaching); minN = 1; nLeachFrac = 0.5; fluxes.drainage = 20; - initState(minN, 0, nLeachFrac); + initNLeachingState(minN, nLeachFrac); if ((fluxes.drainage / params.soilWHC) < 1) { phi = fluxes.drainage / params.soilWHC; } else { @@ -94,7 +206,8 @@ int testNLeaching(void) { } expNLeaching = minN * phi * nLeachFrac; calcNLeachingFlux(); - status |= checkNLeachingFlux(expNLeaching); + //status |= checkNLeachingFlux(expNLeaching); + status |= checkFlux(fluxes.nLeaching, expNLeaching); // Check minN for the last updatePoolsForSoil(); @@ -102,7 +215,7 @@ int testNLeaching(void) { int minStatus = 0; if (!compareDoubles(envi.minN, expMinN)) { minStatus = 1; - logTest("minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); + logTest("Leaching minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); } status |= minStatus; @@ -117,19 +230,21 @@ int testNVolatilization(void) { minN = 2; nVolFrac = 0.1; - initState(minN, nVolFrac, 0); + initNVolatilizationState(minN, nVolFrac); double tEffect = calcTempEffect(climate->tsoil); double mEffect = calcMoistEffect(envi.soilWater, params.soilWHC); expNVolFlux = nVolFrac * minN * tEffect * mEffect; calcNVolatilizationFlux(); - status |= checkVolatilizationFlux(expNVolFlux); + //status |= checkVolatilizationFlux(expNVolFlux); + status |= checkFlux(fluxes.nVolatilization, expNVolFlux); // easy proportionality test - doubling params should double output minN *= 2; expNVolFlux *= 2; - initState(minN, nVolFrac, 0); + initNVolatilizationState(minN, nVolFrac); calcNVolatilizationFlux(); - status |= checkVolatilizationFlux(expNVolFlux); + //status |= checkVolatilizationFlux(expNVolFlux); + status |= checkFlux(fluxes.nVolatilization, expNVolFlux); // Check minN for the last updatePoolsForSoil(); @@ -137,7 +252,7 @@ int testNVolatilization(void) { int minStatus = 0; if (!compareDoubles(envi.minN, expMinN)) { minStatus = 1; - logTest("minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); + logTest("Volatilization minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); } status |= minStatus; @@ -151,7 +266,7 @@ int testFertilization(void) { double expMinN, expEventMinNFlux, expNVolFlux; // init minN 2, nVol 0.1 - initState(initN, nVolFrac, 0); + initNVolatilizationState(initN, nVolFrac); // fert event: 15 5 10 double fertMinN = 10; @@ -201,6 +316,8 @@ int run(void) { status |= testNLeaching(); + status |= testNFixation(); + return status; } diff --git a/tests/smoke/russell_2/sipnet.param b/tests/smoke/russell_2/sipnet.param index 0c84ae1fc..75387038b 100644 --- a/tests/smoke/russell_2/sipnet.param +++ b/tests/smoke/russell_2/sipnet.param @@ -80,3 +80,4 @@ microbePulseEff 0.45 mineralNInit 1.0 nVolatilizationFrac 0.1 nLeachingFrac 0.25 +nFixationFrac 0.11 \ No newline at end of file From da1c528c1e2bab419111b0dcc4fde3ead9bf2b24 Mon Sep 17 00:00:00 2001 From: mswilburn Date: Mon, 24 Nov 2025 21:11:15 -0600 Subject: [PATCH 2/2] Adding N gain from fixation --- src/sipnet/sipnet.c | 15 +++++++------ .../sipnet/test_modeling/testNitrogenCycle.c | 22 ++++++++++--------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/sipnet/sipnet.c b/src/sipnet/sipnet.c index 3358a5079..9a7d0d8f1 100644 --- a/src/sipnet/sipnet.c +++ b/src/sipnet/sipnet.c @@ -393,11 +393,11 @@ void readParamData(ModelParams **modelParamsPtr, const char *paramFile) { initializeOneModelParam(modelParams, "nFixationFrac", &(params.nFixationFrac), ctx.nitrogenCycle); // NOLINTEND - // clang-format on + // clang-format on - readModelParams(modelParams, paramF); + readModelParams(modelParams, paramF); - fclose(paramF); + fclose(paramF); } /*! @@ -411,9 +411,9 @@ void outputHeader(FILE *out) { fprintf(out, "year day time plantWoodC plantLeafC woodCreation "); fprintf(out, "soil microbeC coarseRootC fineRootC "); fprintf(out, "litter soilWater soilWetnessFrac snow "); - fprintf(out, - "npp nee cumNEE gpp rAboveground rSoil rRoot ra rh rtot " - "evapotranspiration fluxestranspiration minN n2oFlux nLeachFlux nFixationFlux\n"); + fprintf(out, "npp nee cumNEE gpp rAboveground rSoil rRoot ra rh rtot " + "evapotranspiration fluxestranspiration minN n2oFlux nLeachFlux " + "nFixationFlux\n"); } /*! @@ -1610,7 +1610,8 @@ void updatePoolsForSoil(void) { climate->length; // Nitrogen Cycle - envi.minN += (fluxes.nFixation - fluxes.nVolatilization - fluxes.nLeaching) * climate->length; + envi.minN += (fluxes.nFixation - fluxes.nVolatilization - fluxes.nLeaching) * + climate->length; } // !!! main runner function !!! diff --git a/tests/sipnet/test_modeling/testNitrogenCycle.c b/tests/sipnet/test_modeling/testNitrogenCycle.c index e94371b1f..2f4c12df4 100644 --- a/tests/sipnet/test_modeling/testNitrogenCycle.c +++ b/tests/sipnet/test_modeling/testNitrogenCycle.c @@ -111,8 +111,7 @@ int checkFlux(double calcFlux, double expFlux) { int status = 0; if (!compareDoubles(calcFlux, expFlux)) { status = 1; - logTest("Calculated N flux is %f, expected %f\n", - calcFlux, expFlux); + logTest("Calculated N flux is %f, expected %f\n", calcFlux, expFlux); } return status; @@ -146,7 +145,7 @@ int testNFixation(void) { double minN; double nFixFrac; double expNFixation; - + minN = 1; nFixFrac = 0.25; initNFixationState(minN, nFixFrac); @@ -167,7 +166,8 @@ int testNFixation(void) { int minStatus = 0; if (!compareDoubles(envi.minN, expMinN)) { minStatus = 1; - logTest("Fixation minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); + logTest("Fixation minN pool is %8.3f, expected %8.3f\n", envi.minN, + expMinN); } status |= minStatus; @@ -192,7 +192,7 @@ int testNLeaching(void) { } expNLeaching = minN * phi * nLeachFrac; calcNLeachingFlux(); - //status |= checkNLeachingFlux(expNLeaching); + // status |= checkNLeachingFlux(expNLeaching); status |= checkFlux(fluxes.nLeaching, expNLeaching); minN = 1; @@ -206,7 +206,7 @@ int testNLeaching(void) { } expNLeaching = minN * phi * nLeachFrac; calcNLeachingFlux(); - //status |= checkNLeachingFlux(expNLeaching); + // status |= checkNLeachingFlux(expNLeaching); status |= checkFlux(fluxes.nLeaching, expNLeaching); // Check minN for the last @@ -215,7 +215,8 @@ int testNLeaching(void) { int minStatus = 0; if (!compareDoubles(envi.minN, expMinN)) { minStatus = 1; - logTest("Leaching minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); + logTest("Leaching minN pool is %8.3f, expected %8.3f\n", envi.minN, + expMinN); } status |= minStatus; @@ -235,7 +236,7 @@ int testNVolatilization(void) { double mEffect = calcMoistEffect(envi.soilWater, params.soilWHC); expNVolFlux = nVolFrac * minN * tEffect * mEffect; calcNVolatilizationFlux(); - //status |= checkVolatilizationFlux(expNVolFlux); + // status |= checkVolatilizationFlux(expNVolFlux); status |= checkFlux(fluxes.nVolatilization, expNVolFlux); // easy proportionality test - doubling params should double output @@ -243,7 +244,7 @@ int testNVolatilization(void) { expNVolFlux *= 2; initNVolatilizationState(minN, nVolFrac); calcNVolatilizationFlux(); - //status |= checkVolatilizationFlux(expNVolFlux); + // status |= checkVolatilizationFlux(expNVolFlux); status |= checkFlux(fluxes.nVolatilization, expNVolFlux); // Check minN for the last @@ -252,7 +253,8 @@ int testNVolatilization(void) { int minStatus = 0; if (!compareDoubles(envi.minN, expMinN)) { minStatus = 1; - logTest("Volatilization minN pool is %8.3f, expected %8.3f\n", envi.minN, expMinN); + logTest("Volatilization minN pool is %8.3f, expected %8.3f\n", envi.minN, + expMinN); } status |= minStatus;