# One Stomatal Model to Rule Them All? Toward Improved Representation of Carbon and Water Exchange in Global Models

Manon Sabot, Martin de Kauwe, Andy Pitman, Belinda Medlyn, David Ellsworth, Nicolas Martin-StPaul, Jin Wu, Brendan Choat, Jean-marc Limousin, Patrick Mitchell, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Manon Sabot, Martin de Kauwe, Andy Pitman, Belinda Medlyn, David Ellsworth, et al. One Stomatal Model to Rule Them All? Toward Improved Representation of Carbon and Water Exchange in Global Models. Journal of Advances in Modeling Earth Systems, 2022, 14 (4), 10.1029/2021MS002761. hal-03842122

HAL Id: hal-03842122

https://hal.inrae.fr/hal-03842122

Submitted on 13 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Sabot, M. E. B., De Kauwe, M. G., Pitman, A. J., Medlyn, B. E., Ellsworth, D. S., Martin-StPaul, N., Wu, J., Choat, B., Limousin, JM., Mitchell, P. J., Rogers, A., & Serbin, S. P. (2022). One Stomatal Model to Rule Them All? Toward Improved Representation of Carbon and Water Exchange in Global Models. *Journal of Advances in Modeling Earth Systems, 14*(4), [e2021MS002761].

https://doi.org/10.1029/2021MS002761

Publisher's PDF, also known as Version of record

License (if available):

CC BY

Link to published version (if available):

10.1029/2021MS002761

Link to publication record in Explore Bristol Research

PDF-document

This is the final published version of the article (version of record). It first appeared online via Wiley at https://doi.org/10.1029/2021MS002761. Please refer to any applicable terms of use of the publisher.

# University of Bristol - Explore Bristol Research

## General rights

This document is made available in accordance with publisher policies. Please cite only the published version using the reference above. Full terms of use are available: http://www.bristol.ac.uk/red/research-policy/pure/user-guides/ebr-terms/
RESEARCH ARTICLE

10.1029/2021MS002761

**Key Points:**

* Parameter identifiability differs among stomatal conductance schemes, implying some are more suitable to global modeling than others

* In some schemes, seemingly good performance can result from misrepresentation of physiological processes and sensitivities to model drivers

* We identify a subset of hydraulic-based stomatal optimization approaches that could improve predictive capacity in novel climate spaces

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

M. E. B. Sabot,
m.e.b.sabot@gmail.com

**Citation:**

Sabot, M. E. B., De Kauwe, M. G., Pitman, A. J., Medlyn, B. E., Ellsworth, D. S., Martin-StPaul, N. K., et al. (2022). One stomatal model to rule them all? Toward improved representation of carbon and water exchange in global models. *Journal of Advances in Modeling Earth Systems, 14*, e2021MS002761.
https://doi.org/10.1029/2021MS002761

Received 14 AUG 2021

Accepted 17 JAN 2022

© 2022 The Authors. Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_3.png" /></figure>

# One Stomatal Model to Rule Them All? Toward Improved Representation of Carbon and Water Exchange in Global Models

Manon E. B. Sabot¹,² ID, Martin G. De Kauwe¹,²,³ ID, Andy J. Pitman¹,² ID, Belinda E. Medlyn⁴ ID, David S. Ellsworth⁴ ID, Nicolas K. Martin-StPaul⁵ ID, Jin Wu⁶,⁷ ID, Brendan Choat⁴ ID, Jean-Marc Limousin⁸ ID, Patrick J. Mitchell⁹ ID, Alistair Rogers¹⁰ ID, and Shawn P. Serbin¹⁰ ID

¹ARC Centre of Excellence for Climate Extremes, Sydney, NSW, Australia, ²Climate Change Research Centre, University of New South Wales, Sydney, NSW, Australia, ³School of Biological Sciences, University of Bristol, Bristol, UK, ⁴Hawkesbury Institute for the Environment, Western Sydney University, Penrith, NSW, Australia, ⁵INRAE, URFM, Domaine Saint Paul, Centre de Recherche PACA, Avignon, France, ⁶School of Biological Sciences, The University of Hong Kong, Hong Kong, China, ⁷State Key Laboratory of Agrobiotechnology, The Chinese University of Hong Kong, Hong Kong, China, ⁸CEFE, Université Montpellier, CNRS, EPHE, IRD, Paul Valéry University Montpellier 3, Montpellier, France, ⁹CSIRO Agriculture and Food, Hobart, TAS, Australia, ¹⁰Department of Environmental and Climate Sciences, Brookhaven National Laboratory, Upton, NY, USA

**Abstract** Stomatal conductance schemes that optimize with respect to photosynthetic and hydraulic functions have been proposed to address biases in land-surface model (LSM) simulations during drought. However, systematic evaluations of both optimality-based and alternative empirical formulations for coupling carbon and water fluxes are lacking. Here, we embed 12 empirical and optimization approaches within a LSM framework. We use theoretical model experiments to explore parameter identifiability and understand how model behaviors differ in response to abiotic changes. We also evaluate the models against leaf-level observations of gas-exchange and hydraulic variables, from xeric to wet forest/woody species spanning a mean annual precipitation range of 361–3,286 mm yr⁻¹. We find that models differ in how easily parameterized they are, due to: (a) poorly constrained optimality criteria (i.e., resulting in multiple solutions), (b) low influence parameters, (c) sensitivities to environmental drivers. In both the idealized experiments and compared to observations, sensitivities to variability in environmental drivers do not agree among models. Marked differences arise in sensitivities to soil moisture (soil water potential) and vapor pressure deficit. For example, stomatal closure rates at high vapor pressure deficit range between −45% and +70% of those observed. Although over half the new generation of stomatal schemes perform to a similar standard compared to observations of leaf-gas exchange, two models do so through large biases in simulated leaf water potential (up to 11 MPa). Our results provide guidance for LSM development, by highlighting key areas in need for additional experimentation and theory, and by constraining currently viable stomatal hypotheses.

**Plain Language Summary** Water availability is critical for plants to maintain normal function, so droughts have considerable impact on natural ecosystems. However, predicting the impact of future drought on ecosystems is hard because current global models make systematic errors in their predictions of plant responses when water is scarce. In turn, uncertainty in the modeled terrestrial water and carbon cycles remains high. Here, we evaluate a range of new modeling approaches that have the capacity to mechanistically capture plant responses to water stress. Both in theoretical experiments and comparisons to observations, we find large differences among these new modeling approaches in response to water availability and atmospheric dryness. Importantly, some approaches achieve what seems like “good” performance through compensatory mechanisms that are not supported by observations and/or through incorrect representation of plant processes. Our results provide important guidance for future model development, by highlighting areas in need of continued research, and by constraining the range of approaches presently able to reduce uncertainty in modeled plant responses and suitable for inclusion in global models.

## 1. Introduction

Central to the coupling between terrestrial carbon, water, and energy exchanges is the physiology of stomatal behavior (Hetherington & Woodward, 2003; Jones, 1998; Sellers et al., 1997). During photosynthesis, plants
open their stomates to allow atmospheric carbon dioxide (CO₂) to diffuse to the sites of carboxylation and, in doing so, transpire water vapor, thus cooling the leaf surface and exerting a feedback on the atmosphere. There is widespread evidence that stomatal conductance (g_s) varies with species and climate (Lin et al., 2015; Oren et al., 1999) and that species differ in their rate and timing of stomatal closure to prevent excessive water losses via transpiration (Martínez-Vilalta & Garcia-Forner, 2017; Tardieu & Simonneau, 1998; Vialet-Chabrand et al., 2017).

Despite the importance of g_s in regulating terrestrial carbon and water fluxes, our understanding of stomatal behavior and capacity to translate it into models remain incomplete (Buckley, 2019; Lawson et al., 2010; Woodward, 1987). Stomatal schemes embedded in global models (Ball et al., 1987; Jarvis, 1976; Leuning, 1995) often ignore widely observed species-level variations in water use and/or rely on empirical formulations that predetermine stomatal sensitivities to a variety of environmental drivers, including atmospheric dryness (i.e., vapor pressure deficit) and [CO₂]. Highlighting the shortcomings associated with the use of empirical g_s schemes, climate models attribute anywhere between 22% and 58% of the annual global terrestrial evapotranspiration (ET) to plant transpiration (Stoy et al., 2019; Wei et al., 2017). By contrast, observations suggest that transpiration is responsible for >60% of the annual ET over land (Good et al., 2015; Li et al., 2019; Schlesinger & Jasechko, 2014).

If climate models under-estimate transpiration, they also likely misrepresent: (a) inter-annual variability in the terrestrial carbon cycle—which is primarily driven by water availability of the vegetation (Ahlstrom et al., 2015; Humphrey et al., 2018; Jung et al., 2017); (b) land-surface feedbacks to the boundary layer (Donat et al., 2018; Miralles et al., 2014) and the amplification of some climate extremes (Miralles et al., 2018; Yunusa et al., 2015); and (c) vegetation changes in water use under future atmospheric [CO₂] (De Kauwe et al., 2013; Mankin et al., 2019; Swann et al., 2016) and associated impacts on, for example, surface water storage (Trancoso et al., 2017; Ukkola, Prentice, et al., 2016). Resolving the knowledge gap around plant contributions to the carbon and water cycles and to local energy budgets is critical to resolving both long term (terrestrial carbon sink) and short term (amplification of extreme events) model uncertainties.

To address model shortcomings, development of stomatal theory, aided by field and experimental synthesis studies (Choat et al., 2012; Lin et al., 2015; Martin-StPaul et al., 2017a; Mencuccini et al., 2019; Miner et al., 2017), has recently become a major focus of global change ecophysiology. First, g_s behavior under well-watered conditions has been shown to vary predictably by species (Lin et al., 2015), which can be used to parameterize stomatal conductance schemes on a plant functional type basis (De Kauwe et al., 2015). Second, g_s can be accurately simulated to decline with decreasing soil moisture empirically (Mäkelä, 1996; Sala & Tenhunen, 1996; Tuzet et al., 2003), or using measured vulnerability curves that represent the progressive impairment of water transport as plant hydraulic conductance is lost (Sperry et al., 2017; Wolf et al., 2016). Third, optimality principles positing that stomatal behavior should balance carbon uptake and water loss (Cowan & Farquhar, 1977) can be used to avoid prescribing the sensitivity of g_s to moisture stress and/or atmospheric drivers like vapor pressure deficit (Prentice et al., 2014; Sperry et al., 2017; Wolf et al., 2016).

Akin to the proliferation of empirical g_s models—see Damour et al. (2010) for a review of >30 empirical g_s formulations—at least 10 new stomatal optimization schemes have been proposed in the past 5 years (Anderegg et al., 2018; Dewar et al., 2017; Eller et al., 2018; Huang et al., 2018; Lu et al., 2020; Mrad et al., 2019; Novick et al., 2016; Sperry et al., 2017; Wolf et al., 2016). Some g_s optimization models have been compared (Anderegg et al., 2018; Bassiouni & Vico, 2021; Dewar et al., 2017; Mrad et al., 2019; Novick et al., 2016; Wang et al., 2020) or evaluated against empirical analogs (Eller et al., 2018, 2020; Sabot et al., 2020; Venturas et al., 2018; Wang et al., 2019). Yet, comparisons rarely control for both univariate model sensitivity (e.g., isolating soil moisture impacts from vapor pressure deficit impacts) and agreement with observations (but see Novick et al., 2016). Importantly, only one g_s optimization model (Sperry et al., 2017) has been evaluated against carbon and water fluxes across leaf, plant, and ecosystem scales (Sabot et al., 2020; Venturas et al., 2018; Wang et al., 2019).

This lack of systematic multi-model evaluation could be due to difficulties ranging from accessing existing code, to coding new schemes into more sophisticated land surface models (LSMs) that account for energy balance and water-carbon feedbacks. In the first extensive comparison of g_s optimization models, Wang et al. (2020) proposed a unified mathematical framework relying on seven theoretical criteria to probe the merits of different optimization models. However, the authors did not assess whether, and to what extent, new optimality-based
<table><caption>**Table 1**
Summary of the Stomatal Conductance (g_s) Schemes</caption><thead><tr><th>Approach</th><th>g<sub>s</sub> form or optimization objective<sup>a</sup></th><th>Model name</th><th>Conductance parameter(s)</th><th>Reduction/ cost parameter(s)<sup>a</sup></th><th>Specific equation(s)</th><th>Main reference(s)</th></tr></thead><tbody><tr><td rowspan="3">Empirical</td><td rowspan="3">g<sub>s</sub> ~ RF</td><td>Medlyn</td><td>g<sub>1,Med</sub></td><td>s<sub>Med</sub></td><td>Equations 4 and 5</td><td>De Kauwe et al. (2015) and Medlyn et al. (2011)</td></tr><tr><td>Tuzet</td><td>k<sub>max</sub>, g<sub>1,Tuz</sub></td><td>Ψ<sub>ref</sub>, s<sub>Tuz</sub><sup>b</sup></td><td>Equations 6 and 7</td><td>Duursma and Medlyn (2012) and Tuzet et al. (2003)</td></tr><tr><td>Eller</td><td>k<sub>max</sub></td><td>–</td><td>Equations 8 and 9</td><td>Eller et al. (2020)</td></tr><tr><td rowspan="12">Optimization</td><td rowspan="3">Avoiding primary impairment to water flow: max(A<sub>n</sub> – CT)</td><td>WUE-f<sub>ψ</sub></td><td>k<sub>max</sub></td><td>λ</td><td>Equation 11</td><td>Wolf et al. (2016)</td></tr><tr><td>CMax</td><td>k<sub>max</sub></td><td>a, b</td><td>Equation 13</td><td>Anderegg et al. (2018) and Wolf et al. (2016)</td></tr><tr><td>ProfitMax</td><td>k<sub>max</sub></td><td>–</td><td>Equation 14</td><td>Sperry et al. (2017)</td></tr><tr><td rowspan="4">Avoiding primary impairment to carbon assimilation: max(A<sub>n</sub> × RF)</td><td>CGain</td><td>k<sub>max</sub></td><td>ω</td><td>Equation 15</td><td>Lu et al. (2020)</td></tr><tr><td>SOX<sub>opt</sub></td><td>k<sub>max</sub></td><td>–</td><td>Equation 16</td><td>Eller et al. (2018)</td></tr><tr><td>ProfitMax2</td><td>k<sub>max</sub></td><td>–</td><td>Equation 17</td><td>Wang et al. (2020)</td></tr><tr><td>LeastCost</td><td>k<sub>max</sub></td><td>η</td><td>Equation 18</td><td>Prentice et al. (2014)</td></tr><tr><td rowspan="2">Non-stomatal optimization: max(A<sub>n</sub>), with internal RF on V<sub>cmax</sub> or C<sub>i</sub></td><td>CAP</td><td>k<sub>max</sub></td><td>Ψ<sub>φ,lim</sub></td><td>Equations 19–21</td><td>Dewar et al. (2017)</td></tr><tr><td>MES</td><td>k<sub>max</sub></td><td>Ψ<sub>φ,lim</sub></td><td>Equations 19, 20, and 22</td><td></td></tr></tbody></table>

developments are advances on existing empirical approaches, nor did they account for leaf-atmosphere feedbacks when assessing model performance.

Here, we compare and evaluate 12 g_s models within a single framework. We ask:

1. How well constrained are g_s scheme parameterizations?

2. Do model behaviors differ in response to abiotic variations?

3. How well do the g_s schemes capture observed leaf-level fluxes and hydraulic status across a large environmental gradient in moisture stress and atmospheric dryness?

We aim to identify the g_s schemes that improve predictability for use in LSMs (e.g., during periods of water stress, at high vapor pressure deficit) but, also, any aspects of those schemes that might limit current application.

## 2. Methods

We addressed our first two questions theoretically, by harmonizing 11 g_s models to a single reference model, the Medlyn et al. (2011) model. The Medlyn model was chosen as a baseline expectation of model behavior because it is widely used in climate models (ACCESS-ESM1.3, Kala et al., 2015; CESM2, Lawrence et al., 2019), and it is well understood (i.e., extensively evaluated in the literature), so we can readily diagnose assumption-driven model-to-reference behavioral departures (De Kauwe et al., 2013; Medlyn et al., 2015). Our experimental setup parallels the “perfect model” framework used in climate predictability studies (Hawkins et al., 2011), except our reference model is not intended as a truth. We address our third question by applying all 12 models to observational records of leaf-level gas exchange from 15 tree species spanning a mean annual precipitation (MAP) range of 363–3,286 mm yr<sup>-1</sup>. A visual summary of the model experiments (Sections 2.2 and 2.3) is presented in Figure S1 of Supporting Information S1.

To distinguish between the empirical and optimization approaches, we refer to the former by author name and the latter by optimality criterion. In this study, the approaches are labeled “empirical” if they explicate a functional form for g_s (i.e., the sensitivity of g_s to both independent and dependent variables is assumed), whereas they are labeled “optimization approaches” if they maximize/minimize function such that the sensitivity of g_s emerges from the optimization criteria. Table 1 provides a summary of the g_s schemes.

<footnote><sup>a</sup>RF stands for reduction factor and CT for cost term. <sup>b</sup>Parameter commonly set to an arbitrary value, that is, s<sub>Tuz</sub> is often set to 2 for the Tuzet model.</footnote>
## 2.1. Modeling Framework

We implemented the $g_s$ schemes within a simplified LSM (Sabot et al., 2020) that broadly emulates the Community Atmosphere-Biosphere Land Exchange LSM (CABLEv2.0; De Kauwe et al., 2015). This framework embeds canopy micrometeorology and coupling between energy balance, transpiration, and photosynthesis. We standardized model responses by prescribing soil moisture, thereby removing the feedback between transpiration and root-zone soil moisture. We also opted to ignore leaf-to-canopy scaling, which allows a direct comparison to leaf-level observations, whilst avoiding LSM-specific scaling assumptions (Rogers et al., 2017).

