Improvement of modeling 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 flux (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 configuration, 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” experiments), when the critical threshold value for inducing soil moisture stress was reduced (“soil14_p0”), and when plants were able to access soil moisture in deeper soil layers (“soil14_dr*2”). For LE, the biases were highest in the default configuration in temperate mixed forests, with overestimation occurring during most of the year. At these sites, reducing soil moisture stress (with the new parameterizations mentioned above) increased LE and increased model biases but improved the simulated seasonal cycle and brought the monthly variance closer to the measured variance of LE. Further evaluation of 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 soil moisture stress. Increasing the soil depth and plant access to deep soil moisture improved many aspects of the simulations, and we recommend these settings in future work using JULES or as a general way to improve land surface carbon and water fluxes in other models. In addition, using soil matric potential presents the opportunity to include plant functional type-specific parameters to further improve modeled fluxes.

1 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; Ukkola et al., 2020). 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; Restrepo-Coupe et al., 2017; 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 (Wolfe et al., 2016). 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 surface-available 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 (Seneviratne et al., 2010). 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 (Ciais et al., 2005; 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 (Sperry et al., 2016; Eller et al., 2018, 2020; De 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₂ 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³ m⁻³). 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 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 (Robinson et al., 2019) and can vary between ecosystem types (Teuling et al., 2010). 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) (Best et al., 2011; 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₂ 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.

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.

Soil moisture stress has been identified as a key driver of variability in JULES projections (Blyth et al., 2011). 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 semi-arid 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 (Blyth et al., 2011) 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₂ 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 (Blyth et al., 2011; 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 understanding soil moisture stress impacts on vegetation, which began in 2016.

2 Methods

2.1 Photosynthesis and stomatal conductance in JULES

The Joint UK Land Environment Simulator (JULES) (Best et al., 2011; 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 (Clark et al., 2011). Potential (non-stressed) photosynthesis is calculated based on three limiting rates: W_c (a RuBisCO-limited rate), W_l (a light-limited rate), and W_e (a transport-limited rate for C₃ plants and a PEP carboxylase limitation for C₄ 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⁻¹) is related to net photosynthesis A (in mol CO₂ m⁻² s⁻¹) through

$$\begin{array}{}\text{(1)}& g_{\text{s}}=-\mathrm{1.6}A{\displaystyle \frac{R{T}^{*}}{c_{\text{i}}-c_{\text{a}}}},\end{array}$$

where c_a and c_i are the atmospheric and intercellular CO₂ concentrations, respectively, in Pa and 1.6 is the molar diffusivity ratio of CO₂ to H₂O in air (Guerrieri et al., 2019). R is the universal gas constant (8.314 J K⁻¹ mol⁻¹), and T^* is the leaf temperature (K). Vapor deficit at the leaf surface (D, kg kg⁻¹) affects stomatal conductance through the gradient between c_a and c_i:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{C_{\text{i}}-{\mathrm{\Gamma }}^{\ast }}{C_{\text{a}}-{\mathrm{\Gamma }}^{\ast }}}=f_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{D}{D_{\text{crit}}}}\right).\end{array}$$

Here, Γ^* is the photorespiration compensation point (Pa) and D_crit and f₀ are PFT-dependent parameters (Cox et al., 1998; Best et al., 2011).

2.2 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 (Bonan et al., 2014; 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:

$$\begin{array}{}\text{(3)}& A=A_{\text{p}}\mathit{\beta }.\end{array}$$

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³ m⁻³) using

