Modelling tropical forest responses to drought and El Niño with a stomatal optimization model based on xylem hydraulics

The current generation of dynamic global vegetation models (DGVMs) lacks a mechanistic representation of vegetation responses to soil drought, impairing their ability to accurately predict Earth system responses to future climate scenarios and climatic anomalies, such as El Niño events. We propose a simple numerical approach to model plant responses to drought coupling stomatal optimality theory and plant hydraulics that can be used in dynamic global vegetation models (DGVMs). The model is validated against stand-scale forest transpiration (E) observations from a long-term soil drought experiment and used to predict the response of three Amazonian forest sites to climatic anomalies during the twentieth century. We show that our stomatal optimization model produces realistic stomatal responses to environmental conditions and can accurately simulate how tropical forest E responds to seasonal, and even long-term soil drought. Our model predicts a stronger cumulative effect of climatic anomalies in Amazon forest sites exposed to soil drought during El Niño years than can be captured by alternative empirical drought representation schemes. The contrasting responses between our model and empirical drought factors highlight the utility of hydraulically-based stomatal optimization models to represent vegetation responses to drought and climatic anomalies in DGVMs. This article is part of a discussion meeting issue ‘The impact of the 2015/2016 El Niño on the terrestrial tropical carbon cycle: patterns, mechanisms and implications’.

CBE, 0000-0002-7795-2574; PM, 0000-0002-2362-0398 The current generation of dynamic global vegetation models (DGVMs) lacks a mechanistic representation of vegetation responses to soil drought, impairing their ability to accurately predict Earth system responses to future climate scenarios and climatic anomalies, such as El Niñ o events. We propose a simple numerical approach to model plant responses to drought coupling stomatal optimality theory and plant hydraulics that can be used in dynamic global vegetation models (DGVMs). The model is validated against stand-scale forest transpiration (E) observations from a long-term soil drought experiment and used to predict the response of three Amazonian forest sites to climatic anomalies during the twentieth century. We show that our stomatal optimization model produces realistic stomatal responses to environmental conditions and can accurately simulate how tropical forest E responds to seasonal, and even long-term soil drought. Our model predicts a stronger cumulative effect of climatic anomalies in Amazon forest sites exposed to soil drought during El Niñ o years than can be captured by alternative empirical drought representation schemes. The contrasting responses between our model and empirical drought factors highlight the utility of hydraulically-based stomatal optimization models to represent vegetation responses to drought and climatic anomalies in DGVMs.
This article is part of a discussion meeting issue 'The impact of the 2015/ 2016 El Niñ o on the terrestrial tropical carbon cycle: patterns, mechanisms and implications'.
Over the last decade, important advances have been made to improve our understanding of the physiological processes determining plant responses to drought [7,8]. Experimental manipulation and field observations have shown that xylem hydraulic conductance loss is an important mechanism triggering drought-induced plant mortality [9][10][11][12]. One of the mechanisms that plants employ to avoid reaching potentially lethal embolism thresholds is the regulation of canopy water potential (C c ) through stomatal control, which creates a coordination between stomatal responses and plant hydraulic conductance losses [13][14][15][16]. While process-based models of stomatal functioning based on plant hydraulics have been proposed recently [17][18][19], most dynamic global vegetation models (DGVMs) rely on empirical drought factors to represent stomatal responses to soil drought [20][21][22][23]. These empirical approaches can perform well under many conditions [24][25][26], but they lack the generality of models that use physiological and ecological theory to predict the responses of vegetation and the global carbon cycle to drier climates [21,27], such as the Amazon climate during El Niñ o events. In this study we describe and test a new model of stomatal response to drought that is numerically simple enough to implement in a DGVM applicable at large spatial scales, without losing recent theoretical advancements made in the field of plant hydraulics and stomatal optimization theory [17,18].
Our model is based on optimality theory, that is, plant structure and functioning have evolved to maximize efficiencies within the limits of genotypic variation and physico-chemical constraints [28][29][30][31][32]. This principle has been widely used to predict stomatal responses to environmental conditions, starting with Cowan [33] and Cowan & Farquhar [34], where stomata are assumed to maximize carbon assimilation (A as carbon mass) while minimizing transpiration (E as water mass) over a given time interval (dt). This concept can be represented by maximizing the function A2lE over dt. The parameter l represents the marginal carbon cost of water (carbon mass per water mass). This E-based optimization approach provides an alternative to empirical models that has been widely used [35][36][37][38][39], such as in Medlyn et al. [40] to derive the unified stomatal optimization model (USO). The USO shows the potential of the E-based optimization theory to predict stomatal conductance (g c ) responses to environmental drivers [40,41]. However, E-based optimization does not account for soil drought effects on g c , which need to be represented empirically as in Zhou et al. [25,26], or with semi-empirical drought factors [36,37].
We represent drought effects on stomatal conductance coupling plant hydraulics with stomatal optimality theory, following the principles outlined in Wolf et al. [19] and Sperry et al. [18] and using an optimization routine similar to Friend [42]. Sperry et al. [18] propose that the costs associated with stomatal opening can be represented as the loss of the plant capacity to transport water, which allows us to replace the need for l with hydraulic traits that determine plant vulnerability to drought-induced embolism. Plant hydraulic traits that determine xylem vulnerability to embolism at the branch-level are currently available for a large number of species of different biomes [43], which makes the hydraulics-based optimization approach particularly attractive for inclusion in ecosystem models.
In this study we validate a stomatal optimization model based on xylem hydraulics (SOX) against scaled-up sap flux observations from an Amazon forest site subject to long-term experimental drought [11,44] and evaluate its predictions against other stomatal models. Subsequently, we investigate how our model predictions differs from empirical drought factors at simulating the response of Amazon forest sites to climatic anomalies during the twentieth century.

