Ecosystem Water‐Saving Timescale Varies Spatially With Typical Drydown Length

Stomatal optimization theory is a commonly used framework for modeling how plants regulate transpiration in response to the environment. Most stomatal optimization models assume that plants instantaneously optimize a reward function such as carbon gain. However, plants are expected to optimize over longer timescales given the rapid environmental variability they encounter. There are currently no observational constraints on these timescales. Here, a new stomatal model is developed and is used to analyze the timescales over which stomatal closure is optimized. The proposed model assumes plants maximize carbon gain subject to the constraint that they cannot draw down soil moisture below a critical value. The reward is integrated over time, after being weighted by a discount factor that represents the timescale (τ) that a plant considers when optimizing stomatal conductance to save water. The model is simple enough to be analytically solvable, which allows the value of τ to be inferred from observations of stomatal behavior under known environmental conditions. The model is fitted to eddy covariance data in a range of ecosystems, finding the value of τ that best predicts the dynamics of evapotranspiration at each site. Across 82 sites, the climate metrics with the strongest correlation to τ are measures of the average number of dry days between rainfall events. Values of τ are similar in magnitude to the longest such dry period encountered in an average year. The results here shed light on which climate characteristics shape spatial variations in ecosystem‐level water use strategy.


Introduction
Terrestrial vegetation is the main link between the water and carbon cycles, and is an important factor in carbonclimate feedbacks (Friedlingstein et al., 2006;Gentine et al., 2019;Lemordant et al., 2018).Plants adjust their stomatal conductance (g) to control how much water diffuses out and how much carbon diffuses in Buckley (2017) in response to environmental factors such as soil moisture, air humidity, air temperature, and photosynthetically active radiation (Ball et al., 1987).However, current ecosystem models suffer from substantial uncertainty in the response of g to soil moisture, with predictions varying widely between models (Sloan et al., 2021;Trugman et al., 2018).
The precise dependence of a plant's stomatal conductance on environmental factors is often modeled based on optimality theory (Buckley, 2017;Medlyn et al., 2011;Wang et al., 2020), an approach pioneered by Cowan and Farquhar (1977) and Hari et al. (1986), among others.Stomatal optimality models assume that plants have evolved to adjust their stomatal conductance to maximize some underlying reward function.This reward function is typically formulated to capture the underlying tradeoff of stomatal opening-carbon gained is subject to water loss costs or constraints.Numerical or analytical optimization of the reward function is then used to predict actual stomatal conductance in response to a specific set of environmental conditions.Typically, these models assume that the plant maximizes the reward instantaneously and independently at each moment in time (Prentice et al., 2014;Wolf et al., 2016).Many instantaneous stomatal optimality models are now in use (e.g., Eller et al., 2020;Medlyn et al., 2011;Sperry et al., 2017;Y. Wang et al., 2020;Wolf et al., 2016).The differences between these models are mostly due to differences in how they express the cost of water use.Another class of optimality models focuses on non-stomatal controls on photosynthesis, such as carboxylation rates and xylem and phloem transport (R. Dewar et al., 2018;Hölttä et al., 2017;Zhou et al., 2013).
Both stomatal and non-stomatal optimization approaches need to represent the stomatal conductance in a compact manner; as such, they optimize the reward function on an instantaneous timescale instead of over an integration time.Moreover, given the stochastic nature of the environment, plants should maximize the reward not just over current conditions, but also over a future time horizon (Blonder et al., 2023).Stomatal behavior that is instantaneously optimal can be suboptimal when integrated over time (Buckley et al., 2023;Feng et al., 2022).Several previous studies have proposed stomatal optimization models that are time-integrated rather than instantaneous.Some of these models maximize a reward function over a timescale that is assumed to be long and could be as long as the plant's entire lifespan.The mathematical and computational complexity of these models is a key limitation (Buckley, 2023), requiring numerical solutions to solve the problem either deterministically (Cowan, 1982;Manzoni et al., 2011;Mrad et al., 2019;Potkay & Feng, 2023b) or stochastically (Cowan, 1986;Lu et al., 2016Lu et al., , 2020)).To simplify their mathematical complexity, these models typically discretize the timespan into the shortest repeated units of time over which their constraints oscillate.This unit of time could be the length of a drydown period when soil moisture is the constraint (Katul et al., 2010;Manzoni et al., 2011;Mrad et al., 2019), or the length of a year when plant carbon storage is a focus (Potkay & Feng, 2023b).These models then maximize their reward functions over one of these units of time, instead of solving for the stomatal conductance that maximizes the reward over the entire timescale.A fundamental issue with this family of models is that the exact length of the optimization timescale remains challenging to identify or independently confirm (Buckley et al., 2017;R. C. Dewar et al., 2009;Mäkelä et al., 2002).While this timescale is thought to be related to the lengths of drydown periods between rainfall events (Mäkelä et al., 1996), that relation has not been extensively tested.The goal of our study is to estimate the stomatal optimization timescale in a variety of ecosystems without invoking a priori assumptions about limitations on the reward function.
Rather than integrating over a fixed time period, an alternative approach to defining an optimality timescale uses the concept of a discount factor.Such a discount factor has also proven useful in other areas of plant physiology.It has been used to understand variation in leaf size and longevity by Westoby et al. (2000), and is reviewed in the context of plant life-history strategy by Castorena et al. (2022).The use of a discount factor was introduced to stomatal modeling by Mäkelä et al. (1996).Their model assumes that the reward function is integrated over an infinite time period while being multiplied by a discount factor that decreases exponentially over time (e -t/τ ), which keeps the integral finite.Within the discount factor, the Mäkelä model defines the quantity 1/τ as the probability of rainfall resetting soil moisture from a dry state to a wet state.Although the Mäkelä model has been tested for at least one ecosystem (Berninger et al., 1996), the model's mathematical complexity requires a numerical solution or an approximate semi-analytical solution.Thus, it is difficult to parametrize and integrate into stomatal predictions in land surface models to be used globally or across many sites.Here, we hypothesize that the relationship between the discount factor timescale and local rainfall climate can be directly inferred from data, potentially enabling spatially variable parametrization of stomatal optimality models based on rainfall statistics.
In this study, we investigate whether the discount factor timescale depends on the local climate.To address that question, we introduce a time-integrated optimality approach that builds on the Mäkelä et al. (1996) model both theoretically and empirically.An analytical equation is obtained for a dynamic feedback strategy of optimal conductance-an equation that is stated only in terms of present conditions but that takes into account planning for future water use.This new equation clarifies the role of τ, the stomatal optimality timescale for water-saving.The equation is then used to estimate the actual values of τ at over 80 sites with eddy covariance data from the FLUXNET program (Pastorello et al., 2020).Finally, we investigate which site characteristics (climate factors and ecosystem types) best predict τ.The questions we seek to answer include: is the timescale of stomatal closure adapted to the dynamics of rainfall?If so, how?Is a particular ecosystem's stomatal behavior calibrated more to save water during a long but rare drought, or more to maximize carbon gain over a shorter but more frequent dry period?Can the proposed model be reconciled with other physiologically based approaches, thus offering a new perspective on the concept of discount factor timescale?