### 2.1.1. Leaf Level Fluxes and Energy Balance

The net rate of carbon assimilation, $A_n$ (µmol m⁻² s⁻¹), is simulated using the Farquhar et al. (1980) photosynthesis model (Text S1 in Supporting Information S1).

Leaf transpiration (E; mmol m⁻² s⁻¹) is assumed to meet the atmospheric demand for water vapor on an instantaneous basis:

$$E = 10^3 \frac{g_s g_b D_t}{(g_s + g_b) P_{atm}} \quad (1)$$

where $g_s$ is expressed on a H₂O basis in mol m⁻² s⁻¹, $g_b$ (mol m⁻² s⁻¹) is the leaf boundary layer conductance to water vapor, $D_t$ (kPa) is the leaf-to-air vapor pressure deficit, $P_{atm}$ (kPa) is atmospheric pressure, and 10³ converts to mmol.

For the empirical schemes (3 models), leaf-level variables (e.g., leaf temperature) are obtained iteratively until the leaf energy balance converges. The optimization schemes (nine models) avoid iteration (except SOXopt; see below) by using two-dimensional (water transport vs. photosynthetic) systems of energy-balanced equations to simultaneously optimize $g_s$, [CO₂] in the leaf intercellular air spaces ($C_i$), and leaf water potential ($Ψ_l$). Note, this implementation difference does not impact the simulations. For more information, refer to the main text and methods S2 of Sabot et al. (2020).

### 2.1.2. Leaf Water Potential

$Ψ_l$ (MPa) is computed by substitution of Equation 1 and the steady-state supply of water from the roots to the leaves (Sperry & Love, 2015):

$$E = \int_{\Psi_s}^{\Psi_l} k_{\Psi} d\Psi \quad (2)$$

where $Ψ$ is the varying water potential between the root-zone soil water potential $Ψ_s$ (MPa) and $Ψ_l$, and $k_Ψ$ (mmol m⁻² s⁻¹ MPa⁻¹) is the associated hydraulic conductance. $Ψ$ cannot drop below $Ψ_{crit}$, a threshold indicative of critical xylem failure at $k_{crit}$, and set to match a 95% loss of hydraulic conductivity (i.e., $P_{95}$).

$k_Ψ$ is downregulated by a cumulative Weibull distribution representing plant vulnerability to cavitation (Neufeld et al., 1992; Sperry et al., 2017):

$$k_{\Psi} = k_{\max} e^{-(\frac{|\Psi|}{s_1})^{s_2}} \quad (3)$$

where $k_{max}$ (mmol m⁻² s⁻¹ MPa⁻¹) is the maximum root-zone-to-leaf hydraulic conductance, and $s_1$ (MPa) and $s_2$ (unitless) are the sensitivity and shape of the vulnerability curve derived from measured percentage conductivity loss (Text S2 in Supporting Information S1).

In segmented representations of plant hydraulics, $k_{max}$ (and $s_1$ and $s_2$) may vary between the rhizosphere, roots, stem, etc. Parameterizing this segmentation presents a challenge for global models, as it requires additional inputs that are not readily available, so we combine all elements into a single hydraulic conductor connecting the root-zone to the leaves.
### 2.1.3. Empirical Schemes (Three Models)

#### The Medlyn Model

The Medlyn $g_s$ model (Medlyn et al., 2011), reworked to account for water stress (De Kauwe et al., 2015), is an analytical approximation of the water use efficiency hypothesis (cf., WUEH; Cowan & Farquhar, 1977):

$$g_s \approx g_0 + 1.57 \left( 1 + \frac{g_{1,Med} \beta}{\sqrt{D_l}} \right) \frac{A_n}{C_s} \quad (4)$$

where $C_s$ (µmol mol⁻¹) is [CO₂] at the leaf surface, 1.57 converts from conductance to CO₂ to conductance to water vapor, $g_0$ (mol m⁻² s⁻¹) is the residual conductance to water vapor (assumed negligible), $g_{1,Med}$ (kPa⁰.⁵) is the slope of the sensitivity of $g_s$ to $A_n$, and $\beta$ (unitless) is an empirical moisture stress factor represented by an exponential dependency on $\Psi_l$ at predawn (i.e., $\Psi_{l,pd}$; Cowan, 1982; Yang et al., 2019; Zhou et al., 2013):

$$\beta = \begin{cases} 1 & \text{if } \Psi_s \ge \Psi_{fc} \\ e^{s_{Med} \Psi_{l,pd}} & \text{if } \Psi_s < \Psi_{fc} \end{cases} \quad (5)$$

where $s_{Med}$ (MPa⁻¹) sets the sensitivity of the stomates to $\Psi_{l,pd}$, $\Psi_{fc}$ is the root-zone water potential at field capacity, and with $\Psi_{l,pd} \sim \Psi_s$ assuming no night-time transpiration.

The choice to attenuate either the ratio of $A_n$; $g_s$ (i.e., the “intrinsic” water use efficiency) or E during periods of water stress follows several LSMs, but we acknowledge that alternative approaches exist (cf., Kennedy et al., 2019).

#### The Tuzet Model

Tuzet et al. (2003) empirically linked $g_s$ and $\Psi_l$, which Duursma and Medlyn (2012) later simplified to:

$$g_s \approx g_0 + 1.57 g_{1,Tuz} \zeta \frac{A_n}{C_s} \quad (6)$$

where $g_{1,Tuz}$ (unitless) is an empirical slope parameter and $\zeta$ (unitless) characterizes stomatal responses to $\Psi_l$:

$$\zeta = \frac{1 + e^{s_{Tuz} \Psi_{ref}}}{1 + e^{s_{Tuz} (\Psi_{ref} \Psi_l)}} \quad (7)$$

where $\Psi_{ref}$ (MPa) is a reference water potential and $s_{Tuz}$ (MPa⁻¹) is the sensitivity of the stomates to variations in $\Psi_l$.

#### The Eller Model

Eller et al. (2020) proposed an analytical approximation to the stomatal optimization based on xylem hydraulics (cf., the SOXₐ₉₅t model below), establishing an empirical link between $g_s$ and plant hydraulic vulnerability:

$$g_s \approx 1.57 \frac{A_n(C_s) - A_n(C_{i,col})}{C_s - C_{i,col}} \left( \sqrt{0.25 + \frac{(C_s - C_{i,col})\xi}{1.57(A_n(C_s) - A_n(C_{i,col}))}} - 0.5 \right) \quad (8)$$

where $C_{i,col}$ (µmol mol⁻¹) is the co-limitation $C_i$—for which biochemical rates of assimilation are equal—and $\xi$ (mol m⁻² s⁻¹) is a reduction factor arising from xylem conductance loss:

$$\xi = \frac{k_{\Psi_{x,pd}}^2 (\Psi_{x,pd} - P_{50}) P_{atm}}{10^3 (k_{\Psi_{x,pd}} - k_{P_{50}}) D_l} \quad (9)$$

where $k_\Psi$ is given by Equation 3, $P_{50}$ (MPa) is the water potential drop at 50% loss of xylem hydraulic conductivity, $10^3$ converts to mmol, and $\Psi_{x,pd}$ is leaf xylem water pressure at predawn, with $\Psi_{x,pd} \sim \Psi_s - 10^{-6} \rho g h$ following Archimede's principle, where $\rho$ (kg m⁻²) is the density of water, $g$ (m s⁻²) is gravitational acceleration, $h$ (m) is canopy height, and $10^{-6}$ converts to MPa.
### 2.1.4. Optimization Schemes (Nine Models)

The optimization models hypothesize that leaf gas exchange downregulation originates from instantaneous $g_s$ controls to avoid primary impairments to water flow (Sperry et al., 2017; Wolf et al., 2016) or assimilation (Eller et al., 2018; Prentice et al., 2014), or that it originates from non-stomatal limitations on photosynthetic function induced by water stress (Dewar et al., 2017). Critically, whether impairment is defined as a function of a dynamic variable or not (e.g., $\Psi_l$ vs. $k_\psi$) and whether it is expressed as a reduction factor (multiplied) or a cost term (subtracted) may affect model responses to changes in the environment.

In the optimization schemes, $C_l$ is obtained by solving a system comprising $A_n$ from the biochemical photosynthetic model (Text S1 in Supporting Information S1) and the diffusive supply of CO$_2$:

$$A_n = \frac{g_s g_b}{1.57 \left( \frac{1.35}{1.57} g_s + g_b \right)} (C_a - C_l) \quad (10)$$

where $C_a$ (µmol mol$^{-1}$) is ambient [CO$_2$], and 1.35 converts from boundary layer conductance to CO$_2$ to boundary layer conductance to water vapor.

Note, our operational implementation uses the optimization criteria forms (as below) instead of the derivative forms (Figure S2 and Text S8 in Supporting Information S1), except for the CMax model (see below).

### 2.1.5. Stomatal Optimizations Avoiding Primary Impairment to Water Flow (Four Models)

#### The WUE-$f_{\psi_l}$ Model

The long-standing water use efficiency hypothesis (WUEH; Cowan & Farquhar, 1977) posits that plants adjust their carbon uptake given a “carbon cost of water loss” ($\lambda$; µmol CO$_2$ mmol$^{-1}$ H$_2$O), by regulation of their stomates to maximize:

$$\max (A_n - \lambda E) \quad (11)$$

However, the WUEH assumes $\lambda$ is constant over time and fails to describe how Equation 11 changes on the timescale over which soil water changes (Manzoni et al., 2011; Wong et al., 1985).

To address this issue, Wolf et al. (2016) linked the WUEH (Equation 11) with hydraulic function through $\Psi_l$, by equating Equations 1 and 2 and Equation 10 with the biochemical photosynthetic assimilation rates (Text S1 in Supporting Information S1), so leaf gas exchange is downregulated as E declines through falling $\Psi_s$ and $\Psi_l$ (Equation 2). This proposition keeps $\lambda$ constant but a priori unspecified, and we henceforth refer to it as the WUE-$f_{\psi_l}$ model.

#### The CMax Model

Wolf et al. (2016) also cast the carbon maximization hypothesis (CMax) as an alternative to the WUEH, replacing $\lambda E$ with a concave-up parabola that increases as $\Psi_l$ drops (Anderegg et al., 2018):

$$\max \left( A_n - \left( \frac{a}{2} |\Psi_l|^2 + b |\Psi_l| + c \right) \right) \quad (12)$$

where $a$ (µmol m$^{-2}$ s$^{-1}$ MPa$^{-2}$) and $b$ (µmol m$^{-2}$ s$^{-1}$ MPa$^{-1}$) and $c$ (µmol m$^{-2}$ s$^{-1}$) are a priori unknown coefficients of a parabola but knowing $c$ is unnecessary, as the maximization criterion can be satisfied by minimizing the derivative of Equation 12 with respect to $\Psi_l$:

$$\frac{\partial A_n}{\partial \Psi_l} - \left( a\Psi_l + b\frac{\Psi_l}{|\Psi_l|} \right) = 0 \quad (13)$$

$a$ can be interpreted as the carbon cost of a falling $\Psi_l$, whilst $lbl$ sets the minimum cost.

#### The ProfitMax Model

To avoid the need for unspecified parameters, Sperry et al. (2017) maximized profit (ProfitMax), that is, the net difference between relative carbon gain and relative hydraulic cost:
$$ \max \left( \frac{A_n}{A_{n,\max}} - \frac{k_{\psi_s} - k_{\psi_l}}{k_{\psi_s} - k_{\text{crit}}} \right) \qquad (14) $$

where $A_{n,\max}$ (µmol m⁻² s⁻¹) is the instantaneous maximum $A_n$ from Equation 10 over the range of possible $Ψ_l$, $g_s$, and $C_i$.

### The CGain Model

Lu et al. (2020) maximized plant net carbon gain (CGain) given a “carbon cost per recovered unit of xylem conductance” (ω; µmol m⁻² s⁻¹). The CGain hypothesis can account for lagged xylem recovery from embolism but we assume that refilling cavitated conduits comes at no extra cost (Text S4 in Supporting Information S1), to remain consistent across models:

$$ \max \left( A_n - \omega \frac{k_{\max} - k_{\psi_l}}{k_{\max}} \right) \qquad (15) $$

Here, ω is constant in time and a priori unspecified, akin to λ in the WUEH.

### 2.1.6. Stomatal Optimizations Avoiding Primary Impairment to Carbon Assimilation (Three Models)

#### The SOX<sub>opt</sub> Model

Inspired by the ProfitMax model, the stomatal optimization based on xylem hydraulics (SOX<sub>opt</sub>; Eller et al., 2018) maximizes carbon uptake whilst $A_n$ is directly impaired by a hydraulic reduction factor. This scheme does not rescale hydraulic cost as $Ψ_s$ drops (i.e., $k_{\max}$ is not evaluated at $Ψ_s$, unlike $k_{\psi_s}$ in ProfitMax), which removes the need for normalization of $A_n$:

$$ \max \left( A_n \frac{k_{\psi_l} - k_{\text{crit}}}{k_{\max} - k_{\text{crit}}} \right) \qquad (16) $$

Eller et al. (2018) made SOX<sub>opt</sub>'s system of equations one-dimensional (instead of two-dimensional as in other optimization schemes), by setting the leaf temperature, directly linking $Ψ_l$ to $A_n$, $C_i$, $g_s$, and $g_b$ through Equations 1, 2, and 10. In this implementation, SOX<sub>opt</sub> thus iterates leaf temperature until the leaf energy balance converges.

#### The ProfitMax2 Model

Instead of a conductance reduction factor, Wang et al. (2020) used a transpiration-based factor to increasingly impair A as E increases. The optimization criterion weights carbon uptake by proximity to hydraulic failure, maximizing carbon profits (ProfitMax2):

$$ \max \left( A_n \frac{E_{\text{crit}} - E}{E_{\text{crit}}} \right) \qquad (17) $$

where $E_{\text{crit}}$ (mmol m⁻² s⁻¹) is calculated from Equation 2 and $Ψ_{\text{crit}}$ triggers a 95% loss of hydraulic conductivity, with $Ψ_{\text{crit}} < Ψ < Ψ_s$.

#### The LeastCost Model

The least cost hypothesis (LeastCost; Prentice et al., 2014) posits that plants minimize the costs of maintaining both transpiration and carboxylation capacity, which is equivalent to maximizing $A_n$ reduced by a factor of the maintenance costs:

$$ \min \left( \frac{\eta E + V_{\text{cmax}}}{A_n} \right) \equiv \max \left( A_n \frac{1}{\eta E + V_{\text{cmax}}} \right) \qquad (18) $$

where η (µmol CO₂ mmol⁻¹ H₂O) is the carbon cost of photosynthetic proteins maintenance through transpiration, and $V_{\text{cmax}}$ (µmol m⁻² s⁻¹) is the maximum rate of carboxylation.

As in the WUE-fψ<sub>l</sub> model, leaf gas exchange is downregulated as E declines through falling $Ψ_s$ and $Ψ_l$ (Equation 2) and we assume η (a priori unspecified) constant over time.
### 2.1.7. Non-Stomatal Optimizations (The CAP and MES Models)

Non-stomatal limitations assume that gas exchange downregulation arise outside of the stomates, from biochemical limitations to photosynthesis or diffusion limitations into the mesophyll cells. Dewar et al. (2017) proposed an optimality framework where hydraulic stress curtails either carboxylation capacity (CAP) or mesophyll conductance (MES), so carbon uptake is maximized at an optimal $\Psi_l$.

$$ \max(A_n) \qquad (19) $$

Hydraulic stress is represented through a unitless linear factor:

$$ \varphi = \max \left( 0, 1 - \frac{\Psi_l}{\Psi_{\varphi, \lim}} \right) \qquad (20) $$

where $\Psi_{\varphi, \lim}$ (MPa) is an a priori unknown leaf water potential under which no photosynthesis can occur.

In CAP, hydraulic stress directly reduces $V_{cmax}$ and $J_{max}$ (µmol m⁻² s⁻¹), the maximum rate of electron-transport, affecting the biochemical photosynthetic assimilation rates (Equations S1.2 and S1.4 in Supporting Information S1):

$$ V_{cmax} \mapsto \varphi V_{cmax} \quad \& \quad J_{max} \mapsto \varphi J_{max} \qquad (21) $$

Instead, MES assumes that downregulation of mesophyll conductance from hydraulic stress yields inequity between [CO₂] in the chloroplast and $C_i$, so the assimilation rates (Equations S1.2 and S1.4 in Supporting Information S1) account for:

$$ C_i \mapsto \varphi(C_i - \Gamma^*) + \Gamma^* \qquad (22) $$

where $\Gamma^*$ (µmol mol⁻¹) is the photorespiration compensation point.

### 2.2. Theoretical Model Comparisons

To perform controlled model comparisons and experiments, we used a weather generator and created 4 weeks of synthetic half-hourly atmospheric forcing during the northern hemisphere's growing season (Figures S3a–S3c in Supporting Information S1; see code). Two synthetic soil moisture profiles were also created (Figure S3d in Supporting Information S1) to contrast well-watered conditions with moderately stressed conditions. Timeseries of $g_s$ associated to the synthetic atmospheric forcing were simulated using the Medlyn model (Figure S4 in Supporting Information S1), depending on soil moisture, and with $g_b$ set to infinity (i.e., assuming perfect leaf-atmosphere coupling).

The other $g_s$ models were calibrated to this idealized reference $g_s$ timeseries only (“Idealized calibrations”) before we perturbed the environmental variables to compare model-specific sensitivities (“Sensitivity analysis”). Calibrating to $g_s$ only (here and in the “Evaluation against observations”), follows standard LSM practice as these models then internally scale fluxes from the leaf to the canopy (De Kauwe et al., 2015; Fisher & Koven, 2020; Franks et al., 2017). Importantly, this approach does not preclude evaluating how representations of $g_s$ impact the prediction of other fluxes, like $A_n$ and $E$.