$$\begin{array}{}\text{(4)}& {\mathit{\beta }}_k=\left\{\begin{array}{ll}\mathrm{1}& {\mathit{\theta }}_k\ge {\mathit{\theta }}_{\text{upp},k}\\ \frac{{\mathit{\theta }}_k-{\mathit{\theta }}_{\text{wilt},k}}{{\mathit{\theta }}_{\text{upp},k}-{\mathit{\theta }}_{\text{wilt},k}}& {\mathit{\theta }}_{\text{wilt},k}\le {\mathit{\theta }}_k\le {\mathit{\theta }}_{\text{upp},k}\\ \mathrm{0}& {\mathit{\theta }}_k\le {\mathit{\theta }}_{\text{wilt},k}\end{array}\right.,\end{array}$$

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₀, a PFT-dependent parameter:

$$\begin{array}{}\text{(5)}& {\mathit{\theta }}_{\text{upp}}={\mathit{\theta }}_{\text{wilt}}+({\mathit{\theta }}_{\text{crit}}-{\mathit{\theta }}_{\text{wilt}})(\mathrm{1}-p_{\mathrm{0}}).\end{array}$$

The parameter p₀ 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₀ 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:

$$\begin{array}{}\text{(6)}& \mathit{\beta }=\sum _k^{n_{\text{soil}}}{r}_k{\mathit{\beta }}_k,\end{array}$$

where

$$\begin{array}{}\text{(7)}& r_k=e^{z/d_{\text{r}}}.\end{array}$$

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₃ and C₄ 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.

Figure 2 Effective root water extraction profiles for JULES. The parameter d_r is the e-folding depth for weighing root water extraction and soil moisture stress. The plant functional types (PFTs) are C₃ and C₄ grasses, evergreen and deciduous shrubs (ESh, DSh), needleleaf evergreen trees (NET), temperate broadleaf evergreen trees (BET-Te), broadleaf deciduous trees (BDT), needleleaf deciduous trees (NDT), and tropical broadleaf evergreen trees (BET-Tr). The horizontal dotted lines show the bottom depth of each layer. Profiles for the experiments in this study are shown: with the default 3 m deep, 4-layer soil (a); with an updated 10.8 m deep, 14-layer soil (b); with the 10.8 m deep soil and doubled d_r (c); with water extraction weighted by layer thickness and 3 m deep soils (d); with water extraction weighted by layer thickness and 10.8 m deep soils (e); and with all PFTs having d_r= 0.5 and the 10.8 m deep soil (f). For comparison, panel (g) shows root fractions from Zeng (2001), where distributions were calculated based on available measurements of root profiles.

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:

$$\begin{array}{}\text{(8)}& {\mathit{ϵ}}_k={\displaystyle \frac{r_k{\mathit{\beta }}_k}{\mathit{\beta }}}.\end{array}$$

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).

Table 1 Summary of the 10 JULES model experiments related to the treatment of soil moisture stress.

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.

2.3 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).

2.3.1 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₃ and C₄ 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₃ and C₄ grasses (d_r=1), and for tropical evergreen trees (d_r=6), 87 % of root water extraction was from the top 7.8 m.

2.3.2 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₀ (Eq. 5), where a value of 0.4–0.5 for p₀ would capture the range of responses found in Verhoef and Egea (2014). In the “p0” experiments, we used p₀=0.4. This was done with both the 4-layer (p0) and 14-layer (soil14_p0) soils.

2.3.3 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; Oleson et al., 2010; Egea et al., 2011) or an S-shaped curve (Tardieu and Davies, 1993; De Kauwe et al., 2015). Nonlinear responses can be represented by a parameter to induce curvature (Egea et al., 2011) or through using the soil matric potential, ψ, rather than θ:

$$\begin{array}{}\text{(9)}& {\mathit{\beta }}_{\mathit{\psi },k}={\displaystyle \frac{{\mathit{\psi }}_k-{\mathit{\psi }}_{\text{close}}}{{\mathit{\psi }}_{\text{open}}-{\mathit{\psi }}_{\text{close}}}},\end{array}$$

Here, ψ_open is the soil matric potential above which β=1, and ψ_open is the soil matric potential below which β=0. We set ψ_open and ψ_close to −0.033 and −1.5 MPa, respectively, which are typical values for field capacity and wilting point. Models that use soil water potential include (Verhoef and Egea, 2014; Fatichi et al., 2012; Manzoni et al., 2013; Lawrence et al., 2019), while other models use leaf water potential (Tuzet et al., 2003; Christina et al., 2017). In the “psi” experiments, we replaced Eq. (4) with Eq. (9). This was done with both the 4-layer (psi) and 14-layer (soil14_psi) soils.

Figure 3 Location of sites used in this study. Details on site characteristics are provided in the Supplement.

2.3.4 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 (Oleson et al., 2010; Christina et al., 2017). Another approach is to use a function of water in the whole root column ($\stackrel{\mathrm{‾}}{\mathit{\theta }}$), 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):

$$\begin{array}{}\text{(10)}& {\mathit{\beta }}_{\text{mod1}}={\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\theta }}-{\mathit{\theta }}_{\text{wilt}}}{{\mathit{\theta }}_{\text{upp}}-{\mathit{\theta }}_{\text{wilt}}}}.\end{array}$$

In this approach, root water extraction per layer is weighted by layer thickness (dz_soil) rather than by beta:

$$\begin{array}{}\text{(11)}& {\mathit{ϵ}}_k^{\mathrm{0}}=d{z}_{\text{soil}}({\mathit{\theta }}_k-{\mathit{\theta }}_{\text{wilt}}).\end{array}$$

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.

2.4 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 pre-processing 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 measurements 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.

Table 3 Summary of model performance for GPP with no prescribed data. The statistics are averages for each biome: Pearson's correlation coefficient (r), normalized absolute annual error (NAE), and variance ratio (VR).

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 50-year 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:

