On the Treatment of Soil Water Stress in GCM Simulations of Vegetation Physiology

Current land surface schemes in weather and climate models make use of the so-called coupled photosynthesis–stomatal conductance (A–g s ) models of plant function to determine the surface fluxes that govern the terrestrial energy, water and carbon budgets. Plant physiology is controlled by many environmental factors, and a number of complex feedbacks are involved, but soil moisture control on root water uptake is primary, particularly in sub-tropical to temperate ecosystems. Land surface models represent plant water stress in different ways, but most implement a water stress factor, β, which ranges linearly (more recently also curvilinearly) between β = 1 for unstressed vegetation and β = 0 at the wilting point, expressed in terms of volumetric water content ( θ ). β is most commonly used to either limit A or g s, and hence carbon and water fluxes, and a pertinent research question is whether these treatments are in fact interchangeable. Following Egea et al. (Agricultural and Forest Meteorology, 2011, 151 (10), 1,370–1,384) and Verhoef et al. (Agricultural and Forest Meteorology, 2014, 191, 22–32), we have implemented new β treatments, reflecting higher levels of biophysical complexity in a state-of-the-art LSM, Joint UK Land Environment Simulator, by allowing root zone soil moisture to limit plant function non-linearly and via individual routes (carbon assimilation, stomatal conductance, or mesophyll conductance) as well as any (non-linear) combinations thereof. The treatment of β does matter to the prediction of water and carbon fluxes: this study demonstrates that it represents a key structural uncertainty in contemporary LSMs, in terms of predictions of gross primary productivity, energy fluxes and soil moisture evolution, both in terms of climate means and response to a number of European droughts, including the 2003 heat wave. Treatments allowing ß to act on vegetation fluxes via stomatal and mesophyll routes are able to simulate the spatiotemporal variability in water use efficiency with higher fidelity during the growing season; they also support a broader range of ecosystem responses, e.g., those observed in regions that are radiation limited or water limited. We conclude that current practice in weather and climate modelling is inconsistent, as well as too simplistic, failing to credibly simulate vegetation response to soil water stress across the typical range of variability that is encountered for current European weather and climate conditions, including extremes of land surface temperature and soil moisture drought. A generalized approach performs better in current climate conditions and promises to be, based on responses to recently observed extremes, more trustworthy for predicting the impacts of climate change.