Material and methods (a) Model description
The SOX model assumes the loss of xylem hydraulic conductance is the main cost associated with stomatal opening. Therefore, we calculate the optimal stomatal conductance for a given set of environmental conditions as the value that maximizes A (mol m 22 s 21 ) given concurrent hydraulic conductance losses, using a numerical routine similar to the PGEN model [42]. A schematic representation of the model is shown in figure 1. The numerical routine we describe here can be coupled to any photosynthesis model that computes A from environmental inputs and the leaf intercellular CO 2 concentration (c i , mol mol 21 ). In this study we use the photosynthesis model from Collatz et al. [45], following Clark et al. [20], described in electronic supplementary material, appendix S1. From an initial value for A, we derive the canopy conductance to CO 2 (g c , mol m 22 leaf s 21 ) and transpiration (E, mol m 22 leaf s 21 ) as: where c a is the CO 2 concentration (mol mol 21 ) in the atmosphere (assumed to be equal to the leaf surface). The leaf-to-air vapour pressure deficit (D, mol mol 21 ) is calculated with the assumption that canopy temperature is close to air temperature. These assumptions are justified on the basis that the model implemented in this study is the proof of concept of a scheme designed to be coupled to larger scale models that often employ more detailed calculations of canopy aerodynamical resistance and energy balance (e.g. Best et al. [46]). The constant 1.6 is the ratio of water vapour to CO 2 diffusivities in the air. The resulting value of E is used to calculate the xylem water potential at the canopy (C c , MPa) using Darcy's Law, assuming steady state conditions (i.e. no contribution of stored water to transpiration): where k rc is the root -canopy hydraulic conductance (mol m 22 leaf s 21 MPa 21 ) and C c,pd is C c at the pre-dawn which, assuming no night-time transpiration, can be approximated as the root C (C r ) adjusted for the canopy height (h, m) induced C gradient: where r is the water density (997 kg m 23 ), g is the Earth's gravitational acceleration (9.8 m s 22 ) and the 10 26 converts Pa to MPa. Stored water can contribute significantly to tropical vegetation transpiration [47,48]. However, this contribution is lower during periods of high water stress, when the internal water reserves are depleted, which makes equation (2.3) a reasonable approximation when C c is more relevant for our model. The k rc in equation (2.3) itself depends on C c for its computation as k rc declines from its maximum value (k rc,max ) as the xylem pressure (C) drops due to cavitation-induced embolism rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 373: 20170315 formation [49]. This process can be described with a function such as the inverse polynomial from Manzoni et al. [50]: where C 50 is C when k rc ¼ 0.5k rc,max and a controls the shape of the function. The empirical relationship between C 50 and a from Christoffersen et al. [51] described in electronic supplementary material, appendix S2, reduces the plant hydraulic parameters needed in SOX to only C 50 , k rc,max and h used in equation (2.4). SOX is designed as a dynamic model that uses the k rc produced at the previous timestep (k rc[t21] ) to compute the current timestep C c and k rc via equation (2.3). The purpose of this choice is to facilitate the incorporation of long-term drought effects in k rc , associated with the incomplete recovery of cavitation [52]. In this study we assume k rc recovers instantaneously as C c and C c,pd increase following a rain event. This might overestimate the fluxes immediately after the dry season, depending on the forest recovery rates from embolism through growth [53,54] or other processes [55,56]. More complex schemes describing partial and gradual k rc recovery processes will be explored in future studies. Sperry & Love [17] and Sperry et al. [18] employ the Kirchhoff transform in equation (2.5) to account for the gradual C drop along the tree, computing k rc as: In SOX we represent the gradual C drop along the tree using the middle value of the root-canopy gradient (C c,mid ): ð2:7Þ Using f (C c,mid ) is numerically simpler than equation (2.6) and provides similar results within a realistic range of C c,pd and C c (electronic supplementary material, figure S1). The k rc produced by f (C c,mid ) normalized as a function of k rc,max , giving k cost , represents the costs of stomatal aperture in SOX: ð2:8Þ photosynthesis model (appendix S1): xylem hydraulic conductance: material, appendix S1) is used to produce the gross carbon assimilation (A, green lines in subpanels a and b) for a set of environmental conditions and initial leaf internal CO 2 concentration (c i ) value. We use equations (2.1)-(2.3) from the main text to calculate the xylem water potential at the canopy (C c ) assoc i ated with the given A value. The midpoint of the root to canopy C c gradient is used to compute the normalized root to canopy hydraulic conductance (k cost ), which represents the cost of stomatal aperture in SOX (red dashed lines in a and b). The k cost and A are used to numerically find the optimum c i (circles on the x-axis of a and b), at the maximum point of the A . k cost function (black lines in a and b). In the subpanels a and b, we represent the SOX optimization routine at increasing soil drought stress (lighter coloured lines represent lower soil water potential, C s ) at two different levels of atmospheric demand. The optimum c i value of the interval evaluated by the SOX routine (0 to atmospheric CO 2 levels) is used to calculate the optimum A, denoted by the circles on the A-c i curve. The equations in the dashed box represent how changes in soil conductance can be incorporated in SOX by modelling the soil to root hydraulic conductance (k rc ) as a function of C s as described in electronic supplementary material, appendix S4. (Online version in colour.) rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 373: 20170315 Lower k cost implies a higher cost associated with a given level of stomatal aperture. These costs are balanced with A using a numerical optimization routine. A detailed discussion of the differences between the cost functions used in SOX and Sperry et al. [18] is given in electronic supplementary material, appendix S3.
The SOX optimization routine is implemented in this paper following similar principles to the PGEN model optimization routine [44], which assumes the optimum g c can be found where the product between A and its unitless drought factor are maximized. In SOX, as A and k cost are functions of c i , the optimum c i , hereafter c i,opt , for a given set of environmental conditions is found at: We use an algorithm (see the SOX model code available as electronic supplementary material) to evaluate c i over the interval (0, c a ) and find the solution to equation (2.9). The c i,opt is used to calculate optimum values of A, g c , E and C c using the photosynthesis model in electronic supplementary material, appendix S1 and equations (2.1) -(2.3).
Changes in soil hydraulic conductance can also be included in SOX by computing C r as a function of soil-to-root conductance as shown in figure 1 and explained in detail in electronic supplementary material, appendix S4. The model evaluations conducted in this study used the simplest version of SOX without the equations from electronic supplementary material, appendix S4 (i.e. assuming C r %C s ), unless noted otherwise.

