Improved model predictions of carbon and water fluxes by including drought legacy effects

Besides simultaneous influences, droughts have lasting impacts on vegetation by impairing hydraulic and photosynthetic capacities, known as the drought legacy effects. The ignorance of legacy effects in numerical simulations, such as lagged xylem recovery, may lead to significant model-observation discrepancies. However, the limited temporal resolution of most observational data makes it challenging to capture the physiological dynamics necessary to improve model accuracy. Here, we investigated the recovery of carbon flux (represented by gross primary productivity, GPP) and water flux (represented by evapotranspiration, ET) following a severe drought in 2012, using half-hourly eddy-covariance flux observations and weekly predawn leaf water potential measurements from a temperate forest in the Central US. We implemented both optimality-based and empirical stomatal models within a land surface model, testing three drought recovery scenarios for each: no recovery, full recovery, and partial recovery of xylem hydraulic conductance and photosynthetic capacity. Before and during the drought, all stomatal models performed similarly for GPP and ET. Post-drought, assuming no recovery led to underestimated ET; assuming full recovery led to overestimated GPP; and assuming partial recovery improved both, indicating persistent biochemical limitations after drought. The observed carbon-water decoupling during and after the event further points to non-stomatal constraints on photosynthesis and unequal stress on carbon and water fluxes. Our work highlights the need to account for delayed recovery of xylem hydraulics and photosynthetic capacity when modeling drought legacy effects. Further research to mechanistically represent dynamic recovery processes, particularly their duration and magnitude, is essential for improving the modeling of global carbon and water fluxes.

1 Introduction

Plants are exposed to increasing frequency and severity of extreme events under climate change, such as heatwaves and droughts (Piao et al., 2019; Xu et al., 2019). During droughts, reduced precipitation and higher evaporation demand driven by rising air temperature usually lead to further depletion of soil moisture (Bastos et al., 2020; Wolf et al., 2016). Although plants partially close their stomata to protect their hydraulic systems from excessive damage (Sperry et al., 2017; Wang et al., 2020), they often suffer irreversible hydraulic impairment due to xylem cavitation. As a result, even after meteorological drought conditions are alleviated, persistent hydraulic damage limits water transport to leaves and a downregulation of photosynthesis capacity may persist (Müller and Bahn, 2022). The reduced photosynthetic capacity during and after drought conditions, caused by insufficient water supply, further delays recovery to pre-drought functionality. This lagged recovery, often referred to as the drought legacy effect, has been documented in drought-induced tree mortality (Matusick et al., 2018), carbon fluxes (Haberstroh et al., 2025), and tree ring growth (Anderegg et al., 2015). Vegetation models, however, often fail to accurately capture drought responses and legacy effects (Anderegg et al., 2015), as they typically oversimplify plant hydraulic processes and stomatal behavior, neglecting the persistent damage to hydraulic and photosynthetic systems. Quantifying the magnitude and timing of hydraulic impairment and recovery is therefore essential particularly for local scale modeling and important also for large-scale land surface and Earth system models to better represent the global carbon cycle.

Typical model simulations simplify the complexity and diversity of ecosystems by using plant functional type (PFT)-specific parameters to represent ecosystem responses. These models are often constrained by ecosystem-level measurements, which reflect the average of potentially divergent responses within the system. To simulate the effects of drought on CO₂ assimilation and transpiration loss, models often use a scaling factor (β) ranging from 0 to 1, to attenuate photosynthesis and/or stomatal conductance. Empirical stomatal models like the Ball-Berry model (Ball et al., 1987), Medlyn model (Medlyn et al., 2011), and Leuning model (Leuning, 1995), typically include a semi-empirical stomatal slope parameter g1, which reflects the efficiency of carbon gain relative to water cost. The accuracy of these models largely depends on the representation of water stress impacts on photosynthesis. Drought responses in photosynthesis are driven by both stomatal and non-stomatal limitations as shown through observations and model studies (Drake et al., 2017; Egea et al., 2011). However, the relative importance of these two components varies across studies (Gourlez de la Motte et al., 2020; Keenan et al., 2009, 2010; Reichstein et al., 2003; Zhou et al., 2015, 2013), and depends on the intensity and duration of water stress (Niinemets and Keenan, 2014). These limitations can be incorporated into models by either reducing the well-watered maximum carboxylation rate (V_cmax) to account for biochemical constraints (Kennedy et al., 2019) or adjusting the stomatal slope parameter (g1) to represent biophysical limitations (De Kauwe et al., 2015).

After drought events, the scaling factor β is often assumed to return to 1 as soil water refills (Fang and Gentine, 2024). Many terrestrial ecosystem models presume that this recovery is both immediate and complete, leading to the full restoration of carboxylation efficiency, stomatal conductance, and associated carbon and water fluxes. However, this assumption is likely invalid, as many studies have documented persistent drought legacy effects lasting from subseasonal to multi-year time scales (Anderegg et al., 2015; Kannenberg et al., 2020; Wu et al., 2018). These effects have been observed through carbon flux measurements (Yu et al., 2022), tree ring width analysis (Anderegg et al., 2015), and vegetation greenness indices (Wu et al., 2018) using both field observations and satellite data. This evidence highlights a critical limitation in the full recovery assumption commonly employed in models. Despite advancements in our understanding of plant physiological responses to drought, there remains a substantial knowledge gap regarding the key factors driving drought legacy effects and their underlying mechanisms. Addressing this gap is essential for improving the accuracy of ecosystem models in representing post-drought recovery processes.