The model experiments aim to isolate model differences, so driving variables do not always co-vary realistically, for example, air temperature ($T_a$) and atmospheric vapor pressure deficit ($D_a$) are coupled, but photosynthetic photon flux density (PPFD) and soil moisture are not.

Table 2 summarizes the ranges of environmental conditions used. We assumed a canopy height of 20 m (this information is only used by the Eller model) and set the hydraulic vulnerability ($P_{50} = -3.1$ MPa, $P_{88} = -4.9$ MPa; Maherali et al., 2006) and photosynthetic ($V_{cmax,25} = 35$ µmol m⁻² s⁻¹, $J_{max,25} = 58.5$ µmol m⁻² s⁻¹; Ellsworth, 2000) traits and $g_{1,Med}$ (2.5 kPa⁰.⁵; Lin et al., 2015) to those of *Pinus taeda*. $s_{Med}$ was set to 2 MPa⁻¹. Tables S1 and S2 in Supporting Information S1 list all other model parameters.
<table><caption>**Table 2**

*Environmental Variable Ranges Used for the Idealized Model Calibrations and Sensitivity Analysis*</caption><thead><tr><th>Variable</th><th>Symbol</th><th>Unit</th><th>Calibrations</th><th>Sensitivity analysis<sup>a</sup></th></tr></thead><tbody><tr><td colspan="5"><strong>Atmospheric forcing</strong></td></tr><tr><td>Photosynthetic photon flux density</td><td>PPFD</td><td>µmol photon m<sup>-2</sup> s<sup>-1</sup></td><td>50.2–1,594</td><td>50–2,500</td></tr><tr><td>Air temperature</td><td>T<sub>a</sub></td><td>°C</td><td>3.1–30</td><td>2–40</td></tr><tr><td>Atmospheric vapor pressure deficit</td><td>D<sub>a</sub></td><td>kPa</td><td>0.5–3</td><td>0.1–10</td></tr><tr><td>Ambient [CO<sub>2</sub>]</td><td>C<sub>a</sub></td><td>ppm</td><td>410</td><td>250–900</td></tr><tr><td colspan="5"><strong>Soil moisture conditions</strong></td></tr><tr><td rowspan="2">Volumetric soil moisture</td><td rowspan="2">θ</td><td rowspan="2">m<sup>3</sup> m<sup>-3</sup></td><td>Wet</td><td>0.28–0.41</td></tr><tr><td>Stressed</td><td>0.12–0.37</td></tr><tr><td rowspan="2">Soil water potential</td><td rowspan="2">Ψ<sub>s</sub></td><td rowspan="2">-MPa</td><td>Wet</td><td>0.006–0.0008</td></tr><tr><td>Stressed</td><td>0.46–0.0014</td></tr></tbody></table>

### 2.2.1. Idealized Calibrations

#### 2.2.2. Calibrated Parameters

Within the LSM framework, we calibrated a maximum of two parameters per scheme (Table 1; see Table S3 in Supporting Information S1 for the sampled parameter spaces) to the reference $g_s$, which avoids overfitting and limits autocorrelation. Apart from the reference $g_s$ model (Medlyn), all the other $g_s$ schemes considered make use of a $k_{max}$ parameter (Equation 3). Although $k_{max}$ should be directly measurable (Mencuccini et al., 2019; see De Kauwe et al. [2020] for an application within a LSM), we chose to calibrate it because models might assume additional/different contributions to $k_{max}$ related to different sensitivities to environmental drivers (cf., Wang et al., 2020).

We assumed that the optimization schemes that avoid primary impairment to water flow (WUE-$f_{\psi_i}$, CMax, ProfitMax, and CGain) would all share similar $k_{max}$, owing to their optimality criteria of the form $max(A_n-CT)$ that implies cost terms (CT; proportional to $k_{max}$) comparable to $A_n$ in magnitude. Thus, upon calibrating the WUE-$f_{\psi_i}$, CMax, and CGain models, we set $k_{max}$ to the value obtained for the ProfitMax model (Equation 14). Note, also, that these three models' additional cost parameters (e.g., $\lambda$ in WUE-$f_{\psi_i}$) buffer behavioral departure from the ProfitMax's $k_{max}$. The Tuzet model adopted the ProfitMax $k_{max}$ here too (but not in the “Evaluation against observations”, see below), as the three-parameter combination of $k_{max} - g_{1,Tuz} - \Psi_{ref}$ is equifinal.

#### 2.2.3. Calibration Method

To answer “how well constrained are $g_s$ scheme parameterizations?,” we calibrated each scheme via seven non-linear least-square minimizers from the LMFIT python package (Newville et al., 2019), ranging from local solvers to global stochastic minimizers capable of handling complex multimodal objective functions (Table S4 in Supporting Information S1). Using a suite of solvers, instead of a single calibration method, avoids favoring specific model objective functions.

We compared calibrations under well-watered conditions (“wet” soil moisture profile) with moderately stressed conditions (“stressed” soil moisture profile); this allows us to test whether parameter predictability and model agreement differ with soil moisture availability. To further characterize solver skill and parameter predictability, we performed “wet” versus “stressed” calibrations on three randomly selected 7-day-long subsets within the weather-generated forcing.

Within data subsets, solver skill was quantified through the Bayesian information criterion (BIC) and given the proximity between the parameter value(s) estimated by each solver and the median parameter value(s) across

<footnote><sup>a</sup>The ranges selected are chosen to yield realistic model estimates.</footnote>
solutions. The best parameterizations (see Figure 1), determined by the most skilled solvers for the complete forcing timeseries, were then used to run idealized simulations (Appendix A and "Sensitivity analysis").

### 2.2.4. Sensitivity Analysis

To understand the physiological responses produced by the models—that is, to answer: “do model behaviors differ in response to abiotic variations?”—we carried out a variance-based global sensitivity analysis. We used Sobol's method (Sobol', 2001; see below), which relies on a probabilistic framework to quantify the sensitivity of a model's outputs over the full domain of variation of its inputs and is particularly robust when dealing with non-linear model responses (Marzban, 2013). Here, the method probes to what degree variability in key environmental variables (PPFD, $T_a$, $D_a$, $C_a$, and $\Psi_s$) affects modeled $g_s$, $\Psi_l$, and $C_i$.

The models were parameterized using the best parameterizations from the “wet” calibrations and forced with 1,008,000 combination sets of PPFD, $T_a$, $D_a$, $C_a$, and $\Psi_s$. The Saltelli cross-sampling method (Saltelli et al., 2007) implemented in the SALib python package (Herman & Usher, 2017) was used to generate forcing samples as uniformly distributed as possible (Table 2). Sobol' sensitivity indices (Saltelli et al., 2010; Sobol', 2001) were calculated for $g_s$, $\Psi_l$, and $C_i$. Variance in modeled variables is estimated via the quasi-Monte Carlo method before being apportioned to each individual driver (first-order index, S1), to synergies between pairs of drivers (second-order index, S2), or to each driver plus all their synergic terms with other drivers (total-order index, ST).

For example, if $g_s'$ S1 by $T_a$ is 30%, its S1 by $D_a$ 40%, and its S2 by $T_a-D_a$ 20%, then 90% of $g_s$ variability is caused by variability in $T_a$ and $D_a$ and their direct interactions. Inaccuracy in these two drivers will be responsible for most of the uncertainty in simulated $g_s$. By contrast, the ST is a summary measure of the overall contributions of one driver to the variation in a modeled variable. For instance, if $T_a$ interacts with other drivers than $D_a$ (if $g_s'$ S2 by $T_a-\Psi_s$ is 5%), then besides the first- and second-order contributions, the ST by $T_a$ will also account for third-order contributions of $T_a$ ($T_a-D_a-\Psi_s$ synergies) to the total variability in $g_s$.

## 2.3. Evaluation Against Observations

To answer our third question, we evaluated the models against leaf-level observations originating from Australia, France, Panama, and the USA. The original data sets include field observations for 26 sites × species combinations (Choat et al., 2006; Limousin et al., 2013; Martin-StPaul et al., 2012; Mitchell et al., 2009; Wu et al., 2020), as well as measurements for four species grown in a common garden experiment (Héroult et al., 2013). After processing the observations (Text S5 in Supporting Information S1), measured $g_s$, $E$, and $A_n$ were available for 16 sites × species, $\Psi_l$ for 15 sites × species, and $C_i$ for 12 sites × species. Observed PPFD, $T_a$, $D_a$, $C_a$, $g_b$, and $\Psi_{l,pd}$ (used as a proxy for $\Psi_s$) were also included as inputs to drive the $g_s$ schemes. All sites × species specifics are listed in Table 3 and Table S5 in Supporting Information S1.

We used the four most skilled solvers from our theoretical “Idealized calibrations” (Figure S5 in Supporting Information S1) to calibrate the 12 schemes to the observations of $g_s$, for each individual site × species. These calibrations were performed within the LSM framework, using measured leaf temperature instead of $T_a$ to avoid energy-balance assumptions at the calibration stage.

Calibrating the Tuzet model, however, required a different approach than that used for the idealized calibrations, as not only was there no unique link between $g_s$ and $\Psi_l$ over all environmental conditions, but the observations were also inconsistently distributed in time, with variable tree age and size classes sampled. The $\Psi_{ref}$ and $s_{TUZ}$ parameters (Equation 7) were directly fitted from the observed $g_s$ and $\Psi_l$ (Text S6, Figures S6a–S6o in Supporting Information S1), and $\Psi_{ref}$ and $s_{TUZ}$ were then pre-set upon calibration of both $g_{1,TUZ}$ and $k_{max}$ within the model framework. Unlike in the “Idealized calibrations,” calibration of a $k_{max}$ specific to the Tuzet model is needed because $\Psi_{ref}$ and $s_{TUZ}$ were not calibrated within the model framework.

Predictive performance was appraised using four metrics (Text S7 in Supporting Information S1). Pearson's correlation coefficient ($r$) evaluates the temporal relationship between modeled and observed variables. The Nash-Sutcliffe Efficiency index (NSE) measures a model's ability to outperform the average observation. The Mean Absolute Scaled Error (MASE) gauges the accuracy of a model's dynamics, by assessing whether its predictions are more skillful than preceding observations. Finally, the ranked Bayesian Information Criterion (rBIC) provides information on parsimonious accuracy within each site × species.
<table><caption>**Table 3**
Site and Observation Data Set Characteristics, Sample Size (N), and Hydraulic and Photosynthetic Traits Used to Parameterize the Models</caption><thead><tr><th rowspan="2">Site</th><th rowspan="2">Lat.</th><th rowspan="2">Lon.</th><th rowspan="2">MAP (mm)</th><th rowspan="2">Species</th><th rowspan="2">PPFD (µmol m<sup>-2</sup> s<sup>-1</sup>)</th><th rowspan="2">D<sub>d</sub> (kPa)</th><th rowspan="2">Ψ<sup>a</sup> (-MPa)</th><th rowspan="2">Data set ref.</th><th rowspan="2">P<sub>50</sub>, P<sub>88</sub> (-MPa)</th><th colspan="2">Hydraulic traits</th><th rowspan="2">V<sub>cmax,25</sub>, J<sub>cmax,25</sub><sup>b</sup> (µmol m<sup>-2</sup> s<sup>-1</sup>)</th><th rowspan="2">Photosynthetic traits ref.</th></tr><tr><th>ref.</th><th>ref.</th></tr></thead><tbody><tr><td>San Lorenzo</td><td>9°16′N</td><td>79°58′W</td><td>3,286</td><td>Carapa guianensis</td><td>36</td><td>159–2,001</td><td>0.5–1.5</td><td>0.9–0.48</td><td>Wu et al. (2020)</td><td>0.8 ± 0.1; 2.3 ± 0.2</td><td>Lopez et al. (2005)</td><td>41.6 ± 1.6; 90.7 ± 2.3</td><td>Rogers et al. (2019)</td></tr><tr><td>Protected Area, Panama<sup>b</sup></td><td></td><td></td><td></td><td>Tachigali versicolor</td><td>30</td><td>121–2,001</td><td>0.4–1.2</td><td>1.8–0.55</td><td></td><td>1.6; 4</td><td>Meinzer et al. (2008)</td><td>53.7 ± 1.8; 99.6 ± 1.6</td><td></td></tr><tr><td>Parque Natural Metro politano, Panama<sup>b</sup></td><td>8°59′N</td><td>79°32′W</td><td>1,826</td><td>Tocoyena pittieri</td><td>32</td><td>100–2,000</td><td>0.5–1.5</td><td>1–0.55</td><td></td><td>3.3 ± 0.1; 4.3 ± 0.3</td><td>Wu et al. (2020)</td><td>50.8 ± 1.6; 83.9 ± 1.2</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>Calyophyllum candidissimum</td><td>32</td><td>51–2,001</td><td>0.5–4.3</td><td>2.2–0.33</td><td></td><td>2.9 ± 0.1; 4.3 ± 0.1</td><td>Brodribb et al. (2003)</td><td>77.4 ± 5.2; 112.1 ± 2.5</td><td></td></tr><tr><td>Many Peaks Range, QLD, Australia<sup>c</sup></td><td>19°11′S</td><td>145°45′E</td><td>1,100</td><td>Alphitonia excelsa</td><td>160</td><td>93–1,854</td><td>1–6.2</td><td>2.2–0.2</td><td>Choat et al. (2006)</td><td>5.6; 9.1</td><td>Choat et al. (2003)</td><td>34.2 ± 3.7; 68.9</td><td>Choat et al. (2006)</td></tr><tr><td></td><td></td><td></td><td></td><td>Austromyrtus bidwillii</td><td>55</td><td>52–1,823</td><td>1.1–4.9</td><td>0.9–0.2</td><td></td><td>5.1; 12</td><td></td><td>23.5 ± 3.8; 47.4</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>Brachychiton australis</td><td>103</td><td>117–2,146</td><td>0.9–7</td><td>0.8–0.18</td><td></td><td>3.2; 5.2</td><td></td><td>48.4 ± 6.8; 97.6</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>Cochlospermum gillivraei</td><td>84</td><td>226–1,967</td><td>1.8–6.8</td><td>1.1–0.28</td><td></td><td>1.4; 2.1</td><td></td><td>42.3 ± 2.5; 85.3</td><td></td></tr><tr><td>Richmond, NSW, Australia<sup>d</sup></td><td>33°19′S</td><td>138°32′E</td><td>462</td><td>Eucalyptus cladocalyx</td><td>90</td><td>398–2,001</td><td>0.6–5.1</td><td>1.4–0.5</td><td>Héroult et al. (2013)</td><td>5.3 ± 0.4; 8.1 ± 0.9</td><td>Bourne et al. (2017)</td><td>98 ± 11; 160.7</td><td>Héroult et al. (2013)</td></tr><tr><td></td><td></td><td></td><td></td><td>Eucalyptus dumii</td><td>92</td><td>398–2,002</td><td>0.5–5.5</td><td>1.2–0.05</td><td></td><td>4.9 ± 0.2; 8.4 ± 1</td><td></td><td>115 ± 11; 188.6</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>Eucalyptus dumii</td><td>78</td><td>193–2,001</td><td>0.5–5.9</td><td>1.3–0.2</td><td></td><td>3.4 ± 0.2; 5.2 ± 0.6</td><td></td><td>96 ± 9; 157.4</td><td></td></tr><tr><td>Corrigin, WA.</td><td>33°03′S</td><td>138°75′E</td><td>361</td><td>Eucalyptus capitillosa</td><td>132</td><td>75–2,202</td><td>1.2–8.9</td><td>4.4–3.2</td><td>Mitchell et al. (2009)</td><td>3.1 ± 0.1; 5.7 ± 0.2</td><td>Mitchell (2009)</td><td>71.7; 117.5</td><td>de Boer et al. (2016)</td></tr><tr><td>Australia<sup>e</sup></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Sevilleta, NM, USA</td><td>34°23′N</td><td>106°31′W</td><td>363</td><td>Juniperus monosperma</td><td>187</td><td>2,000</td><td>0.5–6.7</td><td>5.6–0.85</td><td>Limousin et al. (2013)</td><td>9.4 ± 0.4; 11.8 ± 0.2</td><td>Hudson et al. (2018)</td><td>35.6 ± 1.5; 68.1 ± 1.7</td><td>Limousin et al. (2013)</td></tr><tr><td>Vic la Gar diole, France<sup>b</sup></td><td>43°29′N</td><td>3°45′E</td><td>650</td><td>Pinus edulis</td><td>139</td><td>2,000</td><td>0.7–6.2</td><td>3–0.8</td><td>Martin-StPaul et al. (2012)</td><td>4.4 ± 0.1; 5.4 ± 0.1</td><td>Martin-StPaul et al. (2017b)</td><td>17.4 ± 1.6; 36.1 ± 2.2</td><td>Martin-StPaul et al. (2012)</td></tr><tr><td>Puéchabon, France<sup>b</sup></td><td>43°44′N</td><td>3°35′E</td><td>907</td><td>Quercus ilex</td><td>45</td><td>1,227–1,668</td><td>1–4.2</td><td>3.4–0.35</td><td></td><td>6.9; 9.4</td><td></td><td>56.9 ± 6.7; 98.6 ± 4.2</td><td></td></tr></tbody></table>

