Non-steady-state stomatal conductance modeling and its implications: from leaf to ecosystem

Accurate and efficient modeling of stomatal conductance (g_s) has been a key challenge in vegetation models across scales. Current practice of most land surface models (LSMs) assumes steady-state g_s and predicts stomatal responses to environmental cues as immediate jumps between stationary regimes. However, the response of stomata can be orders of magnitude slower than that of photosynthesis and often cannot reach a steady state before the next model time step, even on half-hourly timescales. Here, we implemented a simple dynamic g_s model in the vegetation module of an LSM developed within the Climate Modeling Alliance and investigated the potential biases caused by the steady-state assumption from leaf to canopy scales. In comparison with steady-state models, the dynamic model better predicted the coupled temporal response of photosynthesis and stomatal conductance to changes in light intensity using leaf measurements. In ecosystem flux simulations, while the impact of g_s hysteresis response may not be substantial in terms of monthly integrated fluxes, our results highlight the importance of considering this effect when quantifying fluxes in the mornings and evenings, as well as interpreting diurnal hysteresis patterns observed in ecosystem fluxes. Simulations also indicate that the biases in the integrated fluxes are more significant when stomata exhibit different speeds for opening and closure. Furthermore, prognostic modeling can bypass the A-C_i iterations required for steady-state simulations and can be robustly run with comparable computational costs. Overall, our study demonstrates the implications of dynamic g_s modeling for improving the accuracy and efficiency of LSMs and for advancing our understanding of plant–environment interactions.

1 Introduction

Modeling stomatal conductance (g_s), the opening and closure of tiny pores on leaves, is one of the key elements and challenges in land surface models (LSMs). Stomata regulate the gas exchange rates of plants, allowing the uptake of CO₂ for photosynthetic assimilation while constraining water loss through transpiration (Berry et al., 2010; Damour et al., 2010). The behavior of stomata, especially their responses to environmental variations, plays a significant role in determining the fluxes of carbon, water, and energy between vegetated surfaces and the atmosphere (Berry et al., 2010; Buckley, 2017). Therefore, accurate and efficient modeling of g_s is important for understanding the current Earth system and projecting future changes.

Stomatal conductance has been traditionally predicted with empirical models, relating g_s to photosynthesis rate and environmental cues with estimated parameters from statistical regressions (Ball et al., 1987; Leuning, 1990, 1995; Medlyn et al., 2011; Damour et al., 2010). Efforts have also been made to constrain stomatal behavior from the principle of optimizing the trade-offs between carbon gain with the related penalty of stomatal opening (Wolf et al., 2016; Venturas et al., 2018; Wang et al., 2020). Additional understanding of stomatal response includes plant hydraulic models that consider the transport of water from soil through plants into the atmosphere (soil–plant–atmosphere continuum, SPAC) (Sperry et al., 1998, 2002; Bonan et al., 2014). However, most existing stomatal models, especially those currently used to scale from leaf to canopy level and implemented in LSMs, assume steady states (Vialet-Chabrand et al., 2017). They predict the opening and closure of stomata in stationary regimes, modeling stomatal response to environmental variations as immediate jumps between states (Vialet-Chabrand et al., 2013).

While steady-state models assume that stomatal conductance changes instantaneously with the environment, the temporal response of stomata in reality can often be an order of magnitude slower than the biochemical response of photosynthesis (Pearcy and Seemann, 1990; Vialet-Chabrand et al., 2013). Plants can experience frequent environmental changes on a timescale of seconds, such as light fluctuations due to cloud cover and canopy shading. Meanwhile, stomatal response times vary from minutes to more than an hour. Thus, a steady state is often not reached when environmental conditions change faster than stomata can respond to (Lawson and Blatt, 2014; Vialet-Chabrand et al., 2017). This slower response of g_s could further impose regulations on assimilation rate via its effects on intercellular CO₂ concentration (C_i), notably under rapidly changing incident radiation in natural environments (Kaiser and Kappen, 2000; Vialet-Chabrand et al., 2017). The mismatch and interaction of photosynthesis and stomatal response could lead to temporal variations in water use efficiency (WUE) as well (Lawson et al., 2011; Venturas et al., 2018). These can all lead to biases, and it is important to consider non-steady-state temporal responses of g_s for more accurate predictions of ecosystem fluxes. Additionally, the inclusion of this factor may also contribute to the hysteresis of plant responses and ecosystem fluxes observed in natural diurnal cycles (Vialet-Chabrand et al., 2013); e.g., evapotranspiration (ET) rates tend to be higher in the afternoon under the same incoming radiation, while canopy conductance overall decreases. These patterns have often been attributed solely to the asymmetry of meteorological variables, especially in temperature and vapor pressure deficit (Zeppel et al., 2004; Bai et al., 2015; Gimenez et al., 2019; Oogathoo et al., 2020; Lin et al., 2019).

The current practice of employing g_s models that assume steady states requires iterations to converge to stable solutions at each simulation step. At the leaf level, this typically involves two nested iteration loops: first to solve the coupled photosynthesis–stomatal conductance (A_n−g_s) model for C_i and then to solve the leaf energy budget for leaf temperature (Collatz et al., 1991; Bonan et al., 2018), as g_s affects latent heat flux through transpiration. This approach can potentially lead to numerical issues (Sun et al., 2012) and increased computational costs, particularly when upscaling with complex canopies, where angular distribution setups are necessary (Wang and Frankenberg, 2022). However, by utilizing prognostic updates of variables, a dynamic model could simplify simulation steps and improve computational efficiency, enabling runs at finer temporal resolutions.