(b) Model evaluation
The model was written in R (v. 3.4.2; [57]), and the code is available as electronic supplementary material; all the subsequent analyses were also conducted in R. The model responses to environmental drivers were evaluated by holding all meteorological inputs constant at their default values (table 1) and varying a single input at a time. Because equation (2.3) depends on k rc[t21] , we run the model at constant environmental conditions for 50 iterations to evaluate SOX instantaneous responses to the environment. This procedure is not necessary when SOX is run as a dynamic model, which is the case when SOX is coupled to a DGVM or in the subsequent model evaluations we conduct in this study.
We evaluated the model capacity to produce realistic predictions of vegetation response to seasonal and experimental soil drought using observations from an evergreen broadleaf tropical forest located in Caxiuanã National Forest in the eastern Brazilian Amazon (electronic supplementary material, figure S3 for site details). We compared the modelled E with the stand-scale sap flux data from two 1 ha plots at the site. One of the plots has been subjected to a throughfall exclusion treatment (TFE) since 2001 [11,58], which provides an ideal scenario to test the capacity of SOX to reproduce vegetation response to severe soil drought. Details on sap flux data collection and procedures to scale the data from tree to stand-level can be found in da Costa et al. [59]. The meteorological forcing data were collected at the top of a 40 m tower at the site, and the soil moisture data were measured with time-domain reflectometry sensors placed at 0.0-0.3, 0.5, 1 and 2.5 m depth. We used the Clapp & Hornberger [60] equation from electronic supplementary material, appendix S4 to obtain C s from the site observations of root mass-weighted soil moisture content (u, m 3 m 23 ), with the soil hydraulic parameters derived from the soil ancillary data used in the Hadley Centre Global Environmental Model Earth System Model (HadGEM2-ES) [44], which is based on the Harmonized World Soil Database (v. 1.2) [61]. The root biomass profile was modelled using the equations from Best et al. [46] assuming soil and root depth were 3 m, which is the default value for broadleaf evergreen tropical trees (BET) used in JULES [20]. We used the site-averaged values of tree hydraulic and physiological data, or the reference JULES values for BET (electronic supplementary material, table S1). Vegetation k rc,max was obtained from branch-level xylem specific conductivity (K x ), h, the ratio between sapwood area and leaf area (i.e. the Huber value, h v ) and a tapering correction factor calculated following Christoffersen et al. [51]; see full description of these calculations in electronic supplementary material, appendix S5. We scale the model predictions from leaf to plot area using the big leaf approach as described in Clark et al. [20], with the light extinction coefficient set to the default BET value (0.5) and leaf area index (LAI) fixed at the mean value observed at the site (4.8 m 2 leaf m 22 soil). We consider the use of a fixed LAI in this study is the most parsimonious choice for the purpose of validating our model, considering the small LAI changes observed at the site (standard deviation of 0.5 m 2 leaf m 22 soil). We compare SOX agreement with observations against a model that uses a drought representation model based on the b-function (b fun ) soil drought factor described by Cox et al. [24]. A description of this model is given in electronic supplementary material figure  S4. We fitted the relationship between A and stomatal conductance of water (g w, mol m 22 s 21 ) predicted by SOX to the unified stomatal optimization model (USO) of Medlyn et al. [40], described in electronic supplementary material, appendix S6. index, which is calculated as the mean sea surface temperature (SST) anomaly from 58N to 58S and 150 -908W [62]. Additionally, we used the site-specific monthly soil moisture product from JULES, applied following the TRENDY protocol [27,63], to drive our simulations. The soil hydraulic parameters for each site were obtained from the HadGEM2-ES soil ancillaries [64]. We used the plant inputs given in electronic supplementary material, table S1 to represent the Caxiuanã site. For the Tapajó s and Manaus sites, we used the mean plant hydraulic [65] and photosynthetic parameters measured at each site to parameterize the models (electronic supplementary material, table S2), while the other parameters were assumed equal to those of the Caxiuanã site. The vegetation hydraulic trait sampling in each site represents approximately 40, 36 and 15% of the forest basal area for the Caxiuanã, Tapajó s and Manaus sites, respectively.
We measured the effects of climatic anomalies on air temperature and atmospheric demand (T a and D) and soil water availability (C s ) by conducting experiments where we drove the models with the 6-hourly data that correspond to an average year based on the historical climate from the CRU-NCEP dataset (1901 to 2016, electronic supplementary material, figure S5). This procedure eliminates climatic anomalies, such as those associated with El Niñ o (electronic supplementary material, figure S3). In total we conducted four simulations for each site: Sim1 is the control run using the unaltered CRU-NCEP dataset; Sim2 is the run without anomalies in T a and D; Sim3 is the run without anomalies in C s ; and Sim4 is without anomalies in any of the previously mentioned variables (T a , D and C s , see electronic supplementary material, table S3 for summary).