<sup>a</sup>Whenever possible, standard errors are reported as per the original source, are inferred from figures in the study of origin, or are calculated off raw data. <sup>b</sup>V<sub>cmax,25</sub> and J<sub>cmax,25</sub> were fitted from A to C<sub>f</sub> curves (see Text S5 in Supporting Information S1). <sup>c</sup>V<sub>cmax,25</sub> was estimated using the “one-point method” for light-saturated rates of photosynthesis (PPFD >500 µmol m<sup>-2</sup> s<sup>-1</sup>) as in De Kauwe et al. (2016) and J<sub>cmax,25</sub> is taken as 2.02 × V<sub>cmax,25</sub> as per Kumarathunge et al. (2019) on average for rainforest species. <sup>d</sup>For species grown in the Richmond common garden experiment, the geographic coordinates and MAP are those of the meteorological station nearest to the seeds' place of origin (Héroult et al., 2013). J<sub>cmax,25</sub> was taken as 1.64 × V<sub>cmax,25</sub> as per Kumarathunge et al. (2019) on average for *Eucalyptus* species. <sup>e</sup>The V<sub>cmax,25</sub> is that of a relative, *Eucalyptus wandoo*, and J<sub>cmax,25</sub> was taken as 1.64 × V<sub>cmax,25</sub> as per Kumarathunge et al. (2019) on average for *Eucalyptus* species.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Spread in calibrated model parameters under well-watered conditions (i.e., “Wet”) and mild soil moisture stress (i.e., “Stressed”). Parameter names appear to the right of the half violin plots and corresponding model names to their left; the models appear from the least (bottom) to the most (top) uncertain in their parameter estimates. Each half violin plot comprises 28 parameter estimates, obtained by calibration to four synthetic data sets (one 4-week data set plus three 7-day data sets) using seven calibration methods. Parameter values are normalized to the median of each calibration data set. Hatched overlays highlight the parameters from the three most skilled calibration methods for each model, when those successfully constrain estimates. The best parameter estimates among those calibrated on the 4-week data set are represented by a star. Note, the scale is non-linear.</figcaption></figure>

### 3. Results

We first evaluate our idealized model calibrations, then explore theoretical model sensitivities, and finally examine whether and/or which theoretical findings are supported by observations.

#### 3.1. How Well Constrained Are $g_s$ Scheme Parameterizations?

Figure 1 shows the spread in parameter values calibrated under “wet” and “stressed” soil moisture conditions, normalized by the median parameter value within each calibration subset. Where the numerical solvers converge, the parameters are well defined and the plots display a single mode (e.g., ProfitMax, $SOX_{opt}$), typically within ±10% of the median parameter value, and with standard errors <10%. Where the solvers do not converge, there can be large spreads in the distributions: these reveal multimodal problems (CAP, MES), low influence parameters ($b$ in CMax), and/or recurring failure from a solver (e.g., LeastCost, $WUE-f\psi_l$) or a class of solvers (e.g., stuck in local minimums).

Regardless of soil moisture conditions, parameters that set the magnitude of $g_s$ and $E$ (e.g., $k_{max}$ in ProfitMax, $SOX_{opt}$, ProfitMax2, Eller; $g_{1,Tuz}$ in Tuzet) are identifiable across solvers, with few outliers to none. In CAP and MES, though, $k_{max}$ is not readily constrained, perhaps because of empirical dependencies (see below) and because reductions in $g_s$ and $E$ are not directly driven by stomatal responses.

The Tuzet model exemplifies a set of well-constrained empirical parameters. Contingent upon a careful solver selection (cf., yellow hatched overlay), the optimization schemes' carbon costs ($\omega$ in CGain, $\eta$ in LeastCost, $\lambda$ in $WUE-f\psi_l$, $a$ in CMax) are also among empirical parameters that can be constrained adequately. By contrast, parameters that empirically link to $\Psi_l$ in optimization schemes—$\Psi_{\phi,lim}$ in CAP and MES, and $b$ (but not $a$) in CMax—are difficult to calibrate, and range between >0.05–20 times the median parameter values.

To address this issue, CMax's $b$ could be fixed or omitted, given its little influence (Zenes et al., 2020). However, this cannot be advocated for CAP and MES whose parameter pairs often show equifinality and high correlation (from covariance matrices: c. 0.8 in >¾ of CAP's parameter estimates and c. 0.7 in ½ of MES"), thus representing trade-offs in optimization strategies (Appendix A; cf., Figure A1 vs. A2).

We gain two important insights from these results. First, combining optimality principles with empirically defined functions risks reducing predictable behavior by adding degrees of freedom. Second, numerical solver choice matters to the calibration of multi-parameter optimization models. Adopting global or Bayesian solvers (as opposed to commonly used local solvers) would be preferable in those cases.

#### 3.2. Do Model Behaviors Differ in Response to Abiotic Variations?

Figure 2 shows the (mostly out-of-sample) influence of atmospheric drivers and soil moisture ($\Psi_s$) on $g_s$, $\Psi_l$, and $C_i$ given by total-order Sobol' sensitivity indices. The total-order Sobol' sensitivity indices align with the first-order indices (Figure S7 in Supporting Information S1) and give a fuller picture of model variance.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Total-order Sobol' sensitivity indices of the stomatal conductance ($g_s$), leaf water potential ($\Psi_l$), and CO$_2$ concentration in the leaf intercellular air spaces ($C_i$) to variability in environmental drivers for the 12 models parameterized under well-watered conditions. The environmental drivers are: atmospheric vapor pressure deficit ($D_a$); soil water potential ($\Psi_s$); ambient CO$_2$ concentration ($C_a$); photosynthetic photon flux density; and air temperature ($T_a$). The concentric circles mark 0.25 increments on a scale of 0-0.75, with 0 signifying no influence and 0.75 high influence; the axes extend up to 1.</figcaption></figure>

In Figure 2 and Figure S7 in Supporting Information S1, where the polynomials have a limited coverage over a given axis, variance is low because model responses are highly deterministic. For example, model behaviors are largely unaffected by variability in PPFD (S1 < 0.04, ST < 0.08), which only yields variability in $C_i$, $g_s$, and $\Psi_l$ when photosynthesis is limited by ribulose-1,5-bisphosphate regeneration (c. 23% of occurrences, on average across models; see Text S1 in Supporting Information S1 for the biochemical limitation rates of assimilation). Even under these conditions, model responses to PPFD are monotonic (Equations S1.4 and S1.5 in Supporting Information S1) and generate little feedback, except immediately near the CO$_2$ compensation point and the co-limitation point where biochemical rates of $A_n$ are equal.

Inclusion of plant hydraulics (i.e., excluding the Medlyn model) confers the models the following features: (a) variability in $g_s$ is primarily driven by variability in $D_a$, with smaller contributions from $T_a$, $\Psi_s$, and $C_a$ that vary by model; and (b) variability in $C_a$, and to a lesser extent in $D_a$ and $T_a$, induces variability in $C_i$. By contrast, using a $\beta$ downregulation factor in Medlyn leads to $\Psi_s$ having a major impact on both $g_s$ and $C_i$.

The picture is more contrasted when looking at the $\Psi_l$ responses (cf., Appendix A). Eller, CMax, and LeastCost are mostly unaffected by variability in $\Psi_s$ but quite sensitive to variability in $D_a$ and $T_a$, and CMax and LeastCost are even affected by variability in $C_a$. Other models are sensitive to $\Psi_s$ in three ways: (a) responses to variability in $\Psi_s$ are smaller than to variability in $D_a$ and/or $T_a$ (Tuzet, WUE-$f_{\Psi_l}$); (b) sensitivity to $\Psi_s$ either slightly dominates
or is on par with sensitivity to other drivers (ProfitMax, CGain, SOX<sub>opt</sub>); and (c) variability is primarily driven by Ψ<sub>s</sub> (ProfitMax2, CAP, MES). We expect that the latter scenario would become prevalent if Ψ<sub>s</sub> were dropping further than the range tested here, and it should happen sooner for scenario (b) than scenario (a).

Sensitivity to Ψ<sub>s</sub> appears as the root-cause for inter-model differences in exerted hydraulic controls (i.e., high sensitivity to Ψ<sub>s</sub> signals tight hydraulic controls, and vice versa; see Appendix A). We also note that contrasting early morning/evening versus afternoon g<sub>s</sub> (Appendix A) most likely arises from high sensitivity to D<sub>a</sub> (and T<sub>a</sub>), with little contribution from PPFD.

### 3.3. How Well Do the g<sub>s</sub> Schemes Capture Observed Leaf-Level Fluxes and Hydraulic Status?

#### 3.3.1. Predictive Performance

We used four metrics to assess the modeled variables against observations. Pooling all sites × species, the Profit-Max model reproduces observed behaviors better than observational averages (highest overall NSE, Figure 3a), possesses accurate dynamics (lowest overall MASE, Figure 3b), and is highly parsimonious (lowest overall rBIC, Figure 3d), whilst the MES model shows high correlation with the observations (highest overall r, Figure 3c). It is worth noting that the demanding MASE benchmark of 1—that marks one-step forecasts of the observations— is met by four models (MES, LeastCost, SOX<sub>opt</sub>, ProfitMax) for g<sub>s</sub> across sites × species. Additionally, half the new models (Eller, ProfitMax, SOX<sub>opt</sub>, ProfitMax2, LeastCost, MES) outperform the well-established Medlyn model (average NSE >0.02 and >0.32, r > 0.42 and >0.53, MASE <1.25 and <1.06, and rBIC < 0.46 and <0.43, over the gas exchange variables and for g<sub>s</sub>, respectively), whilst an additional model displays equivalent performance (CGain).

Given the models were calibrated to observations of g<sub>s</sub> only, evaluating them for their ability to predict other variables degrades performance, even though some individual site × species are highly predictable (e.g., Figure S8 in Supporting Information S1; cf., Figure S9 in Supporting Information S1). Performance degradation can be stark in the tropics, where the models are less skilled at matching the g<sub>s</sub> behavior than in the extra-tropics (Figure 4). Finally, among variables, Ψ<sub>t</sub> is responsible for the largest performance degradation, particularly for the Eller and LeastCost models, which we examine below.

#### 3.3.2. Functional Relationships

Across forest types, Ψ<sub>t</sub> estimates are mostly within observational ranges (Figure 5), but not for the Eller and LeastCost models (and WUE-fψ<sub>t</sub>, to a degree) which simulate unrealistically low Ψ<sub>t</sub> as ecosystems get drier, confirming our theoretical expectation of weak hydraulic controls in these models. Other models are either prone to overestimation of Ψ<sub>t</sub> at the wetter sites (Tuzet, ProfitMax2, CAP, MES) or to underestimation at the xeric sites (CMax, CGain, ProfitMax, SOX<sub>opt</sub>). Ψ<sub>t</sub> bias directly follows from bias in the parameterization of k<sub>max</sub> (Figure S6p in Supporting Information S1), so WUE-fψ<sub>t</sub>, CMax, and CGain share ProfitMax's bias direction (i.e., they use its k<sub>max</sub> and are also biased low at the xeric sites), with WUE-fψ<sub>t</sub> and CMax also displaying more variable errors (Figures 5d and 5e).

Predictions of E and A<sub>n</sub> by sites × species generally fall within observational ranges (Figure 6 and Figure S10 in Supporting Information S1), but the models fail to vary sufficiently, or accurately enough, throughout the whole range of E and A<sub>n</sub>. Indeed, each model occupies its own intrinsic WUE space, depending on its ability to regulate gas exchange with D<sub>a</sub> for example (Medlyn, Tuzet, WUE-fψ<sub>t</sub>, CMax and CGain do not capture the behavior at high D<sub>a</sub> in Figure 6d). Therefore, individual models can be biased with regards to E, A<sub>n</sub>, or both, and may keep one of the two variables within a narrow range whilst mostly varying the other. Although a proportion of the biases can be attributed to the forcing (at the Many Peaks Range, the forcing originates from the nearest weather station), to heterogeneity in the data sets, or to uncertain parameterization (not allowing seasonal variations in V<sub>cmax,25</sub> may cause A<sub>n</sub> estimates to be concentrated around the median observation in Figure 6a [Wu et al., 2019]; leaf morphology and flow rate could affect E), discerning the causes of model-specific biases requires examining the co-variation of g<sub>s</sub>, C<sub>i</sub> and Ψ<sub>t</sub>.

The models generally agree on the magnitude of C<sub>i</sub> per site × species (Figure 7 and Figure S11 in Supporting Information S1), noting that the models that simulate relatively lower C<sub>i</sub> simulate marginally lower A<sub>n</sub> (Eller
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Average model performance across sites × species for: *g*<sub>s</sub>, *C*<sub>i</sub>, *E*, *A*<sub>n</sub>, and *Ψ*<sub>*l*</sub> (in red); all variables excluding *Ψ*<sub>*l*</sub> (in dark blue), and *g*<sub>s</sub> only (in turquoise). Dots show performance averages and horizontal bars show standard errors (cropped for readability in panels [a and b]). The performance metrics are: (a) the Nash-Sutcliffe Efficiency index (NSE); (b) the Mean Absolute Scaled Error (MASE); (c) Pearson's correlation coefficient (*r*); (d) the ranked Bayesian Information Criterion (rBIC). For the NSE, positive values characterize models more skilled than the observed mean, and one perfect forecasts. For the MASE, values <1 identify models more skilled than a one-step forecast of the previous observation, and 0 the minimum forecast error. Yellow dashed lines show the NSE benchmark of 0 and the MASE benchmark of 1. For the rBIC, proximity to 0 indicates the best trade-off between model accuracy and model complexity. The arrows point toward better performance, with higher values meaning higher performance for the NSE and *r*, versus lower values for the MASE and rBIC. The models appear in order of least to most skilled relative to the Medlyn model (for turquoise and dark blue groups of variables), and then relative to each other on each side of the Medlyn model (for all groups of variables).</figcaption></figure>

compared to ProfitMax2 in Figures 6a and 7a), and vice versa. Similarly, relatively lower *E* is associated with relatively lower *g*<sub>s</sub> (LeastCost compared to Tuzet in Figures 6a and 7a), or vice versa. Whereas other models tend toward similar magnitude overestimation of *g*<sub>s</sub> and *C*<sub>i</sub>, or higher *g*<sub>s</sub> than *C*<sub>i</sub>, CAP and MES overestimate *C*<sub>i</sub> more than they do *g*<sub>s</sub>, in fact largely so given their estimates of *A*<sub>n</sub> (Figures 6 and 7; Figures S10 and S11 in Supporting Information S1; see also Appendix A).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Nash-Sutcliffe Efficiency index for the stomatal conductance estimates of 12 models across sites × species. Positive values characterize models more skilled than the observed mean and one perfect forecasts. From left to right, the models appear in order of least (WUE-fψᵢ) to most (LeastCost) skilled.</figcaption></figure>

There are four sets of characteristics that describe $g_s$ to $Ψ_i$ relationships: (a) sensitivities and ranges are consistent with the observations (CAP, MES in Figure 7d and Figures S12l-S12o in Supporting Information S1); (b) sensitivities are consistent with the observations but $Ψ_i$ are offset (ProfitMax, $SOX_{opt}$ in Figure 7d; Figures S12l, S12m, S12o, and S12p in Supporting Information S1); (c) models have the wrong sensitivity (ProfitMax2 in Figure 7d; Figures S12l, S12m, S12o, S12p in Supporting Information S1); (d) models have the wrong sensitivity and are offset (LeastCost in Figures 7c and 7d and Figure S12 in Supporting Information S1). Moreover, for WUE-$fψ_i$ and particularly CMax, $g_s(Ψ_i)$ can split into “branches” at low $g_s$ (below c. 0.5 × $g_{s,norm}$ in Figure 7d; Figures S12e, S12o, and S12p in Supporting Information S1), which suggests diverging $g_s$ closure at low $Ψ_s$ versus high $D_a$.

Average $g_s$ sensitivities to $D_a$ are within observational ranges across models and sites × species (Figure 8a). At low $D_a$, there is a relatively large spread in model behaviors, as per our idealized findings (Appendix A). At high $D_a$, schemes start to diverge from the average observations between 4 and 5 kPa (but CMax diverges from the observations over the whole range of $D_a$). Patterns of stomatal closure (Figure 8b) reveal that the modeled rates of $g_s$ closure are between thrice faster and six times slower than observed. Three models agree with the observed patterns of stomatal closure (CGain, $SOX_{opt}$, ProfitMax), whilst one closes prematurely at high $D_a$ (CMax, with $g_s$ closure c. 45% too high), and the remaining (majority) are prone to largely delayed closure (the stomates are c. 25%—75% too open). Consistency in simulation of $g_s(Ψ_i)$ compared to the observations (i.e., accurate
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Box and whisker plots (line, median; box, interquartile range; whiskers, quartiles ±1.5 times the interquartile ranges; dots, outliers) of the difference between modeled and observed leaf water potential ($Ψ_l$) for 11 models, where positive values indicate overestimation by the models. Site × species are grouped by forest type, from the wettest ecosystem (panel [a]) to the driest one (panel [e]), in the following order: (a) *Carapa guianensis*, *Tachigali versicolor*, *Tocoyena pittieri*, *Calycophyllum candidissimum*; (b) *Alphitonia excelsa*, *Austromyrtus bidwillii*, *Brachychiton australis*, *Cochlospermum gillivraei*; (c) *Eucalyptus dunnii*, *Eucalyptus saligna*, *Eucalyptus cladocalyx*; (d) *Quercus ilex* (Puéchabon), *Quercus ilex* (Vicla Gardiole); (e) *Juniperus monosperma*, *Pinus edulis*. Light gray shadings span ± the maximum range of observations from among individual species within a forest type, and darker gray shadings span ±10% of that range. The scale is symmetrically logarithmic.</figcaption></figure>

functional shape and variability, regardless of bias) appears key to representing stomatal closure as atmospheric aridity increases. For example, CGain, SOXopt, and ProfitMax are characterized by greater $Ψ_l$ biases than the Tuzet and ProfitMax2 models in Figures 5d and 5e, yet their functional forms better match the observed $g_s(Ψ_l)$ (Figure S12 in Supporting Information S1). This implies that Tuzet and ProfitMax2 might be compensating for an overestimated sensitivity to soil moisture by underestimation of the sensitivity to $D_a$. By contrast, CAP and
<figure><img src="image_9.png" /></figure>

