Responses of Marginal and Intrinsic Water‐Use Efficiency to Changing Aridity Using FLUXNET Observations

According to classic stomatal optimization theory, plant stomata are regulated to maximize carbon assimilation for a given water loss. A key component of stomatal optimization models is marginal water‐use efficiency (mWUE), the ratio of the change of transpiration to the change in carbon assimilation. Although the mWUE is often assumed to be constant, variability of mWUE under changing hydrologic conditions has been reported. However, there has yet to be a consensus on the patterns of mWUE variabilities and their relations with atmospheric aridity. We investigate the dynamics of mWUE in response to vapor pressure deficit (VPD) and aridity index using carbon and water fluxes from 115 eddy covariance towers available from the global database FLUXNET. We demonstrate a non‐linear mWUE‐VPD relationship at a sub‐daily scale in general; mWUE varies substantially at both low and high VPD levels. However, mWUE remains relatively constant within the mid‐range of VPD. Despite the highly non‐linear relationship between mWUE and VPD, the relationship can be informed by the strong linear relationship between ecosystem‐level inherent water‐use efficiency (IWUE) and mWUE using the slope, m*. We further identify site‐specific m* and its variability with changing site‐level aridity across six vegetation types. We suggest accurately representing the relationship between IWUE and VPD using Michaelis–Menten or quadratic functions to ensure precise estimation of mWUE variability for individual sites.


Introduction
Terrestrial plants mitigate global warming by sequestering atmospheric carbon dioxide (CO 2 ) through photosynthesis (Beer et al., 2010).However, photosynthesis is inherently linked with plant water loss via transpiration, as CO 2 and water vapor share the same stomatal pathway.Plants risk hydraulic damage during droughts if they maintain high stomatal conductance as soil water availability decreases and atmospheric demand increases, resulting in low leaf water potential and xylem cavitation.Therefore, plants must balance stomatal function to optimize carbon uptake while minimizing transpirational water loss and hydraulic stress (Cowan & Farquhar, 1977;Katul et al., 2010;Sperry et al., 2017;Wang et al., 2020).To predict plant ecophysiological responses to projected changes in atmospheric CO 2 concentration, elevated atmospheric water demand, and more severe and frequent drought events, we need a mechanistic understanding of how different ecosystems regulate the trade-off between photosynthetic carbon assimilation and transpirational water loss (See Table 1 for a glossary of terms related to water-use efficiency used throughout this paper).
Although carbon uptake is usually represented through mechanistic models of photosynthesis (e.g., the Michaelis-Menten equation (Marshall & Biscoe, 1980;Michaelis & Menten, 1913;Thornley, 1976); the Farquhar model (Farquhar, von, et al., 1980;Von Caemmerer, 2000)), water use (i.e., transpiration) is often described based on empirical relationships that prescribe how stomatal conductance responds to environmental drivers and carbon uptakes.For example, the Ball-Berry model (Ball et al., 1987) is one of the most widely used empirical stomatal conductance models (Anderegg et al., 2017;Buckley, 2017;Katul et al., 2010), and has been readily incorporated into many climate models (Bonan et al., 2014).It takes the form: where g s is stomatal conductance (mol m 2 s 1 ), A is carbon assimilation rate (μ mol m 2 s 1 ), c a is atmospheric CO 2 concentration (ppm), RH is relative humidity at the leaf surface, and g 0 and g 1 are empirically fitted parameters.To simulate the non-linear variation in g s with changing humidity, Leuning (1995) modified the Ball-Berry model by replacing relative humidity with a vapor pressure deficit (VPD) response function as follows: where Γ* is CO 2 compensation point for photosynthesis (ppm) and VPD 0 is the empirically determined coefficient, representing the slope of the relationship between g s and VPD.These empirical models are relatively simple, easy to use, and work well for well-watered conditions (Bonan et al., 2014).However, they have an incomplete grounding in physiological theory, leading to uncertainty when they are extrapolated to predict plant  (Franks et al., 2018;Knauer et al., 2015Knauer et al., , 2018;;Medlyn et al., 2012;Sabot et al., 2022).
An alternative way to enable the theoretical interpretation of leaf-level stomatal conductance models is to adopt the principle of stomatal optimization theory (Anderegg et al., 2018;Bonan et al., 2014;Katul et al., 2009Katul et al., , 2010;;Medlyn et al., 2012;Novick, Miniat, & Vose, 2016;Sperry et al., 2017;Wolf et al., 2016).Stomatal optimization theory was originally based on a hypothesis that stomata are regulated to maximize carbon assimilation (A) for a given water loss (transpiration, E).A key parameter in this class of models is the so-called "marginal water-use efficiency (mWUE)," here defined as the ratio of a change in E to a change in A (∂E/∂A) following Cowan and Farquhar (1977), although it is sometimes defined as the inverse form (∂A/∂E) (Katul et al., 2010;Manzoni et al., 2011).The optimality models often maintain the mWUE constant over arbitrary time steps (e.g., daily), assuming abundant water at the canopy (Buckley, 2017;Cowan & Farquhar, 1977;Makela et al., 1996).However, this may not hold true at sub-daily timescales, where high atmospheric demand (i.e., VPD) during midday can decrease water potential at the canopy level even when soil moisture is abundant (Anderegg et al., 2017;Grossiord et al., 2020).
Understanding how mWUE changes under hydrologic stress is necessary for the optimization models in a prognostic sense, yet no consensus on the magnitude or even direction of these changes exists.For instance, Manzoni et al. (2011) and Zhou et al. (2013Zhou et al. ( , 2014) ) performed meta-analyses of leaf gas exchange measurements from previous studies that spanned wide ranges of species and moisture conditions.A major difference in their approaches was the proxy for plant water status; Manzoni et al. (2011) used mid-day leaf water potential, whereas Zhou et al. (2013Zhou et al. ( , 2014) ) used pre-dawn leaf water potential as a proxy for soil moisture availability.Similarly, Lin et al. (2015) compiled a global database of leaf gas exchange measurements spanning diverse plant functional types and estimated a slope parameter (g 1 ) (Medlyn et al., 2012), which is analogous to the slope parameter from empirical models (Equations 1 and 2) and proportional to Medlyn et al., 2012).They further evaluated the relationship between g 1 and a moisture index, defined as the ratio of mean annual precipitation to the equilibrium evapotranspiration.Mäkelä et al. (1996) and Lu et al. (2016) took a theoretical approach to examine shortand long-term optimal stomatal behavior, respectively, in response to the soil moisture availability assuming that plants are adapted to the stochastic rainfall patterns of their environments.More recently, alternative stomatal optimization perspectives have been proposed, which presume stomata function to maximize carbon uptake while minimizing water costs, including those linked to hydraulic damage during droughts (Anderegg et al., 2018;Sperry et al., 2017;Wolf et al., 2016).Although promising, in contrast to the Medlyn et al. (2012) model, these newer formulations have yet to be integrated into land surface model schemes (but refer to Kennedy et al., 2019, for a study implementing plant hydraulics in the Community Land Model).Although theoretical expectation and many studies indicate decreasing mWUE as water stress drives reductions to g s , there is some evidence of increasing mWUE under water stress (Farquhar, Schulze, & Kuppers, 1980;Grieu et al., 1988;Zhou et al., 2013), although reasons for this needed to be clarified.
It is also important to note that canopy water status and water potential are not determined solely by the availability of water supply but by the balance between water supply and demand, with VPD as a major force exerted on the canopy by the atmosphere (Manzoni et al., 2011(Manzoni et al., , 2013;;Novick et al., 2019).Thus, it is reasonable to expect that mWUE needs to be adjusted with changing atmospheric water demand unless other factors limit the plant response (e.g., compromised hydraulic conductivity under water stress, limited soil moisture availability to plants) (Brodribb et al., 2005;Medlyn et al., 2012).Different plants or ecosystems may adjust differently, resulting in divergent responses of mWUE to changing VPD.Understanding the relationship between mWUE and VPD is important given that VPD is expected to keep increasing in the future, which will exert further water stress on plants (Ficklin & Novick, 2017;Grossiord et al., 2020;Novick, Ficklin, et al., 2016;Zhang et al., 2019).Furthermore, while soil moisture is a stochastic variable due to its dependency on intermittent rainfall, VPD is smoother in time and easier to monitor through various meteorological or gas exchange measurement techniques.Although VPD and soil moisture limit plants' carbon uptake and water use independently (Yi et al., 2019), VPD can be used as a proxy of water stress at a sub-daily scale where VPD plays a primary role in regulating stomatal regulation unless severe soil moisture deficiency, as indicated by the models with sub-daily timesteps (e.g., Ball-Berry model and its variations), and in turn influencing the balance between carbon uptake and water loss (i.e., water-use efficiency) at a sub-daily scale (Baldocchi et al., 2022;Grossiord et al., 2020;Novick, Ficklin, et al., 2016).Therefore, examining the association between mWUE and VPD would add insight into the predictability of soil moisture alone.
The objectives of this study are (a) to investigate the variation of mWUE at an hourly timescale in response to changing VPD and (b) to explore approaches for estimating mWUE explicitly from the modeled relationship between intrinsic water-use efficiency (iWUE, carbon assimilation per unit stomatal conductance, representing water-use efficiency at leaf level) and VPD.The Ball-Berry model (Equation 1) reveals that the parameter g 1 , which is proportional to Medlyn et al., 2012), is related to A/g s (= iWUE at leaf level).The iWUE can be more straightforwardly estimated from field measurements across various spatiotemporal scales, including leaf gas exchange (daily to weekly at the leaf level), dendrochronology (seasonal/annual at the tree level), and eddy covariance (hourly at the stand level) (more discussion on iWUE at different scales is available from Beer et al., 2009;Yi et al., 2019).Notably, the inference of iWUE from tree-ring analyses provides an avenue for understanding historical variations in iWUE and, potentially, mWUE.While iWUE has a mathematically simpler form and thus facilitates modeling its response to water stress, the complex mathematical expression of mWUE poses challenges in generalizing its variability at a sub-daily timescale.By elucidating the correlation between iWUE and mWUE, we can gain insights into the response of mWUE to water stress.Additionally, through site comparisons, we further explore whether there is an emerging pattern in the correlation between iWUE and mWUE across different vegetation types and aridity levels.

