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

1. Introduction

El Niño events contribute to major climatic and ecologic impacts over the Amazon basin [1–4]. Climatically, El Niño events are known to make the climate of most of Amazonia drier and warmer, especially affecting the rainfall patterns in northern Amazonia [4]. This drier climate drives a shift in Amazon forest carbon balance towards a net carbon source to the atmosphere [1,5]. The mechanisms involved in this shift are thought to be related to temperature-induced increases in respiration (particularly soil respiration) and drought-induced decreases in gross primary productivity [3,5,6].

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–12]. One of the mechanisms that plants employ to avoid reaching potentially lethal embolism thresholds is the regulation of canopy water potential (Ψ_c) through stomatal control, which creates a coordination between stomatal responses and plant hydraulic conductance losses [13–16]. While process-based models of stomatal functioning based on plant hydraulics have been proposed recently [17–19], most dynamic global vegetation models (DGVMs) rely on empirical drought factors to represent stomatal responses to soil drought [20–23]. These empirical approaches can perform well under many conditions [24–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–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 A−λE over dt. The parameter λ 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–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 λ 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.

2. 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⁻² s⁻¹) 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₂ concentration (c_i, mol mol⁻¹). 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₂ (g_c, mol m⁻² leaf s⁻¹) and transpiration (E, mol m⁻² leaf s⁻¹) as:

and where c_a is the CO₂ concentration (mol mol⁻¹) in the atmosphere (assumed to be equal to the leaf surface). The leaf-to-air vapour pressure deficit (D, mol mol⁻¹) 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₂ diffusivities in the air.

• Download figure
• Open in new tab
• Download PowerPoint

The resulting value of E is used to calculate the xylem water potential at the canopy (Ψ_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⁻² leaf s⁻¹ MPa⁻¹) and Ψ_c,pd is Ψ_c at the pre-dawn which, assuming no night-time transpiration, can be approximated as the root Ψ (Ψ_r) adjusted for the canopy height (h, m) induced Ψ gradient:where ρ is the water density (997 kg m⁻³), g is the Earth's gravitational acceleration (9.8 m s⁻²) and the 10⁻⁶ 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 is more relevant for our model. The k_rc in equation (2.3) itself depends on Ψ_c for its computation as k_rc declines from its maximum value (k_rc,max) as the xylem pressure (Ψ) drops due to cavitation-induced embolism formation [49]. This process can be described with a function such as the inverse polynomial from Manzoni et al. [50]:where Ψ₅₀ is Ψ when k_rc = 0.5k_rc,max and a controls the shape of the function. The empirical relationship between Ψ₅₀ and a from Christoffersen et al. [51] described in electronic supplementary material, appendix S2, reduces the plant hydraulic parameters needed in SOX to only Ψ₅₀, 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[t−1]) to compute the current timestep Ψ_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 and Ψ_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 Ψ drop along the tree, computing k_rc as:

In SOX we represent the gradual Ψ drop along the tree using the middle value of the root–canopy gradient (Ψ_c,mid):

Using f(Ψ_c,mid) is numerically simpler than equation (2.6) and provides similar results within a realistic range of Ψ_c,pd and Ψ_c (electronic supplementary material, figure S1). The k_rc produced by f(Ψ_c,mid) normalized as a function of k_rc,max, giving k_cost, represents the costs of stomatal aperture in SOX:

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 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 Ψ_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 Ψ_r ≈ Ψ_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[t−1], 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 Ψ_s from the site observations of root mass-weighted soil moisture content (θ, m³ m⁻³), 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² leaf m⁻² 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² leaf m⁻² soil).

We compare SOX agreement with observations against a model that uses a drought representation model based on the β-function (β_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⁻² s⁻¹) predicted by SOX to the unified stomatal optimization model (USO) of Medlyn et al. [40], described in electronic supplementary material, appendix S6.

(c) El Niño simulations during the twentieth century