A potential explanation for drought legacy effects is the incomplete recovery of hydraulic conductance, even after soil moisture has refilled, which limits upward water transport to the canopy (Kannenberg et al., 2020). This residual stress can propagate into biochemistry (V_cmax) and impose limitations on intrinsic water use efficiency through g1. However, the incomplete recovery of hydraulic conductance, V_cmax, or g1 has not been explicitly incorporated into most process-based models. Traditionally, water stress in models is represented as a nonlinear function of soil moisture content. Recent advancements in modeling the soil-plant-atmosphere continuum have improved the representation of plant water stress (Eller et al., 2020; Kennedy et al., 2019; Wang et al., 2021; Xu et al., 2016; Yao et al., 2022). These advancements incorporate vertically-resolved water pressure profiles and hydraulic traits through vulnerability curves, enabling the simulation of dynamic hydraulic conductance over time. In this context, imposing limits on hydraulic conductance offers a practical way to represent the incomplete recovery of hydraulic function after soil moisture refilling. These limits reflect the “stress history” experienced by plants, which can be inferred from observational data. Although drought legacy effects have been observed across a variety of metrics and time scales (Anderegg et al., 2015; Yu et al., 2022), higher temporal resolution data often provides better opportunities for model-data assimilation compared to coarser datasets. By integrating high-resolution observational carbon and water flux data with improved modeling frameworks, we can better understand the dominant controls on ecosystem-level drought legacy effects and their long-term impacts on ecosystem functioning.

In the summer of 2012, the Central US experienced a severe drought that rapidly developed and intensified over less than two months (Basara et al., 2019). This study focuses on the Missouri Ozark AmeriFlux (MOFLUX) eddy-covariance site located at the epicenter of the drought, providing an ideal setting to investigate drought effects. Using a land surface model (Climate Modeling Alliance, CliMA, Wang et al., 2021) and flux observations from the MOFLUX site, we aim to test the model's performance by incorporating drought legacy effects assimilated from observations. Our study was driven by three main objectives:

1. simulating the drought response using the CliMA land surface model.

2. identifying the relative strength of biochemical and stomatal factors in limiting carbon and water fluxes after the drought.

3. diagnostically integrating legacy effects to test model performance under different recovery scenarios.

To assess the generalizability of the modeling framework, we extended the analysis by conducting simulations for five additional sites affected by the 2012 drought. The findings of our research will contribute to a better understanding of recovery processes and the factors driving drought legacy effects, enabling more accurate identification of where and how to apply stress factors in models.

2 Materials and methods

2.1 US-MOz site

The MOFLUX eddy-covariance flux tower (AmeriFlux ID: US-MOz) is located in the Central US (latitude: 38.7441, longitude: −92.2000), and characterized by a warm, humid and continental climate zone. The mean annual temperature is 13 °C and the mean annual precipitation is 1,140 mm. Important tree species in the area are white oak (Quercus alba), black oak (Q. velutina), sugar maple (Acer saccharum), shagbark hickory (Carya ovata), and eastern redcedar (Juniperus virginiana) (Gu et al., 2015). We used half-hourly measurements of evapotranspiration (ET), as well as net ecosystem carbon dioxide exchange (NEE) from which gross primary productivity (GPP) was inferred using the night-time partitioning method. The processing of carbon and water fluxes data used here follows the standardized AmeriFlux pipelines (Gu et al., 2016). Site-level leaf area index (LAI) was extracted using the GriddingMachine package (Wang et al., 2022) from MODIS (MODerate resolution Imaging Spectroradiometer) LAI products (Yuan et al., 2011). Predawn leaf water potential (Ψ_leaf,pd) measurements were initiated in 2004 and were carried out weekly to bi-weekly.

The 2012 drought is one of the most extensive drought events over the Central US on record, with rapid development and pronounced area encompassed by persistent dryness and warm temperatures. Specifically, the drought peaked in July, when nearly 64 % of the area experienced at least moderate drought. Wood et al. (2023) divided the 2012 growing season into four stages: peak flux (DOY 125–150), flash stress (DOY 160–180), wilted (DOY 185–240), and recovery (DOY 255–270). We slightly modified these period definitions to eliminate any gaps between them, defining the stages as follows: peak flux (DOY 125–150), flash stress (DOY 150–180), wilted (DOY 180–240), and recovery (DOY 240–280).

2.2 Other flux sites used for model validation

We conducted model simulations at several flux-tower sites that experienced the 2012 drought, including US-MMS, US-GLE, US-UMB, US-NR1, and US-UMd. The US-MMS site is located in Morgan-Monroe State Forest, Indiana, within a temperate deciduous broadleaf forest dominated by a maple–beech–oak–hickory transition community representing secondary successional stands approximately 60–80 years old (Novick and Phillips, 2022). The US-GLE site is located in northern Michigan and consists of a mixed forest composed of aspen, maple, pine, and northern hardwood species (Frank et al., 2022). The US-UMB site, also located in northern lower Michigan, is an aspen-dominated northern hardwood forest (Classen et al., 2023). The US-NR1 site is located in the Rocky Mountains of Colorado within a subalpine evergreen needleleaf forest dominated by Abies lasiocarpa, Picea engelmannii, and Pinus contorta (Blanken et al., 2026). The US-UMd site, located in northern lower Michigan, represents a temperate deciduous broadleaf forest (Gough et al., 2024).

2.3 Vegetation module of CliMA Land

We used the vegetation module of the land surface model developed within the Climate Modeling Alliance (CliMA) (Wang et al., 2021, 2023). The model's performance has been validated in terms of carbon and water fluxes, as well as solar-induced chlorophyll fluorescence (Braghiere et al., 2021; Wang et al., 2021, 2023). The model adopts a modular structure, allowing flexible combinations of radiative transfer schemes (broadband or hyperspectral), photosynthesis models, stomatal models, and plant hydraulic representations.