Overview of Methods
Our approach is divided into two parts.First, a new theoretical model of stomatal optimality is developed, whose solution is explicitly expressed in terms of its timescale factor τ (Section 2.2).Second, the model's τ value is fit to eddy covariance data in a variety of ecosystems (Section 2.3).In the proposed optimization framework, the plant seeks to maximize carbon gain over a drydown (a period with no precipitation), integrated over time with a discount factor, while avoiding drawing soil moisture below a critical value.The plant uses stomatal conductance as a control variable that modulates both the rate of carbon gain and the rate of soil moisture decrease.While other formulations are possible for the control variable (e.g., leaf water potential), the choice here is in keeping with standard optimality schemes.
After the model has been formulated, an empirical approach is developed to estimate τ from eddy-covariance and meteorological data.Since the study uses eddy-covariance data, the empirical τ values obtained are representative of entire ecosystems, aggregating the water use strategies of multiple species.Finally, variation in τ between plant functional types and cross-site correlations between τ and long-term climate factors are assessed.

Optimization Framework
In the proposed optimality model, stomatal conductance depends on the evolution of soil water storage w.The water balance during a drydown is expressed in the following equation: Above, the quantity w refers to the total plant-available water in the soil (with units of mm or mol/m 2 ), rather than the volumetric soil moisture measured at any particular depth.It can be thought of as volumetric soil moisture vertically integrated over the rooting depth.The t is time; here calculated at a daily time step as this time step captures much of the variance in root-zone soil moisture (Ghannam et al., 2016;Nakai et al., 2014).
The g is stomatal conductance to water vapor (units of mm/day), L is leaf area index (LAI, which is used to convert stomatal-scale to canopy-scale conductance), and D is vapor pressure deficit (VPD, defined with respect to atmospheric pressure, so unitless).Note that Equation 1 is only technically accurate at the scale of a leaf rather than a whole canopy, since it neglects aerodynamic resistance and numerous other factors.However, this equation has a long tradition of use in ecohydrological modeling (Feng et al., 2018;Green et al., 2024;Manzoni et al., 2015).
Equation 1 does not contain any loss terms aside from transpiration, so it neglects the effects of bare soil evaporation and soil drainage/subsurface runoff.We assume that in vegetated ecosystems, those fluxes are small compared to transpiration.This assumption is made for simplicity, and to enable a sufficiently simple model that can be tractably optimized for τ.As discussed further in Section 2.3.2,days with rainfall are excluded, as well as the first day following each rainfall event.Drainage is likely to be highest during these days, so removing them minimizes the error associated with this assumption.Bare soil evaporation is expected to be minimal in areas with high LAI, but it may be less negligible in more sparsely vegetated ecosystems; this issue is further discussed in Section 4.3.
The optimization criterion is represented by the following equation, where a stomatal conductance function g(w) that maximizes the reward function L(g,w) is sought:

AGU Advances
10.1029/2023AV001113 The exponential term within the integral is a discount factor that decays over time, with τ being the time constant.This term represents how much a plant values carbon gain in the future relative to gain in the present (or equivalently, saving water for the future relative to using it in the present).For example, if τ has a value of 10 days, the plant values one unit of carbon assimilated 7 days in the future about half as much as the same unit of carbon today, since e 7/10 ≈ 0.497.Equation 2 itself does not presuppose a link between τ and any externally forced environmental variable.However, once Equation 2 is coupled with the water balance of Equation 1, the lack of precipitation input in Equation 1 leads to an interpretation of τ as a timescale between rainfall events.The quantitative details of exactly which rainfall-related timescale τ corresponds to in actual ecosystems are the focus of the present study.
The function A(g(w)) within the integral in Equation 2 represents the carbon gained as a function of opening stomata.We assume that this is an increasing function of g that approaches an asymptote as g increases (i.e., there are diminishing returns to opening stomata due to carboxylation capacity saturating as internal plant carbon dioxide concentration increases).The functional form of A(g(w)) is described further in Section 2.2.2.
The last term within the integral in Equation 2 is an indicator function (I) of w, which is equal to 1 if w is above a critical value w c and 0 if w is equal to or below w c .Thus, the plant cannot gain any reward if soil moisture goes below the critical value.This term is equivalent to enforcing the constraint that the plant cannot draw down the soil water w below w c , and it forces the solution function g(w) to approach 0 as w approaches w c from above.The critical soil water value can be interpreted as an ecosystem-scale wilting point; its value in a given ecosystem will depend on soil texture and plant hydraulics.This wilting point is a proxy for multiple phenomena found at low soil moisture and correspondingly low leaf and xylem water potential that introduce tradeoffs between current and future carbon gain, including but not limited to xylem embolism and other plant hydraulic processes.The variety of such effects is discussed by Feng et al. (2022) and by Potkay and Feng (2023a).In this study, we do not attempt to set w c a priori; rather we estimate τ in a way that does not depend on knowing w c , as our focus is on τ.This estimation strategy is described in Section 2.3.1.
In summary, Equation 2 expresses a plant maximizing photosynthetic carbon gain while simultaneously avoiding the danger of low soil moisture, all while taking into account a discount factor time scale.Equation 2 is similar to the optimization criterion from Mäkelä et al. (1996).However, their model accounts for changes in stomatal conductance at both diurnal and multi-day timescales, and requires numerical integration to fully solve.Here we focus only on multi-day timescales to allow an analytical solution that we use to estimate the value of τ from observations of ecosystem fluxes.As earlier noted, daily time scale variations in w capture well over 90% of the soil moisture spectrum when sampled over half-hour time scales (Ghannam et al., 2016;Nakai et al., 2014).
The optimal stomatal conductance (assumed here to capture the stomatal conductance dynamics of the ecosystems studied) is the function g(w) that maximizes Equation 2 when the dynamics of soil moisture are governed by Equation 1.As shown in Supporting Information S1 via the calculus of variations, this optimal conductance must satisfy the following differential equation: In Equation 3, A′(g) and A″(g) are the first and second derivatives of carbon gain with respect to stomatal conductance.To find an explicit expression for g(w), we therefore need to determine these derivatives.
The proposed optimality model is focused on soil moisture as the state variable that is most directly affected by stomatal closure.We do not account for stochastic day-to-day variations in VPD and light environment, or the effect of stomatal closure on VPD through land-atmosphere coupling.Thus, the solution in Equation 3 is only strictly optimal if VPD and other environmental conditions (aside from soil moisture) are constant over the entire drydown.Nevertheless, we use this constant-conditions assumption as a plausible scenario for the plant's optimization behavior: because the future evolution of VPD over the rest of the drydown is unknown at any given time due to its dependence on external atmospheric conditions.Assuming that VPD will remain similar to its current level for the rest of the drydown provides a reasonable guess for the plant to use in determining how much to close stomata.That is, the model implemented at the daily scale considers the plant as updating its guess of future VPD every day to equal the present VPD and optimizing based on that guess.A similar technique of treating daily-scale environmental conditions as constant within the optimization is used in Mäkelä et al. (1996).

A Photosynthetic Model
In stomatal optimality models, the A(g) function is typically based on a biochemical model of photosynthesis, such as the widely used Farquhar model (Farquhar et al., 1980).However, the parametric form of such a function along with the multiple constraints on photosynthesis (e.g., RuBP, Rubisco) makes it impossible to obtain a closed-form solution for g(w) once the derivatives of A(g) are inserted into Equation 3. Here, a saturating exponential function (Equation 4) is used to model A(g): The form of Equation 4was chosen for several reasons: (a) its concavity guarantees an optimum for finite g, (b) its derivatives allow Equation 3 to be solved analytically, and (c) it preserves the shape of more detailed biophysically based models (Figure S2 in Supporting Information S1).In Equation 4, A max represents the maximum possible carbon gain, and g A represents the stomatal conductance at which A is about 63% (or 1 e 1 ) of its maximum potential value.Equation 4 describes a smooth curve, unlike the full Farquhar model that has an abrupt transition between RuBP limitation and RuBisCO limitation.In this respect Equation 4is closer to the Collatz model which smooths over the Farquhar model's transitions (Collatz et al., 1991;Walker et al., 2021).Besides its mathematical simplicity, the use of a smooth function in photosynthesis modeling can be motivated by the tendency of leaf-to-canopy scaling to smooth out any sharp transitions in the underlying biochemical dynamics (Katul et al., 2007).
While the baseline levels of A max and g A in an ecosystem are largely set by the availability of soil nitrogen and its uptake by plants (Ollinger et al., 2008;Palmroth et al., 2013), the temporal variation of these parameters depends on changing environmental conditions, most importantly incoming photosynthetic radiation (as well as temperature).The exact nature of that dependence should be ecosystem-specific because ecosystems have different photosynthetic parameters and different canopy geometries that affect light attenuation in the canopy (He et al., 2019;Rogers et al., 2017).Here, both A max and g a are parametrized as Michaelis-Menten functions of incoming shortwave radiation R, with each function having a saturated value of V and a half-saturation point of K: This functional form is one of several that has been used for photosynthetic light-response curves in the literature (Hari & Makela, 2003;Marino et al., 2010).

