Improvement of modelling plant responses to low soil moisture in JULESvn4.9 and evaluation against flux tower measurements

. Drought is predicted to increase in the future due to climate change, bringing with it myriad impacts on ecosystems. Plants respond to drier soils by reducing stomatal conductance in order to conserve water and avoid hydraulic damage. Despite the importance of plant drought responses for the global carbon cycle and local and regional climate feedbacks, land surface models are unable to capture observed plant responses to soil moisture stress. We assessed the impact of soil moisture stress on simulated gross primary productivity (GPP) and latent energy ﬂux (LE) in the Joint UK Land Environment Simulator (JULES) vn4.9 on seasonal and annual timescales and evaluated 10 different representations of soil moisture stress in the model. For the default conﬁguration, GPP was more realistic in temperate biome sites than in the tropics or high-latitude (cold-region) sites, while LE was best simulated in temperate and high-latitude (cold) sites. Errors that were not due to soil moisture stress, possibly linked to phenology, contributed to model biases for GPP in tropical savanna and deciduous forest sites. We found that three alternative approaches to calculating soil moisture stress produced more realistic results than the default parameterization for most biomes and climates. All of these involved increasing the number of soil layers from 4 to 14 and the soil depth from 3.0 to 10.8 m. In addition, we found improvements when soil matric potential replaced volumetric water content in the stress equation (the “soil14_psi”


Introduction
Drought has a range of impacts on terrestrial ecosystems (Allen et al., 2010;Choat et al., 2012), plays a role in feedbacks on the weather and climate systems across scales (Seneviratne et al., 2013;Lemordant et al., 2016;Miralles et al., 2019;Lian et al., 2020), and affects the global carbon cycle (Green et al., 2017;Humphrey et al., 2018;Peters et al., 2018). These impacts and feedbacks have the potential to affect society, either directly through moisture availability effects on crops or indirectly by adjusting near-surface temperatures or forcing large-scale variations to the climate system. Roughly 40 % of the vegetated land surface is limited by seasonal water deficits (Nemani et al., 2003;Beer et al., 2010), which are a major control on gross primary productivity (GPP) in sub-humid, semi-arid, and arid regions (Stocker et al., 2018). In the future, soil moisture stress for ecosystems is predicted to increase over large regions (Berg et al., 2016;. In this paper, we define "soil moisture stress" as the physiological stress experienced by vegetation due to its interactions with dry soils. For these reasons, accurate process-based models of plant response to soil moisture stress are needed in coupled land-atmosphere climate models. However, the models used to represent biogeophysical and biogeochemical processes in Earth system models (ESMs) are often unable to properly capture observed responses to soil moisture stress (Beer et al., 2010;Powell et al., 2013;Medlyn et al., 2016;De Kauwe 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. Commonly, plants respond to low SMC by reducing their stomatal aperture to conserve water and protect the xylem from damage (Field and Holbrook, 1989;Sparks and Black, 1999). Embolism is caused by low soil and/or leaf water potential due to dry climatic conditions, and it causes water tension inside the plant to increase enough to drive the formation of air bubbles within the xylem vessels (Lambers et al., 2008;Choat et al., 2012). Embolized xylem is unable to transport water, and for some vegetation types, this is a dominant cause of plant mortality under drought conditions (Brodribb and Cochard, 2009;Choat et al., 2018). To avoid this, many plants limit water loss by reducing their stomatal conductance when soil moisture levels reach a certain threshold (Tyree and Sperry, 1989;Sperry et al., 1998;Choat et al., 2012) or by shedding leaves . High atmospheric vapor pressure deficits (VPD), which sometimes occur in conjunction with meteorological drought, may also result in stomatal closure. The reduced stomatal conductance triggers a cascade of other responses, beginning with reduced rates of photosynthesis (Ball et al., 1987), which reduce carbon uptake and possibly growth and change allocation between above-and below-ground stocks (Merbold et al., 2009b;Doughty et al., 2015). Lower stomatal conductance will reduce transpiration, which causes more surfaceavailable energy to be converted into sensible heat. This transference of latent to sensible heat can contribute to further desiccation of soils, increased land surface temperature, and amplification of heat waves . Over the long term, droughts can lead to changes in plant species composition (Liu et al., 2018) or large-scale forest mortality (Mcdowell et al., 2008), sometimes causing a transient situation where large ecosystems switch from being a sink of carbon dioxide to a source Gatti et al., 2015).
There is a spectrum of mechanisms through which species tolerate or acclimate to drought, meaning a "one-size-fits-all" approach to modeling can be inadequate. Explicit model representations of the xylem hydraulics are complex and difficult to parameterize globally. The emergence of plant trait databases has enabled early models to represent the hydraulic properties of the soil-plant-atmosphere continuum Eller et al., 2018Eller et al., , 2020De Kauwe et al., 2020;Sabot et al., 2020). In addition, new approaches are emerging that focus on "plant profit maximization", where photosynthetic uptake of CO 2 is optimally traded against plant hydraulic function as an alternative to the empirical functions commonly used in models to regulate gas exchange during periods of water stress (Sperry et al., 2017;Sabot et al., 2020).
For now, land surface models (LSMs) more often represent 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 (depending on the model, see Sect. 2). Above a critical SMC, there is no stress (β = 1), and below the critical threshold value, stress increases as SMC decreases until the wilting point is reached (β =0). Alternative, yet related, expressions are available whereby stomatal regulation occurs through changes in the soil matric potential (ψ, expressed in pressure units, such as MPa); θ and matric potential (a measure of how tightly the water is held in the soil pores, thereby affecting water uptake by the roots) are closely related via Figure 1. Comparison of JULES soil moisture stress factor (β) to measurements from various potted experiments from Verhoef and Egea (2014). β is calculated from Eq. (4). Two different values of p0 (Eq. 5) are shown: p0 = 0.4 was used for the "soil14_p0" and "p0" soil moisture stress experiments. the water retention curve. However, using one function for all plant responses to drying soils can result in errors; for example, the parameters describing plant and soil hydraulic responses to soil moisture may change in time  and can vary between ecosystem types . Such variation may be in response to climate change or evolving vegetation and soil properties and their structure.
In this study, we focus on the effects of droughts on vegetation that occur due to low SMC. Although droughts are often associated with changes beyond low precipitation levels, including high air temperatures and VPD, these climate drivers have their own set of impacts on vegetation, adding to the effects of low SMC, which will not be addressed here. We explore different ways in which soil moisture stress can be represented in a widely used model of the terrestrial biosphere, the Joint UK Land Environment Simulator (JULES) Clark et al., 2011). JULES is a community model and is used in coupled or stand-alone mode forced by meteorological variables. Its applications are on timescales ranging from weather forecasting to climate projections, and the model is the terrestrial component of the UK Earth System Model and the HadGEM family of models (The HadGEM2 Development Team, 2011). The spatial scales are similarly diverse. Studies range from single-point modeling of crop yield at one site (Williams et al., 2017), which requires detailed knowledge of one crop variety under carefully controlled conditions, to global predictions of land sources and sinks of CO 2 for the annually updated Global Carbon Project (Friedlingstein et al., 2019), which requires reliable performance for all vegetation types across the globe. The aim of this study is to find an improved general model equation and parameters for global applications of JULES.
Soil moisture stress has been identified as a key driver of variability in JULES projections . Verhoef and Egea (2014) showed that the standard β function in JULES and similar LSMs needs urgent attention as to whether it is the most appropriate functional form and/or if it has been parameterized correctly. For example, JULES calculates β based on θ , but using soil matric potential instead results in a curvilinear increase in stress as soils dry, which may be more realistic ( Fig. 1; Verhoef and Egea, 2014). In an evaluation of the model across 10 flux tower sites, Blyth et al. (2011) showed that the "dry-down" of the sites in semiarid areas was too quick and the seasonal variation of evaporation in the tropics was too great, possibly due to the roots being modeled as too shallow  or due to modeled stress beginning when soils were still relatively wet. Other studies have suggested that the root depths of LSMs were too shallow (Teuling et al., 2006;Wang and Dickinson, 2012). Indeed, some LSMs (CLM, SiB3, TERRA-ML) were able to improve model performance by representing deeper (e.g., 10 m) and more efficient roots (Baker et al., 2008;Akkermans et al., 2012;Liu et al., 2020).
Evaluating the impact of simulated soil moisture stress on vegetation requires that other model errors that also affect CO 2 and water fluxes are minimized. For instance, it is possible that the rapid drying found in Blyth et al. (2011) was due to overestimation of soil evaporation. The fact that land surface models in general overestimate evapotranspiration during wet periods is well documented Mueller and Seneviratne, 2014;Martínez-De La Torre et al., 2019) and leads to unrealistically low soil moisture after long dry periods (Ukkola et al., 2016). The high evaporation (and subsequent low SMC) could be due to errors in factors not being addressed in this study, such as radiation absorption or turbulent exchanges with the atmosphere. Leaf area index (LAI) also strongly affects the magnitude and seasonality of fluxes coming from vegetation and soil (via variations in shading).
This study aimed to evaluate the simulation of GPP and latent energy flux (LE) for a range of biomes and climates, to diagnose sites and seasons when soil moisture stress affects the results, and to evaluate different methods for representing soil moisture stress in JULES as a first step in improving the simulated plant responses to low SMC in global applications of JULES. To do this, we chose a subset of sites in the FLUXNET2015 database and from the Large Scale Biosphere-Atmosphere Experiment in Amazonia (LBA) experiment based on availability of data. Where possible we prescribed soil moisture and LAI from site measurements to differentiate the roles of SMC, β parameterization, or modeled phenology in model biases. We used the GPP calculated before soil moisture stress is applied to understand seasons and locations where the β parameterization was contributing to model errors. We also reviewed other commonly used approaches for modeling soil moisture stress, presented in Sect. 2.2, to motivate the representations evaluated in the remainder of the paper. This work is one of the first published results from a JULES community-wide focus group (called a JULES Process Evaluation Group, or JPEG) on understand-ing soil moisture stress impacts on vegetation, which began in 2016.