Current land surface schemes in weather and climate models make use of the so-called coupled photosynthesis-stomatal conductance (A-g s ) models of plant function to determine the surface fluxes that govern the terrestrial energy, water and carbon budgets. Plant physiology is controlled by many environmental factors, and a number of complex feedbacks are involved, but soil moisture control on root water uptake is primary, particularly in sub-tropical to temperate ecosystems. Land surface models represent plant water stress in different ways, but most implement a water stress factor, β, which ranges linearly (more recently also curvilinearly) between β 1 for unstressed vegetation and β 0 at the wilting point, expressed in terms of volumetric water content (θ). β is most commonly used to either limit A or g s , and hence carbon and water fluxes, and a pertinent research question is whether these treatments are in fact interchangeable. Following Egea et al. (Agricultural and Forest Meteorology, 2011, 151 (10), 1,370-1,384) and Verhoef et al. (Agricultural and Forest Meteorology, 2014, 191, 22-32), we have implemented new β treatments, reflecting higher levels of biophysical complexity in a state-of-the-art LSM, Joint UK Land Environment Simulator, by allowing root zone soil moisture to limit plant function non-linearly and via individual routes (carbon assimilation, stomatal conductance, or mesophyll conductance) as well as any (non-linear) combinations thereof. The treatment of β does matter to the prediction of water and carbon fluxes: this study demonstrates that it represents a key structural uncertainty in contemporary LSMs, in terms of predictions of gross primary productivity, energy fluxes and soil moisture evolution, both in terms of climate means and response to a number of European droughts, including the 2003 heat wave. Treatments allowing ß to act on vegetation fluxes via stomatal and mesophyll routes are able to simulate the spatiotemporal variability in water use efficiency with higher fidelity during the growing season; they also support a broader range of ecosystem responses, e.g., those observed in regions that are radiation limited or water limited. We conclude that current practice in weather and climate modelling is inconsistent, as well as too simplistic, failing to credibly INTRODUCTION Water availability exerts a major control on vegetation gross primary productivity (GPP), as well as on the land surface energy balance. It has been estimated that ∼40% of the global vegetated land surface, particularly in sub-humid, semi-arid, and arid regions (Stocker et al., 2018;O'Sullivan et al., 2020), experiences plant activity/growth limitations caused by seasonal water deficits (Nemani et al., 2003;Beer et al., 2010). In the context of future projections of ecosystem response, soil moisture stress is predicted to increase over large regions (Berg et al., 2016;Ukkola et al., 2020). Consequently, there is a clear requirement to incorporate accurate, process-based models of plant response to soil moisture stress in coupled land-atmosphere climate models. However, the models currently used to represent biogeophysical and biogeochemical processes in Earth System Models, or even simpler GCMs, are often unable to properly capture observed responses to soil moisture stress (e.g., Beer et al., 2010;Powell et al., 2013;Medlyn et al., 2016;De Kauwe et al., 2017;Restrepo-Coupe et al., 2017;Peters et al., 2018;Paschalis et al., 2020).
Plants respond to reductions in soil moisture content (SMC) through a range of drought tolerance and prevention strategies; a thorough review of the state of our knowledge of these processes is provided in Harper et al. (2021). The immediate response is to reduce physiological activity, which has consequences for primary production and for transpiration. Land-atmosphere feedbacks involving anomalously high near-surface vapor pressure deficit and leaf boundary temperature (see e.g., Ball et al., 1987), further exacerbate concurrent soil drought and atmospheric aridity (Zhou et al., 2019). Such conditions are more likely when meteorological drought occurs, often as a result of stagnant atmospheric conditions (e.g., summertime blocking). Such events see a reduction in carbon uptake, with possible consequences for plant growth and below-ground carbon allocation, but also increase the Bowen ratio, raise surface temperature and can lead to further dessication of the soils, at a Clausius-Clapeyron rate (see for instance the review in Seneviratne et al., 2010 andVargas Zeppetello et al., 2019).
Most land surface schemes do incorporate the process of downregulation of photosynthesis, or of stomatal conductance, but this is mostly done in a simplistic way and with a macroscale approach; a review of the range of complexity is available in Verhoef and Egea (2014). A typical LSM represents the regulation of stomatal conductance as a simple generic function of SMC, generally expressed in terms of volumetric water content (θ, m 3 m −3 ). This simple generic function is the so-called "beta" function, where β is a factor between zero and one that limits photosynthesis in some way (this depends on the LSM, see Section 2). Above a critical SMC, θ c , there is no stress (β 1), and below the critical threshold value, stress increases as SMC decreases, until the wilting point, θ w , is reached (β 0). Alternative, yet related, expressions are available whereby stomatal regulation occurs through changes in the soil matric potential, ψ (a measure of how tightly the water is held in the soil pores, thereby affecting water uptake by the roots), expressed in pressure units, such as MPa. θ and ψ are closely related, via the water retention curve, and some models emulate a ψ-type parametrization via a curvilinear dependence on θ (see more details in Verhoef and Egea, 2014). It is important to note, in this context, that the widely adopted linear relationship in most LSMs simulates unrealistically low plants resiliency to water stress in drought conditions (Niu et al., 2020, as well as the discussion in; Verhoef and Egea, 2014).
A recent survey of land surface schemes currently in use in the land surface processes community (available in Peters et al., 2018) has indicated that three strategies, or pathways, exist: 1) those that impose soil moisture stress by regulating stomatal conductance (e.g., LPJ-GUESS, LPJ-C13, CLM); 2) those that impose soil moisture stress by downregulating photosynthesis (e.g., SiB2, CLM, JULES); 3) those that employ some form of control on mesophyll conductance (e.g., SIBCASA, ORCHIDEE) and 4) those that employ all strategies at once (e.g., SIBCASA, ORCHIDEE). Table 1 provides a summary of the above and suggests that there is currently no community consensus on what approach is to be used. At the same time, Table 1 illustrates that individual research groups have chosen the plant water stress strategy that best suited their specific scientific objective (e.g., the simulation of GPP), rather than choosing an approach that considers land surface fluxes and processes in an interconnected way.
A different approach to modelling the stomatal response to drought in LSMs is that based on plant hydraulics modelling (Eller et al., 2020;Sabot et al., 2020). This approach offers a promising and mechanistic alternative to the empirical β function approach, but it still has some limitations for global and long-term modelling related to model parameterization and limited knowledge on plant traits plasticity and acclimation (Anderegg and Venturas, 2020).
The main research question in this paper is what impacts those a-priori decisions on the "β pathway" have on the concurrent prediction of the surface energy balance, of GPP and surface temperature. A secondary question is whether, via triggering feedbacks involving near-surface atmospheric conditions, the choice of β pathway can alter the simulation of weather and climate trajectories and, ultimately, have consequences for ESM projections of climate change, particularly with regards to extremes such as heatwaves.
This study aims therefore to comparatively evaluate the impact of the pathway between soil moisture stress and vegetation function on the simulation of GPP and latent heat flux (LE) for a range of biomes and climates in Europe, under controlled conditions, enabling a direct comparison of the different strategies. A further focus is the ratio of leaf-internal to ambient CO 2 (C i /C a , here denoted as χ), which relates the assimilation rate to stomatal conductance (see for instance Prentice et al., 2014;Dewar et al., 2018), as 1) it is observable, even at the large scale, via the measurement of various isotopes; 2) it can be used as a proxy for water use efficiency (see for instance Peters et al., 2018) to confront model predictions; and 3) it has been shown to be (weakly) dependent on atmospheric CO 2 concentration (e.g. Lavergne et al., 2020), so a more robust indicator of long-term vegetation function. For all these reasons, χ must thus be thoroughly considered as a key variable in the study of anthropogenic climate change, but instead, it has been shown to be poorly represented by current LSMs (see the discussion in Peters et al., 2018), particularly in terms of (continental-scale) interannual variability.
The approach to the investigation is to adopt a generalized, agnostic β approach in a state-of-the-art LSM, driven by observed meteorology over a large region, and to focus on simulation fidelity in terms of climate means and variability, including the response to recent heat waves, particularly the extreme heat wave of 2003. The study is organized as follows: Introduction re-visits the Egea et al. (2011) solutions for aspects that matter to coupling to a full Land Surface Model (LSM); Methods describes Methods and Data used in the modelling study; Results presents results for climate means (3.1), for a composite of European droughts in the 1979-2011 period (3.2) and for the extreme growing season of 2003 (3.3), to then move into Europe-wide water use efficiency, first at the large scale level (3.4), next a mechanistic analysis at the stomatal and PFT level (3.5-3.6), to re-emerge to a scaling up of the results at European level, based on the isotopic fractionation Δ metric from Peters et al. (2018) in 3.7. Discussion discusses the implications of this study for climate studies, revisits Peters et al. (2018) from the perspective of model development, and suggests further avenues for progress, particularly from the point of view of more efficient numerical solutions; Summary and Conclusion provides summary and conclusions.