**Figure 6.** Leaf-level estimates of carbon assimilation ($A_n$) and transpiration ($E$) from 12 models at a subset of wet (panel [a]) to xeric (panel [d]) site × species, compared to observations (light gray crosses). The box and whisker plots aligned with each panel show the observed distributions of $A_n$ and $E$ (line, median; box, interquartile range; whiskers, quartiles ±1.5 times the interquartile ranges). Dashed black lines represent observed behaviors fitted via a generalized additive model and split in three groups: (i) the overall average behavior (unlabeled); (ii) the average behavior when atmospheric vapor pressure deficit ($D_a$) is low; and (iii) the average behavior when $D_a$ is high. Both high and low $D_a$ thresholds were determined relative to the site. For the model estimates, point size is proportional to the number of observations per site × species, and transparency to density. Outliers were excluded by capping the modeled values to 2.5 times the observed maximum at each site.

MES's $g_s(\Psi_t)$ are consistent with the observations, so it is likely they keep $g_s$ high at high $D_a$ to prevent declines in $A_n$ in the absence of water stress.

Our “Evaluation against observations” highlights differences in internal trade-offs. For the non-stomatal models, compensation between overestimated $C_i$ and underestimated $A_n$ can result in high predictive performance. More generally, misrepresentation of the relationship between $g_s$ and $\Psi_t$ does not always degrade the simulation of gas exchange (LeastCost, Eller); however, it may lead to discrepant responses at high $D_a$.

## 4. Discussion

Climate models routinely simulate responses to a climate space without a historical analog (Reu et al., 2014), yet the LSMs included in climate models markedly diverge from observations as soil water availability becomes limiting (Martínez-de la Torre et al., 2019; Ukkola, De Kauwe, et al., 2016). This divergence is, in part, associated with uncertainties in soil hydraulic processes (Van Looy et al., 2017) and processes governing vegetation water stress (Medlyn et al., 2016). Stomatal optimization approaches that account for plant hydraulic functions hold promise for LSMs (Eller et al., 2020; Sabot et al., 2020) because they incorporate more predictive capacity than empirical model parameterizations based on historical behaviors and are simpler to parameterize than detailed representations of plant hydraulics (De Kauwe et al., 2020; Kennedy et al., 2019). In this study, we reviewed the
<figure><img src="image_10.png" /></figure>

**Figure 7.** Functional relationships between stomatal conductance ($g_s$), CO$_2$ concentration in the leaf intercellular air spaces ($C_i$), and leaf water potential ($\Psi_l$) predicted by 12 models at a wet (panels [a and c]) and a dry (panels [b and d]) site × species. Panels (a and b) show the encircled interquartile ranges of simulated to observed $g_s$ ratios against simulated to observed $C_i$ ratios, such that a “perfect” model would be concentrated at the intersection of the 1:1 lines. Panels (c and d) show the modeled decline in $g_s$ with decreasing $\Psi_l$, fitted via a generalized additive model, compared to the observations (light gray crosses). The functional forms were made comparable by normalizing $g_s$ by its model-specific maximum per site × species, and $\Psi_l$ by the critical leaf water potential indicative of total xylem failure ($P_{95}$ in this study). Curves were not fitted if they did not monotonically decrease (e.g., WUE-$f\psi_l$ and CMax in panel [d]), or where the models operate at, or beyond, the $P_{95}$ (e.g., the Eller model in panel [d]).

assumptions and mechanisms underpinning empirical and (non-)stomatal optimization hypotheses, with a view to guiding model development. Our theoretical results were determined using the Medlyn $g_s$ model as reference, but we expect them to be robust to alternative empirical reference choices (e.g., Ball et al., 1987; Leuning, 1995), with subtle differences in the sensitivity to $D_a$ (Medlyn et al., 2011).

## 4.1. Limits to Empirical Approaches

An argument can be made in favor of calibrating soil moisture stress functions in state-of-the-art LSMs, as opposed to changing their complex structure (Leplastrier, 2002; Raupach & Finnigan, 1988). Here, using an exponential dependency on $\Psi_{l,pd}$ (or $\Psi_s$) to down-regulate $g_s$ yielded reasonable skill when the Medlyn et al. (2011) model was evaluated against observations. Nonetheless, empirical moisture stress functions can hamper $g_s$ model performance (see notes S2 and Figure S6 of Sabot et al., 2020) or be implausible when applied outside their calibration sample (Verhoef & Egea, 2014). Finally, as the Medlyn $g_s$ model estimates are not contingent upon estimating $\Psi_l$, its prognostic use beyond the simulation of gas exchange is limited.

An alternative approach to using moisture stress factors is to make $g_s$ an empirical function of leaf water potential, where $\Psi_l$ is obtained by accounting for plant hydraulics (Tuzet et al., 2003). The problem with this approach is
<figure><img src="image_11.png" /><figcaption>**Figure 8.** A comparison of the sensitivity of 12 stomatal conductance ($g_s$) models to vapor deficit ($D_a$) across site × species. The driest two thirds of soil moisture were excluded within each site × species to isolate $D_a$ effects from soil moisture effects. Panel (a) shows the binned average effect of $D_a$ on $g_s$ (normalized by the observed maximum $g_s$ per site × species) over the range of observed $D_a$ (0.4–7 kPa), with each bin spanning (c) 0.4 kPa. The colored lines represent the different models as per (b). The darker gray shading shows the average observed gs, with lighted shading spanning ±1 S.D. around it. The boxed area defines the selection of $g_s$ data shown in (b) and includes a minimum $D_a$ of 2 kPa to account for models that predict relatively early $g_s$ closure. Panel (b) shows ridge plots of the differences between rates of simulated $g_s$ decline (i.e., simulated $g_s$ normalized by its model-specific maximum per site × species and scaled to the maximum observation) and rates of observed $g_s$ decline, relative to observed $g_s$ decline from within the boxed area from (a). Negative spread and peaks signify early stomatal closure compared to the observations, and vice versa.</figcaption></figure>

that there is not always a unique link between $g_s$ and $Ψ_l$ over all environmental conditions (Figures S6a–S6o in Supporting Information S1), so the Tuzet model risks: (a) simulating erroneous $g_s$ (and E) responses to $Ψ_l$ with increasing vapor pressure deficit (this study; Yang et al., 2019); and (b) introducing numerical instability through unrealistic feedbacks between leaf temperature, $D_l$, $Ψ_l$, and $g_s$.

Recently, Eller et al. (2020) proposed an empirical $g_s$ model based on xylem hydraulics which approximates an optimality criterion (i.e., $SOX_{opt}$; Eller et al., 2018). Here, the Eller model showed high skill in the simulation of $g_s$ and other leaf gas-exchange variables (Figure 3) but produced unrealistically low $Ψ_l$ (Figure 5). To approximate an analytical solution to $SOX_{opt}$, Eller et al. (2020) replaced $SOX_{opt}$'s dependency on $Ψ_l$ with a numerically simpler dependency on $0.5(Ψ_{l,pd} + Ψ_l)$, incidentally weakening its hydraulic reduction factor and impacting its analytical solution. Therefore, to achieve similar $Ψ_l$ as predicted by $SOX_{opt}$, all of the Eller model's plant hydraulic parameters (not only $k_{max}$, but also the vulnerability curve traits) need to be refitted (see Figure S3 of Eller et al. [2020]). At Puéchabon (a site also considered in this study), parameterizing the Eller model within the Joint UK LSM (Clark et al., 2011) resulted in a $P_{50}$ of –1.8 MPa (“BET-Te” in Table 2 of Eller et al. [2020]), whereas measured $P_{50}$ is between –7 and –4.5 MPa for *Quercus ilex* (Martin-StPaul et al., 2014, 2017b). Besides, the Eller model was derived from $SOX_{opt}$ assuming a different representation of plant vulnerability to cavitation than that given by Equation 3 (cf., Manzoni et al., 2013), which could have aggravated its tendency to underestimate $Ψ_l$ in our modeling framework.

## 4.2. Added Value and Potential Drawbacks of Optimal Schemes

Optimization approaches avoid reliance on empirical corrections to gas-exchange as water supply becomes limited, instead relying on measurable hydraulic traits (ProfitMax, $SOX_{opt}$, ProfitMax2), or traits plus parameters specific to the optimization hypothesis ($WUE-fΨ_l$, CMax, CGain, LeastCost, CAP, and MES). We found three models (CMax, CAP, and MES) to be difficult to constrain (Figure 1), and one model (CMax) to contain a low
influence parameter that could be fixed or omitted (cf., Zenes et al., 2020). At the ecosystem scale, Bassiouni and Vico (2021) confirmed that the less physiologically meaningful parameters of CMax were the most uncertain among those of three $g_s$ optimization models. From an operational perspective, constraining multi-parameter models can require over three times as many function evaluations than for one-parameter models. This makes “easy to calibrate” models attractive because global parameterizations would come at a considerable computational cost for multimodal optimization schemes and/or schemes sensitive to the calibration conditions (Figure 1).

Optimization was an improvement over the well-established Medlyn model for five optimization models evaluated against leaf-level observations (ProfitMax, SOX<sub>opt</sub>, ProfitMax2, LeastCost, and MES), and the CGain model was largely on par with the Medlyn model (Figure 3). Contrary to previous findings by Anderegg et al. (2018), CMax was not a large improvement over WUE-f<sub>ψ</sub><sub>l</sub> in our modeling framework. CMax was not among the top performing optimization schemes either, whereas the LeastCost model was, at odds with the findings of Wang et al. (2020). Possible explanations for apparent discrepancies between past findings and ours include the following:

1. The calibration methods are different. For instance, Wang et al. (2020) calibrated their models to studentized distributions of $A_n$, E, and $\Psi_l$. Instead, we calibrated our selected models to observed $g_s$ only, which avoids spreading errors across variables and could be critical for models less skilled at capturing $\Psi_l$.

2. Whenever possible, we implemented the optimization criteria forms instead of their derivative forms (Figure S2 and Text S8 in Supporting Information S1).

3. Our model simulations include energy-balance feedbacks between the atmosphere and the leaves, whereas previous studies used leaf-level forcing variables.

4. We evaluated the models' ability to match the observed leaf-level fluxes and water potentials independently, using four metrics that outline various aspects of model performance, as opposed to a single metric averaged across variables.

5. The evaluations use different model selections and data sets, for example, ours includes more *Eucalyptus* species.

Discrepancies aside, our results agree with Wang et al. (2020) in that WUE-f<sub>ψ</sub><sub>l</sub> and (mostly) CMax may behave unrealistically at high $D_a$. We also concur that ProfitMax and SOX<sub>opt</sub> perform to a similar standard, but we find that a slightly more conservative mode of stomatal regulation (i.e., marginally lower gas exchange under well-watered conditions) and less sensitivity to water stress (e.g., higher $\Psi_l$ in Figures 5d and 5e) give ProfitMax the edge over SOX<sub>opt</sub>.

### 4.3. Non-Stomatal Photosynthetic Limitations

The non-stomatal optimization approaches evaluated here (the CAP and MES models) resulted in lower photosynthetic estimates than observed (Figure 6 and Figure S10 in Supporting Information S1), through mechanisms challenged by the observations of co-occurring $g_s$ and $C_i$ (Figure S11 in Supporting Information S1). At the same time, both models proved highly skilled at representing the dynamics of $g_s$, particularly for the *Quercus ilex* sites × species (Figure S9c in Supporting Information S1) for which seasonally varying non-stomatal limitations have been demonstrated (Martin-StPaul et al., 2012). One possible explanation is that, despite clear evidence of apparent $V_{cmax,25}$ downregulation with drought (Keenan et al., 2010; Zhou et al., 2013), the evidence does not unequivocally support instantaneous linear downregulation of photosynthetic capacity with $\Psi_l$ (Wang et al., 2020). Additionally, robust inclusion of non-stomatal limitations in LSMs would require novel fitting of biochemical photosynthetic traits (i.e., in situ rather than apparent $V_{cmax,25}$; Bahar et al., 2018; Crous et al., 2013; Sun et al., 2014) as well as their temperature dependencies (Knauer et al., 2019). We argue that inclusion of non-stomatal limitations in LSMs through optimality principles would be premature, as our understanding of the timescales at which non-stomatal processes interact with plant hydraulic function, let alone how they arise (Yang et al., 2019), remains limited.
4.4. Estimating k_max

$k_{max}$ is a key parameter of hydraulics-enabled models that is measurable yet rarely documented in the plant hydraulics literature. In the absence of measured $k_{max}$, we calibrated model-specific $k_{max}$ parameters from $g_s$ data. These calibrated $k_{max}$ were not always plausible (i.e., outside ranges inferred from observed *E* and $\Psi_l$, as well as vulnerability curves; Figure S6p in Supporting Information S1) for Tuzet, Eller, LeastCost, CAP, and MES. Furthermore, $k_{max}$ appeared: (a) low at the xeric sites but adequate elsewhere (ProfitMax, SOX\_opt); (b) or reasonable at the xeric sites but high at the mesic sites (ProfitMax2, MES); (c) or mostly biased (Tuzet, Eller, LeastCost, CAP).

We opted to use the ProfitMax's $k_{max}$ when parameterizing the WUE-$f_{\psi_l}$, CMax, and CGain models (see “Idealized calibrations” for the rationale). An alternative (but more computationally costly) method would have been to calibrate each of these models' $k_{max}$. Testing (not shown) revealed that calibrated parameter spreads for the WUE-$f_{\psi_l}$ and CGain models would remain analogous to the LeastCost model's in Figure 1, whereas constraining CMax's parameters would be harder still. The predictive performance (Figures 3 and 4) and functional relationships (Figures 5–8) of CGain and CMax would see no substantial change in terms of the gas exchange variables. However, their $\Psi_l$ would risk further departure from the observations through less realistic $k_{max}$ (compared to the ranges shown in Figure S6p of Supporting Information S1). Finally, had we calibrated its $k_{max}$, WUE-$f_{\psi_l}$ would resemble LeastCost in all aspects. That is, WUE-$f_{\psi_l}$'s ability to capture the observations of gas exchange would increase at the expense of $\Psi_l$ and $k_{max}$.

The inclusion of reduction/cost parameters in optimality criteria (e.g., $\lambda$ in WUE-$f_{\psi_l}$, $\eta$ in LeastCost, $\Psi_{\phi,lim}$ in CAP or MES) should be considered carefully, as these additional parameters did not confer a consistent advantage for the simulation of gas exchange and/or degraded the simulation of $\Psi_l$ compared to one-parameter models which required calibrating $k_{max}$ only.

4.5. Future Challenges

Realistic estimates of water potential and embolism are needed for models to represent drought-related mortality (but see De Kauwe et al. (2020) and Venturas et al. (2020) for the associated challenges) and drought legacy effects. Ideally, models should link $\Psi_l$ to $g_s$ through hydraulic reduction factors/cost terms that rely on realistic $k_{max}$. Newly available sap flow measurements from the SAPFLUXNET database (Poyatos et al., 2021), together with existing hydraulic trait databases and aligned with leaf-level observations, offer an opportunity to infer relevant parameterizations of $k_{max}$ and xylem failure thresholds (i.e., $\Psi_{crit}$, $k_{crit}$) that could reduce model bias in $\Psi_l$ estimates. Explicit rhizosphere (Venturas et al., 2018; Wang et al., 2019) and/or symplastic limitations (De Cáceres et al., 2021) under dry conditions could also be explored, noting associated parameterization uncertainties (Xu & Trugman, 2021).

Adding complexity might not suffice where the relationship between $g_s$ and $\Psi_l$ is mischaracterized (e.g., unrealistically shaped) and/or where magnitude biases are very large (Eller, LeastCost, and WUE-$f_{\psi_l}$ to a lesser extent), so reformulation of optimality criteria cost terms/regulation factors might be necessary. Critically, uncertainty in the formulation of plant vulnerability to cavitation (Ogle et al., 2009), in hydraulic vulnerability measurements for long vessel species (Cochard et al., 2013), and in $\Psi_l$ measurements (which can be significantly decorrelated in time from gas-exchange measurements, especially when made in the field) should be accounted for: (a) when selecting parameterizations from the literature (and parameterizations could be constrained using less artifact-prone turgor loss data); (b) when resolving optimality criteria formulations; and (c) in model predictions of $\Psi_l$ and plant water status.

5. Conclusion

Over half the models compared performed well against leaf-level observations of gas exchange. However, care should be taken when analyzing performance metrics, as seemingly good performance can result from problematic internal compensation effects, and/or mischaracterization of model sensitivities to soil moisture and vapor pressure deficit. Only one of the better performing models required calibration of a parameter explicating stomatal sensitivity to the soil moisture conditions (Medlyn), confirming the viability of hydraulics-based stomatal
models for predicting vegetation responses in novel climate spaces. Notably, the better constrained optimization approaches (e.g., ProfitMax, SOX<sub>opt</sub>), both in terms of parameterization and functional representations, embed the potential for acclimation and plasticity in the plant hydraulic and photosynthetic responses (Sabot et al., 2020; Sperry et al., 2019), with a link to plant carbon allocation patterns (Potkay et al., 2021; Trugman et al., 2019) of major significance for the future of global coupled climate-vegetation modeling.

## Appendix A: Harmonized Model Behaviors

Figure A1 shows the average diurnal behavior for a suite of fluxes calibrated under “wet” conditions, highlighting how calibration does not ensure agreement throughout the day, despite producing comparable average estimates. Figure A2 shows the same thing under “stressed” conditions. Average daily g<sub>s</sub> (Figures A1a and A2a) are within ±10% of Medlyn across calibration conditions, but single day average differences can exceed 20% under “wet” conditions and reach up to c. 85% under “stressed” conditions.

<figure><img src="image_12.png" /></figure>

