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

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.

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 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.

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⁻¹. 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.

Table 1 Summary of the Stomatal Conductance (g_s) Schemes

Approach	gs form or optimization objectivea	Model name	Conductance parameter(s)	Reduction/cost parameter(s)a	Specific equation(s)	Main reference(s)
Empirical	gs $$\tilde{\propto }$$ RF	Medlyn	g1,Med	sMed	Equations 4 and 5	De Kauwe et al. (2015) and Medlyn et al. (2011)
		Tuzet	kmax, g1,Tuz	Ψref, sTuz b	Equations 6 and 7	Duursma and Medlyn (2012) and Tuzet et al. (2003)
		Eller	kmax	–	Equations 8 and 9	Eller et al. (2020)
Optimization	Avoiding primary impairment to water flow: max(An − CT)	WUE-$${f}_{{{\Psi}}_{l}}$$	kmax	λ	Equation 11	Wolf et al. (2016)
		CMax	kmax	a, b	Equation 13	Anderegg et al. (2018) and Wolf et al. (2016)
		ProfitMax	kmax	–	Equation 14	Sperry et al. (2017)
		CGain	kmax	ϖ	Equation 15	Lu et al. (2020)
	Avoiding primary impairment to carbon assimilation: max(An × RF)	SOXopt	kmax	–	Equation 16	Eller et al. (2018)
		ProfitMax2	kmax	–	Equation 17	Wang et al. (2020)
		LeastCost	kmax	η	Equation 18	Prentice et al. (2014)
	Non-stomatal optimization: max(An), with internal RF on Vcmax or Ci	CAP	kmax	Ψφ,lim	Equations 19–21	Dewar et al. (2017)
		MES	kmax	Ψφ,lim	Equations 19, 20, and 22

^aRF stands for reduction factor and CT for cost term.

^bParameter commonly set to an arbitrary value, that is, s_Tuz is often set to 2 for the Tuzet model.

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).

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: [Image Omitted. See PDF]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_l (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 SOX_opt; 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).

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): [Image Omitted. See PDF]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₉₅).

k_Ψ is downregulated by a cumulative Weibull distribution representing plant vulnerability to cavitation (Neufeld et al., 1992; Sperry et al., 2017): [Image Omitted. See PDF]where k_max (mmol m⁻² s⁻¹ MPa⁻¹) is the maximum root-zone-to-leaf hydraulic conductance, and s₁ (MPa) and s₂ (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₁ and s₂) 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.

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): [Image Omitted. See PDF]where C_s (μmol mol⁻¹) is [CO₂] at the leaf surface, 1.57 converts from conductance to CO₂ to conductance to water vapor, g₀ (mol m⁻² s⁻¹) is the residual conductance to water vapor (assumed negligible), g_1,Med (kPa^0.5) is the slope of the sensitivity of g_s to A_n, and β (unitless) is an empirical moisture stress factor represented by an exponential dependency on Ψ_l at predawn (i.e., Ψ_l,pd; Cowan, 1982; Yang et al., 2019; Zhou et al., 2013): [Image Omitted. See PDF]where s_Med (MPa⁻¹) sets the sensitivity of the stomates to Ψ_l,pd, Ψ_fc is the root-zone water potential at field capacity, and with Ψ_l,pd ∼ Ψ_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 Ψ_l, which Duursma and Medlyn (2012) later simplified to: [Image Omitted. See PDF]where g_1,Tuz (unitless) is an empirical slope parameter and ζ (unitless) characterizes stomatal responses to Ψ_l: [Image Omitted. See PDF]where Ψ_ref (MPa) is a reference water potential and s_Tuz (MPa⁻¹) is the sensitivity of the stomates to variations in Ψ_l.

The Eller Model

Eller et al. (2020) proposed an analytical approximation to the stomatal optimization based on xylem hydraulics (cf., the SOX_opt model below), establishing an empirical link between g_s and plant hydraulic vulnerability: [Image Omitted. See PDF]where C_i,col (μmol mol⁻¹) is the co-limitation C_i—for which biochemical rates of assimilation are equal—and ξ (mol m⁻² s⁻¹) is a reduction factor arising from xylem conductance loss: [Image Omitted. See PDF]where k_Ψ is given by Equation 3, P₅₀ (MPa) is the water potential drop at 50% loss of xylem hydraulic conductivity, 10³ converts to mmol, and Ψ_x,pd is leaf xylem water pressure at predawn, with Ψ_x,pd ∼ Ψ_s−10⁻⁶ρgh following Archimede's principle, where ρ (kg m⁻²) is the density of water, g (m s⁻²) is gravitational acceleration, h (m) is canopy height, and 10⁻⁶ converts to MPa.

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., Ψ_l vs. k_Ψ) 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_i 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₂: [Image Omitted. See PDF]where C_a (μmol mol⁻¹) is ambient [CO₂], and 1.35 converts from boundary layer conductance to CO₂ 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).

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” (λ; μmol CO₂ mmol⁻¹ H₂O), by regulation of their stomates to maximize: [Image Omitted. See PDF]