Moreover, accurate parameter estimation with steady-state models (e.g., linear fitting for the slope g₁ in empirical models) necessitates measurements to be taken after reaching each equilibrium (Leuning, 1990; Miner et al., 2017). Depending on wide-ranging stomatal response speeds (McAusland et al., 2016), obtaining one accurate response curve could take several hours (Liozon et al., 2000; Duarte et al., 2016). Too short of a time step could result in overestimation, underestimation, or unstable results of parameter estimates with steady-state assumptions (Xu and Baldocchi, 2003), which may often be overlooked (Miner et al., 2017). Alternatively, estimates can be approached with a prognostic model by fitting the entire response curve, where steady-state measurements are not fundamentally necessary.

Limitations of steady-state stomatal modeling have driven efforts to develop dynamic models, primarily at the leaf level (Damour et al., 2010; Vialet-Chabrand et al., 2017). Based on observed variations of g_s, analytical equations of sigmoidal or exponential response have been commonly used (Naumburg and Ellsworth, 2000; Noe and Giersch, 2004; Vialet-Chabrand et al., 2013, 2017; Martins et al., 2016; McAusland et al., 2016); directly adding time-dependent terms to traditional steady-state models has also been proposed (Matthews et al., 2018). While these models have demonstrated effective performance in reproducing leaf-level responses to light intensity in controlled conditions, the impacts of including temporal stomatal dynamics on the simulations of larger-scale fluxes under coupled variations in the natural environment (e.g., transpiration in the coupled diurnal cycles of radiation, temperature, and vapor pressure deficit – VPD) have not been investigated. This may be partly due to the parameterization and complexity of many models optimized for leaf-scale predictions (Kirschbaum et al., 1988; Vialet-Chabrand et al., 2016), which constrains the feasibility of scaling them to the canopy level in LSMs.

In this study, we aim to (1) implement a simplified dynamic stomatal model in the CliMA Land model, i.e., the land component of a new-generation Earth system model within the Climate Modeling Alliance (CliMA); (2) test model performance on leaf-level measurements and demonstrate an alternative method of parameter estimation with the non-steady-state model in a Bayesian nonlinear inversion framework; and (3) compare simulations of the dynamic model with traditional steady-state modeling, primarily focusing on the differences in predictions of canopy fluxes and responses to coupled environmental variations on different timescales.

2 Methods and materials

2.1 Model framework

2.1.1 Dynamic stomatal modeling

The current steady-state modeling approach in LSMs requires convergence of nested iteration loops to solve leaf fluxes at each time step (Fig. 1a) (Bonan et al., 2018). In this study, we proposed replacing the inner loop for steady solutions of the coupled photosynthesis–stomatal conductance (A_n−g_s) model with prognostic updates of g_sw at finer time steps (Fig. 1b).

Figure 1 Comparison of leaf flux calculation flows in (a) steady-state (SS) and (b) non-steady-state (NSS) dynamic modeling. Panel (a) illustrates the two nested loops at each time step in the current practice of a steady-state framework; adapted from Bonan et al. (2018). The inner iteration in the light yellow box represents the flow of solving the coupled photosynthesis–stomatal conductance (A_n−g_s) model for C_i. The outer one solves the leaf energy budget for leaf temperature (T_l). The focus of this study is to implement and compare a dynamic modeling framework to the A_n−g_s model, illustrated in the light blue box in (b), where, instead of iterating for steady solutions, g_sw is updated prognostically at finer time steps based on environmental conditions and a simplified dynamic model (Sect. 2.1.1). This NSS framework of modeling g_sw also allows prognostic updates of T_l. As its implementation is not within the scope of this study, related flows are shown in dashed parts.

At each step, instead of assuming an initial C_i and iterating until convergence, our framework starts with an initial g_sw (e.g., for the first time step of a diurnal simulation from midnight, this can be set as the minimal conductance in the dark), then solves A_n and C_i with biochemical demand and diffusive supply of internal CO₂ (Fig. 1). For instance, when applying the Farquhar photosynthesis model for C₃ plants (Farquhar et al., 1980), with a given g_sw, the RuBisCO-limited rate (A_c) and light-limited rate A_j are calculated using

$$\begin{array}{}\text{(1)}& {\displaystyle }{A}_c=V_{\mathrm{cmax}}\cdot {\displaystyle \frac{C_{\mathrm{i}}-{\mathrm{\Gamma }}^{*}}{C_{\mathrm{i}}+K_{\mathrm{m}}}}=g_{\mathrm{lc}}\cdot (C_{\mathrm{a}}-C_{\mathrm{i}})+R_{\mathrm{d}},\text{(2)}& {\displaystyle }{A}_j=J\cdot {\displaystyle \frac{C_{\mathrm{i}}-{\mathrm{\Gamma }}^{*}}{\mathrm{4}{C}_{\mathrm{i}}+\mathrm{8}{\mathrm{\Gamma }}^{*}}}=g_{\mathrm{lc}}\cdot (C_{\mathrm{a}}-C_{\mathrm{i}})+R_{\mathrm{d}},\end{array}$$

where the middle parts in Eqs. (1) and (2) represent the biochemical demand, and the right part represents the diffusive supply limitation of photosynthesis. V_cmax is the maximum carboxylation rate, C_a is the ambient CO₂ concentration, R_d is the respiration rate, Γ^* is the CO₂ compensation point with the absence of respiration, J is the electron transport rate, K_m is the Michaelis–Menten coefficient, and g_lc is the leaf total conductance to CO₂, which can be calculated using $g_{\mathrm{lc}}^{-\mathrm{1}}=g_{\mathrm{bc}}^{-\mathrm{1}}+\mathrm{1.6}{g}_{\mathrm{sw}}^{-\mathrm{1}}+g_{\mathrm{m}}^{-\mathrm{1}}$, with g_bc the boundary conductance to CO₂ and g_m the mesophyll conductance. Note that computing A_c or A_j requires solving for C_i first. With a known g_lc from g_sw at each time step, rearranging Eqs. (1) and (2) allows for the analytical solution of C_i, A_c, and A_j.