FLUXNET Data
We obtained half-hourly measurements of carbon and energy fluxes, along with ancillary environmental data, from 115 flux towers across FLUXNET sites.These data were collected using the FLUXNET 2015 Tier 1 database (Table S1 in Supporting Information S1) (Pastorello et al., 2020).Eddy covariance records, which have the benefit of providing continuous meteorological and gas exchange data at the high temporal resolution, are very well suited for investigating the relationship between gas exchange dynamics, mWUE, and VPD at the ecosystem scale.
We selected the study sites from six vegetation types (grassland, cropland, shrubland, savanna, broadleaf forest, and needleleaf forest, based on the International Geosphere-Biosphere Programme (IGBP) land cover classification system; Loveland & Belward, 1997) based on the data availability for the variables required for the analysis.For reliable and clear mWUE analysis, we only included the sites that had at least 3 years of data and a strong iWUE-VPD correlation.Specifically, we selected the sites that had a coefficient of determination (R 2 ) > 0.8 with any of the three model fits-linear, quadratic, or Michaelis-Menten, which was the case for more than 70% of the sites over 3 years of data (refer to Section 2.4 for more information about the model fits).In addition, we only used the data where net ecosystem exchange (NEE), latent heat flux (LE), and sensible heat flux (H) were either original measurements (quality control flag = 0) or gap-filled data of good quality (quality control flag = 1) to ensure data quality and make the most of the data.We only used daytime data when net radiation was greater than 0 W m 2 without precipitation.While we acknowledge the potential benefits of excluding more days after rainfall (e.g., Lin et al., 2015), we believe that omitting only the precipitation days is sufficient for our analysis.This is because iWUE had low variability under humid conditions, as evidenced by the low standard deviations of IWUE under low VPD levels in Figure 2. Additionally, we implemented a procedure to remove outliers in soil water content and relative humidity as described in the following paragraph, which would help mitigate the impact of periods after rainy days on our analysis.
We limited our analysis to the growing season, where daily GPP was larger than 10% of the 95th percentiles of daily GPP for each site with >5°C air temperature.We used the GPP partitioned based on the standard daytime method (variable name: GPP_DT_VUT_REF, Lasslop et al., 2010).Additional filtering criteria were applied for some key variables: atmospheric CO 2 concentration between 350 and 420 ppm, friction velocity (u*) greater than 0.1 m s 1 , and canopy conductance calculated by Penman-Monteith equation (Monteith, 1965) greater than 0.05 mol m 2 s 1 .Lastly, we removed outliers of the environmental drivers and biological variables (i.e., air temperature, relative humidity, atmospheric CO 2 concentration, latent heat flux, wind speed, VPD, atmospheric pressure, friction velocity, net radiation, soil water content, canopy conductance, iWUE, and mWUE) by excluding data that were below the fifth or above the 95th percentiles of each variable.Note that the purpose of data filtering was to remove exceptionally low or high values of the variables, which we consider outliers.Our goal was to ensure that the results, especially the variability of mWUE, were not unduly influenced by these outliers.We carefully examined the histograms for the variables for each site to minimize data reduction while retaining useful information.