Results
(a) Theoretical responses to environment SOX predicts that resistance to cavitation produces a stomata behaviour more responsive to changes in incident photosynthetically active radiation (I PAR ), c a , T c and D. However, C s has a much stronger effect in plants more vulnerable to cavitation (figure 2; electronic supplementary material, figure S6).
The asymptotic stomatal response to I PAR (figure 2a) is caused by the light-limitation predicted by the photosynthesis model (electronic supplementary material, appendix S1). The SOX predictions represent a hydraulic effect on the plant light response, as plants more resistant to cavitation can sustain light-saturated g c 2-3 times higher than the more vulnerable plants. The SOX response to c a (figure 2c,d) is driven by equation (2.2) producing lower g c for a given A as c a increases. The lower intrinsic water use efficiency (i.e. A/g c ) at low c a is partially compensated in cavitationresistant plants, which can maintain stomatal aperture with low cavitation costs, reducing the c a and c i gradient (electronic supplementary material, figure S6). The T c response in figure 2e,f results from the V cmax -T c relationship present in the photosynthesis model (electronic supplementary material, appendix S1), which is more pronounced in plants vulnerable to cavitation. These plants maintain a greater distance to their potential maximum g c due to premature, hydraulicallyinduced stomatal closure.
The SOX response to atmospheric demand (figure 2g,h) results from the increased k rc loss associated with the lower C c necessary to maintain carbon assimilation when the atmosphere is drier (i.e. higher D, see equations (2.2) and (2.3)). The lower C c induces an exponential decline in g c as the c i,opt shifts closer to 0 (figure 1). This pattern is commonly observed [66][67][68], and more accentuated in plants with a less negative C 50 owing to their increased cavitation costs as the atmosphere dries ( figure 2g,h). The b fun predicts a more gradual stomatal closure in response to D, which approximates SOX predictions for C 50 ¼ 22.5 MPa when D . 0.02 mol mol 21 (figure 2g).
The SOX response to a drying soil emerges from the same mechanism as its responses to D, that is, increased cavitation costs due to the lower C c necessary to maintain A when C s is low. SOX predicts a more gradual decline in g c in response to drying soil when compared with the b fun model (figure 2i,j ). The b model predicts g c ¼ 0 when C s ¼ 21.5 MPa and u ¼ u w (electronic supplementary material, figure S4; figure 2i,j), whereas SOX predicts that even plants very vulnerable to cavitation (C 50 ¼ 21 MPa) still have 30% of their maximum g c at the same C s. Accounting for changes in k sr with equations in the electronic supplementary material, appendix S3 makes g c more responsive to C s , affecting particularly plants more resistant to cavitation. At C s ¼ 25 MPa, even plants with C 50 ¼ 25 MPa will have dropped to 1% of their maximum g c , whereas in the model that assumes C r % C s the plant would still have 58% of its maximum g c . Using a steeper vulnerability curve (higher a from equation (2.5)) also greatly increases the plant sensitivity to soil drought ( figure 2k,l ). The a also affects the C c response to C s , which is linear when a is low, but higher a produces a more stable C c at high soil moisture (figure 2j,l ).