However, the WUEH assumes λ 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 Ψ_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 Ψ_s and Ψ_l (Equation 2). This proposition keeps λ 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 λE with a concave-up parabola that increases as Ψ_l drops (Anderegg et al., 2018): [Image Omitted. See PDF]where a (μmol m⁻² s⁻¹ MPa⁻²) and b (μmol m⁻² s⁻¹ MPa⁻¹) and c (μmol m⁻² s⁻¹) 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 Ψ_l: [Image Omitted. See PDF]a can be interpreted as the carbon cost of a falling Ψ_l, whilst |b| 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: [Image Omitted. See PDF]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: [Image Omitted. See PDF]

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

Stomatal Optimizations Avoiding Primary Impairment to Carbon Assimilation (Three Models) The SOX_opt Model

Inspired by the ProfitMax model, the stomatal optimization based on xylem hydraulics (SOX_opt; 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: [Image Omitted. See PDF]

Eller et al. (2018) made SOX_opt'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_opt 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): [Image Omitted. See PDF]where E_crit (mmol m⁻² s⁻¹) is calculated from Equation 2 and Ψ_crit triggers a 95% loss of hydraulic conductivity, with Ψ_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: [Image Omitted. See PDF]where η (μmol CO₂ mmol⁻¹ H₂O) is the carbon cost of photosynthetic proteins maintenance through transpiration, and V_cmax (μmol m⁻² s⁻¹) is the maximum rate of carboxylation.

As in the WUE-$${f}_{{{\Psi}}_{l}}$$ 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.

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 Ψ_l. [Image Omitted. See PDF]

Hydraulic stress is represented through a unitless linear factor: [Image Omitted. See PDF]where Ψ_φ,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): [Image Omitted. See PDF]

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: [Image Omitted. See PDF]where Γ^∗ (μmol mol⁻¹) is the photorespiration compensation point.

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₅₀ = −3.1 MPa, P₈₈ = −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^0.5; 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 2 Environmental Variable Ranges Used for the Idealized Model Calibrations and Sensitivity Analysis

[Table omitted. See PDF]

^aThe ranges selected are chosen to yield realistic model estimates.