Two Different Approaches Describing mWUE
We used two different approaches for describing the mWUE: two optimization-theory-driven mWUE, the solution of "∂E/∂A" suggested by Katul et al. (2010) and the "g 1 " parameter proposed by Medlyn et al. (2012).The difference between the optimization-theory-driven mWUE is based on their interpretation of stomatal optimization.Katul et al. (2010) assumed that stomata are optimizing for photosynthesis limited by Rubisco activity (i.e., carbon-limited), and plant stomatal optimality is subject to change (i.e., mWUE is not constant).On the other hand, Medlyn et al. (2012) assumed that stomata are optimized for photosynthesis limited by Ribulose-1,5bisphosphate (RuBP) regeneration (i.e., light-limited).In either case, the optimization objective should result in constant mWUE values at short timescales- Katul et al. (2010) suggested approximately 10 min, whereas Medlyn et al. (2012) suggested daily or longer-although it may change at longer timescales as hydrologic conditions evolve.
Following Katul et al. (2010), the ∂E/ ∂A emerges from an optimality condition determined with a linearized variant of the Farquhar, Schulze, and Kuppers (1980) photosynthesis model, defined as: where iWUE is defined as a ratio of A to g s at the leaf-scale (Beer et al., 2009).
The other perspective on optimality proposed by Medlyn et al. (2012) takes an analogous form to an empirical model proposed by Leuning (1995) (Equation 2): This approach indicates that the parameter g 1 represents a slope between g s and A/ c a ̅̅̅̅̅̅̅̅̅̅ VPD √ and is proportional to Lin et al., 2015;Medlyn et al., 2012).Therefore, to facilitate comparison between the two approaches, we compare ∂E/∂A with squared g 1 (i.e., g 1 2 ) throughout the results.Equation 4 was rearranged with an assumption that g 0 , which represents cuticular conductance to water vapor, is negligible (but refer to Manzoni et al. (2011) and Lanning et al. (2020) for discussion of the role of cuticle conductance under drier conditions): Consequently, two different mWUE parameters, ∂E/∂A (mol H 2 O • kPa • mol 1 of dry air) and g 1 (mol H 2 O • kPa 0.5 • mol 1 of dry air), were expressed as functions of iWUE, c a , and VPD.Assuming c a is relatively stable over a short period, we focus on how iWUE (as a biological factor) and VPD (as an indicator of water stress governing plant response at a short temporal scale, e.g., sub-daily) affect both mWUE parameters.We applied an approximation of iWUE at the ecosystem level, inherent WUE (IWUE), defined by Beer et al. (2009).IWUE (μ mol C mol 1 H 2 O) was particularly suitable for our study because IWUE can be calculated from the measurements of carbon and water fluxes by eddy covariance technique and ancillary meteorological data, that is, gross primary productivity (GPP; μ mol m 2 s 1 ) from NEE representing canopy-level carbon assimilation, evapotranspiration rate (ET, mol m 2 s 1 ) from latent heat flux, VPD under the assumption of equal temperatures of leaves and atmosphere, and atmospheric pressure (P a , kPa): Several important assumptions for the definition of IWUE include (a) small and invariant soil evaporation (E) compared to plant transpiration (T ) over the course of the day (hence ∆ET ∼ ∆T ) especially during days without rainfall (conditions we used for our analysis), (b) thermal equilibrium between leaf and air, which influences VPD, and (c) disregarding of aerodynamic resistance through the boundary layer that can change depending on the vegetation structure (refer to Beer et al. (2009) for more discussion about IWUE as a proxy of ecosystem-level iWUE).We confirmed the robustness of IWUE as a proxy of iWUE at the ecosystem level by comparing it with a few other definitions of iWUE (the comparison results are available in the Supporting Information S1; Figures S1  and S2).Note that IWUE and mWUE were computed using half-hourly FLUXNET data; hence, their variabilities discussed here represent plant physiological response at a sub-hourly scale.