Idealised Analytical Solutions
The study by Egea et al. 2011, (E11 hereafter) developed and exploited analytical solutions for the simultaneous equation set that governs the fast, dynamic evolution of assimilation (A), stomatal conductance (g s ) and leaf internal CO 2 concentration (C i ), the more dynamic component of χ, observed during the diurnal cycle. Numerical solutions were based on an approach by Baldocchi (1994), and eliminated the need for iteration, while enabling complete freedom in the treatment of β. In order to understand the results in this study, which moves beyond E11 in implementing the framework into a full land surface simulation scheme, offered by the JULES LSM, the E11 equations are first revisited in a highly idealized manner. This is a useful exercise, because in coupling E11 to a land surface scheme such as JULES, some feedbacks with the canopy environment become possible, at least in part (when JULES is used offline, as driven by observed meteorology) or more completely (when JULES is within the parent HadGEM3-GC31 GCM), clouding our interpretation of the primary chain of mechanisms and of the feedback responses. This preliminary step is thus important, because it enables the reader to better appreciate the motivation of the study, as well as to build expectations for what understanding can be achieved at this stage, prior to coupling to a complex land surface scheme, or even to a GCM. Figure 1 provides a conceptual-level prediction of the shape of the functional relationships between soil moisture availability and land surface prognostic variables, when applying the E11 framework to the coupled A-g s model contained in an LSM such as JULES. The relative stomatal conductance, g s /g s (β 1), is shown on the x-axis, while relative assimilation, A/A(β 1), is shown on the y axis, where g s (β 1) and A(β 1) are the unstressed stomatal conductance and assimilation, respectively. If no other stress, or any feedbacks, are imposed, the relationship should be linear, as in the typical Ball-Berry plot, and that linear relationship is indicated by the black straight line. The three panels in Figure 1 show, instead, the responses obtained in the E11 framework: by imposing soil moisture stress (β) on the stomatal conductance pathway (left panel, blue curve), a curvilinear relationship is revealed, and this shape generally implies a higher water use efficiency (WUE) at intermediate β levels; the other two pathways (biochemical and mesophyll) show no deviation from the linear relationship. Adding an idealized vapor pressure deficit (VPD) feedback term to E11 (center panel, this is done by editing the leaf-level atmospheric input file to the E11 model that increases VPD in proportion to 1/β), causes no change in the stomatal pathway, but now causes the two non-stomatal pathways to also become curvilinear, and to partially emulate the stomatal pathway, with a higher WUE, maximized near the β mid-range. The addition of a term that emulates (leaf) temperature feedback (right panel, this is done by increasing leaf temperature linearly by a factor proportional to 1/β) further increases WUE (right panel) for all three pathways. The rightmost panel is the one that most resembles field observations under conditions of vegetation stress, e.g., what is shown in E11 for some agricultural crops in southern Spain. Figure 2 shows the response of χ to relative stomatal aperture (as directly caused by the soil moisture stress) for the three pathways, under the same conditions imposed for Figure 1: when no feedback with the atmosphere is allowed (left panel), any of the non-stomatal routes show a very high value of χ under stressed conditions, which results from the suppression of assimilation, as directly caused by β. This is not observed in the stomatal pathway, because stomata close, but assimilation carries on unhindered by β, using CO 2 inside the stomatal chamber. The three pathways behave identically for unstressed vegetation. When a VPD feedback is allowed (center panel), the two non-stomatal pathways show a reduction of χ in the lower range of stomatal conductance, while the stomatal pathway solution resembles the one in the leftmost panel. The introduction of a leaf temperature feedback (right panel) has different effects on the biochemical and stomatal pathways: the former returns to higher levels of χ for low stomatal conductance conditions (strong soil moisture stress), while the latter shows a reversal at low stomatal conductance levels, with χ increasing. The latter is the response that most strongly resembles the field   Egea et al. (2011) model: relative stomatal conductance (g s /g s (β 1), x axis) versus internal CO 2 concentration, C i (y axis, expressed as Eq. 5) for the three pathways and for zero feedback with the canopy air environment (left panel), VPD feedback with the canopy air environment (center panel) and VPD plus leaf temperature feedback with the canopy air environment (right panel).
Frontiers in Environmental Science | www.frontiersin.org August 2021 | Volume 9 | Article 689301 observations in E11, as well as those presented in Yan et al. (2017). These exploratory analyses strengthen the need for an investigation of the E11 under conditions that more closely resemble those of a GCM, albeit without incurring the additional complexity of GCM errors in the surface energy balance, which would hinder our understanding.

METHODS
The JULES Land Surface Component of HadGEM3-GC31 JULES (the Joint United Kingdom Land Environment Simulator) is based on the Hadley Centre land surface scheme MOSES2. Full descriptions of JULES are provided in Best et al., 2011;Clark et al., 2011. JULES simulates the exchange of water, momentum and energy between the soil, land surface, and atmosphere. It is driven by sub-daily (typically 3-hourly) time series of radiation, precipitation, temperature, humidity, wind speed and surface pressure. The soil is divided into layers, with the thermal and hydraulic characteristics defined for all layers. It is important for the interpretation of the experiments presented in this study to note that the soil column in each grid box is vertically resolved by four layers of increasing thickness, down to a total depth of 3 m. Land surface tiles share a single soil column, and soil moisture is thus shared across tiles in each grid box.
In its standard configuration, JULES divides the land-surface into nine surface types: broadleaf trees, needle leaf trees, C3 (temperate) grass, C4 (tropical) grass, shrubs, urban, inland water, bare soil, and ice. Crops are treated as grasses. Sub-grid heterogeneity is represented by tiling of land-surface types (for example, Essery et al., 2003). All land grid boxes can be made up of any mixture of the nine land-surface types, except ice. The surface fluxes of moisture and heat are computed individually for each tile, and the state of the grid box is prognosed via the aggregation of tiles fluxes.
The biophysical state of each vegetation tile is defined by its leaf area index (LAI), canopy height and rooting depth. In the JULES experiments presented in this study a seasonal vegetation phenology, based on a MODIS climatology, is imposed at each grid point, so that no aspect of dynamic vegetation is enabled; this conscious decision imposes a buffering effect on the perturbation experiments, but makes it easier to interpret the sensitivity to β pathway, also because the dynamic vegetation aspects of JULES have been tuned in the past for the biochemical β route.
The surface fluxes of moisture and heat in JULES depend on the atmospheric boundary conditions and the characteristics of the land surface. The fluxes are computed by a network of parameterizations including soil surface, stomatal and aerodynamic resistances. Stomatal resistance also controls the intake of CO 2 , and is thus the link between the carbon, water and energy cycles. The stomatal resistance parameterization (Cox et al., 1998) depends on environmental conditions, including the ambient concentration of CO 2 .
JULES calculates the photosynthesis rate using the method described in Collatz et al. (1991) for the C3 pathway and Collatz et al. (1992) for the C4 pathway. Potential photosynthesis takes place at the minimum of three limiting rates: light; enzyme kinematics (Rubisco); and the transport of photosynthetic products. The potential photosynthesis is reduced under water-stressed conditions by β, the soil moisture availability factor. Photosynthesis is scaled from the leaf to the grid box scale under the assumption that the rate of photosynthesis is a function of the LAI.