We ran CliMA Land in hyperspectral mode for radiative transfer with a vertically resolved canopy, each layer of which is further split into sunlit and shaded fractions based on uniform leaf inclination and azimuthal distributions (Wang and Frankenberg, 2022; Zhao et al., 2019). At each time step, sunlit fractions and absorbed photosynthetically active radiation (APAR) are calculated, and leaf-level photosynthesis is solved by coupling the diffusion model with the Farquhar model (Farquhar et al., 1980) (Eq. 1). LAI is prescribed from 8 d MODIS observations and linearly interpolated to daily values, thereby capturing seasonal phenology and structural changes during the 2012 drought (Fig. S1 in the Supplement). Soil respiration is modeled using the Arrhenius equation with parameterized temperature sensitivity of soil respiration.

$$\begin{array}{}\text{(1)}& A_{\mathrm{net}}=V_{\mathrm{cmax}}{\displaystyle \frac{C_{\mathrm{i}}-{\mathrm{\Gamma }}^{\ast }}{C_{\mathrm{i}}+K_{\mathrm{m}}}}-R_{\mathrm{d}}=g_{\mathrm{s}}\cdot (C_{\mathrm{a}}-C_{\mathrm{i}})\end{array}$$

where V_cmax is the maximum rate of Rubisco activity and K_m is the Michaelis-Menten coefficient for Rubisco kinetics. Γ^∗ is the CO₂ compensation point in the absence of dark respiration and R_d is the dark respiration rate. C_a is ambient CO₂ concentration. C_i is the intercellular CO₂ concentration. The model details of stomatal conductance (g_s) can be found in Sect. 2.4.

Canopy transpiration is computed using the diffusion model as Eq. (2).

$$\begin{array}{}\text{(2)}& T=\sum _{i=\mathrm{1}}^{n_{\mathrm{canopy}}}\mathrm{1.6}\cdot g_{\mathrm{s}}\cdot D\cdot {\mathrm{LAI}}_i\end{array}$$

LAI_i is the leaf area index of ith layer. D is leaf-to-air vapour pressure deficit. 1.6 is the unitless relative diffusivity of water vapor with respect to CO₂.

The hydraulic system is represented as a vertically solved soil-plant-atmosphere continuum, including two root layers and 19 canopy layers. A fixed root : stem : leaf resistance ratio of $\mathrm{2}:\mathrm{1}:\mathrm{1}$ is assumed (Sperry et al., 2017). Leaf water potential (Ψ_leaf) is updated following Darcy's law based on vertical water flow and hydraulic conductance (Eq. 3).

$$\begin{array}{}\text{(3)}& {\mathit{\psi }}_{\mathrm{leaf}}={\mathit{\psi }}_{\mathrm{soil}}-F\times \sum _{i\in \mathrm{leaf},\mathrm{stem},\mathrm{root}}{\displaystyle \frac{\mathrm{1}}{K_i}}\end{array}$$

F is transpiration-driven vertical water flow and K_i is hydraulic conductance of each segment. Ψ_soil is soil water potential.

Wood capacitance and water storage dynamics were not included. Our model does not simulate soil evaporation and re-evaporation of rainfall interception, so we used a fixed $T/\text{ET}$ ratio (0.67, falling in range reported in Fatichi and Pappas, 2017) to link the modeled plant transpiration to ET observation. We also tested varying T/ET ratio with observation-based partitioned transpiration (Zhou et al., 2016).

2.4 Stomatal conductance models and water stress formulation

The CliMA Land model incorporates several stomatal optimality models (e.g., Eller et al., 2020; Sperry et al., 2017; Wang et al., 2020) and three empirical stomatal models: Medlyn (Medlyn et al., 2011), Ball-Berry (Ball et al., 1987) and Leuning models (Leuning, 1995). In this study, we used the Medlyn stomatal model and an optimality stomatal model (Wang et al., 2021). In empirical formulations, drought responses are introduced through a water stress factor β. Two pathways of stress propagation are considered. The first one operates by adjusting parameters of the empirical stomatal models, such as the slope parameter (g1) (Eq. 4). The second pathway involves the modulation of photosynthetic capacity through down-regulation of V_cmax (Eq. 5), which also reduces stomatal conductance as photosynthesis is reduced. A schematic is shown in Fig. 1.

$$\begin{array}{}\text{(4)}& {\displaystyle }{g}_{\mathrm{s}}=g_{\mathrm{0}}+\mathit{\beta }\cdot \left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}}}{\sqrt{D}}}\right)\cdot {\displaystyle \frac{A_{\mathrm{net}}}{C_{\mathrm{a}}}}\text{(5)}& {\displaystyle }{g}_{\mathrm{s}}=g_{\mathrm{0}}+\left(\mathrm{1}+{\displaystyle \frac{g_{\mathrm{1}}}{\sqrt{D}}}\right)\cdot {\displaystyle \frac{A_{\mathrm{net}}\left(\mathit{\beta }\right)}{C_{\mathrm{a}}}}\end{array}$$

where D is leaf-to-air vapour pressure deficit. g₀ and g₁ are empirical parameters.

Figure 1 Schematic of the modeling framework. Abbreviations: gs, stomatal conductance; ET, evapotranspiration; GPP, gross primary productivity; Ψ63: the water potential when 63 % of hydraulic conductance is lost, b: shape parameter in Weibull function; K_max, maximum hydraulic conductance. (a) Workflow of simulations and data constraints. (b) Conceptual updating of leaf water potential (Ψ); in empirical stomatal models, β is applied either to g1 or to V_cmax to represent water-stress effects, and β also scales hydraulic conductance ($K=K_{\max }\times \mathit{\beta }$), thereby affecting Ψ. (c) Time series of percentage loss of hydraulic conductance under the three recovery scenarios (full, partial, no recovery).

