Forest productivity recovery or collapse? Model‐data integration insights on drought‐induced tipping points

More frequent and severe droughts are driving increased forest mortality around the globe. We urgently need to describe and predict how drought affects forest carbon cycling and identify thresholds of environmental stress that trigger ecosystem collapse. Quantifying the effects of drought at an ecosystem level is complex because dynamic climate–plant relationships can cause rapid and/or prolonged shifts in carbon balance. We employ the CARbon DAta MOdel fraMework (CARDAMOM) to investigate legacy effects of drought on forest carbon pools and fluxes. Our Bayesian model‐data fusion approach uses tower observed meteorological forcing and carbon fluxes to determine the response and sensitivity of aboveground and belowground ecological processes associated with the 2012–2015 California drought. Our study area is a mid‐montane mixed conifer forest in the Southern Sierras. CARDAMOM constrained with gross primary productivity (GPP) estimates covering 2011–2017 show a ~75% reduction in GPP, compared to negligible GPP change when constrained with 2011 only. Precipitation across 2012–2015 was 45% (474 mm) lower than the historical average and drove a cascading depletion in soil moisture and carbon pools (foliar, labile, roots, and litter). Adding 157 mm during an especially stressful year (2014, annual rainfall = 293 mm) led to a smaller depletion of water and carbon pools, steering the ecosystem away from a state of GPP tipping‐point collapse to recovery. We present novel process‐driven insights that demonstrate the sensitivity of GPP collapse to ecosystem foliar carbon and soil moisture states—showing that the full extent of GPP response takes several years to arise. Thus, long‐term changes in soil moisture and carbon pools can provide a mechanistic link between drought and forest mortality. Our study provides an example for how key precipitation threshold ranges can influence forest productivity, making them useful for monitoring and predicting forest mortality events.