Optimality Solution With the Simple Photosynthetic Model
When substituting the first and second derivatives of Equation 4 into Equation 3, the A max term cancels.After solving the resulting ODE, an explicit equation is found for the optimal stomatal conductance: This equation predicts that stomatal conductance should scale with the inverse square root of VPD, a relationship that has been found in optimality models since the work of Hari et al. (1986), and which was analyzed in detail in Katul et al. (2009) in terms of its similarity to other empirical VPD-g relationships such as those discussed in Oren et al. (1999).It also predicts a positive (and thus realistic) response of stomatal conductance to increasing soil moisture.Additionally, the intrinsic water use efficiency (iWUE, that is, carbon gained per unit of conductance) can be calculated by combining Equations 4 and 7.As shown in Figure S1 in Supporting Information S1, the
proposed model predicts iWUE to increase with decreasing plant available water, as found in previous studies (Battipaglia et al., 2014;Vickers et al., 2012;Yi et al., 2019).
The presence of g A in Equation 7 reflects how, based on Equation 4, g A sets the scale of conductance at which diminishing returns to stomatal opening become important.For example, if g A is 0.05 mol/m 2 /s, then an increase in stomatal conductance from 0.1 mol/m 2 /s to 0.15 mol/m 2 /s causes carbon assimilation to increase by only 9.9%, while if g A is 0.1 mol/m 2 /s the same change in stomatal conductance results in a 23% increase in assimilation.Thus, under the larger g A value, there is more incentive for the plant to open stomata further, which results in a higher optimal conductance.Since g A is an increasing function of incoming light (Equation 6), so is the predicted stomatal conductance.Furthermore, it can be shown that Equation 7 can be made equivalent to general expressions for optimal conductance under Rubisco limitations presented in Katul et al. (2010), if g A is suitably redefined.This equivalence is arrived at when the Katul et al. (2010) general stomatal conductance model is integrated over the length of a single drydown (or unknown duration) without any hydrologic balance necessary beyond initial and final soil moisture.The equivalence between the two models is completed by enforcing the initial optimal conductance at the start of a dry down in the earlier model to g A .An explanation of this equivalence is presented in Supporting Information S1.
To demonstrate the effect of τ in the overall model, Figure 1 plots the predictions of Equation 7 coupled with the water balance from Equation 1, with two example values of τ, compared to eddy covariance observations from two sites.For each site, a single example drydown is shown.To isolate the effect of changing τ in this illustrative figure, the values of initial soil moisture above w c are set so that the initial modeled ET is equal to the observed value.For both sites shown, when τ = 20 days, modeled ET decreases more quickly (i.e., ET changes more relative to the same change in soil moisture) than when τ = 50 days.Figure 1 also shows that the observed ET in this example drydown for US-Santa Rita Mesquite (SRM), a savanna in Arizona, exhibits behavior similar to that of the model with τ of 20 days.By contrast, US-Me5, a coniferous forest in Oregon, matches the behavior associated with larger τ values.It can also be shown (see Supporting Information S1) that the time it takes to deplete a fixed initial soil moisture down to w c is proportional to the square root of τ (in addition to being related to VPD and the other quantities in Equation 5).This result illustrates the role of τ as the timescale of vegetation water-saving.
A comparison to the instantaneous optimization Medlyn model-which is well validated and commonly usedillustrates that the response of the proposed optimization model to VPD and soil moisture is realistic.If the photosynthetic model used here (Equation 4) is substituted into the Medlyn model, the dependence of optimal g on VPD (holding plant-available water constant) is similar to that of the proposed model (Figure 2a).The Medlyn model does not prescribe a specific stomatal response to soil moisture, but such a response is often included by multiplying the model's g 1 parameter by a so-called "β function" of soil moisture (Zhou et al., 2013).We can fit a β(w) function that, when it is used in the Medlyn model, produces a response to soil moisture (holding VPD constant) that is similar to the square-root dependence of the proposed model in Equation 5 (Figure 2b).In this example, β(w) is parametrized as a power function of linearly transformed soil moisture, as in several other studies (Drake et al., 2017;Egea et al., 2011;Keenan et al., 2010).This power function reflects the fact that soil water tension, not moisture, regulates the driving force for the delivery of water from the soil to roots.The combined dependence on soil moisture and VPD is shown in Figures 2c-2f, again revealing the similarity between the predictions of the time-integrated model and the Medlyn model, the latter of which has been empirically tested on numerous data sets.