**Figure A1.** Average diurnal behavior of leaf-level water (panels a, c, and e) and carbon (panels b and d) fluxes predicted by the 12 models under well-watered conditions for 26 simulated data points per day. The Medlyn model's simulation of stomatal conductance (g<sub>s</sub>, panel a) was used as a reference to calibrate the other models. Other simulated fluxes shown are: (b) the CO<sub>2</sub> concentration in the leaf intercellular air spaces (C<sub>i</sub>) as a proportion of the ambient CO<sub>2</sub> concentration (C<sub>a</sub>); (c) the leaf water potential (Ψ<sub>l</sub>); (d) the net rate of carbon assimilation (A<sub>n</sub>); and (e) the transpiration (E).
<figure><img src="image_13.png" /></figure>

**Figure A2.** Average diurnal behavior of leaf-level water (panels a, c, and e) and carbon (panels b and d) fluxes predicted by the 12 models under mild soil moisture stress for 26 simulated data points per day. The Medlyn model's simulation of stomatal conductance ($g_s$, panel a) was used as a reference to calibrate the other models. Other simulated fluxes shown are: (b) the CO₂ concentration in the leaf intercellular air spaces ($C_i$) as a proportion of the ambient CO₂ concentration ($C_a$); (c) the leaf water potential ($Ψ_l$); (d) the net rate of carbon assimilation ($A_n$); and (e) the transpiration ($E$).

Compared to Medlyn in Figure A1a, five optimal models (CMax, ProfitMax, CGain, SOXopt, and MES) are characterized by higher and earlier $g_s$ peaks in the morning when $T_a$, $D_a$, and PPFD are low (see Figure S3 in Supporting Information S1), resulting in asymmetric biases between a few simulated points in the early morning/ evening and the rest of the day (>15% greater $g_s$ in SOXopt than in Medlyn before 10 a.m. and after 5 p.m. vs. <10% smaller $g_s$ otherwise; >5% greater $g_s$ in CMax before 10 a.m. and after 5 p.m. vs. >14% smaller $g_s$ otherwise). For WUE-$f\Psi_l$ and LeastCost, $g_s$ peaks occur later in the morning than for Medlyn. Throughout the day, the Eller model is the stediest, possibly because it assumes biochemical co-limitation (Equation 8), which smooths photosynthetic transition effects (see below), and because its downregulation factor depends on $Ψ_{x, pd}$ and $P_{50}$ rather than $Ψ_l$ (Equation 9), which limits $D_a$ effects.

CAP and MES excepted (see below), well-watered estimates of diurnal $C_i$ (Figure A1b) vary together with $g_s$, so higher $g_s$ relative to Medlyn also means higher $C_i$. Under “stressed” conditions, marginal adjustments in water use may result in higher $C_i$ versus lower $g_s$ than Medlyn, at various points of the day (Eller, ProfitMax2 in Figures A2a and A2b). Marked transitions in $C_i$:$C_a$ result from transitions between RuBP-regeneration limited photosynthesis (morning and evening) and Rubisco-limited photosynthesis (afternoon; see Method S1). Schemes that simulate
higher morning and evening $C_i:C_a$ often transition half an hour later than Medlyn—particularly in the morning— and those that simulate lower $C_i:C_a$ transition about half an hour earlier—particularly in the evening.

In the early morning and evening, relative differences in $g_s$ and $C_i$ make most optimization schemes more prof- ligate water spenders than Medlyn (up to >20% increase in $E$; Figure A1e) for negligible carbon gains (<3% increase in $A_n$; Figure A1d), and vice versa in the afternoon. Having a somewhat riskier strategy, WUE-$f_{\psi_l}$ (both under “wet” and “stressed” conditions) and LeastCost (under “stressed” conditions) can transpire more in the afternoon (up to c. 10%) for modest carbon gains (up to <4%) that risk being offset by photorespiration when investments in $E$ are insufficient. Note, neither WUE-$f_{\psi_l}$ nor LeastCost explicitly account for hydraulics in their optimization criteria formulation. CAP and MES stand out, always simulating lower $A_n$ than Medlyn, because their hydraulic controls directly downregulate $V_{cmax}$ and $J_{max}$ (Equation 21) or mesophyll conductance (Equa- tion 22). Consequently, both models must maintain high $C_i$ (Figure A1b) to prevent stark reductions in $A_n$.

Simulation of $\Psi_l$ (Figure A1c) is where schemes most diverge. LeastCost and Eller simulate lower $\Psi_l$ than all other schemes, even beyond the $P_{88}$. These models' calibrated $k_{max}$ is smaller than the ProfitMax's (Table S6 in Supporting Information S1), so LeastCost and Eller have lower $\Psi_l$ at the same magnitude $g_s$ than other models, pointing to almost no hydraulic controls. Conversely, ProfitMax2, CAP, and MES exert tight hydraulic controls, relying on higher calibrated $k_{max}$ than ProfitMax, and so, yielding higher $\Psi_l$. Both CGain and SOX$_{opt}$ behave simi- larly to ProfitMax, with ever so marginally lower $\Psi_l$. Under “stressed” conditions, CAP and MES (Figure A2c) operate in a contrasting mode of optimization compared to the other models, where $k_{max}$ is over three times that of the “wet” calibrations, so $\Psi_l$ barely drops throughout the day. Finally, CMax simulates stomatal closure in the afternoon, accompanied by a reversal of the diurnal leaf water potential drop (most visible in Figure A2c), when $T_a$ and $D_a$ rise.

## Data Availability Statement

All model, analysis code, and data are freely available from https://github.com/ManonSabot/One_gs_model_to- rule_them_all and https://doi.org/10.5281/zenodo.5932661 (Sabot, 2022).

## Acknowledgments

M. E. B. Sabot, M. G. De Kauwe, A. J. Pitman acknowledge support from the Australian Research Council Centre of Excellence for Climate Extremes (CE170100023) and M. E. B. Sabot was also supported by the UNSW Scien- tia PhD Scholarship Scheme. M. G. De Kauwe and A. J. Pitman acknowl- edge support from the ARC Discovery Grant (DP190101823) and M. G. De Kauwe also acknowledges the NSW Research Attraction and Acceleration Program. B. E. Medlyn acknowledges support from ARC Laureate Fellowship FL190100003. A. Rogers and S. P. Serbin were supported by the Next-Generation Ecosystem Experiments Tropics project that is supported by the Office of Biolog- ical and Environmental Research in the Department of Energy, Office of Science, and through the United States Department of Energy contract No. DE-SC0012704 to Brookhaven National Laboratory. J. Wu was in part supported by the Innovation and Technology Fund (funding support to State Key Laboratories in Hong Kong of Agrobiotechnology) of the HKSAR, China. The authors are grateful to John Finnegan and Ying-Ping Wang for helpful discussion of leaf-level boundary layer feedbacks. The authors also thank anon- ymous reviewers for their constructive comments.

## References

Ahlstrom, A., Raupach, M. R., Schurgers, G., Smith, B., Arneth, A., Jung, M., et al. (2015). The dominant role of semi-arid ecosystems in the trend and variability of the land CO$_2$ sink. *Science*, **348**(6237), 895–899. https://doi.org/10.1126/science.aaa1668

Anderegg, W. R. L., Wolf, A., Arango-Velez, A., Choat, B., Chmura, D. J., Jansen, S., et al. (2018). Woody plants optimise stomatal behaviour relative to hydraulic risk. *Ecology Letters*, **21**(7), 968–977. https://doi.org/10.1111/ele.12962

Bahar, N. H. A., Hayes, L., Scafaro, A. P., Atkin, O. K., & Evans, J. R. (2018). Mesophyll conductance does not contribute to greater photosyn- thetic rate per unit nitrogen in temperate compared with tropical evergreen wet-forest tree leaves. *New Phytologist*, **218**(2), 492–505. https:// doi.org/10.1111/nph.15031

Ball, J. T., Woodrow, I. E., & Berry, J. A. (1987). A model predicting stomatal conductance and its contribution to the control of photosynthesis under different environmental conditions. In J. Biggins (Ed.), *Progress in photosynthesis research* (pp. 221–224). Springer Netherlands. https:// doi.org/10.1007/978-94-017-0519-6_48

Bassiouni, M., & Vico, G. (2021). Parsimony vs predictive and functional performance of three stomatal optimization principles in a big-leaf framework. *New Phytologist*, **231**(2), 586–600. https://doi.org/10.1111/nph.17392

Bourne, A. E., Creek, D., Peters, J. M. R., Ellsworth, D. S., & Choat, B. (2017). Species climate range influences hydraulic and stomatal traits in *Eucalyptus* species. *Annals of Botany*, **120**(1), 123–133. https://doi.org/10.1093/aob/mcx020

Brodribb, T. J., Holbrook, N. M., Edwards, E. J., & Gutierrez, M. V. (2003). Relations between stomatal closure, leaf turgor and xylem vulnera- bility in eight tropical dry forest trees. *Plant, Cell & Environment*, **26**(3), 443–450. https://doi.org/10.1046/j.1365-3040.2003.00975.x

Buckley, T. N. (2019). How do stomata respond to water status? *New Phytologist*, **224**(1), 21–36. https://doi.org/10.1111/nph.15899

Choat, B., Ball, M., Luly, J., & Holtum, J. (2003). Pit membrane porosity and water stress-induced cavitation in four co-existing dry rainforest tree species. *Plant Physiology*, **131**(1), 41–48. https://doi.org/10.1104/pp.014100

Choat, B., Ball, M. C., Luly, J. G., Donnelly, C. F., & Holtum, J. A. M. (2006). Seasonal patterns of leaf gas exchange and water relations in dry rain forest trees of contrasting leaf phenology. *Tree Physiology*, **26**(5), 657–664. https://doi.org/10.1093/treephys/26.5.657

Choat, B., Jansen, S., Brodribb, T. J., Cochard, H., Delzon, S., Bhaskar, R., et al. (2012). Global convergence in the vulnerability of forests to drought. *Nature*, **491**, 752–755. https://doi.org/10.1038/nature11688

Clark, D. B., Mercado, L. M., Sitch, S., Jones, C. D., Gedney, N., Best, M. J., et al. (2011). The Joint UK Land Environment Simulator (JULES), model description—Part 2: Carbon fluxes and vegetation dynamics. *Geoscientific Model Development*, **4**(3), 701–722. https://doi.org/10.5194/ gmd-4-701-2011

Cochard, H., Badel, E., Herbette, S., Delzon, S., Choat, B., & Jansen, S. (2013). Methods for measuring plant vulnerability to cavitation: A critical review. *Journal of Experimental Botany*, **64**(15), 4779–4791. https://doi.org/10.1093/jxb/ert193

Cowan, I. R. (1982). Regulation of water use in relation to carbon gain in higher plants. In O. L. Lange, P. S. Nobel, C. B. Osmond, & H. Ziegler (Eds.), *Physiological plant ecology II* (pp. 589–613). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-68150-9_18

Cowan, I. R., & Farquhar, G. D. (1977). Stomatal function in relation to leaf metabolism and environment. *Symposia of the Society for Experi- mental Biology*, **31**, 471–505.
Crous, K. Y., Quentin, A. G., Lin, Y.-S., Medlyn, B. E., Williams, D. G., Barton, C. V. M., & Ellsworth, D. S. (2013). Photosynthesis of temperate *Eucalyptus globulus* trees outside their native range has limited adjustment to elevated CO₂ and climate warming. *Global Change Biology*, *19*(12), 3790–3807. https://doi.org/10.1111/gcb.12314

Damour, G., Simonneau, T., Cochard, H., & Urban, L. (2010). An overview of models of stomatal conductance at the leaf level. *Plant, Cell & Environment*, *33*(9), 1419–1438. https://doi.org/10.1111/j.1365-3040.2010.02181.x

de Boer, H. J., Drake, P. L., Wendt, E., Price, C. A., Schulze, E.-D., Turner, N. C., et al. (2016). Apparent overinvestment in leaf venation relaxes leaf morphological constraints on photosynthesis in arid habitats. *Plant Physiology*, *172*(4), 2286–2299. https://doi.org/10.1104/pp.16.01313

De Cáceres, M., Mencuccini, M., Martin-StPaul, N., Limousin, J.-M., Coll, L., Poyatos, R., et al. (2021). Unravelling the effect of species mixing on water use and drought stress in Mediterranean forests: A modelling approach. *Agricultural and Forest Meteorology*, *296*, 108233. https://doi.org/10.1016/j.agrformet.2020.108233

De Kauwe, M. G., Kala, J., Lin, Y.-S., Pitman, A. J., Medlyn, B. E., Duursma, R. A., et al. (2015). A test of an optimal stomatal conductance scheme within the CABLE land surface model. *Geoscientific Model Development*, *8*(2), 431–452. https://doi.org/10.5194/gmd-8-431-2015

De Kauwe, M. G., Lin, Y.-S., Wright, I. J., Medlyn, B. E., Crous, K. Y., Ellsworth, D. S., et al. (2016). A test of the 'one-point method' for estimating maximum carboxylation capacity from field-measured, light-saturated photosynthesis. *New Phytologist*, *210*(3), 1130–1144. https://doi.org/10.1111/nph.13815

De Kauwe, M. G., Medlyn, B. E., Ukkola, A. M., Mu, M., Sabot, M. E. B., Pitman, A. J., et al. (2020). Identifying areas at risk of drought-induced tree mortality across South-Eastern Australia. *Global Change Biology*, *26*(10), 5716–5733. https://doi.org/10.1111/gcb.15215

De Kauwe, M. G., Medlyn, B. E., Zaehle, S., Walker, A. P., Dietze, M. C., Hickler, T., et al. (2013). Forest water use and water use efficiency at elevated CO₂: A model-data intercomparison at two contrasting temperate forest FACE sites. *Global Change Biology*, *19*(6), 1759–1779. https://doi.org/10.1111/gcb.12164

Dewar, R., Mauranen, A., Mäkelä, A., Hölttä, T., Medlyn, B., & Vesala, T. (2017). New insights into the covariation of stomatal, mesophyll and hydraulic conductances from optimization models incorporating nonstomatal limitations to photosynthesis. *New Phytologist*, *217*, 571–585. https://doi.org/10.1111/nph.14848

Donat, M. G., Pitman, A. J., & Angélil, O. (2018). Understanding and reducing future uncertainty in midlatitude daily heat extremes via land surface feedback constraints. *Geophysical Research Letters*, *45*(19), 10627–10636. https://doi.org/10.1029/2018GL079128

Duursma, R. A., & Medlyn, B. E. (2012). MAESPA: A model to study interactions between water limitation, environmental drivers and vegetation function at tree and stand levels, with an example application to [CO₂] × drought interactions. *Geoscientific Model Development*, *5*(4), 919–940. https://doi.org/10.5194/gmd-5-919-2012

Eller, C. B., Rowland, L., Mencuccini, M., Rosas, T., Williams, K., Harper, A., et al. (2020). Stomatal optimization based on xylem hydraulics (SOX) improves land surface model simulation of vegetation responses to climate. *New Phytologist*, *226*, 1622–1637. https://doi.org/10.1111/nph.16419

Eller, C. B., Rowland, L., Oliveira, R. S., Bittencourt, P. R. L., Barros, F. V., da Costa, A. C. L., et al. (2018). Modelling tropical forest responses to drought and El Niño with a stomatal optimization model based on xylem hydraulics. *Philosophical Transactions of the Royal Society B: Biological Sciences*, *373*(1760), 20170315. https://doi.org/10.1098/rstb.2017.0315

Ellsworth, D. S. (2000). Seasonal CO₂ assimilation and stomatal limitations in a *Pinus taeda* canopy. *Tree Physiology*, *20*(7), 435–445. https://doi.org/10.1093/treephys/20.7.435

Farquhar, G. D., Caemmerer, S. von, & Berry, J. A. (1980). A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species. *Planta*, *149*(1), 78–90. https://doi.org/10.1007/BF00386231

Fisher, R. A., & Koven, C. D. (2020). Perspectives on the future of land surface models and the challenges of representing complex terrestrial systems. *Journal of Advances in Modeling Earth Systems*, *12*(4). https://doi.org/10.1029/2018MS001453

Franks, P. J., Berry, J. A., Lombardozzi, D. L., & Bonan, G. B. (2017). Stomatal function across temporal and spatial scales: Deep-time trends, land-atmosphere coupling and global models. *Plant Physiology*, *174*(2), 583–602. https://doi.org/10.1104/pp.17.00287

Good, S. P., Noone, D., & Bowen, G. (2015). Hydrologic connectivity constrains partitioning of global terrestrial water fluxes. *Science*, *349*(6244), 175–177. https://doi.org/10.1126/science.aaa5931

Hawkins, E., Robson, J., Sutton, R., Smith, D., & Keenlyside, N. (2011). Evaluating the potential for statistical decadal predictions of sea surface temperatures with a perfect model approach. *Climate Dynamics*, *37*(11–12), 2495–2509. https://doi.org/10.1007/s00382-011-1023-3

Herman, J., & Usher, W. (2017). SALib: An open-source Python library for sensitivity analysis. *The Journal of Open Source Software*, *2*(9), 97. https://doi.org/10.21105/joss.00097

Héroult, A., Lin, Y.-S., Bourne, A., Medlyn, B. E., & Ellsworth, D. S. (2013). Optimal stomatal conductance in relation to photosynthesis in climatically contrasting *Eucalyptus* species under drought: Stomatal responses of eucalyptus under drought. *Plant, Cell & Environment*, *36*(2), 262–274. https://doi.org/10.1111/j.1365-3040.2012.02570.x

Hetherington, A. M., & Woodward, F. I. (2003). The role of stomata in sensing and driving environmental change. *Nature*, *424*(6951), 901–908. https://doi.org/10.1038/nature01843