addition, the consequences of drought are further compounded as risks of forest fires and insect outbreaks increase (McDowell et al., 2020). These processes can limit both the forest's resilience to stress and its recovery (Anderegg, Flint, et al., 2015;Anderegg, Schwalm, et al., 2015).
Here we examine ecosystem resilience as the capacity to assimilate C at pre-drought conditions following sustained environmental stress-measured as gross primary productivity (GPP). Extended periods of drought stress and extreme temperatures from climate change are likely to alter photosynthetic capacity and overall ecosystem function (Brodribb et al., 2020;Hammond et al., 2022). For example, drought stress from the severe 2003 European drought was linked to 30% reduction in GPP across Europe (Ciais et al., 2005).
Based on this example, a more resilient forest can either tolerate or recover from a water stress event, eventually recovering GPP to predrought conditions. If instead the forest is stressed beyond recovery, post-drought GPP may be insufficient for C accumulation and allocation to plant tissues-in particular, foliar biomass-which, in turn, will further reduce GPP. Leaf area index (LAI) is a primary determinant of photosynthesis and C allocation (Running & Coughlan, 1988) and can explain approximately 33%-80% of the variation in GPP in different ecosystems (Flack-Prain et al., 2019;Street et al., 2007). The feedback between declining GPP and declining foliar biomass can in principle lead to a partial or complete GPP collapse, indicating a shift in the equilibrium state (foliar biomass) and function (GPP) of the ecosystem ( Figure 1).
The equilibrium state of an ecosystem depends on both external forces and its internal processes (Luo & Weng, 2011) which impact its ability to photosynthesize. The dynamics of the ecosystem can result in systems where a gradual change in the environment can cause a drastic, non-recoverable change in the equilibrium state of F I G U R E 1 Potential consequences of drought on the state of a forest ecosystem. (a) GPP decreases in a healthy forest ecosystem (black line) during (b) sustained drought. Depending on the resilience of the ecosystem, ecosystem GPP recovers under pre-drought precipitation levels (green line) or ecosystem GPP falls below the point of recovery and the ecosystem collapses despite precipitation conditions improving (red line). (c) The ecosystem state has a direct effect on GPP, and in this context, is primarily controlled by precipitation deficit. There is a stable equilibrium state (solid line) between ecosystem state and precipitation deficit until the deficit reaches a tipping point (Tp) at which the upper stable equilibrium disappears. The dashed line indicates an unstable equilibrium state and collapse of the ecosystem. The arrows represent the trajectory of a return to steady state. GPP, gross primary productivity. the ecosystem (Figure 1c). The critical threshold at which an ecosystem can abruptly shift to an alternate state is known as a "tipping point" (Angeli et al., 2004). Tipping points can occur across a range of spatial and temporal scales where disturbances like climate change and deforestation can cause dramatic shifts in vegetation type and/ or ecosystem function (Reyer et al., 2015). Despite many theoretical models that describe critical thresholds (Luo & Weng, 2011;Scheffer et al., 2001), and numerous observations of ecosystem change and collapse (Allen et al., 2010;Anderegg, Flint, et al., 2015;Anderegg, Schwalm, et al., 2015;Van Mantgem et al., 2009), quantifying where tipping points lie through well-defined metrics remains extremely difficult (Allen et al., 2015;Forzieri et al., 2022;Trugman et al., 2021).
Californian forests have suffered many severe and prolonged droughts, making it an ideal study area for assessing drought-induced tipping points. A remarkable drought in 2012-2015 is considered the most severe in the last 1200 years (Griffin & Anchukaitis, 2014) and led to the mortality of an estimated 30.8 million trees (Goulden & Bales, 2019;USFS, 2016). Here, we use an eddy covariance site in the southern Sierra Nevada mountains-located in one of the most severely affected forests (Goulden et al., 2012). This mid-montane, ponderosa pine forest experienced 80% canopy mortality from 2011 to 2015 (Goulden & Bales, 2019). Severe soil moisture overdraft from 5 to 15 m soil depth drove tree death at this study site and across the Sierra Nevada (Goulden & Bales, 2019). This dataset provides a unique opportunity for identifying where precipitation induced tipping points for forests could lie and how droughts can irreversibly affect forest C balance.
Multi-year in situ observations provide critical evidence for the processes regulating drought-recovery cycles in disturbed forests.
Ecosystem structural and C cycle observations, including leaf area estimates, eddy covariance flux tower measurements of net ecosystem exchange, and derived estimates of GPP, can all help constrain and resolve the temporal evolution of ecosystem C cycling and their responses to environmental forcings .
Furthermore, reconciling land surface observations with mechanistic C cycle models can provide additional constraints on unknown process parameters, including biomass C states, uncertain allocation and turnover rate parameters, and the overall sensitivity of ecosystem fluxes to climatic variability (Bloom & Williams, 2015;Clark et al., 2011). To this end, terrestrial biosphere models that optimize dynamic constraints on ecosystem parameters against diverse observations tend to show better predictive performance (Famiglietti et al., 2020).
The CARbon DAta MOdel fraMework (CARDAMOM) provides a mechanistic framework for reconciling ecosystem C cycle dynamics with observations (Bloom & Williams, 2015). A key advantage of CARDAMOM compared to other models is that it does not rely on plant functional type information or prior assumptions of plant traits. CARDAMOM has successfully revealed the effect of fire decline on the tropical land C sink , the concurrent effects of disturbance history on tropical carbon balance , the state of global terrestrial C balance (Bloom et al., 2016), and helped diagnose climate variability on evergreen GPP (Stettz et al., 2021). CARDAMOM is a Bayesian model-data fusion approach that uses meteorological forcings and land-surface observations to estimate time-invariant parameters and initial C and water states Bloom & Williams, 2015). This provides a continuous and self-consistent mechanistic representation of climate-forced C and water cycling, permitting a more mechanistic diagnosis of processes such as GPP collapse. The Data Assimilation Linked Ecosystem Carbon model (DALEC; Williams et al., 2005) is a key component of CARDAMOM, representing the key processes regulating ecosystem carbon storage and allocation. These include internal C processes, such as the transfer across pools (e.g., foliar, woody biomass, soil organic carbon), and external forces, such as climate, that together drive the C balance of an ecosystem (Luo & Weng, 2011). Thus, the ecosystem carbon states and their interactions with climate affect an ecosystem's ability to gain carbon through GPP, allowing inference into the potential collapse in forest productivity by tracking these pools over time.
In this study, we use multi-year in situ GPP observations to resolve and assess forest ecosystem drought-induced tipping points; specifically, we address four questions: (1) How well can we capture forest productivity collapse using a terrestrial biosphere model that is constrained with GPP observations? (2) How does drought affect intra-and inter-annual forest carbon dynamics? (3) Did the precipitation deficit exceed an ecosystem critical tipping point and how can we use CARDAMOM to help inform the mechanistic cause of this?
And (4) Would GPP collapse have been avoidable under higher predrought foliar C or wetter soil states? To address these questions, we integrate tower-based net ecosystem C exchanges and meteorological drivers within the CARDAMOM Bayesian model-data fusion framework. Furthermore, we use the emergent ecosystem C and water cycle dynamics to quantitatively resolve the mechanistic links between precipitation deficits on observed GPP declines. In the face of uncertainty on ecosystem dynamics, CARDAMOM provides a unique opportunity to reconcile observed C fluxes and states with underlying model parameters to resolve the processes at play during drought-induced GPP collapse.

