Research Stomatal optimization based on xylem hydraulics (SOX) improves land surface model simulation of vegetation responses to climate Cleiton B. Eller1,2 , Lucy Rowland1 , Maurizio Mencuccini3,4 , Teresa Rosas3,4 , Karina Williams5 , Anna Harper6 , Belinda E. Medlyn7 , Yael Wagner8 , Tamir Klein8 , Grazielle S. Teodoro9 , Rafael S. Oliveira2 , Ilaine S. Matos10 , Bruno H. P. Rosado10 , Kathrin Fuchs11 , Georg Wohlfahrt12 , Leonardo Montagnani13 , Patrick Meir14,15 , Stephen Sitch1 and Peter M. Cox6 1 College of Life and Environmental Sciences, University of Exeter, Exeter, EX4 4QF, UK; 2Department of Plant Biology, University of Campinas, Campinas 13083-862, Brazil; 3CREAF, Bellaterra 08193, Barcelona, Spain; 4ICREA, Pg. Lluıs Companys 23, 08010 Barcelona, Spain; 5Met Office Hadley Centre, Fitzroy Road, Exeter, EX1 3PB, UK; 6College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK; 7Hawkesbury Institute for the Environment, Western Sydney University, Locked Bag 1797, Penrith, NSW 2751,Australia; 8Department of Plant & Environmental Sciences, Weizmann Institute of Science, 76100, Rehovot, Israel; 9Institute of Biological Sciences, Federal University of Para, Belem 66075-110, Brazil; 10Department of Ecology – IBRAG, Rio de Janeiro State University (UERJ), Rio de Janeiro 20550-013, Brazil; 11Department of Environmental Systems Science, ETH Zurich, Universit€atstrasse 2, 8092 Zurich, Switzerland; 12Department of Ecology, University of Innsbruck, Innsbruck 6020, Austria; 13Forest Services, Autonomous Province of Bolzano, Via Brennero 6, 39100 Bolzano, Italy; 14Research School of Biology, The Australian National University, Acton, ACT 2601, Australia; 15School of Geosciences, University of Edinburgh, Edinburgh, EH9 3FF, UK Summary Author for correspondence:  Land surface models (LSMs) typically use empirical functions to represent vegetation Cleiton B. Eller responses to soil drought. These functions largely neglect recent advances in plant ecophysiol- Tel: +44 07864 755668 ogy that link xylem hydraulic functioning with stomatal responses to climate. Email: c.breder-eller@exeter.ac.uk  We developed an analytical stomatal optimization model based on xylem hydraulics (SOX) Received: 5 November 2019 to predict plant responses to drought. Coupling SOX to the Joint UK Land Environment Simu- Accepted: 3 January 2020 lator (JULES) LSM, we conducted a global evaluation of SOX against leaf- and ecosystem- level observations. New Phytologist (2020) 226: 1622–1637  SOX simulates leaf stomatal conductance responses to climate for woody plants more accu- doi: 10.1111/nph.16419 rately and parsimoniously than the existing JULES stomatal conductance model. An ecosys- tem-level evaluation at 70 eddy flux sites shows that SOX decreases the sensitivity of gross Key words: climate change, drought, eddy primary productivity (GPP) to soil moisture, which improves the model agreement with obser- covariance, land-surface models, stomatal vations and increases the predicted annual GPP by 30% in relation to JULES. SOX decreases optimization, xylem hydraulics. JULES root-mean-square error in GPP by up to 45% in evergreen tropical forests, and can sim- ulate realistic patterns of canopy water potential and soil water dynamics at the studied sites.  SOX provides a parsimonious way to incorporate recent advances in plant hydraulics and optimality theory into LSMs, and an alternative to empirical stress factors. 2018b), thereby restricting our capability to predict the effect of Introduction increased aridity on vegetation distribution and its feedbacks on Large areas of the globe will be exposed to increased aridity in the the global carbon cycle and climate. Many LSMs represent the near future (Sheffield & Wood, 2008; Duffy et al., 2015; effects of reduced soil moisture on canopy carbon assimilation Marengo et al., 2018). As drought events become more intense (A) using an empirical drought factor commonly referred as b- and frequent, accurately representing vegetation–climate feed- factor (Cox et al., 1998). The b-factor approach has been shown backs in Earth system models (ESMs) is increasingly important, to overestimate plant responses to seasonal and experimentally as these interactions can drastically influence model projections induced drought (Ukkola et al., 2016; Restrepo-Coupe et al., of global climate change (Cox et al., 2000). The current genera- 2017; Eller et al., 2018b). The b-factor has a large impact on the tion of land surface models (LSMs) does not accurately simulate modelled global carbon budget, supressing 30–40% of the vegetation carbon dynamics during drought (Sitch et al., 2008; annual gross primary productivity (GPP) in large areas of arid Powell et al., 2013; Medlyn et al., 2016; Ukkola et al., 2016; and semiarid ecosystems (Trugman et al., 2018). Despite its Restrepo-Coupe et al., 2017; Rogers et al., 2017; Eller et al., importance, there is scarce empirical support for the drought 1622 New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust 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. New Phytologist Research 1623 functions used in most LSMs (Medlyn et al., 2016). The lack of (Kattge et al., 2011; Choat et al., 2012) to simulate gs responses a theoretical or empirical basis for the b-factor implies an urgent to environmental conditions, producing a realistic gs decline in need for new modelling approaches to replace this important response to elevated atmospheric CO2 and soil drought (Sperry component of LSMs so as to improve our capacity to predict veg- et al., 2017; Venturas et al., 2018; Eller et al., 2018b; Wang etation–climate interactions. et al., 2019). This approach predicts a tight coordination between Stomatal responses of plants to soil drought involve complex stomatal and xylem functioning which is widely corroborated by chemical signalling and hydrodynamic processes in leaf cells, observations (Hubbard et al., 2001; Meinzer et al., 2009; Klein, some of which have not been entirely elucidated (Buckley, 2017, 2014). Another advantage of this approach is its capacity to simu- 2019; Qu et al., 2019). Stomatal optimization models are a use- late a diversity of contrasting stomatal behaviours, from iso- to ful approach to model stomatal behaviour that circumvents the anisohydric (Martinez-Vilalta et al., 2014; Klein, 2014). need to explicitly represent the physiological processes involved Sperry et al. (2017) proposes a model that assumes that, as in stomatal regulation. Optimization models employ a ‘goal-ori- xylem hydraulic conductance declines, the increased risk of ented’ approach, assuming that plant stomata behaviour has been hydraulic failure is the main fitness cost associated with low Ψ. selected through plant evolutionary history to maximize a given Eller et al. (2018b) adapted the Sperry et al. (2017) model into objective function (Cowan, 2002; Dewar et al., 2009; Prentice the stomatal optimization model based on xylem hydraulics et al., 2014; Buckley, 2017). The traditional approach to model (SOX), which differs from the Sperry et al. (2017) model princi- optimal stomatal behaviour is derived from the seminal work of pally by using a different optimization target. The SOX opti- Cowan & Farquhar (1977). This approach proposes that optimal mization target is based on the PGEN model (Friend, 1995), stomatal behaviour maximizes A minus the carbon cost of water which assumes that stomata optimize plant dry matter produc- lost (kE) over a given time interval, where E is transpiration and tion, represented by the product of photosynthesis and a linear k is the Lagrange multiplier that represents the carbon cost of a function of Ψ. The SOX model in Eller et al. (2018b) uses a unit of water lost. This model, hereafter labelled CF, after Cowan numerical routine to find the optimum gs. However, the PGEN and Farquhar, is capable of simulating many patterns of stomatal optimization target can also be found analytically (Friend & Cox, responses to climate over short timescales (Farquhar et al., 1980; 1995; Dewar et al., 2018). A parsimonious analytical formulation Berninger & Hari, 1993), and has provided the theoretical basis for SOX would facilitate its incorporation into existing LSMs for several widely used semi-empirical stomatal models (Jacobs, and provide a practical alternative to the b-function for mod- 1994; Leuning, 1995; Medlyn et al., 2011). However, CF pre- elling stomatal responses to drought at global scales. dicts that stomatal conductance (gs) increases in response to ele- In this study we develop an analytical approximation for the vated CO2 when A is Rubisco-limited, which contradicts most numerical SOX model presented in Eller et al. (2018b). We then observations (Mott, 1988; Medlyn et al., 2001). Other limita- create a new configuration for the Joint UK Land Environment tions are related to the k, as the CF hypothesis does not link k to Simulator (JULES; Best et al., 2011; Clark et al., 2011) that uses measurable plant traits or environmental quantities (Buckley, SOX to compute vegetation gs from environmental and plant 2017), and assumes k is constant over the period of reference hydraulic data. Using a global dataset of xylem hydraulic traits, (Cowan & Farquhar, 1977), which makes the original CF unable together with an extensive leaf gas-exchange and eddy covariance to predict long-term gs decline in response to soil moisture deple- dataset, we calibrate the SOX parameters and compare the tion. JULES-SOX performance to the default JULES using the b- Since the original CF work, many attempts have been made to function, across all major global biomes. Our goals in this paper incorporate the effects of declining soil moisture in the CF stom- are twofold: to test SOX agreement with global observations of gs atal optimization framework (Cowan, 1986; M€akel€a et al., 1996; to assess the generality of the underlying hypothesis in SOX, that Williams et al., 1996; Manzoni et al., 2013). Some of these is, that plant stomata evolved to balance carbon assimilation with attempts, such as the soil–plant–atmosphere (SPA) model of Wil- the loss of hydraulic conductance; and to evaluate the effect of liams et al. (1996), employ principles of plant hydraulics to con- SOX on JULES ecosystem-scale predictions of carbon and water strain stomatal optimization and have been successfully fluxes, and their agreement with observations. incorporated into LSMs (Bonan et al., 2014). The numerical approach used by SPA employs a hydraulic threshold to set a Materials and Methods lower water potential limit (Ψmin) for gs, which simulates a strict isohydric stomatal regulation (Fisher et al., 2006). Despite using Analytical SOX description plant hydraulics, SPA still relies on a water-use efficiency opti- mization similar to CF to model stomatal behaviour when The SOX central hypothesis can be summarized as ‘stomatal con- Ψ > Ψmin (Williams et al., 1996; Bonan et al., 2014). ductance (gs) is such as to maximize the product of leaf photosyn- Alternative routes to model plant optimal stomatal behaviour thesis and xylem hydraulic conductance’ and is given by: have been proposed recently (for a review, see Mencuccini et al., 2019a). These approaches circumvent the CF limitations by A½ci ðgs ÞK ½Wm ðgs Þ; Eqn 1 assuming that plant optimal stomatal behaviour minimizes the instantaneous fitness costs associated with low Ψ. These new where A is leaf net CO2 assimilation (mol CO2 m2 s1), which optimization models use widely available plant hydraulic traits is a function of leaf internal CO2 partial pressure (ci, in Pa), Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1624 Research Phytologist which is itself a function of stomatal conductance to CO2 (gs, We evaluated SOX as a stand-alone leaf-level model, and cou- mol m2 s1). The K is the normalized (0–1) xylem hydraulic pled to JULES (hereafter JULES-SOX). The leaf-level model was conductance computed as: evaluated against leaf gas exchange data as an ‘assumption cen- tred’ (sensu Medlyn et al., 2015) test of the hypothesis underlying 1 SOX. The JULES-SOX was then evaluated against ecosystem- K ðWÞ ¼ ; Eqn 2 ½1 þ ðW=W50 Þa  level eddy flux data, which constituted the first practical test of the utility of SOX for LSMs. where Ψ50 is Ψ when K = 0.5 and the parameter a gives the shape of the curve, with a higher a producing a steeper response to Ψ. JULES b-function description We use the mean (Ψm, MPa) of the canopy water potential at the predawn (Ψpd, MPa) and the canopy water potential (Ψc, MPa) The JULES model (Best et al., 2011; Clark et al., 2011) uses the to compute K with Eqn 2 to account for the gradual decline in Ψ Collatz et al. (1991, 1992) photosynthesis model for C3 and C4 along the soil to canopy hydraulic pathway (see details in Sup- plants (Notes S3) to produce unstressed rates of A based on the porting Information Notes S1). The gs value that maximizes Eqn colimitation of light, Rubisco carboxylation capacity, and the 1 is found at: transport of photoassimilates (for C3 plants) and PEPcarboxylase limitation (for C4 plants). The effect of soil moisture in A in the @AK default JULES is given by multiplying A by the b factor, com- ¼ 0: Eqn 3 puted using the b-function from Cox et al. (1998): @gs 8 < 1 for h [ hc The gs value that satisfies Eqn 3 was found numerically in Eller b ¼ hhh w for h w \h  hc ; Eqn 7 et al. (2018b), but a computationally efficient analytical solution : c hw 0 for h  hw is preferable for application in dynamic global vegetation models (DGVMs) and ESMs. We developed an analytical approximation for the optimal SOX gs using the partial derivatives of A with where h is the mean soil moisture in the root zone (m3 m3), and respect to ci and K with respect to Ψm. All steps of the model hc and hw, are the critical and wilting points, which are defined derivation are described in Notes S1. The resulting SOX equa- by Cox et al. (1998) as the h when soil Ψ is  0.033 tion for the optimal gs is: and  1.5 MPa, respectively. The default JULES formulation employs the Jacobs (1994) equation to predict ci from D, ca and sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! the CO2 compensation point, Γ (Pa): @A 4n gs ¼ 0:5 þ11 ; Eqn 4   @ci @A=@ci D c i ¼ f0 1  ðca  CÞ; Eqn 8 Dcrit The benefit of stomatal opening is represented here by the sen- sitivity of leaf photosynthesis to the internal CO2 concentration where f0 and Dcrit are empirical parameters (Jacobs, 1994; Cox (@A=@ci ). By contrast, the parameter ξ represents the cost of et al., 1998). stomatal opening in terms of loss of xylem conductivity under The JULES-SOX configuration replaces Eqns 7 and 8, com- low Ψpd and/or higher leaf-to-air vapour pressure (D, puting gs from environmental data and plant hydraulic inputs mol mol1): with Eqns 4 and 5. To compute A from the gs predicted by Eqn 4, we solved the limiting photosynthetic rates from the Collatz 2 et al. (1991, 1992) model as functions of ca and gs, as described n¼ : Eqn 5 1=K @K =@Wm rp 1:6D in Notes S3. Low ξ indicates high hydraulic costs occurring during drought Leaf-level SOX evaluation (i.e. lower Ψpd and higher D; Fig. S1). SOX simulates dynamic changes on the plant hydraulic resistance (rp), computing rp as a We used a global compilation of leaf gas exchange data to evalu- function of Ψpd and the plant minimum hydraulic resistance ate the SOX capacity to reproduce leaf stomatal responses of a (rpmin, m2 s MPa mol1 H2O): wide range of woody plants. This dataset contains observations compiled by Lin et al. (2015), complemented with other pub- rp;min lished and unpublished data (see Table S1 and Fig. S2 for addi- rp ¼  : Eqn 6 K Wpd tional information). In total, there are 3597 measurements of gs and Ψpd together with environmental variables used for driving Solving SOX main equations (Eqns 4, 5) requires computing the model, that is, incident photosynthetic active radiation (Ipar), the partial derivatives of A and K, @A/@ci and @K/@Ψm, respec- air temperature (Ta), ca and D. These data come from 30 woody tively. These derivatives were estimated numerically in this study plant species collected in 15 sites around the world (Fig. S2b). as described in Notes S2. The Ψpd was measured on the same day as gs, and the New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1625 environmental data was measured simultaneously with gs. The the photosynthesis model used in SOX (Notes S3) were set equal dataset included field and glasshouse observations, with environ- to those in Harper et al. (2016). mental conditions varying from well-watered to extreme drought The model predictive skill was evaluated using the model root- (Ψpd = 7 MPa). These observations were grouped into the mean-square errors (RMSE) and the Nash & Sutcliffe (1970) global plant functional type (PFT) categories from Harper et al. model efficiency index (NSE). The NSE varies from ∞ to 1, (2016) (Table 1). Harper et al. (2016) divides angiosperm tree with 1 indicating perfect agreement between model and observa- species into broadleaf evergreen tropical trees (BET-Tr), tions, while NSE < 0 indicates that the mean value of the obser- broadleaf evergreen temperate trees (BET-Te) and broadleaf vations is a better predictor than the model. The model deciduous trees (BDT), while gymnosperms tree species are parsimony was evaluated using the Akaike information criterion divided into needleleaf evergreen trees (NET) and needleleaf (AIC), which penalizes model overparameterization (Bozdogan, deciduous trees (NDT). Shrub species were divided into ever- 1987). We compared SOX AIC score with the b-function (Eqns green shrubs (ESh) and deciduous shrubs (DSh), and two grass 7, 8). The parameters f0 and Dcrit, (Eqn 8) were fitted to the PFT PFTs defined by their photosynthetic pathway (C3 and C4). The gs data, while hc and hw were held at their default values (0.033 grass PFTs and the NDT were excluded from the leaf-level evalu- and  1.5 MPa, respectively). ation because no stomatal conductance data were available for The uncertainty in plant hydraulic parameters caused by these PFTs in the dataset used in this study. within-PFT hydraulic variability was propagated to the model The plant hydraulic parameters used in SOX (i.e. Ψ50, a, and predictions using bootstrapped 95% confidence intervals. We rpmin) were fitted to the gs data using an algorithm that minimizes created the interval based on 1000 model runs with parameters the model residual sum of squares within the constraints of the resampled from the hydraulic trait data for each PFT. observed Ψ50, a and rpmin. We compiled hydraulic data for each PFT from the literature to constrain the leaf-level model fit. The Eddy-covariance based JULES-SOX evaluation Ψ50 for woody plants was obtained from a version of the Choat et al. (2012) dataset updated recently by Mencuccini et al. We evaluated default JULES and JULES-SOX against daily GPP (2019b). The shape parameter a of the xylem vulnerability func- estimates derived from eddy flux tower data at 62 FLUXNET tion (Eqn 2) was estimated from the linear gradient between Ψ50 sites (http://fluxnet.fluxdata.org, Baldocchi et al., 2001) and eight and the Ψ when the plant loses 88% of its maximum hydraulic LBA sites (https://daac.ornl.gov/LBA, Saleska et al., 2013), cover- conductance. The rpmin was estimated from branch-level ing all the major biomes of the world (Fig. S2; Table S2). In 10 hydraulic conductivity measurements scaled from branch to of these sites we also had data for surface (5–15 cm) soil moisture whole plant, taking into account plant height, Huber value and content, which was used to evaluate the model soil moisture xylem tapering using the calculations described in Christoffersen dynamics predictions. We classified the land cover on each site et al. (2016) and Savage et al. (2010) (Notes S4). All the data using the International Geosphere-Biosphere Programme (IGBP) used for these calculations were obtained from the hydraulic classification (Loveland et al., 2000). Each site was classified as dataset from Mencuccini et al. (2019b). We note that scaling one of the following categories according to its prescribed PFT branch to whole tree rpmin requires several assumptions about tree cover (Table S2): cropland (CRO), deciduous broadleaf forests hydraulic architecture (Notes S4). Therefore, the presented values (DBF), deciduous needleleaf forests (DNF), temperate evergreen of rpmin must be considered as a reference useful only to assess if broadleaf forests (EBF-Te), tropical evergreen broadleaf forests the rpmin input values used in the model are within the same (EBF-Tr), evergreen needleleaf forest (ENF), grassland (GRA), order of magnitude of the observations. The other parameters of mixed forest (MF), savannah (SAV), shrubland (SHR) and wet- lands (WET). We grouped the IGBP categories open and closed Table 1 Residual sum of squares (RSS), number of leaf-level stomatal shrublands into SHR, as we only had a single closed shrubland conductance observations (N) used to fit n parameters to the data, and the site. Similarly, woody savannah was grouped with SAV, as we resulting Akaike information criterion differences (DAIC) between stomatal only had two woody savannah sites. We divided the evergreen optimization based on xylem hydraulics (SOX) and the b-function. broad leaf forests category into EBF-Te and EBF-Tr, as these sites were dominated by distinct PFTs (BET-Te and BET-Tr, respec- SOX b-function tively). PFT N RSS n RSS n DAIC We evaluated JULES-SOX using the SOX hydraulic parame- ters (i.e. Ψ50, a, and rpmin) that minimized the residual sum of BET-Tr 434 4.83 3 6.53 2 126.1 BET-Te 1334 19.68 3 37.37 2 853.2 squares between SOX predictions and the eddy flux GPP obser- BDT 71 3.48 3 3.04 2 11.6 vations from a subset of the sites used for model evaluation NET 1571 0.65 3 2.29 2 1926.4 (Fig. S2; Table S2). Each site was used to calibrate the hydraulic ESh 133 3.37 3 7.94 2 112 parameters for its dominant PFT (i.e. the PFT covering > 50% DSh 64 2.76 3 8.03 2 66.4 of the site area), except for DSh, which was not dominant in any PFT, plant functional type; BET-Tr, broadleaf evergreen tropical tree; BET- of the available sites. We used a site with DSh cover of 35% (US- Te, broadleaf evergreen temperate tree; BDT, broadleaf deciduous tree; SRM) to calibrate the hydraulic parameters of this PFT. The NET, needleleaf evergreen tree; ESh, evergreen shrubs; DSh, deciduous hydraulic parameters of the others PFTs (if any) present on the shrubs. site were kept constant during the model runs for parameter Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1626 Research Phytologist calibration. Similar to the leaf-level evaluation, the parameter cal- complex: lower a increased gs sensitivity to less negative Ψpd, but ibration in JULES-SOX was constrained within the range of the decreased gs sensitivity to very negative Ψpd values (Fig. 1c). observed values of Ψ50, a, and rpmin for all PFTs, except for The patterns produced by the analytical SOX were similar to NDT, which did not have enough observations to satisfactorily the numerical version from Eller et al. (2018b), with a correlation constrain the model parameters. The Ψ50 for grasses was obtained coefficient ranging from 0.92 to 1 (Fig. S3). However, the use of from the Lens et al. (2016) dataset updated with data from linear gradients in Eqns 4 and 5 (Notes S2) can cause discrepan- Ocheltree et al. (2016). cies between the different model versions under certain ranges of environmental conditions. The analytical version of SOX under- estimated gs at low D (Fig. S3), overestimated gs at low ca, and gs Model setup increased faster in response to light (Fig. S3) than in the numeri- The JULES and JULES-SOX configurations used in this study cal model. employed the 10-layer canopy scheme with sunlit and shaded leaves in each layer as described in Clark et al. (2011). The SOX leaf-level evaluation canopy radiation profile was given by the two-stream approach from Sellers (1985), with the sun-fleck penetration scheme from Stomatal optimization based on xylem hydraulics simulated the Mercado et al. (2009), and an exponential decrease of photosyn- observed leaf-level gs responses to soil drought better than the b- thetic capacity through the canopy (Mercado et al., 2007). All function in all the studied woody PFTs, except BDT (Fig. 2). the model runs used in this study were site-level simulations The b-function predicted that all PFTs will reach gs = 0 at driven with hourly local meteorological data. Vegetation dynam- Ψpd > 2 MPa, whereas SOX predicted gs > 0 even when ics (Cox, 2001) was turned off and the site PFT coverage by site Ψpd < 4 MPa in some PFTs (Fig. 2b,e). The less conservative was prescribed based on the site vegetation description obtained stomatal behaviour predicted by SOX produced a NSE that was, from the site principal investigators (Table S2) or information on average, 0.65 higher and a RMSE that was 26% lower than from the site available on the FLUXNET website (https://fluxne the b-function. Most of the observed gs was within SOX 95% t.fluxdata.org/sites/site-list-and-pages/). Site soil hydraulic prop- confidence bounds derived from the hydraulic parameters’ uncer- erties were parameterized using Brooks & Corey (1964) relations. tainty (shaded region in Fig. 2). The only values outside SOX These properties were derived from data collected at each site or, uncertainty boundaries were the highest gs values in BET-Tr and when local data were not available, calculated from the sand/silt/ BET-Te (Fig. 2a,b), and the lowest NET gs values when clay fractions in the nearest gridbox in the high-resolution input Ψpd > 3.5 MPa (Fig. 2d). file to the Met Office Central Ancillary Program (Dharssi et al., Stomatal optimization based on xylem hydraulics produced a 2009), using approximations from Cosby et al. (1984). The better fit to the gs data, which resulted in a lower AIC than the b- model was spun up by recycling the meteorological data at each function for all PFTs, except BDT (Table 1). Fitting the two site for up to 50 yr. empirical parameters of the Jacobs (1994) equation (f0 and Dcrit; Eqn 8) to the gs data results in a b-function AIC score that is 512.1 higher than SOX (Table 1). For the BDT observations, Results the b-function results in an AIC score that is 11.6 lower than SOX. Our BDT observations were restricted to relatively well- SOX sensitivity to environmental and hydraulic drivers watered conditions (lowest Ψpd was  1.2 MPa), which limits the The SOX analytical approximation (Eqns 4, 5) has gs responses utility of this dataset to evaluate the model responses to soil to climate which are consistent with the patterns commonly drought. reported in the literature (Mott, 1988; Leuning, 1995; Dewar et al., 2018). The gs responses to Ipar and ca in SOX (Fig. 1a) are JULES-SOX site-level calibration given by the @A/@ci gradient decreasing at low light because of the changes in the light response curve, as A starts being limited The hydraulic parameters that maximized the JULES-SOX fit to by light (Notes S3), or at high ca (Notes S2). SOX correctly pre- the GPP data at the calibration sites (Table S2; Fig. S2) were dicted stomatal closure in response to increased ca under within 1 SD of the mean observed hydraulic parameters for most Rubisco-limited conditions (Mott, 1988; Fig. 1a). The classical PFTs (Table 2). The gymnosperm PFTs (NDT and NET) exponential gs responses to D (Leuning, 1995) was reproduced in required Ψ50 values 1.6 MPa less negative than their observed SOX (Fig. 1a) through the D effect on ξ (Eqn 5; Fig. S1a). An Ψ50 means to fit the GPP data, which is lower than the observed exponential gs decline was also predicted by SOX in response to SD range but still within the range of Ψ50 observations for NET decreasing Ψpd (Fig. 1a), which summarizes both the responses to (Ψ50 ranges from  2.3 to  7.5 MPa in NET). The NDT and the soil water availability in the root zone and the hydraulic stress BET-Tr calibrated a were also slightly lower than the SD range of transporting water to the top of the canopy (Eqn S1.2 in Notes (Table 2), but within the observed a range for BET-Tr (a ranges S1). The plant hydraulic parameters modulated the model sensi- from 1.8 to 7.8 in BET-Tr). The only PFT with a calibrated tivity to D or Ψpd (Fig. 1b–d), with a less negative Ψ50 or a higher rpmin outside the SD range of the mean rpmin was ESh (Table 2). rpmin increasing the gs sensitivity to Ψpd and D (Fig. 1c,d). The The monthly GPP modelled by JULES-SOX fitted the eddy effect of the vulnerability curve shape parameter a was more covariance GPP data better than the default JULES in eight out New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1627 (a) (b) 1.0 D a pd 50 r pmi n 0.8 I par ca 0.6 Fig. 1 (a, b) Stomatal conductance (gs) sensitivity to environmental drivers (a) and plant hydraulic traits (b) as modelled by 0.4 stomatal optimization based on xylem hydraulics (SOX) (D, vapour pressure deficit; 0.2 Ψpd, predawn water potential; Ipar, incident photosynthetically active radiation; ca, atmospheric CO2 partial pressure; Ψ50, Ψ 0.0 when plant loses 50% of its maximum conductance; a, shape of vulnerability 0 1 2 3 4 0 1 2 3 4 function; rpmin, minimum plant hydraulic D, pd , I par and c a relative to their default values 1 resistance). Variables were changed gs mol m 2 s individually while the others were held (c) (d) constant at their reference values 1.0 (D = 0.5 kPa, Ψpd = 0.5 MPa, Reference Ipar = 600 µmol m2 s1, ca = 36 Pa, a 2.5 0.8 Ψ50 = 2 MPa, a = 3, 50 1.5 MPa rpmin = 1 m2 s MPa mmol1). For (c) and (d) r pmi n 2 mmol 1 m2 s MPa the reference lines (dashed black) represent 0.6 values of Ψ50 = 3 MPa, a = 5, rpmin = 1 mmol1 m2 s MPa, and the coloured 0.4 lines show how changing each hydraulic parameter affects gs response to Ψpd and D. In (c) and (d), Ipar was set to 0.2 2000 µmol m2 s1. The Rubisco maximum carboxylation rate at 25°C (Vcmax25) was set 0.0 to 100 µmol m2 s and the rest of the photosynthetic parameters follow the 0 −1 −2 −3 −4 0 1 2 3 broadleaf evergreen tropical tree (BET-Tr) parameterization from Harper et al. (2016). pd MPa D kPa of the nine sites used for parameter calibration (Table S2; Figs 3, modelled Ψc at the calibration sites is within the interquartile S2). The default JULES NSE was 0.01 higher in the DSh site range of the observed minimum Ψc at midday for all woody (Fig. 3i), whereas in all the other sites JULES-SOX had a better PFTs, except NDT (Fig. 4). fit. The difference between JULES-SOX and default JULES NSE ranged from 0.03 for C3 grasses (Fig. 3f) to 11.44 for BET-Tr Biome-level JULES-SOX evaluation (Fig. 3a). The large improvement in the BET-Tr site was caused by the lower GPP decline predicted by SOX during January– Using JULES-SOX with calibrated SOX hydraulic parameters March and September–December. The decline in BET-Tr GPP produced a better fit to the GPP data than default JULES for 50 in default JULES can be attributed to the b-factor overestimating out of the 70 eddy flux evaluation sites (Tables 3, S2; Fig. 5). the effects of soil moisture on the vegetation carbon assimilation Across all biomes the JULES-SOX median NSE was 0.19 higher during drier periods (Fig. S4a). On average, JULES-SOX NSE than default JULES, and its RMSE was 19% lower (Table 3). for GPP was 1.59 higher than default JULES, while its RMSE The difference between JULES-SOX and JULES skill was highest was 38% lower than JULES. at EBF-Tr sites, which have a median NSE 3.18 higher and The less conservative stomatal behaviour predicted by SOX RMSE 45% lower in JULES-SOX (Table 3; Fig 5a). The fit of resulted in higher evapotranspiration rates throughout the year EBT-Te to data was also improved substantially by JULES-SOX, (Figs S5, S6), which depleted soil moisture to lower values than with JULES-SOX having a median NSE 1.01 higher and RMSE the b-function in default JULES during drier periods (Figs S4, 18% lower (Fig. 5a; Table 3). Default JULES only outperformed S7). The soil moisture dynamics from JULES-SOX are more JULES-SOX at CRO, which had a median NSE 0.08 lower in closely aligned with the monthly soil moisture observations in JULES-SOX, and GRA, where the RMSE 5% was higher in eight out of the 10 sites where soil moisture data were available JULES-SOX (Fig. 5a; Table 3). (Fig. S7). JULES-SOX NSE for monthly soil moisture was 1.67 Default JULES significantly underestimated the observed higher and RMSE was 19% lower than default JULES. JULES- mean annual GPP by 143.3 g C m2 across all biomes, which SOX also simulates realistic Ψc for most PFTs (Figs 4, S4). The corresponds to 13.6% of the observed mean annual GPP Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1628 Research Phytologist (a) (b) 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.0 0.0 Fig. 2 Predicted and observed (grey points) 0 −1 −2 −3 −4 −5 0 −1 −2 −3 −4 −5 stomatal conductance (gs) response to (c) (d) changes in leaf predawn water potential 0.25 1.0 (Ψpd) for the woody plant functional types (PFT) from Harper et al. (2016), except for 0.20 0.8 needleleaf deciduous trees, which were not present in the dataset used in this study. The 1 0.15 red and blue lines are the best fits from the 0.6 gs mol m 2 s stomatal optimization based on xylem hydraulics (SOX) and b-function (Eqns 7, 8), 0.10 0.4 respectively. The shaded regions are nonparametric 95% confidence boundaries 0.05 0.2 derived from 1000 bootstrapping replications of the SOX hydraulic inputs. All environmental conditions except Ψpd were 0.00 0.0 held constant at their median values when 0 −1 −2 −3 −4 −5 0 −1 −2 −3 −4 −5 the gs measurements were taken. The Ψpd (e) (f) was converted in soil volumetric water 1.0 1.0 content to drive the b-function using the Brooks & Corey (1964) equations 0.8 0.8 parameterized with soil physical properties derived from the Met Office Central Ancillary 0.6 0.6 Program (Dharssi et al., 2009). The model fit to data is shown as the root-mean-square errors (RMSE) and Nash-Sutcliffe (1970) 0.4 0.4 model efficiency index (NSE). The PFT abbreviations in each panel are as follows: (a) 0.2 0.2 broadleaf evergreen tropical tree (BET-Tr); (b) broadleaf evergreen temperate tree (BET- Te); (c) broadleaf deciduous tree (BDT); (d) 0.0 0.0 0 −1 −2 −3 −4 −5 0 −1 −2 −3 −4 −5 needleleaf evergreen tree (NET); (e) evergreen shrubs (ESh); and (f) deciduous pd MPa shrubs (DSh). (Fig. 5b). JULES-SOX deviation from the observed mean annual moisture and light availability controlling the vegetation GPP was considerably smaller (71.6 g C m2; Fig. 5b). The sig- GPP. The JULES LUE declined as soil dried out, with a nificantly lower annual GPP predicted by default JULES can be mean linear slope of 1.21 ( 0.1) across all biomes. By con- attributed to b-function-induced GPP declines, which also pro- trast, the JULES-SOX LUE–soil moisture relationship had a duced a stronger GPP seasonality than is present in the data mean slope of 0.73 ( 0.21), with some biomes, such as (Fig. 5c). JULES overestimated the median observed GPP sea- DBF, reaching a slope as low as 0.22 (Fig. 6b). The conse- sonality by 70%, compared with a 13% overestimation by quence of sustaining higher LUE at low soil moisture in JULES-SOX (Fig. 5c). This difference means JULES predicts JULES-SOX is a greater depletion of soil moisture, as indi- 17% of the sites have a markedly seasonal GPP with a Seasonality cated by the more left-skewed soil moisture probability dis- Index (SI; Walsh & Lawler, 1981) higher than 0.8, while just 4% tribution predicted by JULES-SOX (lower panels in Fig. 6). of the sites actually have SI > 0.8. JULES-SOX predicts only 8% The mean moisture content of the top 1 m of soil predicted of the sites would have SI > 0.8. by JULES-SOX was, on average, 10% lower than default The light-use efficiency (LUE; Fig. 6) is the ratio between JULES. In JULES-SOX some biomes, such as ENF, could GPP and the Ipar absorbed by the canopy (Stocker et al., reach a soil moisture, on average, 17% lower than JULES 2018), and can be used to disentangle the effects of soil (Fig. 6f). New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1629 Table 2 Observed (Obs) mean ( SD) hydraulic parameters compiled from the literature for each plant functional type (PFT) from JULES (Harper et al., 2016). Ψ50 (MPa) a (unitless) rpmin (103 mol1 m2 s MPa) PFT N Obs Cala Nb Obs Cal N Obs Cal BET-Tr 77 1.9 ( 1.3) 1.7 20 4.4 ( 2.1) 2.1 40 2.2 ( 3.4) 0.6 BET-Te 44 2.7 ( 1.5) 1.8 17 3.7 ( 1.8) 2.8 40 3.1 ( 8) 5 BDT 87 2.6 ( 1.4) 1.6 43 5.5 ( 3.8) 3.5 31 5.3 ( 5.6) 0.5 NET 48 4.2 ( 1.2) 2.6 25 8.7 ( 4.9) 4.9 20 2.4 ( 1.8) 4.2 NDT 5 3.4 ( 0.6) 1.8 2 7.4 ( 5) 1.8 2 8 ( 4.3) 9 C3 45 3.1 ( 1.6) 2.4 – – 2.2 – – 3.2 C4 15 2.7 ( 1.7) 1.5 – – 1.8 – – 9.5 ESh 61 4 ( 2.2) 2.1 53 4.1 ( 3.3) 2.5 49 1.5 ( 1.8) 9.5 DSh 26 4 ( 2.3) 1.8 3 3.4 ( 2.2) 2.1 4 2.6 ( 2.4) 5 BET-Tr, broadleaf evergreen tropical tree; BET-Te, broadleaf evergreen temperate tree; BDT, broadleaf deciduous tree; NET, needleleaf evergreen tree; NDT, needleleaf deciduous tree; C3, C3 grasses; C4, C4 grasses; ESh, evergreen shrubs; DSh, deciduous shrubs. a The calibrated (Cal) columns are the parameter values that maximize the fit of the Joint UK Land Environment Simulator–stomatal optimization based on xylem hydraulics (JULES-SOX) to observed gross primary productivity (GPP) in the calibration sites (see Supporting Information Table S2; Fig. S2). b The N column is the number of species compiled for the correspondent parameter. recently, Wang et al. (2019) shows that a similar hydraulics-based Discussion model can predict stomatal responses to increased CO2 better We report the first evaluation of a LSM using a stomatal opti- than the Ball–Berry–Leuning empirical model (Leuning, 1995). mization model fully based on xylem hydraulics to drive the vege- These results show the potential of using plant hydraulics to tation stomatal responses to climate. Our results provide support model the stomatal behaviour of plants across contrasting envi- for the SOX underlying hypothesis that stomata evolved to bal- ronmental conditions, and supports its use in ESMs to project ance carbon assimilation with instantaneous hydraulic conduc- the evolution of global climate. tance loss. The risk of mortality through hydraulic failure (Choat The analytical formulation developed for SOX facilitates its et al., 2012; Rowland et al., 2015; Anderegg et al., 2016; Adams coupling to LSMs, allowing the host LSM to constrain its predic- et al., 2017) should drive the evolution of mechanisms to prevent tions using plant hydraulic information. We show that inclusion the plant from reaching lethal embolism thresholds (Sperry, of plant hydraulics in JULES through SOX improves its capabili- 2004). There is abundant evidence that stomata controls xylem ties to simulate GPP and soil moisture dynamics in most of the tension, and consequently embolism (Hubbard et al., 2001; Bro- studied biomes (Figs 3–5). In addition, SOX opens new possibili- dribb et al., 2003; Meinzer et al., 2009; Klein, 2014). Our model ties to evaluate LSM predictions and expands the range of represents this xylem–stomata coordination through the assump- hypotheses that can be tested with JULES. Using JULES-SOX tion of optimization by natural selection (Wolf et al., 2016). within ESMs will allow us to understand how hydraulic processes Whereas our model fits the observations of most PFTs better affect climatic and biogeochemical cycles at the global scale, as than its empirical alternative, there is still a considerable amount well as to investigate the role of plant hydraulics on vegetation of unexplained variance in the data (Fig. 2). This can be partially distribution and its response to climate change. attributed to the large hydraulic heterogeneity within each PFT, but we must also acknowledge that many processes not directly SOX parametrization and parsimony related to xylem hydraulics are important to plant life history and stomatal evolution. Processes related to nutrient use and acquisi- Other LSMs and DGVMs have already successfully employed tion, carbohydrate allocation and storage, the maintenance of tis- principles of plant hydraulics (Hickler et al., 2006; Bonan et al., sues and biochemical apparatus, and protection from pathogens 2014; Kennedy et al., 2019), but JULES-SOX is the first LSM to and herbivores (Melotto et al., 2008; Cramer et al., 2009; Pren- use the new generation of hydraulically based stomatal optimiza- tice et al., 2014) could all explain part of our model residual vari- tion models (Wolf et al., 2016; Sperry et al., 2017; Anderegg ance. It is extremely important to explore the relevance of these et al., 2018b; Eller et al., 2018b) to predict stomatal responses to processes in future research on stomatal optimality. However, the climate. The SPA (Williams et al., 1996) adaptation to the com- SOX model as we propose it already provides a parsimonious munity land model (CLM) by Bonan et al. (2014) was one of the alternative to the empirical models commonly used in LSMs. first approaches to link plant stomatal function to plant hydraulic Our findings that xylem hydraulics-based models can ade- processes in a LSM. Despite SPA being an extremely useful quately simulate stomatal behaviour agree with other recent stud- model, SOX has an advantage in circumstances where assuming a ies. For example, Anderegg et al. (2018b) shows that a strict isohydric behaviour is not appropriate (Klein, 2014; hydraulics-based optimization model can simulate the stomatal Martinez-Vilalta et al., 2014). In relation to SOX, SPA does not behaviour of woody plants better than the CF model. More represent dynamic changes in the plant hydraulic conductance or Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1630 Research Phytologist LBA_K67(BETTr) (a) FR_Pue(BETTe) (b) US_MMS(BDT) (c) SOX RMSE 43.74 NSE −1.89 SOX RMSE 22.4 NSE 0.64 SOX RMSE 31.08 NSE 0.89 RMSE 97.39 NSE −13.33 RMSE 43.64 NSE −0.36 RMSE 41.27 NSE 0.81 300 200 100 0 IT_Ren(NET) (d) RU_SkP(NDT) (e) CN_HaM(C3) (f) SOX RMSE 19.35 NSE 0.94 SOX RMSE 30.28 NSE 0.47 SOX RMSE 33.18 NSE 0.76 RMSE 60.27 NSE 0.38 RMSE 38.37 NSE 0.14 RMSE 34.66 NSE 0.73 300 GPP g C m 2 per month 200 100 0 SD_Dem(C4) (g) IT_Noe(ESh) (h) US_SRM(DSh) (i) SOX RMSE 52.48 NSE 0.09 SOX RMSE 16.11 NSE 0.79 SOX RMSE 11.57 NSE 0.82 RMSE 69.04 NSE −0.57 RMSE 23.68 NSE 0.55 RMSE 11.26 NSE 0.83 300 200 100 0 Jan Feb Mar Apr Jun Jun Jul May Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr Jun May Jul Aug Sep Oct Nov Dec Jan Feb Mar Apr May Aug Sep Oct Nov Dec Month Fig. 3 Monthly mean gross primary production (GPP) modelled by default Joint UK Land Environment Simulator (JULES, blue line) and JULES-stomatal optimization based on xylem hydraulics (JULES-SOX, red line) vs observations (grey points are means and bars are 2 9 SE) at each eddy flux site used for calibrating the SOX hydraulic parameters (plant functional type (PFT); Supporting Information Table S2; Fig. S3). The model fit to data is shown as the root-mean-square errors (RMSE) and Nash-Sutcliffe (1970) model efficiency index (NSE). The PFT abbreviations in each panel are as follows: (a) broadleaf evergreen tropical tree (BET-Tr); (b) broadleaf evergreen temperate tree (BET-Te); (c) broadleaf deciduous tree (BDT); (d) needleleaf evergreen tree (NET); (e) needleleaf deciduous tree (NDT); (f) C3 grasses (C3); (g) C4 grasses (C4); (h) evergreen shrubs (ESh); (i) deciduous shrubs (DSh). New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1631 0 −2 Fig. 4 Minimum observed midday leaf water potential (Ψmidday) from 279 woody plant species compiled from the literature grouped using the Harper et al. (2016) plant −4 MPa functional type (PFT) categories. The Ψmidday for each of the calibration sites as modelled midday by stomatal optimization based on xylem hydraulics (SOX) (see Supporting −6 Information Table S2; Fig. S2) is plotted in red. The circle is the mean Ψmidday and the arrows indicate the minimum and maximum Ψmidday. The data for the deciduous PFT were −8 restricted to the growing season. The PFT abbreviations in each panel are as follows: broadleaf evergreen tropical tree (BET-Tr); broadleaf evergreen temperate tree (BET- −10 Observations Te); broadleaf deciduous tree (BDT); SOX needleleaf evergreen tree (NET); needleleaf BETTr BETTe BDT NET NDT ESh DSh deciduous tree (NDT); C3 grasses (C3); (g) C4 C3 C4 grasses (C4); evergreen shrubs (ESh); deciduous shrubs (DSh). Table 3 Median Nash-Sutcliffe (1970) model efficiency index (NSE) and changes in hydraulic conductance in different compartments root-mean-square error (RMSE) for the biomes used for evaluating the along the soil–atmosphere continuum, providing a more detailed Joint UK Land Environment Simulator–stomatal optimization based on representation than SOX of hydraulic processes occurring along xylem hydraulics (JULES-SOX) and the default JULES. the soil–plant hydraulic pathway. However, PHS still requires JULES-SOX JULES empirical parameters to represent stomatal responses to soil drought and D (Kennedy et al., 2019), namely the g0 and g1 a b Biome N NSE RMSE NSE RMSE parameters from the Medlyn et al. (2011) model, and the critical CRO 3 0.49 123.12 0.57 141.1 and wilting points used in the empirical stress factor. The main DBF 7 0.89 37.32 0.83 47.19 advantage of SOX is providing an alternative to the b-function DNF 1 0.58 25.93 0.37 31.97 and empirical stomatal parameters by linking plant hydraulic EBF-Te 3 0.23 45.22 1.24 66.36 processes directly to stomatal functioning. As we treat the soil– EBF-Tr 6 0.41 40.36 2.77 73.53 plant–atmosphere pathway as a single hydraulic compartment, ENF 5 0.9 34.14 0.59 40.58 GRA 12 0.22 32.31 0.01 30.62 SOX only requires the hydraulic parameters rpmin, Ψ50 and a to MF 3 0.85 47.87 0.59 79.29 predict stomatal responses to climate. This makes SOX even SAV 5 0.4 59.72 2.12 89.69 more parsimonious than default JULES, which requires four SHR 4 0.78 14.90 0.64 15.92 empirical parameters to simulate stomatal responses to climate WET 21 0.68 32.23 0.46 38.67 (Eqns 7, 8) and does not simulate any aspect of vegetation a Biome abbreviations are as follows: CRO, cropland; DBF, deciduous hydraulic functioning (Clark et al., 2011). broadleaf forests; DNF, deciduous needleleaf forests; EBF-Te, temperate We show that the SOX hydraulic parameters in most PFTs evergreen broadleaf forests; EBF-Tr, tropical evergreen broadleaf forests; can be constrained with plant branch-level hydraulic observations ENF, evergreen needleleaf forest; GRA, grassland; MF, mixed forest; SAV, (Table 2), which is an advantage over models that employ empir- savannah; SHR, shrubland; WET, wetlands. b The N column is the number of sites representing the biome in the eddy ical parameters difficult to constrain and interpret biologically. covariance dataset. However, we observed discrepancies between the SOX-calibrated parameters and the observed hydraulic traits in certain PFTs (Table 2). In some cases, such as NDT, the parameter discrep- an anisohydric mode of stomatal regulation (Williams et al., ancy may have been a result of a very restricted observational 1996; Fisher et al., 2006). However, SPA accounts for plant sampling of hydraulic parameters in this group. The NDT only hydraulic capacitance, which can be important for plant func- had Ψ50 data for five species and a and rpmin for two species tioning, especially during the early morning (Goldstein et al., (Table 2). Considering that the observations used in this study 1998), and is currently not implemented in SOX. were not collected in the same FLUXNET sites used to evaluate Recently, Kennedy et al. (2019) implemented a plant SOX, some of the observed discrepancies between calibrated and hydraulic scheme (PHS) in a CLM. The PHS simulates dynamic measured parameters might reflect hydraulic differences between Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1632 Research Phytologist 1 − CRO 5 − EBF−Tr 9 − SAV 2 − DBF 6 − ENF 10 − SHR 3 − DNF 7 − GRA 11 − WET SOX (a) 0.1 0.2 4 − EBF−Te 8 − MF (b) 1.5 0.3 1.2 0.4 Obs GPP 0.8 0.5 9 C 1.5 0.4 0.6 ( kg m 2 yr 1( or re Standard deviation (normalized) la Model GPP 0.7 tio 9 n −1.2 −0.8 −0.4 1.0 1.0 4 5 0.8 11 (c) 7 0.9 1.0 11 3 Obs GPP SI 7 2 0.5 5 2 6 0.95 10 6 0.5 0.5 3 4 8 10 RM Model GPP SI 1 8 SE 1 0.99 0.0 0.0 −0.5 ALL CRO DBF DNF ENF GRA MF SHR WET SAV EBF−Te EBF−Tr 0.0 0.5 1.0 1.5 Fig. 5 (a) The Taylor diagram shows the difference in Joint UK Land Environment Simulator (JULES) and JULES-stomatal optimization based on xylem hydraulics (JULES-SOX) skill to predict the monthly gross primary productivity (GPP) in each studied biome. Green lines are the model-centred root-mean- square errors (RMSE), and points closer to the reference circle on the x-axis indicate higher model skill. The two arrows highlight the improvement in model skill for tropical evergreen broadleaf forests (EBF-Tr) and temperate evergreen broadleaf forests (EBF-Te). The boxplot panels show the differences between models (default JULES in blue, JULES-SOX in red) and observations (Obs) in the annual GPP (b) and the GPP seasonality (GPP SI) (c). Data gaps were excluded from the annual GPP calculations for both models and observations, and therefore the differences can be used to evaluate the model skill, but the absolute values do not represent the total annual GPP in each biome. The GPP SI was computed using the approach from Walsh & Lawler (1981). Boxes filled with lines are different (at a = 0.05) from 0 in a one-sample t-test. The biomes are abbreviated as follows: cropland (CRO); deciduous broadleaf forests (DBF); deciduous needleleaf forests (DNF); temperate evergreen broadleaf forests (EBF-Te); tropical evergreen broadleaf forests (EBF-Tr); evergreen needleleaf forest (ENF); grassland (GRA); mixed forest (MF); savannah (SAV); shrubland (SHR); and wetlands (WET). populations treated as the same PFT in this study. For example, hydraulic conductance and, ultimately, control its water use the deciduous angiosperms species present in the XFT dataset (Fisher et al., 2007). These bottlenecks could bias the SOX-cali- used in this study contain mostly hydraulic data from cold-decid- brated hydraulic parameters towards the limiting component and uous temperate species (Mencuccini et al., 2019b), which might explain its departure from the branch-level hydraulic data. In this not be adequate to describe the hydraulic system of tropical and case, SOX parameterization would benefit from the use of more subtropical drought-deciduous. Our hydraulic scheme opens up integrative methodologies to estimate hydraulic parameters that possibilities of improving the representation of different global represent the entire soil–plant hydraulic vulnerability (Eller et al., vegetation types in JULES with different hydraulic and pheno- 2018a). Alternatively, the SOX structure (i.e. the K function in logical strategies. Capturing the large diversity of ecological Eqn 2) would need to explicitly represent the variability between strategies in plants is important to simulate species-rich ecosys- different hydraulic compartments along the soil–plant–atmo- tems such as tropical forests (Xu et al., 2016). sphere pathway, similar to SPA or other models (Eller et al., Anderegg et al. (2018a) computed the community-weighted 2018b; Kennedy et al., 2019; Mencuccini et al., 2019b). average values for Ψ50 in two of the FLUXNET sites used in this study (US-MMS and IT-Ren) and obtained values closer to the Ecosystem-level implications of SOX calibrated values for BDT and NET (2.1 and 3.6 MPa, respectively) than the means from our compiled hydraulic dataset Stomatal optimization based on xylem hydraulics improved (Table 2). In Eller et al. (2018b) a numerical version of SOX out- JULES GPP simulation in over 70% of the 70 studied sites, and performed the b-function approach when parameterized with soil moisture dynamics in 80% of the 10 sites where soil moisture locally measured branch-level hydraulic data from EBF-Tr. These data were available. This improved fit was achieved using findings suggest that SOX can be constrained with in situ hydraulic parameters calibrated against the GPP data of a small hydraulic measurements when these are available. However, we subset of eddy flux sites (the sites in Fig. S2), which suggests that must also consider the possibility that there are intrinsic limita- the calibrated parameters are generic enough to be used in global tions in using branch-level hydraulic data to parameterize the simulations. The lower sensitivity of SOX to soil moisture model. Roots and leaves can be more vulnerable to embolism improved the simulations of annual GPP (Fig. 5) and predicted than branches (Bartlett et al., 2016; Wolfe et al., 2016), which terrestrial biomes to assimilate on average 2.58 Mg C ha1 yr1 can make these tissues bottleneck plant hydraulic conductance or 30% more than predicted by default JULES. This increased during drought. The soil outside the roots can also limit plant carbon assimilation could affect Earth’s atmospheric CO2 New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1633 SOX β 1 (a) (b) (c) 0.5 y −17.1+25.5x y 3.1+2.4x y −1.5+7.9x 0.5 y −6.4+9.8x y −1+4.2x y −2.4+5.8x 0 0 2 4 CRO DBF DNF 6 1 (d) (e) (f) 0.5 y 0.5+5.2x y 0.4+7.3x y 1.7+5.8x Light use efficiency normalised 0.5 y −4.2+7.7x y −4.1+8.2x y −2.7+6.6x 0 0 Probability density 2 4 EBF−Te EBF−Tr ENF 6 1 Fig. 6 Model predictions of the normalized (g) (h) (i) 0.5 light-use efficiency responses to soil y 0.2+7.2x y −0.8+7.9x y 1.9+5.1x moisture, expressed as a fraction of the soil moisture saturation point at the top 1 m of 0.5 soil. The light use efficiency is the ratio y −2.3+6.2x y −4.8+8.6x y −2.1+5.7x between gross primary productivity and the 0 0 photosynthetic active radiation absorbed by 2 the canopy. The default JULES predictions 4 are in blue and JULES-SOX predictions in red. GRA MF SAV 6 The lines in the scatter plot panels are linear 1 (j) (k) 0.2 0.4 0.6 0.8 1.0 regressions fit to the data. The histograms on 0.5 the bottom panels are the soil moisture y −0.3+7.7x y −0.8+7.8x probability density predicted by each model. The biomes are abbreviated as follows: (a) 0.5 cropland (CRO); (b) deciduous broadleaf y −2.5+6.2x y −2.5+5.8x forests (DBF); (c) deciduous needleleaf 0 0 forests (DNF); (d) temperate evergreen 2 broadleaf forests (EBF-Te); (e) tropical 4 evergreen broadleaf forests (EBF-Tr); (f) SHR WET 6 evergreen needleleaf forest (ENF); (g) 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 grassland (GRA); (h) mixed forest (MF); (i) savannah (SAV); (j) shrubland (SHR); and (k) wetlands (WET). Fraction of soil moisture evolution and climate change projections (Cox et al., 2000; Win- et al., 2010; Meir et al., 2018). da Costa et al. (2018) showed that kler et al., 2019). even after 15 yr of partial rainfall exclusion, Amazon trees can The JULES-SOX model particularly improved the fit of EBF- maintain or even increase their transpiration rates (albeit follow- Tr sites to the observations (Fig. 5; Table 3), using hydraulic ing significant mortality). Whereas tropical forest resistance to parameters very similar to those observed in BET-Tr (Table 2). drought has previously been attributed only to deep roots pos- Considering that SOX is also able to capture the response of sessed by the vegetation (Nepstad et al., 1994), our results indi- EBF-Tr even to extreme experimental drought (Eller et al., cate that plants more resistant to embolism could maintain their 2018b), JULES-SOX may contribute to decrease the large uncer- carbon assimilation during drought even without a deeper root tainty in how these important ecosystems will respond to climate system. change (Sitch et al., 2008). Tropical forest productivity estimated The unavoidable consequence of maintaining stomatal gas by SOX is less sensitive to seasonal soil drought (Figs 3, S4), exchange during soil drought is a greater depletion of soil mois- which is consistent with the little seasonality often observed in ture reserves (Figs 6, S4, S7). This behaviour is a direct conse- tropical forest–atmosphere CO2 exchange (Grace et al., 1995; quence of the main assumption in SOX, which reflects a ‘use or Carswell et al., 2002; Alden et al., 2016), as well as to forest lose it’ stomatal regulation strategy with respect to soil moisture responses to experimental drought (Meir et al., 2009; da Costa (Sperry et al., 2017). SOX assumes plants with a more Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1634 Research Phytologist conservative water-use strategy will be outcompeted by neigh- Maurizio Mencuccini https://orcid.org/0000-0003-0840- bouring plants with a less conservative stomatal behaviour (Wolf 1477 et al., 2016). The demographic consequences of the stomatal reg- Rafael S. Oliveira https://orcid.org/0000-0002-6392-2526 ulation strategy embedded in SOX should be explored in future Bruno H. P. Rosado https://orcid.org/0000-0002-8924-8672 studies using the dynamic vegetation component of JULES (Cox, Teresa Rosas https://orcid.org/0000-0002-8734-9752 2001; Moore et al., 2018). The more competitive soil moisture Lucy Rowland https://orcid.org/0000-0002-0774-3216 dynamics predicted by SOX, together with a more accurate repre- Stephen Sitch https://orcid.org/0000-0003-1821-8561 sentation of vegetation drought-induced mortality, which also Grazielle S. Teodoro https://orcid.org/0000-0002-5528-8828 can be developed from SOX, might be the key to predicting sud- Yael Wagner https://orcid.org/0000-0002-2588-9278 den and widespread vegetation die-off during droughts that have Karina Williams https://orcid.org/0000-0002-1185-535X been increasingly reported in ecosystems around the globe (Allen Georg Wohlfahrt https://orcid.org/0000-0003-3080-6702 et al., 2010; Worrall et al., 2010; Meir et al., 2015). Acknowledgements References Adams HD, Zeppel MJB, Anderegg WRL, Hartmann H, Landh€a usser SM, This research was supported by the Newton Fund through the Tissue DT, Huxman TE, Hudson PJ, Franz TE, Allen CD et al. 2017. A Met Office Climate Science for Service Partnership Brazil (CSSP multi-species synthesis of physiological mechanisms in drought-induced tree Brazil) and an NERC independent fellowship grant NE/ mortality. Nature Ecology and Evolution 1: 1285–1291. N014022/1 to LR. This work used eddy covariance data acquired Alden CB, Miller JB, Gatti LV, Gloor MM, Guan K, Michalak AM, van der and shared by the FLUXNET community, including the follow- Laan-Luijkx IT, Touma D, Andrews A, Basso LS et al. 2016. Regional atmospheric CO2 inversion reveals seasonal and geographic differences in ing networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, Amazon net biome exchange. Global change biology 22: 3427–3443. CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet- Allen CD, Macalady AK, Chenchouni H, Bachelet D, McDowell N, Vennetier Canada, GreenGrass, ICOS, KoFlux, LBA, NECC, OzFlux- M, Kitzberger T, Rigling A, Breshears DD, Hogg EH et al. 2010. A global TERN, TCOS-Siberia, and USCCC. The ERA-Interim reanaly- overview of drought and heat-induced tree mortality reveals emerging climate sis data were provided by ECMWF and processed by LSCE. The change risks for forests. Forest Ecology and Management 259: 660–684. Anderegg WRL, Klein T, Bartlett M, Sack L, Pellegrini AFA, Choat B, Jansen FLUXNET eddy covariance data processing and harmonization S. 2016. Meta-analysis reveals that hydraulic traits explain cross-species were carried out by the European Fluxes Database Cluster, patterns of drought-induced tree mortality across the globe. Proceedings of the AmeriFlux Management Project, and Fluxdata project of National Academy of Sciences, USA 113: 5024–5029. FLUXNET, with the support of CDIAC and ICOS Ecosystem Anderegg WRL, Konings AG, Trugman AT, Yu K, Bowling DR, Gabbitas R, Thematic Center, and the OzFlux, ChinaFlux and AsiaFlux Karp DS, Pacala S, Sperry JS, Sulman BN et al. 2018a. Hydraulic diversity of forests regulates ecosystem resilience during drought. Nature 561: 538. offices. Anderegg WRL, Wolf A, Arango-Velez A, Choat B, Chmura DJ, Jansen S, Kolb T, Li S, Meinzer FC, Pita P et al. 2018b. Woody plants optimise stomatal behaviour relative to hydraulic risk. Ecology Letters 21: 968–977. Author contributions Baldocchi D, Falge E, Gu L, Olson R, Hollinger D, Running S, Anthoni P, Bernhofer C, Davis K, Evans R et al. 2001. FLUXNET: a new tool to study CBE, LR, MM, SS and PMC led the scientific development of the temporal and spatial variability of ecosystem-scale carbon dioxide, water SOX. PMC and CBE derived the analytical solution. CBE evalu- vapor, and energy flux densities. Bulletin of the American Meteorological Society ated leaf-level SOX using data provided by LR, PM, MM, TR, 82: 2415–2434. BEM, YW, TK, GST, RSO, ISM, BHPR. CBE and KW coded Bartlett MK, Klein T, Jansen S, Choat B, Sack L. 2016. The correlations and sequence of plant stomatal, hydraulic, and wilting responses to drought. SOX into JULES. KW and AH created a JULES suite used by Proceedings of the National Academy of Sciences, USA 113: 13098–13103. CBE to evaluate JULES-SOX against eddy covariance data col- Berninger F, Hari P. 1993. Optimal regulation of gas exchange: evidence from lected by KF, GW, LM, among other FLUXNET and LBA PIs. field data. Annals of Botany 71: 135–140. All authors contributed to writing the manuscript. Best MJ, Pryor M, Clark DB, Rooney GG, Essery R, Menard CB, Edwards JM, Hendry MA, Porson A, Gedney N. 2011. The Joint UK Land Environment Simulator (JULES), model description–Part 1: energy and water fluxes. ORCID Geoscientific Model Development 4: 677–699. Bonan GB, Williams M, Fisher RA, Oleson KW. 2014. Modeling stomatal Peter M. Cox https://orcid.org/0000-0002-0679-2219 conductance in the earth system: linking leaf water-use efficiency and water Cleiton B. Eller https://orcid.org/0000-0002-7795-2574 transport along the soil-plant-atmosphere continuum. Geoscientific Model Development 7: 2193–2222. Kathrin Fuchs https://orcid.org/0000-0003-1776-283X Bozdogan H. 1987. Model selection and Akaike’s Information Criterion (AIC): Anna Harper https://orcid.org/0000-0001-7294-6039 the general theory and its analytical extensions. Psychometrika 52: 345–370. Tamir Klein https://orcid.org/0000-0002-3882-8845 Brodribb TJ, Holbrook NM, Edwards EJ, Gutierrez MV. 2003. Relations Ilaine S. Matos https://orcid.org/0000-0001-5557-5133 between stomatal closure, leaf turgor and xylem vulnerability in eight tropical Belinda E. Medlyn https://orcid.org/0000-0001-5728-9827 dry forest trees. Plant, Cell & Environment 26: 443–450. Brooks R, Corey T. 1964. Hydraulic properties of porous media. Hydrology Patrick Meir https://orcid.org/0000-0002-2362-0398 Papers, Colorado State University 24: 37. Leonardo Montagnani https://orcid.org/0000-0003-4412- Buckley TN. 2017. Modeling stomatal conductance. Plant Physiology 174: 572– 4243 582. New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1635 Buckley TN. 2019. How do stomata respond to water status? New Phytologist Eller CB, Bittencourt PRL, Oliveira RS. 2018a. Using sap flow to measure 224: 21–36. whole-tree hydraulic conductance loss in response to drought. Acta Carswell FE, Costa AL, Palheta M, Malhi Y, Meir P, Costa JDPR, Ruivo MDL, Horticulturae 1222: 75–84. Leal LDSM, Costa JMN, Clement RJ et al. 2002. Seasonality in CO2 and Eller CB, Rowland L, Oliveira RS, Bittencourt PRL, Barros FV, da Costa ACL, H2O flux at an eastern Amazonian rain forest. Journal of Geophysical Research Meir P, Friend AD, Mencuccini M, Sitch S et al. 2018b. Modelling tropical D: Atmospheres 107: LBA-43. forest responses to drought and El Ni~ no with a stomatal optimization model Choat B, Jansen S, Brodribb TJ, Cochard H, Delzon S, Bhaskar R, Bucci SJ, based on xylem hydraulics. Philosophical transactions of the Royal Society of Feild TS, Gleason SM, Hacke UG et al. 2012. Global convergence in the London. Series B: Biological sciences 373: 20170315. vulnerability of forests to drought. Nature 491: 752–755. Farquhar G, Schulze E, Kuppers M. 1980. Responses to humidity by stomata of Christoffersen BO, Gloor M, Fauset S, Fyllas NM, Galbraith DR, Baker TR, Nicotiana glauca L. and Corylus avellana L. are consistent with the optimization Kruijt B, Rowland L, Fisher RA, Binks OJ et al. 2016. Linking hydraulic traits of carbon dioxide uptake with respect to water loss. Australian Journal of Plant to tropical forest function in a size-structured and trait-driven model (TFS v.1- Physiology 7: 315–327. Hydro). Geoscientific Model Development 9: 1–29. Fisher RA, Williams M, Da Costa AL, Malhi Y, Da Costa RF, Almeida S, Meir Clark DB, Mercado LM, Sitch S, Jones CD, Gedney N, Best MJ, Pryor M, P. 2007. The response of an Eastern Amazonian rain forest to drought stress: Rooney GG, Essery RLH, Blyth E et al. 2011. The Joint UK Land results and modelling analyses from a throughfall exclusion experiment. Global Environment Simulator (JULES), model description – Part 2: carbon fluxes Change Biology 13: 2361–2378. and vegetation. Geoscientific Model Development Discussions 4: 641–688. Fisher RA, Williams M, Do Vale RL, Da Costa AL, Meir P. 2006. Evidence Collatz GJ, Ball JT, Grivet C, Berry JA. 1991. Physiological and environmental from Amazonian forests is consistent with isohydric control of leaf water regulation of stomatal conductance, photosynthesis and transpiration: a model potential. Plant, Cell & Environment 20: 151–165. that includes a laminar boundary layer. Agricultural and Forest meteorology 54: Friend AD. 1995. PGEN: an integrated model of leaf photosynthesis, and 107–136. conductance. Ecological Modelling 77: 233–255. Collatz GJ, Ribas-Carbo M, Berry JA. 1992. Coupled photosynthesis-stomatal Friend AD, Cox PM. 1995. Modelling the effects of atmospheric CO2 on conductance model for leaves of C4 plants. Functional Plant Biology. 19: 519– vegetation–atmosphere interactions. Agricultural and Forest Meteorology 73: 538. 295. Cosby BJ, Hornberger GM, Clapp RB, Ginn TR. 1984. A statistical exploration Goldstein G, Andrade JL, Meinzer FC, Holbrook NM, Cavelier J, Jackson P, of the relationships of soil moisture characteristics to the physical properties of Celis A. 1998. Stem water storage and diurnal patterns of water use in tropical soils. Water Resources Research 20: 682–690. forest canopy trees. Plant, Cell & Environment 21: 397–406. Cowan IR. 1986. Economics of carbon fixation in higher plants. In: Givnish TJ, Grace J, Lloyd J, McIntyre J, Miranda AC, Meir P, Miranda HS, Nobre C, ed. On the economy of plant form and function. Cambridge, UK: Cambridge Moncrieff J, Massheder J, Malhi Y et al. 1995. Carbon dioxide uptake by an University Press, 133–170. undisturbed tropical rain forest in Southwest Amazonia, 1992 to 1993. Science Cowan I. 2002. Fit, fitter, fittest; where does optimisation fit in? Silva Fennica 36: 270: 778–780. 745–754. Harper AB, Cox PM, Friedlingstein P, Wiltshire AJ, Jones CD, Sitch S, Cowan IR, Farquhar GD. 1977. Stomatal function in relation to leaf metabolism Mercado LM, Groenendijk M, Robertson E, Kattge J et al. 2016. Improved and environment. Symposia of the Society for Experimental Biology 31: 471–505. representation of plant functional types and physiology in the Joint UK Land Cox PM. 2001. Description of the "TRIFFID" Dynamic Global Vegetation Model. Environment Simulator (JULES v4.2) using plant trait information. Hadley Centre technical note 24. Bracknell, UK: Met Office. Geoscientific Model Development 9: 2415–2440. Cox PM, Betts RA, Jones CD, Spall SA, Totterdell IJ. 2000. Acceleration of Hickler T, Prentice IC, Smith B, Sykes MT, Zaehle S. 2006. Implementing global warming due to carbon-cycle feedbacks in a coupled climate model. plant hydraulic architecture within the LPJ Dynamic Global Vegetation Nature 408: 184. Model. Global Ecology and Biogeography 15: 567–577. Cox PM, Huntingford C, Harding RJ. 1998. A canopy conductance and Hubbard RM, Ryan MG, Stiller V, Sperry JS. 2001. Stomatal conductance and photosynthesis model for use in a GCM land surface scheme. Journal of photosynthesis vary linearly with plant hydraulic conductance in ponderosa Hydrology 212–213: 79–94. pine. Plant, Cell & Environment 24: 113–121. Cramer MD, Hawkins HJ, Verboom GA. 2009. The importance of nutritional Jacobs C. 1994. Direct impact of atmospheric CO2 enrichment on regional regulation of plant water flux. Oecologia 161: 15–24. transpiration. PhD thesis, Wageningen Agricultural University, Wageningen, da Costa ACL, Galbraith D, Almeida S, Portela BTT, da Costa M, de Athaydes the Netherlands. Silva Junior J, Braga AP, de Goncßalves PHL, de Oliveira AA, Fisher R et al. Kattge J, Dıaz S, Lavorel S, Prentice IC, Leadley P, B€onisch G, Garnier E, 2010. Effect of 7 yr of experimental drought on vegetation dynamics and Westoby M, Reich PB, Wright IJ et al. 2011. TRY – a global database of plant biomass storage of an eastern Amazonian rainforest. New Phytologist 187: 579– traits. Global Change Biology 17: 2905–2935. 591. Kennedy D, Swenson S, Oleson KW, Lawrence DM, Fisher R, da Costa ACL, da Costa ACL, Rowland L, Oliveira RS, Oliveira AAR, Binks OJ, Salmon Y, Gentine P. 2019. Implementing plant hydraulics in the Community Land Vasconcelos SS, Junior JAS, Ferreira LV, Poyatos R et al. 2018. Stand Model, version 5. Journal of Advances in Modeling Earth Systems 11: 485–513. dynamics modulate water cycling and mortality risk in droughted tropical Klein T. 2014. The variability of stomatal sensitivity to leaf water potential across forest. Global Change Biology 24: 249–258. tree species indicates a continuum between isohydric and anisohydric Dewar RC, Franklin O, M€a kel€a A, McMurtrie RE, Valentine HT. 2009. behaviours. Functional Ecology 28: 1313–1320. Optimal function explains forest responses to global change. BioScience 59: Lens F, Picon-Cochard C, Delmas CEL, Signarbieux C, Buttler A, Cochard H, 127–139. Jansen S, Chauvin T, Doria LC, Del Arco M et al. 2016. Herbaceous Dewar R, Mauranen A, M€a kel€a A, H€oltt€a T, Medlyn B, Vesala T. 2018. New angiosperms are not more vulnerable to drought-induced embolism than insights into the covariation of stomatal, mesophyll and hydraulic conductances angiosperm trees. Plant Physiology 172: 661–667. from optimization models incorporating nonstomatal limitations to Leuning R. 1995. A critical appraisal of a combined stomatal - photosynthesis photosynthesis. New Phytologist 217: 571–585. model for C3 plants. Plant, Cell & Environment 18: 339–355. Dharssi I, Vidale PL, Verhoef A, Macpherson B, Jones C, Best M. (2009). New Lin Y, Medlyn B, Duursma R. 2015. Optimal stomatal behaviour around the soil physical properties implemented in the Unified Model at PS18. Met Office world. Nature Climate Change 5: 459. technical report 528. Exeter, UK: Met Office. Loveland TR, Reed BC, Ohlen DO, Brown JF, Zhu Z, Yang L, Merchant JW. Duffy PB, Brando P, Asner GP, Field CB. 2015. Projections of future 2000. Development of a global land cover characteristics database and IGBP meteorological drought and wet periods in the Amazon. Proceedings of the DISCover from 1 km AVHRR data. International Journal of Remote Sensing 21: National Academy of Sciences, USA 112: 13172–13177. 1303–1330. Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com New 1636 Research Phytologist M€akel€a A, Berninger F, Hari P. 1996. Optimal control of gas exchange during Nepstad DC, De Carvalho CR, Davidson EA, Jipp PH, Lefebvre PA, Negreiros drought: theoretical analysis. Annals of Botany 77: 461–468. GH, Da Silva ED, Stone TA, Trumbore SE, Vieira S. 1994. The role of deep Manzoni S, Vico G, Palmroth S, Porporato A, Katul G. 2013. Optimization of roots in the hydrological and carbon cycles of Amazonian forests and pastures. stomatal conductance for maximum carbon gain under dynamic soil moisture. Nature 372: 666–669. Advances in Water Resources 62: 90–105. Ocheltree TW, Nippert JB, Prasad PVV. 2016. A safety vs efficiency trade-off Marengo JA, Souza CM, Thonicke K, Burton C, Halladay K, Betts RA, Alves identified in the hydraulic pathway of grass leaves is decoupled from LM, Soares WR. 2018. Changes in climate and land use over the Amazon photosynthesis, stomatal conductance and precipitation. New Phytologist 210: region: current and future variability and trends. Frontiers in Earth Science 6: 97–107. 228. Powell TL, Galbraith DR, Christoffersen BO, Harper A, Imbuzeiro HMA, Martinez-Vilalta J, Poyatos R, Aguade D, Retana J, Mencuccini M. 2014. A new Rowland L, Almeida S, Brando PM, da Costa ACL, Costa MH et al. 2013. look at water transport regulation in plants. New Phytologist 204: 105–115. Confronting model predictions of carbon fluxes with measurements of Amazon Medlyn BE, Barton CVM, Broadmeadow MSJ, Ceulemans R, De Angelis P, forests subjected to experimental drought. New Phytologist 200: 350–365. Forstreuter M, Freeman M, Jackson SB, Kellom€aki S, Laitat E et al. 2001. Prentice IC, Dong N, Gleason SM, Maire V, Wright IJ. 2014. Balancing the Stomatal conductance of forest species after long-term exposure to elevated costs of carbon gain and water transport: testing a new theoretical framework CO2 concentration: a synthesis. New Phytologist 149: 247–264. for plant functional ecology. Ecology letters 17: 82–91. Medlyn BE, De Kauwe MG, Zaehle S, Walker AP, Duursma RA, Luus K, Qu X, Cao B, Kang J, Wang X, Han X, Jiang W, Shi X, Zhang L, Cui L, Hu Z Mishurov M, Pak B, Smith B, Wang YP et al. 2016. Using models to guide et al. 2019. Fine-tuning stomatal movement through small signaling peptides. field experiments: a priori predictions for the CO2 response of a nutrient- and Frontiers in Plant Science 10: 69. water-limited native Eucalypt woodland. Global Change Biology 22: 2834– Restrepo-Coupe N, Levine NM, Christoffersen BO, Albert LP, Wu J, Costa 2851. MH, Galbraith D, Imbuzeiro H, Martins G, da Araujo AC et al. 2017. Do Medlyn BE, Duursma RA, Eamus D, Ellsworth DS, Prentice IC, Barton CVM, dynamic global vegetation models capture the seasonality of carbon fluxes in Crous KY, De Angelis P, Freeman M, Wingate L. 2011. Reconciling the the Amazon basin? A data-model intercomparison. Global Change Biology 23: optimal and empirical approaches to modelling stomatal conductance. Global 191–208. Change Biology 17: 2134–2144. Rogers A, Medlyn BE, Dukes JS, Bonan G, von Caemmerer S, Dietze MC, Medlyn BE, Zaehle S, De Kauwe MG, Walker AP, Dietze MC, Hanson PJ, Kattge J, Leakey ADB, Mercado LM, Niinemets U et al. 2017. A roadmap for Hickler T, Jain AK, Luo Y, Parton W et al. 2015. Using ecosystem improving the representation of photosynthesis in Earth system models. New experiments to improve vegetation models. Nature Climate Change 5: 528. Phytologist 213: 22–42. Meinzer FC, Johnson DM, Lachenbruch B, McCulloh KA, Woodruff DR. Rowland L, Da Costa ACL, Galbraith DR, Oliveira RS, Binks OJ, Oliveira 2009. Xylem hydraulic safety margins in woody plants: coordination of AAR, Pullen AM, Doughty CE, Metcalfe DB, Vasconcelos SS et al. 2015. stomatal control of xylem tension with hydraulic capacitance. Functional Death from drought in tropical forests is triggered by hydraulics not carbon Ecology 23: 922–930. starvation. Nature 528: 119–122. Meir P, Brando PM, Nepstad D, Vasconcelos S, Costa ACL, Davidson E, Saleska SR, da Rocha HR, Huete AR, Nobre AD, Artaxo P, Shimabukuro YE. Almeida S, Fisher RA, Sotta ED, Zarin D et al. 2009. The effect of drought 2013. LBA-ECO CD-32 flux tower network data compilation, Brazilian Amazon: on Amazonian rain forests. In: Gash J, Keller M, Bustamante M, Silva Dias P, 1999–2006. Data set. Oak Ridge, TN, USA: Oak Ridge National Laboratory eds. Amazonia and global change. Washigton, DC, USA: American Geophysical Distributed Active Archive Center. Union, 429–449. Savage VM, Bentley LP, Enquist BJ, Sperry JS, Smith DD, Reich PB, von Meir P, Mencuccini M, Binks O, da Costa AL, Ferreira L, Rowland L. 2018. Allmen EI. 2010. Hydraulic trade-offs and space filling enable better Short-term effects of drought on tropical forest do not fully predict impacts of predictions of vascular structure and function in plants. Proceedings of the repeated or long-term drought: gas exchange versus growth. Philosophical National Academy of Sciences, USA 107: 22722–22727. Transactions of the Royal Society of London. Series B: Biological Sciences 373: Sellers PJ. 1985. Canopy reflectance, photosynthesis and transpiration. 20170311. International Journal of Remote Sensing 6: 1335–1372. Meir P, Wood TE, Galbraith DR, Brando PM, Da Costa ACL, Rowland L, Sheffield J, Wood EF. 2008. Projected changes in drought occurrence under Ferreira LV. 2015. Threshold responses to soil moisture deficit by trees and soil future global warming from multi-model, multi-scenario, IPCC AR4 in tropical rain forests: insights from field experiments. BioScience 65: 882–892. simulations. Climate Dynamics 31: 79–105. Melotto M, Underwood W, He SY. 2008. Role of stomata in plant innate Sitch S, Huntingford C, Gedney N, Levy PE, Lomas M, Piao SL, Betts R, Ciais immunity and foliar bacterial diseases. Annual Review of Phytopathology 46: P, Cox P, Friedlingstein P et al. 2008. Evaluation of the terrestrial carbon 101–122. cycle, future plant geography and climate-carbon cycle feedbacks using five Mencuccini M, Manzoni S, Christoffersen B. 2019a. Modelling water fluxes in Dynamic Global Vegetation Models (DGVMs). Global Change Biology 14: plants: from tissues to biosphere. New Phytologist 222: 1207–1222. 2015–2039. Mencuccini M, Rosas T, Rowland L, Choat B, Cornelissen JHC, Jansen S, Sperry JS. 2004. Coordinating stomatal and xylem functioning–an evolutionary Kramer K, Lepenas A, Manzoni S, Niinemets U et al. 2019b. Leaf economics perspective. New Phytologist 162: 568–570. and xylem hydraulics drive leaf/wood area ratios. New Phytologist 224: 1544– Sperry JS, Venturas MD, Anderegg WRL, Mencuccini M, Mackay DS, Wang Y, 1556. Love DM. 2017. Predicting stomatal responses to the environment from the Mercado LM, Bellouin N, Sitch S, Boucher O, Huntingford C, Wild M, Cox optimization of photosynthetic gain and hydraulic cost. Plant Cell and PM. 2009. Impact of changes in diffuse radiation on the global land carbon Environment 40: 816–830. sink. Nature 458: 1014. Stocker BD, Zscheischler J, Keenan TF, Prentice IC, Pe~ nuelas J, Seneviratne SI. Mercado LM, Huntingford C, Gash JHC, Cox PM, Jogireddy V. 2007. 2018. Quantifying soil moisture impacts on light use efficiency across biomes. Improving the representation of radiation interception and photosynthesis for New Phytologist 218: 1430–1449. climate model applications. Tellus, Series B: Chemical and Physical Meteorology Trugman AT, Medvigy D, Mankin JS, Anderegg WRL. 2018. Soil moisture 59: 553–565. stress as a major driver of carbon cycle uncertainty. Geophysical Research Letters Moore JR, Zhu K, Huntingford C, Cox PM. 2018. Equilibrium forest 45: 6495–6503. demography explains the distribution of tree sizes across North America. Ukkola AM, De Kauwe MG, Pitman AJ, Best MJ, Abramowitz G, Haverd V, Environmental Research Letters 13: 084019. Decker M, Haughton N. 2016. Land surface models systematically Mott KA. 1988. Do stomata respond to CO2 concentrations other than overestimate the intensity, duration and magnitude of seasonal-scale intercellular? Plant Physiology 86: 200–203. evaporative droughts. Environmental Research Letters 11: 104012. Nash JE, Sutcliffe JV. 1970. River flow forecasting through conceptual models Venturas MD, Sperry JS, Love DM, Frehner EH, Allred MG, Wang Y, part I – a discussion of principles. Journal of Hydrology 10: 282–290. Anderegg WRL. 2018. A stomatal control model based on optimization of New Phytologist (2020) 226: 1622–1637 Ó 2020 The Authors www.newphytologist.com New Phytologist Ó 2020 New Phytologist Trust New Phytologist Research 1637 carbon gain versus hydraulic risk predicts aspen sapling responses to drought. Fig. S4 Daily drought evolution modelled by JULES and New Phytologist 220: 835–850. JULES-SOX. Walsh RPD, Lawler DM. 1981. Rainfall seasonality: description, spatial patterns and change through time. Weather 36: 201–208. Wang Y, Sperry JS, Venturas MD, Trugman AT, Love DM, Anderegg WR. Fig. S5 Differences between evapotranspiration predicted by 2019. The stomatal response to rising CO2 concentration and drought is JULES and JULES-SOX. predicted by a hydraulic trait-based optimization model. Tree Physiology 39: 1416–1427. Fig. S6 Seasonal variation in modelled and observed evapotran- Williams M, Rastetter EB, Fernandes DN, Goulden ML, Wofsy SC, Shaver GR, Melillo JM, Munger JW, Fan SM, Nadelhoffer KJ. 1996. Modelling the spiration. soil-plant-atmosphere continuum in a Quercus-acer stand at Harvard forest: the regulation of stomatal conductance by light, nitrogen and soil/plant Fig. S7 Seasonal variation in modelled and observed soil mois- hydraulic properties. Plant, Cell & Environment 19: 911–927. ture. Winkler AJ, Myneni RB, Alexandrov GA, Brovkin V. 2019. Earth system models underestimate carbon fixation by plants in the high latitudes. Nature Communications 10: 885. Notes S1 Analytical SOX derivation. Wolf A, Anderegg WRL, Pacala SW. 2016. Optimal stomatal behavior with competition for water and risk of hydraulic impairment. Proceedings of the Notes S2 Computing A and V numerical derivatives. National Academy of Sciences, USA 113: E7222–E7230. Wolfe BT, Sperry JS, Kursar TA. 2016. Does leaf shedding protect stems from Notes S3 Leaf photosynthesis model solved for stomatal conduc- cavitation during seasonal droughts? A test of the hydraulic fuse hypothesis. New Phytologist 212: 1007–1018. tance. Worrall JJ, Marchetti SB, Egeland L, Mask RA, Eager T, Howell B. 2010. Effects and etiology of sudden aspen decline in southwestern Colorado, USA. Notes S4 Whole-tree hydraulic conductance and xylem tapering Forest Ecology and Management 260: 638–648. calculations. Xu X, Medvigy D, Powers JS, Becknell JM, Guan K. 2016. Diversity in plant hydraulic traits explains seasonal and inter-annual variations of vegetation dynamics in seasonally dry tropical forests. New Phytologist 212: 80–95. Table S1 Details for the data used in the SOX leaf-level evalua- tion. Supporting Information Table S2 Sites used for the ecosystem-level evaluation of JULES- Additional Supporting Information may be found online in the SOX. Supporting Information section at the end of the article. Fig. S1 Responses of Eqns 4 and 5 to environmental drivers. Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the Fig. S2 Maps of the observation sites used on this study. authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office. Fig. S3 Agreement between numerical and analytical SOX. See also the Commentary on this article by Anderegg & Venturas, 226: 1535–1538. Ó 2020 The Authors New Phytologist (2020) 226: 1622–1637 New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com