$$\begin{array}{}\text{(12)}& \text{NAE}={\displaystyle \frac{\stackrel{\mathrm{‾}}{X_{\text{obs}}}-\stackrel{\mathrm{‾}}{X_{\text{mod}}}}{\stackrel{\mathrm{‾}}{X_{\text{obs}}}}},\end{array}$$

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.

3 Results

3.1 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 grasslands 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).

Figure 4 Correlation coefficient for simulated monthly mean 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.

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).

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.

Table 4 Summary of model performance for LE with no prescribed data. The statistics are averages for each biome: Pearson's correlation coefficient (r), normalized absolute annual error (NAE), and variance ratio (VR).

Figure 6 Average seasonal cycle of GPP (g C m⁻² d⁻¹) for representative sites in biomes with large biases. Full dates of simulations are provided in the Supplement; here we give the years that are included in each simulation: AU-Fog (2006–2008), BE-Vie (1996–2006), CA-Oas (1996–2010), DE-Tha (1996–2014), IT-Col (1996–2014), LBA-BAN (2004–2006), LBA-K67 (2002–2003), RU-Che (2002–2005), US-Ne1 (2001–2012), US-Ton (2001–2014), ZA-Kru (2000–2013).

3.2 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-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 semi-arid 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 removing 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).

4. 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.

3.3 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 seasonal cycles of GPP, simulated β, and LE are shown in Figs. S6 and S8.

Figure 7 Example of impacts of soil moisture stress representations on GPP model skill for two soil-moisture-stressed sites: GF-Guy (tropical evergreen forest) and DE-Tha (temperate evergreen needleleaf forest). The GF-Guy simulations included the years 2007–2009, and the DE-Tha simulations included the years 1996–2014. Details of the simulations are provided in Sect. 2.3 and Tables 1–2. The solid lines indicate the centered rms error (based on the normalized standard deviations).

3.3.1 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₀ 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⁻² d⁻¹ on average to 1.59, 1.54, and 1.73 g C m⁻² d⁻¹, 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⁻², respectively, for soil14_p0, soil14_psi, and soil14_dr*2) compared to the default experiment (19.78 W m⁻²), and the correlation coefficient was > 0.81 (Table S4).

Table 5 Average results of soil moisture stress experiments for GPP at the soil-moisture-stressed sites GF-Guy, BE-Vie, DE-Tha, and CA-Oas.

Figure 8 Example of impacts of changes to the representation of soil moisture stress (see Table 1) on simulated seasonal cycle of GPP at two soil-moisture-stressed sites (see Sect. 3.3; similar examples for BE-Vie and CA-Oas are given in Fig. S9a). The GF-Guy simulations included the years 2007–2009, and the DE-Tha simulations included the years 1996–2014. GF-Guy is a tropical evergreen broadleaf forest and DE-Tha is an evergreen needleleaf forest. Details of the simulations are provided in Sect. 2.3 and Tables 1–2.

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 stress experienced by the plants. Another site in the soil-moisture-stressed category was DE-Tha, where most simulations yielded reasonable results (r> 0.9) (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.

3.3.2 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⁻² for the p0, psi, soil14_p0 and soil14_psi experiments, respectively, compared to 19.67 W m⁻² 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.

Table 6 Average results of soil moisture stress experiments for GPP at the sites with Mediterranean climate (IT-Col, US-Var, and US-Ton).

Table 7 Average results of soil moisture stress experiments for GPP at sites with soil moisture stress plus other errors (LBA-K67, LBA-BAN, RU-Che, ZA-Kru).

3.3.3 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⁻² d⁻¹ in the default simulation to 2.08, 2.14, and 2.17 g C m⁻² d⁻¹ 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⁻² to < 18 W m⁻² 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⁻² and correlation coefficient decreased to 0.35.

3.3.4 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.

4 Discussion and conclusions

4.1 Default model configuration

Tables 3–4 summarize some of the key findings from this study pertaining to the default configuration. JULES-simulated 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; Von Randow 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 (Restrepo-Coupe et al., 2013). The bedrock is relatively shallow at RJA (2–3 m) (Christoffersen et al., 2014); 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 year-round. 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 simulations and in reality at this site (Merbold et al., 2009a). The C₃ 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.

4.2 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.

4.3 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 (Gentine et al., 2019). 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₀ (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 (Bonan et al., 2014; Christoffersen et al., 2016; Kennedy et al., 2019), including JULES (Eller et al., 2020). Instead, they provide a practical, alternative way to represent some aspects of the soil–plant 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; Oleson et al., 2010; 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₂ 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 (Harper et al., 2013; 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.