Regression Strategy to Estimate τ
Estimating τ from Equation 7 requires quantifying the soil water storage s (as well as the photosynthetic potential factor g A , which we estimate from eddy covariance data as described in Section 2.3.2).Although many FLUXNET sites measure soil moisture, these measurements are typically made only close to the surface, so they are not representative of the root zone and cannot be used to estimate s.Instead, we took an approach to estimate τ that does not require measured soil moisture but rather infers it on a drydown-specific basis.We defined a drydown as a period of at least 5 consecutive days with zero recorded precipitation.We excluded the first day of each drydown (the day immediately following rainfall) to avoid periods where runoff and/or canopy evaporation might represent a substantial part of total ecosystem water fluxes.
To estimate τ, we began by noting that the right-hand side of Equation 1 equals ET, and by combining Equations 1 and 7 to derive an equation for ET: Additionally, by integrating both sides of Equation 1, soil water storage during each drydown can be related to the cumulative sum of ET (ignoring other losses of soil moisture, an assumption already made in Equation 1): where i indexes time in days during the drydown.It is important to note that, to predict ET on a given day, only the sum of ET up to the previous day is used, thus avoiding circularity in the optimization.
Rearranging Equation 8 and plugging in Equation 9, the following equation is obtained: The left-hand side of this equation is effectively a version of ET that has been normalized to account for the effect of VPD, LAI, and photosynthetic potential, leaving soil water as the only environmental driver whose effect remains on the right-hand side.
To derive a version of Equation 10 that can be used to estimate τ, we start by defining the normalized ET, ET norm , and accumulated ET, ET acc : Then Equation 9can be re-written as follows: Equation 13 only applies to drydowns where ET is water-limited.To determine which drydowns are waterlimited, we calculate the correlation between ET norm and ET acc within each drydown.If this correlation is negative, then as ET draws down soil moisture, ecosystem conductance controlling for VPD and photosynthetic capacity also decreases (consistent with a positive τ).For drydowns that had a positive correlation between ET norm and ET acc , ET was inferred to be insensitive to soil water availability.Those non-water-limited drydowns were excluded from the rest of our analysis.
To estimate τ, Equation 13 can be interpreted as a linear model, with a slope a = 2/τ and a different intercept term for each drydown, b dd = 2 τ (w 0 w c ) : ET norm (t,dd) = a ET acc (t,dd) + b dd (dd) (14) In addition to estimating τ as 2/a, the coefficients of this regression model also yield estimates of the quantity w 0 w c = τ 2 b dd for every drydown.The right-hand side of Equation 14 then provides estimates of w throughout each drydown.As a reasonableness check, we then calculate the correlation between these estimates of soil water storage and the measured surface soil moisture from the FLUXNET data, over all the drydowns that were used to fit Equation 14 at a site.The observational soil moisture data was not used anywhere else in the analysis, so it can be treated as a completely independent data set for comparison.Finally, a 95% confidence interval for τ was calculated using the standard errors of the slope a and the inverse relation between a and τ.
The fitted regression model provides predictions of ET norm based on ET acc .These predictions can be transformed to predictions of ET itself using Equation 11.Although the goal of our study is to characterize τ rather than to develop a full predictive model of ET, checking the proposed model's ET predictive skill is still useful to ensure the model's realism.However, simply comparing predicted ET to observed ET may give overly optimistic accuracy metrics due to the structure of Equation 14-each drydown has its own fitting parameter in the form of its initial soil moisture, so there are many degrees of freedom.To reduce this problem, we instead divide both observed and predicted ET within each drydown by the drydown's initial observed value (ET 0 ) before calculating an R 2 of ET/ET 0 for all drydowns together for a site.Thus, the R 2 values we present represent the model's skill at predicting changes in ET over the course of multiple drydowns, not for predicting absolute values of ET.

Use of Eddy Covariance Data
Data on meteorological conditions and ecosystem fluxes from the FLUXNET2015 compilation of eddycovariance sites (Pastorello et al., 2020) was used in the present study.We excluded sites that were labeled as wetlands or as snow/ice under the IGBP biome classification, since the proposed optimization model is not valid in sites where water input is dominated by sources other than rainfall.Due to the small number of sites classified as "woody savanna" and "closed shrubland," we grouped those IGBP classes with "savanna" and "open shrubland" respectively, when analyzing the results by biome.For evapotranspiration (ET), we used the latent heat flux data field (not the energy balance corrected version), and for GPP we used the average of the nighttime and daytime estimates provided with the FLUXNET2015 data.Additionally, we used LAI data derived from the MODIS satellite: the 500 m Terra/Aqua 8-day composite MCD15A2H product (Myneni et al., 2021).For each site, LAI data was temporally smoothed according to the methods in Ukkola et al. (2022).It should be noted that our method for estimating τ does not depend on the absolute magnitude of LAI at a site, but only on the temporal dynamics of LAI.
We used flux data during the growing season only, to avoid time periods when rapidly changing LAI was the main driver of changing ET, as opposed to water limitation being the driver.This assumption also reduced errors associated with possible imperfections in the LAI data set.We defined the growing season based on climatological average LAI.After finding the maximum value of the LAI climatology and the "peak day of year" at which that value is achieved, we defined the start of season as the last day before the peak day that had an LAI of less than 75% of the maximum value.Similarly, we defined the end of season as the first day after the peak day that had an LAI of less than 75% of the maximum value.Thus, the growing season is defined as a contiguous range of days with high LAI, roughly centered around the peak day.To ensure data quality, in the remaining analyses we further limited ET and GPP to days on which daytime VPD was higher than 0.5 kPa, daytime shortwave radiation was at least 100 W/m 2 , and when no rainfall was recorded on that day or on the previous day (for examples of similar quality filters on eddy covariance data, see Knauer et al., 2018;Novick, Ficklin, et al., 2016).
To estimate g A -a prerequisite to estimating τ-we fit Equations 4-6 to eddy covariance data at each site separately.We estimated hourly canopy conductance by dividing eddy covariance ET by VPD, for consistency with Equation 1.We then divided by LAI to convert the estimated stomatal conductance from canopy to stomatal scale.The coefficients V A , K A , V g , and K g in Equations 5 and 6 were estimated simultaneously by using nonlinear

10.1029/2023AV001113
HOLTZMAN ET AL. optimization to minimize the squared errors in predicted hourly GPP during the daytime.Finally, predicted hourly g A based on Equation 6was averaged to a daily timescale as an input for estimating τ.
Out of 120 FLUXNET sites (excluding snow/ice and wetlands), there were 91 sites where there was enough data to fit the model of growing-season GPP using Equations 4-6 (given observed canopy conductance, shortwave radiation, and LAI), and for which that model had an R 2 value above 0. Of those sites, 82 sites satisfied the following criteria to produce a reliable estimate of τ: • There were at least three drydowns that had a negative correlation between ET norm and ET acc (water-limited drydowns) • The site's 95% confidence interval for τ, from Equation 14, did not overlap 0 and had a width of less than 100 days These 82 sites were used in the rest of the analyses, and they are shown on the map in Figure 3.