Idealized Calibrations 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}}_{l}}$$, 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}}_{l}}$$, 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., λ in WUE-$${f}_{{{\Psi}}_{l}}$$) 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 − Ψ_ref is equifinal.

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 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”).

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 Ψ_s) affects modeled g_s, Ψ_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 Ψ_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, Ψ_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-Ψ_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-Ψ_s synergies) to the total variability in g_s.

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, Ψ_l for 15 sites × species, and C_i for 12 sites × species. Observed PPFD, T_a, D_a, C_a, g_b, and Ψ_l,pd (used as a proxy for Ψ_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.

Table 3 Site and Observation Data Set Characteristics, Sample Size (N), and Hydraulic and Photosynthetic Traits Used to Parameterize the Models

Site	Lat.	Lon.	MAP (mm)	Species	N	PPFD (μmol m−2 s−1)	Da (kPa)	Ψs (-MPa)	Data set ref.	P50; P88 a (-MPa)	Hydraulic traits ref.	Vcmax,25; Jmax,25 a (μmol m−2 s−1)	Photosynthetic traits ref.
San Lorenzo Protected Area, Panamab	9°16′N	79°58′W	3,286	Carapa guianensis	36	159–2,001	0.5–1.5	0.9–0.48	Wu et al. (2020)	0.8 ± 0.1; 2.3 ± 0.2	Lopez et al. (2005)	41.6 ± 1.6; 90.7 ± 2.3	Rogers et al. (2019)
				Tachigali versicolor	30	121–2,001	0.4–1.2	1.8–0.55		1.6; 4	Meinzer et al. (2008)	53.7 ± 1.8; 99.6 ± 1.6
				Tocoyena pittieri	32	100–2,000	0.5–1.5	1–0.55		3.3 ± 0.1; 4.3 ± 0.3	Wu et al. (2020)	50.8 ± 1.6; 83.9 ± 1.2
Parque Natural Metropolitano, Panamab	8°59′N	79°32′W	1,826	Calycophyllum candidissimum	32	51–2,001	0.5–4.3	2.2–0.33		2.9 ± 0.1; 4.3 ± 0.1	Brodribb et al. (2003)	77.4 ± 5.2; 112.1 ± 2.5
Many Peaks Range, QLD, Australiac	19°11′S	145°45′E	1,100	Alphitonia excelsa	160	93–1,854	1–6.2	2.2–0.2	Choat et al. (2006)	5.6; 9.1	Choat et al. (2003)	34.2 ± 3.7; 68.9	Choat et al. (2006)
				Austromyrtus bidwillii	55	52–1,823	1.1–4.9	0.9–0.2		5.1; 12		23.5 ± 3.8; 47.4
				Brachychiton australis	103	117–2,146	0.9–7	0.8–0.18		3.2; 5.2		48.4 ± 6.8; 97.6
				Cochlospermum gillivraei	84	226–1,967	1.8–6.8	1.1–0.28		1.4; 2.1		42.3 ± 2.5; 85.3
Richmond, NSW, Australiad	33°19′S	138°32′E	462	Eucalyptus cladocalyx	90	398–2,001	0.6–5.1	1.4–0.5	Héroult et al. (2013)	5.3 ± 0.4; 8.1 ± 0.9	Bourne et al. (2017)	98 ± 11; 160.7	Héroult et al. (2013)
	28°47′S	153°09′E	1,622	Eucalyptus dunnii	92	399–2,002	0.5–5.5	1.2–0.05		4.9 ± 0.2; 8.4 ± 1		115 ± 11; 188.6
	30°58′S	152°11′E	945	Eucalyptus saligna	78	193–2,001	0.5–5.9	1.3–0.2		3.4 ± 0.2; 5.2 ± 0.6		96 ± 9; 157.4
Corrigin, WA, Australiae	33°03′S	138°75′E	361	Eucalyptus capillosa	132	75–2,202	1.2–8.9	4.4–3.2	Mitchell et al. (2009)	3.1 ± 0.1; 5.7 ± 0.2	Mitchell (2009)	71.7; 117.5	de Boer et al. (2016)
Sevilleta, NM, USA	34°23′N	106°31′W	363	Juniperus monosperma	187	2,000	0.5–6.7	5.6–0.85	Limousin et al. (2013)	9.4 ± 0.4; 11.8 ± 0.2	Hudson et al. (2018)	35.6 ± 1.5; 68.1 ± 1.7	Limousin et al. (2013)
				Pinus edulis	139	2,000	0.7–6.2	3–0.8		4.4 ± 0.1; 5.4 ± 0.1		17.4 ± 1.6; 36.1 ± 2.2
Vic la Gardiole, Franceb	43°29′N	3°45′E	650	Quercus ilex	45	1,227–1,668	1–4.2	3.4–0.35	Martin-StPaul et al. (2012)	6.9; 9.4	Martin-StPaul et al. (2017b)	56.9 ± 6.7; 98.6 ± 4.2	Martin-StPaul et al. (2012)
Puéchabon, Franceb	43°44′N	3°35′E	907		32	1,494–1,515	1–5.2	3.2–0.4				60.9 ± 3.4; 113.8 ± 3.3

^aWhenever 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.

^bV_cmax,25 and J_max,25 were fitted from A to C_i curves (see Text S5 in Supporting Information S1).

^cV_cmax,25 was estimated using the “one-point method” for light-saturated rates of photosynthesis (PPFD >500 μmol m⁻² s⁻¹) as in De Kauwe et al. (2016) and J_max,25 is taken as 2.02 × V_cmax,25 as per Kumarathunge et al. (2019) on average for rainforest species.

^dFor 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_max,25 was taken as 1.64 × V_cmax,25 as per Kumarathunge et al. (2019) on average for Eucalyptus species.

^eThe V_cmax,25 is that of a relative, Eucalyptus wandoo, and J_max,25 was taken as 1.64 × V_cmax,25 as per Kumarathunge et al. (2019) on average for Eucalyptus species.

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 Ψ_l over all environmental conditions, but the observations were also inconsistently distributed in time, with variable tree age and size classes sampled. The Ψ_ref and s_TUZ parameters (Equation 7) were directly fitted from the observed g_s and Ψ_l (Text S6, Figures S6a–S6o in Supporting Information S1), and Ψ_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 Ψ_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.

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.

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).

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.

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 (ϖ in CGain, η in LeastCost, λ 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 Ψ_l in optimization schemes—Ψ_φ,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.

Do Model Behaviors Differ in Response to Abiotic Variations?

Figure 2 shows the (mostly out-of-sample) influence of atmospheric drivers and soil moisture (Ψ_s) on g_s, Ψ_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 2. Total-order Sobol’ sensitivity indices of the stomatal conductance (gs), leaf water potential (Ψl), and CO2 concentration in the leaf intercellular air spaces (Ci) to variability in environmental drivers for the 12 models parameterized under well-watered conditions. The environmental drivers are: atmospheric vapor pressure deficit (Da); soil water potential (Ψs); ambient CO2 concentration (Ca); photosynthetic photon flux density; and air temperature (Ta). 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.

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 Ψ_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₂ 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, Ψ_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 β downregulation factor in Medlyn leads to Ψ_s having a major impact on both g_s and C_i.

The picture is more contrasted when looking at the Ψ_l responses (cf., Appendix A). Eller, CMax, and LeastCost are mostly unaffected by variability in Ψ_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 Ψ_s in three ways: (a) responses to variability in Ψ_s are smaller than to variability in D_a and/or T_a (Tuzet, WUE-$${f}_{{{\Psi}}_{l}}$$); (b) sensitivity to Ψ_s either slightly dominates or is on par with sensitivity to other drivers (ProfitMax, CGain, SOX_opt); and (c) variability is primarily driven by Ψ_s (ProfitMax2, CAP, MES). We expect that the latter scenario would become prevalent if Ψ_s were dropping further than the range tested here, and it should happen sooner for scenario (b) than scenario (a).

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

How Well Do the g_s Schemes Capture Observed Leaf-Level Fluxes and Hydraulic Status? Predictive Performance

We used four metrics to assess the modeled variables against observations. Pooling all sites × species, the ProfitMax 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_opt, ProfitMax) for g_s across sites × species. Additionally, half the new models (Eller, ProfitMax, SOX_opt, 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_s, respectively), whilst an additional model displays equivalent performance (CGain).