(b) Model evaluation
The SOX predictions agree with the E observed at both plots in Caxiuanã more consistently than the alternative b fun models (figure 3; electronic supplementary material, figure S7). We were able to approximate the sap flux at both the control and the TFE plot, even though we made the simplifying assumption that the vegetation of both plots was identical (electronic supplementary material, table S1) and used the relationship from Christoffersen et al. The relationship between g w and A predicted by SOX agrees with the Medlyn et al. [40] USO model under high I PAR (figure 4), and produces estimates of g 1 and its response to C c,pd (electronic supplementary material, table S4; figure 4) within the range observed for tropical trees in other studies [26,69]. Deviations from the 1 : 1 line occur at low g w and are associated with low I PAR periods. These deviations are present even if we set the minimum conductance parameter from USO, g 0 , to 0. Therefore, the SOX A-g w relationship implies a dependency of the water marginal carbon costs (related with the USO parameter g 1 , see electronic supplementary material, appendix S6) on the light regime that is not present in USO. The g 1 predicted by USO is lower at the TFE plot than in the control plot (electronic supplementary material, table S4), indicating that SOX predicts a higher water carbon cost at the TFE. This pattern cannot be observed  Figure 2. Response curves of stomatal conductance to CO 2 (g c , left panels) and canopy water potential (C c , right panels) to changes in incident photosynthetically active radiation (I PAR ), atmospheric carbon dioxide partial pressure (c a ), canopy temperature (T c ), vapour pressure deficit (D) and soil water potential (C s ). All the other environmental inputs and plant inputs were held constant at their default values (table 1)