Relating τ to Climate Characteristics
Across FLUXNET sites, the estimated values of τ were compared to various site climate characteristics to investigate what factors might drive ecosystems to evolve different τ values.The following climate characteristics were used as possible predictors of τ: • Mean annual precipitation Aridity index was calculated as net radiation times the latent heat of evaporation of water divided by precipitation.
For the last two variables in the above list, we treat the end of each year's growing season as the end of a drydown even if there is no rain that day.This step is necessary to exclude drydowns that start near the end of the growing season and last long after the growing season is over.All site climate metrics were calculated using the FLUXNET 2015 data.

Estimates of τ at FLUXNET Sites
Across 82 sites, the estimated τ values ranged from 5.0 to 122 days, with a mean of 21.6 days.For comparison, Berninger et al. (1996) obtained a value of 24.6 days for the timescale parameter of a similar stomatal model that was fit to data from a pine forest in southern Finland; our analysis found a τ value of 20.0 days for the same site.Supporting Information S2 contains a full table of τ values and other calculated quantities for all 82 sites.The number of water-limited drydowns per site ranged from 4 to 83, with a mean of 25.1.The number of unique years per site with water-limited drydowns ranged from 1 to 14 with a mean of 6.6 years.The proportion of all drydowns that were water-limited varied by site from 44% to 100% with an average of 74%; there were 78 out of 82 sites where more than half of all drydowns were water-limited.The 16 sites that were broadleaf forest or mixed forest had a lower proportion of water-limited drydowns (63% vs. 76%) compared to the sites in other biomes.
The relative standard error of the regression slope a (Equation 14) ranged from 3.2% to 44%, with a mean of 15%.Calculated from these standard errors, the width of the 95% confidence interval for τ ranged from 1.6 to 69 days, with a mean of 14 days.Sites with more water-limited drydowns in the analysis had more precise estimates of a, and thus of τ.For example, the 34 sites with at least 25 water-limited drydowns had a mean relative standard error for a of 10.0%, compared to 18.4% for the other 48 sites with fewer than 25 drydowns to analyze.This result suggests that longer data records, if they were available, would help reduce the uncertainty in our τ estimates.
As a plausibility check on the underlying assumptions used in the estimation of τ, we assessed the fit of the GPP model, as well as the fit of the proposed optimality model's predicted ET (normalized by drydown-initial ET) and soil moisture.Across sites, the GPP model had R 2 values with an interquartile range of 0.29-0.69and a mean of 0.48 (Figure S3 in Supporting Information S1).The obtained R 2 values were interpreted as acceptably high for use of the simplified GPP model in Section 2.2.2 to be reasonable.The R 2 values for ET/ET 0 (calculated according to Section 2.3.1) had an interquartile range of 0.45-0.82and a mean of 0.53; 74 out of 82 sites had a positive R 2 (Figure S3 in Supporting Information S1).The correlations between plant-available water inferred from Equation 14 and measured surface soil moisture at each site had an interquartile range between 0.36 and 0.73, with a mean of 0.51 (Figure S3 in Supporting Information S1).These correlations have a similar range to correlations between surface and 50 cm-deep soil moisture at field sites that measure both quantities (Feldman et al., 2023).This result suggests that using Equation 14 to simultaneously estimate soil water storage and τ is valid.

Typical Drydown Length Explains Spatial Variability of τ
The sites we used in the final analysis span a range of annual average temperatures from 2.7 to 27 degrees Celsius, and average annual precipitation ranging from 17 to 200 cm (Figure S4 in Supporting Information S1).
We found no significant differences in τ values between biomes, according to an ANOVA with p > 0.5.By contrast, τ is significantly and positively correlated with growing season D max , with R 2 = 0.59 (Figure 4).Additionally, as shown in Figure 5, growing season D mean was approximately as highly correlated with τ (R 2 = 0.54) as growing season D max was, likely because of the high correlation between the two drydown length metrics (R 2 = 0.94).The other climate metrics tested all had an R 2 of less than 0.2 for predicting τ (Figure 5).The magnitude of τ is much closer to the growing season D max than growing season D mean .Based on linear regression with the intercept fixed at 0, the slope between τ and growing season D max is close to 1 (slope = 1.08), whereas it is much greater for D mean (slope = 2.00).This result suggests that ecosystem stomatal behavior tends to be adapted to the longest dry period in a typical year.To test whether the correlations in Figures 4 and 5 are strongly influenced by outliers, we repeated that analysis while excluding sites with the two largest τ values and the two largest values of each climate variable (Figure S5 in Supporting Information S1).Although the correlations with large values removed were lower than they were with the full data set, the results were qualitatively similar to Figures 4 and 5.