Figure 3. Average model performance across sites × species for: gs, Ci, E, An, and Ψl (in red); all variables excluding Ψl (in dark blue), and gs 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 [less than]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).

Given the models were calibrated to observations of g_s 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_s behavior than in the extra-tropics (Figure 4). Finally, among variables, Ψ_l is responsible for the largest performance degradation, particularly for the Eller and LeastCost models, which we examine below.

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Ψl ${f}_{{{\Psi}}_{l}}$) to most (LeastCost) skilled.

Functional Relationships

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

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 (Vic la 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.

Predictions of E and A_n 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_n. Indeed, each model occupies its own intrinsic WUE space, depending on its ability to regulate gas exchange with D_a for example (Medlyn, Tuzet, WUE-$${f}_{{{\Psi}}_{l}}$$, CMax and CGain do not capture the behavior at high D_a in Figure 6d). Therefore, individual models can be biased with regards to E, A_n, 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_cmax,25 may cause A_n 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_s, C_i and Ψ_l.

Figure 6. Leaf-level estimates of carbon assimilation (An) 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 An 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 (Da) is low; and (iii) the average behavior when Da is high. Both high and low Da 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.

The models generally agree on the magnitude of C_i per site × species (Figure 7 and Figure S11 in Supporting Information S1), noting that the models that simulate relatively lower C_i simulate marginally lower A_n (Eller compared to ProfitMax2 in Figures 6a and 7a), and vice versa. Similarly, relatively lower E is associated with relatively lower g_s (LeastCost compared to Tuzet in Figures 6a and 7a), or vice versa. Whereas other models tend toward similar magnitude overestimation of g_s and C_i, or higher g_s than C_i, CAP and MES overestimate C_i more than they do g_s, in fact largely so given their estimates of A_n (Figures 6 and 7; Figures S10 and S11 in Supporting Information S1; see also Appendix A).

Figure 7. Functional relationships between stomatal conductance (gs), CO2 concentration in the leaf intercellular air spaces (Ci), and leaf water potential (Ψ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 gs ratios against simulated to observed Ci 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 gs with decreasing Ψl, fitted via a generalized additive model, compared to the observations (light gray crosses). The functional forms were made comparable by normalizing gs by its model-specific maximum per site × species, and Ψl by the critical leaf water potential indicative of total xylem failure (P95 in this study). Curves were not fitted if they did not monotonically decrease (e.g., WUE-fΨl ${f}_{{{\Psi}}_{l}}$ and CMax in panel [d]), or where the models operate at, or beyond, the P95 (e.g., the Eller model in panel [d]).

There are four sets of characteristics that describe g_s to Ψ_l 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 Ψ_l 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}_{{{\Psi}}_{l}}$$ and particularly CMax, g_s(Ψ_l) 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(Ψ_l) compared to the observations (i.e., accurate functional shape and variability, regardless of bias) appears key to representing stomatal closure as atmospheric aridity increases. For example, CGain, SOX_opt, 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 MES's g_s(Ψ_l) 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.

Figure 8. A comparison of the sensitivity of 12 stomatal conductance (gs) models to vapor deficit (Da) across site × species. The driest two thirds of soil moisture were excluded within each site × species to isolate Da effects from soil moisture effects. Panel (a) shows the binned average effect of Da on gs (normalized by the observed maximum gs per site × species) over the range of observed Da (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 gs data shown in (b) and includes a minimum Da of 2 kPa to account for models that predict relatively early gs closure. Panel (b) shows ridge plots of the differences between rates of simulated gs decline (i.e., simulated gs normalized by its model-specific maximum per site × species and scaled to the maximum observation) and rates of observed gs decline, relative to observed gs decline from within the boxed area from (a). Negative spread and peaks signify early stomatal closure compared to the observations, and vice versa.

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 Ψ_l does not always degrade the simulation of gas exchange (LeastCost, Eller); however, it may lead to discrepant responses at high D_a.

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 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).

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 Ψ_l,pd (or Ψ_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 Ψ_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 Ψ_l is obtained by accounting for plant hydraulics (Tuzet et al., 2003). The problem with this approach is 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₅₀ of −1.8 MPa (“BET-Te” in Table 2 of Eller et al. [2020]), whereas measured P₅₀ 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.

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}_{{{\Psi}}_{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_opt, 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}_{{{\Psi}}_{l}}$$ 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 Ψ_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 Ψ_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}_{{{\Psi}}_{l}}$$ and (mostly) CMax may behave unrealistically at high D_a. We also concur that ProfitMax and SOX_opt 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 Ψ_l in Figures 5d and 5e) give ProfitMax the edge over SOX_opt.

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 Ψ_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.

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 Ψ_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 Ψ_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 Ψ_l and k_max.