| ME THODS
To address the questions outlined above, we (1) analyze the trajectory of ecosystem carbon fluxes and states from in situ and modeled estimates, (2) integrate data into the CARDAMOM model-data fusion framework and (3) assess the sensitivity of GPP to changes in C and water ecosystem states. We describe datasets and CARDAMOM framework in Sections 2.1-2.3 and simulations in Section 2.4.

| Study area and data collection
The study area in the Southern Sierras was severely affected by the 2012-2015 California drought (Goulden & Bales, 2019). We focus on the Southern Sierra Critical Zone Observatory eddy covariance tower denoted "CZO2" (37.03106, −119.2566). The site is at an elevation of 1160 m and dominated by ponderosa pine forest (Pinus ponderosa) with some Californian black oak (Quercus kelloggii) and incense cedar and (Calocedrus decurrens). Soils at the site consist of fine, coarse, and sandy loams derived from granitic parent material, with a weathered bedrock-hard bedrock boundary at 4.5 m (O'Geen et al., 2018). A 40 m eddy covariance flux tower collected data from 2011 to 2017 including CO 2 exchange, evapotranspiration (ET), and meteorology data averaged at half-hourly intervals (Goulden et al., 2012). Half-hourly gross ecosystem exchange was calculated as the difference between observed net ecosystem exchange and respiration, using the nighttime partitioning method over 30-day periods (Goulden et al., 2012) and we used these values to derive a monthly GPP for the site. Historical precipitation data for the area were determined from the ERA5 reanalysis dataset (resolution is 9 km on a reduced Gaussian grid at 1 h; Hersbach et al., 2018). Annual precipitation is calculated using the water year, October to September. From 1982 to 2010, the historical annual average precipitation was 1045 mm year −1 , with a standard deviation of 352 mm year −1 . To quantify precipitation deficits, we calculate the anomaly using z scores by (annual precipitation -historical average)/ standard deviation.

| CARDAMOM framework
CARDAMOM is built around the DALEC model ( Figure S1). DALEC v2 is a carbon cycle model that uses meteorological inputs to drive 33 parameters that assess aboveground and belowground carbon and water ecological processes such as photosynthesis and turnover rates (Table S1, Williams et al., 2005). Some examples of these time-invariant parameters with their respective prior ranges include fraction of NPP to autotrophic respiration (0.2%-0.8%), fraction of litter decomposition per day (0.0001-0.01 per day), and initial plant available water (PAW; 1.0-10,000.0 mm). There is no assumption of plant functional type in CARDAMOM, eliminating the need to rely on prior information for model initialization. Parameter optimization is flexible relative to the meteorological drivers and observational constraints but is restricted to their respective prior ranges. Thus, CARDAMOM is advantageous when modeling multiple ecosystems or in ecosystems with multiple plant species.
DALEC estimates the size and residence times of seven ecosystem pools . There are six C pools (g C m −2 day −1 ): nonstructural carbohydrates (C NSC ), leaves (C F ), litter (C L ), roots (C R ), soil organic matter (C SOM ), woody biomass (C W ), and one PAW pool (mm). Major flux estimates include GPP allocation to the four live pools (g C m −2 day −1 ), C NSC to C F (g m −2 day −1 ) and wood to soil carbon (m −2 day −1 , Table S2). The DALEC model uses PAW to determine soil water deficit stress on photosynthesis. The PAW model is similar to a hydrological "bucket model" and is dynamically derived as the sum of precipitation input, and ET and runoff outputs. The ET semi-empirical parametrization (Boese et al., 2017) is derived using GPP, PAW, and vapor pressure deficit; the runoff parametrization-representing both surface and subsurface lateral water flows .
PAW effectively represents the amount of water available to the rootzone and does not resolve deep or static stores of water unavailable to plants. The PAW state directly influences GPP through a linear scaling function below a "water stress onset" PAW threshold value. We expect water parameters (such as the runoff parametrization constant) to be indirectly informed by GPP, since GPP observations narrow the water balance configurations (i.e., water-related parameter combinations) that best represent the observed GPP responses to soil water deficits. Plant biomass losses are represented by fractional losses from three live biomass C pools (C F , C R , and C W ) to two dead organic C pools (C L and C SOM ). Heterotrophic respiration (Rh) is a function of temperature and precipitation-rather than PAW-given that Rh is expected to be highly correlated with nearsurface moisture, which is more closely related to precipitation.
CARDAMOM is a Bayesian model-data fusion framework that constrains initial conditions and parameters in DALEC using adaptive

Metropolis-Hastings Markov Chain Monte Carlo (AP-MHMCMC).
The starting values for each of the 33 parameters are selected randomly from a uniform prior distribution. To optimize the values for each parameter, we use the AP-MHMCMC algorithm to modify the step size. Likelihood for a suite of parameters for a given iteration is generated and compared using Bayesian inference across 10 8 iterations, over four chains. The CARDAMOM framework checks against 12 ecological dynamical constraints that are simple, defined relationships between certain parameters and help inform the model of more realistic values between the priors (Table S3, Bloom & Williams, 2015). For example, turnover rate of soil organic matter should be slower than that of litter. To reduce parameter uncertainty and to guide more efficient model development, we can include multiple observations such as soil organic carbon, abovegroundbelowground biomass or productivity (Famiglietti et al., 2020).
Posterior distributions provide uncertainty estimates for all parameters, fluxes, and pools. Model output can then be assessed by comparing posterior distributions to priors and against independent observations that were not included during model development.

| Initial-GPP-constrained versus annual-GPP-constrained models for monitoring GPP collapse
We perform two CARDAMOM runs using monthly forcings for DALEC. Inputs from the flux tower are precipitation, temperature, vapor pressure deficit, and solar irradiance (Goulden et al., 2012).
We collect inputs of monthly CO 2 atmospheric concentrations from the Mauna Loa Observatory which vary over time (Keeling & Keeling, 2017). To constrain initial pools sizes, we use a single estimate of aboveground and belowground biomass from a tree inven- CARDAMOM U (initial-GPP constrained) uses monthly GPP EC for 2011 only to constrain initial conditions and then estimates 2012-2017 GPP solely based on observed meteorological forcing and associated climate anomalies. Thus, CARDAMOM U does not leverage empirical knowledge of the drought to make predictions and mostly reflects prior parameter distributions in DALEC. In contrast, CARDAMOM C (annual-GPP-constrained) empirically learns from data as it uses monthly GPP EC for the entire period (2011)(2012)(2013)(2014)(2015)(2016)(2017). This provides parametric constraints on GPP during wet and dry periods and as a result, assimilates GPP by accounting for time variations of GPP EC and meteorological forcing. We assess the performance of CARDAMOM to capture GPP collapse via the correlation coefficient between assimilated GPP from each CARDAMOM run and GPP EC .
We also include independent validation of site-specific soil moisture, eddy covariance ET, and a moderate resolution imaging spectroradiometer (MODIS) estimate of LAI (MOD15A2H, Myneni et al., 2015).
Using four soil moisture probes measuring the top 30 cm depths, we compute interannual change of soil moisture by taking a rolling 12-month mean to reduce the amplitude variation that is associated with drying/wetting of topsoil.

| Resolving the mechanistic links between drought and GPP collapse
We conduct a series of sensitivity scenarios to provide insight on whether GPP collapse is due to (1) precipitation variability and/or (2) initial ecosystem C states and water availability using CARDAMOM C .
To test the sensitivity of GPP to changes in precipitation (experiment 1), we conduct a series of perturbed precipitation runs using CARDAMOM C parameters and predrought (initial) conditions. We modify the monthly precipitation in 2014 (an especially low water year where annual precipitation equals 293 mm) by adding 5%-70% (in 5% increments) of annual precipitation (1045 mm). For example, in the 15% perturbation, we compute the sum of 15% annual precipitation (157 mm) and 2014 precipitation (293 mm) which equals 450 mm. The simulated increase from 293 to 450 mm is equivalent to 1.536 and we use this factor to multiply the precipitation of each month. We plot the change in GPP across 14 simulations ( Figure S4).
For each simulation, we regress GPP against time (all 48 months between 2014 and 2017) and extract the slope. We define "recovery" if the slope is above zero (indicating positive slope) and simulated GPP in 2017 is approximately 7 g C m −2 day −1 . We define the critical threshold range for precipitation tipping points as the minimum annual precipitation between GPP collapse and recovery. This approach allows us to characterize and identify the unstable equilibrium ( Figure 1c).
From this, we assess how C and water pools can drive changes in GPP under different precipitation conditions. After determining the minimum amount of added precipitation required for GPP recovery (from the above sensitivity analysis), we simulate how C and water pools change with additional precipitation. We add precipitation to 2014, the year just prior to GPP collapse. We define the critical threshold range for ecosystem pool tipping points as the range in summer peak values from 2015 between collapse and recovery rainfall simulations.
We then conduct a dedicated sensitivity simulation to investigate if pre-drought (in this case, 2011 conditions) water or C states affect the vulnerability of the forest to GPP collapse (experiment 2). We first augment pre-drought PAW to 2250 mm which is equal to 100% soil moisture in this study and should be seen as a theoretical maximum. A maximum PAW of 2250 mm is reasonable given that other studies have reported 1200 mm PAW in southern Sierra Nevada mid-montane forests (Stocker et al., 2023) and 2017 was a particularly wet year with 1703 mm annual precipitation. We assume that the rooting depth maximum for the site is 4.5 m given that there is weathered bedrock-hard bedrock boundary at 4.5 m and that porosity for the study area is 0.5. Second, we augment predrought C F in five simulations. We set pre-drought C F to 200 g C m −2 which reflects the maximum values observed in this study. We then include four additional simulations by decreasing pre-drought C F (44 g C m −2 day −1 ) from 80% to 20%, in 20% increments. Although the drought breaks from 2016, summer GPP EC maximum further declines to 2.0-2.1 g C m −2 day −1 and winter GPP EC remains at 0 g C m −2 day −1 . From 2011 to 2017, there is an approximate 71%-75% reduction in summer GPP EC and 100% loss of winter GPP EC .

| Precipitation deficit can lead to GPP collapse
These data are also presented in Goulden and Bales (2019).

| Annual-GPP-constrained models can detect drought-induced GPP collapse
Using annual GPP from 2011 to 2017 to optimize CARDAMOM carbon-water cycles and associated processes provide a robust representation of GPP variability in response to climate. CARDAMOM C successfully captures the observed GPP collapse (Figure 2b).  Figure 2b). Furthermore, the CARDAMOM U GPP (GPP U ) seasonal cycle shows peaks that are smaller in amplitude and later in the season. Overall, we find lower agreement of GPP U with GPP EC (Figure 2c, R 2 = .31, RMSE = 2.3 g C m −2 day −1 , and bias = −1.16). To further support the robustness of CARDAMOM C , we find that PAW estimated from CARDAMOM C shows strong correlations with an independent validation of soil moisture (R = .7) compared with CARDAMOM U (R = .6, Figure S2). CARDAMOM C ET has a stronger correlation with eddy covariance ET (R = .61) compared with CARDAMOM U (R = .34, Figure S2). In addition, we find weak correlations with MODIS LAI and CARDAMOM C LAI (R = .29, Figure S2).
There are major differences in parameter estimates between CARDAMOM runs when only the first year of GPP is assimilated (CARDAMOM U ) and when wet and dry periods are