For prognostic updates of g_sw, we implemented a simplified dynamic model, adapted from previous studies on leaf-level prognostic modeling (Kirschbaum et al., 1988; Rayment et al., 2000; Noe and Giersch, 2004; Vialet-Chabrand et al., 2016):

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{\Delta }{g}_t}{\mathrm{\Delta }t}}={\displaystyle \frac{(g_{\mathrm{ss}}-g_t)}{\mathit{\tau }}},\end{array}$$

where Δt is the time step of the simulation, g_t represents the conductance at the current time step, g_ss is the target conductance calculated with steady-state models at the current conditions, and τ is the time constant, representing the timescale of stomatal responses. In this study, we used the Ball–Berry model (Ball et al., 1987) to compute the g_ss for leaf-level simulations for simplicity and the Medlyn model (Medlyn et al., 2011) for the canopy-scale simulations, as the vegetation trait dataset (De Kauwe et al., 2015) we employed for our study region is only available for the Medlyn parameters. We should note that the selection of the empirical stomatal model is of minor relevance to our primary findings, as our study focuses on the differences between steady-state and prognostic schemes only.

As indicated in the flowchart (Fig. 1), our dynamic modeling avoids nested iterations for steady solutions while requiring updates of variables at finer time steps (e.g., 5–10 min, compared to 30 or 60 min time steps of current LSMs) for the stability of simulations, which we tested and discuss in Sects. 2.3.2 and 3.3. The prognostic updates of leaf temperature can be implemented accordingly, but as it is not within the scope of this study, we prescribed the leaf temperature updates with measurements in our simulations.

2.1.2 Implementation in LSM