Inferred τ Values in the Context of Plant Hydraulic Strategy
The spatial variability in τ, the timescale with which carbon gain is discounted for stomatal closure, was found to correlate positively with D max , the typical duration of a year's longest drydown.That is, in climates with longer maximum dry spells during the growing season, plants value future carbon assimilation more-and save water over a longer term-than in climates where rain is frequent and the longest dry spells are short.The similarity between our model and that of Katul et al. (2010) emphasizes the physiological connection between τ and typical drydown length.Although the Katul et al. (2010) model was originally derived without time integration, the two model solutions become equivalent when the Katul et al. (2010) model is integrated over the length of a drydown.This equivalence is found despite the Katul et al. (2010) model taking into account detailed leaf physiology, while the model in the current study is based on a simplified canopy-level parametrization of photosynthesis.
We found that τ is also significantly correlated with D mean .On average, τ is equal in magnitude to D max , and considerably larger than D mean , suggesting a potential (though not certain) adaptation to D max rather than to D mean .This suggests relatively conservative stomatal behavior, as has been previously hypothesized to be evolutionarily optimal for preventing dryness-induced damage to plant tissues, as discussed in detail in Bartlett et al. (2016) and Potkay and Feng (2023a).Recall that dryness-induced limitations-which impact the plant through water transport limitations and embolism-are accounted for in our model via the wilting point parameter w c as an indirect proxy.The conservative stomatal behavior that is suggested by our results might be favored because vegetation recovery from drought is not instantaneous.Xylem that is damaged by embolism takes time to repair (if it is repairable at all), during which the plant's capability for photosynthesis is reduced (Brodersen & McElrone, 2013;Kannenberg et al., 2019;McDowell et al., 2019).In addition to xylem embolism, other mechanisms of long-term damage to plants from water stress include leaf shedding and reduced growth rates (McDowell et al., 2008).It may be advantageous for plants to save water for a longer drydown than the mean drydown length, because doing so reduces the likelihood of damage that would impair future carbon assimilation.We also note that in Equation 2, τ is an e-folding timescale and not an absolute cutoff.Thus, our finding that τ is on average twice the magnitude of D mean could also be interpreted without reference to D max , in terms of plants valuing carbon an average drydown-length in the future 14% as much as carbon in the present.
The plant's water-saving timescale (τ) studied here can be thought of as a functional trait that is part of a larger water-use strategy (Kannenberg et al., 2022).Indeed, there is extensive evidence of evolutionary coordination between traits such as plant height, rooting depth, xylem hydraulic conductivity, xylem vulnerability curve parameters, and biochemical photosynthetic capacity (X.Li et al., 2018;H. Liu et al., 2019;Y. Liu et al., 2021;Sanchez-Martinez et al., 2020;Xu et al., 2021).Analogously to our findings with τ, some of these traits have been observed to vary systematically with local climate factors.Among seedlings of one eucalyptus species taken from different locations in Australia and grown in the same greenhouse, C. Li et al. (2000) found that hydraulic architecture traits and transpiration levels were correlated with the average amounts of dry-season rainfall in the seedlings' locations of origin.Correlations between climate dryness metrics and hydraulic trait values have also been found among European beech trees intraspecifically (Schuldt et al., 2016), across multiple tree species over aridity gradients in Australia (X.Li et al., 2018;Pfautsch et al., 2016), and across 150 woody legume species worldwide (H.Liu et al., 2024).
Future work could investigate whether τ coordinates with other hydraulic traits as part of a whole-plant strategy or whether it varies independently from them.An implication of our results is that, in the absence of changes to stomatal behavior, if a location's value of D max changes due to changing precipitation patterns, plants at that location will no longer be functioning optimally.Since drydown durations are lengthening across much of the globe (Fischer et al., 2013), this may put plants under greater water stress than their stomata have evolved to cope with, leaving plants vulnerable to hydraulic damage.Such dis-equilibrium, however, might be partially mitigated by changing species cover (García-Valdés et al., 2021;Trugman et al., 2020) or by plasticity in water use traits (Blackman et al., 2017;Challis et al., 2022), although the effect of these processes on τ is unclear.
Although our analysis is only focused on the growing season, there is evidence that phenology is interrelated with stomatal behavior.For example, Manzoni et al. (2015) found a relationship between plant isohydricity and leaf phenology in contrasting drought regimes.Additionally, climate seasonality can affect which ecohydrologic strategies are most productive (Vico et al., 2015), much like rainfall intermittency affects the optimal τ here.We can revisit the two illustrative sites shown in Figure 1 in light of their climates.US-SRM has a slightly drier climate than US-Me5 overall: annual mean precipitation values of 0.91 mm/day and 1.1 mm/day respectively, and annual D max values of 67 and 34 days respectively.Thus, it may be surprising that US-SRM has a much lower value of τ than US-Me5 (15 vs. 37 days).The site in the wetter climate seems to be adapted to longer periods without rain.However, the intra-annual patterns of precipitation and phenology at the two sites are different.At US-SRM, the main growing season is in phase with the main rainy season (the summer monsoon), while at US-Me5, the growing season is during the dry summer after the winter precipitation has already fallen.During each site's growing season, US-SRM has more frequent rainfall than US-Me5, which explains the difference in τ values from a climate adaptation perspective.

Implications for Stomatal Modeling
The best way to prognostically integrate the stomatal model developed here into land surface models or other large-scale models is unclear, as well as whether such an effort would even be valuable.Several simplifications are made in our treatment of water balance and photosynthesis (see Section 4.3) which are not acceptable in a full land surface model.Thus, in the context of a model with dynamically updating soil moisture and without these simplifications, the conductance in Equation 7 may no longer be the optimal solution for the optimization criterion in Equation 2. Indeed, the primary aim of our study was not prognostic but rather diagnostic: not to develop a new model for stomatal conductance, but to understand spatial variations in the optimality timescale.Nevertheless, the value of accounting for non-instantaneous considerations in stomatal closure may outweigh the simplified aspects of the proposed optimality model, such that progress toward incorporating a similar modeling approach in a land surface model could prove beneficial.Thus, as a next step, we recommend that the proposed time-integrated stomatal model be compared to existing stomatal models in terms of predictive accuracy and sensitivity to environmental drivers over longer time scales, as in Sabot et al. (2022) or Bassiouni and Vico (2021).
The strong relationship between τ and D max (Figure 4) also motivates the parametrization of stomatal closure in models as a function of local climate, not just using plant functional types, as is commonly done.This approach is known as environmental filtering (e.g., Famiglietti et al., 2023;Walker et al., 2017).Our findings are consistent with previous findings of instantaneous stomatal parameters being related to aridity (Bassiouni et al., 2023;Lin et al., 2015;Y. Liu et al., 2022) and rainfall intensity and frequency (Bassiouni et al., 2023).The use of environmental filtering in stomatal closure model parametrization may improve ET estimation, as previously shown across the globe using a simplified Shuttleworth-Wallace model (Wu et al., 2020).Efforts to incorporate the effect of local climate on stomatal parameters would likely benefit from an integrated study of the coordination in climate sensitivity between multiple different stomatal traits.