The Implementation of a Flexible Treatment of β Into JULES
The three equations that must be solved simultaneously in order to produce land surface prognostics for temperature, moisture, GPP, and relative fluxes to the atmosphere are: where V cmax is the maximum carboxylation rate (μmol m −2 s −1 ) and J is the electron transport rate (μmol m −2 s −1 ), the K c,o are the Michaelis-Menton kinetic factors for carboxylation and oxygenation (Pa), Γ* is the chloroplastic CO 2 photocompensation point in the absence of mitochondrial respiration (μmol mol −1 ) and C c is the chloroplastic CO 2 concentration (μmol mol −1 ).
In Eq. 2 g sc (mol m −2 s −1 ) is the stomatal conductance to CO 2 , g 0 (mol m −2 s −1 ) is the cuticular conductance, C s (μmol mol −1 ) is the CO 2 concentration at the leaf surface, Γ is the CO 2 compensation point (μmol mol −1 ), D s is the vapour pressure deficit at the leaf surface, m is equal to 1/(1-f 0 ) and D* is D max /(m-1); f 0 and D max are the parameters defined by Jacobs et al. (1996).
In order to solve the simultaneous set above, the three equations that describe A, g sc and g m as a function of C c were coupled together and solved for C c , resulting in a cubic relationship, which was solved using the Baldocchi (1994) method. This approach is very flexible and removes the requirement for iteration, which is computationally inefficient. The mathematical details are fully explained in the Supplementary Material section. Similar to the approach in E11, the scheme was first tested at a small selection of FLUXNET sites (not shown), e.g., for accuracy and for numerical stability, but then applied to the regional distributed simulations shown here.
In order to find solutions for the A-g s sub-model, enabling soil moisture stress to impact the biochemical pathway (assimilation) the conductance pathways (stomatal, mesophyll), or any combinations thereof, the functional relationships, as in E11, are: August 2021 | Volume 9 | Article 689301 where β, as in E11, can be applied to any single pathway, or combinations, or to all, with different shapes (linear or curvilinear) and weights (see E11 for a large selection of tests). The subscripts i,j B, S and M correspond to biochemical, stomatal, and mesophyll limitations, respectively (see Table 2).
The general expression for β, also enabling non-linear dependence, is, as in E11: where β i is a soil moisture dependent function ranging from 1 to 0. The subscript j for the q j exponent signifies that independent functional shapes can be adopted in the definition of each β i , enabling complete freedom for each pathway (see Table 2); for further details, see the discussion, after Eq. 16, in E11. Following E11, in order to test the relative influence of soil moisture stress on canopy-atmosphere exchange processes, we have set up Europe-wide JULES model experiments with the following four model configurations, constituting four sensitivity experiments.
Configuration QB fully corresponds in its design (application to photosynthetic capacity) to that of the original JULES, and could in principle be considered our control experiment. We have chosen, however, to include an unmodified JULES integration (CTL) in our analysis, as it was deemed useful in order to understand the impact of using the analytic solution method in E11 (see the full development in Additional Materials) and to quantify any uncertainty. Configuration QS applies soil moisture limitation stress to stomatal closure exclusively; configuration QM applies soil moisture limitation stress to the mesophyll conductance exclusively; C6, finally, applies soil moisture limitation stress to all three pathways, albeit with differing weights and functional shapes (see the motivation for these choices in E11).

Experimental Design
We have configured a Europe-wide (35N-65N, 10W to 40E) JULES (a full set was run with v2.2 initially, but for all results presented here we have re-run with v4.4, Clark et al., 2011), driven by meteorological data (net radiation, precipitation, wind, reference height temperature and humidity; WATCH, Weedon et al., 2011). The model grid configuration, in the horizontal and vertical, is identical to that in Clark et al. (2011); model parameters, in particular those controlling photosynthesis, are also standard. Each grid point comprises 9 "tiles", five of which are plant functional types (PFTs), the others are ice, bare soil, urban. The model is run continuously, at half-hourly time steps, for 33 years, from 1979 to 2011. We apply an initial spin-up to each simulation, to create 1979 balanced initial conditions for soil prognostics; spin-up completion is governed by convergence of soil moisture and soil temperature solution (to within 1% between spin-up cycles). Typically, each model goes through 30 spin-up cycles of 10 years each, until full convergence, for a maximum total of 300 years.
The E11 scheme was implemented as an alternative to the standard A-g s scheme (see Cox et al., 1998), albeit only for C3 plants at this initial (developmental) stage. The abundance of C4 plants in the European domain chosen for this study is minimal (mostly less than 10%, with only a single grid point found to be dominant in Turkey) and any grid points containing a fraction of C4 plants above 15% are discarded from the analysis presented in this paper. This approach is fully consistent with our forcing of the model using observed meteorology, by which individual grid points are fully independent from each other (and tiles partially so, because of the shared within-soil solutions), so that any tiles potentially "contaminated" by C4 biophysics cannot contribute to the evolution of fluxes in any neighboring grid cells. This is because in the off-line setup there is no horizontal advection of air between grid points. Future applications with the coupled GCM, in which air can move freely between grid points, will require development of an equivalent β pathway treatment for C4 plants.
Atmospheric CO 2 , thus C a C s were fixed at 360 μmol mol −1 for the period.

Special High-Frequency Stomatal-Level Diagnostics
New stomatal level diagnostics for C i (also for leaf-level VPD, not shown here), for each PFT, were introduced in JULES, in order to study the evolution of χ in high detail (time-step level), and to enable intrinsic Water Use Efficiency (iWUE) and isotopic discrimination metric (Δ, ‰) diagnostics, see 2.6, to be used for comparison with FLUXNET observations in the diurnal and seasonal cycle of C i (χ). The C i equation is: The equation used to estimate Δ from C i is, as in Peters et al. (2018): where Δ p (27‰) and Δ d (4.4‰) are the isotopic discriminations during assimilation catalysed by the enzyme Rubisco in C3 photosynthesis and molecular diffusion of CO 2 through the stomata, respectively.
iWUE ≈ c a · 1 − Ci Ca 1.6 (7) These time-step level diagnostics of C i and Δ, thus iWUE, were next resampled to 3-hourly frequency for each PFT (5 in total in this study), GPP-weighted as in observations, and for each grid box in the domain. As such, the output is large, and more detailed than what can be achieved in a GCM setting. The same highfrequency diagnostics, for simulation QB, were used in Peters et al. (2018). The ecosystem-level water use efficiency, eWUE, was computed using the ratio of NPP and evapotranspiration from grid-box daytime averages, as represented by 6-hourly model statistics, then projected onto monthly means.