The inclusion of reduction/cost parameters in optimality criteria (e.g., λ in WUE-$${f}_{{{\Psi}}_{l}}$$, η in LeastCost, Ψ_φ,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 Ψ_l compared to one-parameter models which required calibrating k_max only.

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 Ψ_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., Ψ_crit, k_crit) that could reduce model bias in Ψ_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 Ψ_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 Ψ_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 Ψ_l and plant water status.

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_opt), 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.

A AppendixHarmonized 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_s (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 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 (gs, panel a) was used as a reference to calibrate the other models. Other simulated fluxes shown are: (b) the CO2 concentration in the leaf intercellular air spaces (Ci) as a proportion of the ambient CO2 concentration (Ca); (c) the leaf water potential (Ψl); (d) the net rate of carbon assimilation (An); and (e) the transpiration (E).

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 (gs, panel a) was used as a reference to calibrate the other models. Other simulated fluxes shown are: (b) the CO2 concentration in the leaf intercellular air spaces (Ci) as a proportion of the ambient CO2 concentration (Ca); (c) the leaf water potential (Ψl); (d) the net rate of carbon assimilation (An); and (e) the transpiration (E).

Compared to Medlyn in Figure A1a, five optimal models (CMax, ProfitMax, CGain, SOX_opt, 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 SOX_opt 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 steadiest, 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₅₀ 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 profligate 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 (Equation 22). Consequently, both models must maintain high C_i (Figure A1b) to prevent stark reductions in A_n.

Simulation of Ψ_l (Figure A1c) is where schemes most diverge. LeastCost and Eller simulate lower Ψ_l than all other schemes, even beyond the P₈₈. These models' calibrated k_max is smaller than the ProfitMax's (Table S6 in Supporting Information S1), so LeastCost and Eller have lower Ψ_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 Ψ_l. Both CGain and SOX_opt behave similarly to ProfitMax, with ever so marginally lower Ψ_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 Ψ_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.

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 Scientia PhD Scholarship Scheme. M. G. De Kauwe and A. J. Pitman acknowledge 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 Biological 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 anonymous reviewers for their constructive comments.

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).