Sensitivity of mWUE Parameters to Moisture Conditions
Variations of mWUE parameters in response to moisture conditions (i.e., atmospheric water demand and sitelevel aridity) were evaluated at the individual site level and across sites.For the individual sites, mWUE parameters were partitioned into discrete bins spanning a range of VPD.To avoid biases from unevenly distributed data points across the range of VPD (i.e., sample sizes at low and high VPD are smaller than those for the intermediate level of VPD), data binning was performed in a way that the sample sizes were evenly distributed into 30 bins across the range of VPD at each site.Then, mWUE-VPD relationships were produced based on the mean mWUE values generated for the different VPD bins.
To compare across the sites, the relationships between site-specific mWUE and aridity index (AI) were evaluated (refer to Figure S3 in the Supporting Information S1 for AI at all the study sites).AI was defined as the ratio of annual precipitation (P) to annual potential evapotranspiration (PET) and averaged over the observation period for each site (UNEP, 1992): The annual PET was determined by summing up the half-hourly PET values over the course of a year, using the United Nations Food and Agriculture Organization (FAO) Penman-Monteith method as outlined by Allen et al. (1998): where Δ is the slope of vapor pressure curve (kPa °C 1 ), R n is the net radiation (MJ m 2 hr 1 ), G is the soil heat flux density (MJ m 2 hr 1 ), γ is the psychrometric constant (kPa °C 1 ), T a is the air temperature (°C), u is the wind speed (m s 1 ), e s is the saturation vapor pressure (kPa), and e a is the actual vapor pressure (kPa).The estimation of AI is sensitive to gaps in precipitation data.Therefore, we used long-term mean annual precipitation provided on the site information page at the FLUXNET website (https://fluxnet.org/sites/site-list-and-pages/)rather than calculating mean annual precipitation from the FLUXNET2015 data set.For the sites where annual precipitation records were not provided, the high-frequency precipitation record in the FLUXNET2015 data set was used.

Assessing the Relationship Between mWUE and IWUE
As a first step to conceptually understand the relationship between mWUE and IWUE, the relationship between IWUE and VPD was modeled by three hypothetical functions-linear, quadratic, and the Michaelis-Menten functions-based on the observations across the sites.
where m is the slope of IWUE L , n is IWUE L at VPD = 0, IWUE max is the maximum potential IWUE, k is the VPD at which IWUE proceeds at half IWUE max , a represents the curvature of IWUE Q , b is the vertex, c is the maximum IWUE Q at the vertex.
The expected dynamics of mWUE across the FLUXNET sites in response to changing VPD were simulated based on an empirically driven IWUE-VPD model to understand how the mWUE metrics would respond to changing VPD and IWUE.To generate possible patterns of mWUE-VPD, the range of coefficients in the IWUE models was determined empirically from the data across the sites.To facilitate interpretation, the patterns were simulated by changing the curvature of the quadratic equation (Equation 11), assuming the intercept is equal to zero.For the simulation of mWUE, a constant c a was applied by calculating its average across the sites to focus on the interactions among VPD, IWUE, and mWUE (Equations 3 and 5).
Lastly, we investigated how IWUE (as a biological factor) and AI (as an environmental driver) influence the variability of mWUE.Based on Equations 3 and 5, we hypothesized that a simple relationship between mWUE and the inverse of IWUE (IWUE 1 ) might emerge and would be affected by changing moisture conditions.Therefore, we identified a relationship between mWUE and IWUE 1 for each study site and examined whether the relationship can be generalized across the sites based on the site-specific AI.

Empirical Response of IWUE to Changing VPD or AI
To test the robustness of IWUE as a proxy of iWUE at the ecosystem level, we first compared the two different definitions of intrinsic water-use efficiencies at stand level, GPP divided by surface conductance (G s ) (i.e., iWUE = GPP/G s ) and inherent WUE (i.e., IWUE = GPP/ET × VPD/P a ).The two WUE definitions were linearly correlated across the study sites (Figure 1), and most sites had R 2 values larger than 0.95 (Figure 1b), indicating the robustness of IWUE as a proxy of iWUE at the ecosystem level (refer to the Supporting Information S1 for an additional comparison of multiple definitions of iWUE; Figures S1 and S2 in Supporting Information S1).We also performed the entire analysis using these two WUE definitions and observed similar results, which led to the same conclusion.Therefore, we only show the results from using IWUE hereafter.
In most cases, the Michaelis-Menten model and the quadratic model explained empirical IWUE patterns across the range of VPD better than the linear model (Figure 2 and Figure S3 in the Supporting Information S1).Specifically, the Michaelis-Menten model worked better for the sites where the increase of IWUE plateaued at high VPD, and the quadratic model worked better for the sites where IWUE started decreasing at very high VPD.
On the other hand, the linear model often overestimated IWUE at low and high VPD, except the sites where IWUE-VPD was highly linear.
When the site-specific IWUE-VPD slope values derived from the linear model (i.e., m in Equation 9) were combined, we found increasing m with rising AI (P < 0.001, Figure 3a).However, site-level aridity did not influence the intercept of IWUE-VPD relationship (P > 0.05, not shown here).When the sites were divided by their vegetation types, m increased with a rising AI in all vegetation types.However, the trend was only significant in grasslands, croplands, and shrublands (P < 0.05, Figure 3).VPD (Figure 4).When the iWUE-VPD is strongly linear, mWUE decreased exponentially and became less variable as VPD increased (Brighter curves in Figures 4b and 4c).However, as the iWUE-VPD relationship became more non-linear, mWUE declined at lower VPD and then increased at higher VPD (i.e., concave-up), rendering the mWUE-VPD relationship non-monotonic (darker curves in Figures 4b and 4c).

Response of mWUE to
The simulated patterns of mWUE-VPD agreed well with the patterns from the empirical observation when the appropriate function for the IWUE-VPD relationship was applied.We show mWUE-VPD relationships from three study sites as examples (Figure 5), of which IWUE-VPD was represented best by linear, the Michaelis-Menten, and quadratic functions, respectively (refer to Figure 2 for their corresponding IWUE-VPD relationships.Also, refer to Figure S5 in the Supporting Information S1 for the results of all study sites).As indicated by the simulation, the site with highly linear IWUE-VPD (IT-BCi) showed exponentially decreasing mWUE with rising VPD.In contrast, the other two sites with highly non-linear IWUE-VPD relationships had a concave-up pattern of mWUE-VPD.Notably, the mWUE-VPD relationship generated using a less optimal IWUE-VPD model can differ substantially from the empirical pattern.For example, application of linear IWUE-VPD function to the CA-NS2 and US-Ton, the sites represented best by the Michaelis-Menten and quadratic functions, respectively, generated concave-down mWUE-VPD pattern that is opposite to the actual pattern (Figure 5).The disagreements between models and observations increased as VPD approached very high and very low extremes.
The variability of mWUE to changing VPD was substantial in most cases (Figure 6).Out of the total of 115 study sites, the percent increase of ∂E/∂A (i.e., growth in ∂E/∂A from the lowest to the largest value at a site) was larger than 50% in 43 sites, and larger than 100% in 22 sites.Note that the reported percent increase was determined by excluding the upper and lower 10% of values.This step was taken to prevent exaggeration caused by extremely high ∂E/∂A at low VPD, which is commonly observed across the study sites (refer to Figure S5 in the Supporting Information S1 for the variability of ∂E/∂A with VPD at all the study sites).As a result, the reported percent increase represents a conservative estimate overall.

Correlation Between mWUE and IWUE
Although the trend of mWUE-VPD seems hard to generalize, the simulated mWUE had a clear linear relationship with IWUE 1 for the majority of IWUE's range regardless of the linearity of the IWUE-VPD relationship except when IWUE is very high (i.e., under high VPD, Figure 7).Although limited to a small portion of the entire range, a sharp directional change in the variation of mWUE was near a point where IWUE 1 was smallest, and strong linearities between and IWUE 1 were found before and after the transitional point.Substantial hysteresis became more evident as the IWUE-VPD pattern became more curved (darker curves in Figure 4).
As predicted by the simulated mWUE-IWUE 1 relationships (Figure 7), the empirical mWUE-IWUE 1 relationship was strongly linear (P < 0.001 at all sites, Figure 8).A sign of hysteresis was noticeable for the site that showed decreasing iWUE under very high VPD (US-Ton, refer to Figure 2 for its IWUE-VPD relationship).In contrast, hysteresis was less evident at the other sites.When the relationship was drawn by grouping data by different levels of IWUE (black dots in Figure 8), hysteresis was not observed, and the mWUE-IWUE 1 relationship was strongly linear.
We investigated whether the relationship between mWUE and IWUE 1 could be generalized across the sites based on the site-specific AI.Specifically, the linear IWUE 1 -mWUE slopes (hereafter m*) from all study sites were merged, and their variability in response to changing AI was evaluated.We found a significant linear relationship between m* and AI when both are scaled by log 10 (P < 0.001, Figure 9).The m* was larger at the drier sites (i.e., sites of lower AI) than at the wetter sites (i.e., sites of larger AI).However, we did not find a significant relationship between the IWUE 1 -mWUE intercept and AI (P > 0.05, not shown here).
We further tested whether we could find the similar relationship when the sites were grouped by the vegetation type.We found decreasing m* with rising AI in grasslands, croplands, and shrublands (P < 0.01, Figure 10).On the other hand, m* was relatively constant across the range of AI in savannas, deciduous broadleaf forests, and evergreen needleleaf forests (P > 0.05, Figure 10).

Discussion
Stomatal optimization theory, which originated with the work of Cowan and Farquhar (1977), has experienced a recent surge in popularity as the modeling community continually seeks ways to inject more theoretical rigor into Earth system models (Anderegg et al., 2018;Bassiouni & Vico, 2021;Bonan et al., 2014;Feng et al., 2022;Katul et al., 2009Katul et al., , 2010;;Lin et al., 2015Lin et al., , 2018;;Lu et al., 2016Lu et al., , 2020;;Medlyn et al., 2012Medlyn et al., , 2017;;Novick, Miniat, & Vose, 2016;Sabot et al., 2022;Sperry et al., 2017;Wolf et al., 2016).The mWUE is a key parameter in this type of model, but we still need a clear understanding of how mWUE is regulated biologically and environmentally.Lin et al. (2018) previously suggested suboptimal mWUE in response to VPD at a sub-daily scale by estimating site-specific, best-fitted exponent for VPD based on the variation model of optimality theory (Medlyn model), which inspired our study.In comparison, our study is unique in analyzing the dynamics of mWUE observed at the half-hourly timescale in response to changing VPD owing to the long-term continuous carbon and water flux data from the network of eddy covariance towers.
Another motivation for our study was the conflicting arguments over the dynamics of mWUE in response to water stress.Although mWUE is in general considered to be nearly constant during a day under stable soil moisture conditions (Berninger & Hari, 1993;Fites & Teskey, 1988;Hall & Schulze, 1980;Hari et al., 2000), several studies showed that mWUE changed with different levels of water stress.For example, theoretical considerations indicate a monotonic decrease of mWUE with higher water stress when the stochasticity of rainfall depths is neglected (Cowan, 1982;Makela et al., 1996), although some empirical estimates showed that mWUE increases under severe water stress (Farquhar, Schulze, & Kuppers, 1980;Grieu et al., 1988).On the other hand, Manzoni et al. ( 2011) performed a meta-analysis of 50 species to estimate mWUE from gas exchange observations along gradients of soil moisture and showed that mWUE decreases under mild water stress but increases under severe water stress (note that they defined λ ∂A/∂E, which is inverse of the definition used by Cowan and Farquhar (1977) and this study).

Relationship Between IWUE and VPD
Based on the two equations of stomatal optimization theory (Equations 3 and 5), we first characterized the variability of mWUE based on the relationship between IWUE and VPD, representing biological and environmental factors, respectively.We show that the variability of IWUE needs to be modeled accurately to emulate the variability of mWUE in response to water stress correctly.For example, as demonstrated in Figure 5 (CA-NS2 and US-Ton), oversimplifying the IWUE-VPD relationship by modeling it with a linear function can incorrectly interpret mWUE variability.
The non-linear IWUE-VPD relationship is presumably driven by different rates of carbon assimilation and water loss in response to changing VPD at an hourly scale, reflecting the balance between stomatal and non-stomatal limitations to photosynthesis at the leaf level (Farquhar, 1978;Jones, 2014).Under low to moderate VPD conditions, photosynthesis is less sensitive to changing intercellular CO 2 concentration because stomatal conductance is high enough to maintain the high intercellular CO 2 when VPD is low to moderate.Therefore, the quantity Figure 6.Sorted percent increase of ∂E/∂A (from the lowest ∂E/∂A) (GRA: grassland, CRO: cropland, SH: shrubland, SAV: savanna, BF: broadleaf forest, NF: needleleaf forest).Embedded plots in GRA and SH are zoomed in for those sites where percent increases are lower than 100%.Note that the percent increases were calculated after removing values of the highest 10% and lowest 10% to avoid exaggeration due to very high ∂E/∂A at low vapor pressure deficit at some sites.Therefore, the reported percent increase values are conservative estimates for most sites.
of reduced water loss by stomatal closure (ET at an ecosystem level) outweighs the quantity of reduced carbon assimilation (GPP at an ecosystem level) when VPD rises (i.e., |∆GPP| < |∆ET|), resulting in the increasing phase of As VPD keeps increasing, photosynthesis can be limited when the reduction of stomatal conductance under high VPD conditions substantially reduces intercellular CO 2 concentration (i.e., |∆GPP| ≈ |∆ET|), resulting  in the steady phase of IWUE.As VPD becomes excessively high, photosynthesis will be further suppressed by high temperature (Yamori et al., 2014) and low leaf water potential (Lawlor & Tezara, 2009) that are associated with dry conditions (i.e., |∆A| > |∆g |), leading to the decreasing phase of IWUE.
Therefore, assuming a linear IWUE-VPD relationship may not only fail to emulate observations but also oversimplify the physiological responses to water stress.Our analysis recommends using the Michaelis-Menten function for most sites while utilizing a quadratic function for sites exhibiting extreme cases where IWUE declines under high VPD conditions.The Michaelis-Menten function can be beneficial to characterize the IWUE-VPD relationship because the maximum potential IWUE and the rate of IWUE increase can be identified during parameterization (Equation 10).Although the quadratic function can emulate IWUE-VPD relationships very well or performs even better than the Michaelis-Menten function in some cases where IWUE decreases at high VPD, it is parameterized empirically and as a result, lacks mechanistic information.Nevertheless, the quadratic function is preferable to the linear function.
It is also important to consider the definition of water-use efficiency for accuracy.We used IWUE as a proxy of iWUE at the ecosystem level as suggested by Beer et al. (2009), which can be estimated by GPP and ET inferred  from the flux tower observations.This approximation is particularly useful over a more common ecosystem-level iWUE = GPP/G s because IWUE requires fewer variables and is easier to extrapolate to a larger scale.However, it is important to note that ET/VPD in the equation of IWUE (Equation 6) is a proxy of canopy conductance, assuming the canopy is well coupled to the atmosphere, boundary layer resistance is small, and thermal equilibrium between leaf and air is achieved (Beer et al., 2009).In other words, IWUE under non-equilibrium between canopies and atmosphere can be overestimated due to the higher VPD than the vapor pressure gradient near the canopy surface (i.e., the difference between intercellular vapor pressure (e i ) and atmospheric vapor pressure (e a ), e i -e a ).Difference between leaf and air can also influence the e i -e a ; if leaf temperature is higher than air temperature (as it often is, e.g., Novick & Barnes, 2023;Yi et al., 2020), e i will increase while e a remains constant, resulting in larger e i -e a than measured VPD and consequently underestimate IWUE.Therefore, it is important to address this potential bias to quantify iWUE accurately.According to our results, the correlation between the two ecosystem-level iWUE proxies was strong at the site level (Figure 1), implying that the choice of ecosystem-level iWUE definition is unlikely to influence our interpretation of the iWUE and mWUE variabilities substantially.Furthermore, our comparison of multiple definitions of iWUE using a mechanistic model, CAN-VEG (refer to the Supporting Information S1 for more details), indicated that IWUE is a good proxy of leaf-level iWUE and meets the general assumptions to address scaling issues.Thus, we conclude that eddy covariance observation of carbon and water fluxes is suitable to model the variability of iWUE in response to changing VPD.
Of note, the linear relationship between the slope of IWUE-VPD and AI (Figure 4) was stronger in the ecosystems characterized by lower vegetation types (e.g., grasslands, croplands, and shrubland).In contrast, ecosystems with higher vegetation (e.g., savannas, broadleaf forests, and needleleaf forests) exhibited a weaker relationship.This observation implies a potential link between water-use efficiency and the vertical structure of vegetation, although the exact underlying mechanism remains uncertain.

Modeling the Variability of mWUE
We compared two solutions of mWUE by Katul et al. (2010) (∂E/∂A) and Medlyn et al. (2012) (g 1 ) developed based on different assumptions on stomatal optimality (carbon-limited vs. light-limited) for more robust conclusion.Despite the difference in the assumption, both solutions yielded very similar results throughout our analysis, confirming that the optimality assumption had little effect on evaluating the variability of mWUE in response to changing moisture conditions.
We characterized the trend of mWUE by using VPD as an environmental driver (Figures 4 and 5), where its variability in response to VPD was unique and not necessarily unidirectional, thus making it hard to generalize with commonly available functions.Specifically, the variability of mWUE was simpler and decreased exponentially with rising VPD when the IWUE-VPD relationship was more linear, making it easy to model the mWUE-VPD relationship (Figures 4 and 5).However, the variability of mWUE was not unidirectional when the IWUE-VPD relationship was non-linear, as observed in most cases (Figure S5 in the Supporting Information S1); high variability in mWUE is usually observed at low-and high-ends of VPD.On the other hand, when mWUE was calculated under conditions of moderate VPD level only, the variability of mWUE can be overlooked and considered constant.This complex pattern signifies the importance of a comprehensive view of IWUE and mWUE across the full potential range of VPD.Observation under conditions of a partial range of environmental factors is common in many types of field measurements that have coarser time resolution (hourly vs. daily to weekly, e.g., eddy covariance vs. leaf gas exchange measurements) unless they are performed frequently over a long period to cover non-typical conditions.We were able to estimate precise variability of mWUE matching with the hypothetical models owing to the large amount of data (FLUXNET2015) collected every half-hour over the long period throughout the network of flux towers (total 1,036 site years with many sites offering data collected over more than a decade), highlighting the value of long-term, continuous measurements.Overall, our result of the mWUE-VPD relationship supports the results of Manzoni et al. (2011) among the various conflicting results over the response of mWUE in response to water stress, which found decreasing mWUE under mild water stress and increasing mWUE under severe water stress from a meta-analysis of gas exchange observations.As a solution to model unique patterns of mWUE, we attempt to address its variability with information that can be obtained easily from various types of field measurements (e.g., eddy covariance, gas exchange, and tree-ring cores) and modeled empirically-IWUE.The relationship between mWUE and IWUE was inferred from the two equations of the optimization theory (Equations 3 and 5).We found a strong linear correlation between IWUE 1 and mWUE from both empirical data (Figure 8) and modeling exercise (Figure 7).In other words, the variability of mWUE in response to changing VPD can be characterized by (a) the function of IWUE-VPD relationship and (b) the slope between IWUE 1 and mWUE.The relationship between IWUE-VPD is relatively simple and can be identified with various field measurements.This raises the question of whether a simple way exists to identify the slope between IWUE 1 and mWUE.By synthesizing the IWUE 1 -mWUE slopes across the sites, we found that the IWUE-mWUE slope is highly correlated with the site-specific AI that can be characterized for different vegetation types (Figure 9).The correlation is conceivable from the equations of mWUE (Equations 3 and 5).If, for instance, Equation 3 is rearranged, ∂E/∂A IWUE 2 ∝ VPD (12) indicating that the slope between mWUE and the inverse of IWUE is proportional to VPD, which is commensurate with AI at a site-level.The correlation between the IWUE 1 -mWUE slope and the AI at a site level implies that the AI is a good surrogate for the site-specific IWUE 1 -mWUE slope.
We further illustrated how the correlations between the IWUE 1 -mWUE slope (m*) and AI vary across vegetation types (Figure 10).Among the vegetation types, GRA, CRO, and SH had strong correlations between m* and AI, which indicated that using different m* depending on the site-level dryness would be appropriate.On the other hand, the low variability of m* observed in SAV, BF, and NF indicates that constant m* can generate a reasonably accurate mWUE-VPD relationship regardless of the site-level dryness.Although the reasons for this difference are not entirely clear, this empirical relationship can help more accurately model the variability of mWUE in response to changing VPD across the sites and biomes.Growth in data availability across the flux tower network helps ensure the coverage of the full potential range of environmental factors.More data availability can be achieved by consistently collecting good-quality data from existing study sites and establishing new sites in underrepresented areas.Furthermore, additional data would also help the development of m* in detail, for instance, based on the plant water-use strategies, with the aid of conjoined field measurements such as water potential (Ψ ) of soil and plant.

Implications for Future Research
It is important to note that plant response to water stress is not only determined by the water demand (i.e., atmospheric dryness or VPD) but also by the availability of water sources (i.e., soil moisture).Although volumetric soil moisture content (θ) is often considered as a metric of soil water available to plants, soil water potential (Ψ S ) is the driving force of water flows that becomes available to plants by moving along gradients of water potential through the plant (stem and leaf) and eventually to the air.Moreover, Ψ S is not only determined by the θ but also by soil physical properties, and thus can differ even under conditions of the same θ (Campbell, 1974;van Genuchten, 1980).Unlike Ψ S , θ is widely measured and usually available with flux data, and carbon and water fluxes are often explained as a function of θ (Green et al., 2019;Novick, Ficklin, et al., 2016).However, θ may not characterize soil moisture availability to plants properly, and its relationship with carbon and water fluxes is usually nonlinear and threshold-driven (Feldman et al., 2019;Novick et al., 2022;Stocker et al., 2018), making the modeling of the relationship between IWUE and soil moisture availability challenging.Therefore, enhanced accessibility to Ψ S data by improving the ease and reliability of Ψ S observations, for example, by building a centralized and standardized network of Ψ (Novick et al., 2022) would be a necessary step to better characterize the effect of soil moisture availability on plant responses such as IWUE and mWUE.
In this study, we tested the two stomatal optimization models (Katul et al., 2010;Medlyn et al., 2012) that are elaborations of the original Cowan and Farquhar (1977) model with few modifications because our goal was to characterize variability of mWUE in response to dryness (VPD and AI) using IWUE that can be calculated from the extensive, long-term continuous data from the network of eddy covariance.Meanwhile, more recent optimization models are incorporating additional physiological penalties than the water loss, for instance, damage to the vascular system induced by water stress (Anderegg et al., 2018;Sperry et al., 2017;Wolf et al., 2016), which may enhance prediction of long-term plant responses to climate change.Although monitoring the integrity of the vascular system, which can be informed by the dynamics of hydraulic conductivity, has not been widely conducted, recent advances in psychrometric approaches allowing continuous measurements of plant Ψ (e.g., PSY1 manufactured by ICT International) and Ψ S (e.g., TEROS 21 manufactured by Meter Group) are now enabling the monitoring the dynamics of hydraulic conductivity.Moreover, the relationship between plant and soil Ψ can be used to identify plant water-use strategies (e.g., isohydry framework; Martinez-Vilalta et al., 2014), which can help develop m* based on plant water-use strategies.The measurements of carbon and water fluxes using the eddy covariance technique with the aid of the centralized and standardized deployment of Ψ sensors (Novick et al., 2022) will have a great potential to test models and theories of stomatal optimization and advance our knowledge of it.
Changing VPDBoth of the mWUE indices, ∂E/∂A and squared g 1 (g 1 2 ), showed a very similar response to changing VPD and indicated that the directional change of mWUE can be interpreted differently depending on the pattern of IWUE-

Figure 1 .
Figure 1.Comparison of two different definitions of water-use efficiencies at all sites (a) and at three sample sites (c, d, e): inherent water-use efficiency (IWUE) at the ecosystem level (=GPP/ET × VPD/P a ), and intrinsic water-use efficiency (iWUE) at the ecosystem level (=GPP/G s ).Refer to Beer et al. (2009) for a comparison of different definitions of water-use efficiencies at the leaf and ecosystem levels.Individual dots in Panels (a, c, d, and e) indicate WUE partitioned into discrete bins spanning a range of VPD.Solid red lines indicate significant linear regressions (P < 0.05), and dashed red lines indicate 95% confidence interval.Dashed gray lines represent 1:1 lines.Panel b shows the histogram of coefficients of determination (R 2 ) of the linear fits between IWUE and iWUE across the study sites.

Figure 2 .
Figure 2. Examples of empirical (black dots) and modeled (linear: blue, Michaelis-Menten: green, quadratic: red) responses of IWUE to changing VPD.The examples include three sites best represented by the linear model (IT-BCi, cropland), the Michaelis-Menten function (CA-NS2, needleleaf forest), and the quadratic model (US-Ton, savanna), respectively.Each error bar (light gray) represents the standard deviation of IWUE for each VPD bin (95% confidence).Refer to Figure S4 in the Supporting Information S1 for the IWUE-VPD relationships of all the study sites (n = 115).

Figure 3 .
Figure 3. Relationship between the site-level aridity index and the regression slope of IWUE-VPD from individual sites (i.e., m in Equation 9).Panel a shows the relationship when all sites were consolidated.The relationship is also illustrated separately for six different vegetation types in Panels (b)-(g) (GRA: grassland, CRO: cropland, SH: shrubland, SAV: savanna, BF: broadleaf forest, NF: needleleaf forest).Each circle represents m from an individual site.Error bars represent standard errors of linear regressions.Solid lines indicate significant linear relationships (P < 0.05) and dashed lines indicate 95% confidence intervals.

Figure 4 .
Figure 4. Hypothetical models of IWUE-VPD relationship (a), simulated ∂E/∂A-VPD (b) and g 1 2 -VPD (c) relationships based on typical cases, and their corresponding patterns illustrated using observations from all study sites (d, e, and f).The marginal water-use efficiency curves are the results of using the IWUE-VPD relationships of the corresponding colors.Note that IWUE-VPD relationships become more linear with lighter colors.

Figure 5 .
Figure 5. Examples of empirical (black dots) and modeled (linear: blue, Michaelis-Menten: green, quadratic: red) relationships between ∂E/∂A (analytical solution by Katul et al., 2010) and VPD, and between g 1 2 (Medlyn et al., 2012) and VPD.The examples include three sites best represented by the linear IWUE-VPD model (IT-BCi, cropland), the Michaelis-Menten function (CA-NS2, needleleaf forest), and the quadratic model (US-Ton, savanna), respectively.Note that the terms "linear", "Michaelis-Menten", and "quadratic" denote the regression fits for the IWUE-VPD relationships (Refer to Figure 2 for the IWUE-VPD relationships at the corresponding sites).Each error bar (light gray) represents the standard error of the mean inherent water-use efficiency for each VPD bin (95% confidence).Refer to Figure S5 in the Supporting Information S1 for the ∂E/∂A -VPD relationships at the 115 study sites.

Figure 7 .
Figure 7. Simulated relationship between marginal water-use efficiency metrics (∂E/∂A (a) and g 1 2 (b)) and IWUE 1 (based on the hypothetical IWUE-VPD model in Figure 4).The colors of the curves correspond to those used in Figure 4: IWUE-VPD relationships become more linear with lighter colors.Dashed arrows in panel a represent the directional change of VPD from low to high VPD.

Figure 8 .
Figure 8. Examples of empirical relationship between marginal water-use efficiency metrics (∂E/∂A and g 1 2 ) and IWUE 1 .The examples include three sites best represented by the linear IWUE-VPD model (IT-BCi, cropland), the Michaelis-Menten function (CA-NS2, needleleaf forest), and the quadratic model (US-Ton, savanna), respectively.Refer to Figure 2 for the IWUE-VPD relationships at the corresponding sites.Colorful dots represent hourly data points shaded based on the level of VPD (refer to color bars for the scale of VPD).Black dots represent data binned by IWUE 1 : Data binning was performed to distribute sample sizes evenly across bins (∼30 samples per bin).Error bars represent standard deviations.The red and black solid lines indicate linear fits for hourly and binned data, respectively.Dashed red lines represent confidence intervals for the slopes of linear regressions.Note that red and black linear regressions and their confidence intervals overlap.Refer to Figure S6 in the Supporting Information S1 for the ∂E/∂A-IWUE 1 relationships at the 115 study sites.

Figure 9 .
Figure 9. Relationships between IWUE 1 -mWUE slope (∂E/∂A (a) and g 1 2 (b)) and aridity index (=P/PET) derived from all the study sites (n = 115).Each circle represents the slope obtained from an individual site.Both the x and y axes are scaled by log 10 .The numbers in parentheses next to the x-axis tick labels represent the aridity indices before the log transformation.The solid lines indicate linear regressions, and the dashed lines indicate confidence intervals (95% confidence interval).

Figure 10 .
Figure 10.Relationships between log-transformed IWUE 1 -mWUE slope and aridity index in different vegetation types (GRA: grassland, CRO: cropland, SH: shrubland, SAV: savanna, BF: broadleaf forest, NF: needleleaf forest).Each circle represents the log-transformed slope obtained from an individual site.The numbers in parentheses next to the x-axis tick labels represent the aridity indices before the log transformation.Solid lines indicate significant linear relationships (P < 0.05), and dashed lines indicate 95% confidence intervals.

Table 1 A
Glossary of Terms Related to Water-Use Efficiency iWUE Intrinsic water-use efficiency; leaf-level water-use efficiency (=A/g s ) IWUE Inherent water-use efficiency; a proxy of intrinsic water-use efficiency at the ecosystem level (=GPP/ET × VPD/P a , Beer et al., 2009) m* The slope of the linear relationship between IWUE 1 and mWUE mWUE Marginal water-use efficiency, the ratio of a change in E to a change in A (=∂E/∂A) The quadratic model of IWUE-VPD (hereafter IWUE Q ) depicts the case where IWUE increases with VPD until it reaches a maximum and then decreases afterward.describing the actual response of IWUE to VPD, while the Michaelis-Menten function is a more intermediate case.Mathematically, the IWUE L , IWUE M , and IWUE Q take the forms: In other words, when VPD is low, increasing IWUE with increasing VPD reflects a faster decrease of g s than A (due to the high intercellular CO 2 concentration, c i ), whereas decreasing IWUE with increasing VPD at high VPD reflects a faster decrease of A than g s (low g s at high VPD reduces c i and eventually causes the steep decline of A).The linear model (hereafter IWUE L ), on the other hand, represents a simplified IWUE-VPD relationship where IWUE would keep increasing with rising VPD assuming IWUE is only limited by g s but not by photosynthetic capacity.The Michaelis-Menten function (hereafter IWUE M ) represents the saturating IWUE under high VPD but does not account for IWUE reduction.Thus, the linear and quadratic functions are considered plausible "end members"