Discussion
Our results show that a xylem hydraulics-based stomatal optimization scheme can produce realistic stomatal responses  This finding complements recent studies that have established the theoretical basis for a hydraulically-based model of plant stomatal responses to drought [17,18], and supports the recent findings of Anderegg et al. [70], showing the potential of xylem hydraulics-based optimization approaches to simulate the responses of tropical forests to drought. The SOX predictions agree with other models based on the optimality theory, such as the USO, under most circumstances ( figure 4). However, SOX predictions are considerably different from the drought factor approach, represented here by the b fun model (figures 4 and 5; electronic supplementary material, figure S4). The drastic differences that emerge from long-term simulations between SOX and the b fun model ( figure 5) highlight the importance of using a more mechanistic plant hydraulic representation to simulate the effects of climatic anomalies, such as El Niñ o, on forest carbon and water fluxes.
The drastic divergence between the b fun predictions and observations found in our study (figure 3) could be partly explained by the choice of using soil moisture data to drive our simulations. The b fun model and other empirical drought factors used in DGVMs are often coupled to a soil hydrology scheme [20,46,71,72]. The influence of plant transpiration on soil moisture dynamics could attenuate the extreme soil drought responses we observed ( figure 5). However, other studies show that even when soil hydrology is accounted for, b fun might still overestimate soil drought responses [73]. The approach we adopted can be considered a conservative test of the model capability to predict forest transpiration, as no assumptions were made modelling the soil water dynamics.