We compared SOX's sensitivity to drought events with the β_fun model (electronic supplementary material, figures S4) using meteorological and vegetation hydraulic observations coupled to the modelled soil moisture dynamics of three Amazonian sites (electronic supplementary material, figure S3). We used the CRU-NCEP (v. 4, see supplementary figure S5; N. Viovy, Laboratoire des Sciences du Climat et de l’Environnement, (LSCE), France, 2016, personal communication) 6-hourly meteorological data from 1901 to 2016 to drive our models (see electronic supplementary material, figure S5). These forest sites possess distinct climatic responses to El Niño events (electronic supplementary material, figure S3), represented by the Niño-3 index, which is calculated as the mean sea surface temperature (SST) anomaly from 5°N to 5°S and 150–90°W [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 (Ψ_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 Ψ_s; and Sim4 is without anomalies in any of the previously mentioned variables (T_a, D and Ψ_s, see electronic supplementary material, table S3 for summary).

3. 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, Ψ_s has a much stronger effect in plants more vulnerable to cavitation (figure 2; electronic supplementary material, figure S6).

• Download figure
• Open in new tab
• Download PowerPoint

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 cavitation-resistant 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, hydraulically-induced stomatal closure.

The SOX response to atmospheric demand (figure 2g,h) results from the increased k_rc loss associated with the lower Ψ_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 induces an exponential decline in g_c as the c_i,opt shifts closer to 0 (figure 1). This pattern is commonly observed [66–68], and more accentuated in plants with a less negative Ψ₅₀ owing to their increased cavitation costs as the atmosphere dries (figure 2g,h). The β_fun predicts a more gradual stomatal closure in response to D, which approximates SOX predictions for Ψ₅₀ = −2.5 MPa when D > 0.02 mol mol⁻¹(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 necessary to maintain A when Ψ_s is low. SOX predicts a more gradual decline in g_c in response to drying soil when compared with the β_fun model (figure 2i,j). The β model predicts g_c = 0 when Ψ_s = −1.5 MPa and θ = θ_w (electronic supplementary material, figure S4; figure 2i,j), whereas SOX predicts that even plants very vulnerable to cavitation (Ψ₅₀ = −1 MPa) still have 30% of their maximum g_c at the same Ψ_s. Accounting for changes in k_sr with equations in the electronic supplementary material, appendix S3 makes g_c more responsive to Ψ_s, affecting particularly plants more resistant to cavitation. At Ψ_s = −5 MPa, even plants with Ψ₅₀ = −5 MPa will have dropped to 1% of their maximum g_c, whereas in the model that assumes Ψ_r ≈ Ψ_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 response to Ψ_s, which is linear when a is low, but higher a produces a more stable Ψ_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 β_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. [51] to estimate the shape parameter of the vulnerability curve (a), which predicts a = 2.1. Optimizing a to the observations of each plot produces a very high agreement on the control plot (a = 2.4; R = 0.94) and a strong agreement on the TFE plot (a = 1.1; R = 0.44). Accounting for changes in k_sr (electronic supplementary material, appendix S4) allows us to improve even further the agreement between SOX predictions and observations in the TFE plot (R = 0.5). SOX can also reproduce well the observed seasonal fluctuations in Ψ_c (electronic supplementary material, figure S8). The β_fun model greatly overestimates the soil moisture effects, leading to excessive stomatal regulation in the control treatment (R = −0.3) and almost complete stomatal closure in the TFE (figure 3b; electronic supplementary material, figure S8). A model that ignores soil moisture effects (β_off) can fit the control plot data (R = 0.94) but cannot capture the seasonality in the TFE plot (R = 0.11).

• Download figure
• Open in new tab
• Download PowerPoint

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₁ and its response to Ψ_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₀, to 0. Therefore, the SOX A–g_w relationship implies a dependency of the water marginal carbon costs (related with the USO parameter g₁, see electronic supplementary material, appendix S6) on the light regime that is not present in USO. The g₁ 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 with the USO parameters estimated from the β model's output (electronic supplementary material, table S4).

• Download figure
• Open in new tab
• Download PowerPoint

(c) El Niño predictions during the twentieth century

Tapajós was the only site where JULES predicted significant soil drought, which could be particularly intense in El Niño years (electronic supplementary material, figures S3, S5 and S9). At this site, the β_fun model is oversensitive to soil drought, strongly limiting A (electronic supplementary material, figure S10) in a similar way to what is observed in figure 3. The β_fun model is also more sensitive to soil drought anomalies, as shown by the greater interannual variability between Sim1 and Sim3 in Tapajós (figure 5). Both models produce a similar magnitude of negative responses to soil drought anomalies of ca −0.7 kg C m⁻², but β_fun predicts that A can rise by up to 0.52 kg C m⁻² in years when the soil is more humid than usual, while SOX predicts a maximum increase of 0.16 kg C m⁻² yr⁻¹ (figure 5e). This divergence amplifies over the years, leading to the cumulative effect of soil drought anomalies in Tapajós predicted by the β_fun model being −0.49 kg C m⁻² after 115 years, while SOX predicts a strong negative cumulative effect of −5.37 kg C m⁻².

• Download figure
• Open in new tab
• Download PowerPoint

The effect of atmospheric anomalies is comparatively small in Tapajos (figure 5d,e), but is the dominant effect in Caxiuanã and Manaus (figure 5a–c,g–i). Atmospheric anomalies tended to increase A until ca 1950, with β_fun predicting a maximum effect of accumulated anomalies of 2.1 kg C m⁻² in 1947 at Manaus, whereas SOX predicts only 0.76 kg C m⁻² at the same year (figure 5g). The increase of frequency and magnitude of positive climatic anomalies in the second half of the twentieth century (electronic supplementary material, figure S9) had a detrimental effect on forest A, particularly strong in Manaus. The β_fun model predicts that at the end of the 115 years climatic anomalies would reduce forest A by 0.85 kg C m⁻², while SOX predicts a reduction of 0.92 kg C m⁻². The responses of Caxiuanã to climatic anomalies are similar to Manaus but less pronounced, with an overall cumulative effect of climatic anomalies of −0.62 kg C m⁻² according to β_fun and −0.15 kg C m⁻² by SOX (figure 5c).

4. Discussion

Our results show that a xylem hydraulics-based stomatal optimization scheme can produce realistic stomatal responses to environmental variables (figure 2), being able to predict the observed responses of a tropical forest to seasonal, and even severe experimentally-induced soil drought (figure 3). 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 β_fun model (figures 4 and 5; electronic supplementary material, figure S4). The drastic differences that emerge from long-term simulations between SOX and the β_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 β_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 β_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, β_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 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 hydraulics-based 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–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 time-scales 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–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–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. A−Eλ), 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 Ψ_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/semi-empirical functions [20,24,72] is highlighted in the long-term 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 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.5° × 0.5° 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 Ψ_s variability (electronic supplementary material, figure S3), and it was the site where the divergences between the β_fun model and SOX were largest (figure 5; electronic supplementary material, figures S9 and S10). The β_fun excessive soil moisture response and highly variable response to climatic anomalies reflect the steep gradient between the critical and wilting points of the β_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 β_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.

5. 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. ([email protected]). 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.

Authors' contributions

C.B.E., L.R., S.S., P.C., M.M., R.S.O. and A.D.F. designed the model. L.R., R.S.O., P.R.L.B., F.V.B., A.C.L.d.C. and P.M. collected the data used to parameterize and evaluate the model. C.B.E. wrote the model and the manuscript with help from all authors.

Competing interests

We declare we have no competing interests.

Funding

This study was funded by the Newton Fund through the Met. Office Climate Science for Service Partnership Brazil (CSSP Brazil), UK NERC independent fellowship grant no. NE/N014022/1 to L.R., UK NERC grant no. NE/J010154/1 to S.S., UK NERC grant no. NE/J011002 to S.S., P.M. and M.M., ARC grant DP170104091 to P.M., and CNPQ grant no. 457914/2013-0/MCTI/CNPq/FNDCT/LBA/ESECAFLOR to A.C.L.d.C.

Footnotes

One contribution of 22 to a discussion meeting issue ‘The impact of the 2015/2016 El Niño on the terrestrial tropical carbon cycle: patterns, mechanisms and implications’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4224362.