Observational Datasets Used for Model Assessment
The FLUXCOM data were obtained from Jung et al. (2019) and are used as monthly means over 0.5 o grid boxes. While we realize that FLUXCOM are not actually observations, it is often assumed, for evaluation purposes, to be a suitable extrapolation of the measured net ecosystem exchange from many eddy-covariance sites across Europe. This also enables the study to be verified with a relatively independent dataset (the version chosen upscales the FLUXNET observations by using neural networks to combine it with remote sensing and meteorological data, in this case also WFDEI). This is the best compromise achievable in bridging the field scale (FLUXNET towers) to the large scale used in this study, without direct use of re-analysis products, which are strongly model-dependent in terms of surface fluxes. The measurements used as reference for iWUE are derived from FLUXNET data, as described in Peters et al. (2018), Section S5, and summarized in their Table S3. To compare these two better, relative changes in ‰ from the baseline (years 2002-2006) are shown in Results. Averaging of C i , thus iWUE, over time was also GPP-weighted, such that high-photosynthesis hours are represented the most.
Definition of JULES Soil Moisture Droughts for Masking of Regional Flux Averages JULES soil moisture droughts are defined by the seasonal average (AMJJAS) of the linearly calculated β anomaly (Δβ) at each grid point, as deviations from the 1979-2011 climatology in CTL (enabling a common comparison). A regional mask is constructed from these Δβ for each year through a three-step process, using a quantile approach. At first, drought thresholds were calculated for each grid cell and each month as the 20th percentile of daily β values using the CTL simulation. Subsequently, differences between monthly averages of β and monthly drought thresholds were summed up over the period April to September for each grid cell and each year, and grid cells were marked as dry for a particular year when this sum was negative. Integrating over the period April to September, which was selected based on GPP values, allows to capture droughts of different intensities and durations, i.e., short, intense droughts as well as prolonged, weaker droughts. Finally, drought masks were created by selecting the largest connected area identified as dry within the model domain for each year. These drought masks were applied to all simulations, to guarantee consistency in spatial averages across simulations.

RESULTS
All experiments (CTL, QB, QS, QM, C6) were spun-up individually and run for the full period. However, in terms of presentation, in order to reduce the number of plots and subplots, CTL is only shown in terms of long-term soil moisture memory, as QB is virtually identical to CTL in terms of landatmosphere fluxes. First the climatological biases of some variables of interest will be shown, then the 2003 heat wave response will be used as a surrogate for future climate conditions, as suggested for instance in Schär et al. (2004), where it was shown how the 2003 event represents an outlier for current climate, but is, instead, compatible with a temperature probability distribution function (PDF) extracted from climate projections forced by a strong climate change scenario. For completeness, we have also included a more general assessment of all major droughts in the period.

Climatological Biases
All biases have been computed relative to FLUXCOM products (Jung et al., 2019). Similar plots have also been produced with alternative products in the ILAMB verification toolkit, but are not shown here, as they would not bring extra information. The focus for the fluxes is on the growing season (MAM + JJA), as autumn and wintertime responses to the perturbation experiments are of much smaller amplitude, but soil moisture will be shown over a more extended period, because of soil moisture memory effects, which are of interest. Figure 3 shows the GPP climatological biases for experiments QB, QS, and C6, which are the most salient. Overall, the figure indicates that the QB model (nearly identical to CTL, not shown) tends to overestimate European GPP in spring, particularly in central Europe, which is made worse, and more widespread, by QS and C6. In summer, following an overactive use of soil moisture (mostly by vegetation) at the center of the domain, the GPP tends to be underestimated by QB, as well as, in the North and South of the domain, by QS and C6; there is a slight overestimation at the center of the domain, and over the Alps. The summer biases are sufficient to overwhelm the annual mean. Figure 4 shows the LE response, which, to a large extent, reflects the transpiration response, given the high vegetation fraction for these large-scale experiments. For the three models, QB, QS, and C6, latent heat tends to be overestimated in Spring nearly everywhere in Europe, while for Summer a dipole emerges, with a slight overestimation in North and Central Europe, but a pronounced underestimation in Southern Europe, very likely due to soil dessication. No significant differences are seen, however, between QB and C6 in any of the seasons, which is an important result in view of what has been found for the GPP response. Model QS shows some indication of a bias reduction in the NE portion of the domain, but a worsening in the SE (over Greece and Turkey). Figure 5 shows the sensible heat flux (SH) response, which is nearly insensitive to the model formulation in spring, and shows small regional responses in summer, which mirror the responses seen in the latent heat flux maps (Figure 4): this small response to model formulation is again rather expected (because the model was run with observed meteorology and is thus constrained) and Frontiers in Environmental Science | www.frontiersin.org August 2021 | Volume 9 | Article 689301 rather reassuring, as in most years European plants are not water stressed for periods of time long enough to impact climate means. Figure 6 shows the available soil moisture response in terms of the β factor in the soil: since it is a relative measure, and weighted by layer thickness, it reflects the effect on the bulk of the roots for the JULES vegetation (mostly forests and C3 grasses in this European domain), with a possible buffering effect from the deepest layer, which is more accessible to trees than grasses. The response of each model is quite strikingly different: while QB is virtually identical to CTL (some differences are expected, since these are long-term spun-up experiments, and any numerical residual would show strongly in a cumulative variable such as soil moisture), the QM response is mostly a strong drying, seen in all seasons, but with a strong North-West to South-East gradient. The response of QS is, instead, a moistening, again with a strong NW to SE gradient. The C6 response is, remarkably, nearly identical to the QB (CTL) response, indicating that this combined pathway for soil moisture impact on the surface    Figure 7, it is next possible to produce a synthesis plot that presents the response to the treatment of β, for each year in which a significant drought has been identified. Figure 8 provides a synthesis of model behavior for 20 soil moisture drought years, indicated by dots of different magnitude, according to spatial extent and intensity, expanding from the focus on the largest disturbance, 2003, which remains a prominent benchmark in terms of area extent and intensity. Figure 8 shows the summary of model flux errors, averaged over each drought year, extracted from the data used for Figures 3-5 (and defined as difference from identically processed FLUXCOM observations), against the annual AMMJAS Δβ. Figure 8A shows that the majority of drought years correspond to years of (overly) depressed GPP; however, the depressed vegetation production response is strongly dependent on the treatment of β, which can be appreciated by how clustered the coloured dots are, and less dependent on the intensity of the drought, as the dots tend to be organized in horizontal bands, rather than on a diagonal. Simulations QS and C6 correspond to the data clusters with the least bias in GPP for those soil moisture drought years, independent of location. Figure 8B shows the errors in the evaporative fraction (EF), which is included here to account for differences in the sum of SH and LE between simulations and FLUXCOM products. The errors are now aligned on a diagonal, with smaller values of Δβ mostly corresponding to a negative bias (caused by a positive error in SH and a negative error in LE, see Figure 8C, and Figure 8D), and larger values of Δβ leading to a positive bias in EF. Simulation QS is closest to EF observations for 2003, with C6 second-best, albeit similar to QB, and with QM displaying the largest bias. The small amplitude Δβ (negative EF bias) seem to correspond to points in the Northern portion of the domain.
For LE flux, Figure 8C shows that the error metric is now aligned on a diagonal (unlike GPP) and tends to be dominated by the southern portion of the domain, where JULES strongly underestimates LE in the climatology. As the drought severity increases, the LE flux bias is reduced, suggesting that the LE response to Δβ starts to plateau, and the model has an opportunity to "catch up" to observations. There are indications that QM, the most liberal treatment of β in terms of soil water usage, has the least error overall, while QS, the most conservative, is exaggerating the LE response, particularly at moderate Δβ levels. Relative advantages are not so clear for the other treatments (QB, C6), even for the largest events. Figure 8D shows the error in SH as a function of Δβ, which is aligned on a diagonal with slope opposite to the one for LE: SH is strongly underestimated in severe drought cases, meaning that the observed heating from the anomalous 2003 soil dessication is not Frontiers in Environmental Science | www.frontiersin.org August 2021 | Volume 9 | Article 689301 simulated faithfully, and slightly overestimated in moderate drought cases. For this variable, simulation QM seems to be slightly worse than the others for the most pronounced droughts (better skill for low intensity droughts), while QS seems to have a small advantage overall at medium Δβ levels; C6 is quite similar to QB in this respect. The overall lesson is that, while GPP errors are affected systematically by the choice of β treatment, the errors in energy fluxes are shifted to different equilibria, reflecting climatological shifts in soil moisture, and could be interpreted as a near-neutral response.