| Reduced PAW can lead to loss of live biomass pools and ecosystem GPP collapse
GPP collapse in 2017 could have been avoided if at least 157 mm, 15% of the annual average, was added to 2014 (Figure 4). Within this simulation, we find the 2017 summer GPP C maximum increases to 7.2 g C m −2 day −1 , closely resembling pre-drought GPP C for the forest (7.4-7.9 g C m −2 day −1 ) and winter GPP C maximum increase to  To highlight the effect of water ecosystem status on C pools, we reduce the precipitation deficit in 2014 from −2.1σ to −1.6σ by adding 15% (157 mm) of the historical average. Growing season peak GPP C increases from ~5 g C day −1 in 2014 to ~7.1 g C day −1 in 2017.
The simulation results show recovery across all six C pools from 2014 ( Figure 5). The small differences between the 2015 summer peaks of actual and simulated pools indicate a potential minimum size of C and water state for the growing season; approximately 825 PAW mm and 18 g C F m −2 .
Pre-drought pool sizes for 2011 water and C ecosystem states also affect GPP resilience. A maximized pre-drought PAW state of 2250 mm supports GPP recovery by 2017 (Figure 6a). When predrought C F is maximized to 200 g C m −2 , we instead see earlier and more severe GPP collapse (Figure 6b). Alternatively, if pre-drought C F reduces to 18 g C F m −2 , 40% of pre-drought foliar C (44 g C F m −2 ), we see GPP recovery from 2017 ( Figure 6b). Overall, we find that the change in pre-drought conditions (in 2011) leads to noticeable changes in GPP in the following years (2012-2017).