Photosynthesis and stomatal conductance in JULES
The Joint UK Land Environment Simulator (JULES) Clark et al., 2011) is a process-based model that simulates the fluxes of carbon, water, energy, and momentum between the land surface and the atmosphere. JULES treats each vegetation type as existing on a separate tile within a grid box. Energy and carbon flux calculations are performed separately for each tile, depending on plant functional type (PFT)-dependent parameters. The tiles share a common soil column. Leaf-level net photosynthesis is integrated over the canopy, according to the canopy radiation scheme specified.
In the present study, we used 10 canopy layers of equal LAI (in JULES this is "canopy radiation model 6"), although another option in JULES is to use a "big leaf" approach . Potential (non-stressed) photosynthesis is calculated based on three limiting rates: W c (a RuBisCOlimited rate), W l (a light-limited rate), and W e (a transportlimited rate for C 3 plants and a PEP carboxylase limitation for C 4 plants). For full details on the photosynthesis scheme in JULES, see Clark et al. (2011) and Harper et al. (2016). Stomatal conductance to water vapor g s (in m s −1 ) is related to net photosynthesis A (in mol CO 2 m −2 s −1 ) through where c a and c i are the atmospheric and intercellular CO 2 concentrations, respectively, in Pa and 1.6 is the molar diffusivity ratio of CO 2 to H 2 O in air (Guerrieri et al., 2019). R is the universal gas constant (8.314 J K −1 mol −1 ), and T * is the leaf temperature (K). Vapor deficit at the leaf surface (D, kg kg −1 ) affects stomatal conductance through the gradient between c a and c i : Here, * is the photorespiration compensation point (Pa) and D crit and f 0 are PFT-dependent parameters (Cox et al., 1998;.

Soil moisture stress in JULES and other terrestrial biosphere models
Many land surface, terrestrial biosphere, and crop models include a β function to represent the effect of soil moisture stress on vegetation. The implementation of the stress factor can generally be split into two categories: stomatal and biochemical limitation De Kauwe et al., 2015). JULES falls under the latter category, with potential leaf-level carbon assimilation, A p , being converted to the water-limited net leaf photosynthesis through multiplication with the stress factor: Other land surface models apply biochemical limitation through reducing RuBisCO or reducing electron transport (e.g., ORCHIDEE, Krinner et al., 2005). CABLE applies limits to both the stomata (via reducing g s ) and A (De Kauwe et al., 2015). In JULES, soil moisture stress (β, unitless) for each soil layer k is a function of volumetric water content (θ ) in each layer (θ k , m 3 m −3 ) using where θ wilt and θ upp are the water contents at the wilting point and at which the plant starts to become water stressed, respectively (Cox et al., 1998). θ upp is a function of θ crit , the critical water content (usually defined as the field capacity), and p 0 , a PFT-dependent parameter: The parameter p 0 was added to JULES in version 4.6 to allow β = 1 for θ < θ crit , in other words delaying the critical threshold value for inducing stress as soils dry below the field capacity. In the default configuration, p 0 is set to 0 (meaning θ upp = θ crit ), and θ wilt and θ crit correspond to soil matric potentials of −1.5 and −0.033 MPa, respectively. Equation (4) means that, for each soil layer, soil moisture stress completely limits root water extraction from that layer if θ k is at or below the wilting point (β k = 0), while there is no soil moisture stress (β k = 1) if θ k is at or above θ upp,k . In between these points, there is a linear increase in stress (decrease in β k ) as water content decreases (blue line in Fig. 1). An effective root fraction per layer (r k ) is used to calculate the overall soil moisture stress factor: where r k = e z/d r .
In Eq. (7), z is the depth of each soil layer, and d r is a PFT-specific parameter that weighs the effective root fraction within each layer. r k is an effective root fraction and is not the same as the actual root mass distribution, as it accounts for other traits and processes not present in JULES, such as the surface area of roots, conductivity, and hydraulic redistribution. JULES has four soil layers (n soil = 4) that together extend to 3 m depth (Fig. 2a). The smaller the d r , the more emphasis is given to shallow layers, while deeper layers are emphasized with a larger d r . As a specific example, with JULES default soil depth of 3 m, 87 % of the root water extraction is from the top 1 m for C 3 and C 4 grasses (d r = 0.5), compared to 45 % in the top 1 m for tropical broadleaf evergreen trees (d r = 3.0). As Fig. 2 shows, d r is not the root depth because roots are present in every soil layer, even though the fraction of roots is very small towards the bottom of the column for small values of d r .
The stress factor is also applied to leaf maintenance respiration (and optionally to stem and root maintenance respiration). The effective root distribution and stress factor also affect the fraction of total plant transpiration extracted from each soil layer, k : Although not used in this study, it is worth noting that many land surface and terrestrial biosphere models apply soil moisture stress through limiting stomatal conductance (the "stomatal" grouping from Bonan et al., 2014) (Egea et al., 2011;Fatichi et al., 2012;De Kauwe et al., 2015). These include JSBACH and DLEM (Raddatz et al., 2007;Tian et al., 2010). For example, CABLE uses β to modify the slope of the relationship between stomatal conductance and net photosynthesis (De Kauwe et al., 2015). In other models (e.g., crop model WOFOST), they interact through allowing the actual or potential evapotranspiration to impact the soil moisture threshold for unstressed vegetation (Tardieu and Davies, 1993). Models that limit stomatal conductance from soil moisture stress can include the explicit consideration of the plant or soil hydraulics (Williams et al., 1996;Zhou et al., 2013;Bonan et al., 2014;Mirfenderesgi et al., 2016;Eller et al., 2018;Kennedy et al., 2019;De Kauwe et al., 2020) and/or chemical signalling, such as the abscisic acid (ABA) concentration in the xylem sap (Tardieu and Davies, 1993;Dewar, 2002;Verhoef and Egea, 2014;Huntingford et al., 2015;Takahashi et al., 2018). In other models, β can affect root growth and leaf senescence (Arora and Boer, 2005;Song et al., 2013;Wang et al., 2016) or reduce mesophyll conductance (Keenan et al., 2010).