Response to Extreme Atmospheric Forcing:
The European Heat Wave of 2003 Figure 9 shows the evolution of GPP during the 2003 summer heat wave, as monthly means, from April to August, for simulations QB, QS, and C6, the most interesting ones, analyzed in terms of the evolution of model error (against FLUXCOM). While it is well-known that Spring 2003 saw clear skies and warmer conditions, leading to enhanced vegetation activity while there was still ample soil moisture availability, both QB and C6 exaggerate this response, particularly over France and Germany, with C6 showing the largest response in April, which implies a larger soil moisture usage. This deficit can then lead, via soil memory, to a suppressed GPP in each of the JJA months; in this respect, C6 is an improvement over QB, as it is apparently able to sustain a larger JJA GPP. Simulation QS is nearly identical to C6, albeit with a slightly smaller amplitude for the JJA signal for the Western Europe mean, therefore the smallest error overall, while simulation QM (not shown) has a far smaller response.
For latent heat flux during the Summer 2003 heat wave (Figure 10), simulation C6 hardly shows any difference to CTL (better evidenced in the soil moisture maps, Figures 6, 11), while simulation QS displays a more interesting seasonal behavior, with a regional worsening of the dry bias in Spain and Eastern Europe (June), then in a region that covers France, Germany and Poland by August. This indicates that the enhanced GPP response in Spring has indeed consequences for the water cycle, via soil moisture memory, so that the errors in JJA are larger than for QB and C6. Simulation QM (not shown) shows a similar, albeit even more pronounced enhanced spring behavior, with a strong overestimation throughout, particularly from May to July, in most of central Europe, but then a relative reduction of the dry bias in August. These behaviors might amplify and lead to far worse summer conditions if the heat wave had been repeated in 2004 in this region. Figure 11 reveals a soil moisture evolution for 2003 that mostly explains the above results for GPP and latent heat. Simulation QS is increasingly accumulating soil moisture in the soil, in a rather linear way, from April to August, but this is then reflected in too strong a limitation of transpiration, thus latent heat, as will be shown in Figure 12. Simulation C6 is closest to CTL, and demonstrates a neutral response to the soil moisture stress, despite having achieved the best performance in terms of the simulation of GPP in Figure 9. Simulation QM (not shown) responds as an enhanced QS: starting from already dry conditions, it develops a seasonally increasing drying in the soil, which is compatible with too liberal a use of soil moisture, e.g., for Eastern Europe in spring and summer. Figure 12 shows a comparison of ecosystem water use efficiency (eWUE, calculated here from NPP and evapotranspiration at grid-box level) for the year 2003, one of the key years identified in Peters et al. (2018), as a seasonal evolution. It is clear that, as early as April, a large area of Europe experienced high eWUE, as the atmosphere was clear and net radiation abundant, but springtime soil moisture was already starting to become depressed, as was seen in Figure 11. As the season evolved, and assimilation (A) dropped due to increased plant water stress, particularly for the southern part of the domain, differences start to become obvious for the treatments involving a stomatal conductance route: QS and C6 show a clear advantage over the entire 45-60 latitude band, from May to August. This is quite easy to explain for simulation QS, which showed a large 2003 soil moisture surplus when compared to the other simulations, but more subtle for simulation C6, which shows virtually no change in soil moisture when compared to CTL, and even to QB. The C6 response must then be explained by changes in carbon assimilation, which must have become more efficient. It remains to be seen whether this advantage represents a realistic response, by further verification against observations.  comparing the different responses, from left to right, pathway QB sees a strong reduction in stomatal conductance in 2003, as compared to 2002 and later years, in which a vegetation recovery occurred. An even stronger response can be seen in the second column, corresponding to pathway QS, which is easy to understand from the point of view of the β (Eq. 4). The mesophyll pathway (QM) shows the least reduction of stomatal conduction in 2003, while C6 exhibits an interim response in comparison to QB and QM.

Comparison of Ecosystem Water Use Efficiency
The Diurnal and Seasonal Cycles of C i (χ) Figure 14 shows the diurnal and seasonal evolution of C i for the same grid point shown in Figure 12, which is fully indicative of χ evolution, since reference level CO 2 is prescribed in these simulations, due to the atmosphere being forced. The difference between the four different simulations is striking, particularly for years in which plants experienced strong environmental stress, such as 2003. Simulation QS exhibits the strongest response: by closing the stomata during periods of increased drying, such as early 2003, the C i level drops dramatically, as had been anticipated in Figure 2. However, at the peak of the 2003 summer heat wave, during the late stage of the diurnal cycle, when the soil moisture stress (and the mounting temperature stress, which compounds the effect) starts to reach its peak, C i returns towards the C a value (in the early morning and late evening), because photosynthesis stops, also as seen previously in Figure 2, due to a compound stress effect. Simulation C6 shows a response that is qualitatively similar, albeit with smaller magnitude, while simulations QB and QM hardly drop the C i level at all, also consistent with what had been seen in Figure 2.