| DISCUSS ION
Carbon and water ecosystem states are highly sensitive to precipitation and can impact GPP response of the forest. Severe periods of multi-year droughts can cause a sustained PAW deficit, leading to a decline in new production, and subsequent GPP collapse. A droughtinduced tipping point is reached when a gradual decline in C pools

F I G U R E 4
The change in predicted GPP following a one-time increase in monthly precipitation in the year prior to GPP collapse (2014). The blue box indicates that precipitation was added in 2014 only. Added precipitation ranges from 0% to 70% (of annual precipitation = 1045 mm year −1 ) in 5% intervals (color bar). For comparisons, observed GPP from eddy covariance measurements (grey line) and median assimilated GPP from model. GPP, gross primary productivity. leads to an irreversible shift in the ecosystem; in this case, GPP collapse. Our study highlights the sensitivity of a mid-montane, ponderosa pine forest to accumulated water deficit, with significant lags in GPP recovery attributable to previous years' drought conditions. The tipping point for GPP collapse at the ecosystem scale is complex and challenging to predict. At this pine site with soils composed of fine, coarse, and sandy loams, the critical threshold range for precipitation deficit tipping points could lie from −2.1 to −1.6σ where a deficit greater than −2.1σ (following 2 years of sustained drought) leads to GPP collapse (Figure 7). This ultimately led to significant canopy and tree mortality, as was observed in this region (Goulden & Bales, 2019;USFS, 2016). Significant deviations from historical precipitation averages have driven global tree mortality events. For example, monthly precipitation z scores were significantly lower (−0.21σ ± 0.03) than the long-term average during a tree mortality year (Hammond et al., 2022) and every 35 mm decrease from the historical average leads to an increase in tree mortality risk (1.07 times) and rate (0.85% year −1 , Stovall et al., 2019). The 752 mm deficit in 2014 (which could equate to an increase of 23 times mortality risk and 18.3% mortality rate) leads to a PAW estimate of 446 mm, which is equivalent to 20% soil moisture content. 20% soil moisture levels are typically associated with largely negative soil water potentials in weathered granitic soils (Hubbert et al., 2001), which can induce stomatal closure in drought-stressed plants (Choat et al., 2018).
CARDAMOM C highlights the catastrophic feedback between the C and water cycles that can occur during prolonged drought.
From 2011 to 2015, there is a 75% decline in C F from maximum (186 g C F m −2 in 2013) matching the observed 80% canopy mortality event in the area (Goulden & Bales, 2019) and decline in LAI observed from MODIS data ( Figure S2). Canopy loss is a common plant response to conserve water under drought (Choat et al., 2018) but has detrimental effects on GPP and plant C pools. The reinforcing feedback exacerbates the drought and hinders GPP recovery postdrought. CARDAMOM C suggests that canopy-level C starvation is a result of reduced GPP and could explain the coincident decline in C NSC and C F ; however, low NSCs alone are not necessarily a driver of mortality (Adams et al., 2017). The carbon starvation hypothesis proposes that as trees exhaust their carbon reserves, they become more susceptible to further stress and may ultimately die (McDowell et al., 2008). The decline in available carbohydrates (as a result of decreases in photosynthesis) can weaken the tree's ability to grow new tissue, repair damaged structure, or produce defense compounds, making them more vulnerable to secondary stressors such as insects or disease. Hydraulic failure is another likely cause of forest mortality during drought (Choat et al., 2018;Nolan et al., 2021) and is a possible cause of long-term reduction of forest GPP postdrought (Skelton et al., 2017). While CARDAMOM does not explicitly account for hydraulic failure, attributing reductions in C pools from reduced GPP with reduced PAW suggests hydraulic failure on ecosystem C states. Notably, a model like CARDAMOM is unable to explicitly attribute if carbon starvation or hydraulic failure is or is not at play here, but including information on the coupled carbon-water cycle does offer insight.
If ecosystem resilience relies on its ability to gain C, we can assess vulnerability as a function of C and water resources available increase forest drought resilience as soil dryness is a long-term driving factor of widespread forest mortality (Goulden & Bales, 2019).
In contrast, too dense, or too sparse pre-drought forest canopy can also affect drought resilience. We simulate pre-drought C F perturbations to show the damaging effects of frequent drought cycles on a Californian forest. Reducing pre-drought forest canopy to 18 g C m −2 at this site led to GPP recovery, whereas increasing it to 200 g C m −2 showed the opposite response-GPP collapse was instead earlier and more severe (Figure 6b). Our findings support the negative effects of "structural overshoot" where drought-stressed forests immediately respond to favorable conditions by quickly expanding their canopy but are then at serious risk of losing too much water due to the increased leaf area in the event of subsequent drought (Jump et al., 2017). Indeed, antecedent structural overshoot from 2009 to 2012, along with warmer temperatures F I G U R E 6 The change in predicted GPP following a simulated change in predrought (a) PAW or (b) C F states. Scenarios included a maximum of 2250 PAW (mm), 200 g C F m −2 , and decreases of the predrought C F (44 g C F m −2 ) from 80% to 20%, in 20% increments. For comparisons, observed GPP from eddy covariance measurements (grey line) is shown. GPP, gross primary productivity; PAW, plant available water. and precipitation deficit, drove the severe mortality event in the Sierra mountains (Goulden & Bales, 2019). We note if C F is too low (9.2 C F m −2 ), forest productivity remains below 1 g C m −2 day −1 (Figure 6b), suggesting a sweet-spot between the leaf area required to gain C and to reduce water loss.
It can take years for an ecosystem's response to changing climate to manifest into a measurable change in ecosystem GPP. Trees are long-lived and respond to the variability in current and past climate and soil conditions, leading to potential temporal mismatch between conditions and ecosystem GPP response (i.e. lagged effects). For example, deciduous forests can store enough C to replace the canopy up to four times (Hoch et al., 2003). There is also a time-lag for forest recovery (Figures 4-6). We show that the sensitivity of GPP recovery increases with time, where the full extent of recovery is observed 3 years following added precipitation (Figures 4 and 5). Recovery from water stress can be slow as NSC pools can continue to decline following drought to cope with initial costs of canopy recovery (Lloret et al., 2018) and C is needed to repair hydraulic damage due to water stress before C gain (Ruehr et al., 2019;Trugman et al., 2018). Furthermore, mortality effects can occur years following a drought event in a variety of forest types (Phillips et al., 2010;Trugman et al., 2018). Thus, GPP responses are not necessarily reflective of current precipitation conditions but of the previous years (i.e. legacy effects) and can explain why nonlinear relationships can exist between productivity and precipitation (Felton et al., 2021). Given that lagged effects can contribute 64% of the interannual variability of net biosphere exchange in the tropics , not accounting for this can lead to large uncertainties in terrestrial C balance and productivity (Anderegg, Flint, et al., 2015;Anderegg, Schwalm, et al., 2015).
While the mismatch between drought and ecosystem response has been observed in many studies (Allen et al., 2010;Anderegg, Flint, et al., 2015;Anderegg, Schwalm, et al., 2015;Lloret et al., 2004;Phillips et al., 2009;Trugman et al., 2018), this statistical challenge is difficult to resolve without the guidance of long-term data input from stressful events.
Given that CARDAMOM does not use plant species or ecosystemspecific information, it enables flexibility for ecosystem parameters to better reflect times of stress. GPP observations that cover wet and drought periods provide critical information that constrains the complex feedback between C and water cycles. Greater water sensitivity of photosynthesis is a key reason why CARDAMOM C predicts GPP well across the whole sequence, whereas CARDAMOM U predicts the drought period poorly ( Figure 2) and provides broad estimates of many parameters (Figure 3; Figure S3). For example, CARDAMOM C F I G U R E 7 The relationship between precipitation deficit in 2014 on pre-drought (a) PAW and (b) C F ecosystem state and how it can lead to GPP growth or collapse from 2014 to 2017. With increasing precipitation deficit, all ecosystem states gradually decline to an irreversible state of (A) GPP collapse*, i.e., "tipping point" (Tp), irrespective of ecosystem state value. (a) Availability of water has a positive effect on GPP where (A) actual drought conditions lead to GPP collapse > , (B) reduce precipitation deficit to 595-648 leads to GPP growth, (C, D) increase predrought PAW to 2000 or 2250 leads to GPP growth. (b) Predrought C F has a non-linear effect on GPP where (A, B) actual drought caused severe precipitation deficit that could not support a C F state of 44 or 200, leading to GPP collapse^, (C) reduce precipitation deficit to 595-648 leads to GPP growth, (D) reducing predrought C F to 18 allows GPP growth, (E) very low C F of 09 leads to GPP collapse # . GPP, gross primary productivity; PAW, plant available water.
shows that this pine forest has leaves that can do higher rates of photosynthesis but are more sensitive to water compared with predictions from CARDAMOM U . Estimated LMCA are also better constrained in CARDAMOM C (Figure 3, range 25-75 = 35-72 g C m −2 ) and the median estimate (52 g C m −2 ) is similar to that found in other North American ponderosa pine stands (52-64 g C m −2 , O' Hara & Nagel, 2006). This contrasts with a less constrained and overestimate from CARDAMOM U (90 g C m −2 , range 25-75 = 51-136). Plant trait information from databases like TRY (Kattge et al., 2011) could be incorporated in future work to further assess model parameterizations and provide insight into prior ranges for future models.
CARDAMOM C estimates short leaf longevities, a bimodal budburst day and highly seasonal foliar C pools, which are unusual for an evergreen forest. Our estimate of leaf life span of ~1 year for an evergreen forest would contradict assumptions of this plant functional type but possibly reflects what is happening (lots of needles dropping) in a particularly drought-stressed forest. Indeed, premature leaf abscission is a mechanism that both evergreen and deciduous trees can employ to avoid water stress (Dallstream & Piper, 2021) and reflects the study area's 80% canopy loss from 2012 to 2015 (Goulden & Bales, 2019). The underestimation of GPP in 2015 via CARDAMOM C also suggests that our framework may be overly reliant on leaf fall and flushing to capture seasonal changes in GPP at this site. We find that MODIS LAI shows seasonal changes in LAI that mirror C F from CARDAMOM C but are not as severe ( Figure S2).

Models that allow forests to adjust leaf longevity in response to
severe drought may warrant consideration in ecosystem models used for carbon cycle projections and ecological forecasting and will be important for modeling non-typical ecosystem response during times of stress. The flexibility during parameterization is a key strength of the CARDAMOM framework and future work may need to consider structural changes to the model to better reflect seasonal dynamics of GPP.
CARDAMOM has a static IWUE, which could be an issue when covering wet and drought periods, given that stomatal closure under water stress increases IWUE (Beer et al., 2009). This may explain the limited performance of CARDAMOM to capture ET ( Figure S2).
In this study, we used ET values for validating the model; however, including ET observations during CARDAMOM development could help us better understand if IWUE stays constant or changes after a major mortality event. CARDAMOM C also predicts only modest changes in woody biomass pre-and post-drought. This is because there is no partitioning of living and non-living woody biomass in CARDAMOM-all respiration occurs as a proportion of GPP-and there is no mechanism for woody biomass mortality. This is unlikely to impact our results, however, because woody C does not influence living carbon-water interactions in CARDAMOM. Future models could benefit from better representations of tree mortality; such as whole-plant turnover being a function of C availability, and capping the maximum amount of leaf growth possible per unit area as a function of live woody biomass. However, this would rely on more detailed quantitative data that can track changes in living and nonliving woody biomass.
We provide a framework to assess tipping points leading to ecosystem collapse in other forests, enabling better monitoring of forest vulnerability. While our study is limited to a single ponderosa pine forest in California, applying the CARDAMOM framework to eddy covariance observations of a forest experiencing multi-year drought allows us to calibrate and track ecosystem processes in situ.
We suggest critical threshold ranges for precipitation deficit tipping points, which can be tested against other studies that assess the role of water stress on forests and under future climate scenarios (Forzieri et al., 2022). Quantifying how far an ecosystem is from a precipitation deficit tipping points will also provide information on its resilience to other risk factors such as warming temperatures, forest fires, and insect attacks Stovall et al., 2019;Van Nieuwstadt & Sheil, 2005).

ACK N OWLED G M ENTS
We are grateful to M. Williams, T.L. Smallman, and D.M. Milodowski and to three reviewers whose comments greatly improved the manuscript. We acknowledge D. Malgarini Perez as the graphical artist for Figure 1.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this study are publicly available on the Ameriflux