Alternative representations of soil moisture stress
In this study, we evaluated JULES GPP and LE using alternative parameterizations for β based on a review of methods found in the literature and supported by measurements. The 10 experiments are summarized in Tables 1 and 2, including settings in the default configuration. To summarize, these experiments aim to capture the impact of the following variables: 1. using deeper soils and roots ("soil14" and "soil14_dr*2" experiments, Sect. 2.3.1); 2. reducing the critical soil moisture content below which stress begins to increase ("p0" experiments, Sect. 2.3.2); 3. using soil matric potential rather than θ to calculate soil moisture stress ("psi" experiments, Sect. 2.3.3); 4. emphasizing deep roots that may have small fraction of total root biomass but can extract large amounts of soil water ("mod1" experiments, Sect. 2.3.4); 5. assuming a strong decay rate of root functioning for all PFTs ("soil14_dr0.5", Sect. 2.3.5). psi Use soil matric potential (Eq. 8) rather than volumetric water content (Eq. 4) to calculate β. Induces a curvilinear response. p0 Reduce the critical volumetric water content where stress begins. p0 in Eq. (5) is changed from 0 to 0.4 (dashed green line in Fig. 1).

mod1
Allow plants to access all soil moisture in the column. Equation (9) replaces Eq. (4), and Eq. (10) replaces Eq. (7). Double default d r (max value 3). d r is the maximum depth of roots instead of e-folding depth. Root profile in Fig. 2d soil14 Increase soil layers to 14 to 10.8 m depth, but d r remains unchanged. Root profile in Fig. 2b

Deeper soil column and roots (soil14 and soil14_dr*2)
Several studies have found that deep roots are an essential part of modeling plant drought responses (Canadell et al., 1996;Teuling et al., 2006;Baker et al., 2008;Akkermans et al., 2012;Wang and Dickinson, 2012). Canadell et al. (1996) found that the global average maximum root depth is 7 ± 1.2 m for trees and 2.6 ± 0.1 m for herbaceous plants, although maximum rooting depth is difficult to ascertain. For example, one study found that only 9 % of 475 rooting profiles extended to depths where roots were no longer present (Schenk and Jackson, 2005). We evaluated the impact of deeper soils by using a 14-layer soil, extending to 10.8 m depth. The 14-layer soil is being evaluated for use in future global configurations of JULES both offline and coupled in the UK Earth System Model. For example, it has been used for studying freeze-thaw dynamics in permafrost regions (Chadburn et al., 2015), but the impacts on surface fluxes in the middle and low latitudes have not yet been evaluated. In the "soil14" experiments, n soil increased from 4 to 14 and the thickness of each soil layer (dz soil ) was changed as in Table 1 to give a total depth of 10.8 m. This also increased the vertical resolution of layers in the top 2.8 m of soil, which is more accurate for solving the nonlinear Richards' equation (Mu et al., 2021). The parameter d r remained unchanged, resulting in the effective root profiles shown in Fig. 2b. As a result, for C 3 and C 4 grasses (d r = 0.5), 99 % of root water extraction was from the top 2.4 m, while for tropical broadleaf evergreen trees (d r = 3), 95 % of root water extraction was from the top 7.8 m (the remaining 5 % was from the bottom soil layer, which extended from 7.8 to 10.8 m). These numbers compare well to the observed maximum rooting depths (Canadell et al., 1996).
To evaluate the impact of placing more emphasis on deeper soil layers (in Eqs. 6 and 8), we doubled d r in an additional experiment ("soil14_dr*2") ( Fig. 2c). In these experiments, 99 % of root water extraction was from the top 4.8 m for C 3 and C 4 grasses (d r = 1), and for tropical evergreen trees (d r = 6), 87 % of root water extraction was from the top 7.8 m.

Delayed onset of stress (p0 and soil14_p0)
Measurements of transpiration rates show that plants do not limit transpiration until intermediate levels of soil dryness occur ( Fig. 1) (Verhoef and Egea, 2014). In JULES, having no stress until soils dry below field capacity can be represented with the parameter p 0 (Eq. 5), where a value of 0.4-0.5 for p 0 would capture the range of responses found in Verhoef and Egea (2014). In the "p0" experiments, we used p 0 = 0.4. This was done with both the 4-layer (p0) and 14layer (soil14_p0) soils.