For the optimality-based model, stomatal conductance maximizes the difference between leaf-level carbon gain and hydraulic risk, with risk represented through leaf hydraulics and their coupling to photosynthesis (Eq. 6).

$$\begin{array}{}\text{(6)}& g_{\mathrm{s}}=\mathrm{argmax}\left(A_{\mathrm{net}}-A_{\mathrm{net}}{\displaystyle \frac{E}{E_{\mathrm{crit}}}}\right)\end{array}$$

where E is the leaf-level transpiration rate, and E_crit is the critical transpiration rate at which the leaf xylem hydraulic conductance decreases to 0.1 % of the maximum value (namely at 99.9 % loss of hydraulic conductance). A_net is the net photosynthesis rate. β does not directly scale stomatal conductance; instead, it is applied to hydraulic conductance, thereby influencing transpiration-driven water flow.

2.5 Model parameters

We utilized vegetation traits such as the vulnerability curve (VC), V_cmax at 25 °C, semi-empirical slope parameter (g1) that characterizes the efficiency of carbon gain relative to water cost, and maximum hydraulic conductance (K_max), from published references wherever possible. The well-watered V_cmax value was obtained from $A/C_{\mathrm{i}}$ curves measurements made at the site in May 2012, when the plants had not been hit by the severe water stress and Ψ_leaf,pd was close to 0 MPa (Gu et al., 2010). We took the average of the extracted V_cmax from these measurements (65.8±18.8 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). We used the “fitaci” function in package “plantecophys” version 1.4-6 (Duursma, 2015) to derive the fitted V_cmax at 25 °C. The root mean squared error from fitting is between 0.48–2.68 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. We used g1 value for deciduous broadleaf forests from De Kauwe et al. (2015).

For the VC of white oak, we used the sigmoid format parameters from Benson et al. (2022) and transformed it to Weibull function format, which better fits the observations. VC was used to compute the stress factor β (Eq. 7). Ψ63 indicates the xylem pressure at ∼63 % loss of conductivity, and b is the shape parameter (Neufeld et al., 1992). This β formula is the same for stomatal and non-stomatal limitations.

$$\begin{array}{}\text{(7)}& \mathit{\beta }=e^{-(\frac{\mathit{\psi }}{\mathit{\psi }\mathrm{63}}{)}^b}\end{array}$$

where $\mathrm{\Psi }\mathrm{63}=-\mathrm{2.94}$ MPa, b=2.52.

We did not find measurements for K_max, therefore, this parameter was fitted by minimizing the error of ET between the model simulations and eddy-covariance observations over the pre-recovery periods (Eq. 8). Different K_max values were obtained depending on the chosen stomatal models, the location of β, and recovery assumptions (Table S1 in the Supplement). We also optimized K_max over longer pre-drought periods, extending to earlier years, which revealed notable inter-annual variability. For instance, in years with less water stress, such as 2010, the optimized K_max was significantly higher (2.81 mol m⁻² s⁻¹ MPa⁻¹) compared to 2012, a drought year (0.89 mol m⁻² s⁻¹ MPa⁻¹). Consequently, we optimized K_max using 2012 data alone.

$$\begin{array}{}\text{(8)}& \text{For }{K}_{\max }:\text{ error function }=\sum {\displaystyle \frac{\mathrm{|}{\text{ET}}_{\mathrm{mod}}-{\text{ET}}_{\mathrm{obs}}\mathrm{|}}{\text{SD}\left({\text{ET}}_{\mathrm{obs}}\right)}}\end{array}$$

Subscripts of “obs” and “mod” denote the flux observation and modeled term. The operator ∑ takes the sum of the absolute difference between flux observation and model output. The operator “SD” takes the standard deviation of flux observation.

2.6 Hydraulic conductance recovery scenarios 

Three scenarios are considered regarding the recovery of hydraulic conductance for each stomatal model setup, including optimality, Medlyn-V_cmax, and Medlyn-g1. In total, the “partial recovery scenario” falls between the “full recovery scenario” and the “no recovery scenario” (Fig. 1c).

• Full recovery scenario (FR): Hydraulic conductance can either increase or decrease and the shape of the vulnerability curve remains constant, i.e. is not affected by the previous severe water stress condition. In other words, the water potential and hydraulic conductance can be updated without specific boundaries (Eq. 9). Only transient (or immediate) responses are considered.

  K^FR(t) is updated according to the water stress factor β (Eq. 9).

  $$\begin{array}{}\text{(9)}& K^{\mathrm{FR}}\left(t\right)=K_{\max }\cdot \mathit{\beta }\left(\mathit{\psi }\right(t\left)\right)\end{array}$$

  ψ(t) is leaf water potential at time t. K^FR(t) is hydraulic conductance under FR at time t.