Limitations of This Study
The optimality model proposed here does not account for competition between plants for soil water, which can have substantial effects on optimal plant water use (Cowan, 1986;Lu et al., 2020;Mrad et al., 2019).More competitive ecosystems would be expected to have smaller τ values, as individual plants would use soil water quickly to prevent their neighbors from using it first.Future work could study whether Equation 7 represents an evolutionary stable strategy in the context of competing plant populations.Also, the model does not account for diurnal changes in plant water status or differences between soil and leaf water potentials, so it cannot represent all aspects of plant hydraulic strategy.This exclusion of internal plant hydraulics could lead to inaccuracy in the relative strengths of stomatal conductance responses to VPD and soil moisture (Y.Liu et al., 2020).Another limitation of our study is our treatment of all ET as transpiration.To simplify the optimality model, we ignored other sources of evaporation, such as soils, water bodies, or canopy interception.This assumption is most reasonable in areas where dense canopy cover minimizes the fraction of ET that is bare soil evaporation.Our results for more sparsely vegetated ecosystems (savannas, shrublands, and grasslands) should be taken with caution.Even in forests, soil evaporation can make up as much as 25% of total ET (X.Li et al., 2019;Paul-Limoges et al., 2020;Qubaja et al., 2020).Indeed, the R 2 values for predicted ET/ET 0 (values shown in Figure S2 in Supporting Information S1), were positively correlated with site average LAI (r = 0.37, p < 0.01, excluding sites with negative R 2 values), suggesting that the proposed optimality model performs better in more densely vegetated sites.
In our comparison between climate factors and water-saving timescale, several sites show values of τ greater than the growing season D max .In these ecosystems, plants may typically "prepare for" (i.e., they close stomata in a manner consistent with the possibility of) a more extreme drought that might occur only once every few years.However, some of the scatter in Figure 4 is simply due to noise in our estimates of τ and D max , so it is unclear if these sites reflect different plant behavior or simply noise.Sources of uncertainty in our τ estimates include the limited number of drydowns at sites with wet climates, noise in the ET observations and GPP estimates, and processes such as runoff that are not captured by the proposed model.A possible source of bias in our τ estimates is that we use VPD measurements from the tops of eddy covariance towers located above the canopy, rather than within the canopy where stomatal closure and photosynthesis are occurring.Uncertainty in our D max estimates is largely due to the limited number of years in each site's FLUXNET data.

Conclusions
We introduced a new stomatal optimality model that describes the dynamic feedback between soil moisture and stomatal conductance, and parametrizes integrated carbon gain over time using a discount factor with a timescale of τ.The model treats photosynthesis and water balance in a way that maintains realism while enabling a closedform equation for optimal stomatal conductance as a function of soil moisture, VPD, and incoming shortwave radiation, with τ as a parameter.We then fit that equation to FLUXNET data to estimate values of τ in a range of ecosystems.Across FLUXNET sites, τ was found to be similar in magnitude to a particular climate-related timescale: the length of the longest period without rain in an average year.This result helps characterize what level of drought plant water use tends to be adapted to.Our findings could be used to parametrize land surface models, and combined with climate change projections to understand where terrestrial carbon and water cycling may be particularly vulnerable or resilient to changing drought occurrence.

Figure 1 .
Figure 1.Optimality model predictions of ET during example drydowns from two FLUXNET sites under two different illustrative values of τ (blue and orange lines), compared to observations (black line) at each site.The sites are US-Santa Rita Mesquite in southeastern Arizona, USA (top of panel); and US-Me5 (Metolius first young-aged pine) in central Oregon, USA (bottom).The x-axis starts at day 2 of each drydown because we exclude the day immediately after rainfall.

Figure 2 .
Figure 2. The new time-integrated optimality model and the Medlyn model show similar responses to vapor pressure deficit (VPD) and soil water if the Medlyn model g 1 parameter and soil β function are calibrated to match the optimality model.The stomatal conductance values predicted by time-integrated and Medlyn models are denoted by g INT and g MED , respectively.(a) Shows stomatal conductance predicted by both models as a function of VPD, with soil moisture held constant.(b) Shows stomatal conductance predicted by both models as a function of soil moisture, with VPD held constant.(c) Shows the timeintegrated model predicted conductance as a function of both VPD and soil moisture.(d) Shows the Medlyn model predicted conductance as a function of both VPD and soil moisture.(e) Shows the difference between (c) and (d).(f) Shows a scatter plot between the values in (c) and the values in (d).

Figure 3 .
Figure 3. Locations of FLUXNET sites used in the final analysis.

Figure 4 .
Figure 4. Across ecosystems, the water-saving timescale τ is similar in magnitude to the growing season D max (average annual-maximum dry period length).

Figure 5 .
Figure 5. Relationships between τ and climate characteristics.Aridity index is abbreviated as AI, and growing season is abbreviated GrowSeas.