Curvilinear response (psi and soil14_psi)
While Eq. (4) assumes a linear increase in stress as water content decreases, some models assume a curvilinear increase in stress (Sinclair, 2005;Egea et al., 2011) or an S-shaped curve (Tardieu and Davies, 1993;3276 A. B. Harper et al.: Improvement of modeling plant responses to low soil moisture in JULESvn4.9 Table 2. Default parameter settings (changed in experiments summarized in Table 1). In the JULES code, p0 is called fsmc_p0, n soil is called sm_levels, d r is called rootd, ψ open is psi_open, and ψ close is psi_close.

JULES parameter Explanation
Default setting Change in experiments fsmc_shape Switch that controls whether β decreases linearly with VWC θ or with soil matric potential ψ. . Nonlinear responses can be represented by a parameter to induce curvature (Egea et al., 2011) or through using the soil matric potential, ψ, rather than θ : Here, ψ open is the soil matric potential above which β = 1, and ψ open is the soil matric potential below which β = 0.

Remove root-weighted access to soil moisture (mod1 and soil14_mod1)
The measure of water availability for β can be a function of each layer's water content (Eq. 6), water in the wettest layer (Martens et al., 2017), or the contribution of the water in each layer can be weighted by the root density or plant and soil hydraulics Christina et al., 2017). Another approach is to use a function of water in the whole root column (θ ), rather than layer-by-layer, which is equivalent to assuming that plants can access water anywhere in the soil column if there are roots present (Baker et al., 2008;Harper et al., 2013): In this approach, root water extraction per layer is weighted by layer thickness (dz soil ) rather than by beta: In the "mod1" experiments, Eqs. (4) and (6) were replaced with Eq. (10) and Eq. (8) was replaced with Eq. (11). In addition, d r was implemented as the maximum root depth instead of the e-folding depth and was double its default value (with a maximum depth of 3 m). The effective root fraction in each soil layer was set equal to the proportional thickness of each layer, up to the maximum depth of roots (Fig. 2d).
In "soil14_mod1", d r was double its default value (Table 2) but without enforcing a maximum depth of 3 m (Fig. 2e).
With the default interpretation of d r , roots are present in every layer, but in these experiments plants could not access water at depths below the parameter d r . Therefore, this approach should benefit deep-rooted PFTs, as they could access more of the soil column than shallow-rooted grasses and shrubs.
2.3.5 Exponential decline of roots with depth (soil14_dr0.5) The effective root profile from grasses with n soil = 14 and depth of 10.8 m more closely resembles the observed rapid decay of root biomass with depth than the profiles for other PFTs (Zeng, 2001) (Fig. 2g). We evaluated the impact of using more realistic root distributions by setting d r to 0.5 for all PFTs in the "soil14_dr0.5" experiment ( Fig. 2f). Essentially, this gave more emphasis to shallow layers in calculating root water extraction and β and was an opposite approach of the "mod1" experiments, which gave more emphasis to the thickest soil layers.

Model set up and evaluation
We evaluated JULES at 40 sites covering eight general biome types from the tropics to the Arctic (Fig. 3, Table S1). Each JULES simulation was run with meteorological measurements taken at each site (i.e., point-scale runs rather than simulating the entire grid box). The meteorological and flux tower observations were obtained from the LBA Model Intercomparison Project (sites with "LBA" in the name) or FLUXNET2015 dataset (Pastorello et al., 2020). We selected sites with soil moisture measurements at the time of our original data request (26 July 2016). At each site, we extracted temperature, precipitation, wind speed, surface pressure, specific humidity, and longwave and shortwave radiation for running JULES at either half-hourly or hourly resolution, depending on the data available. We then used measured LE and calculated GPP as supplied in both datasets (for the FLUXNET2015 data, these are variables LE_F_MDS and GPP_NT_VUT_REF, respectively). Details of the data preprocessing are provided in the Supplement. We individually contacted site principal investigators (PIs) to gather details on LAI; the depth of soil moisture measurements (where available); and other details on soil texture, physical properties, and root depth. Based on the responses, this resulted in a subset of 21 sites with soil moisture measurements plus the additional information necessary for prescribing soil moisture in JULES. Of these sites, 14 also had the information necessary for prescribing LAI (these are listed in Table S1). Often the time period of LAI/SM mea-surements was shorter than the full record, and we only ran JULES for the time periods with the most data to avoid the need for gap-filling. The time periods of the simulations and soil layers for prescribing data are provided in Table S1.
The default plant parameter set was taken from Harper et al. (2016). When LAI was not prescribed, we used the JULES phenology scheme to predict LAI. This scheme predicts leaf growth and senescence based on temperature alone. Fractions of each PFT (or bare soil) present at the site were determined from the vegetation class (Tables S1, S2). We calculated soil properties from information supplied by site PIs where possible; otherwise, we used the grid box sand, silt, clay fractions of the Met Office Central Ancillary Program (CAP) high-resolution input file (Dharssi et al., 2009) to derive the Brooks and Corey (1964) parameters, along with the approximations of the parameters (via pedotransfer function) required for the soil hydraulic properties as detailed in Cosby et al. (1984) (Table S3). Each simulation began with a 50year spin-up of the soil moisture using recycled meteorology.
This evaluation focused on seasonal and annual timescales of fluxes. We started with daily measurements from the sites, then masked any modeled outputs on days when measurements were not available and calculated monthly means when > 50 % of the data was present. To evaluate the model performance, we used four metrics: normalized absolute error (NAE), variance ratio (VR), correlation coefficient (r), and root-mean-squared error (RMSE). The NAE gives an indication of the average model-data mismatch: where X obs is the observed flux, X mod is the modeled flux, and the overbar denotes an average taken over the entire simulation period. The other metrics were calculated from monthly mean fluxes. The VR is the ratio of variance in the simulations to the observations. For a perfect fit, the VR would be 1: lower values mean the model variance is too low and vice versa (Carvalhais et al., 2008). R is the Pearson's correlation coefficient and it gives an indication of model-data agreement on both a seasonal and year-to-year timescale. For the soil moisture stress experiments, we used Taylor diagrams based on monthly mean fluxes to evaluate the best fit, along with RMSE from fluxes averaged over daily and monthly periods and VR and correlation calculated from monthly fluxes.

Simulated GPP and ET
On average, JULES matched the pattern of observed seasonal cycle of GPP well for sites in non-agricultural biomes in temperate and cold climates (mean r > 0.79) (Fig. 4, Table 3). The correlation was fairly good for sites in tropical grass-lands and savannas (mean r > 0.70) and cropland (r = 0.67). However, the seasonal cycle was not well represented for sites in tropical dry forests (mean r = 0.43) or tropical evergreen forests (mean r = −0.10). In terms of model biases, the NAE was lowest (mean < 0.2) for GPP at tropical evergreen forest and temperate woody savanna sites, while NAE was highest in tropical grassland, tropical savanna, and cold grassland sites (mean > 0.50) (Fig. 5). The variance ratio (VR) indicates the amount of simulated variability in comparison to observations, a perfect simulation would have a VR of 1.0. A low VR indicates that simulated variability (either magnitude of seasonal cycle or interannual variability) was too low -this was the case for sites in cold grasslands and cropland (average VR of 0.35 and 0.21, respectively). On average, VR was between 0.55-0.92 for sites in tropical savannas, temperate non-agricultural biomes, and boreal forest. Conversely, a high VR indicates that simulated variability was higher than observed. Sites in tropical dry and evergreen forests and tropical grasslands had an average VR of 4.8, 5.5, and 4.8, respectively, due to an overestimated seasonal cycle (ie LBA-K67 in Fig. 6).
The model tended to perform best in temperate midlatitude climates. The average NAE and correlation (r) for temperate forest sites was 0.15 and 0.92, compared to 0.51 and 0.75 for the three sites in a Mediterranean climate (IT-CA1, IT-Ren, and IT-Col). Sites in temperate grasslands had an average NAE of 0.35 and were better simulated than those in cold and tropical grasslands (NAE = 0.50 and 0.99, respectively). NAE also was significantly higher for sites in tropical savannas (NAE = 0.79) compared to those in temperate savannas in the US (NAE = 0.14).
The model performance was also more related to climate than biome for LE. On average, the seasonal cycle of LE was well simulated for sites outside of the tropics (mean r per biome > 0.84) and for sites in tropical savannas (r = 0.79) ( Table 4, Fig. S1). However, in tropical dry and evergreen forests and tropical grasslands, the seasonal cycle was overestimated, as indicated by low correlations (mean r = 0.52, 0.29, 0.35, respectively) and high variance ratios (mean VR = 1.9, 3.8, 5.2, respectively). Model variance was close to observed for the tropical savanna sites (VR = 0.99). Unlike for GPP, the highest NAE occurred in temperate mixed forests (NAE = 0.55) (Fig. S2). The NAE was lowest for the cropland sites (NAE = 0.03), followed by tropical evergreen and dry forest sites (NAE = 0.13 for both).

Role of soil moisture stress in GPP errors
Based on the above analysis, on average the model performance is poorest for evergreen broadleaf sites, Mediterranean climates, cold and tropical grasslands, and tropical savannas. We compared the GPP that JULES would calculate if there was no soil moisture stress to the actual simulated GPP (Figs. 6, S3) to elucidate the role of soil moisture stress in generating model bias from 3.1. This was possible through a new diagnostic added to the model, which output GPP prior to multiplication by β. At the tropical evergreen forest sites (GF-Guy, LBA-K34, LBA-K67, LBA-K83, and LBA-BAN), simulated GPP decreased during the dry season, while the unstressed GPP and observed GPP remained high or even increased during dry seasons (Fig. S3), which indicates that the model was overestimating soil moisture stress during the dry season. At the tropical grassland and savanna sites (AU- Figure 5. Normalized absolute errors for simulated GPP at Fluxnet sites for the following 10 biomes: TrEF stands for tropical evergreen forests, TrDF stands for tropical deciduous forests, TrG stands for tropical grasslands, TrS stands for tropical savannas, TeMF stands for temperate mixed forests, TeG stands for temperate grasslands, TeS stands for temperate savannas, Cr stands for cropland, CoG stands for continental or high-altitude grasslands, and BoF stands for boreal forests. The sites that fall into each category are listed in the Supplement. Fog, CG-Tch, LBA-PDG, LBA-K77, and LBA-FNS), the modeled GPP was often too high, and the unstressed GPP was even higher. An exception was ZA-Kru, where the observed GPP was somewhere in between simulated GPP and unstressed GPP. There were mixed results for the sites with a Mediterranean climate (IT-CA1 deciduous broadleaf forest, US-Ton woody savanna, and US-Var grassland): stress was impacting the GPP but other processes were also affecting the simulation. For example, at IT-CA1 the modeled GPP  was very close to measured values when observed soil moisture and LAI were used, indicating that errors in soil hydrology and phenology were important at this site. At other semiarid sites (IT-Col deciduous broadleaf forest, US-Ton, and US-Var), the bias occurred during the peak growing season, when JULES GPP was lower than observed but unstressed GPP was closer to observations, indicating that soil moisture stress was impacting results at these sites. In the cold grassland sites, soil moisture stress sometimes resulted in underestimated GPP (e.g., RU-Che), possibly due to JULES not simulating enough unfrozen soil moisture at these sites. Conversely, at two temperate climate grasslands (AT-Neu and CH-Cha), the simulated GPP was too low even with soil moisture stress removed. Other sites where JULES showed a large improvement with the unstressed GPP were the aspen site in Canada (CA-Oas), Tharandt evergreen needleleaf forest in Germany (DE-Tha), the deciduous broadleaf forest in Belgium (BE-Vie), and the cropland site (US-Ne1). This analysis gives a list of sites that are useful for further exploring the role of soil moisture status in vegetation functioning: all sites with a Mediterranean climate or in tropical evergreen forests, as well as ZA-Kru, RU-Che, CA-Oas, DE-Tha, BE-Vie, and US-Ne1. These sites are further evaluated in Sect. 3.3. When prescribing soil moisture and LAI (see Sect. 2.4), the general trends in model performance were similar to prior simulations, although often the simulated GPP was less realistic with more prescribed data. This could be due to other errors within the soil physical parameterizations related to infiltration or soil evaporation (Van Den Hoof et al., 2013). The simulations at the tropical evergreen forest sites still did not resemble the measured GPP (as indicated by very low or negative correlations), even with prescribed LAI and soil moisture. It is possible that soil layers below those typically measured are influencing the forests soil water balance and canopy exchange processes, so more data are needed to accurately prescribe the full soil moisture profile. Only 14 sites had enough data to prescribe both soil moisture and LAI from site observations (Sect. 2.4), and often the time resolution of data was monthly, which for soil moisture could miss impact of extremely wet or dry periods. However, most often adding the LAI data resulted in an improved simulation of GPP, indicating biases resulting from the JULES phenology scheme. The improvements with incorporation of prescribed LAI were particularly large for the cropland sites and at LBA-RJA, which is a seasonally dry tropical forest.
We categorize the sites depending on the impact of soil moisture stress on their simulation of GPP with the most available prescribed data (for example, in the simulation with soil moisture and LAI prescribed at LBA-BAN, and for the simulation with soil moisture only at CN-HaM). The four categories are as follows.
1. Sites with underestimated GPP. Simulated GPP was lower than observed. However, β was often 1 and re-moving soil moisture stress had a small effect on the simulation, indicating the importance of other processes in regulating GPP at these sites. Two tropical (LBA-K34, LBA-RJA) and two temperate grasslands (AT-Neu, CH-Cha) sites fall into this category 2. Sites with overestimated GPP. Simulated GPP was higher than observed, so removing soil moisture stress increased GPP and made the simulation worse. This category includes one tundra site (CN-HaM), a Mediterranean woodland (IT-CA1), two coniferous evergreen forests in Finland and Italy (FI-Hyy and IT-Ren), an arid grassland (US-SRG), and two tropical savanna sites (CG-Tch, SD-Dem).
3. Soil-moisture-stressed sites. As in the first set of sites, there was a low bias in GPP but removing soil moisture stress improved the simulation. The "stressed" sites include three temperate mixed forests (BE-Vie, DE-Tha, and US-UMB), a Mediterranean deciduous forest (IT-Col), a boreal aspen forest (CA-Oas), a tropical evergreen forest (GF-Guy), and a cropland site (US-Ne1).

Stressed sites plus other errors.
At several sites, removing soil moisture stress made the simulation slightly better, but other missing processes also apparently affect the simulation. The difference between this category and the soil-moisture-stressed sites is the fact that there would still be a large bias even without soil moisture stress. Sites in this category include tropical forests (LBA-Ban, LBA-K83, LBA-K67), cropland (US-Ne2, US-Ne3), two savanna sites (ZA-Kru and US-SRM), and a tundra site (RU-Che).
The challenge is to determine a representation of soil moisture stress which improves the simulations at sites falling into categories 3 and 4 without degrading the simulation at the other sites. Clearly, we do not want to completely remove soil moisture stress as this plays an important role in regulating seasonal cycles in many ecosystems. In the remainder of the paper, we will focus on examples of changes at some of these sites.

New treatments of soil moisture stress
We ran the 10 experiments (Sect. 2.3, Table 1) at a subset of 11 sites that span the categories listed in Sect. 3.2. This included four sites where soil moisture stress was the main contributor to model biases (soil-moisture-stressed sites: GF-Guy, BE-Vie, DE-Tha, and CA-Oas), sites with a Mediterranean climate (IT-Col, US-Var, US-Ton), and sites with soil moisture stress plus other errors (LBA-K67, LBA-BAN, ZA-Kru, and RU-Che). Because some experiments focused on extending the soils far below the deepest soil moisture measurements available, we were unable to use prescribed data for these experiments. Taylor diagrams for GPP and LE for all sites are shown in Figs. S5 and S7, respectively, and

Soil-moisture-stressed sites
At these sites, there was an improvement when the 14-layer soil was combined with model settings p0, psi, or dr*2 (representing, respectively, setting p 0 in Eq. (3) to 0.4; using Eq. (9), which depends on the soil matric potential, to represent β; and doubling the parameter d r ). Monthly RMSE decreased from 2.30 g C m −2 d −1 on average to 1.59, 1.54, and 1.73 g C m −2 d −1 , respectively, in the soil14_p0, soil14_psi, and soil14_dr*2 experiments, averaged across the four sites.
There was also an improvement in the VR and the correlation coefficient ( Table 5). The VR reduced from 2.15 in the default simulation to nearly 1 in the soil14, soil14_p0, and soil14_mod1 experiments. For LE, the RMSE was slightly higher in these experiments (22.57, 22.49, and 20.77 W m −2 , respectively, for soil14_p0, soil14_psi, and soil14_dr*2) compared to the default experiment (19.78 W m −2 ), and the correlation coefficient was > 0.81 (Table S4). At the tropical forest site (GF-Guy), experiments with default 3 m soil depth had correlation coefficients r < 0.4 and an exaggerated seasonal cycle, as indicated by the high normalized standard deviation in the Taylor diagrams (Fig. 7). In the soil14_p0, soil14_psi, and soil14_dr*2 experiments, the correlation r was > 0.7 (compared to 0.2 in the default configuration), and the standard deviation was closer to observed. The GF-Guy site experienced the lowest amount of soil moisture stress in the soil14_p0 and soil14_psi experiments, which led to a more realistic simulation of GPP at this site (Fig. 8). Using a shallower effective root profile (setting d r to 0.5) produced the worst results, and β was very low during the dry season at the tropical forest sites in the "soil14_dr0.5" experiments (Fig. 8). In the "soil14_dr0.5" simulation, β was still weighted by root distribution, so the dry top soil layers had a relatively large impact on the  (Fig. 7). Only the default and soil14_dr0.5 simulations produced results outside the standard deviation of measured GPP (Fig. 8). Variability (denoted by standard deviation in the Taylor diagram as well as VR close to 1) was best in the soil14_p0, p0, soil14_psi, and psi simulations.

Mediterranean climate sites
At the sites with a Mediterranean climate (IT-Col, US-Var, US-Ton), soil14_psi and soil14_p0 removed the most stress, but p0 and psi with the default soil depth also produced a good fit for GPP (Figs. S5b, S6b, Table 6). However, the RMSE for LE was significantly higher in these four experiments (RMSE = 22.55,23.59,25.52,and 26.09 W m −2 for the p0, psi, soil14_p0 and soil14_psi experiments, respectively, compared to 19.67 W m −2 in the default simulation), while the correlation coefficient was high (r = 0.85-0.87 compared to 0.88 in the default) (Fig. S7b, Table S5). US-Var and US-Ton are dominated by grass and shrubs, which have an effective root depth d r of 0.5 and 1 m, respectively. At these sites, the "soil14_mod1" experiments had β < 0.5, and GPP was underestimated during the growing season (Fig. S6b). In these experiments, access to soil moisture was not weighted by effective root fractions, and d r was double its default value and was interpreted as the maximum root depth. This meant that grasses and shrubs could not access water below 1 and 2 m depth, respectively, resulting in the strong soil moisture stress seen at the US-Ton and US-Var sites.

Sites with soil moisture stress and other errors
At the sites with soil moisture stress plus other errors, there were fewer improvements; however, RMSE decreased from 2.81 g C m −2 d −1 in the default simulation to 2.08,   Table 1) on simulated seasonal cycle of GPP at two soil-moisture-stressed sites (see Sect 2.14, and 2.17 g C m −2 d −1 in the soil14_psi, soil14_p0, and soil14_dr*2 simulations, respectively (Figs. S5c, S6c, Table 7). These sites are LBA-K67, LBA-BAN, ZA-Kru, and RU-Che. The VR was best captured in the soil14_dr*2 simulations, while the correlation coefficient was highest in the default simulation and in the soil14_dr0.5 simulation. At LBA-K67 (a tropical forest site), soil14_psi and soil14_p0 had the lowest RMSE and seasonal variation in GPP, although for all experiments the correlation coefficient was negative (Fig. S5c). When d r = 0.5 m (as in "soil14_dr0.5"), there were proportionally more roots in the top soil layers, and as these dried out, there was a sharp decline in β. This is further illustrated in Fig. S9 at the LBA-K67 site, which plots β against soil moisture in the top 1 m. In comparison, with d r = 3 m (the default value) the trees were able to access water from deeper layers, so β did not decline as rapidly. At ZA-Kru, all results were within the range of the measurements, although the growing season GPP was underestimated ( Fig. S6c). At LBA-BAN, soil14_dr*2, soil14_psi, and soil14_p0 gave the lowest RMSE, but VR was very high (> 3), and the correlation coefficient was low (r < 0.4) for all simulations. There was very little difference between any of the simulations at RU-Che, and β was < 0.25 year-round for all experiments. For LE, there was a significant reduction in RMSE from 22.54 W m −2 to < 18 W m −2 for all experiments with 14-layer soil at these sites (Table S6). The correlation coefficient was also significantly improved in these experiments (from 0.48 in the default simulation to > 0.67). The exception to these improvements was the "soil14_dr0.5" experiment, where the RMSE increased to 25.17 W m −2 and correlation coefficient decreased to 0.35.

Average response across sites
Averaging across the 11 sites where we performed additional experiments, the lowest RMSE for GPP occurred in the soil14_p0, soil14_psi, and soil14_dr*2 experiments (on both daily and monthly timescales). The variability was best captured by the soil14, soil14_p0, and soil14_psi experiments (as denoted by VR of 1.06, 1.06, and 0.98, respectively). The mean correlation coefficient was similar across all experiments (0.50-0.57). All of the experiments were an improvement compared to the default configuration, except for the p0, mod1, and soil14_dr0.5 experiments. For LE, averaged across all sites, the daily and monthly RMSE was lowest for the soil14 experiment, and this was the only experiment with RMSE lower than the default configuration. There was an improvement in the VR for the soil14, soil14_p0, soil14_psi, soil14_mod1, and soil14_dr*2 experiments compared to the default (with VR between 1.26-1.44 compared to 1.58 in the default). The correlation was highest (r ∼ 0.74-0.76 compared to default r = 0.70) for all experiments with 14-layer soil except for soil14_dr0.5.

Default model configuration
Tables 3-4 summarize some of the key findings from this study pertaining to the default configuration. JULESsimulated GPP was more realistic in temperate biome sites than in the tropics or high-latitude (cold-region) sites, as indicated by three statistics to measure annual biases (NAE), seasonal cycles (r), and variability (VR). LE was best simulated in temperate and high-latitude (cold) sites based on the same statistics (except for temperate mixed forests). For sites in the tropics, the default β parameterization contributed to an exaggerated seasonal cycle of GPP compared to the measurements, especially in tropical evergreen forests. Although the NAE was low in tropical evergreen forest sites (e.g., LBA sites K34, K83, K67, and BAN), the seasonal cycle was overestimated (despite LAI being nearly constant all year), as indicated by high VR and low correlation coefficients. A similar result was observed with LE in most tropical sites: the seasonal cycle was incorrect and the VR was high. For example, at LBA-K67, the measurements show an increasing trend in GPP from August to October (coinciding with the dry season), while JULES predicted a decreasing trend during this time. Even with soil moisture and LAI prescribed for the four tropical evergreen forest sites, the correlation coefficients were negative. At these sites, it is possible that including a seasonally varying photosynthetic capacity would improve the results, as in Wu et al. (2017). The dry season is often accompanied by enhanced carbon uptake in Amazon forests, due to a combination of fewer clouds and increased incoming solar radiation (Saleska et al., 2003;Restrepo-Coupe et al., 2013;Zeri et al., 2014) and seasonal leaf flushing (Wu et al., 2016). The observed seasonality in GPP is enabled by deep roots that can access ample soil moisture and by the relatively high photosynthetic capacity of new leaves (Wu et al., 2017), a process not yet represented in JULES.
Other errors, possibly linked to phenology, also contributed to model biases in tropical savanna and deciduous forest sites. The improvements seen when LAI was prescribed at LBA-RJA (a seasonally dry tropical forest site) further suggest that JULES' lack of a moisture-driven phenology scheme could be affecting the results at this site. LBA-RJA serves as interesting comparison to LBA-K67: RJA receives a similar amount of annual rainfall, but the dry season is more intense, with about half as much rainfall during the dry season compared to K67 . The bedrock is relatively shallow at RJA (2-3 m) ; therefore, deep soil moisture is not present. At this site, measured GPP drops steadily from January until reaching a minimum in the middle of the dry season. JULES captured this seasonal cycle very well, although the amplitude was slightly dampened with predicted GPP being higher than observed during most of the year (with prescribed LAI and soil moisture).
In cold grassland sites, JULES underpredicted the variability of GPP and had high annual biases. The biases were due to very little GPP being simulated, with β being low yearround. At RU-Che, giving more emphasis to deeper layers (with "soil14_dr*2") did not increase GPP, which is not unexpected due to the presence of frozen soils both in the simu-lations and in reality at this site (Merbold et al., 2009a). The C 3 grass PFT at this site has the most roots in the top 0.5 m, which indicates that evaporation or sublimation could be drying the soils too much in the layers with the most roots and unfrozen soil moisture content.

Overview of alternative approaches for
representing soil moisture stress We found that three alternative approaches to calculating soil moisture stress produced more realistic results than the default parameterization for most biomes and climates: 14-layer soil with a curvilinear stress response function ("soil14_psi", Eq. 9), 14-layer soil with delayed induction of stress ("soil14_p0", Eq. 3), and 14-layer soil with deeper roots ("soil14_dr*2"). Within the default configuration, LE biases were greatest in temperate mixed forests, with overestimation occurring during Spring-Autumn. At these sites, reducing soil moisture stress (i.e., with soil14_psi, soil14_p0, and soil14_dr*2) increased LE and increased RMSE, but improved the simulated seasonal cycle and variance. Further evaluation into the reason for the high bias in LE at many of the sites would enable improvements in both carbon and energy fluxes with new parameterizations for β. There is ample justification for having deeper soils and roots in JULES. Total soil column depth and root distribution determine the total amount of water and nutrients available to plants. Deep roots can access soil moisture at depth (Christina et al., 2017) and potentially the water table and hence contribute to tree transpiration during dry periods, e.g., for GF-Guy where many canopy trees are not impacted by dry season droughts (Stahl et al., 2013a, b). Deep roots have been found to be important for many vegetation types and ecosystems (Canadell et al., 1996;Pierret et al., 2016;Germon et al., 2020): for multiple tree species in tropical forests (Nepstad et al., 1994;Jipp et al., 1998;Strey et al., 2017;Brum et al., 2019), for Acacias in semi-arid savannas such as SD-Dem (Ardö et al., 2008), and for fast-growing Eucalypt and Acacia mangium plantations in Brazil (Christina et al., 2011;Laclau et al., 2013;Germon et al., 2018), to name a few examples. In particular, in tropical forests, the global average maximum rooting depth is approximately 7 m (Canadell et al., 1996). Although estimates of maximum rooting depth are uncertain (Schenk and Jackson, 2005;Pierret et al., 2016), these examples contrast with the shallow soils (3 m) in the default JULES simulations. In addition, weighting root water uptake or soil moisture stress by fraction of roots in each layer could produce too much stress if the shallow layers (with the most roots) dry out too quickly. Deep roots are very efficient at moving water, for example, specific hydraulic conductivities (K s ) of deep roots can be as much as 15 times higher than K s of superficial roots for Banksia sp. (Pate et al., 1995), and deep roots can redistribute water from deep to shallow layers (Caldwell et al., 1998;Burgess et al., 2001;Oliveira et al., 2005). However, not all plants rely on deep roots during a drought (Prechsl et al., 2015;Brinkmann et al., 2019), and at sites dominated by grasses and shrubs there were high biases in the "soil14_mod1" experiments (weighting the contribution of each layer's β i by the thickness of that layer rather than by the effective root fraction in that layer). Studies with other land surface models have drawn similar conclusions. Increasing the soil column from 3.5 to 10 m and allowing roots to access this entire reservoir improved the fit of the SiB3 model to observations at the LBA-K83 site (Baker et al., 2008). Similarly, the ability of the G'DAY model to accurately simulate wood production in fast-growing sub-tropical plantations was considerably improved by accounting for tree ability to uptake water in deep soil layers (Marsden et al., 2013). On the other hand, using the default calculation for β with an e-folding depth d r = 0.5 m emphasized shallow layers, and the overall soil moisture stress increased at most sites, resulting in a poor fit to measured GPP and LE in the "soil14_dr0.5" experiments.

Outlook for modeling soil moisture stress in JULES and other land surface models
In this study, we used flux tower observations and detailed site information when possible. Working with site researchers enabled us to narrow down reasons for model biases by prescribing soil moisture and LAI at some sites and to better understand mechanisms of drought responses at others. These are invaluable benefits of working with site-level data. Future studies could benefit from incorporating more sites (the full FLUXNET2015 dataset includes 212 sites), particularly if the focus is reducing biome-scale model biases. There is potential to extract even more information from available datasets to improve the representation of soil moisture-vegetation interactions . This includes better utilization of satellite data, and one particular opportunity is to consider soil moisture measurements in parallel with those of solar-induced fluorescence, which is used to estimate photosynthesis (Lee et al., 2013). Satellite records have large spatial coverage, and modern machine learning algorithms could be used to characterize Earth observation datasets of drought conditions (Huntingford et al., 2019). Such methods could address the difficulty in modeling the high complexity and geographical diversity of plant adaptive responses to soil moisture deficits that exist in nature. Future work should build upon these results to further evaluate JULES response with these parameterizations, focusing on deeper soils and either using a non-zero p 0 (we used 0.4 in this study) or using the soil matric potential (ψ) rather than volumetric water content for calculating β. We note that such alternative parameterizations are not a replacement for improved representations of the soil-plant hydraulic system that have been developed for many models Kennedy et al., 2019), including JULES (Eller et al., 2020). Instead, they provide a practi-cal, alternative way to represent some aspects of the soilplant hydraulic system, including hydraulic differences between PFTs through the parameters ψ open and ψ close (Eq. 9), which can be adopted by any model that uses the β function to represent vegetation responses to soil moisture. Several other land surface models use soil water potential (e.g., CLM;Lawrence et al., 2019) for calculating soil moisture stress, and a further benefit of this approach is the ability to set PFT-specific values for ψ open and ψ close (Eq. 9), with measured turgor loss points serving as a starting point for ψ close (Bartlett et al., 2012). Whereas our new parametrization generally improves JULES skill to simulate GPP and LE, it remains to be tested if similar results would be achieved by other models, including models that apply the β function at different parts of their photosynthesis and stomatal conductance schemes (e.g., Keenan et al., 2010;De Kauwe et al., 2015).
Currently, the land partially offsets anthropogenic CO 2 emissions by photosynthetic drawdown, but this could be reversed if droughts increase in frequency or intensity in the future. Feedbacks from the land surface can amplify and lock-in existing drought conditions (Morillas et al., 2017), and land surface responses to regional drought can affect precipitation and circulation in other regions Lian et al., 2020). Improving responses of vegetation to drought in land surface models such as JULES would have far-reaching implications for global climate modeling and are therefore of utmost importance.
Code availability. Both the model code and the files for running it are available from the Met Office Science Repository Service: https://code.metoffice.gov.uk/ (last access: 13 April 2021). Registration is required and code is freely available subject to completion of a software license. The results presented in this paper were obtained from running JULES branch https://code.metoffice.gov.uk/trac/jules/browser/main/branches/ dev/karinawilliams/r9227_add_gpp_unstressed_diagnostic (last access: 13 April 2021,  which is a branch of JULESv4.9 with the additional unstressed GPP diagnostic added. The runs were completed with the Rose suite https://code. metoffice.gov.uk/trac/roses-u/browser/a/l/7/5/2/u-al752-jpegpaper (last access: 13 April 2021, , which also includes Python scripts for creating the plots. The Taylor diagrams (Figs. 7, S5 and S7) were made with Python scripts from Yannick Copin (https://gist.github.com/ycopin/3342888, last access: 13 April 2021, Copin, 2018).
Author contributions. This study is the result of a large effort within the JULES community to better understand soil moisture stress and simulated responses to soil moisture deficits. All coauthors contributed at some point to writing or improving the manuscript. Flux tower researchers provided particular insight into their sites: LeM (IT-Ren); IM (FI-Hyy); DB (GF-Guy); AG and YN (CG-Tch); GW (AT-Neu); NB, LuM, and KF (CH-Cha); and LuM (RU-Che, Za-Kru).
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. Flux tower measurements used in this study are from FLUXNET2015 and the LBA project. Several site PIs contributed data for the soil moisture and LAI prescribed runs, we gratefully acknowledge their contribution here: Simone Sabbatini (  The GF-Guy site is supported by an Investissement d'Avenir grant from the Agence Nationale de la Recherche (CEBA: ANR-10-LABX-0025; ARBRE: ANR-11-LABX-0002-01). CA-Oas is part of the Fluxnet Canada network, supported by the Natural Science and Engineering Research Council of Canada (NSERC) and the Canadian Foundation for Climate and Atmospheric Science (CF-CAS).
Review statement. This paper was edited by Jatin Kala and reviewed by two anonymous referees.