• Partial recovery scenario (PR): In reality, the full recovery of the hydraulic system may not necessarily occur due to irreversible processes such as a broken hydraulic transport system (e.g., embolism). Recovery can vary depending on the time scale; for instance, the construction of new xylem over time may facilitate greater recovery, while immediate recovery following a drought may be limited. The extent of potential plant recovery can be inferred from the response of carbon and water fluxes to post-drought climate conditions by applying a similar concept of optimizing ecosystem models through assimilating eddy-covariance measurements in Wolf et al. (2006). By minimizing the difference between modeled and observed carbon and water fluxes, we diagnostically derived the water stress time series and updated modeled leaf water potential at each time step. β is estimated by minimizing the model-observation mismatch (Eq. 10). K^PR(t) is then updated using the inverted β (Eq. 11).

  $$\begin{array}{}\text{(10)}& {\displaystyle }\begin{array}{rl}{\mathit{\beta }}^{\ast }\left(t\right)& =\arg \sum _t\left\{{\left[{\displaystyle \frac{{\text{GPP}}_{\mathrm{mod}}(t;\mathit{\beta })-{\text{GPP}}_{\mathrm{obs}}\left(t\right)}{\text{SD}\left({\text{GPP}}_{\mathrm{obs}}\right)}}\right]}^{\mathrm{2}}\right.\\ & \left.+{\left[{\displaystyle \frac{{\text{ET}}_{\mathrm{mod}}(t;\mathit{\beta })-{\text{ET}}_{\mathrm{obs}}\left(t\right)}{\text{SD}\left({\text{ET}}_{\mathrm{obs}}\right)}}\right]}^{\mathrm{2}}\right\}\end{array}\text{(11)}& {\displaystyle }{K}^{\mathrm{PR}}\left(t\right)=K_{\max }\cdot {\mathit{\beta }}^{\ast }\left(t\right)\end{array}$$

  K^PR(t) is hydraulic conductance under PR at time t.

  In empirical stomatal models, β serves as an integrated proxy for the recovery state of both hydraulic and biochemical or stomatal processes, and their dynamics jointly contribute to the simulated drought legacy effects. Comparing Medlyn-V_cmax and Medlyn-g1 configurations reveals which process, biochemical or nonbiochemical, better explains the recovery of carbon and water fluxes.

• No recovery scenario (NR): Hydraulic conductance can't increase (recover) over time. During severe water stress events, hydraulic conductance decreases due to a drop in water potential. Subsequently, during the post-drought period, hydraulic conductance cannot increase or return to its pre-drought state. K^NR(t) is constrained to be non-increasing (Eqs. 12–13). When drought causes ψ(t) to become more negative, instantaneous hydraulic conductance K^inst drops and K^NR drops. After drought, ψ(t) becomes less negative and K^inst increases. However, the minimum operator in Eq. (13) prevents recovery of K^NR.

  $$\begin{array}{}\text{(12)}& {\displaystyle }{K}^{\mathrm{inst}}\left(t\right)=K_{\max }\cdot \mathit{\beta }\left(\mathit{\psi }\right(t\left)\right)\text{(13)}& {\displaystyle }{K}^{\mathrm{NR}}\left(t\right)=min\left(K^{\mathrm{NR}}\right(t-\mathrm{1}),K^{\mathrm{inst}}(t\left)\right)\end{array}$$

  where K^NR(t) is hydraulic conductance under NR at time t.

All three recovery scenarios share the same initial condition corresponding to the pre-drought state (Eq. 14).

$$\begin{array}{}\text{(14)}& K^{\mathrm{FR}}\left(t_{\mathrm{0}}\right)=K_{\max }\cdot \mathit{\beta }\left(\mathit{\psi }\right(t_{\mathrm{0}}\left)\right)\end{array}$$

where t₀ represents time before recovery.

Figure 2 Drought response of meteorology, ecosystem fluxes, and plant water status. Time series of monthly mean air temperature (T_air), monthly maximum air temperature (T_max), monthly mean precipitation rate, vapor pressure deficit (VPD), near-surface soil moisture, GPP, NEE, ET, and community predawn leaf water potential (Ψ_leaf,pd). The 2012 event is shown as a bold black line; the 2004–2019 record (excluding 2012) appears as gray shading indicating mean ±1σ.

3 Results

3.1 2012 drought conditions

In 2012, the Central US experienced a severe drought that impacted various climatological variables, carbon and water fluxes as shown in Fig. 2 for the MOFLUX site. During the earlier growing season, the air temperature exceeded records from other years, and precipitation reached the upper limit of recorded range. These conditions resulted in an overshoot (or peak flux) detected in the fluxes of GPP and ET. However, in the summer months of June and July, the monthly mean and maximum air temperature were 3 °C higher than those of other years spanning from 2007 to 2019, while summer precipitation only reached 30 % of normal levels. Furthermore, the vapour pressure deficit (VPD) in July was almost double that of other years. Periodic measurements of community Ψ_leaf,pd also showed extremely negative values as low as −4 MPa near the end of August, which was a rare occurrence. This meteorological drought led to a significant reduction in LAI, GPP, NEE and ET in the middle of 2012 (Figs. 2, S1). However, in September, soil water levels began to replenish, and carbon and water fluxes gradually started to recover. While climatological stress eased at that time, fluxes only reached the lower end of the expected range from multi-year records.

Figure 3 Model performance for GPP and ET across four periods in 2012. Curves show the optimality-based stomatal model (cyan) and the Medlyn model with the water-stress factor applied to Vcmax (red) or to g1 (blue). FR, full recovery; NR, no recovery. Error bars indicate variability within each period.

3.2 Model performance