(a) Generality and limitations of SOX
Our model is designed to be coupled to large-scale ecosystem models such as DGVMs, and therefore its performance depends on the coupled routines representing vegetation  1900 1940 1980 2020 1900 1940 1980 2020 1900 1940 1980 2020 1900 1940 1980 2020 1900 1940 1980 2020 1900 1940 1980 2020 1900 1940 1980  processes (e.g. photosynthesis, canopy energy balance, and phenology), soil hydrology and atmospheric processes. For this study we assumed constant leaf area over time when scaling from leaf-level to plot-level in figure 3, as the LAI variation at this site is relatively small. However, phenology schemes [20,74] should be easily integrated with SOX; in addition, our model opens the possibility for plant hydraulics-driven phenology schemes. Linking vegetation phenology to drought responses is a much-needed functionality in many ecosystem models [21,75], and could further improve how SOX represents vegetation responses to extreme drought (figure 3b). The hydraulic processes represented by SOX also open up the possibility for a more explicit representation of drought-induced mortality in DGVMs. The thresholds of hydraulic conductance loss associated with increased risks of plant mortality, thought to be close to 0.5k rc,max for gymnosperms [12] and 0.12k rc,max for angiosperms [9,10], can be linked from the SOX output into a DGVM module that controls vegetation demographic processes, such as the TRIFFID module currently used in JULES [74]. The good performance of the simplified SOX implementation we show in this study, which is comparable even to that of more detailed models previously used on the site [51,76], illustrates the parsimony of the xylem hydraulicsbased optimization approach. Our model evaluation at the Caxiuanã TFE plot shows that accounting for soil hydraulic conductance loss is an important step for reproducing long-term drought effects (electronic supplementary material, figure S7). These results complement previous work made at the site [76], showing that even after over a decade of experimental drought, soil hydraulic conductance loss remains an important driver for forest response to drought. Even accounting for changes in soil conductance, the performance of SOX in figure 3b shows that there is room for improvement in how we model long-term drought in SOX. Together with phenological responses to soil drought mentioned above, legacy effects of cavitation [52] could be an important mechanism driving the TFE plot responses. The SOX treatment of k rc in equation (2.3) makes it simple to incorporate the processes determining the recovery of k rc by the plant.
The accuracy of our model predictions requires further testing against observations from other ecosystems and plant functional types (PFTs). The agreement of our model predictions with data depends on the two main theoretical assumptions of optimality theory being satisfied: (1) that it is physiologically possible for plant stomata to operate close to the SOX definition of optimum, and (2) the optimization criterion used in SOX can be strongly linked to plant fitness [29 -32]. Plant stomata have been often observed to function close to a theoretical optimum [34,[38][39][40][41]77], but deviations from this behaviour have also been observed [78,79]. These departures can be interpreted as consequences of physical and biochemical limitations on stomatal reaction times [80]. These effects should be more conspicuous at short timescales and in PFTs with slower stomatal responses, such as gymnosperms [79]. Other mechanisms that have been proposed to cause stomatal departure from a theoretical optimum include non-stomatal limitations to A, such as a reduction of Rubisco activity [26,35] and mesophyll conductance [81].
The second SOX assumption concerns our optimization criterion as the maximization of the cost-regulated carbon assimilation product (A . k cost ). The optimality theory replaces the need for detailed physiological parameterization, with evolutionary assumptions that depend on the impact of specific processes and structures on the fitness of organisms [29 -32]. The link between A and plant fitness is clear, as the reproductive success of a plant depends on its energetic investment in reproductive tissues over its lifespan [82], and in tissues necessary for survival and acquisition of resources other than carbon. The cost term in SOX, represented by xylem hydraulics dysfunction, implies that the complete loss of hydraulic conductance (i.e. k cost ¼ 0) would be associated with plant mortality, which represents the ultimate fitness cost [18,19,82]. There is substantial evidence that high levels of xylem cavitation-induced embolism are in fact associated with plant mortality [9][10][11], which corroborates this assumption. Even non-lethal loss of hydraulic conductivity should be detrimental to plant fitness, as recovery of hydraulic conductivity through construction of new vessels [53,54,83], or through active refilling of embolized vessels [84][85][86][87], requires carbon investment, which would necessarily detract from plant tissue growth and reproductive investments. Differences in plant capabilities of recovering hydraulic conductivity, be it through refilling or through the construction of new vessels, imply that a given level of hydraulic damage predicted by the xylem vulnerability function might not fully represent the costs of stomatal opening, as the long-term carbon balance impact of embolism are not explicitly represented. Even though the normalized xylem vulnerability-based cost function we use here represents a satisfactory first approximation, an appropriated weighting of the carbon costs associated with the recovery of hydraulic conductivity [54,56] might be a necessary theoretical development to improve the generality and accuracy of xylem hydraulics-based optimization models.