Huang, C.-W., Domec, J.-C., Palmroth, S., Pockman, W. T., Litvak, M. E., & Katul, G. G. (2018). Transport in a coordinated soil-root-xylem-phloem leaf system. *Advances in Water Resources*, *119*, 1–16. https://doi.org/10.1016/j.advwatres.2018.06.002

Hudson, P. J., Limousin, J. M., Krofcheck, D. J., Boutz, A. L., Pangle, R. E., Gehres, N., et al. (2018). Impacts of long-term precipitation manipulation on hydraulic architecture and xylem anatomy of piñon and juniper in Southwest USA: Xylem acclimation in piñon-juniper woodland. *Plant, Cell & Environment*, *41*(2), 421–435. https://doi.org/10.1111/pce.13109

Humphrey, V., Zscheischler, J., Ciais, P., Gudmundsson, L., Sitch, S., & Seneviratne, S. I. (2018). Sensitivity of atmospheric CO₂ growth rate to observed changes in terrestrial water storage. *Nature*, *560*(7720), 628–631. https://doi.org/10.1038/s41586-018-0424-4

Jarvis, P. G. (1976). The interpretation of the variations in leaf water potential and stomatal conductance found in canopies in the field. *Philosophical Transactions of the Royal Society B: Biological Sciences*, *273*(927), 593–610. https://doi.org/10.1098/rstb.1976.0035

Jones, H. G. (1998). Stomatal control of photosynthesis and transpiration. *Journal of Experimental Botany*, *49*(Special), 387–398. https://doi.org/10.1093/jxb/49.Special_Issue.387

Jung, M., Reichstein, M., Schwalm, C. R., Huntingford, C., Sitch, S., Ahlström, A., et al. (2017). Compensatory water effects link yearly global land CO₂ sink changes to temperature. *Nature*, *541*(7638), 516–520. https://doi.org/10.1038/nature20780

Kala, J., De Kauwe, M. G., Pitman, A. J., Lorenz, R., Medlyn, B. E., Wang, Y.-P., et al. (2015). Implementation of an optimal stomatal conductance scheme in the Australian Community Climate Earth Systems Simulator (ACCESS1.3b). *Geoscientific Model Development*, *8*(12), 3877–3889. https://doi.org/10.5194/gmd-8-3877-2015

Keenan, T., Sabate, S., & Gracia, C. (2010). The importance of mesophyll conductance in regulating forest ecosystem productivity during drought periods. *Global Change Biology*, *16*(3), 1019–1034. https://doi.org/10.1111/j.1365-2486.2009.02017.x
Kennedy, D., Swenson, S., Oleson, K. W., Lawrence, D. M., Fisher, R., Lola da Costa, A. C., & Gentine, P. (2019). Implementing plant hydraulics in the community land model, version 5. *Journal of Advances in Modeling Earth Systems*, 11(2), 485–513. https://doi.org/10.1029/2018MS001500

Knauer, J., Zaehle, S., De Kauwe, M. G., Bahar, N. H. A., Evans, J. R., Medlyn, B. E., et al. (2019). Effects of mesophyll conductance on vegetation responses to elevated CO₂ concentrations in a land surface model. *Global Change Biology*, 25(5), 1820–1838. https://doi.org/10.1111/gcb.14604

Kumarathunge, D. P., Medlyn, B. E., Drake, J. E., Tjoelker, M. G., Aspinwall, M. J., Battaglia, M., et al. (2019). Acclimation and adaptation components of the temperature dependence of plant photosynthesis at the global scale. *New Phytologist*, 222(2), 768–784. https://doi.org/10.1111/nph.15668

Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G., et al. (2019). The Community Land Model version 5: Description of new features, benchmarking, and impact of forcing uncertainty. *Journal of Advances in Modeling Earth Systems*, 11(12), 4245–4287. https://doi.org/10.1029/2018MS001583

Lawson, T., von Caemmerer, S., & Baroli, I. (2010). Photosynthesis and stomatal behaviour. In U. E. Lüttge, W. Beyschlag, B. Büdel, & D. Francis (Eds.), *Progress in Botany* (Vol. 72, pp. 265–304). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-13145-5_11

Leplastrier, M. (2002). Exploring the relationship between complexity and performance in a land surface model using the multicriteria method. *Journal of Geophysical Research*, 107(D20), 4443. https://doi.org/10.1029/2001JD000931

Leuning, R. (1995). A critical appraisal of a combined stomatal-photosynthesis model for C3 plants. *Plant, Cell and Environment*, 18(4), 339–355. https://doi.org/10.1111/j.1365-3040.1995.tb00370.x

Li, X., Gentine, P., Lin, C., Zhou, S., Sun, Z., Zheng, Y., et al. (2019). A simple and objective method to partition evapotranspiration into transpiration and evaporation at eddy-covariance sites. *Agricultural and Forest Meteorology*, 265, 171–182. https://doi.org/10.1016/j.agrformet.2018.11.017

Limousin, J.-M., Bickford, C. P., Dickman, L. T., Pangle, R. E., Hudson, P. J., Boutz, A. L., et al. (2013). Regulation and acclimation of leaf gas exchange in a piñon-juniper woodland exposed to three different precipitation regimes: Rainfall manipulation in piñon-juniper woodland. *Plant, Cell & Environment*, 36(10), 1812–1825. https://doi.org/10.1111/pce.12089

Lin, Y.-S., Medlyn, B. E., Duursma, R. A., Prentice, I. C., Wang, H., Baig, S., et al. (2015). Optimal stomatal behaviour around the world. *Nature Climate Change*, 5(5), 459–464. https://doi.org/10.1038/nclimate2550

Lopez, O. R., Kursar, T. A., Cochard, H., & Tyree, M. T. (2005). Interspecific variation in xylem vulnerability to cavitation among tropical tree and shrub species. *Tree Physiology*, 25(12), 1553–1562. https://doi.org/10.1093/treephys/25.12.1553

Lu, Y., Duursma, R. A., Farrior, C. E., Medlyn, B. E., & Feng, X. (2020). Optimal stomatal drought response shaped by competition for water and hydraulic risk can explain plant trait covariation. *New Phytologist*, 225(3), 1206–1217. https://doi.org/10.1111/nph.16207

Maherali, H., Moura, C. F., Caldeira, M. C., Willson, C. J., & Jackson, R. B. (2006). Functional coordination between leaf gas exchange and vulnerability to xylem cavitation in temperate forest trees. *Plant, Cell and Environment*, 29(4), 571–583. https://doi.org/10.1111/j.1365-3040.2005.01433.x

Mäkelä, A. (1996). Optimal control of gas exchange during drought: Theoretical analysis. *Annals of Botany*, 77(5), 461–468. https://doi.org/10.1006/anbo.1996.0056

Mankin, J. S., Seager, R., Smerdon, J. E., Cook, B. I., & Williams, A. P. (2019). Mid-latitude freshwater availability reduced by projected vegetation responses to climate change. *Nature Geoscience*, 12(12), 983–988. https://doi.org/10.1038/s41561-019-0480-x

Manzoni, S., Vico, G., Katul, G., Fay, P. A., Polley, W., Palmroth, S., & Porporato, A. (2011). Optimizing stomatal conductance for maximum carbon gain under water stress: A meta-analysis across plant functional types and climates: Optimal leaf gas exchange under water stress. *Functional Ecology*, 25(3), 456–467. https://doi.org/10.1111/j.1365-2435.2010.01822.x

Manzoni, S., Vico, G., Katul, G., Palmroth, S., Jackson, R. B., & Porporato, A. (2013). Hydraulic limits on maximum plant transpiration and the emergence of the safety-efficiency trade-off. *New Phytologist*, 198(1), 169–178. https://doi.org/10.1111/nph.12126

Martínez-de la Torre, A., Blyth, E., & Robinson, E. (2019). Evaluation of drydown processes in global land surface and hydrological models using flux tower evapotranspiration. *Water*, 11(2), 356. https://doi.org/10.3390/w11020356

Martínez-Vilalta, J., & Garcia-Forner, N. (2017). Water potential regulation, stomatal behaviour and hydraulic transport under drought: Deconstructing the iso/anisohydric concept: Deconstructing the iso/anisohydric concept. *Plant, Cell & Environment*, 40(6), 962–976. https://doi.org/10.1111/pce.12846

Martin-StPaul, N., Delzon, S., & Cochard, H. (2017a). Plant resistance to drought depends on timely stomatal closure. *Ecology Letters*, 20(11), 1437–1447. https://doi.org/10.1111/ele.12851

Martin-StPaul, N., Delzon, S., & Cochard, H. (2017b). SurEau database: A database of hydraulic and stomatal traits for modelling drought resistance in plants [Data set]. Zenodo. https://doi.org/10.5281/zenodo.854700

Martin-StPaul, N., Longepierre, D., Huc, R., Delzon, S., Burlett, R., Joffre, R., et al. (2014). How reliable are methods to assess xylem vulnerability to cavitation? The issue of “open vessel” artifact in oaks. *Tree Physiology*, 34(8), 894–905. https://doi.org/10.1093/treephys/tpu059

Martin-StPaul, N. K., Limousin, J.-M., Rodríguez-Calcerrada, J., Ruffault, J., Rambal, S., Letts, M. G., & Misson, L. (2012). Photosynthetic sensitivity to drought varies among populations of Quercus ilex along a rainfall gradient. *Functional Plant Biology*, 39(1), 25. https://doi.org/10.1071/FP11090

Marzban, C. (2013). Variance-based sensitivity analysis: An illustration on the Lorenz'63 model. *Monthly Weather Review*, 141(11), 4069–4079. https://doi.org/10.1175/MWR-D-13-00032.1

Medlyn, B. E., De Kauwe, M. G., Zaehle, S., Walker, A. P., Duursma, R. A., Luus, K., et al. (2016). Using models to guide field experiments: A priori predictions for the CO₂ response of a nutrient- and water-limited native eucalypt woodland. *Global Change Biology*, 22(8), 2834–2851. https://doi.org/10.1111/gcb.13268

Medlyn, B. E., Duursma, R. A., Eamus, D., Ellsworth, D. S., Prentice, I. C., Barton, C. V. M., et al. (2011). Reconciling the optimal and empirical approaches to modelling stomatal conductance. *Global Change Biology*, 17(6), 2134–2144. https://doi.org/10.1111/j.1365-2486.2010.02375.x

Medlyn, B. E., Zaehle, S., De Kauwe, M. G., Walker, A. P., Dietze, M. C., Hanson, P. J., et al. (2015). Using ecosystem experiments to improve vegetation models. *Nature Climate Change*, 5(6), 528–534. https://doi.org/10.1038/nclimate2621

Meinzer, F. C., Woodruff, D. R., Domec, J.-C., Goldstein, G., Campanello, P. I., Gatti, M. G., & Villalobos-Vega, R. (2008). Coordination of leaf and stem water transport properties in tropical forest trees. *Oecologia*, 156(1), 31–41. https://doi.org/10.1007/s00442-008-0974-5

Mencuccini, M., Manzoni, S., & Christoffersen, B. (2019). Modelling water fluxes in plants: From tissues to biosphere. *New Phytologist*, 222(3), 1207–1222. https://doi.org/10.1111/nph.15681

Miner, G. L., Bauerle, W. L., & Baldocchi, D. D. (2017). Estimating the sensitivity of stomatal conductance to photosynthesis: A review. *Plant, Cell & Environment*, 40(7), 1214–1238. https://doi.org/10.1111/pce.12871

Miralles, D. G., Gentine, P., Seneviratne, S. I., & Teuling, A. J. (2018). Land-atmospheric feedbacks during droughts and heatwaves: State of the science and current challenges. *Annals of the New York Academy of Sciences*, 1436(1), 19–35. https://doi.org/10.1111/nyas.13912
Miralles, D. G., Teuling, A. J., van Heerwaarden, C. C., & Vilà-Guerau de Arellano, J. (2014). Mega-heatwave temperatures due to combined soil desiccation and atmospheric heat accumulation. *Nature Geoscience*, 7(5), 345–349. https://doi.org/10.1038/ngeo2141

Mitchell, P. J. (2009). *From conduits to communities: Plant water use strategies and evapotranspiration in a semi-arid ecosystem in south-western Australia*. The University of Western Australia. Retrieved from https://research-repository.uwa.edu.au/en/publications/from-conduits-to-communities-plant-water-use-strategies-and-evapo-

Mitchell, P. J., Veneklaas, E., Lambers, H., & Burgess, S. S. O. (2009). Partitioning of evapotranspiration in a semi-arid eucalypt woodland in south-western Australia. *Agricultural and Forest Meteorology*, 149(1), 25–37. https://doi.org/10.1016/j.agrformet.2008.07.008

Mrad, A., Sevanto, S., Domec, J.-C., Liu, Y., Nakad, M., & Katul, G. (2019). A dynamic optimality principle for water use strategies explains isohydric to anisohydric plant responses to drought. *Frontiers in Forests and Global Change*, 2. https://doi.org/10.3389/ffgc.2019.00049

Neufeld, H. S., Grantz, D. A., Meinzer, F. C., Goldstein, G., Crisosto, G. M., & Crisosto, C. (1992). Genotypic variability in vulnerability of leaf xylem to cavitation in water-stressed and well-irrigated sugarcane. *Plant Physiology*, 100(2), 1020–1028. https://doi.org/10.1104/ pp.100.2.1020

Newville, M., Otten, R., Nelson, A., Ingargiola, A., Stensitzki, T., Allan, D., et al. (2019). Imfit/lmfit-py 1.0.0. Zenodo. https://doi.org/10.5281/ ZENODO.3588521

Novick, K. A., Miniat, C. F., & Vose, J. M. (2016). Drought limitations to leaf-level gas exchange: Results from a model linking stomatal optimi- zation and cohesion-tension theory. *Plant, Cell & Environment*, 39(3), 583–596. https://doi.org/10.1111/pce.12657

Ogle, K., Barber, J. J., Willson, C., & Thompson, B. (2009). Hierarchical statistical modeling of xylem vulnerability to cavitation. *New Phytolo- gist*, 182(2), 541–554. https://doi.org/10.1111/j.1469-8137.2008.02760.x

Oren, R., Sperry, J. S., Katul, G. G., Pataki, D. E., Ewers, B. E., Phillips, N., & Schäfer, K. V. R. (1999). Survey and synthesis of intra- and inter- specific variation in stomatal sensitivity to vapour pressure deficit: Intra- and interspecific variation in stomatal sensitivity to vapour pressure deficit. *Plant, Cell & Environment*, 22(12), 1515–1526. https://doi.org/10.1046/j.1365-3040.1999.00513.x

Potkay, A., Trugman, A. T., Wang, Y., Venturas, M. D., Anderegg, W. R. L., Mattos, C. R. C., & Fan, Y. (2021). Coupled whole-tree optimality and xylem-hydraulics explain dynamic biomass partitioning. *New Phytologist*, 230, 2226–2245. https://doi.org/10.1111/nph.17242

Poyatos, R., Granda, V., Flo, V., Adams, M. A., Adorján, B., Aguadé, D., et al. (2021). Global transpiration data from sap flow measurements: The SAPFLUXNET database. *Earth System Science Data*, 13(6), 2607–2649. https://doi.org/10.5194/essd-13-2607-2021

Prentice, I. C., Dong, N., Gleason, S. M., Maire, V., & Wright, I. J. (2014). Balancing the costs of carbon gain and water transport: Testing a new theoretical framework for plant functional ecology. *Ecology Letters*, 17(1), 82–91. https://doi.org/10.1111/ele.12211

Raupach, M. R., & Finnigan, J. J. (1988). “Single-layer models of evaporation from plant canopies are incorrect but useful, whereas multilayer models are correct but useless”: Discuss. *Australian Journal of Plant Physiology*, 15(6), 705–716. https://doi.org/10.1071/pp9880705

Reu, B., Zaehle, S., Bohn, K., Pavlick, R., Schmidtlein, S., Williams, J. W., & Kleidon, A. (2014). Future no-analogue vegetation produced by no-analogue combinations of temperature and insolation: No-analogue vegetation, temperature and insolation. *Global Ecology and Biogeog- raphy*, 23(2), 156–167. https://doi.org/10.1111/geb.12110

Rogers, A., Medlyn, B. E., Dukes, J. S., Bonan, G., von Caemmerer, S., Dietze, M. C., et al. (2017). A roadmap for improving the representation of photosynthesis in Earth system models. *New Phytologist*, 213(1), 22–42. https://doi.org/10.1111/nph.14283

Rogers, A., Serbin, S., Ely, K., Wu, J., Wolfe, B., Dickman, T., et al. (2019). CO2 response (ACI) gas exchange, calculated Vcmax & Jmax parameters, February 2016–May 2016, PA-SLZ, PA-PNM: Panama [Data set]. Next-Generation Ecosystem Experiments Tropics; Brookhaven National Laboratory; Los Alamos National Laboratory; STRI. https://doi.org/10.15486/NGT/1411867

Sabot, M. E. B. (2022). One GS_model_to_rule_them_all (v1.0). Zenodo. https://doi.org/10.5281/zenodo.5932661

Sabot, M. E. B., De Kauwe, M. G., Pitman, A. J., Medlyn, B. E., Verhoef, A., Ukkola, A. M., & Abramowitz, G. (2020). Plant profit maximization improves predictions of European forest responses to drought. *New Phytologist*, 226, 1638–1655. https://doi.org/10.1111/nph.16376