Comparison With Intrinsic Water Use Efficiency From FLUXNET Observations
It is expected that in years of drought vegetation stress will also be revealed by measurements of isotopic discrminination, which was the main topic of Peters et al. (2018). The mechanisms involved in the process of photosynthesis favor the lighter 12 CO 2 molecule over the heavier 13 CO 2 molecule, at several stages in the mechanistic chain. This leads to 'discrimination': stressed vegetation will present lower δ 13 C 13 C/ 12 C values, relative to the atmosphere, and the difference is denoted by the symbol Δ δ 13 C a −δ 13 C v (units of ‰). In JULES, Δ can be estimated approximately, by using in Eq. 6, Methods. Supplementary Figure S1 shows the domain average behavior of the four C3 vegetation types in the JULES experiments, as annual anomalies for the year 2003. Because of drought stress, it is expected that NPP will be reduced compared to other years, and that Δ will be smaller, which is: see Figure S1 of Peters et al. for the curve of C i /C a vs Δ, what is shown by observations and by the model intercomparison in Peters et al. (2018).
Supplementary Figure S1 reveals that grasses suffer the most during the severe drought year of 2003, with NPP reductions of over 100 Tg C, while needleleaf trees are quite resilient (also because they are mostly located outside the area of drought, but a similar calculation with a drought mask shows the same qualitative behavior). Simulation C6 shows the least drop in NPP and the most realistic response in terms of the combined NPP/Δ metric, as compared to the data in Peters et al. (2018), with the anomaly in Δ reaching values of -0.15 for grasses, while simulation QB tends to have a Δ closer to 0 or above. It is next possible to plot the combined gridbox Δ data on a map and to compare with observations. Figure 15 shows a superposition of iWUE in observations (Peters et al., 2018), and as produced by JULES, computed by using Eq. 7, and using 3-hourly C i at PFT level. The four panels show the response (marked by dots and color-coded for change in iWUE, as originally presented in Peters et al., 2018) Peters et al., 2018), indicated that most LSMs struggle to simulate the shift to higher iWUE (greener colors), which is summarized for the observation sites in the PDFs in each inset. The original JULES-CTL had in fact failed to capture both the magnitude and spatial extent of the impact, similar to the QB simulation developed for this study, which clearly missed the high level of iWUE reached by plants in Western Europe during the 2003 heat wave, so that it was replaced with configuration C6. iWUE and eWUE originate at different scales in the model, the first stomatal, via the C i route, and the second involving area averaging of carbon and water fluxes over each grid box, weighted by PFT abundance, and are thus compared against different types of observations. The consistency of the improvement, from iWUE to eWUE, over a coarser scale and over a larger domain, for well-observed events like 2003, gives us further confidence that our proposed pathways strongly improve the JULES drought response, not just for reductions in GPP and LE (achieved in all schemes) but especially for their relative ratio (captured best by the C6 configuration).

DISCUSSION
The E11 model was implemented in a state-of-the-art LSM, JULES, now with complete freedom in enabling soil moisture stress to limit photosynthesis (QB), stomatal conductance (QS) or mesophyll conductance (QM), independent of any a-priori assumptions on what are the controlling mechanisms and pathways. A combined treatment (C6) was also tested, taken from the E11 recommendations. A conscious decision was made, at this stage, to choose a framework in which GCM biases in surface energy balance would not cloud the investigation, thus observed meteorology was used as a forcing. Additionally, the parsimonious requirements for computational resources enabled a broad set of experiments, including a lengthy LSM spin-up (cycles up to 300 years) for each, could be made feasible, unlike in a full GCM set-up. It has also to be remembered that, at this stage, the solution algorithm has yet to be made efficient, and these simulations, based on polynomial solutions, are slow, unsuitable for a GCM environment. The primary purpose of this LSM development was indeed to answer the question of whether or not the a-priori decision of implementing the β limitation on just one of the model components, e.g., photosynthesis, as in models based on the Collatz et al. (1991) formulation, has consequences for our ability to credibly predict the future of land surface dynamics. There is strong demand for predicting the fate of feedbacks between the energy and water cycles (but including carbon over longer periods), as triggered by soil moisture deficits, but such an investigation must be carried out incrementally, with many of the feedbacks purposefully disabled, in order to enable clearer understanding.
The overall model sensitivity seems to be revealed more systematically on the GPP side than on the energy fluxes, as adjustments occur, as revealed by EF, also in response to longterm soil moisture memory. This flux adjustment can be seen as a reassuring result in this particularly constrained experimental setup. Because of non-linearities in the system, models in which β is applied to stomatal conductance tend to limit transpiration, thus to conserve soil moisture, while only partially affecting GPP, except perhaps in parts of Southern Europe at the end of summer. Models in which β is applied to mesophyll conductance tend to overuse soil moisture, creating a deficit, particularly evident in semi-arid regions, and causing stress (from drought and heating) at the end of summer, but this does not seem to have any effect on GPP, implying that water use efficiency is increased in Southern Europe.
Analysis focusing on the 2003 summer heat wave revealed that the environmental stress, as an earlier and more intensive use of soil moisture by evaporation across Europe, causes an even stronger sensitivity response in the four models that have been analyzed. While a negative bias in GPP is present in all models, the 2003 heat wave is most credibly simulated by applying the combined pathway (C6) for the months June, July and August, despite an earlier overestimation of GPP. The heat wave can, however, develop more severely (making the dry bias even worse) when using models with β applied Frontiers in Environmental Science | www.frontiersin.org August 2021 | Volume 9 | Article 689301 directly to stomatal control, particularly in Southern Europe, or be classified as less severe when using any models that involve mesophyll conductance (e.g., QM, or even the more moderate C6). It is unclear whether the more intensive water use incurred by models including the mesophyll pathway would be sustainable if the anomalously dry and warm conditions were to last for more than 1 year, unlike what happened in post-2003: for Europe winter precipitation is normally sufficient to recharge soil moisture on an interannual basis, but this is not true for other regions, e.g. for the Southern United States, which has been experiencing multi-year drought conditions in recent years. For χ, the diurnal and seasonal behavior uncovered in the latter part of the study are fully consistent with the expectations raised by the idealized study in the introduction: during the heat wave, enhanced stress conditions (e.g., at noon and/or at the peak of summer) result in stomatal closure, which then results in a strong reduction of χ, but does not necessarily need imply a direct reduction of GPP, as caused instead by making the decisions, conscious or not, of imposing β on carbon assimilation, as is done in many second-generation LSMs. Comparisons provided in E11, and many references therein, indicate that low levels of χ, made possible in JULES with the new generalized analytical scheme, are credible, and indicative of a correct chain of mechanisms.
A strong confirmation of this interpretation is the fact that results of all experiments were now used to compute iWUE from C i , and showed that the C6 formulation performs comparatively better against the data collected and analysed in Peters et al. (2018) study, revealing a very realistic response of the model to the 2003 perturbation, which seems governed by mechanisms at the stomatal level, and primarily by the fast evolution of C i . The comparison of eWUE (performed on grid-scale variables) and iWUE (performed at stomatal and PFT level C i every 3 h) show complete consistency and indicates that applying the β stress to photosynthesis alone is unable to reproduce a credible vegetation response to the large 2003 event, as well as, from the new analysis in this study, to the other European droughts in the simulated period.
It has also to be remembered that C i can also be used to compute the isotopic discrimination metric, Δ, and that in the Peters et al. (2018) study JULES (in the C6 configuration) performed best when assessed against independent isotopic measurements, and against other models, some of which had been specifically developed to simulate isotopic discrimination. This further confirmation (and more analysis, see Supplementary Material), via the C i chain of mechanisms, indicates that C i dynamics is more credibly represented by allowing soil moisture stress to affect all plant function at once, which has There are, however, shortcomings in the adopted methodology: it is impossible, with the current set of experiments, driven by observed meteorology, to estimate the full impact of the surface temperature feedback, because the experimental setup only allows feedback to the leaf-level temperature, but this is then nudged towards the value of the atmospheric near-surface temperature, for each JULES model time step in the 3-hourly meteorological forcing data intervals. If anything, because of the use of observed meteorology as drivers, the present results are indicative of a muted version of possible land-atmosphere feedbacks that can be unleashed by choosing one model configuration over another.
The impact of feedbacks would very likely be significantly larger in a free-running GCM experiment. This may result in more dramatic conclusions with respect to suitability for the CMIP-type climate predictions, particularly if any form of dynamic vegetation (e.g., dynamic phenology or phenology plus competition) is taken into account, as the present results indicate that soil moisture response is over large scales and long time periods, thus with the potential for altering important components of the terrestrial carbon cycle.
However, while the experiments in this study have demonstrated that the generalized β scheme is fully functional (and numerically stable), the solution scheme is more expensive (at least a factor 3x) than the original JULES scheme, because of the use of complex polynomials (see Supplementary Material), which are onerous and hard to optimize, even with modern compilers.