CliMA Land (https://github.com/CliMA/Land, last access: 23 March 2023), a new-generation LSM, is highly modularized and offers flexible model schemes (Wang et al., 2021, 2023), enabling easy implementation and assessment of the dynamic stomatal model across scales. To test the model performance and compare simulations with different stomatal modeling schemes, we implemented our non-steady-state framework in CliMA Land. More information on the CliMA Land configuration can be found in the Supplement (Sect. S1.1) and Wang et al. (2021, 2023).

2.2 Performance on leaf-level measurements

2.2.1 Leaf gas exchange

To test our model and determine key parameters, we recorded light response curves of grape (Vitis vinifera) and walnut (Juglans regia cv.) leaves using an LI-6800 portable photosynthesis system (LI-COR, Inc., Lincoln, NE, USA). Saplings of Vitis vinifera and Juglans regia cv. were planted in 5-gallon (19 L) pots with UC soil mix. 44.4 mL of Osmocote^® Smart-Release^® Plant Food Plus fertilizer was added to each pot. The plants were grown in a UC Davis lath house. The plants were watered to maintain around 75 % of completely saturated soil by weight (details in Meeker et al., 2021). The youngest, fully expanded, intact leaf was chosen and dark-adapted for 30 min. During the measurements, the photosynthetic photon flux density (PPFD) was sequentially increased following the gradient of 50, 100, 200, 400, 600, 900, 1200, 1500, and 1800 µmol m⁻² s⁻¹, with a time step of 30 min at each light level. The chamber air temperature was set at 25 °C; CO₂ partial pressure was controlled at 400 ppm, and the relative humidity in the chamber was maintained around 50 %.

2.2.2 Parameter optimization

We applied a Bayesian nonlinear inversion framework (Rodgers, 2000; Dutta et al., 2019) to jointly fit the response curves of the net photosynthetic assimilation (A_n) and stomatal conductance (g_s) for each leaf with the non-steady-state model. The forward problem in this case can be represented as follows:

$$\begin{array}{}\text{(4)}& y=\mathcal{F}(X;b)+\mathit{ϵ},\end{array}$$

where y represents the measurements, i.e., the light response curves of both A_n and g_s (see Sect. 2.2.1); ℱ represents the forward model, CliMA Land with the dynamic g_s model (see Sect. 2.1); and X is the state vector of parameters to be retrieved, which in our case includes the maximum carboxylation rate (V_cmax), the slope (g₁), and the minimum conductance (g₀) of the BB model, the mesophyll conductance (g_m) (Sun et al., 2014), and the time constant (τ). We also included a scaling factor for A_n to account for variations in the respiration rate and the ratios between CO₂ and H₂O fluxes; b is the vector of other parameters that have influences on the measurements and are known to some accuracy but not intended to be retrieved, e.g., the ratio between J_max (the maximum electron transport rate) and V_cmax, which is assumed to be 1.6 in this study but may vary across conditions (Medlyn et al., 2002), and ϵ is the error term.

The Levenberg–Marquardt (LM) iterations (Levenberg, 1944; Marquardt, 1963; Rodgers, 2000) were utilized to solve the nonlinear inversion problem and find the best estimate of key parameters:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{x}}_{i+\mathrm{1}}& ={\mathit{x}}_i+{\left((\mathrm{1}+\mathit{\gamma }){\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_i^{\mathrm{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}_i\right)}^{-\mathrm{1}}\\ & \cdot \left({\mathbf{K}}_i^{\mathbf{T}}{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}[y-F({\mathit{x}}_i\left)\right]-{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}[{\mathit{x}}_i-{\mathit{x}}_{\mathrm{a}}]\right),\end{array}\end{array}$$

where x_a is the prior estimate of the state (in this study, V_cmax: 70 µmol m⁻² s⁻¹, g₁: 9, g₀: 0.03 mol H₂O m⁻² s⁻¹, g_m: 0.4 mol CO₂ m⁻² s⁻¹, τ: 600 s, scaling_A: 1), and S_a is the prior covariance matrix, assumed to be purely diagonal, with Gaussian uncertainties in the prior state (the assumed prior standard deviation of V_cmax: 30, g₁: 3, g₀: 0.005, g_m: 0.02, τ: 100, scaling_A, 0.01). K_i is the Jacobian matrix at the ith iteration. γ is adjusted at each step, ensuring that each update of the state vector moves towards minimizing the cost function. S_ϵ is the error covariance matrix; in our case, errors were assumed to be mainly from measurement uncertainties and calculated based on the standard deviation and mean of the ΔCO₂ and ΔH₂O in LI-6800 measurements.

2.2.3 Uncertainties in traditional parameter estimation

To illustrate the influence of time steps on parameter estimation in the traditional method, which assumes steady states, we used the NSS model to generate leaf response curves to the same PPFD sequence but with different time intervals. For example, in the 5 min time step simulation, light intensity input jumped every 5 min, and measurements were assumed to be taken right before the next jump, following the traditional method. We then employed these curves to calculate the estimated g₁ and g₀ values using the traditional linear fitting method for the Ball–Berry model. The potential biases were assessed by comparing fitted parameters with different applied time steps.

2.3 Comparison of models in diurnal cycles

To assess the potential bias of the current steady-state modeling in LSMs, we compared the predictions of surface fluxes from models with different assumptions under natural environmental variations. We evaluated and compared the simulation results on both the leaf and canopy flux scales.

2.3.1 Environmental drivers and plant traits

As the light intensity tends to be the most rapidly changing environmental condition that stomata response to, we employed high-temporal-resolution radiation measurements in the field as the incoming irradiance inputs and ran CliMA Land with both setups. Photosynthetically active radiation (PAR) in a crop field (42.481677° N, 93.523521° W) was recorded with an LI-190R quantum sensor (LI-COR, Inc., Lincoln, NE, USA) at 1 s temporal resolution during August 2017.

In addition to the fluctuations of total incoming photon density that the PAR sensor can provide, canopies in natural environments also experience variations in the fraction of direct and diffuse components in the total radiation. This variation affects the distribution of PAR received by individual leaves across different layers of the canopy structure (Durand et al., 2021). To account for this effect, we employed an empirical fitting with the hourly radiation data from ERA5 (Fig. S2 in the Supplement) to estimate the partitioning between the direct and diffuse radiation (Boland et al., 2001). The empirical relationship was then applied to high-temporal-resolution PAR measurements to obtain the direct and diffuse components in the recorded total radiation, which were used as inputs for simulations at the canopy scale.

Other meteorological variables (e.g., air temperature, dew-point temperature, volumetric soil water, wind speed) were extracted from the ERA5 hourly reanalysis dataset (Hersbach et al., 2018) and input as environmental drivers for the simulations on the canopy scale. Linear interpolations were applied for runs at sub-hourly time steps. Key plant traits (e.g., V_cmax, g₁, leaf area index – LAI) were extracted from several globally gridded datasets using GriddingMachine (Wang et al., 2022; Croft et al., 2020; Butler et al., 2017; Luo et al., 2021; De Kauwe et al., 2015; Yuan et al., 2011; He et al., 2012; see also Wang et al., 2023, for detailed information on global-scale datasets used in CliMA Land).

2.3.2 Model simulations

We ran the CliMA Land surface flux simulations with different stomatal modeling schemes to assess the effects of g_s temporal response on model predictions. In the SS runs, iterations were employed to converge to steady-state solutions at each time step. For NSS mode, previous studies have observed different time constants for stomatal opening (τ_op) and closure (τ_cl) in various species, as well as a positive correlation between τ_op and τ_cl, with the ${\mathit{\tau }}_{\mathrm{op}}/{\mathit{\tau }}_{\mathrm{cl}}$ ratio varying from around $\mathrm{1}/\mathrm{3}$ to 3 (McAusland et al., 2016; Vico et al., 2011; Ozeki et al., 2022). Based on the average time constant retrieved from the leaf response curves in Sect. 2.2.1 and 2.2.2 as well as previous estimates of the time constant variations (Vialet-Chabrand et al., 2013; McAusland et al., 2016; Vialet-Chabrand et al., 2017; Vico et al., 2011; Ozeki et al., 2022), we tested several sets of τ_op and τ_cl varying from 300 to 900 s, including (a) ${\mathit{\tau }}_{\mathrm{op}}={\mathit{\tau }}_{\mathrm{cl}}=\mathrm{900}$ s, as the base comparison; (b) τ_op>τ_cl, with ${\mathit{\tau }}_{\mathrm{op}}/{\mathit{\tau }}_{\mathrm{cl}}$ ratios varying from 1.2 to 3, e.g., τ_op = 900 s, τ_cl=300 s; and (c) τ_cl>τ_op, with ${\mathit{\tau }}_{\mathrm{cl}}/{\mathit{\tau }}_{\mathrm{op}}$ ratios from 1.2 to 3, e.g., τ_op=300 s, τ_cl=900 s.

For the leaf-scale runs, we used the key parameters retrieved in previous sections and tested model predictions for an ideal clear-sky day. To investigate the differences in ecosystem fluxes, we ran and assessed the NSS and SS simulations using a time step of 1 min, with the inputs of meteorological drivers and plant traits at the location of the PAR measurement for the month of August 2017. Meteorological variables are updated at each time step. In order to further evaluate the potential contribution of g_s hysteresis to the observed diurnal hysteresis of ecosystem fluxes, we compared the standard runs with the model predictions where environmental variables (e.g., temperature, VPD, soil water content – SWC) were held constant over the daytime (as the mean of daytime values in each day). This approach allowed us to isolate the effect of hysteresis in g_s response and assess its potential contribution to the observed diurnal hysteresis of canopy and ecosystem fluxes.

Furthermore, to test the stability of prognostic modeling and assess the computational cost, we compared NSS simulations using different time steps, as well as the SS simulation run at a time step of 30 min, which is commonly used in current LSMs. This enabled us to evaluate the sensitivity of NSS predictions to the time step used, as well as to compare the computational cost for stable NSS runs and standard SS simulations. We resampled the environmental drivers from ERA5 and the PAR sensor to match the temporal resolution of the simulations, while maintaining constant average values for each diurnal cycle across simulations with different time steps.

3 Results

3.1 Model performance and parameter estimates on leaf measurements

With the parameters estimated from the LM inversion framework, the non-steady-state model predicted the temporal responses of g_sw and A_n well (Fig. 2). The model was able to capture the gradual increases in g_sw and A_n after each step change in APAR, and the reproduced curves were close to the measurements, with all R² higher than 0.98. Fitted time constant τ showed a variation between the two example leaves (292 and 2028 s for the Vitis vinifera leaf and the Juglans regia cv. leaf, respectively). The relative difference in the time constant matched the variations of response speed observed in the measured response curves (Fig. 2). Compared to the SS model, the dynamic model provided more accurate prediction of the temporal responses. The improvements in R² were more prominent in the predictions of the Juglans regia cv. leaf responses, which have a larger time constant, than in those of the Vitis vinifera leaf. Other parameters estimated for the Vitis vinifera leaf include V_cmax of 71 µmol m⁻² s⁻¹, g₁ of 11.3, g₀ of 0.023 mol H₂O m⁻² s⁻¹, g_m of 0.18 mol CO₂ m⁻² s⁻¹, and scaling_A of 1.1. For the Juglans regia cv. leaf, V_cmax is 152 µmol m⁻² s⁻¹, g₁ is 3.9, g₀ is 0.052 mol H₂O m⁻² s⁻¹, g_m is 0.34 mol CO₂ m⁻² s⁻¹, and scaling_A is 1.0.

Figure 2 Modeled and measured temporal responses of the stomatal conductance (g_sw) and net photosynthesis rate (A_n) to the step changes in APAR for different leaves. The shaded area indicates the difference between the prediction of the steady-state (SS) model and the non-steady-state (NSS) dynamic model. (a–b) The temporal responses of the Vitis vinifera leaf and (c–d) the Juglans regia cv. leaf.

The dynamic model was also able to better capture the temporal variations of internal CO₂ concentration (Fig. 3). Particularly, the NSS model reproduced the undershooting of the intercellular CO₂ concentration (C_i) after each step change in light intensity, which resulted from the differences in the speed of g_s and A_n responses and their interactions. As shown in the measured time series (Fig. 2), after each increase in the incident light, photosynthesis was able to respond almost instantaneously, leading to a rapid decrease in C_i, while stomata opened gradually, slowly bringing up C_i over time. This then led to a gradual rise of A_n after the initial rapid response, indicating the regulation of g_s on A_n through its impacts on the internal CO₂ supply. In the meantime, the increasing A_n further promoted the opening of stomata with a higher internal CO₂ demand, demonstrating their coupled responses to environmental variations.

Figure 3 Modeled and measured temporal responses of intercellular CO₂ concentration (C_i) for (a) the Vitis vinifera leaf and (b) the Juglans regia cv. leaf. As indicated by the labels, measured curves were shifted to illustrate the comparison of modeled and measured response patterns, as the absolute values are not directly comparable due to different assumptions of LI-6800 and CliMA Land in calculating the internal CO₂. The shaded area indicates the difference between the prediction of the steady-state (SS) model and the non-steady-state (NSS) dynamic model.

Figure 4 Parameter estimates for the Ball–Berry model with the traditional linear fitting method using model-reproduced response curves with different time steps (5, 10, 30, 60, 120 min). (a) Fitting results for the Vitis vinifera leaf and (b) the Juglans regia cv. leaf. Corresponding Ball–Berry index and g_sw are plotted, along with the fitted lines and parameters (i.e., the Ball–Berry slope, g₁, and the intercept, namely, the minimum conductance, g₀).

With the dynamic model and optimized parameters that accurately reproduced the measured leaf responses, we investigated the influence of time steps used in light response curves (i.e., the length of intervals between step changes in light intensity) on parameter estimates obtained with traditional methods (Fig. 4). The results showed that, particularly for the Juglans regia cv. leaf that has a long time constant over 2000 s, the values and relationship between the Ball–Berry index and g_sw varied significantly depending on the time step used, resulting in notable uncertainties in fitted g₁ and g₀ with too short of a interval to reach equilibrium. This also suggested that obtaining reliable estimates for this leaf using the traditional method could require more than an hour for stable readings at each step.

3.2 Model comparison in diurnal cycles

3.2.1 Leaf responses

To compare NSS and SS models over the course of a day, we evaluated the differences in their predictions of leaf responses to an ideal diurnal cycle of light with other environmental conditions (e.g., temperature, VPD, CO₂) held constant (Fig. 5). Results showed that compared to NSS, the SS model predicted a higher A_n and g_s in the morning, as it assumed the stomata could respond immediately to an increase in light, while in the more realistic NSS simulation, the gradual opening of stomata limited the CO₂ supply for photosynthesis with a lower C_i. The opposite was true for the afternoon, but the overestimation of A_n and g_s in SS modeling in the morning was more significant than the underestimation in the afternoon, leading to slightly higher diurnally integrated predictions than those of the NSS model. This was due to the fact that in the course of sunset, the major limiting factor on productivity was the decreasing light, in contrast to the sunrise where it was the available C_i regulated by g_s responses that mainly constrained A_n increases. The relative differences (RDs) in integrated g_s in the morning and afternoon were both higher than those of the photosynthesis, reflecting the differences in the response speed.

The differences in predictions of A_n and g_sw responses also led to RDs in the intrinsic water use efficiency (WUE, i.e., the ratio between A_n and g_sw). Although the mean instantaneous WUE during the daytime was higher in the NSS simulation, diurnal WUE calculated from the integrated A_n and g_sw was lower. This was because the gradual opening of stomata during sunrise limited assimilation in the morning, whereas during sunset, delayed closure led to unnecessary water loss when carbon gain was constrained by low light.

Figure 5 Predictions of the leaf diurnal course of (a) net photosynthesis rate (A_n), (b) stomatal conductance to water vapor (g_sw), (c) intercellular CO₂ concentration (C_i), and (d) intrinsic water use efficiency (WUE) for a leaf with a stomatal time constant of 900 s on an ideal clear-sky day (${\mathit{\tau }}_{\mathrm{op}}={\mathit{\tau }}_{\mathrm{cl}}=\mathrm{900}$ s). Other environmental conditions (e.g., leaf temperature, VPD) were held constant. The shaded areas indicate the differences between the NSS and SS simulations (blue: SS > NSS; orange: SS < NSS) in both absolute and relative terms. Relative differences (RD, NSS − SS) in the temporal integrals are also presented for morning (05:00–12:00 local solar time), afternoon (12:00–19:00), and daytime (05:00–19:00). RD^* of WUE represents the ratio between integrated A_n and g_sw, differing from the RD, the integral of the instantaneous WUE during the daytime.

3.2.2 Canopy fluxes

To quantify the impacts of the inclusion of g_s temporal response, we analyzed the simulated canopy fluxes under natural radiation variations and coupled dynamics of environmental conditions. As shown in the examples of diurnal cycle simulations (Figs. 6 and 7), the SS model predicted higher variations in instantaneous fluxes in response to rapid fluctuations in radiation, particularly in transpiration rates.

Figure 6 Comparison of the predicted diurnal cycles of ecosystem fluxes, as well as (a) the net productivity and (b) the transpiration rate for DOY 219 in 2017. The shaded areas in blue and orange indicate the differences between the NSS and SS predictions (blue: SS > NSS; orange: SS < NSS). The shaded areas in gray under the radiation curves represent the diffuse component of the total radiation. Relative differences (RD, NSS − SS) in the temporal integrated fluxes are also presented for morning (05:00–12:00), afternoon (12:00–19:00), and daytime (05:00–19:00).

Figure 7 Comparison of the predicted diurnal cycles of ecosystem fluxes, as well as (a) the net productivity and (b) the transpiration rate for DOY 224 in 2017. The shaded areas in blue and orange indicate the differences between the NSS and SS predictions (blue: SS > NSS; orange: SS < NSS). The shaded areas in gray under the radiation curves represent the diffuse component of the total radiation. Relative differences (RD, NSS − SS) in the temporal integrated fluxes are also presented for morning (05:00–12:00), afternoon (12:00–19:00), and daytime (05:00–19:00).

When τ_op=τ_cl, the overestimate of morning g_sw in SS predictions was mostly compensated for by an underestimate of afternoon g_sw, resulting in relatively minor differences in daily average g_sw (Fig. 8c; the mean RD of daytime mean g_sw was 0.5 %). However, when the time constants of stomatal opening and closure were not equal, there was overall underestimation or overestimation in both mornings and afternoons, and the daily mean RDs could be notable (Figs. 8g and k, S3, S4). For example, when τ_op=3τ_cl, the faster closure than opening of stomata led to overall lower conductance over the diurnal cycles compared to SS runs (Fig. 8k; the mean RD of daytime mean g_sw was −6.1 %).

Figure 8 Relative differences (NSS − SS; RD) in the predicted daytime mean fluxes of the NSS (time step: 1 min) and SS (1 min) simulations for August 2017. The solid line in each box indicates the median, and the dashed line represents the mean. The results for the transpiration rate (H₂O flux), net productivity (CO₂ flux), canopy-averaged stomatal conductance to water (g_sw), and water use efficiency (WUE) are shown in the respective columns from left to right. (a–d) ${\mathit{\tau }}_{\mathrm{op}}={\mathit{\tau }}_{\mathrm{cl}}=\mathrm{900}$ s, (e–h) τ_op=300 s, τ_cl=900 s, (i–l) τ_op=900 s, τ_cl=300 s. Diurnal: 05:00–19:00, morning: 05:00–12:00, afternoon: 12:00–19:00.

In the simulations with the same time constants of stomata opening and closure, the differences in fluxes between the NSS and SS predictions were not significant when integrated over monthly periods (e.g., the mean RD of transpiration in August 2017 was 0.87 %, and the median was 1.0 %) but could be notable at sub-diurnal scales depending on the environmental conditions (e.g., the variation of afternoon RDs ranged from −7.4 % to 6.1 %). When there were differences in τ_op and τ_cl, the divergences between NSS and SS predictions could be more significant (e.g., when ${\mathit{\tau }}_{\mathrm{op}}=\mathrm{1}/\mathrm{3}{\mathit{\tau }}_{\mathrm{cl}}$, the mean RD of transpiration in August 2017 was 4.9 %, and the maximum daily mean RD of transpiration was 9.0 %).

The overall tendency to overestimate productivity with traditional SS models was also observed at the canopy scale, when τ_op was equal to or larger than τ_cl, as the regulation of g_s hysteresis on the supply of CO₂ for photosynthesis was not considered (Fig. 8b and j). For example, in Fig. 6, when rapid spikes of radiation occurred in the afternoon, the speed of g_s response constrained the increases in photosynthesis in the NSS simulation. However, when τ_op was smaller than τ_cl, predicted daily mean photosynthesis was slightly higher in the NSS simulation (Fig. 8f; the mean RD of productivity in August 2017 was 0.2 %). This resulted from the overall higher g_sw due to faster opening than closure (Fig. 8g), as higher conductance resulted in higher C_i (Fig. S5), leading to generally higher rates of photosynthesis.

In contrast to the leaf-scale results, when accounting for other covarying environmental drivers (e.g., temperature, VPD, soil water content), the SS model tended to underestimate canopy transpiration rates when τ_op=τ_cl (Figs. 6b, 7b, 8a). This could be because the transpiration rates were determined by both g_sw and VPD. During the daytime, VPD typically increased following air temperature and peaked in the afternoon, when the slow response of stomata to the increasing VPD and decreasing radiation could result in excess water loss (Figs. 7b, 8a and c). The overestimation of productivity and underestimation of transpiration in SS simulations both contributed to the overestimation of WUE. When τ_op<τ_cl, slower stomatal closure led to increased water loss and thus a more significant underestimation of transpiration in the SS predictions (Fig. 8e; the mean RD of transpiration in August 2017 was 4.9 %), resulting in further overestimation of the WUE. When τ_op>τ_cl, the NSS model predicted lower transpiration rates and higher WUE compared to the SS model because of the overall lower g_sw (Fig. 8i and k). As lower g_sw yielded overall lower C_i (Fig. S5), larger gradients of CO₂ concentration across stomata contributed to higher WUE on daily and monthly timescales.

3.2.3 Diurnal hysteresis

To investigate the relative contributions of g_s hysteresis and environmental variables to the hysteresis observed in plant behaviors and ecosystem fluxes, we separated the effects of these two factors by comparing predicted response curves in NSS and SS simulations with and without diurnal environmental variations (e.g., temperature, VPD, soil water content). While the asymmetry of environmental variables in the diurnal cycle could lead to a modeled hysteresis of g_s in response to radiation, where g_s tended to be lower in the afternoon mainly due to higher VPD and temperature, our results (Fig. 9) showed that the kinetic lag of g_s could partially offset this effect (Fig. 9b and d), even presenting an opposite tendency at low radiation. Additionally, only the NSS model simulations predicted a hysteresis of canopy transpiration, with or without the consideration of coupled environmental variations (Fig. 9g and h), in which canopy H₂O fluxes tended to be higher in the afternoon. Differences between τ_op and τ_cl affected the magnitudes of hysteresis, but the overall patterns remained similar (Figs. S6, S7).

Figure 9 Hysteresis of the canopy mean stomatal conductance (g_sw) and canopy transpiration rate (T) in response to radiation during an ideal clear-sky day, when ${\mathit{\tau }}_{\mathrm{op}}={\mathit{\tau }}_{\mathrm{cl}}=\mathrm{900}$ s. (a, e) SS model and (b, f) SS model with coupled diurnal variations of environmental conditions (Env, e.g., air temperature, VPD). (c, g) NSS model and (d, h) NSS model with Env. (a–d) Normalized g_sw responses and (e–h) normalized T responses. In simulations without Env variations, except for the radiation, all the other environmental drivers were kept at the daytime means. g_sw and T are normalized with the values at noon (12:00). Arrows indicate the increasing and decreasing parts of the diurnal courses.

3.3 Stability of the dynamic model

We further assessed the sensitivity of the dynamic modeling to the time step of simulation. Figure 10 shows that the NSS model was be able to run at a time step of 10 min stably and still demonstrated the impacts of gradual g_s responses compared to the traditional practice of SS modeling at a time step of 30 min.

Figure 10 Simulations of dynamic g_sw using different time steps (2, 6, 10 min) and comparison with the traditional steady-state modeling (30 min resolution) predictions. (a) Using high-temporal-resolution PAR as radiation input, values are resampled accordingly to match the time step used; (b) using ERA5 hourly radiation as input, values are linearly interpolated to 30 min resolution. The shaded areas in gray under the radiation curves represent the diffuse component of the total radiation.

4 Discussion

In this study, we demonstrated the feasibility and benefits of implementing a non-steady-state stomatal conductance modeling framework from the leaf to canopy scale in a new-generation LSM, CliMA Land. Our results suggested that compared to traditional steady-state models, the dynamic model was able to provide more realistic and accurate predictions of leaf temporal responses to the changes in light intensity (Sect. 3.1). In the meantime, modeling g_s with prognostic updates – similar to how plants control their stomata movements gradually in natural environments – increased neither computational cost nor model complexity, as simulations were simplified with iterations to solve for steady states avoided. Sun et al. (2012) pointed out that the default three-step fix-point iteration for solving the coupled A_n−g_s model in CLM4 (the Community Land Model version 4) does not always converge, leading to uncertainties in flux predictions. In our simulations at the canopy scale (Sect. 3.3), the dynamic model could be stably run at a temporal resolution (10 min) that presented comparable efficiency to the current practice of 30 min resolution SS simulations commonly used in LSMs (three-step prognostic updates for each SS default three-step iteration). This also indicates that the dynamic model can enable predictions of canopy flux dynamics at a finer time resolution with higher efficiency and accuracy.

With the non-steady-state model, we were able to apply a Bayesian nonlinear inversion framework to jointly fit the light response curves of both A_n and g_s and obtain estimates for key parameters (Sect. 3.1). As suggested in our results (Fig. 4) and previous studies (Xu and Baldocchi, 2003; Miner et al., 2017), the time step of light response curves can notably influence the estimated parameters obtained from the traditional linear fitting method for steady-state empirical models. Our framework with the dynamic model can help reduce the time required for accurate parameter estimations, particularly for leaves with long time constants, as equilibrium is not required. Although the retrieval setups presented in this study may not be optimal for estimating V_cmax, which is typically derived from A−C_i response curves (Medlyn et al., 2002; Miao et al., 2009; Duarte et al., 2016), a similar framework can be applied to other scenarios for estimation of various parameters, including an A−C_i curve for V_cmax.

Furthermore, we evaluated how the inclusion of g_s temporal responses could affect model predictions of leaf and canopy fluxes in diurnal cycles with natural environmental variations (Sect. 3.2.2). The comparison of NSS and SS simulations indicated that the differences in fluxes depended on the integration timescales, relative speed of stomatal opening and closure, and environmental variations. In terms of instantaneous effects, slow opening of stomata tended to limit productivity responses to rapid radiation increases, and delayed closure of g_s following decreases in radiation or increases in environmental stress (e.g., increasing VPD) resulted in unnecessary water loss. The divergence of NSS and SS schemes was less significant when considering the monthly integrated canopy fluxes compared to daily or sub-diurnal-scale results. The monthly differences were more notable when the speeds of stomata opening and closure differed. The overall effects on WUE also depended on the relative speed of opening and closure. In the simulations where stomata open at a speed similar to or faster than they close, excessive water loss in the afternoons, when VPD was high, led to a lower WUE. This also suggested that traditional steady-state simulations may overestimate WUE. Similar impacts have been noted in studies on leaf-scale response to PPFD fluctuations (Lawson et al., 2011; Lawson and Blatt, 2014; McAusland et al., 2016). Meanwhile, when stomata opened more slowly than they closed, plants exhibited both a lower maximum g_s during diurnal cycles and a lower average g_s compared to the SS runs. This resulted in reduced transpiration and increased WUE, even though productivity was also suppressed. These results suggest that the temporal hysteresis of g_s can have impacts on integrated canopy fluxes, and further studies on variations of stomata opening and closure speeds across plants can be helpful to assess these effects more comprehensively on larger scales.

In addition, the hysteresis of leaf-level g_s response can contribute to the hysteresis patterns at the ecosystem scale, which have often been solely attributed to the asymmetry of environmental variables during the daytime. For instance, higher evapotranspiration (ET) fluxes and sap velocity (i.e., an indicator of plant transpiration rate) have been observed in the field, with the explanation often focused on higher VPD in the afternoon following increased air temperature (Zeppel et al., 2004; Gimenez et al., 2019; Oogathoo et al., 2020; Lin et al., 2019). Our simulations showed that the SS model with diurnal environmental variations was unable to reproduce this hysteresis pattern, while it was captured in NSS runs. This indicated the significance of considering g_s temporal dynamics when interpreting diurnal hysteresis in transpiration (Sect. 3.2.3). Moreover, observed patterns of lower g_s in the afternoon have also been commonly explained by similar environmental asymmetry (Bai et al., 2015; Lin et al., 2019), while our results suggested the kinetic lag of g_s could partially offset this effect and should thus be taken into account in understanding the hysteresis patterns.

Our study mainly focused on taking the first step to implement prognostic stomatal modeling in an LSM, including the impacts on canopy flux simulations. Further improvements can be made in assessing other effects of g_s temporal responses in LSM projections, as well as validating the comparisons with site-level observations. For example, while daily effects on canopy productivity were minor, they may add up to significant differences in long-term vegetation growth trajectories. As plant traits were prescribed in our simulations, the accumulative effects were not included in our analysis of the short-term predictions. The transient limitation on photosynthesis from the slow temporal response of g_s could also cause potential photoinhibitory damage to the photosystem II reaction centers. Future studies can focus on the parameterization of these impacts in LSMs and the evaluation of cumulative effects on plant growth and hydraulics in the long term.

The dynamic g_s model enables predictions of temporal changes in latent heat flux through transpiration in leaf energy balance, which allows a similar prognostic framework to be employed for the modeling of leaf temperature. Coupled dynamic modeling of stomatal conductance and leaf temperature will enhance our ability to evaluate the influences of g_s hysteresis on the feedback between leaf transpiration and thermal condition. This is out of the scope of our current study but can be a valuable direction for future research efforts. Bonan et al. (2018) implemented a non-steady-state framework for leaf temperature modeling, but as steady-state g_s models were employed, iterations for stable solutions were still required. With the dynamic g_s model presented in this study, the traditional nested iteration loops in leaf flux calculations, which can take up to 40 iterations to solve for a single simulation step in CLM4.5 (Bonan et al., 2018), can be replaced by more efficient and accurate prognostic updates of variables with ordinary differential equations (ODEs). Such an approach can also facilitate better couplings of LSMs with other components in Earth system models (ESMs), where ODE systems are commonly used.

5 Conclusions

We implemented a simplified dynamic stomatal conductance model in CliMA Land and evaluated its impacts on model simulations across scales. In comparison with the traditional steady-state model, the dynamic model better predicted the coupled temporal responses of A_n, g_s, and C_i observed in leaf measurements. We also found uncertainties in parameter estimation for steady-state g_s models with the traditional linear fitting method, when too short of a time step used resulted in unstable estimates. We proposed an alternative approach using a Bayesian nonlinear inversion framework with a dynamic model, which could help reduce the time investment for estimation, particularly for leaves with long time constants. Our results on canopy-scale simulations suggested that the effects of temporal g_s responses on ecosystem fluxes depend on the timescales of integration, relative speed of stomatal opening and closure, and environmental variations. Although the differences in monthly integrated fluxes between NSS and SS simulations were relatively minor, the hysteresis of g_s should be taken into account when predicting diurnal courses and quantifying sub-diurnal-scale fluxes, as well as explaining the hysteresis patterns observed in diurnal cycles. In addition, we also show that these divergences become notable when stomata open and close at different speeds.

We demonstrated that the more realistic prognostic modeling of gradual g_s response simplified the simulation as iteration loops for solving steady states at each time step were avoided, and the dynamic model can be run at a finer time resolution that presents comparable computational costs to the current practice of steady-state leaf flux calculation. A similar framework can be extended to leaf temperature modeling, which will enable prognostic updates of leaf-level variables with higher efficiency and accuracy, towards better couplings of LSMs with other components in Earth system models (ESMs).