(b) Agreement with alternative drought-representation schemes
The relationship between g w and A predicted by SOX agrees well with that of the USO model from Medlyn et al. [40], which reflects the agreement between the different optimization principles underlying each model. The USO assumes stomata maximize the mass of carbon gain per mass of water lost (i.e. A2El), while SOX maximizes the fraction of xylem lost per mass of carbon assimilated (A . k cost ). The association between these principles can be interpreted as a result of the dependency between E and k loss (equations (2.3) and (2.5)). As high E has no direct detrimental effect on plant fitness, its association with plant hydraulics provides the necessary theoretical link between E and plant fitness to satisfy the fundamental assumption of optimality theory [28 -32]. Xylem hydraulics-based optimization models have the advantage of combining stomatal responses to D and C s using a few hydraulic parameters (table 1) that are currently widely available [43]. The difference between our integrated drought representation and approaches usually employed in DGVMs that rely on combining two empirical/semiempirical functions [20,24,72] is highlighted in the longterm simulations and their responses to climatic anomalies (figure 5). The carbon assimilation in Caxiuanã and, especially, in Manaus was dominated by atmospheric anomalies, as there was little soil drought in the driving rstb.royalsocietypublishing.org Phil. Trans. R. Soc. B 373: 20170315 data used for this experiment. The soil moisture data used to drive these simulations were the product of large-scale JULES simulations and meteorological datasets (0.58 Â 0.58 resolution), which explains their contrast with the environmental data collected at the site that was used to drive the model in the evaluation against sap flux data from Caxiuanã (figure 3). Atmospheric demand is an important driver of vegetation carbon and water fluxes [88], and a more likely mechanism driving Amazon forest responses to climatic anomalies than soil water stress, as the latter often requires multiple years of sustained rainfall reduction to produce a significant response in tropical forests [11,44,59].
Tapajó s was the only site with a significant interannual C s variability (electronic supplementary material, figure S3), and it was the site where the divergences between the b fun model and SOX were largest (figure 5; electronic supplementary material, figures S9 and S10). The b fun excessive soil moisture response and highly variable response to climatic anomalies reflect the steep gradient between the critical and wilting points of the b fun equation (electronic supplementary material, figure S4), producing a stronger decline in g c in response to soil drought than SOX, especially for plants more resistant to cavitation and with lower a value (figure 2i,k). Other studies have also shown that the excessive stomatal regulation produced by the b fun produces divergences between model predictions and seasonal GPP patterns in Tapajó s [73]. The large discrepancy between the two models, especially over the last 50 years, indicates that tropical forest sites exposed to soil water limitations during El Niñ o years might have stronger responses to climatic anomalies than can be captured by models based on empirical drought factor schemes.

Conclusion
Our stomatal optimization model, SOX, provides a simple but theoretically robust approach to simulate tropical forest responses to drought, capable of reproducing the effects of even severe experimental droughts. A process-based representation of atmospheric and soil drought responses is essential for the unbiased simulation of tropical forest responses to El Niñ o-style climatic anomalies. Improving the representation of plant hydrodynamics is a priority for the current generation of ecosystem models [21 -23,27]. The flexibility, relative simplicity and small number of parameters required by SOX make it an attractive candidate to be used in large-scale modelling of tropical forest responses to climate change and extreme climatic anomalies. More studies are necessary to assess the generality of our approach in distinct PFTs and environments, and there is a potential need to incorporate additional mechanisms, such as processes involved in the recovery of hydraulic conductance, hydraulically-driven phenological changes, and mortality.
Data accessibility. The JULES soil moisture output used in this study as well as the meteorological driving data were obtained from the TRENDY project. The full TRENDY dataset (http://dgvm.ceh.ac.uk/index.html) is available, subject to the individual modelling group approval, via a request to S.S. (s.a.sitch@exeter.ac.uk). The sap flux data used for model validation is published in da Costa et al. [59]. The R code for the models used in this paper and the plant input data for each site used in this study are available as electronic supplementary material.