From 52beabda49a8eb1a8f683896cfc737c1698b9e40 Mon Sep 17 00:00:00 2001 From: ktrizzo <95317969+ktrizzo@users.noreply.github.com> Date: Tue, 10 Nov 2026 21:16:57 -0800 Subject: [PATCH] Added new stomatal model from Buckley et al 2026 --- plugins/photosynthesis/doc/Photosynthesis.dox | 2 + .../include/PhotosynthesisModel.h | 17 +- .../src/PhotosynthesisModel.cpp | 85 +++-- plugins/photosynthesis/tests/selfTest.cpp | 74 +++++ .../planthydraulics/doc/PlantHydraulics.dox | 6 +- .../include/PlantHydraulicsModel.h | 20 ++ .../src/PlantHydraulicsModel.cpp | 52 ++++ plugins/planthydraulics/tests/selfTest.cpp | 79 +++++ .../doc/StomatalConductance.dox | 163 +++++++++- .../include/StomatalConductanceModel.h | 98 ++++++ .../src/StomatalConductanceModel.cpp | 292 ++++++++++++++++++ .../stomatalconductance/tests/selfTest.cpp | 237 ++++++++++++++ 12 files changed, 1098 insertions(+), 27 deletions(-) diff --git a/plugins/photosynthesis/doc/Photosynthesis.dox b/plugins/photosynthesis/doc/Photosynthesis.dox index da9cafd30..3719180b9 100755 --- a/plugins/photosynthesis/doc/Photosynthesis.dox +++ b/plugins/photosynthesis/doc/Photosynthesis.dox @@ -57,6 +57,8 @@ Gamma\_CO2\f$\Gamma\f$\f$\mu\f$mol CO2/mol\htmlonlyfloat\endhtmlonlyFarquhar (C3)Photosynthetic CO2 compensation point (including "dark respiration"). limitation\_stateN/AN/A\htmlonlyint\endhtmlonlyFarquhar (C3), C4Photosynthesis limitation state. For the C3 Farquhar model: 0 = Rubisco-limited, 1 = electron-transport-limited. For the C4 model: 1 = enzyme (\f$A_c\f$) limited, 2 = electron-transport (\f$A_j\f$) limited. electron\_transport\_ratio\f$J/J_{max}\f$unitless\htmlonlyfloat\endhtmlonlyFarquhar (C3)Ratio of actual electron transport rate to maximum electron transport rate. Ranges from 0 (dark) to ~1 (light-saturated). Useful for computing fluorescence quantum yield. + electron\_transport\_rate\f$J\f$\f$\mu\f$mol/m2-sec\htmlonlyfloat\endhtmlonlyFarquhar (C3)Potential electron transport rate at the leaf's light and temperature, the smaller root of \f$\theta J^2-(\alpha Q+J_{max})J+\alpha QJ_{max}=0\f$. Does not depend on \f$C_i\f$. Written automatically when the Buckley et al. (2026) stomatal conductance model is in use. + dark\_respiration\f$R_d\f$\f$\mu\f$mol CO2/m2-sec\htmlonlyfloat\endhtmlonlyFarquhar (C3)Day respiration at the leaf temperature. Written automatically when the Buckley et al. (2026) stomatal conductance model is in use. Cm\f$C_m\f$\f$\mu\f$bar\htmlonlyfloat\endhtmlonlyC4Mesophyll cytosolic CO2 partial pressure, \f$C_m = C_i - A/g_m\f$. Vp\f$V_p\f$\f$\mu\f$mol CO2/m2/s\htmlonlyfloat\endhtmlonlyC4PEP carboxylation rate. diff --git a/plugins/photosynthesis/include/PhotosynthesisModel.h b/plugins/photosynthesis/include/PhotosynthesisModel.h index 54fbd681a..31474b544 100644 --- a/plugins/photosynthesis/include/PhotosynthesisModel.h +++ b/plugins/photosynthesis/include/PhotosynthesisModel.h @@ -911,7 +911,11 @@ class PhotosynthesisModel { //! Add optional output primitive data values to the Context /** * \param[in] label Name of primitive data. Available labels for the Farquhar model: "Ci" (intercellular CO2), "limitation_state" (0=Rubisco-limited, 1=electron transport-limited), "Gamma_CO2" (CO2 compensation point), "electron_transport_ratio" - * (relative light saturation Ja/Je of van der Tol et al. 2014: electron transport used by carbon metabolism divided by its light-limited potential; 1 in the dark and falling as light saturates photosynthesis; drives the RadiationModel SIF fluorescence yield). Additional labels for the C4 model: "Ci", "limitation_state" (1=enzyme-limited, 2=electron-transport-limited), "Cm" (mesophyll cytosolic CO2), "Vp" (PEP carboxylation rate), "electron_transport_ratio". + * (relative light saturation Ja/Je of van der Tol et al. 2014: electron transport used by carbon metabolism divided by its light-limited potential; 1 in the dark and falling as light saturates photosynthesis; drives the RadiationModel SIF + * fluorescence yield), "electron_transport_rate" (potential electron transport rate J at the leaf's light and temperature, umol/m^2/s), "dark_respiration" (day respiration Rd at the leaf temperature, umol/m^2/s). Additional labels for the C4 + * model: "Ci", "limitation_state" (1=enzyme-limited, 2=electron-transport-limited), "Cm" (mesophyll cytosolic CO2), "Vp" (PEP carboxylation rate), "electron_transport_ratio". + * \note "electron_transport_rate" and "dark_respiration" are also written automatically for the Farquhar model when the StomatalConductanceModel Buckley et al. (2026) model is in use (global data "buckley2026_stomatal_model_enabled" is 1), since + * that model needs them. */ void optionalOutputPrimitiveData(const char *label); @@ -992,6 +996,17 @@ class PhotosynthesisModel { static float respondToTemperature(const PhotosyntheticTemperatureResponseParameters *temperatureResponseParameters, float T); + //! Potential electron transport rate and day respiration of the Farquhar model at a given light and leaf temperature + /** + * Neither depends on the intercellular CO2 concentration, so they are evaluated directly rather than by the Ci solver. + * \param[in] params Farquhar model coefficients + * \param[in] i_PAR Incident PAR photon flux density (umol/m^2/s) + * \param[in] TL Leaf temperature (K) + * \param[out] J Potential electron transport rate (umol/m^2/s): the smaller root of theta*J^2 - (alpha*Q + Jmax)*J + alpha*Q*Jmax = 0 + * \param[out] Rd Day respiration (umol/m^2/s) + */ + static void evaluateElectronTransportAndRespiration(const FarquharModelCoefficients ¶ms, float i_PAR, float TL, double &J, float &Rd); + float i_PAR_default; float TL_default; float CO2_default; diff --git a/plugins/photosynthesis/src/PhotosynthesisModel.cpp b/plugins/photosynthesis/src/PhotosynthesisModel.cpp index a97d926b2..0bb340f2f 100644 --- a/plugins/photosynthesis/src/PhotosynthesisModel.cpp +++ b/plugins/photosynthesis/src/PhotosynthesisModel.cpp @@ -949,6 +949,17 @@ void PhotosynthesisModel::run(const std::vector &lUUIDs) { WarningAggregator warnings; warnings.setEnabled(message_flag); // Respect existing message flag + // The electron transport rate and day respiration are written when requested as optional output, and also whenever the StomatalConductanceModel Buckley et al. (2026) model is in use, + // which signals this through global data and computes the CO2-saturated electron-transport-limited rate J/4 - Rd from them. + bool buckley_inputs_requested = false; + if (context->doesGlobalDataExist("buckley2026_stomatal_model_enabled") && context->getGlobalDataType("buckley2026_stomatal_model_enabled") == HELIOS_TYPE_INT) { + int buckley_enabled = 0; + context->getGlobalData("buckley2026_stomatal_model_enabled", buckley_enabled); + buckley_inputs_requested = (buckley_enabled == 1); + } + const bool write_electron_transport_rate = buckley_inputs_requested || std::find(output_prim_data.begin(), output_prim_data.end(), "electron_transport_rate") != output_prim_data.end(); + const bool write_dark_respiration = buckley_inputs_requested || std::find(output_prim_data.begin(), output_prim_data.end(), "dark_respiration") != output_prim_data.end(); + for (uint UUID: lUUIDs) { float i_PAR; @@ -1164,6 +1175,18 @@ void PhotosynthesisModel::run(const std::vector &lUUIDs) { context->setPrimitiveData(UUID, "Vp", Vp_c4); } } + + if (model == "farquhar" && (write_electron_transport_rate || write_dark_respiration)) { + double J; + float Rd; + evaluateElectronTransportAndRespiration(getCoefficientsForPrimitive_Farquhar(UUID), i_PAR, TL, J, Rd); + if (write_electron_transport_rate) { + context->setPrimitiveData(UUID, "electron_transport_rate", static_cast(J)); + } + if (write_dark_respiration) { + context->setPrimitiveData(UUID, "dark_respiration", Rd); + } + } } warnings.report(std::cerr); @@ -1293,7 +1316,7 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector &var const FarquharModelCoefficients modelcoeffs = *reinterpret_cast(parameters); - float Vcmax, Jmax, TPU, Rd, alpha, theta; + float Vcmax, TPU, Rd, alpha; float CO2 = variables[0]; float Q = variables[1]; @@ -1318,17 +1341,9 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector &var Vcmax = respondToTemperature(&modelcoeffs.VcmaxTempResponse, TL); } - if (modelcoeffs.Jmax > 0) { - Jmax = modelcoeffs.Jmax * expf(modelcoeffs.dH_Jmax * (invDiffRT)); - } else { - Jmax = respondToTemperature(&modelcoeffs.JmaxTempResponse, TL); - } - - if (modelcoeffs.Rd > 0) { - Rd = modelcoeffs.Rd * expf(modelcoeffs.dH_Rd * (invDiffRT)); - } else { - Rd = respondToTemperature(&modelcoeffs.RdTempResponse, TL); - } + // Electron transport rate J and day respiration Rd do not depend on Ci + double J; + evaluateElectronTransportAndRespiration(modelcoeffs, Q, TL, J, Rd); if (modelcoeffs.alpha > 0) { alpha = modelcoeffs.alpha; @@ -1338,7 +1353,6 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector &var TPU = respondToTemperature(&modelcoeffs.TPUTempResponse, TL); - theta = respondToTemperature(&modelcoeffs.thetaTempResponse, TL); // Mesophyll conductance gm (mol CO2 / m^2 / s / bar). If gm is infinite/huge, // Cc ≡ Ci and we take the original legacy path below. Otherwise we solve a @@ -1353,14 +1367,6 @@ float PhotosynthesisModel::evaluateCi_Farquhar(float Ci, std::vector &var float Ko = exp(c_Ko - dH_Ko / (R * TL)); float Kco = Kc * (1.f + O / Ko); - // thetaJ^2 + -(alphaQ+Jmax)J + alphaQJmax = 0 - double a = std::max(theta, 0.0001f); - double ia = 1.000f / a; - double b = -(alpha * Q + Jmax); - double c = alpha * Q * Jmax; - double J = (-b - sqrt(pow(b, 2.000f) - 4.f * a * c)) * 0.5f * ia; - // J = Jmax * alpha * Q / (alpha * Q + Jmax); - float A; float limitation_state; @@ -1693,6 +1699,40 @@ float PhotosynthesisModel::evaluateC4Model(const C4ModelCoefficients ¶ms, fl } +void PhotosynthesisModel::evaluateElectronTransportAndRespiration(const FarquharModelCoefficients ¶ms, float i_PAR, float TL, double &J, float &Rd) { + + const float R = 0.0083144598; // molar gas constant (kJ/K/mol) + const float invDiffRT = (1.f / 298.15f - 1.f / TL) / R; + + float Jmax; + if (params.Jmax > 0) { + Jmax = params.Jmax * expf(params.dH_Jmax * (invDiffRT)); + } else { + Jmax = respondToTemperature(¶ms.JmaxTempResponse, TL); + } + + if (params.Rd > 0) { + Rd = params.Rd * expf(params.dH_Rd * (invDiffRT)); + } else { + Rd = respondToTemperature(¶ms.RdTempResponse, TL); + } + + float alpha; + if (params.alpha > 0) { + alpha = params.alpha; + } else { + alpha = respondToTemperature(¶ms.alphaTempResponse, TL); + } + const float theta = respondToTemperature(¶ms.thetaTempResponse, TL); + + // thetaJ^2 + -(alphaQ+Jmax)J + alphaQJmax = 0 + double a = std::max(theta, 0.0001f); + double ia = 1.000f / a; + double b = -(alpha * i_PAR + Jmax); + double c = alpha * i_PAR * Jmax; + J = (-b - sqrt(pow(b, 2.000f) - 4.f * a * c)) * 0.5f * ia; +} + float PhotosynthesisModel::respondToTemperature(const PhotosyntheticTemperatureResponseParameters *params, float T_in_Kelvin) { float T = T_in_Kelvin; float R = 0.0083144598f; @@ -1779,7 +1819,8 @@ void PhotosynthesisModel::enableMessages() { void PhotosynthesisModel::optionalOutputPrimitiveData(const char *label) { - if (strcmp(label, "Ci") == 0 || strcmp(label, "limitation_state") == 0 || strcmp(label, "Gamma_CO2") == 0 || strcmp(label, "electron_transport_ratio") == 0 || strcmp(label, "Cm") == 0 || strcmp(label, "Vp") == 0) { + if (strcmp(label, "Ci") == 0 || strcmp(label, "limitation_state") == 0 || strcmp(label, "Gamma_CO2") == 0 || strcmp(label, "electron_transport_ratio") == 0 || strcmp(label, "Cm") == 0 || strcmp(label, "Vp") == 0 || + strcmp(label, "electron_transport_rate") == 0 || strcmp(label, "dark_respiration") == 0) { output_prim_data.emplace_back(label); } else { if (message_flag) { diff --git a/plugins/photosynthesis/tests/selfTest.cpp b/plugins/photosynthesis/tests/selfTest.cpp index 08d3ec2d6..b379d3608 100644 --- a/plugins/photosynthesis/tests/selfTest.cpp +++ b/plugins/photosynthesis/tests/selfTest.cpp @@ -2357,6 +2357,80 @@ DOCTEST_TEST_CASE("PhotosynthesisModel Farquhar library values match documentati DOCTEST_CHECK(pist_male.getVcmaxTempResponse().value_at_25C != doctest::Approx(pist_female.getVcmaxTempResponse().value_at_25C).epsilon(1e-4)); } +DOCTEST_TEST_CASE("PhotosynthesisModel - electron_transport_rate and dark_respiration optional outputs") { + Context context; + PhotosynthesisModel photosynthesis(&context); + photosynthesis.disableMessages(); + + FarquharModelCoefficients coeffs; + coeffs.setVcmax(100.f); + coeffs.setJmax(150.f); + coeffs.setRd(1.5f); + coeffs.setQuantumEfficiency_alpha(0.3f); + coeffs.setLightResponseCurvature_theta(0.7f); + photosynthesis.setModelCoefficients(coeffs); + + const float PAR = 200.f; // W/m^2 + uint leaf = context.addPatch(); + context.setPrimitiveData(leaf, "radiation_flux_PAR", PAR); + context.setPrimitiveData(leaf, "temperature", 298.15f); + context.setPrimitiveData(leaf, "moisture_conductance", 0.3f); + + // Not written unless requested + photosynthesis.run(); + DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaf, "electron_transport_rate")); + DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaf, "dark_respiration")); + + photosynthesis.optionalOutputPrimitiveData("electron_transport_rate"); + photosynthesis.optionalOutputPrimitiveData("dark_respiration"); + photosynthesis.run(); + + // J is the smaller root of theta J^2 - (alpha Q + Jmax) J + alpha Q Jmax = 0, with the temperature-independent coefficients above + const double aQ = 0.3 * PAR * 4.57, Jmax = 150.0, theta = 0.7; + const double J_expected = ((aQ + Jmax) - std::sqrt((aQ + Jmax) * (aQ + Jmax) - 4.0 * theta * aQ * Jmax)) / (2.0 * theta); + float J, Rd; + context.getPrimitiveData(leaf, "electron_transport_rate", J); + context.getPrimitiveData(leaf, "dark_respiration", Rd); + DOCTEST_CHECK(J == doctest::Approx(J_expected).epsilon(1e-4)); + DOCTEST_CHECK(Rd == doctest::Approx(1.5f).epsilon(1e-6)); + + // J and Rd do not depend on Ci, so J/4 - Rd bounds the electron-transport-limited net assimilation from above + float A; + context.getPrimitiveData(leaf, "net_photosynthesis", A); + DOCTEST_CHECK(A < 0.25f * J - Rd); +} + +DOCTEST_TEST_CASE("PhotosynthesisModel - electron_transport_rate and dark_respiration are written when the Buckley et al. (2026) stomatal model requests them") { + Context context; + PhotosynthesisModel photosynthesis(&context); + photosynthesis.disableMessages(); + photosynthesis.setModelCoefficients(photosynthesis.getFarquharCoefficientsFromLibrary("Almond")); + + uint leaf = context.addPatch(); + context.setPrimitiveData(leaf, "radiation_flux_PAR", 300.f); + context.setPrimitiveData(leaf, "temperature", 303.15f); + + context.setGlobalData("buckley2026_stomatal_model_enabled", 0); + photosynthesis.run(); + DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaf, "electron_transport_rate")); + + context.setGlobalData("buckley2026_stomatal_model_enabled", 1); + photosynthesis.run(); + float J, Rd; + DOCTEST_REQUIRE(context.doesPrimitiveDataExist(leaf, "electron_transport_rate")); + DOCTEST_REQUIRE(context.doesPrimitiveDataExist(leaf, "dark_respiration")); + context.getPrimitiveData(leaf, "electron_transport_rate", J); + context.getPrimitiveData(leaf, "dark_respiration", Rd); + DOCTEST_CHECK(J > 0.f); + DOCTEST_CHECK(Rd > 0.2f); // almond Rd25 = 0.2, higher at 30 C + + // In the dark J is zero + context.setPrimitiveData(leaf, "radiation_flux_PAR", 0.f); + photosynthesis.run(); + context.getPrimitiveData(leaf, "electron_transport_rate", J); + DOCTEST_CHECK(J == doctest::Approx(0.f).epsilon(1e-6)); +} + int PhotosynthesisModel::selfTest(int argc, char **argv) { return helios::runDoctestWithValidation(argc, argv); } diff --git a/plugins/planthydraulics/doc/PlantHydraulics.dox b/plugins/planthydraulics/doc/PlantHydraulics.dox index c298aa4b0..88a53ee22 100644 --- a/plugins/planthydraulics/doc/PlantHydraulics.dox +++ b/plugins/planthydraulics/doc/PlantHydraulics.dox @@ -41,8 +41,10 @@ | Primitive Data Label | Symbol | Units | Data Type | Description | | -------------------- | ------ | ------ | --------- | ----------- | -| hydraulic_conductance | \f$K\f$ | mol m\f$^{-2}\f$ s\f$^{-1}\f$ MPa\f$^{-1}\f$ | \htmlonlyfloat\endhtmlonly | Hydraulic conductance | -| hydraulic_capacitance | \f$C\f$ | mol m\f$^{-2}\f$ MPa\f$^{-1}\f$ | \htmlonlyfloat\endhtmlonly | Hydraulic capacitance | +| hydraulic_conductance | \f$K_{leaf}\f$ | mol m\f$^{-2}\f$ s\f$^{-1}\f$ MPa\f$^{-1}\f$ | \htmlonlyfloat\endhtmlonly | Leaf hydraulic conductance. Enabled with PlantHydraulicsModel::outputConductancePrimitiveData(). | +| hydraulic_capacitance | \f$C\f$ | mol m\f$^{-2}\f$ MPa\f$^{-1}\f$ | \htmlonlyfloat\endhtmlonly | Hydraulic capacitance. Enabled with PlantHydraulicsModel::outputCapacitancePrimitiveData(). | +| plant_hydraulic_conductance | \f$K_{plant}\f$ | mol m\f$^{-2}\f$ s\f$^{-1}\f$ MPa\f$^{-1}\f$ | \htmlonlyfloat\endhtmlonly | Whole-plant (soil-to-leaf) hydraulic conductance per unit leaf area, \f$K_{plant}=\left(K_{root}^{-1}+K_{stem}^{-1}+K_{leaf}^{-1}\right)^{-1}\f$. Enabled with PlantHydraulicsModel::outputPlantConductancePrimitiveData(), or automatically when the Buckley et al. (2026) stomatal conductance model is in use. | +| soil_water_potential | \f$\psi_{soil}\f$ | MPa | \htmlonlyfloat\endhtmlonly | Soil water potential of the leaf's plant. Written together with plant_hydraulic_conductance, but only for plants whose soil water potential was set with PlantHydraulicsModel::setSoilWaterPotentialOfPlant(), so that a value set directly as primitive data is not replaced by the default placeholder. | \section PH_Intro Introduction The Plant Hydraulics plugin models the water potentials, water contents, hydraulic conductances, and diff --git a/plugins/planthydraulics/include/PlantHydraulicsModel.h b/plugins/planthydraulics/include/PlantHydraulicsModel.h index 7ef552dd6..a669ca88a 100644 --- a/plugins/planthydraulics/include/PlantHydraulicsModel.h +++ b/plugins/planthydraulics/include/PlantHydraulicsModel.h @@ -20,6 +20,7 @@ SPDX-License-Identifier: LGPL-2.1-or-later #define PLANTHYDRAULICS_MODEL #include "Context.h" +#include struct HydraulicConductance { HydraulicConductance() { @@ -611,6 +612,16 @@ class PlantHydraulicsModel { * */ void outputCapacitancePrimitiveData(bool toggle); + //! Toggle the writing out of optional primitive data - whole-plant hydraulic conductance and soil water potential + /** + * Writes "plant_hydraulic_conductance" (\f$K_{plant}\f$, the soil-to-leaf conductance of the root, stem, and leaf conductances in series, per unit leaf area, mol/m^2/s/MPa) and "soil_water_potential" (MPa, of the + * leaf's plant) to each leaf primitive. This differs from the optional output "hydraulic_conductance", which is the leaf conductance \f$K_{leaf}\f$ only. "soil_water_potential" is written only for plants whose + * soil water potential was set with setSoilWaterPotentialOfPlant(), so that a value set directly as primitive data is not replaced by the default placeholder. + * \param[in] toggle + * \note Both are also written automatically when the StomatalConductanceModel Buckley et al. (2026) model is in use (global data "buckley2026_stomatal_model_enabled" is 1), since that model uses them. + * */ + void outputPlantConductancePrimitiveData(bool toggle); + //! Group leaf primitives into a plant that share stem, root, and soil hydraulics /** * \param[in] UUIDs Leaf primitives @@ -736,11 +747,20 @@ class PlantHydraulicsModel { std::map stemWaterPotentialsByPlantID; std::map rootWaterPotentialsByPlantID; std::map soilWaterPotentialsByPlantID; + //! Plants whose soil water potential was set with setSoilWaterPotentialOfPlant() (others use a placeholder value, which is not written as "soil_water_potential" primitive data) + std::set soilWaterPotentialSetByUser; std::vector plantIDs; bool message_flag = false; bool OutputConductance = false; bool OutputCapacitance = false; + bool OutputPlantConductance = false; + + //! True if "plant_hydraulic_conductance" and "soil_water_potential" should be written: requested with outputPlantConductancePrimitiveData(), or the Buckley et al. (2026) stomatal model is in use + bool plantConductanceOutputRequested() const; + + //! Whole-plant (soil-to-leaf) hydraulic conductance K_plant per unit leaf area: the root, stem, and leaf conductances in series (mol/m^2/s/MPa); 0 if any of them is not positive + static float computePlantConductance(float K_root, float K_stem, float K_leaf); void updateRootAndStemWaterPotentialsOfPlant(const std::vector &UUIDs, int timespan, int timestep); diff --git a/plugins/planthydraulics/src/PlantHydraulicsModel.cpp b/plugins/planthydraulics/src/PlantHydraulicsModel.cpp index 0d17427ae..9d85b7d4a 100644 --- a/plugins/planthydraulics/src/PlantHydraulicsModel.cpp +++ b/plugins/planthydraulics/src/PlantHydraulicsModel.cpp @@ -348,6 +348,21 @@ void PlantHydraulicsModel::updateLeafWaterPotentialsOfPlant(const std::vector 0; + } + if (steadystate) { for (uint i = 0; i < UUIDs.size(); i++) { auto it = modelcoeffs_map.find(UUIDs.at(i)); @@ -370,6 +385,12 @@ void PlantHydraulicsModel::updateLeafWaterPotentialsOfPlant(const std::vectorsetPrimitiveData(UUIDs.at(i), "turgor_pressure", leaf_turgor_pressure); context->setPrimitiveData(UUIDs.at(i), "osmotic_potential", leaf_osmotic_potential); context->setPrimitiveData(UUIDs.at(i), "relative_water_content", leaf_relative_water_content); + if (write_plant_conductance) { + context->setPrimitiveData(UUIDs.at(i), "plant_hydraulic_conductance", computePlantConductance(K_root_plant, K_stem_plant, K_leaf)); + if (write_soil_water_potential) { + context->setPrimitiveData(UUIDs.at(i), "soil_water_potential", soil_water_potential_plant); + } + } } } else { // Non-steady-state @@ -420,6 +441,13 @@ void PlantHydraulicsModel::updateLeafWaterPotentialsOfPlant(const std::vectorsetPrimitiveData(UUIDs.at(i), "hydraulic_capacitance", C_leaf); } + if (write_plant_conductance) { + K_leaf = computeConductance(coeffs.LeafHydraulicConductance, leaf_water_potential, leaf_temperature); // at the updated leaf water potential + context->setPrimitiveData(UUIDs.at(i), "plant_hydraulic_conductance", computePlantConductance(K_root_plant, K_stem_plant, K_leaf)); + if (write_soil_water_potential) { + context->setPrimitiveData(UUIDs.at(i), "soil_water_potential", soil_water_potential_plant); + } + } } } } @@ -428,6 +456,29 @@ void PlantHydraulicsModel::outputConductancePrimitiveData(bool toggle) { OutputConductance = toggle; } +void PlantHydraulicsModel::outputPlantConductancePrimitiveData(bool toggle) { + OutputPlantConductance = toggle; +} + +bool PlantHydraulicsModel::plantConductanceOutputRequested() const { + if (OutputPlantConductance) { + return true; + } + if (context->doesGlobalDataExist("buckley2026_stomatal_model_enabled") && context->getGlobalDataType("buckley2026_stomatal_model_enabled") == helios::HELIOS_TYPE_INT) { + int buckley_enabled = 0; + context->getGlobalData("buckley2026_stomatal_model_enabled", buckley_enabled); + return buckley_enabled == 1; + } + return false; +} + +float PlantHydraulicsModel::computePlantConductance(float K_root, float K_stem, float K_leaf) { + if (!(K_root > 0.f) || !(K_stem > 0.f) || !(K_leaf > 0.f)) { + return 0.f; + } + return 1.f / (1.f / K_root + 1.f / K_stem + 1.f / K_leaf); +} + void PlantHydraulicsModel::outputCapacitancePrimitiveData(bool toggle) { OutputCapacitance = toggle; } @@ -638,6 +689,7 @@ float PlantHydraulicsModel::getSoilWaterPotentialOfPlant(uint plantID) { void PlantHydraulicsModel::setSoilWaterPotentialOfPlant(uint plantID, float soil_water_potential) { soilWaterPotentialsByPlantID[plantID] = soil_water_potential; + soilWaterPotentialSetByUser.insert(plantID); } std::vector PlantHydraulicsModel::getPrimitivesWithoutPlantID(std::vector UUIDs) { diff --git a/plugins/planthydraulics/tests/selfTest.cpp b/plugins/planthydraulics/tests/selfTest.cpp index 9d957f336..f01eda002 100644 --- a/plugins/planthydraulics/tests/selfTest.cpp +++ b/plugins/planthydraulics/tests/selfTest.cpp @@ -621,6 +621,85 @@ DOCTEST_TEST_CASE("PlantHydraulicsModel - Documentation Example 2") { } } +DOCTEST_TEST_CASE("PlantHydraulicsModel - whole-plant hydraulic conductance K_plant and soil water potential outputs") { + // K_plant is the root, stem, and leaf conductances in series (per unit leaf area); it is written on request, or when the Buckley et al. (2026) stomatal model asks for it + const float K_root = 0.08f, K_stem = 0.05f, K_leaf = 0.02f; // mol/m^2/s/MPa, independent of water potential and temperature + const float K_plant_expected = 1.f / (1.f / K_root + 1.f / K_stem + 1.f / K_leaf); + const float psi_soil = -0.3f; + + auto setup = [&](Context &context, PlantHydraulicsModel &hydraulics) { + uint objID = context.addTileObject(nullorigin, make_vec2(1, 1), nullrotation, make_int2(2, 2)); + std::vector leaves = context.getObjectPrimitiveUUIDs(objID); + context.setObjectData(objID, "plantID", 1); + context.setPrimitiveData(leaves, "latent_flux", 100.f); + PlantHydraulicsModelCoefficients coeffs; + coeffs.setRootHydraulicConductance(K_root); + coeffs.setStemHydraulicConductance(K_stem); + coeffs.setLeafHydraulicConductance(K_leaf); + coeffs.setRootHydraulicCapacitance(1.f); // fixed capacitances, needed by the time-stepped solution + coeffs.setStemHydraulicCapacitance(1.f); + hydraulics.setModelCoefficients(coeffs); + hydraulics.setSoilWaterPotentialOfPlant(1, psi_soil); + return leaves; + }; + + SUBCASE("not written by default") { + Context context; + PlantHydraulicsModel hydraulics(&context); + std::vector leaves = setup(context, hydraulics); + hydraulics.run(leaves); + DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaves.front(), "plant_hydraulic_conductance")); + DOCTEST_CHECK(!context.doesPrimitiveDataExist(leaves.front(), "soil_water_potential")); + } + + SUBCASE("requested with outputPlantConductancePrimitiveData, steady state and time-stepped") { + for (int timespan: {0, 1000}) { + Context context; + PlantHydraulicsModel hydraulics(&context); + std::vector leaves = setup(context, hydraulics); + hydraulics.outputPlantConductancePrimitiveData(true); + hydraulics.run(leaves, timespan, 100); + for (uint UUID: leaves) { + float K_plant, psi; + DOCTEST_REQUIRE(context.doesPrimitiveDataExist(UUID, "plant_hydraulic_conductance")); + context.getPrimitiveData(UUID, "plant_hydraulic_conductance", K_plant); + context.getPrimitiveData(UUID, "soil_water_potential", psi); + DOCTEST_CHECK(K_plant == doctest::Approx(K_plant_expected).epsilon(1e-5)); + DOCTEST_CHECK(psi == doctest::Approx(psi_soil).epsilon(1e-6)); + } + } + } + + SUBCASE("a soil water potential never set with setSoilWaterPotentialOfPlant() does not replace one set directly") { + Context context; + PlantHydraulicsModel hydraulics(&context); + uint objID = context.addTileObject(nullorigin, make_vec2(1, 1), nullrotation, make_int2(2, 2)); + std::vector leaves = context.getObjectPrimitiveUUIDs(objID); + context.setObjectData(objID, "plantID", 2); + context.setPrimitiveData(leaves, "latent_flux", 100.f); + context.setPrimitiveData(leaves, "soil_water_potential", -1.5f); // set directly, e.g. from a soil model + hydraulics.outputPlantConductancePrimitiveData(true); + hydraulics.run(leaves); + float psi; + context.getPrimitiveData(leaves.front(), "soil_water_potential", psi); + DOCTEST_CHECK(psi == doctest::Approx(-1.5f).epsilon(1e-6)); + DOCTEST_CHECK(context.doesPrimitiveDataExist(leaves.front(), "plant_hydraulic_conductance")); // K_plant is still written + } + + SUBCASE("requested by the Buckley et al. (2026) stomatal model through global data") { + Context context; + PlantHydraulicsModel hydraulics(&context); + std::vector leaves = setup(context, hydraulics); + context.setGlobalData("buckley2026_stomatal_model_enabled", 1); + hydraulics.run(leaves); + float K_plant; + DOCTEST_REQUIRE(context.doesPrimitiveDataExist(leaves.front(), "plant_hydraulic_conductance")); + context.getPrimitiveData(leaves.front(), "plant_hydraulic_conductance", K_plant); + DOCTEST_CHECK(K_plant == doctest::Approx(K_plant_expected).epsilon(1e-5)); + DOCTEST_CHECK(K_plant < K_leaf); // K_plant is not the leaf conductance + } +} + int PlantHydraulicsModel::selfTest(int argc, char **argv) { return helios::runDoctestWithValidation(argc, argv); } diff --git a/plugins/stomatalconductance/doc/StomatalConductance.dox b/plugins/stomatalconductance/doc/StomatalConductance.dox index 88583af02..295584ccf 100755 --- a/plugins/stomatalconductance/doc/StomatalConductance.dox +++ b/plugins/stomatalconductance/doc/StomatalConductance.dox @@ -145,7 +145,7 @@ \f$\mu\f$mol/m2-s \htmlonlyfloat\endhtmlonly Net CO2 (photosynthetic) flux. - BWB, BBL, Mopt + BWB, BBL, Mopt, Buckley2026 Can be computed by \ref PhotosynthesisModel plug-in. 0 @@ -179,6 +179,46 @@ N/A -0.1 MPa + + electron_transport_rate + \f$J\f$ + \f$\mu\f$mol/m2-s + \htmlonlyfloat\endhtmlonly + Potential electron transport rate at the leaf's light and temperature. + Buckley2026 + Computed by the \ref PhotosynthesisModel plug-in (Farquhar model), automatically once the Buckley2026 model is selected. + None (required) + + + dark_respiration + \f$R_d\f$ + \f$\mu\f$mol/m2-s + \htmlonlyfloat\endhtmlonly + Day respiration at the leaf temperature. + Buckley2026 + Computed by the \ref PhotosynthesisModel plug-in (Farquhar model), automatically once the Buckley2026 model is selected. + None (required) + + + soil_water_potential + \f$\psi_{soil}\f$ + MPa + \htmlonlyfloat\endhtmlonly + Water potential of the soil (source) supplying the plant. + Buckley2026 + Written by the \ref PlantHydraulicsModel plug-in (for plants whose soil water potential was set with setSoilWaterPotentialOfPlant()), automatically once the Buckley2026 model is selected, or set directly. + -0.1 MPa + + + plant_hydraulic_conductance + \f$K_{plant}\f$ + mol/m2-s-MPa + \htmlonlyfloat\endhtmlonly + Whole-plant (soil-to-leaf) hydraulic conductance per unit leaf area, at the leaf temperature. Not the leaf conductance \f$K_{leaf}\f$. + Buckley2026 + Written by the \ref PlantHydraulicsModel plug-in, automatically once the Buckley2026 model is selected. + Model coefficient \f$K_{plant,25}\f$ scaled to leaf temperature + @@ -215,7 +255,7 @@ - +
Primitive Data LabelUnitsData TypeDescription
vapor\_pressure\_deficitmmol/mol\htmlonlyfloat\endhtmlonlyVapor pressure deficit between the stomatal cavity and leaf surface. Only applicable for BMF model.
vapor\_pressure\_deficitmmol/mol\htmlonlyfloat\endhtmlonlyVapor pressure deficit between the stomatal cavity and leaf surface. Only applicable for the BMF and Buckley2026 models.
model\_parametersN/A\htmlonlystd::string\endhtmlonlyModel parameters for each primitive. See table below for the actual primitive data labels that are created.
@@ -258,6 +298,18 @@ \f$b\f$b_BMF + + Buckley et al. (2026)\f$s\f$s_Buckley2026 + + + \f$K_{plant}\f$ (value used, at leaf temperature)K_plant_Buckley2026 + + + \f$b_o\f$b_o_Buckley2026 + + + \f$\psi_c\f$psi_c_Buckley2026 + \section StomatalIntro Introduction @@ -713,6 +765,105 @@ Because \f$g_s\f$ appears on both sides of these expressions but enters only lin As with the other models, \f$D_s\f$ is obtained by iteratively solving the water vapor flux balance at the leaf surface, so the Bailey model responds to the boundary-layer conductance. This model requires the xylem water potential (primitive data "xylem_water_potential") rather than the net photosynthetic flux. +\section Buckley2026Theory Buckley et al. (2026) Model + +The model of Buckley et al. (2026) (Buckley, T. N., J. Lamour, D. M. Barnard, D. W. Buckley, A. J. Jarvis, A. Diaz-Espejo, A. Rogers, and L. Sack. 2026. "A Simple Stomatal Model That Unifies the Metabolic and Hydraulic Control of Carbon and Water Flux." Global Change Biology 32: e70999) simplifies the hydromechanical/biochemical model of Buckley, Mott, and Farquhar (2003). Stomatal conductance is driven by the residual photosynthetic capacity, the excess of the CO2-saturated assimilation rate \f$A_m\f$ over the actual net assimilation rate, and limited by the water supply from the soil through the plant (their Equation 3): + +\f[g_s = \frac{\beta\, s\, b\,(\psi_{soil}-\psi_c)}{1+s\,b\,D_s/K_{plant}},\qquad b = \max\left(0,\ \frac{J}{4}-R_d-A+b_o\right),\f] + +where \f$J\f$ is the potential electron transport rate at the leaf's light and temperature, \f$R_d\f$ is day respiration, \f$A\f$ is net photosynthesis, \f$\psi_{soil}\f$ is the soil (source) water potential, \f$D_s\f$ is the vapor pressure deficit between the stomatal cavity and the leaf surface (mol/mol; calculated as described for the other models above), and \f$K_{plant}\f$ is the whole-plant (soil-to-leaf) hydraulic conductance per unit leaf area. \f$J/4-R_d\f$ is the limit of the Farquhar electron-transport-limited rate \f$J(C_c-\Gamma^*)/(4C_c+8\Gamma^*)-R_d\f$ as \f$C_c\to\infty\f$. The denominator accounts for the fall in leaf water potential with transpiration, \f$\psi_{leaf}=\psi_{soil}-g_sD_s/K_{plant}\f$. The parameters are \f$s\f$ (mol \f$\mu\f$mol-1 MPa-1), the closure water potential \f$\psi_c\f$ (MPa), \f$b_o\f$ (\f$\mu\f$mol/m2-s), and \f$K_{plant,25}\f$ (mol/m2-s-MPa), the whole-plant hydraulic conductance at 25 \f$^\circ\f$C. When primitive data "plant_hydraulic_conductance" is not available, \f$K_{plant}=K_{plant,25}\,\phi(T_L)\f$, where \f$\phi(T)=4.891\times10^{10}\exp(-4209/T-0.04527T+3.376\times10^{-5}T^2)\f$ is the fluidity of water relative to 25 \f$^\circ\f$C. Conductance is zero when \f$\psi_{soil}\le\psi_c\f$. + +Note that \f$A\f$, \f$J\f$, and \f$R_d\f$ must be computed by the \ref PhotosynthesisModel plug-in (Farquhar model); once this model is selected, the photosynthesis plug-in writes \f$J\f$ and \f$R_d\f$ automatically. \f$\psi_{soil}\f$ and \f$K_{plant}\f$ can be computed by the \ref PlantHydraulicsModel plug-in, or \f$\psi_{soil}\f$ can be set directly as primitive data, in which case \f$K_{plant}\f$ is given by the model coefficient \f$K_{plant,25}\f$. + +\subsection Buckley2026Library Species library + +The coefficients below were fitted to diurnal LI-600 porometer surveys of field and orchard trees (gsdb, Davis, CA; conductances corrected per Rizzo & Bailey) with soft-L1 least squares. \f$\psi_{soil}\f$ was the first (predawn or early-morning) stem water potential measured on each tree and day. \f$\psi_c\f$ was fixed at the pressure-volume turgor loss point \f$\psi_{TLP}\f$ of the species (measured, or from the literature where marked), as recommended by Buckley et al. (2026), and fitted for eastern redbud, which has no pressure-volume data. With a single \f$\psi_{soil}\f$ per day only \f$s(\psi_{soil}-\psi_c)\f$ and \f$K_{plant}(\psi_{soil}-\psi_c)\f$ are identified, so \f$s\f$ and \f$K_{plant,25}\f$ depend on the value used for \f$\psi_c\f$. \f$A\f$, \f$J\f$, and \f$R_d\f$ were computed with the \ref PhotosynthesisModel Farquhar library entry listed for each species, and the fit depends on those photosynthesis parameters through \f$b\f$: use each entry together with that Farquhar library entry. \f$K_{plant,25}\f$ was fitted, so it is the value consistent with \f$s\f$ and \f$b_o\f$. \f$R^2\f$ is for the means of each date and hour, with the leave-one-hour-out cross-validated value in parentheses. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Species\f$s\f$ (mol \f$\mu\f$mol-1 MPa-1)\f$K_{plant,25}\f$ (mol/m2-s-MPa)\f$b_o\f$ (\f$\mu\f$mol/m2-s)\f$\psi_c\f$ (MPa)Farquhar library entryReadings / days\f$R^2\f$ (CV)
Almond0.015950.032661.416-2.329 (turgor loss point, literature)Almond147 / 10.75 (0.56)
Elderberry0.02710.0074911.27-2.752 (turgor loss point)Elderberry389 / 20.83 (0.77)
Cowpea0.26260.057141.074-0.846 (turgor loss point)Cowpea258 / 10.76 (0.65)
EasternRedbud0.0041950.0017441.774-3.667 (fitted)EasternRedbudSunlit, EasternRedbudShaded2548 / 40.63 (0.57)
Grape0.033410.0074711.268-1.718 (turgor loss point, literature)Grape231 / 10.55 (0.16)
Olive0.021560.0045625.033-2.043 (turgor loss point)Olive57 / 10.88 (0.73)
Pistachio0.087660.0074111.11-4.103 (turgor loss point)PistachioFemale (or Pistachio)824 / 10.76 (0.60)
+ +The default \ref Buckley2026coefficients are the pistachio entry. The cowpea entry is a single genotype measured on one day. + \section StomatalErrors Error Conditions The plug-in follows a fail-fast philosophy: physically meaningless inputs produce an error rather than a silently substituted value. \ref StomatalConductanceModel::run() "run()" will throw an error if any of the following conditions is encountered. @@ -770,6 +921,14 @@ Each model has a data structure containing member variables for each of the mode Buckley, Mott, Farquhar \ref BMFcoefficients "BMFcoefficients" + + Bailey + \ref BBcoefficients "BBcoefficients" + + + Buckley et al. (2026) + \ref Buckley2026coefficients "Buckley2026coefficients" + In order to modify model coefficients from the default values, one should 1) declare and instance of the coefficient structure, 2) modify the appropriate data elements of the data structure, and 3) pass the data structure to the \ref StomatalConductanceModel::setModelCoefficients() function. Model coefficients can be set to the same value for all primitives, or differently for a subset of primitives based on a vector of UUIDs. Example code is given below diff --git a/plugins/stomatalconductance/include/StomatalConductanceModel.h b/plugins/stomatalconductance/include/StomatalConductanceModel.h index 79239e201..3dc4f46c1 100644 --- a/plugins/stomatalconductance/include/StomatalConductanceModel.h +++ b/plugins/stomatalconductance/include/StomatalConductanceModel.h @@ -94,6 +94,30 @@ struct BBcoefficients { float chi; // mol/m^2-s-MPa }; +//! Coefficients for the Buckley et al. (2026) stomatal conductance model +/** + * Buckley, T. N., J. Lamour, D. M. Barnard, D. W. Buckley, A. J. Jarvis, A. Diaz-Espejo, A. Rogers, and L. Sack. 2026. A Simple Stomatal Model That Unifies the Metabolic and Hydraulic Control of Carbon and + * Water Flux. Global Change Biology 32: e70999. https://doi.org/10.1111/gcb.70999 + * + * \f$g_s = \beta\,s\,b\,(\psi_{soil}-\psi_c)/(1+s\,b\,D_s/K_{plant})\f$ with \f$b = J/4 - R_d - A + b_o\f$ (floored at zero), where \f$J/4-R_d\f$ is the CO2-saturated electron-transport-limited + * photosynthetic rate. \f$J\f$, \f$R_d\f$, and \f$A\f$ come from the Farquhar photosynthesis model, and \f$K_{plant}\f$ is the whole-plant (soil-to-leaf) hydraulic conductance per unit leaf area. + */ +struct Buckley2026coefficients { + + Buckley2026coefficients() { + // Pistachio (see the species library) + s = 0.08766f; + K_plant_at_25C = 0.007411f; + b_o = 1.11f; + psi_c = -4.103f; + } + + float s; // mol/umol/MPa, stomatal scaling and sensitivity coefficient (mol m^-2 s^-1 MPa^-1 of stomatal conductance per umol m^-2 s^-1 of b) + float K_plant_at_25C; // mol/m^2/s/MPa, whole-plant (soil-to-leaf) hydraulic conductance per unit leaf area at 25 C; scaled by the fluidity of water at leaf temperature. Not the leaf conductance K_leaf. + float b_o; // umol/m^2/s + float psi_c; // MPa, leaf water potential at which stomata close (Buckley et al. 2026 recommend the turgor loss point when it is not measured) +}; + class StomatalConductanceModel { public: //! Default constructor @@ -170,6 +194,57 @@ class StomatalConductanceModel { */ void setModelCoefficients(const BBcoefficients &coeffs, const std::vector &UUIDs); + //! Set the model coefficient values for all primitives - Buckley et al. (2026) model + /** + * \param[in] coeffs Model coefficient values + * \note Selecting this model sets global data "buckley2026_stomatal_model_enabled" to 1, which makes the PhotosynthesisModel write "electron_transport_rate" and "dark_respiration", and the + * PlantHydraulicsModel write "plant_hydraulic_conductance" and "soil_water_potential", on their next run. The global data is not reset when another model is selected. + */ + void setModelCoefficients(const Buckley2026coefficients &coeffs); + + //! Set the model coefficient values for a subset of primitives based on their UUIDs - Buckley et al. (2026) model + /** + * \param[in] coeffs Model coefficient values + * \param[in] UUIDs Universal unique identifiers for primitives to set + */ + void setModelCoefficients(const Buckley2026coefficients &coeffs, const std::vector &UUIDs); + + //! Set Buckley et al. (2026) model coefficients for a material by label + /** + * \param[in] material_label String identifier for the material + * \param[in] coeffs Model coefficient values + */ + void setModelCoefficients(const std::string &material_label, const Buckley2026coefficients &coeffs); + + //! Set the Buckley et al. (2026) model coefficients for all primitives using one of the available species in the library + /** + * \param[in] species Name of species + * \note Each entry was fitted together with the PhotosynthesisModel Farquhar library entry of the same species, so use it with PhotosynthesisModel::setFarquharCoefficientsFromLibrary() for that species. + * \note Throws a runtime error if the species is not in the library; see the stomatal conductance plug-in documentation for the list of valid names. + */ + void setBuckley2026CoefficientsFromLibrary(const std::string &species); + + //! Set the Buckley et al. (2026) model coefficients for a subset of primitives using one of the available species in the library + /** + * \param[in] species Name of species + * \param[in] UUIDs Universal unique identifiers for primitives to set + */ + void setBuckley2026CoefficientsFromLibrary(const std::string &species, const std::vector &UUIDs); + + //! Set the Buckley et al. (2026) model coefficients for a material using one of the available species in the library + /** + * \param[in] species Name of species + * \param[in] material_label String identifier for the material + */ + void setBuckley2026CoefficientsFromLibrary(const std::string &species, const std::string &material_label); + + //! Get the Buckley et al. (2026) model coefficients of one of the available species in the library + /** + * \param[in] species Name of species + * \note Throws a runtime error if the species is not in the library; see the stomatal conductance plug-in documentation for the list of valid names. + */ + Buckley2026coefficients getBuckley2026CoefficientsFromLibrary(const std::string &species); + //! Set the model coefficient values for a subset of primitives based on their UUIDs /** * \param[in] coeffs Vecotr of sets of model coefficient values @@ -317,6 +392,7 @@ class StomatalConductanceModel { MOPTcoefficients MOPTcoeffs; BMFcoefficients BMFcoeffs; BBcoefficients BBcoeffs; + Buckley2026coefficients Buckley2026coeffs; float i_default; float TL_default; @@ -329,12 +405,14 @@ class StomatalConductanceModel { float air_CO2_default; float blconductance_default; float beta_default; + float soil_water_potential_default; std::map BWBmodel_coefficients; std::map BBLmodel_coefficients; std::map MOPTmodel_coefficients; std::map BMFmodel_coefficients; std::map BBmodel_coefficients; + std::map Buckley2026model_coefficients; // Cache to avoid repeated Context lookups during run() mutable std::map material_coefficient_cache_BWB; @@ -342,6 +420,7 @@ class StomatalConductanceModel { mutable std::map material_coefficient_cache_MOPT; mutable std::map material_coefficient_cache_BMF; mutable std::map material_coefficient_cache_BB; + mutable std::map material_coefficient_cache_Buckley2026; // Helper to retrieve coefficients with caching BWBcoefficients getCoefficientsForPrimitive_BWB(uint UUID) const; @@ -349,6 +428,10 @@ class StomatalConductanceModel { MOPTcoefficients getCoefficientsForPrimitive_MOPT(uint UUID) const; BMFcoefficients getCoefficientsForPrimitive_BMF(uint UUID) const; BBcoefficients getCoefficientsForPrimitive_BB(uint UUID) const; + Buckley2026coefficients getCoefficientsForPrimitive_Buckley2026(uint UUID) const; + + //! Select the Buckley et al. (2026) model and signal the photosynthesis and plant hydraulics plug-ins (global data "buckley2026_stomatal_model_enabled") to write the inputs it needs + void selectBuckley2026Model(); std::map dynamic_time_constants; // .x is tau_open, .y is tau_close @@ -357,6 +440,21 @@ class StomatalConductanceModel { static float evaluate_MOPTmodel(float esurf, std::vector &variables, const void *parameters); static float evaluate_BMFmodel(float esurf, std::vector &variables, const void *parameters); static float evaluate_BBmodel(float esurf, std::vector &variables, const void *parameters); + static float evaluate_Buckley2026model(float esurf, std::vector &variables, const void *parameters); + + //! Buckley et al. (2026) conductance for a given surface vapor pressure deficit + /** + * \param[in] Ds Vapor pressure deficit between the stomatal cavity and the leaf surface (mol/mol) + * \param[in] b \f$J/4-R_d-A+b_o\f$ (umol/m^2/s), floored at zero by the caller + * \param[in] psi_soil Soil (source) water potential (MPa) + * \param[in] K_plant Whole-plant hydraulic conductance per unit leaf area at the leaf temperature (mol/m^2/s/MPa) + * \param[in] beta Soil moisture factor + * \param[in] coeffs Model coefficients + */ + static float evaluate_Buckley2026conductance(float Ds, float b, float psi_soil, float K_plant, float beta, const Buckley2026coefficients &coeffs); + + //! Fluidity of liquid water at temperature T relative to 25 C (the temperature scaling of K_plant_at_25C) + static float waterFluidityRelativeTo25C(float T_in_Kelvin); //! Closed-form solution of the implicit Bailey guard-cell relation for a given surface vapor pressure deficit static float evaluate_BBconductance(float Ds, float Psix, float i, float beta, const BBcoefficients &coeffs); diff --git a/plugins/stomatalconductance/src/StomatalConductanceModel.cpp b/plugins/stomatalconductance/src/StomatalConductanceModel.cpp index d1d14788f..75b33cbbe 100644 --- a/plugins/stomatalconductance/src/StomatalConductanceModel.cpp +++ b/plugins/stomatalconductance/src/StomatalConductanceModel.cpp @@ -48,6 +48,15 @@ namespace { constexpr const char *LABEL_BB_sigma = "gs_bb_sigma"; constexpr const char *LABEL_BB_chi = "gs_bb_chi"; + // Material data labels for Buckley et al. (2026) model + constexpr const char *LABEL_Buckley2026_s = "gs_buckley2026_s"; + constexpr const char *LABEL_Buckley2026_K_plant_at_25C = "gs_buckley2026_K_plant_at_25C"; + constexpr const char *LABEL_Buckley2026_b_o = "gs_buckley2026_b_o"; + constexpr const char *LABEL_Buckley2026_psi_c = "gs_buckley2026_psi_c"; + + // Global data through which the Buckley et al. (2026) model asks the photosynthesis and plant hydraulics plug-ins for its inputs + constexpr const char *GLOBAL_Buckley2026_enabled = "buckley2026_stomatal_model_enabled"; + //! Boundary-layer conductance that, in series with the total stomatal conductance, carries the same flux as the stomata of each leaf face in series with the boundary layer of that face /** * The stomatal conductance is split between the front and back faces in the proportions \f$\zeta\f$ and \f$1-\zeta\f$, and each face exchanges gas only through its own boundary layer. The surface @@ -144,6 +153,8 @@ StomatalConductanceModel::StomatalConductanceModel(helios::Context *m_context) { beta_default = 1.f; // unitless + soil_water_potential_default = -0.1; // MPa + model = "BMF"; // default model - Buckley, Mott, Farquhar } @@ -414,6 +425,109 @@ void StomatalConductanceModel::setModelCoefficients(const BBcoefficients &coeffs model = "BB"; } +void StomatalConductanceModel::selectBuckley2026Model() { + model = "Buckley2026"; + // The request is never withdrawn automatically: another StomatalConductanceModel instance may still use the Buckley et al. (2026) model, and an unneeded request only costs extra + // outputs from the other plug-ins, whereas a missing one would leave this model without its inputs. + context->setGlobalData(GLOBAL_Buckley2026_enabled, 1); +} + +void StomatalConductanceModel::setModelCoefficients(const Buckley2026coefficients &coeffs) { + Buckley2026coeffs = coeffs; + Buckley2026model_coefficients.clear(); + selectBuckley2026Model(); +} + +void StomatalConductanceModel::setModelCoefficients(const Buckley2026coefficients &coeffs, const std::vector &UUIDs) { + for (uint UUID: UUIDs) { + Buckley2026model_coefficients[UUID] = coeffs; + } + selectBuckley2026Model(); +} + +void StomatalConductanceModel::setModelCoefficients(const std::string &material_label, const Buckley2026coefficients &coeffs) { + context->setMaterialData(material_label, LABEL_Buckley2026_s, coeffs.s); + context->setMaterialData(material_label, LABEL_Buckley2026_K_plant_at_25C, coeffs.K_plant_at_25C); + context->setMaterialData(material_label, LABEL_Buckley2026_b_o, coeffs.b_o); + context->setMaterialData(material_label, LABEL_Buckley2026_psi_c, coeffs.psi_c); + + // Clear cache for this material + uint matID = context->getMaterialIDFromLabel(material_label); + material_coefficient_cache_Buckley2026.erase(matID); + + selectBuckley2026Model(); +} + +void StomatalConductanceModel::setBuckley2026CoefficientsFromLibrary(const std::string &species) { + setModelCoefficients(getBuckley2026CoefficientsFromLibrary(species)); +} + +void StomatalConductanceModel::setBuckley2026CoefficientsFromLibrary(const std::string &species, const std::vector &UUIDs) { + setModelCoefficients(getBuckley2026CoefficientsFromLibrary(species), UUIDs); +} + +void StomatalConductanceModel::setBuckley2026CoefficientsFromLibrary(const std::string &species, const std::string &material_label) { + setModelCoefficients(material_label, getBuckley2026CoefficientsFromLibrary(species)); +} + +Buckley2026coefficients StomatalConductanceModel::getBuckley2026CoefficientsFromLibrary(const std::string &species) { + // Fitted to LI-600 porometer diurnals (gsw corrected per Rizzo & Bailey) with the measured predawn or first stem water potential of each day as psi_soil, and psi_c from the + // pressure-volume turgor loss point (Buckley et al. 2026, Section 2.2) where available (fitted for eastern redbud), with A, J, and Rd from the PhotosynthesisModel Farquhar library entry named + // in each comment. Use each entry with that Farquhar entry: s and b_o depend on the photosynthesis parameters. + Buckley2026coefficients coeffs; + const std::string &s = species; + if (s == "Almond" || s == "almond") { + // fitted with Farquhar library entry Almond; psi_c turgor loss point, literature; R2 of date x hour means 0.75 (leave-one-hour-out 0.56) + coeffs.s = 0.01595f; + coeffs.K_plant_at_25C = 0.03266f; + coeffs.b_o = 1.416f; + coeffs.psi_c = -2.329f; + } else if (s == "Elderberry" || s == "elderberry" || s == "blue_elderberry") { + // fitted with Farquhar library entry Elderberry; psi_c turgor loss point; R2 of date x hour means 0.83 (leave-one-hour-out 0.77) + coeffs.s = 0.0271f; + coeffs.K_plant_at_25C = 0.007491f; + coeffs.b_o = 1.27f; + coeffs.psi_c = -2.752f; + } else if (s == "Cowpea" || s == "cowpea") { + // fitted with Farquhar library entry Cowpea; psi_c turgor loss point; R2 of date x hour means 0.76 (leave-one-hour-out 0.65) + coeffs.s = 0.2626f; + coeffs.K_plant_at_25C = 0.05714f; + coeffs.b_o = 1.074f; + coeffs.psi_c = -0.846f; + } else if (s == "EasternRedbud" || s == "easternredbud") { + // fitted with Farquhar library entry EasternRedbudSunlit and EasternRedbudShaded; psi_c fitted; R2 of date x hour means 0.63 (leave-one-hour-out 0.57) + coeffs.s = 0.004195f; + coeffs.K_plant_at_25C = 0.001744f; + coeffs.b_o = 1.774f; + coeffs.psi_c = -3.667f; + } else if (s == "Grape" || s == "grape") { + // fitted with Farquhar library entry Grape; psi_c turgor loss point, literature; R2 of date x hour means 0.55 (leave-one-hour-out 0.16) + coeffs.s = 0.03341f; + coeffs.K_plant_at_25C = 0.007471f; + coeffs.b_o = 1.268f; + coeffs.psi_c = -1.718f; + } else if (s == "Olive" || s == "olive") { + // fitted with Farquhar library entry Olive; psi_c turgor loss point; R2 of date x hour means 0.88 (leave-one-hour-out 0.73) + coeffs.s = 0.02156f; + coeffs.K_plant_at_25C = 0.004562f; + coeffs.b_o = 5.033f; + coeffs.psi_c = -2.043f; + } else if (s == "Pistachio" || s == "pistachio" || s == "PistachioFemale" || s == "pistachiofemale" || s == "pistachio_female" || s == "Pistachio_Female" || s == "Pistachio_female") { + // fitted with Farquhar library entry PistachioFemale (or Pistachio); psi_c turgor loss point; R2 of date x hour means 0.76 (leave-one-hour-out 0.60) + coeffs.s = 0.08766f; + coeffs.K_plant_at_25C = 0.007411f; + coeffs.b_o = 1.11f; + coeffs.psi_c = -4.103f; + } else { + helios_runtime_error("ERROR (StomatalConductanceModel::getBuckley2026CoefficientsFromLibrary): Unknown species '" + s + + "'. Valid species are: Almond, Elderberry (or blue_elderberry), Cowpea, EasternRedbud, Grape, Olive, Pistachio (or PistachioFemale). Names are also accepted in all lower case."); + } + if (message_flag) { + std::cout << "Returning Buckley et al. (2026) stomatal model coefficients for " << s << std::endl; + } + return coeffs; +} + void StomatalConductanceModel::validateDynamicTimeConstants(float tau_open, float tau_close) { // The time constants appear in the denominator of the forward Euler update gs = gs_old + (gs_ss-gs_old)*dt/tau. A zero value gives a non-finite conductance, and a negative value inverts the // direction of the relaxation so that stomata diverge away from the steady-state value. @@ -616,6 +730,41 @@ BBcoefficients StomatalConductanceModel::getCoefficientsForPrimitive_BB(uint UUI return BBcoeffs; } +Buckley2026coefficients StomatalConductanceModel::getCoefficientsForPrimitive_Buckley2026(uint UUID) const { + uint materialID = context->getPrimitiveMaterialID(UUID); + + if (material_coefficient_cache_Buckley2026.find(materialID) != material_coefficient_cache_Buckley2026.end()) { + return material_coefficient_cache_Buckley2026.at(materialID); + } + + Buckley2026coefficients coeffs; + bool found = true; + + try { + const Material &mat = context->getMaterial(materialID); + if (mat.doesMaterialDataExist(LABEL_Buckley2026_s) && mat.doesMaterialDataExist(LABEL_Buckley2026_K_plant_at_25C) && mat.doesMaterialDataExist(LABEL_Buckley2026_b_o) && mat.doesMaterialDataExist(LABEL_Buckley2026_psi_c)) { + + mat.getMaterialData(LABEL_Buckley2026_s, coeffs.s); + mat.getMaterialData(LABEL_Buckley2026_K_plant_at_25C, coeffs.K_plant_at_25C); + mat.getMaterialData(LABEL_Buckley2026_b_o, coeffs.b_o); + mat.getMaterialData(LABEL_Buckley2026_psi_c, coeffs.psi_c); + + material_coefficient_cache_Buckley2026[materialID] = coeffs; + return coeffs; + } else { + found = false; + } + } catch (const std::exception &) { + found = false; + } + + if (!found && Buckley2026model_coefficients.find(UUID) != Buckley2026model_coefficients.end()) { + return Buckley2026model_coefficients.at(UUID); + } + + return Buckley2026coeffs; +} + void StomatalConductanceModel::run() { run(context->getAllUUIDs()); } @@ -636,6 +785,13 @@ void StomatalConductanceModel::run(const std::vector &UUIDs, float dt) { size_t assumed_default_An = 0; size_t assumed_default_Gamma = 0; + size_t assumed_default_psi_soil = 0; + + // Request the Buckley et al. (2026) inputs from the photosynthesis and plant hydraulics plug-ins (also done when the model is selected; repeated in case the global data was cleared) + if (model == "Buckley2026") { + context->setGlobalData(GLOBAL_Buckley2026_enabled, 1); + } + bool warn_dt_too_large = false; bool warn_tau_unspecified = false; bool warn_old_gs_unspecified = false; @@ -819,6 +975,38 @@ void StomatalConductanceModel::run(const std::vector &UUIDs, float dt) { } } + // Buckley et al. (2026) model inputs: net photosynthesis, electron transport rate, and day respiration from the photosynthesis model, and the soil water potential and whole-plant + // hydraulic conductance K_plant from the plant hydraulics model (the coefficient K_plant_at_25C, scaled to leaf temperature, is used when K_plant is not available) + float J = 0.f; + float Rd = 0.f; + float psi_soil = soil_water_potential_default; + float K_plant = -1.f; // < 0: not provided as primitive data + if (model == "Buckley2026") { + // A, J, and Rd have no meaningful defaults (b = J/4 - Rd - A + b_o drives the model), so they are required, as "electron_transport_ratio" is for SIF in the radiation model + for (const char *label: {"net_photosynthesis", "electron_transport_rate", "dark_respiration"}) { + if (!context->doesPrimitiveDataExist(UUID, label) || context->getPrimitiveDataType(label) != HELIOS_TYPE_FLOAT) { + helios_runtime_error("ERROR (StomatalConductanceModel::run): The Buckley et al. (2026) model requires primitive data \"" + std::string(label) + "\", which does not exist for primitive " + std::to_string(UUID) + + ". Run the PhotosynthesisModel with the Farquhar (C3) model before the stomatal conductance model: once the Buckley et al. (2026) model is selected it writes " + "\"net_photosynthesis\", \"electron_transport_rate\", and \"dark_respiration\" automatically. The C4 and empirical photosynthesis models do not provide them."); + } + } + context->getPrimitiveData(UUID, "net_photosynthesis", An); + context->getPrimitiveData(UUID, "electron_transport_rate", J); + context->getPrimitiveData(UUID, "dark_respiration", Rd); + if (context->doesPrimitiveDataExist(UUID, "soil_water_potential") && context->getPrimitiveDataType("soil_water_potential") == HELIOS_TYPE_FLOAT) { + context->getPrimitiveData(UUID, "soil_water_potential", psi_soil); + } else { + assumed_default_psi_soil++; + } + if (context->doesPrimitiveDataExist(UUID, "plant_hydraulic_conductance") && context->getPrimitiveDataType("plant_hydraulic_conductance") == HELIOS_TYPE_FLOAT) { + context->getPrimitiveData(UUID, "plant_hydraulic_conductance", K_plant); + if (!std::isfinite(K_plant) || K_plant < 0.f) { + helios_runtime_error("ERROR (StomatalConductanceModel::run): The whole-plant hydraulic conductance for primitive " + std::to_string(UUID) + " is " + std::to_string(K_plant) + + " mol/m^2/s/MPa, but it must be a finite non-negative value. Check the primitive data \"plant_hydraulic_conductance\", or remove it to use the model coefficient K_plant_at_25C."); + } + } + } + // Solve for the steady-state conductance. The leaf-surface balance and the conductance model are evaluated for given effective boundary-layer conductances to water vapor and CO2 (see // effectiveBoundaryLayerConductance()). Returns false, without evaluating gs, if the resulting leaf-surface CO2 concentration is not physically meaningful for the model: Cs <= 0, which would flip @@ -910,6 +1098,35 @@ void StomatalConductanceModel::run(const std::vector &UUIDs, float dt) { context->setPrimitiveData(UUID, "chi_BB", coeffs.chi); } + } else if (model == "Buckley2026") { + + // model coefficients + Buckley2026coefficients coeffs = getCoefficientsForPrimitive_Buckley2026(UUID); + + // whole-plant hydraulic conductance at the leaf temperature (mol/m^2/s/MPa) + const float K_plant_leaf = (K_plant >= 0.f) ? K_plant : coeffs.K_plant_at_25C * waterFluidityRelativeTo25C(TL); + + // b = J/4 - Rd - A + b_o, where J/4 - Rd is the CO2-saturated electron-transport-limited photosynthetic rate (umol/m^2/s) + const float b = max(0.f, 0.25f * J - Rd - An + coeffs.b_o); + + std::vector variables{b, es, ea, gbw_effective, press, beta, psi_soil, K_plant_leaf}; + + float esurf = solveSurfaceVaporPressure(evaluate_Buckley2026model, variables, &coeffs, ea, es, warnings); + float Ds = max(0.f, (es - esurf) / press); // mol/mol + + gs_result = evaluate_Buckley2026conductance(Ds, b, psi_soil, K_plant_leaf, beta, coeffs); + + if (write_outputs && std::find(output_prim_data.begin(), output_prim_data.end(), "model_parameters") != output_prim_data.end()) { + context->setPrimitiveData(UUID, "s_Buckley2026", coeffs.s); + context->setPrimitiveData(UUID, "K_plant_Buckley2026", K_plant_leaf); + context->setPrimitiveData(UUID, "b_o_Buckley2026", coeffs.b_o); + context->setPrimitiveData(UUID, "psi_c_Buckley2026", coeffs.psi_c); + } + + if (write_outputs && std::find(output_prim_data.begin(), output_prim_data.end(), "vapor_pressure_deficit") != output_prim_data.end()) { + context->setPrimitiveData(UUID, "vapor_pressure_deficit", Ds * 1000.f); // mmol/mol + } + } else { // model coefficients @@ -1089,6 +1306,11 @@ void StomatalConductanceModel::run(const std::vector &UUIDs, float dt) { std::cout << "WARNING (StomatalConductanceModel::run): The Ball-Berry-Leuning stomatal conductance model requires net photosynthesis, but primitive data net_photosynthesis could not be found for " << assumed_default_An << " primitives. Did you forget to run the photosynthesis model?" << std::endl; } + if (model == "Buckley2026" && assumed_default_psi_soil > 0) { + std::cout << "WARNING (StomatalConductanceModel::run): The Buckley et al. (2026) stomatal conductance model requires the soil water potential, but primitive data soil_water_potential could not be found for " << assumed_default_psi_soil + << " primitives, and the default value of " << soil_water_potential_default + << " MPa was used. Set it directly, or run the PlantHydraulicsModel after setting the soil water potential of each plant with setSoilWaterPotentialOfPlant()." << std::endl; + } if (model == "BBL" && assumed_default_Gamma > 0) { std::cout << "WARNING (StomatalConductanceModel::run): The Ball-Berry-Leuning stomatal conductance model requires the CO2 compensation point Gamma, but primitive data Gamma_CO2 could not be found for " << assumed_default_Gamma << " primitives. Did you forget to set optional output primitive data Gamma_CO2 in the photosynthesis model?" << std::endl; @@ -1203,6 +1425,47 @@ float StomatalConductanceModel::evaluate_BMFmodel(float esurf, std::vector 0.f) || !(dpsi > 0.f) || !(K_plant > 0.f)) { + return 0.f; // no guard-cell advantage, soil at or below the closure potential, or no water supply + } + const float gs = sb * dpsi / (1.f + sb * Ds / K_plant); + + // The soil moisture factor scales the result in the same way as for the other models. + return max(0.f, beta * gs); +} + +float StomatalConductanceModel::evaluate_Buckley2026model(float esurf, std::vector &variables, const void *parameters) { + + // We want to find the vapor pressure at the surface, esurf, that balances the equation gs(esurf)*(es-esurf) = gbw*(esurf-ea). This function returns the residual of this equation. + + const auto *coeffs = reinterpret_cast(parameters); + + float b = variables[0]; + float es = variables[1]; + float ea = variables[2]; + float gbw = variables[3]; + float press = variables[4]; + float beta = variables[5]; + float psi_soil = variables[6]; + float K_plant = variables[7]; + + float Ds = max(0.f, (es - esurf) / press); // mol/mol + float gs = evaluate_Buckley2026conductance(Ds, b, psi_soil, K_plant, beta, *coeffs); + + return gs * (es - esurf) - gbw * (esurf - ea); +} + float StomatalConductanceModel::evaluate_BBconductance(float Ds, float Psix, float i, float beta, const BBcoefficients &coeffs) { // Fraction of the epidermal water potential transmitted to the guard cells, and the xylem-to-epidermis hydraulic resistance (MPa m^2 s / mol). These are treated as fixed in this simplified form of @@ -1296,6 +1559,10 @@ void StomatalConductanceModel::printDefaultValueReport(const std::vector & size_t assumed_default_CO2 = 0; size_t assumed_default_gbw = 0; size_t assumed_default_beta = 0; + size_t assumed_default_J = 0; + size_t assumed_default_Rd = 0; + size_t assumed_default_psi_soil = 0; + size_t assumed_coefficient_K_plant = 0; size_t Nprimitives = UUIDs.size(); @@ -1358,6 +1625,22 @@ void StomatalConductanceModel::printDefaultValueReport(const std::vector & if (!context->doesPrimitiveDataExist(UUID, "xylem_water_potential") || context->getPrimitiveDataType("xylem_water_potential") != HELIOS_TYPE_FLOAT) { assumed_default_Psix++; } + } else if (model == "Buckley2026") { + if (!context->doesPrimitiveDataExist(UUID, "net_photosynthesis") || context->getPrimitiveDataType("net_photosynthesis") != HELIOS_TYPE_FLOAT) { + assumed_default_An++; + } + if (!context->doesPrimitiveDataExist(UUID, "electron_transport_rate") || context->getPrimitiveDataType("electron_transport_rate") != HELIOS_TYPE_FLOAT) { + assumed_default_J++; + } + if (!context->doesPrimitiveDataExist(UUID, "dark_respiration") || context->getPrimitiveDataType("dark_respiration") != HELIOS_TYPE_FLOAT) { + assumed_default_Rd++; + } + if (!context->doesPrimitiveDataExist(UUID, "soil_water_potential") || context->getPrimitiveDataType("soil_water_potential") != HELIOS_TYPE_FLOAT) { + assumed_default_psi_soil++; + } + if (!context->doesPrimitiveDataExist(UUID, "plant_hydraulic_conductance") || context->getPrimitiveDataType("plant_hydraulic_conductance") != HELIOS_TYPE_FLOAT) { + assumed_coefficient_K_plant++; + } } } @@ -1423,6 +1706,15 @@ void StomatalConductanceModel::printDefaultValueReport(const std::vector & << std::endl; } + if (model == "Buckley2026") { + std::cout << "net photosynthesis: " << assumed_default_An << " of " << Nprimitives << " used default value of " << An_default << " because net_photosynthesis primitive data did not exist" << std::endl; + std::cout << "electron transport rate: " << assumed_default_J << " of " << Nprimitives << " used default value of 0 because electron_transport_rate primitive data did not exist" << std::endl; + std::cout << "day respiration: " << assumed_default_Rd << " of " << Nprimitives << " used default value of 0 because dark_respiration primitive data did not exist" << std::endl; + std::cout << "soil water potential: " << assumed_default_psi_soil << " of " << Nprimitives << " used default value of " << soil_water_potential_default << " because soil_water_potential primitive data did not exist" << std::endl; + std::cout << "whole-plant hydraulic conductance: " << assumed_coefficient_K_plant << " of " << Nprimitives + << " used the model coefficient K_plant_at_25C (scaled to leaf temperature) because plant_hydraulic_conductance primitive data did not exist" << std::endl; + } + std::cout << "------------------------------------------------------" << std::endl; } diff --git a/plugins/stomatalconductance/tests/selfTest.cpp b/plugins/stomatalconductance/tests/selfTest.cpp index c9101cc50..43fde685b 100644 --- a/plugins/stomatalconductance/tests/selfTest.cpp +++ b/plugins/stomatalconductance/tests/selfTest.cpp @@ -1499,3 +1499,240 @@ DOCTEST_TEST_CASE("StomatalConductanceModel - Intermediate stomatal sidedness do // An assimilation rate that cannot be supplied even with both faces' full boundary-layer conductance is still rejected. DOCTEST_CHECK_THROWS(gs = runModel(60.f)); } + +// --------------------------------------------------------------------------------------------------------------------- +// Buckley et al. (2026) model +// --------------------------------------------------------------------------------------------------------------------- + +namespace { + // A leaf at air temperature with a very large boundary-layer conductance, so that the leaf-surface vapor pressure is that of the air and Ds = (es - ea)/P + uint addBuckleyTestLeaf(Context &context, float TL, float rh, float A, float J, float Rd, float psi_soil) { + uint UUID = context.addPatch(); + context.setPrimitiveData(UUID, "temperature", TL); + context.setPrimitiveData(UUID, "air_temperature", TL); + context.setPrimitiveData(UUID, "air_humidity", rh); + context.setPrimitiveData(UUID, "air_pressure", 101000.f); + context.setPrimitiveData(UUID, "boundarylayer_conductance", 1.0e5f); + context.setPrimitiveData(UUID, "net_photosynthesis", A); + context.setPrimitiveData(UUID, "electron_transport_rate", J); + context.setPrimitiveData(UUID, "dark_respiration", Rd); + context.setPrimitiveData(UUID, "soil_water_potential", psi_soil); + return UUID; + } + + float buckleyTestVPD(float TL, float rh) { // mol/mol, with the plug-in's saturation vapor pressure + const float es = 611.f * std::exp(17.502f * (TL - 273.f) / ((TL - 273.f) + 240.97f)); + return (es - rh * es) / 101000.f; + } + + float buckleyClosedForm(const Buckley2026coefficients &c, float A, float J, float Rd, float psi_soil, float K_plant, float Ds) { + const float sb = c.s * std::max(0.f, 0.25f * J - Rd - A + c.b_o); + return std::max(0.f, sb * (psi_soil - c.psi_c) / (1.f + sb * Ds / K_plant)); + } +} // namespace + +DOCTEST_TEST_CASE("StomatalConductanceModel - Buckley2026 reproduces the coupled steady state of the fitting code") { + // Reference values: steady state of Buckley et al. (2026) Eqn 3a coupled to the Helios Farquhar model (PhotosynthesisModel library entry "PistachioFemale"), solved independently + // in Python with the code used to fit the species library. Given that solution's A, J, and Rd, the stomatal model must return its conductance. + Context context; + StomatalConductanceModel gsm(&context); + gsm.disableMessages(); + gsm.setBuckley2026CoefficientsFromLibrary("Pistachio"); + + uint leaf1 = addBuckleyTestLeaf(context, 303.15f, 0.4f, 28.16217f, 183.7890f, 2.04238f, -0.5f); + uint leaf2 = addBuckleyTestLeaf(context, 298.15f, 0.6f, 9.08960f, 57.7394f, 1.5f, -1.5f); + gsm.run(); + + float gs1, gs2; + context.getPrimitiveData(leaf1, "moisture_conductance", gs1); + context.getPrimitiveData(leaf2, "moisture_conductance", gs2); + DOCTEST_CHECK(gs1 == doctest::Approx(0.959348f).epsilon(1e-3)); + DOCTEST_CHECK(gs2 == doctest::Approx(0.649252f).epsilon(1e-3)); +} + +DOCTEST_TEST_CASE("StomatalConductanceModel - Buckley2026 matches its closed form and uses K_plant from primitive data") { + Context context; + StomatalConductanceModel gsm(&context); + gsm.disableMessages(); + + Buckley2026coefficients coeffs; + coeffs.s = 0.05f; + coeffs.K_plant_at_25C = 0.02f; + coeffs.b_o = 1.0f; + coeffs.psi_c = -2.0f; + gsm.setModelCoefficients(coeffs); + + const float TL = 298.15f, rh = 0.5f, A = 15.f, J = 120.f, Rd = 1.2f, psi_soil = -0.4f; + uint leaf_coefficient = addBuckleyTestLeaf(context, TL, rh, A, J, Rd, psi_soil); + uint leaf_hydraulics = addBuckleyTestLeaf(context, TL, rh, A, J, Rd, psi_soil); + const float K_plant_hydraulics = 0.005f; // mol/m^2/s/MPa, as written by the plant hydraulics plug-in + context.setPrimitiveData(leaf_hydraulics, "plant_hydraulic_conductance", K_plant_hydraulics); + gsm.run(); + + const float Ds = buckleyTestVPD(TL, rh); + const float fluidity_25C = 4.891e10f * std::exp(-4209.f / TL - 0.04527f * TL + 3.376e-5f * TL * TL); + float gs_coefficient, gs_hydraulics; + context.getPrimitiveData(leaf_coefficient, "moisture_conductance", gs_coefficient); + context.getPrimitiveData(leaf_hydraulics, "moisture_conductance", gs_hydraulics); + DOCTEST_CHECK(gs_coefficient == doctest::Approx(buckleyClosedForm(coeffs, A, J, Rd, psi_soil, coeffs.K_plant_at_25C * fluidity_25C, Ds)).epsilon(1e-3)); + DOCTEST_CHECK(gs_hydraulics == doctest::Approx(buckleyClosedForm(coeffs, A, J, Rd, psi_soil, K_plant_hydraulics, Ds)).epsilon(1e-3)); + DOCTEST_CHECK(gs_hydraulics < gs_coefficient); // the smaller K_plant limits the conductance +} + +DOCTEST_TEST_CASE("StomatalConductanceModel - Buckley2026 responds to soil water potential, VPD, and photosynthesis") { + Context context; + StomatalConductanceModel gsm(&context); + gsm.disableMessages(); + gsm.setBuckley2026CoefficientsFromLibrary("Almond"); + const Buckley2026coefficients coeffs = gsm.getBuckley2026CoefficientsFromLibrary("Almond"); + + uint wet = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + uint dry = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -1.2f); + uint closed = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, coeffs.psi_c - 0.1f); // soil below the closure potential + uint humid = addBuckleyTestLeaf(context, 298.15f, 0.9f, 10.f, 100.f, 1.f, -0.2f); + uint high_A = addBuckleyTestLeaf(context, 298.15f, 0.5f, 20.f, 100.f, 1.f, -0.2f); // same capacity J/4 - Rd, smaller excess b + gsm.run(); + + auto gs = [&](uint UUID) { + float value; + context.getPrimitiveData(UUID, "moisture_conductance", value); + return value; + }; + DOCTEST_CHECK(gs(wet) > 0.f); + DOCTEST_CHECK(gs(dry) < gs(wet)); + DOCTEST_CHECK(gs(closed) == 0.f); + DOCTEST_CHECK(gs(humid) > gs(wet)); + DOCTEST_CHECK(gs(high_A) < gs(wet)); +} + +DOCTEST_TEST_CASE("StomatalConductanceModel - Buckley2026 requests its inputs through global data") { + Context context; + StomatalConductanceModel gsm(&context); + gsm.disableMessages(); + uint leaf = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + + auto flag = [&]() { + int value = -1; + if (context.doesGlobalDataExist("buckley2026_stomatal_model_enabled")) { + context.getGlobalData("buckley2026_stomatal_model_enabled", value); + } + return value; + }; + DOCTEST_CHECK(flag() == -1); // not requested by the default model + + gsm.setBuckley2026CoefficientsFromLibrary("Olive"); + DOCTEST_CHECK(flag() == 1); + + // A second instance using another model does not withdraw the request + StomatalConductanceModel gsm_other(&context); + gsm_other.disableMessages(); + gsm_other.setModelCoefficients(BMFcoefficients()); + gsm_other.run(std::vector{leaf}); + DOCTEST_CHECK(flag() == 1); + + // The request is not withdrawn automatically, even when the instance that made it selects another model, since another instance may still use the Buckley et al. (2026) model + gsm.setModelCoefficients(BMFcoefficients()); + gsm.run(std::vector{leaf}); + DOCTEST_CHECK(flag() == 1); + + // UUID- and material-based coefficients also make the request, and run() restores it if the global data was cleared + for (int setter = 0; setter < 3; setter++) { + context.setGlobalData("buckley2026_stomatal_model_enabled", 0); + StomatalConductanceModel gsm_new(&context); + gsm_new.disableMessages(); + if (setter == 0) { + gsm_new.setModelCoefficients(Buckley2026coefficients(), std::vector{leaf}); + } else if (setter == 1) { + context.addMaterial("buckley_leaf_material"); + gsm_new.setBuckley2026CoefficientsFromLibrary("Grape", "buckley_leaf_material"); + } else { + gsm_new.setBuckley2026CoefficientsFromLibrary("Grape"); + context.setGlobalData("buckley2026_stomatal_model_enabled", 0); + gsm_new.run(std::vector{leaf}); + } + DOCTEST_CHECK(flag() == 1); + } +} + +DOCTEST_TEST_CASE("StomatalConductanceModel - Buckley2026 library, UUID, and material coefficients") { + Context context; + StomatalConductanceModel gsm(&context); + gsm.disableMessages(); + + for (const std::string species: {"Almond", "Elderberry", "blue_elderberry", "Cowpea", "EasternRedbud", "Grape", "Olive", "Pistachio", "PistachioFemale", "pistachio"}) { + Buckley2026coefficients coeffs; + DOCTEST_CHECK_NOTHROW(coeffs = gsm.getBuckley2026CoefficientsFromLibrary(species)); + DOCTEST_CHECK(coeffs.s > 0.f); + DOCTEST_CHECK(coeffs.K_plant_at_25C > 0.f); + DOCTEST_CHECK(coeffs.b_o >= 0.f); + DOCTEST_CHECK(coeffs.psi_c < 0.f); + } + DOCTEST_CHECK_THROWS_AS(gsm.getBuckley2026CoefficientsFromLibrary("UnknownSpecies"), std::runtime_error); + DOCTEST_CHECK_THROWS_AS(gsm.setBuckley2026CoefficientsFromLibrary("UnknownSpecies"), std::runtime_error); + + // The default coefficients are the pistachio entry + const Buckley2026coefficients pistachio = gsm.getBuckley2026CoefficientsFromLibrary("Pistachio"); + const Buckley2026coefficients defaults; + DOCTEST_CHECK(defaults.s == pistachio.s); + DOCTEST_CHECK(defaults.K_plant_at_25C == pistachio.K_plant_at_25C); + DOCTEST_CHECK(defaults.b_o == pistachio.b_o); + DOCTEST_CHECK(defaults.psi_c == pistachio.psi_c); + + // Coefficients by UUID and by material, reported through the optional model_parameters output + uint leaf_global = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + uint leaf_uuid = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + uint leaf_material = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + context.addMaterial("olive_leaves"); + context.assignMaterialToPrimitive(leaf_material, "olive_leaves"); + + gsm.optionalOutputPrimitiveData("model_parameters"); + gsm.optionalOutputPrimitiveData("vapor_pressure_deficit"); + gsm.setBuckley2026CoefficientsFromLibrary("Almond"); + gsm.setBuckley2026CoefficientsFromLibrary("Cowpea", std::vector{leaf_uuid}); + gsm.setBuckley2026CoefficientsFromLibrary("Olive", "olive_leaves"); + gsm.run(); + + auto s_of = [&](uint UUID) { + float value; + context.getPrimitiveData(UUID, "s_Buckley2026", value); + return value; + }; + DOCTEST_CHECK(s_of(leaf_global) == gsm.getBuckley2026CoefficientsFromLibrary("Almond").s); + DOCTEST_CHECK(s_of(leaf_uuid) == gsm.getBuckley2026CoefficientsFromLibrary("Cowpea").s); + DOCTEST_CHECK(s_of(leaf_material) == gsm.getBuckley2026CoefficientsFromLibrary("Olive").s); + for (const char *label: {"K_plant_Buckley2026", "b_o_Buckley2026", "psi_c_Buckley2026", "vapor_pressure_deficit"}) { + DOCTEST_CHECK(context.doesPrimitiveDataExist(leaf_global, label)); + } + float Ds_out; + context.getPrimitiveData(leaf_global, "vapor_pressure_deficit", Ds_out); + DOCTEST_CHECK(Ds_out == doctest::Approx(buckleyTestVPD(298.15f, 0.5f) * 1000.f).epsilon(1e-3)); // mmol/mol +} + +DOCTEST_TEST_CASE("StomatalConductanceModel - Buckley2026 without photosynthesis outputs, and invalid K_plant") { + Context context; + StomatalConductanceModel gsm(&context); + gsm.disableMessages(); + gsm.setBuckley2026CoefficientsFromLibrary("Pistachio"); + + // A, J, and Rd are required (no meaningful defaults): missing any of them is an error, e.g. photosynthesis not run, or a C4 or empirical photosynthesis model + uint bare = context.addPatch(); + DOCTEST_CHECK_THROWS_AS(gsm.run(std::vector{bare}), std::runtime_error); + for (const char *missing: {"net_photosynthesis", "electron_transport_rate", "dark_respiration"}) { + uint leaf_missing = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + context.clearPrimitiveData(leaf_missing, missing); + DOCTEST_CHECK_THROWS_AS(gsm.run(std::vector{leaf_missing}), std::runtime_error); + } + DOCTEST_CHECK_NOTHROW(gsm.printDefaultValueReport(std::vector{bare})); + + // A missing soil water potential falls back to the default (with a warning), and K_plant to the coefficient + uint no_soil = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + context.clearPrimitiveData(no_soil, "soil_water_potential"); + DOCTEST_CHECK_NOTHROW(gsm.run(std::vector{no_soil})); + float gs; + context.getPrimitiveData(no_soil, "moisture_conductance", gs); + DOCTEST_CHECK(gs > 0.f); + + uint leaf = addBuckleyTestLeaf(context, 298.15f, 0.5f, 10.f, 100.f, 1.f, -0.2f); + context.setPrimitiveData(leaf, "plant_hydraulic_conductance", -0.01f); + DOCTEST_CHECK_THROWS_AS(gsm.run(std::vector{leaf}), std::runtime_error); +}