All tested stomatal schemes from CliMA Land demonstrated their ability to capture seasonal variations in carbon and water fluxes during the growing season of 2012 under the assumption of full recovery (Figs. 3, S2). Prior to the recovery period, model simulations performed similarly (Fig. 3). The performance of the optimality model and empirical stomatal models did not show distinct differences when incorporating a stress factor calculated from the loss of hydraulic conductance. Using the Medlyn stomatal model, we observed no significant differences in modeled ET between simulations where the stress factor was applied to V_cmax and g1 over the entire growing season (Wilcoxon rank test: P=0.13). However, applying the stress factor to g1 produced a higher GPP simulation compared to applying the stress factor to V_cmax (Wilcoxon rank test: P<0.001). Additionally, simulations using the stress factor on V_cmax exhibited lower RMSE for GPP at the daily scale compared to those applying the stress factor to g1 (1.40 vs. 2.67 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

Analyzing different time periods, GPP simulations showed higher RMSE during the recovery period compared to the peak, flash stress, and wilted periods, when the stress factor was applied to V_cmax (Fig. S3). During the recovery period, applying the stress factor to g1 resulted in a larger positive bias in GPP simulations relative to observations. For ET, simulations showed comparable lower biases when applying the stress factor to either V_cmax or g1 across all time periods, with no significant differences between ET simulations and flux observations over the entire growing season (Wilcoxon rank test: P>0.10).

Overall, under the full recovery assumption, model simulations tended to overestimate carbon fluxes after the drought event, especially when the stress factor was applied to g1, while accurately capturing ET dynamics (Fig. 3). In contrast, under the no-recovery assumption, modeled GPP aligned more closely with observations during the recovery period when using the optimality-based stomatal model and the Medlyn-g1 setup, though ET simulations showed an underestimation (Fig. 3). These comparisons between full and no-recovery scenarios highlight a lagged recovery process, underscoring the need to incorporate delayed recovery mechanisms in modeling carbon and water fluxes. Using the hydraulic trait parameters of Acer produced similar modeling outcomes (Fig. S4), although the underestimation of photosynthesis was relatively smaller.

Figure 4 Carbon–water coupling strength across drought stages. Pearson correlation between daily GPP⋅VPD^0.5 and ET across four periods (peak, flash-stress, wilted, recovery). Points show flux-tower observations (black) and model simulations: optimality-based stomatal model (cyan) and Medlyn with the water-stress factor applied to V_cmax (red) or g1 (blue). Correlations are computed from daily values and averaged within each period. Error bars on observations denote interannual variability; the multi-year mean is calculated over 2004–2019 excluding 2012.

3.3 Delayed recovery of carbon fluxes

The delayed recovery of carbon fluxes relative to water fluxes suggests a lagged response driven by either stomatal or biochemical limitations. To identify the presence of non-stomatal limitations in the eddy-covariance data, we computed the daily correlation between GPP and ET (Fig. 4). Our findings indicate that this correlation coefficient, representing the strength of carbon-water coupling, was relatively high during the peak flux and flash stress periods, suggesting a tightly coupled carbon and water cycle. However, during the wilted period, this correlation dropped significantly to 0.54 (compared to a multi-year average of 0.89). During the recovery period, the correlation rose to around 0.74 – higher than in the wilted period but still below the multi-year average. The weaker GPP-ET coupling derived from flux observations indicates the emergence of non-stomatal limitations during drought, which can persist into the recovery phase. In contrast, model simulations tend to overestimate this coupling, likely because GPP and transpiration are tightly linked through stomatal conductance.

Given the persistence of non-stomatal limitations after drought, we investigated whether our model could capture this lagged recovery behavior by imposing additional non-stomatal limitations. It should be noted that the water stress is derived from loss of hydraulic conductance in our modeling framework, and further applied to either V_cmax or g1, such legacy effect includes also the delayed recovery of hydraulic conductance.

Figure 5 Improvement in model performance after accounting for drought legacy effects (recovery period only). Scatterplots compare RMSE under the partial-recovery scenario (RMSE_PR) against the full-recovery baseline (RMSE_FR) for GPP and ET. Symbols denote simulations assuming full recovery (black circles) and partial recovery using the optimality-based model (cyan) or the Medlyn model with the water-stress factor applied to V_cmax (red) or g1 (blue). RMSEs are computed against flux observations; lower values indicate better fit.

We found that introducing a partial recovery setup with a fixed bound is inadequate to capture the post-drought responses of carbon and water fluxes, necessitating the data assimilation of using observations as constraint (Figs. S5, S6). After assimilating flux observations as constraints, we observed a slower alleviation of water stress compared to the full recovery assumption and a similarly delayed increase in carbon fluxes (Figs. S7, S8). When comparing improvements in RMSE for GPP and ET between the model simulation and flux observation (shown as “PR” relative to “FR”), we found that incorporating lagged recovery of V_cmax enhanced model performance for both carbon and water fluxes (Fig. 5). In contrast, applying the stress factor to g1 reduced RMSE for carbon fluxes but introduced a greater ET bias, revealing that this model setup fails to yield consistent improvements for both fluxes. These contrasting performances can be explained by the representations in the two model setups. Applying β to V_cmax represents a biochemical limitation, directly reducing photosynthetic capacity and indirectly constraining stomatal conductance via reduced assimilation. This leads to a stronger and more persistent reduction in GPP, with comparatively weaker impacts on transpiration. In contrast, applying β to g1 represents a non-biochemical (stomatal) limitation, directly reducing stomatal conductance and CO? diffusion, thereby suppressing transpiration and secondarily reducing photosynthesis. As a result, when model performance of transpiration is similar between the two setups, as in the FR scenario, where transpiration is used to constrain K_max, photosynthesis simulated by the Medlyn-g1 tends to be slightly higher than that from the Medlyn-V_cmax configuration (Fig. 5b, c). Similarly, under the PR scenario, reducing the overestimation of carbon fluxes is accompanied by an underestimation of transpiration (Fig. 5c, f). In the optimality-based stomatal model, β acts on hydraulic conductance only, and thus affects stomatal conductance indirectly, resulting in a response similar to the Medlyn-g1 formulation and a comparable underestimation of transpiration. Overall, the flux data and our modeling experiments suggest that an impairment in V_cmax delays the recovery of carbon fluxes, with the recovery process spanning days to weeks rather than occurring immediately. This finding highlights the crucial role of biochemical limitations in shaping the post-drought recovery trajectory of terrestrial ecosystems and the asymmetric recoveries of GPP and ET. Instead of relying solely on a fixed $T/\text{ET}$ ratio, we applied the ET partitioning approach of Zhou et al. (2016) to directly compare modeled transpiration with observation-based partitioned transpiration. This method allows the $T/\text{ET}$ ratio to vary dynamically. We find that using either a fixed $T/\text{ET}$ ratio or observation-based partitioned transpiration does not affect our conclusions, namely, that biochemical limitations delay the recovery of both carbon and water fluxes (Figs. S9, S10).

Figure 6 Diurnal cycles of GPP and ET across drought stages. Panels show (a–d) GPP and (e–h) ET for Peak flux, Flash stress, Wilted, and Recovery. Curves compare the optimality-based stomatal model (cyan) and the Medlyn model with the water-stress factor applied to V_cmax (red) or g1 (blue). Scenario keys: FR, full recovery; PR, partial recovery; NR, no recovery.

When analyzing the diurnal cycle, the model simulations successfully captured the overall shape of the diurnal changes in carbon and water fluxes (Fig. 6). During the peak flux period, both carbon and water fluxes reached their maximum values at noon. As conditions transitioned into the flash stress period, the peak flux shifted earlier to around 10:00 am (LT). With increasing water stress, the peak flux occurred even earlier in the morning. During the recovery periods, the diurnal shape became flatter, with the peak value returning to around noon and remaining relatively stable between 10:00 am and 03:00 pm. Under the full recovery assumption, the model predicted GPP at noon to be twice as high as observed when the stress factor was applied directly to g1, while applying the stress factor to V_cmax led to a 35 % overestimation relative to observations during the 10:00 am to 02:00 pm window. However, the partial recovery approach helped reduce the model-observation bias, particularly during daytime hours. Model simulations using both the optimality-based stomatal model and empirical stomatal models with the stress factor applied to V_cmax exhibited less bias across the four time periods.

It is important to note, however, that modeled peak ET around noon during the recovery periods was slightly dampened. From a modeling perspective, a decrease in V_cmax reduces carbon assimilation, which in turn suppresses stomatal conductance and transpiration. Thus, the stress applied to V_cmax effectively encompasses lumped effects from changes in Rubisco enzyme content (or activation) and environmental stress. The slightly weaker ET performance could also be due to potential shifts in the $T/\text{ET}$ ratio between periods, highlighting the need for further exploration of inconsistent drought recovery in carbon and water fluxes. Since $T/\text{ET}$ typically increases under dry conditions (Nie et al., 2021), we tested the model by increasing $T/\text{ET}$ during the wilted and recovery periods and recalibrated it accordingly. This recalibration required slightly higher K_max values to align with the elevated $T/\text{ET}$ during the drought. Nevertheless, this adjustment does not alter our conclusion: the delayed recovery of carbon fluxes can be effectively addressed by incorporating biochemical limitations, thereby reducing the model-observation gap.

4 Discussion

This study examines the recovery of ecosystem carbon and water fluxes following the extreme drought in the Central US in 2012, using in-situ eddy-covariance flux measurements and the land surface model CliMA. Our findings show that CliMA accurately captured ecosystem behavior during the peak flux, flash stress, and wilted periods (Fig. 3). However, during the recovery period, the model tended to overestimate carbon fluxes when assuming a complete recovery of hydraulic conductance immediately after water stress subsided. This indicates that despite the resumption of water fluxes, lingering stomatal or non-stomatal limitations may impede the recovery of carbon fluxes. Our findings align with the analysis by He et al. (2018), which found that GPP takes longer to recover than ET. The quick recovery behavior of ET should be treated with caution as it contains both plant transpiration and soil evaporation, with the latter recovering rapidly after rewetting and potentially masking the slower recovery of transpiration. Supporting this interpretation, observation-based partitioned transpiration shows an intermediate recovery between the simulated full- and no-recovery scenarios, consistent with the presence of drought legacy effects in plant water use. In contrast, this pattern is less apparent in ET, suggesting that enhanced soil evaporation following rainfall may compensate for reduced transpiration and obscure the persistence of vegetation drought legacy. However, unlike He et al. (2018), who attributed the recovery of carbon and water fluxes primarily to canopy conductance, our findings indicate that introducing a progressive biochemical recovery mechanism into empirical stomatal conductance models improves model performance for both carbon and water fluxes. Specifically, we conceptualized the post-drought flux response as a composite of transient and lagged responses to environmental conditions. In this framework, the “full recovery” scenario considers only transient responses, while the “partial recovery” scenario accounts for lagged responses as well. This approach involves constructing a transient baseline for the post-drought period and quantifying deviations from the ecosystem's actual response, thereby capturing the legacy effects. Our methodology is similar to the legacy detection framework used by Yu et al. (2022), which uses a machine learning model to separate transient and lagged responses, with the residuals representing the lagged component. Based on our tests of legacy detection methods, this study underscores the significance of V_cmax impairment in shaping delayed recovery of carbon fluxes, with hydraulic stress being jointly considered in drought legacy effect. The importance of V_cmax downregulation has been noted in previous studies, such as Yang et al. (2019), which demonstrated that imposing only hydraulic limitations can capture drought responses but often in a physically unrealistic manner. Other studies have also proposed adjustments to V_cmax that incorporate climate effects, yielding a better model-observation agreement (Galmés et al., 2013; Gimeno et al., 2019). Therefore, accurately modeling vegetation responses to drought requires explicitly incorporating and quantifying both stomatal and non-stomatal constraints.

While our findings highlight the significant role of lagged responses of V_cmax at the MOFLUX site, this does not imply that this phenomenon fully explains all post-drought behaviors. Model simulations at another site, US-MMS, similarly showed marked declines in both carbon and water fluxes during the drought, along with overestimations of GPP and ET in the post-drought period (Fig. S11). Under the partial recovery scenario, applying the β factor to either V_cmax or g1 improved the model performance in reducing these overestimations (Fig. S12). Additionally, simulations at four other sites affected by the 2012 drought successfully reproduced the seasonal dynamics of carbon and water fluxes without notable overestimations (Figs. S13, S14). These results suggest that challenges remain in parameterized stomatal and non-stomatal limitations, as the dominant limitation may vary depending on environmental conditions. The choice of whether to prioritize V_cmax or g1 should depend on the severity and duration of drought stress, underscoring the need for flexible, context-specific modeling approaches. Short-term mild drought events may not significantly impair Rubisco activity, with stomatal closure serving as the primary limitation and non-stomatal effects playing a minor role (Wilson et al., 2000). Under prolonged drought conditions, both stomatal and non-stomatal limitations likely work together to regulate physiological responses (Galmés et al., 2007). For example, observations at the MOFLUX site revealed carbon-water decoupling during drought and recovery periods, implying the important role of biochemical limitations in regulating these processes. Stomatal limitations are observed across a wide range of soil moisture levels, while non-stomatal limitations become more prominent under moderate to extreme water deficits (Drake et al., 2017). Reductions in the stomatal slope parameter often precede reductions in V_cmax, although Rubisco content may increase under prolonged drought due to acclimation (Zhou et al., 2018). Similarly, Flexas et al. (2006) noted a strong connection between stomatal conductance and V_cmax, with Rubisco activity declining only when stomatal conductance drops below a critical threshold. The transition between these two types of limitations could depend on factors such as leaf position and water content (Song et al., 2020), while Drake et al. (2017) argued that these limitations often act concurrently rather than in distinct phases. As decreased Rubisco activity can hinder rapid recovery even after rehydration (Varone et al., 2012), understanding the thresholds beyond which non-stomatal limitations dominate requires further sensitivity experiments to mechanistically unravel the interconnected dynamics of carbon and water fluxes during recovery from extreme climate events.

Although the model simulation in this study parameterized a single set of plant hydraulics parameters based on the site's dominant tree species of oak, known for its drought-resistant characteristics, we recognize that drought responses can vary significantly among tree species with diverse hydraulic regulation strategies (Drake et al., 2017; Zhou et al., 2014). For instance, Rehschuh and Ruehr (2022) observed no necessary connection between hydraulic impairment and lagged photosynthetic recovery in Scots pine. Xeric tree species, on the other hand, tend to experience more negative leaf water potentials at which photosynthetic efficiency declines sharply (Zhou et al., 2014). Species-specific differences also emerge in the relative contributions of stomatal and biochemical limitations (Zhou et al., 2015). These findings underscore the need to extend our analysis framework to encompass a broader range of tree species, considering variations in life stage, the drought timing, intensity, and duration of stress period, as well as species-specific resistance and resilience to drought.

Our findings emphasize the significance of accounting for drought-induced damage to V_cmax and its subsequent delayed impacts on carbon and water fluxes. This study represents a preliminary application of our modeling framework at a representative eddy-covariance site impacted by the 2012 mega drought, aimed at identifying lagged responses in carbon and water flux as well as stomatal and non-stomatal limitations. Our analysis revealed that, beyond the transient responses of carbon flux, a lagged biochemical recovery helps bridge the model-observation gap. This site-level research lays the groundwork for potential regional expansions to delve further into the factors influencing drought recovery and to validate the generalizability of the identified responses. The present framework should therefore be viewed as a diagnostic tool for identifying drought recovery processes at observation-rich sites. Extending it to regional and global applications will require predictive relationships between inferred recovery characteristics and observable ecosystem properties, rather than prescribing a universal recovery limit. For example, key recovery characteristics, such as the time required for the release of accumulated water deficit and the upper bound of hydraulic conductance recovery, can be inferred from flux tower data. These features can subsequently be linked to commonly available environmental variables using nonlinear regression or machine-learning approaches, enabling their prediction across broader spatial domains without eddy-covariance coverage. To enhance predictions of ecosystem resilience to future climate challenges, further studies should evaluate the asynchronous responses of biochemical and stomatal limitations to drought and explore their underlying drivers. A critical implication of our research is that most current Earth system models lack the inclusion of this vital mechanism. As a result, their projections of carbon sequestration may be overly optimistic, raising concerns about the reliability of these models in informing mitigation policies. By incorporating these crucial processes into Earth system models, we can develop more reliable projections and help make informed decisions to effectively address the climate change challenges. Additionally, considering the potential asynchronous recovery rates of carbon and water fluxes will provide valuable insights into the possible carbon-water decoupling behavior during extreme events (De Kauwe et al., 2019; Krich et al., 2022).

To reduce the long-term impacts of drought legacy effects on forest ecosystems, management strategies should focus on enhancing resilience and adaptive capacity. Selective thinning can lessen competition for water, enabling surviving trees to recover more effectively after drought. Increasing species diversity – particularly through the inclusion of drought-tolerant and deep-rooted species – can help buffer stands against future water stress. Regular monitoring of physiological indicators, such as canopy water content and stomatal conductance, can reveal delayed drought responses and guide timely interventions. Integrating drought legacy effects into forest growth models will be critical for sustaining carbon sequestration and maintaining ecosystem services under a changing climate.

5 Conclusion

In this study, we conducted a comprehensive investigation of the drought response during the 2012 drought in the Central US, utilizing in-situ eddy-covariance flux measurements and the CliMA model. Our drought response simulations showed that assuming full recovery of hydraulic conductance leads to an overestimation of carbon flux recovery. To address this limitation, we introduced a progressive biochemical recovery process into the model diagnostically, accounting for both delayed hydraulic recovery and impairment of Rubisco activity. This modification led to enhanced model performance in simulating both carbon and water fluxes, offering more accurate representation of the ecosystem's recovery dynamics. The implications of this finding are particularly relevant in the context of increasing frequency and intensity of extreme climate events. Ignoring legacy effects and relying on overly optimistic estimations of future land carbon uptake could be misleading, emphasizing the need for more in-depth studies on understanding ecosystem resilience to extreme climate events.