A Practical Proposal for Future Solution Schemes for the E11 Model
As discussed above, the new implementation of the E11 model, now inside a state-of-the-art land surface scheme, JULES, was successful, and the prognostics that are produced seem sound in all respects. The scheme is, however, more expensive than the original (iterative) scheme in JULES; this is currently untenable for use in GCMs, but this study suggests that the same set of experiments should now be run in a full GCM setting. In order to find viable configurations, Supplementary Figure S2 in Supplementary Materials suggests possible future avenues for a more computationally parsimonious implementation. Starting from the top, the standard JULES scheme is summarized, then the E11 scheme implemented in this paper. Two further solution schemes are proposed. The first one is inspired by a classic forward-backward in time (FTBT) scheme, used in a variety of ocean and atmospheric models (e.g. in WRF), either A or g s are diagnosed at the start of each time step, and then a quadratic equation for C c is solved at each time step (depending on g s or A on alternate time steps), which is cheaper than solving a cubic equation. A second additional scheme is proposed, in which g s is initially estimated as a simple prognostic (as in Sellers et al., 1996, but also reminiscent of the modified Matsuno, 1966, scheme, in which an initial Euler step precedes an imitation of an implicit time scheme), leading again to a quadratic equation for the other two variables. This implementation was successful in SiB2 and also reduced fast time oscillations, which can potentially be incurred by FTBT-type schemes. This same approach could in the future be applied to the prediction of χ, or to C c .
Once a suitable numerical implementation is found, with computational costs comparable to the standard JULES scheme, and once a C4 scheme for the β pathway is developed, it will be possible to run the same type of experiment in the HadGEM3-GC31 GCM, albeit only for the most promising configuration, C6, given the large computational costs involved.

SUMMARY AND CONCLUSION
We have implemented the generalized, analytical and simultaneous soil water stress scheme of Egea et al. (2011) into a state-of-the-art land surface model, JULES, which can be used at multiple scales, including global offline and coupled to the HadGEM3 GCM. In this study, we have chosen to focus on the European region, and to drive JULES with observed meteorology, in order to retain a degree of control on feedbacks and to enable comparison with the standard (CTL) JULES configuration. The implementation required the development of a new solution set for the three simultaneous equations that are used to prognose A, g s , and C c . A further development in Egea et al. (2011) that was implemented in JULES is the inclusion of mesophyll conductance. The new prognostic scheme was shown to be successful and to enable complete freedom in imposing soil water stress, β, via any pathways (stomatal, biochemical, mesophyll) or any combinations thereof, including non-linear relationships between soil moisture values and β.
Results show that the treatment of soil moisture stress matters to the simulation of land surface climate in Europe, particularly for summer, but even extending into Spring and Autumn. All land surface prognostics are affected by the choice of β pathway, but particularly those controlling water and carbon fluxes between the land and the atmosphere, as well as soil moisture dynamics, which shows a cumulative effect, with strong depletion if a mesophyll conductance pathway is chosen, and strong surplus (particularly in the Southern part of the domain) if a stomatal pathway is chosen. Seasonal and climatological feedbacks between soil moisture levels and land surface fluxes are then triggered.
Responses to the drought involved in the exceptional 2003 summer heat wave, also expanded to other droughts in the last 30 years, demonstrate how important vegetation function feedbacks, from intra-seasonal to interannual, can be initiated even in a modelling framework in which the atmospheric signal is imposed. Under such conditions, it has been shown that models including a stomatal pathway can reach very low levels of stomatal conductance, accompanied by very low, and realistic, χ levels, which result in more realistic simulation of intrinsic and ecosystem Water Use Efficiency, as well as, via a specific metric of isotopic fractionation (Δ), large scale anomalies in vegetation activity. This skill in representing the chain of mechanisms involved in drought response is not normally seen in LSMs, suggesting that the realistic prediction of the feedbacks involved in changes in the carbon cycle requires re-visiting some of the most fundamental assumptions in LSMs used for climate prediction.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
PV, GE and AV designed the experiments. GE and PV implemented the E11 model into JULES, which required the development of new numerical solution methods (see Supplementary Materials), mostly accomplished by GE. PV and BB-S ran all experiments with JULES 2.2; PV then reported the E11 model to the more recent JULES 4.4 and ran initial quality assurance and stability experiments, but PM finally completed the full set of experiments presented in this article. MT and OM performed the analysis of climatological biases and variability; MT produced the analysis of drought years and the synthesis of annual GPP responses, as well as revised energy flux figures. PV conducted the idealised WUE and χ analyses at the start of the article and wrote the diagnostic code for C i , VPD and T* at stomatal and PFT levels. WP helped with the water-use efficiency analysis and provided comments on the text.