Sala, A., & Tenhunen, J. D. (1996). Simulations of canopy net photosynthesis and transpiration in Quercus ilex L. under the influence of seasonal drought. *Agricultural and Forest Meteorology*, 78(3–4), 203–222. https://doi.org/10.1016/0168-1923(95)02250-3

Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., & Tarantola, S. (2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. *Computer Physics Communications*, 181(2), 259–270. https://doi.org/10.1016/j.cpc.2009.09.018

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., et al. (2007). *Global sensitivity analysis. The primer*. John Wiley & Sons, Ltd. https://doi.org/10.1002/9780470725184

Schlesinger, W. H., & Jasechko, S. (2014). Transpiration in the global water cycle. *Agricultural and Forest Meteorology*, 189–190, 115–117. https://doi.org/10.1016/j.agrformet.2014.01.011

Sellers, P. J., Dickinson, R. E., Randall, D. A., Betts, A., Hall, F. G., Berry, J. A., et al. (1997). Modeling the exchanges of energy, water, and carbon between continents and the atmosphere. *Science*, 275(5299), 502–509. https://doi.org/10.1126/science.275.5299.502

Sobol', I. M. (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. *Mathematics and Computers in Simulation*, 55(1–3), 271–280. https://doi.org/10.1016/S0378-4754(00)00270-6

Sperry, J. S., & Love, D. M. (2015). What plant hydraulics can tell us about responses to climate-change droughts. *New Phytologist*, 207(1), 14–27. https://doi.org/10.1111/nph.13354

Sperry, J. S., Venturas, M. D., Anderegg, W. R. L., Mencuccini, M., Mackay, D. S., Wang, Y., & Love, D. M. (2017). Predicting stomatal responses to the environment from the optimization of photosynthetic gain and hydraulic cost: A stomatal optimization model. *Plant, Cell & Environment*, 40(6), 816–830. https://doi.org/10.1111/pce.12852

Sperry, J. S., Venturas, M. D., Todd, H. N., Trugman, A. T., Anderegg, W. R. L., Wang, Y., & Tai, X. (2019). The impact of rising CO₂ and acclimation on the response of US forests to global warming. *Proceedings of the National Academy of Sciences*, 116(51), 25734–25744. https:// doi.org/10.1073/pnas.1913072116

Stoy, P. C., El-Madany, T. S., Fisher, J. B., Gentine, P., Gerken, T., Good, S. P., et al. (2019). Reviews and syntheses: Turning the challenges of partitioning ecosystem evaporation and transpiration into opportunities. *Biogeosciences*, 16(19), 3747–3775. https://doi.org/10.5194/ bg-16-3747-2019

Sun, Y., Gu, L., Dickinson, R. E., Pallardy, S. G., Baker, J., Cao, Y., et al. (2014). Asymmetrical effects of mesophyll conductance on fundamental photosynthetic parameters and their relationships estimated from leaf gas exchange measurements: Asymmetrical mesophyll conductance effects. *Plant, Cell & Environment*, 37(4), 978–994. https://doi.org/10.1111/pce.12213

Swann, A. L. S., Hoffman, F. M., Koven, C. D., & Randerson, J. T. (2016). Plant responses to increasing CO₂ reduce estimates of climate impacts on drought severity. *Proceedings of the National Academy of Sciences*, 113(36), 10019–10024. https://doi.org/10.1073/pnas.1604581113

Tardieu, F., & Simonneau, T. (1998). Variability among species of stomatal control under fluctuating soil water status and evaporative demand: Modelling isohydric and anisohydric behaviours. *Journal of Experimental Botany*, 49, 419–432. https://doi.org/10.1093/jxb/49.special issue.419

Trancoso, R., Larsen, J. R., McVicar, T. R., Phinn, S. R., & McAlpine, C. A. (2017). CO₂-vegetation feedbacks and other climate changes impli- cated in reducing base flow. *Geophysical Research Letters*, 44(5), 2310–2318. https://doi.org/10.1002/2017GL072759
Trugman, A. T., Anderegg, L. D. L., Sperry, J. S., Wang, Y., Venturas, M., & Anderegg, W. R. L. (2019). Leveraging plant hydraulics to yield predictive and dynamic plant leaf allocation in vegetation models with climate change. *Global Change Biology*, 25(12), 4008–4021. https://doi.org/10.1111/gcb.14814

Tuzet, A., Perrier, A., & Leuning, R. (2003). A coupled model of stomatal conductance, photosynthesis and transpiration. *Plant, Cell and Environment*, 26(7), 1097–1116. https://doi.org/10.1046/j.1365-3040.2003.01035.x

Ukkola, A. M., De Kauwe, M. G., Pitman, A. J., Best, M. J., Abramowitz, G., Haverd, V., et al. (2016). Land surface models systematically overestimate the intensity, duration and magnitude of seasonal-scale evaporative droughts. *Environmental Research Letters*, 11(10), 104012. https://doi.org/10.1088/1748-9326/11/10/104012

Ukkola, A. M., Prentice, I. C., Keenan, T. F., van Dijk, A. I. J. M., Viney, N. R., Myneni, R. B., & Bi, J. (2016). Reduced streamflow in water-stressed climates consistent with CO₂ effects on vegetation. *Nature Climate Change*, 6(1), 75–78. https://doi.org/10.1038/nclimate2831

Van Looy, K., Bouma, J., Herbst, M., Koestel, J., Minasny, B., Mishra, U., et al. (2017). Pedotransfer functions in Earth system science: Challenges and perspectives. *Reviews of Geophysics*, 55(4), 1199–1256. https://doi.org/10.1002/2017RG000581

Venturas, M. D., Sperry, J. S., Love, D. M., Frehner, E. H., Allred, M. G., Wang, Y., & Anderegg, W. R. L. (2018). A stomatal control model based on optimization of carbon gain versus hydraulic risk predicts aspen sapling responses to drought. *New Phytologist*, 220(3), 836–850. https://doi.org/10.1111/nph.15333

Venturas, M. D., Todd, H. N., Trugman, A. T., & Anderegg, W. R. L. (2020). Understanding and predicting forest mortality in the western United States using long-term forest inventory data and modeled hydraulic damage. *New Phytologist*, 230, 1896–1910, 17043. https://doi.org/10.1111/nph.17043

Verhoef, A., & Egea, G. (2014). Modeling plant transpiration under limited soil water: Comparison of different plant and soil hydraulic parameterizations and preliminary implications for their use in land surface models. *Agricultural and Forest Meteorology*, 191, 22–32. https://doi.org/10.1016/j.agrformet.2014.02.009

Vialet-Chabrand, S. R. M., Matthews, J. S. A., McAusland, L., Blatt, M. R., Griffiths, H., & Lawson, T. (2017). Temporal dynamics of stomatal behavior: Modeling and implications for photosynthesis and water use. *Plant Physiology*, 174(2), 603–613. https://doi.org/10.1104/pp.17.00125

Wang, Y., Sperry, J. S., Anderegg, W. R. L., Venturas, M. D., & Trugman, A. T. (2020). A theoretical and empirical assessment of stomatal optimization modeling. *New Phytologist*, 227, 311–325. https://doi.org/10.1111/nph.16572

Wang, Y., Sperry, J. S., Venturas, M. D., Trugman, A. T., Love, D. M., & Anderegg, W. R. L. (2019). The stomatal response to rising CO₂ concentration and drought is predicted by a hydraulic trait-based optimization model. *Tree Physiology*, 39(8), 1416–1427. https://doi.org/10.1093/treephys/tpz038

Wei, Z., Yoshimura, K., Wang, L., Miralles, D. G., Jasechko, S., & Lee, X. (2017). Revisiting the contribution of transpiration to global terrestrial evapotranspiration. *Geophysical Research Letters*, 44(6), 2792–2801. https://doi.org/10.1002/2016GL072235

Wolf, A., Anderegg, W. R. L., & Pacala, S. W. (2016). Optimal stomatal behavior with competition for water and risk of hydraulic impairment. *Proceedings of the National Academy of Sciences*, 113(46), E7222–E7230. https://doi.org/10.1073/pnas.1615144113

Wong, S.-C., Cowan, I. R., & Farquhar, G. D. (1985). Leaf conductance in relation to rate of CO₂ assimilation: III. Influences of water stress and photoinhibition. *Plant Physiology*, 78(4), 830–834. https://doi.org/10.1104/pp.78.4.830

Woodward, F. I. (1987). Stomatal numbers are sensitive to increases in CO₂ from pre-industrial levels. *Nature*, 327(6123), 617–618. https://doi.org/10.1038/327617a0

Wu, J., Rogers, A., Albert, L. P., Ely, K., Prohaska, N., Wolfe, B. T., et al. (2019). Leaf reflectance spectroscopy captures variation in carboxylation capacity across species, canopy environment and leaf age in lowland moist tropical forests. *New Phytologist*, 224(2), 663–674. https://doi.org/10.1111/nph.16029

Wu, J., Serbin, S. P., Ely, K. S., Wolfe, B. T., Dickman, L. T., Grossiord, C., et al. (2020). The response of stomatal conductance to seasonal drought in tropical forests. *Global Change Biology*, 26(2), 823–839. https://doi.org/10.1111/gcb.14820

Xu, X., & Trugman, A. T. (2021). Trait-based modeling of terrestrial ecosystems: Advances and challenges under global change. *Current Climate Change Reports*, 7(1), 1–13. https://doi.org/10.1007/s40641-020-00168-6

Yang, J., Duursma, R. A., De Kauwe, M. G., Kumarathunge, D., Jiang, M., Mahmud, K., et al. (2019). Incorporating non-stomatal limitation improves the performance of leaf and canopy models at high vapour pressure deficit. *Tree Physiology*. https://doi.org/10.1093/treephys/tpz103

Yunusa, I. A. M., Eamus, D., Taylor, D., Whitley, R., Gwenzi, W., Palmer, A. R., & Li, Z. (2015). Partitioning of turbulent flux reveals contrasting cooling potential for woody vegetation and grassland during heat waves. *Quarterly Journal of the Royal Meteorological Society*, 141(692), 2528–2537. https://doi.org/10.1002/qj.2539

Zenes, N., Kerr, K. L., Trugman, A. T., & Anderegg, W. R. L. (2020). Competition and drought alter optimal stomatal strategy in tree seedlings. *Frontiers in Plant Science*, 11, 478. https://doi.org/10.3389/fpls.2020.00478

Zhou, S., Duursma, R. A., Medlyn, B. E., Kelly, J. W. G., & Prentice, I. C. (2013). How should we model plant responses to drought? An analysis of stomatal and non-stomatal responses to water stress. *Agricultural and Forest Meteorology*, 182–183, 204–214. https://doi.org/10.1016/j.agrformet.2013.05.009

## References From the Supporting Information

Anderegg, W. R. L., Wolf, A., Arango-Velez, A., Choat, B., Chmura, D. J., Jansen, S., et al. (2018). Woody plants optimise stomatal behaviour relative to hydraulic risk. *Ecology Letters*, 21(7), 968–977. https://doi.org/10.1111/ele.12962

Choat, B., Ball, M. C., Luly, J. G., Donnelly, C. F., & Holtum, J. A. M. (2006). Seasonal patterns of leaf gas exchange and water relations in dry rain forest trees of contrasting leaf phenology. *Tree Physiology*, 26(5), 657–664. https://doi.org/10.1093/treephys/26.5.657

Collatz, G. J., Ball, J. T., Grivet, C., & Berry, J. A. (1991). Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: A model that includes a laminar boundary layer. *Agricultural and Forest Meteorology*, 54(2–4), 107–136. https://doi.org/10.1016/0168-1923(91)90002-8

De Kauwe, M. G., Lin, Y.-S., & Medlyn, B. E. (2015). FitFarquharmodel: V<sub>cmax</sub> one-point method. *Zenodo*. https://doi.org/10.5281/ZENODO.30954

De Pury, D. G. G., & Farquhar, G. D. (1997). Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models. *Plant, Cell and Environment*, 20(5), 537–557. https://doi.org/10.1111/j.1365-3040.1997.00094.x

Farquhar, G. D., von Caemmerer, S., & Berry, J. A. (1980). A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species. *Planta*, 149(1), 78–90. https://doi.org/10.1007/BF00386231
Héroult, A., Lin, Y.-S., Bourne, A., Medlyn, B. E., & Ellsworth, D. S. (2013). Optimal stomatal conductance in relation to photosynthesis in climatically contrasting *Eucalyptus* species under drought: Stomatal responses of eucalyptus under drought. *Plant, Cell & Environment, 36*(2), 262–274. https://doi.org/10.1111/j.1365-3040.2012.02570.x

Hyndman, R. J., & Koehler, A. B. (2006). Another look at measures of forecast accuracy. *International Journal of Forecasting*, 22(4), 679–688.
https://doi.org/10.1016/j.ijforecast.2006.03.001

Kirschbaum, M., & Farquhar, G. (1984). Temperature dependence of whole-leaf photosynthesis in *Eucalyptus pauciflora* Sieb. Ex Spreng. *Functional Plant Biology*, 11(6), 519. https://doi.org/10.1071/PP9840519

Lasdon, L., Duarte, A., Glover, F., Laguna, M., & Martí, R. (2010). Adaptive memory programming for constrained global optimization. *Computers & Operations Research*, 37(8), 1500–1509. https://doi.org/10.1016/j.cor.2009.11.006

Limousin, J.-M., Bickford, C. P., Dickman, L. T., Pangle, R. E., Hudson, P. J., Boutz, A. L., et al. (2013). Regulation and acclimation of leaf gas exchange in a piñon-juniper woodland exposed to three different precipitation regimes: Rainfall manipulation in piñon-juniper woodland. *Plant, Cell & Environment*, 36(10), 1812–1825. https://doi.org/10.1111/pce.12089

Lu, Y., Duursma, R. A., Farrior, C. E., Medlyn, B. E., & Feng, X. (2020). Optimal stomatal drought response shaped by competition for water and hydraulic risk can explain plant trait covariation. *New Phytologist*, 225(3), 1206–1217. https://doi.org/10.1111/nph.16207

Martin-StPaul, N. K., Limousin, J.-M., Rodríguez-Calcerrada, J., Ruffault, J., Rambal, S., Letts, M. G., & Misson, L. (2012). Photosynthetic sensitivity to drought varies among populations of Quercus ilex along a rainfall gradient. *Functional Plant Biology*, 39(1), 25. https://doi.org/10.1071/FP11090

Medlyn, B. E., Dreyer, E., Ellsworth, D., Forstreuter, M., Harley, P. C., Kirschbaum, M. U. F., et al. (2002). Temperature response of parameters of a biochemically based model of photosynthesis. II. A review of experimental data. *Plant, Cell and Environment*, 25(9), 1167–1179. https://doi.org/10.1046/j.1365-3040.2002.00891.x

Mitchell, P. J., Veneklaas, E., Lambers, H., & Burgess, S. S. O. (2009). Partitioning of evapotranspiration in a semi-arid eucalypt woodland in south-western Australia. *Agricultural and Forest Meteorology*, 149(1), 25–37. https://doi.org/10.1016/j.agrformet.2008.07.008

Monteith, J. L., & Unsworth, M. H. (1990). *Principles of environmental physics* (2nd ed.). E. Arnold; Distributed in the USA by Routledge, Chapman and Hall.

Nash, J. E., & Sutcliffe, J. V. (1970). River flow forecasting through conceptual models part I—A discussion of principles. *Journal of Hydrology*, 10(3), 282–290. https://doi.org/10.1016/0022-1694(70)90255-6

Nelder, J. A., & Mead, R. (1965). A simplex method for function minimization. *The Computer Journal*, 7(4), 308–313. https://doi.org/10.1093/comjnl/7.4.308

Neufeld, H. S., Grantz, D. A., Meinzer, F. C., Goldstein, G., Crisosto, G. M., & Crisosto, C. (1992). Genotypic variability in vulnerability of leaf xylem to cavitation in water-stressed and well-irrigated sugarcane. *Plant Physiology*, 100(2), 1020–1028. https://doi.org/10.1104/pp.100.2.1020

Powell, M. J. D. (1964). An efficient method for finding the minimum of a function of several variables without calculating derivatives. *The Computer Journal*, 7(2), 155–162. https://doi.org/10.1093/comjnl/7.2.155

Powell, M. J. D. (1994). A direct search optimization method that models the objective and constraint functions by linear interpolation. In S. Gomez & J.-P. Hennart (Eds.), *Advances in optimization and numerical analysis* (pp. 51–67). Springer Netherlands. https://doi.org/10.1007/978-94-015-8330-5_4

Prentice, I. C., Dong, N., Gleason, S. M., Maire, V., & Wright, I. J. (2014). Balancing the costs of carbon gain and water transport: Testing a new theoretical framework for plant functional ecology. *Ecology Letters*, 17(1), 82–91. https://doi.org/10.1111/ele.12211

Storn, R., & Price, K. (1997). Differential Evolution – A simple and efficient heuristic for global optimization over continuous spaces. *Journal of Global Optimization*, 11(4), 341–359. https://doi.org/10.1023/A:1008202821328

Venables, W. N., & Ripley, B. D. (2010). *Modern applied statistics with S* (4th ed., [Nachdr.]). Springer.

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. *Nature Methods*, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2

Wales, D. J., & Doye, J. P. K. (1997). Global optimization by basin-hopping and the lowest energy structures of Lennard-Jones clusters containing up to 110 atoms. *The Journal of Physical Chemistry A*, 101(28), 5111–5116. https://doi.org/10.1021/jp970984n

Wang, Y., Sperry, J. S., Anderegg, W. R. L., Venturas, M. D., & Trugman, A. T. (2020). A theoretical and empirical assessment of stomatal optimization modeling. *New Phytologist*, 227(2), 311–325. https://doi.org/10.1111/nph.16572

Wu, J., Serbin, S. P., Ely, K. S., Wolfe, B. T., Dickman, L. T., Grossiord, C., et al. (2020). The response of stomatal conductance to seasonal drought in tropical forests. *Global Change Biology*, 26(2), 823–839. https://doi.org/10.1111/gcb.14820

Xiang, Y., & Gong, X. G. (2000). Efficiency of generalized simulated annealing. *Physical Review E*, 62(3), 4473–4476. https://doi.org/10.1103/PhysRevE.62.4473