Improved representation of plant physiology in the JULES-vn5.6 land surface model: photosynthesis, stomatal conductance and thermal acclimation

Carbon and water cycle dynamics of vegetation are controlled primarily by photosynthesis and stomatal conductance (g_s). Our goal is to improve the representation of these key physiological processes within the JULES land surface model, with a particular focus on refining the temperature sensitivity of photosynthesis, impacting modelled carbon, energy and water fluxes. We test (1) an implementation of the Farquhar et al. (1980) photosynthesis scheme and associated plant functional type-dependent photosynthetic temperature response functions, (2) the optimality-based g_s scheme from Medlyn et al. (2011) and (3) the Kattge and Knorr (2007) photosynthetic capacity thermal acclimation scheme. New parameters for each model configuration are adopted from recent large observational datasets that synthesise global experimental data. These developments to JULES incorporate current physiological understanding of vegetation behaviour into the model and enable users to derive direct links between model parameters and ongoing measurement campaigns that refine such parameter values. Replacement of the original Collatz et al. (1991) C₃ photosynthesis model with the Farquhar scheme results in large changes in GPP for the current day, with ∼ 10 % reduction in seasonal (June–August, JJA, and December–February, DJF) mean GPP in tropical forests and ∼ 20 % increase in the northern high-latitude forests in JJA. The optimality-based g_s model decreases the latent heat flux for the present day (∼ 10 %, with an associated increase in sensible heat flux) across regions dominated by needleleaf evergreen forest in the Northern Hemisphere summer. Thermal acclimation of photosynthesis coupled with the Medlyn g_s scheme reduced tropical forest GPP by up to 5 % and increased GPP in the high-northern-latitude forests by between 2 % and 5 %. Evaluation of simulated carbon and water fluxes by each model configuration against global data products shows this latter configuration generates improvements in these key areas. Thermal acclimation of photosynthesis coupled with the Medlyn g_s scheme improved modelled carbon fluxes in tropical and high-northern-latitude forests in JJA and improved the simulation of evapotranspiration across much of the Northern Hemisphere in JJA. Having established good model performance for the contemporary period, we force this new version of JULES offline with a future climate scenario corresponding to rising atmospheric greenhouse gases (Shared Socioeconomic Pathway (SSP5), Representative Concentration Pathway 8.5 (RCP8.5)). In particular, these calculations allow for understanding of the effects of long-term warming. We find that the impact of thermal acclimation coupled with the optimality-based g_s model on simulated fluxes increases latent heat flux (+50 %) by the year 2050 compared to the JULES model configuration without acclimation. This new JULES configuration also projects increased GPP across tropical (+10 %) and northern-latitude regions (+30 %) by 2050. We conclude that thermal acclimation of photosynthesis with the Farquhar photosynthesis scheme and the new optimality-based g_s scheme together improve the simulation of carbon and water fluxes for the current day and have a large impact on modelled future carbon cycle dynamics in a warming world.

1 Introduction

Photosynthesis and stomatal conductance (g_s) together exert a strong control over the exchange of carbon, water and energy between the land surface and the atmosphere. The behaviour of stomatal pores on the leaf surface link these processes, controlling the amount of carbon dioxide (CO₂) entering and water leaving each leaf. Photosynthesis represents the largest exchange of carbon between the land and atmosphere (Friedlingstein et al., 2020), being more substantial than respiration loss. This imbalance is central to the global carbon cycle because it slows the rate of accumulation of CO₂ in the atmosphere caused by fossil fuel burning and therefore also lowers the rate of atmospheric temperature increase. As stomata open to take up CO₂ for photosynthesis, plants also lose water through transpiration, and this flux has been estimated to account for 60 %–80 % of evapotranspiration (ET) across the land surface (Jasechko et al., 2013; Schlesinger and Jasechko, 2014). Hence, for vegetated surfaces, transpiration is the primary driver of the latent heat flux (LE), the latter describing the overall transfer of water vapour to the atmosphere. The partitioning of available net radiation between LE and sensible heat (H) is also a key determinant of land surface temperature, therefore having a feedback on photosynthesis and other key metabolic processes that influence the global carbon cycle such as plant respiration.

Land surface models (LSMs) simulate the exchange of carbon, water and energy between the land surface and the atmosphere, providing the lower boundary conditions for the atmospheric component of Earth system models (ESMs) when run in a coupled configuration. ESM projections form the main tool to predict future climate change and underpin much of the regular United Nations Intergovernmental Panel on Climate Change (IPCC) reports that inform policymakers. However, ESM predictions of the global carbon sink are fraught with large uncertainties surrounding projections of future carbon uptake (Friedlingstein et al., 2014), causing uncertainty in any translation from CO₂ emissions to atmospheric CO₂ trajectory. A lack of knowledge in how the global carbon cycle operates creates uncertainties in translating from emissions to global warming, and these uncertainties are a sizeable fraction of those associated with unknowns of physical climate processes (Huntingford et al., 2009). Therefore, given the critical role of both photosynthesis and g_s in determining land–atmosphere exchanges, their accurate representation and parameterisation in LSMs are of paramount importance. Booth et al. (2012) show that a significant uncertainty is the temperature sensitivity of photosynthesis and suggest that thermal acclimation of photosynthesis – where plants adjust their optimum temperature for photosynthesis to growth conditions experienced over the timescale of days to weeks – might reduce the spread in modelled carbon exchange. Yet despite strong evidence of the thermal acclimation capability of plant photosynthesis (Dusenge et al., 2020; Slot et al., 2021; Way et al., 2017; Way and Yamori, 2014; Yamaguchi et al., 2016), incorporation of this process in large-scale LSMs is limited to only a few, for example, TEM (Chen and Zhuang, 2013), CLM4.5 (Lombardozzi et al., 2015), LM3 (Smith et al., 2016), JULES (Mercado et al., 2018), ORCHIDEE (Krinner et al., 2005) and BETHY (Ziehn et al., 2011), and is not yet commonly represented in ESMs. Currently, the majority of LSMs and ESMs use simple fixed (i.e. non-acclimating) temperature response functions for photosynthetic capacity parameters (Smith and Dukes, 2013), which, in general, cause the rate of leaf photosynthesis to increase with temperature to an optimum and then decrease under higher temperatures. These functional forms are either generic for all C₃ and C₄ species and fixed in time and space or are dependent on a small number of plant functional types (PFTs) but again fixed in time and space. Consequently, climate–carbon feedbacks in ESMs are sensitive to the assumed value of the fixed optimum temperature for photosynthetic capacity (T_opt) because the amount of carbon assimilated depends on whether leaf temperature is dominantly above or below T_opt. Improved process representation of g_s, photosynthesis and its temperature sensitivity in LSMs is necessary to support robust predictions of global climate change via their coupling into ESMs. Modelling studies have shown how photosynthesis and g_s impact climate feedbacks, play a critical role in how climate will change, and strongly influence climate-induced impacts such as water resources (Betts et al., 2007; Cruz et al., 2010; De Arellano et al., 2012; Gedney et al., 2006; Kooperman et al., 2018; Zeng et al., 2017).

This study, therefore, updates the plant physiology routines in the Joint UK Land Surface Environment Simulator (JULES-vn5.6) LSM, the land surface component of the UK Hadley Centre ESM (Sellar et al., 2019). To date, JULES has employed the mechanistic C₃ photosynthesis scheme of Collatz et al. (1991) (“Collatz”). However, the Farquhar et al. (1980) (“Farquhar”) scheme is more generally adopted by those modelling photosynthetic response and by researchers analysing data from empirical studies. The Farquhar scheme has been recently implemented in JULES by Mercado et al. (2018) for C₃ plant types, though a big leaf canopy scaling approach was used, and it was not parameterised and evaluated for global applications. Here we build on that previous study by using a data-driven approach incorporating data from multiple biomes to parameterise the Farquhar model photosynthetic capacity parameters and their temperature sensitivity so it is amenable for use in global studies. Our specific rationale for including the Farquhar photosynthesis scheme is twofold. Firstly, studies by Rogers et al. (2017) and Walker et al. (2021) demonstrate that despite only the Collatz or Farquhar descriptions of leaf photosynthesis being in general use, simulated photosynthesis varies significantly between LSMs. This variation is attributed to several factors, including (1) differences in prescribed Rubisco kinetic constants and their temperature responses (Rogers et al., 2017), (2) structural differences, namely the method used to determine the transition point between the limiting rates of photosynthesis which has a disproportionate impact on estimates (Huntingford and Oliver, 2021; Walker et al., 2021), and (3) the sensitivity of photosynthesis to temperature, in terms of the under-representation of parameters from different biomes to describe the short-term instantaneous response of photosynthesis to temperature (Rogers et al., 2017). In particular, these differences imply that parameter values derived calibrating the Collatz model against data will differ to those derived using Farquhar against the same set of measurements. Parameter values are not transferable between models; hence such differences will lead to inconsistencies and projection errors if parameters are fitted to data but then applied within the alternative model. Building in the capacity of an LSM to run with either photosynthesis scheme greatly enhances flexibility in modelling. Importantly, this flexibility allows for consistency between parameters used by empiricists to derive leaf-level photosynthetic parameters from observations and those used in large-scale modelling. Additionally, our re-parameterisation of the photosynthetic capacity and temperature sensitivity parameters is based on recent global datasets that are more extensive, including species from a range of different biomes, further enhancing the capacity for global modelling applications. Our second rationale is that the Farquhar photosynthesis scheme is required as the underlying model to implement the Kattge and Knorr (2007) thermal acclimation scheme.

Leaf-level g_s response to water vapour is commonly represented in LSMs empirically (Jarvis et al., 1976) or with a semi-empirical model (Ball et al., 1987; Damour et al., 2010; Leuning, 1995). Values of g_s are subsequently scaled, yielding an estimate of canopy conductance for vegetation in different ecosystems. De Kauwe et al. (2013) showed that 10 of the 11 ecosystem models studied in their inter-comparison used a form of the Ball–Berry–Leuning approximation. This model form links g_s to changes in environmental conditions and directly to photosynthetic rate. However, there is increasing interest in using models based on optimisation theory (Franks et al., 2017, 2018), using evidence that stomata may behave to maximise CO₂ gain whilst minimising water loss. The major advantage of optimality theory is that the optimisation criterion will apply under any environmental conditions, past or future. Hence the derived equations can replace uncertain mechanistic formulations and may also have more predictive capability corresponding to future climate regimes. JULES traditionally uses the empirically based Jacobs (1994) g_s scheme (“Jacobs”), and in this study we compare the behaviour of this scheme against the Medlyn et al. (2011) g_s scheme (“Medlyn”), which is based on optimisation theory. The Medlyn g_s model has been previously implemented in JULES by Oliver et al. (2018). However, in this study, we advance on that previous work by calibrating for the increased number of plant functional types now in JULES (nine PFTs, as opposed to five in the original study), and we parameterise using data from a global synthesis of experimental observations.

There is increasing evidence that the short-term vegetation temperature responses are themselves sensitive to temperatures experienced over longer timescales (days to weeks to seasons) and, in particular, have the capability to acclimate to growth temperature (T_growth) (Kattge and Knorr, 2007). Observational evidence of thermal acclimation of photosynthesis has been widely reported, primarily for temperate and boreal ecosystems (Atkin et al., 2006; Gunderson et al., 2000, 2010; Hikosaka et al., 2007; Way and Yamori, 2014; Yamori et al., 2014). The effect is defined as the fast temporal adjustment of the temperature response of photosynthesis driven by a change in T_growth. Thermal acclimation of photosynthesis typically results in a shift in the optimum temperature (T_opt) for photosynthesis towards the new growth temperature, which can result in an increase or maintenance of the photosynthetic rate respective to T_growth (Yamori et al., 2014). In this study, we implement thermal acclimation of photosynthetic capacity in JULES using the scheme from Kattge and Knorr (2007). The scheme attributes all changes in the photosynthetic response to changing T_growth, without specifically separating adaptation from acclimation processes. Of those LSMs that do account for thermal acclimation of photosynthesis (e.g. TEM, CLM4.5, LM3, JULES) (Chen and Zhuang, 2013; Lombardozzi et al., 2015; Mercado et al., 2018; Smith et al., 2016), all similarly use this numerical algorithm from Kattge and Knorr (2007). Mercado et al. (2018) investigated the impacts of thermal acclimation on the future land carbon sink using an implementation of the Kattge and Knorr (2007) in JULES though using a simple big leaf scaling approach. In this study we apply the thermal acclimation scheme in the updated JULES model (i.e. newly parameterised Farquhar scheme, running with a multi-layer canopy and nine PFTs) and updated with the Medlyn g_s scheme and related parameters.

This paper therefore brings together these three key recent developments of the JULES plant physiology routines, (1) implementation of the Farquhar photosynthesis scheme, (2) the optimisation-based Medlyn model of stomatal opening and (3) thermal acclimation of photosynthesis, along with updated parameters and an evaluation of model behaviour. We make incremental additions of the different processes to the JULES model in a set of factorial simulations and run the model with current-day (1979 to 2013) near-surface meteorological forcing and CO₂ levels. First, we present the different factorial simulations in the context of a thorough evaluation of simulated contemporary carbon and energy fluxes. Such evaluation includes comparison against individual eddy covariance sites and at spatial scales up to the global scale against satellite products. Timescales analysed are both seasonal and annual. Secondly, we apply the new model configurations within a past-to-future climate change simulation based on a high-end emissions scenario (Shared Socioeconomic Pathway (SSP5) Representative Concentration Pathway 8.5 (RCP8.5)). We use output from HadGEM3-GC3.1 spanning the years 1960 to 2050 to explore sensitivity of global vegetation to future climate change. This choice of scenario is to allow for eventual comparison between these offline simulations and the equivalent in the coupled global climate model to investigate land–atmosphere feedbacks resulting from these changes to the plant physiology routines. This is currently work being undertaken. This updated version of the JULES model is now available in official JULES releases for use by the community (see “Code and data availability” section). It is therefore also readily available for full coupling into the UK community ESM (UKESM), a process that is just starting.

2 Model description

2.1 JULES land surface model

Our modelling framework is JULES (https://jules.jchmr.org, last access: 5 April 2022), the land surface component of the Hadley Centre climate models, which includes the new UK community Earth System Model (UKESM1) (Sellar et al., 2019). JULES can be run offline, as in this study, forced with observed meteorology, at different spatial scales (from a single location to global). A full description of JULES is provided in Best et al. (2011), Clark et al. (2011) and Harper et al. (2016). Of particular relevance for this study is the plant physiological representation in JULES. JULES uses a leaf-level coupled model of photosynthesis and g_s (Cox et al., 1998) based on Collatz et al. (1991, 1992) (for C₃ and C₄ plants) and Jacobs (1994) respectively. Photosynthesis and g_s are modelled to respond to changes in environmental drivers of temperature, humidity deficit, light, CO₂ concentration and water availability. Soil moisture content is modelled using a dimensionless soil water stress factor, which is related to the mean soil water concentration in the root zone, and the soil water contents at the critical and wilting point (Best et al., 2011). The critical and wilting point soil moisture concentrations vary by soil type in these simulations. In this study, JULES uses a multilayer canopy radiation interception and photosynthesis scheme (i.e. 10 layers) that accounts for vertical variation of incoming direct and diffuse radiation, sunfleck penetration through the canopy, change in photosynthetic capacity with depth into the canopy and inhibition of leaf respiration in the light and differentiates calculation of sunlit and shaded photosynthesis at each layer (Clark et al., 2011; Mercado et al., 2009). The implementation of a multilayer canopy for light interception in JULES was shown to improve modelled canopy-scale photosynthetic fluxes at eddy covariance sites compared to the “big leaf approach” (Blyth et al., 2011; Jogireddy et al., 2006; Mercado et al., 2007). Specifically, the multi-layer approach better captured the light response and diurnal cycles of canopy photosynthesis. While light inhibition of leaf respiration and changing photosynthetic capacity with canopy depth are supported by observations (Atkin et al., 2000, 1998; Meir et al., 2002). Sunfleck penetration through the canopy and the differential effects of direct and diffuse beam radiation on modelled carbon and water exchange in JULES were studied by Mercado et al. (2009). This enabled JULES to reproduce the different light-response curves of GPP under diffuse and direct radiation conditions at both a broadleaf and needleleaf temperate forest.

2.2 Physiology developments

2.2.1 Farquhar photosynthesis for C₃ plants and parameterisation

We implement the Farquhar photosynthesis scheme (Farquhar et al., 1980) to describe the leaf-level biochemistry of photosynthesis for C₃ vegetation following the approach of Mercado et al. (2018). Here the leaf-level photosynthesis is calculated as the minimum (note no smoothing) of two potentially limiting rates (Eq. 1a). These two rates are (i) Rubisco-limited photosynthesis (Eq. 2) and (ii) light-limited photosynthesis with a dependence on the incident photosynthetically active photon flux density and the potential electron transport rate (Eqs. 3 and 4). Note, as in the original Farquhar formulation, we do not include a TPU-limited (triose phosphate utilisation) rate. Further, recent empirical studies suggest that TPU limitation rarely limits photosynthesis under present-day CO₂ concentrations and is also unlikely to limit photosynthesis at elevated CO₂ (Kumarathunge et al., 2019a). This, and the current uncertainty in the formulation of TPU limitation of photosynthesis, led Rogers et al. (2021) to conclude it is an unnecessary complication in LSMs. Hence,

$$\begin{array}{}\text{(1a)}& {\displaystyle }{A}_{\mathrm{p}}=min\left\{A_v,A_j\right\}-R_{\mathrm{d}}\text{(1b)}& {\displaystyle }{A}_{\mathrm{n}}=A_{\mathrm{p}}\mathit{\beta },\end{array}$$

where A_p is the net potential (i.e. unstressed) leaf photosynthetic carbon uptake (mol m² s⁻¹), R_d is the rate of leaf respiration in the dark (mol m² s⁻¹) and A_n is the net photosynthetic rate (mol m² s⁻¹), which accounts for the impact of soil moisture stress on photosynthetic rate by multiplying A_p by the soil water stress factor β. Rubisco-limited photosynthesis (A_v, mol m² s⁻¹) is calculated as in Eq. (2). The maximum rate of carboxylation of Rubisco is determined by V_cmax (mol m² s⁻¹), c_i and o_a are the intercellular concentrations of CO₂ and O₂ (both Pa), K_c and K_o (both units of Pa) are the Michaelis Menten coefficients for Rubisco carboxylation and oxygenation respectively, and Γ (Pa) is the CO₂ compensation point in the absence of mitochondrial respiration.

$$\begin{array}{}\text{(2)}& A_v={\displaystyle \frac{V_{\mathrm{cmax}}\left(c_{\mathrm{i}}-\mathrm{\Gamma }\right)}{\left[c_{\mathrm{i}}+K_{\mathrm{c}}\left(\mathrm{1}+\frac{o_{\mathrm{a}}}{K_{\mathrm{o}}}\right)\right]}}\end{array}$$

The light-limited rate of photosynthesis (A_j, mol m² s⁻¹) (Eq. 3) is a function of the rate of electron transport J (mol m² s⁻¹) which is represented in Eq. (4). J depends on the incident photosynthetically active photon flux density Q (mol quanta m² s⁻¹), the potential rate of electron transport J_max (mol m² s⁻¹), the apparent quantum yield of electron transport α (mol electrons mol⁻¹ photon) fixed at 0.3 (mol electrons mol⁻¹ photon) following Medlyn et al. (2002) and θ a non-rectangular hyperbola smoothing parameter, which takes a value of 0.9 (unitless) following Medlyn et al. (2002). The factor of 4 used in the Farquhar model in Eq. (3) accounts for four electrons being required per carboxylation/oxygenation reaction.

$$\begin{array}{}\text{(3)}& {\displaystyle }{A}_j=\left({\displaystyle \frac{J}{\mathrm{4}}}\right){\displaystyle \frac{\left(c_{\mathrm{i}}-\mathrm{\Gamma }\right)}{\left(c_{\mathrm{i}}+\mathrm{2}\mathrm{\Gamma }\right)}}\text{(4)}& {\displaystyle }\mathit{\theta }{J}^{\mathrm{2}}-\left(\mathit{\alpha }Q+J_{max}\right)J+\mathit{\alpha }Q{J}_{max}=\mathrm{0}\end{array}$$

JULES currently uses Q₁₀ functions in the Collatz scheme to describe the temperature dependency of V_cmax, K_c, K_o and Γ (see Sect. S1 in the Supplement). In our implementation of the Farquhar scheme, temperature sensitivities for the K_c, K_o and Γ are taken from Bernacchi et al. (2001) as described in Medlyn et al. (2002). These are the same temperature sensitivities used by experimentalist to derive estimates of photosynthetic capacity parameters (Rogers et al., 2017). Of particular importance to our analysis here are the temperature responses of V_cmax and J_max. Equation (5) describes the temperature response of both parameters:

$$\begin{array}{}\text{(5)}& k_T=k_{\mathrm{25}}\exp \left[H_{\mathrm{a}}{\displaystyle \frac{\left(T_{\mathrm{l}}-T_{\mathrm{ref}}\right)}{T_{\mathrm{ref}}R{T}_{\mathrm{l}}}}\right]{\displaystyle \frac{\mathrm{1}+\exp \left[\frac{T_{\mathrm{ref}}\mathrm{\Delta }S-H_{\mathrm{d}}}{T_{\mathrm{ref}}R}\right]}{\mathrm{1}+\exp \left[\frac{T_{\mathrm{l}}\mathrm{\Delta }S-H_{\mathrm{d}}}{T_{\mathrm{l}}R}\right]}}.\end{array}$$

Here, k_T (µmol m² s⁻¹) is either V_cmax or J_max at leaf temperature T_l (K), k₂₅ (µmol m² s⁻¹) is the rate of V_cmax or J_max at the reference temperature T_ref of 25 ^∘C (298.15 K), R is the universal gas constant (8.314 J mol⁻¹ K⁻¹), H_a and H_d (J mol⁻¹) are the activation and deactivation energies respectively, and ΔS (J mol⁻¹ K⁻¹) is an entropy term (see Table 1 for PFT-specific parameter values). Broadly, H_a describes the rate of exponential increase in the function below the optimum temperature (T_opt), and H_d describes the rate of decrease above the T_opt. ΔS and T_opt are related by Eq. (6), which is used to calculate the T_opt of V_cmax and J_max (Table 1):

$$\begin{array}{}\text{(6)}& T_{\mathrm{opt}}={\displaystyle \frac{H_{\mathrm{d}}}{\mathrm{\Delta }S-R\ln \left[\frac{H_{\mathrm{a}}}{H_{\mathrm{d}}-H_{\mathrm{a}}}\right]}}.\end{array}$$

Table 1 PFT-specific parameters for the required temperature dependency of V_cmax and J_max in the Collatz and Farquhar photosynthesis schemes. PFT codes (left column) are BET-tr – broadleaf evergreen tropical tree, BET-te – broadleaf evergreen temperate tree, BDT – broadleaf deciduous tree, NET – needle leaf evergreen tree, NDT – needle leaf deciduous tree, C₃ – C₃ grass, C₄ – C₄ grass, ESH – evergreen shrub, DSH – deciduous shrub.

To find new estimates for V_cmax and the $J_{max}:V_{\mathrm{cmax}}$ ratio at T_ref of 25 ^∘C for use with the Farquhar model for the nine PFTs in JULES, we used the global dataset from Walker et al. (2014), which includes data from 356 species. For V_cmax and J_max, Walker et al. (2014) re-analysed the data to remove the variation in these two parameters across studies caused by different parametric assumptions used in their derivation from A−C_i curves (e.g. using a common set of kinetic parameters and reporting values at 25 ^∘C). We calculated the mean V_cmax and J_max across studies conducted at ambient CO₂ concentration for each of the JULES PFTs (Table 2). To parameterise the deciduous needleleaf tree (NDT) PFT, we use the values for the evergreen needleleaf tree (NET) PFT because the data for NDT were from a single study on one juvenile (3 years old) species. An exception was the tropical broadleaf evergreen tree (BET-tr) PFT, where we use V_cmax and J_max from the dataset collated in the more recent compilation by Kumarathunge et al. (2019b), as this study includes many more tropical tree species than any previous meta-analysis.

Parameter values for the temperature response functions for V_cmax and J_max (Eq. 5) in the Farquhar scheme were taken from a global dataset of photosynthetic CO₂ response curves, which entrained data from 141 C₃ species, ranging from the tropical rainforest to Arctic tundra (Kumarathunge et al., 2019b). The study provides parameter values for tree PFTs that match those in JULES, for example, tropical broadleaf evergreen trees (BET-tr PFT in JULES), temperate broadleaf evergreen trees (BET-te), broadleaf deciduous trees (BDTs) and needleleaf evergreen trees (NETs). For the remaining JULES PFTs, BDT values are used for NDT and deciduous shrubs (DSH), and BET-te values are used for evergreen shrubs (ESH). Kumarathunge et al. (2019b) do not include data for C₃ grasses; therefore to parameterise the temperature dependency of V_cmax and J_max for this PFT, we fitted both to the existing V_cmaxtemperature response function in the Collatz scheme for C₃ grasses because of a scarcity of data in the literature. Figure S1 in the Supplement shows the temperature dependency of V_cmax, J_max and gross photosynthesis for Collatz and Farquhar using the PFT-specific parameters in Tables 1 and 2.

Table 2 PFT-specific parameters for the Collatz and Farquhar photosynthesis schemes.

2.2.2 Medlyn model of g_s and parameterisation

In JULES, g_s (m s⁻¹) is represented in Eq. (7).

$$\begin{array}{}\text{(7)}& g_{\mathrm{s}}=\mathrm{1.6}R{T}_{\mathrm{l}}{\displaystyle \frac{A_{\mathrm{n}}}{c_{\mathrm{a}}-c_{\mathrm{i}}}},\end{array}$$

where the factor 1.6 accounts for g_s being the conductance for water vapour rather than CO₂, R is the universal gas constant (J mol⁻¹ K⁻¹), T_l is the leaf surface temperature (K), c_a and c_i (both Pa) are the leaf surface and internal CO₂ partial pressures respectively, and A_n is the net photosynthetic rate. Here, c_i is unknown and is calculated in JULES using the Jacobs scheme as in Eq. (8) and relates the ratio of ambient (c_a) to leaf intercellular (c_i) partial pressure of CO₂ ($c_{\mathrm{i}}/c_{\mathrm{a}}$), to leaf humidity deficit:

$$\begin{array}{}\text{(8)}& c_{\mathrm{i}}=\left(c_{\mathrm{a}}-\mathrm{\Gamma }\right)f_{\mathrm{0}}\left(\mathrm{1}-{\displaystyle \frac{d_q}{d_{q_{\mathrm{crit}}}}}\right)+\mathrm{\Gamma },\end{array}$$

where Γ (Pa) is the CO₂ photorespiration compensation point, d_q is the specific humidity deficit at the leaf surface (kg kg⁻¹), and $d_{q_{\mathrm{crit}}}$ (kg kg⁻¹) and f₀ are PFT-specific parameters representing the critical humidity deficit at the leaf surface and the leaf internal-to-atmospheric CO₂ ratio ($c_{\mathrm{i}}/c_{\mathrm{a}}$) at the leaf-specific humidity deficit (Best et al., 2011). To implement the Medlyn model, Eq. (9) is used to calculate c_i, retaining Eq. (7) to calculate g_s. In Eq. (9), g₁ (kPa^0.5) is a PFT-specific model parameter, and d_q is the vapour pressure deficit at the leaf surface expressed in kilopascals (kPa). The Medlyn scheme is based on optimisation theory and so assumes that stomatal aperture is regulated to maximise carbon gain while simultaneously minimising water loss:

$$\begin{array}{}\text{(9)}& c_{\mathrm{i}}=c_{\mathrm{a}}\left({\displaystyle \frac{g_{\mathrm{1}}}{g_{\mathrm{1}}+\sqrt{d_q}}}\right).\end{array}$$

PFT-specific values of the g₁ parameter were derived for the nine JULES PFTs from the global database of Lin et al. (2015) (Table 3). The g₁ parameter represents the sensitivity of g_s to the assimilation rate, i.e. plant water-use efficiency, and was derived as in Lin et al. (2015), by fitting the Medlyn et al. (2011) model to observations of g_s, photosynthesis and vapour pressure deficit (VPD), assuming an intercept of zero. A non-linear mixed-effect model was used to estimate the model slope coefficient, g₁, for each PFT, where individual species were assumed to be the random effect to account for the differences in the g₁ slope among species within the same group, following Lin et al. (2015).

Table 3 PFT-specific parameters required for the Jacobs and Medlyn g_s schemes.

2.2.3 Thermal acclimation of photosynthetic capacity

The Kattge and Knorr (2007) acclimation algorithm (“AcKK”) is based on the parameters of the Farquhar photosynthesis scheme; hence acclimation is implemented in the Farquhar model. The AcKK algorithm uses empirical relationships to describe the response of V_cmax, J_max and the $J_{max}:V_{\mathrm{cmax}}$ ratio to changes in T_growth (defined in AcKK as the average temperature (day and night) of the previous 30 d), and importantly it represents combined acclimation and adaptation processes. Kattge and Knorr (2007) found that ΔS_v, ΔS_j and the $J_{max}:V_{\mathrm{cmax}}$ ratio decrease linearly with increasing T_growth following Eq. (10). This means according to these relationships, the optimum temperatures (T_opt) of V_cmax and J_max ($T_{{\mathrm{opt}}_v}$ and $T_{{\mathrm{opt}}_j}$) increase by 0.44 and 0.33 ^∘C per degree increase in T_growth respectively, and the $J_{max}:V_{\mathrm{cmax}}$ ratio at 25 ^∘C decreases by 0.035 ^∘C per degree increase in T_growth.

$$\begin{array}{}\text{(10)}& x_i=a_i+b_i{T}_{\mathrm{growth}}\end{array}$$

The x is either ΔS_v, ΔS_j or the $J_{max}:V_{\mathrm{cmax}}$ ratio, and the sub-index i refers to the parameter values (a and b shown in Table 4) for V_cmax, J_max or the $J_{max}:V_{\mathrm{cmax}}$ ratio. T_growth is the growth temperature (calculated online as the mean temperature of the previous 30 d).

Table 4 Parameter values derived by Kattge and Knorr (2007) and used in this study in Eq. (10) to model thermal acclimation of photosynthesis using the AcKK scheme.

3 Model evaluation and application

3.1 Site-level simulations

JULES was applied using four model configurations (Table 5) with observed meteorology and evaluated against data from 17 eddy covariance sites (Table S1, Fig. S2). This collection of eddy covariance measurements represents a range of climates and land cover types (Table S1, Fig. S2). In all simulations the vegetation cover was prescribed, removing any biases that the modelled competition may introduce through self-diagnosis of PFT extents. The prescribed leaf area index (LAI) was used where site data were available, otherwise the JULES phenology scheme was switched on, allowing the LAI to evolve in the model. Model output was evaluated against fluxes of gross primary productivity (GPP) and evaporative fraction (EF). We used EF rather than latent heat flux to minimise issues with incomplete closure of the energy balance (that can typically range from 5 % to 30  % at some eddy covariance sites; Liu et al., 2006). For analysis we used daytime values only (i.e. where the shortwave radiation was >10 W m²) from days with no missing data and compare mean seasonal diurnal cycles of modelled GPP and EF against the observed fluxes. The mean seasonal cycle calculated over the entire measurement period is used in order to assess the mean model behaviour.

We evaluate the site-level simulations with RMSE (root mean square error) for the seasonal diurnal cycle of simulated (daytime) fluxes (GPP and EF). For each site, the time period of the simulation and therefore evaluation period is stated in Table S1. We summarise the changes in RMSE using the relative improvement for each model configuration (i) compared to the current standard JULES configuration of Collatz with Jacobs (Clz.Jac). The statistic is calculated so that positive values show an improvement compared to Clz.Jac and therefore a better comparison to the observations:

$$\begin{array}{}\text{(11)}& {\mathrm{RMSE}}_{{\mathrm{rel}}_i}={\displaystyle \frac{{\mathrm{RMSE}}_{\mathrm{Clz}.\mathrm{Jac}}-{\mathrm{RMSE}}_i}{{\mathrm{RMSE}}_{\mathrm{Clz}.\mathrm{Jac}}}}.\end{array}$$

Table 5 Description of the four model experiments performed both at site level and globally, with the JULES land surface model.

n/a – not applicable.

3.2 Global-scale simulations

Four JULES simulations were performed globally for the period 1979–2013 as outlined in Table 5. These global present-day simulations were run at 0.5^∘ × 0.5^∘ spatial resolution. The WFDEI meteorological dataset was used to drive the model (Weedon et al., 2014). This has a 3 h temporal resolution that JULES interpolated down to an hourly model time step. To focus on the direct effects of the model changes on GPP and surface energy fluxes, the land surface properties of the model were prescribed. We use a static map of land cover (in terms of different PFT extents) derived from the European Space Agency's Land Cover Climate Change Initiative (ESA LC_CCI) global vegetation distribution version 1.6 for the 2010 epoch (Poulter et al., 2015) (Fig. S3) following that used in Harper et al. (2016). Seasonally varying LAI values were derived from the Global LAnd Surface Satellite (GLASS) dataset (Xiao et al., 2016). Prescribed parameters were used for the hydraulic and thermal properties of the soil from a modified version of the H1 lookup table from Zhang and Schaap (2017) that depends upon the soil textural type from SoilGrids (Hengl et al., 2014). We also prescribe transient atmospheric CO₂ concentrations based on annual mean observations from Mauna Loa (Tans and Keeling, 2014). A spin-up of 80 years was performed (recycling through the period 1979 to 1999), which is sufficient to equilibrate soil temperature and soil moisture.

The global offline present-day simulations were compared against the global evaluation products, and for both model output and observations, we calculate seasonal means over the period 2002 to 2012. We used the global FluxCom product to evaluate modelled GPP, LE, H and ET (Jung et al., 2020; Tramontana et al., 2016). We compare our simulations against the FluxCom ensemble product (RS + MET) driven with the same forcing (WFDEI), as is recommended by Jung et al. (2019) to minimise deviations due to different climate input data. To convert LE to ET, we assume a constant latent heat of vaporisation of 2.5 MJ mm⁻¹. We also use the model-derived product from GLEAM-v3.3a to evaluate ET and additionally use the MODIS GPP product (Zhao et al., 2005; Zhao and Running, 2010; Zhao et al., 2006) to evaluate simulated global GPP.

Global future climate simulations were performed forced with meteorological output (1960 to 2050) from the HadGEM3-GC3.1 model atmosphere-only simulations at 3 h temporal resolution and N512 spatial resolution (Roberts et al., 2019; Williams et al., 2018). These projections follow the CMIP6 HighResMIP protocol (Haarsma et al., 2016). This choice of forcing to drive JULES is to allow for comparison of the offline runs performed in this study with the equivalent simulations currently being undertaken in the coupled HadGEM3-GC3.1 model to explore land–atmosphere feedbacks arising from changes implemented in the plant physiology routines in this work. The factorial set of offline simulations in this work provide a systematic sensitivity study that is less computationally expensive with which to help understand behaviour seen in the coupled model. The output at N512 was re-gridded to 0.5^∘ × 0.5^∘ using conservative interpolation, which ensures the physical conservation of each variable. Figure S4 shows the mean temperature and precipitation change by region over the study period and the atmospheric CO₂ concentration. Atmospheric CO₂ concentrations were prescribed based on observations up to 2014 as described in historical CMIP6 simulations (Eyring et al., 2016). From 2015 onwards, atmospheric CO₂ concentrations were based on a high-end emission scenario of the Shared Socioeconomic Pathways (SSP5) with the Representative Concentration Pathway 8.5 (RCP8.5) (Haarsma et al., 2016). As for the current-day simulations, LAI, land cover and soil properties were prescribed using the same datasets. A spin-up period of 80 years (recycling through the period 1960 to 1980) was again used to equilibrate soil temperature and soil moisture.

We analyse the future global simulations using the “difference of difference” approach. This method explicitly targets the change in the variable of interest over the study period resulting from the change in process alone and negates differences that may arise from different initial starting points of each simulation (different initial conditions):

$$\begin{array}{}\text{(12)}& \mathrm{Effect}=\left({\stackrel{\mathrm{‾}}{X}}_{\mathrm{2050}}-{\stackrel{\mathrm{‾}}{X}}_{\mathrm{1980}}\right)-\left({\stackrel{\mathrm{‾}}{Y}}_{\mathrm{2050}}-{\stackrel{\mathrm{‾}}{Y}}_{\mathrm{1980}}\right),\end{array}$$

where X represents the simulation with the process of interest, and Y represents the simulation with the alternative representation, and 2050 and 1980 represent the end and start of the simulation analysis period respectively (calculated as the mean over 2040 to 2050 and 1980 to 1990 respectively). For example, to look at the impact of changing photosynthesis schemes, $X=\mathrm{Fq}.\mathrm{Jac}$ and $Y=\mathrm{Clz}.\mathrm{Jac}$. In this case, both configurations are using the Jacobs g_s scheme; only the photosynthesis scheme changes from Collatz to Farquhar. The impact of changing g_s scheme is assessed, where $X=\mathrm{Fq}.\mathrm{Med}$ and $Y=\mathrm{Fq}.\mathrm{Jac}$. The impact of thermal acclimation is assessed where $X=\mathrm{AcKK}.\mathrm{Med}$ and $Y=\mathrm{Fq}.\mathrm{Med}$; here both simulations use the Farquhar photosynthesis scheme and the Medlyn g_s scheme, but X has the addition of thermal acclimation of photosynthesis.

4 Results

4.1 Site-level evaluation

Results from the FLUXNET sites comparing the mean seasonal diurnal cycles of GPP and EF against observed fluxes are summarised in Fig. 1, where reds and yellows indicate reduced RMSE relative to the “standard” JULES configuration of Collatz with Jacobs (Clz.Jac) and therefore closer agreement to site-level FLUXNET observations. Results are variable by site and season (Figs. 1, S5 and S6), some of which will be due to other site-specific characteristics that are not simulated well by the model, such as LAI for those sites that rely on model-derived estimates. On the other hand, soil properties are prescribed by parameters that describe the thermal and hydraulic characteristics of the soil; uncertainties in these parameterisations have consequences for the simulated soil moisture content at each site, for example, which impacts simulated carbon and water fluxes. We first consider results for the five tropical sites. Results are mixed for the simulated seasonal diurnal cycle of GPP at the tropical (EBF $/$ BET-tr) sites; GPP is improved (reduced) with the new JULES model configurations at three out of the five tropical sites in March–April–May (MAM; Figs. 1a, S5), with thermal acclimation leading to the greatest improvements. However in June–July–August (JJA; Figs. 1b, S5), this improvement is only found at two of the tropical sites. At the EBF sites, implementing the Farquhar photosynthesis model means V_cmax is lower (BET-tr, Table 2), and this in addition to the change in temperature sensitivity (Table 1; Fig. S1a–c) and model structural changes from Collatz to Farquhar results in lower simulated GPP compared to Collatz. Thermal acclimation allows for further adjustments of the $T_{{\mathrm{opt}}_v},T_{{\mathrm{opt}}_j}$ and the $J_{max}:V_{\mathrm{cmax}}$ ratio which results in lower simulated photosynthesis and therefore GPP compared to Farquhar (Fig. S5). The change from Jacobs g_s model to Medlyn has minimal impact on simulated GPP for the tropical tree PFT because in both schemes the modelled c_i has a similar sensitivity to humidity deficit at the leaf surface, with the exception at very low humidity deficit (Figs. S7; S5). The simulated seasonal diurnal cycle of EF is improved (reduced) at four out of the five tropical sites in both MAM and JJA, again with some of the largest improvements seen with thermal acclimation (Figs. 1c, d; S6).

Figure 1 Relative changes in RMSE for each JULES model configuration compared to Collatz with Jacobs (Clz.Jac) for hourly daytime (a) GPP (March–April–May), (b) GPP (June–July–August), (c) EF (March–April–May) and (d) EF (June–July–August). Calculated according to Eq. (11), positive values (reds and yellows) are where RMSE is lower compared to the Clz.Jac configuration and therefore indicates an improvement compared to the Clz.Jac baseline, and the FLUXNET observations. EBF – broadleaf evergreen tropical tree, GRA – C₃ grassland, BDT – broadleaf deciduous tree, NET – needle leaf evergreen tree. The fit of each model configuration to observations and the RMSE are shown in Fig. S5 (GPP) and Fig. S6 (EF).

At the C₃ grassland sites (GRA), improved simulated GPP (higher GPP) is seen across all sites in JJA with the Medlyn g_s scheme and thermal acclimation (Figs. 1b, S5). This is matched by improvements in simulated EF (higher EF) across all grassland sites in both seasons, with the exception of US_var in JJA (Figs. 1c, d; S6). The change from Collatz to Farquhar at the GRA sites means a lower V_cmax is used (C₃, Table 2), although the temperature sensitivity is similar (Table 1, Fig. S1p, q); this results in lower GPP simulated by Farquhar compared to Collatz, which compares worse to the observations (GPP and EF, Figs. 1, S5). In contrast to using Farquhar with the Jacobs g_s scheme, using Farquhar with the Medlyn scheme improves simulated GPP and EF; both are increased because for the C₃ grass PFT, as the humidity deficit at the leaf surface increases, c_i simulated by Medlyn is less sensitive compared to Jacobs (Figs. S7; S5), leading to higher c_i, higher net canopy photosynthesis and GPP, and higher transpiration and LE. These results suggest the Medlyn scheme has a large impact on simulated carbon and water fluxes for the C₃ grass PFT in the JULES model. In JJA, the adjustment of the temperature sensitivity of photosynthesis to the T_growth by the thermal acclimation scheme tends to increase GPP compared to Farquhar with no acclimation, and this compares better to the observations (Figs. 1, S5).

At the broadleaf deciduous tree (BDT) sites, simulated GPP is improved with all JULES model configurations in MAM (higher GPP) at three out of the four sites (Fig. 1a). However in JJA, improvements are mainly seen with thermal acclimation (lower GPP compared to Fq.Med; Fig. 1b). Medlyn g_s performs worse at all sites in JJA, suggesting either the model formulation or parameters are not suitable to correctly capture stomatal behaviour in this season for this PFT (Figs. 1b, S5). Compared to Collatz, the Farquhar model for the BDT PFT uses a lower V_cmax (Table 2) and has a considerably lower $T_{{\mathrm{opt}}_v}$ (Table 1; Fig. S1h), which means that at leaf temperatures below ∼ 22 ^∘C, photosynthesis is higher with the Farquhar model and above this photosynthesis is lower than Collatz (Fig. S1g). Consequently, warmer temperatures in JJA lead to lower GPP simulated by Farquhar compared to Collatz, and cooler temperatures in MAM result in slightly higher GPP with Farquhar compared to Collatz (Fig. S5). Using the Medlyn model means simulated c_i is more sensitive to increasing leaf humidity deficit for the BDT PFT (Fig. S7). Medlyn simulates a lower c_i as humidity deficit increases compared to Jacobs, which leads to lower GPP and LE, the magnitude of which depends on the local site humidity conditions. In JJA the Medlyn g_s model performs worse at all sites for GPP (Fig. 1b), although improvements in simulated EF are seen in JJA, where both Medlyn and thermal acclimation improve model performance at three out of four BDT sites (Figs. 1d, S6).

At the evergreen needleleaf (NET) sites, the most consistent improvements to simulated GPP are seen with the Farquhar model, where simulated GPP in JJA is substantially improved (GPP reduced) at three out of four sites (Figs. 1b, S5); in this season both Medlyn and thermal acclimation generate larger improvements in the simulated GPP (reducing GPP further), but this is just at two out of the four sites. In our implementation of the Farquhar model, the NET PFT has a lower V_cmax compared to Collatz (Table 2) and a slightly higher $T_{{\mathrm{opt}}_v}$ (Table 1, Fig. S1k). The resulting shape of the temperature response curve for photosynthesis (Fig. S1j) means that at leaf temperatures below ∼ 10 ^∘C, Farquhar photosynthesis is higher. However, above 10 ^∘C, Farquhar photosynthesis is lower compared to Collatz, resulting in simulated GPP in MAM that tends to be higher with Farquhar than Collatz, and in JJA the opposite occurs (Fig. S5). In MAM and JJA, the Medlyn g_s model simulates some large improvements in EF; c_i simulated by Medlyn is more sensitive to increasing leaf humidity deficit compared to Jacobs (Fig. S7), which results in lower transpiration and EF, and this compares better to the observations (Figs. 1, S6).

4.2 Global evaluation

4.2.1 Spatial differences between model configurations

The impact of changes in the photosynthesis scheme, g_s scheme, adding thermal acclimation of photosynthetic capacity and the overall change on simulated GPP, LE and H, is shown in Fig. 2 by comparing each of the new JULES configurations to the configuration with the alternative process representation. For GPP, the biggest change is moving from the Collatz photosynthesis scheme to the Farquhar photosynthesis scheme (Fig. 2a). Most notably, this change results in decreased GPP in the tropical region in JJA of up to 1.5 g C m² d⁻¹ (up to 10 % reduction), whilst in the high northern latitudes, GPP is increased by up to 1.5 g C m² d⁻¹ (up to 20 % increase). This is consistent with results from the site-level simulations where GPP was reduced with implementation of the Farquhar model at tropical sites and increased in cooler months (MAM) at the evergreen needleleaf forest sites (here increased GPP in NET-dominated areas are in the forests of the high northern latitudes, which is consistent with cooler temperatures). Impacts on LE and H resulting from the move from Collatz to Farquhar are not as extensive as those seen with GPP (Fig. 2b and c). The change from Jacobs g_s scheme to Medlyn impacts LE and H most, resulting in a pronounced pattern of decreased LE in northern latitudes (up to 10 W m², equivalent to a 10 % reduction) and corresponding increase in H in JJA (Fig. 2e and f). In these JULES simulations, this region is dominated by NET forest, and the high-latitude changes are consistent with results from the site-level simulations, where using the Medlyn g_s scheme at NET sites resulted in some of the biggest improvements in simulated EF (lower LE and therefore lower EF). Including thermal acclimation of photosynthesis has the most extensive impacts on simulated GPP in contrast to LE and H. In the tropical forests, GPP is reduced by up to 1 g C m² d⁻¹ (between 2 % and 5 % reduction) in JJA (Fig. 2g). The impact of acclimation is spatially variable in the temperate region in JJA, with GPP decreased in Europe (between 2 % and 5 %) but increased in eastern United States (up to 20 %). Some areas of the boreal region see increased GPP (between 2 % to 5 %). This GPP response demonstrates the impact of thermal acclimation, which allows for the parameters of the temperature sensitivity functions for photosynthetic capacity (V_cmax, J_max and $J_{max}:V_{\mathrm{cmax}}$) to move in response to the temperature of the growth environment, leading to spatially and temporally different values of the T_opt for photosynthesis for each C₃ PFT. Thermal acclimation impacts LE and H to a lesser extent, but where changes are seen, acclimation increases LE with a corresponding decrease in H (Fig. 2h and i). Figure 2j, k and l show the overall change that results from moving from the traditional JULES set-up of Collatz with Jacobs (Clz.Jac) to Farquhar with thermal acclimation and Medlyn g_s (AcKK.Med), and the impacts on simulated GPP, LE and H can clearly be seen as the trade-off between the dominating effects from each model configuration. For LE and H the response of the simulated energy fluxes is dominated by the change in the representation of g_s, and for GPP the response of simulated carbon fluxes is dominated by the change in the representation of photosynthesis and its response to temperature (i.e. thermal acclimation).

Figure 2 Absolute difference between JULES-modelled GPP, latent (LE) and sensible heat (H) for the different JULES model configurations in June–July–August (JJA) to show the impact of (a, b, c) changing photosynthesis scheme (Fq.Jac–Clz.Jac), (d, e, f) changing g_s scheme (Fq.Med–Fq.Jac), (g, h, i) accounting for thermal acclimation of photosynthesis (AcKK.Med–Fq.Med) and (j, k, l) the overall change (AcKK.Med–Clz.Jac), under present-day meteorological conditions. For each variable, the mean over the period 2002 to 2012 is used. The absolute mean value simulated by each model configuration (JJA) is shown in Fig. S8. DJF is shown in Fig. S9 (mean absolute values) and Fig. S10 (absolute difference).

4.2.2 Comparison to global estimates: seasonal mean GPP and ET

Evaluation of simulated global mean GPP by season using FluxCom and MOD17 global GPP products is presented in Fig. 3a and using global ET from both FluxCom and GLEAM is shown in Fig. 3b. The seasonal means show thermal acclimation compares best to observations (FluxCom) in JJA (AcKK.Med underestimates GPP by just 4 %, whereas Clz.Jac underestimates GPP by 6 %; Fig. 3a and Table S2) and MAM (AcKK.Med underestimates GPP by just 5 %, whereas Clz.Jac underestimates GPP by 11 %; Fig. 3a and Table S2) and is in reasonable agreement with FluxCom in DJF (AcKK.Med overestimates GPP by just 2 %, whereas Clz.Jac underestimates GPP by 4 %; Fig. 3a and Table S2). All JULES model configurations have a high GPP bias in SON compared to FluxCom, and in all seasons GPP is overestimated by all model configurations compared to MOD17; similarly this is largest in SON. For simulated ET, seasonally the model performance is very similar between the different JULES configurations but in both SON and DJF Medlyn (Fq.Med) compares better to both FluxCom and GLEAM, but the differences are very small (Fig. 3b and Table S3).

Figure 3 Seasonal mean global (a) GPP and (b) ET for each JULES model configuration compared to FluxCom (closed symbols) and MOD17 (GPP) or GLEAM (ET) (open symbols).

4.2.3 Comparison to global estimates: latitudinal mean GPP and ET

Figures 4 and 5 present comparisons of seasonal zonal-mean GPP and ET respectively. Firstly, Figs. 4 and 5 highlight the differences between global products used to evaluate GPP and ET (see, for example, Spafford and MacDougall, 2021). For example, FluxCom generally predicts higher GPP in the tropics compared to MOD17, especially in DJF and MAM, and in JJA the different distribution of GPP by latitude means in the tropics MOD17 GPP is higher than FluxCom in the southern latitudes, and FluxCom GPP is higher in the northern tropics. Comparison of the two ET products shows that GLEAM tends to give higher ET in the tropics, particularly in DJF and MAM. Bearing in mind uncertainties in observation-based estimates of fluxes at this scale, we now consider how the different model configurations compare. Notably, all the JULES model configurations in this study simulate comparable global carbon and water fluxes for the recent contemporary period and are in reasonable agreement with the global products used for evaluation. Differences in RMSE between the different model configurations are small for both GPP and ET. Importantly, the most consistent change is the improvement (lowest RMSE) of modelled GPP in the tropics with the Farquhar model (Fq.Jac). This improvement is evident in all seasons and holds when comparing to both FluxCom and MOD17 (Fig. 4). Similarly, estimates of ET are improved in the tropics (lowest RMSE) with the Farquhar model (Fq.Jac) in DJF and JJA and with the Medlyn model (Fq.Med) in MAM and SON, and again this result is not dependent on the choice of observation-based product (Fig. 5). Another notable change is the improvement of simulated GPP in the temperate north and boreal regions in MAM with thermal acclimation (AcKK.Med). Deficiencies in the model stand out, but these biases are common to all model configurations. For example, all configurations simulate an over-prediction of GPP and ET in SON in the temperate north and boreal regions, overestimated GPP in MAM in tropical southern latitudes (0 to −20^∘ S), under-predicted GPP and ET in MAM in temperate north and boreal regions, and an over-prediction of ET in MAM in the temperate and tropical south.

Figure 4 Mean (2002 to 2012) GPP (g C m² d⁻¹) by latitude band and season for each JULES model configuration compared to the FluxCom and MOD17 global GPP products. The bars along the side indicate which model configuration gives the lowest RMSE and therefore better comparison to FluxCom (righthand bar) and MOD17 (left-hand bar) derived GPP for each region. RMSE values are shown in Tables S4 (FluxCom) and S5 (MOD17). The grey shaded area shows the uncertainty in the FluxCom GPP product, provided as the median absolute deviation of ensemble members, this is scaled to a robust estimate of the standard deviation of a normal distribution by multiplying by 1.4826 according to Jung et al. (2019).

Figure 5 Mean (2002 to 2012) evapotranspiration (ET mm d⁻¹) by latitude band and season for each JULES model configuration compared to the FluxCom and GLEAM global ET products. The bars along the side indicate which model configuration gives the lowest RMSE and therefore better comparison to FluxCom (righthand bar) and GLEAM (left-hand bar) derived ET for each region. RMSE values are shown in Table S6 (FluxCom) and Table S7 (GLEAM). The grey shaded area shows the uncertainty in the FluxCom ET product, provided as the median absolute deviation of ensemble members, this is scaled to a robust estimate of the standard deviation of a normal distribution by multiplying by 1.4826 according to Jung et al. (2019).

4.2.4 Comparison to global estimates: spatial variability of mean GPP and ET

The spatial variability of simulated GPP and ET is shown in Fig. 6 during JJA (Fig. S11 for DJF). We show which of the JULES model configurations gives the lowest RMSE compared to observation-based estimates of GPP and ET from FluxCom, MODIS and GLEAM (actual RMSE in Figs. S12 and S13). The differences in RMSE are typically small between the different JULES model configurations; however some clear patterns emerge. Figure 6a and b show that in the tropical forests of the Amazon basin, central Africa and Southeast Asia (Indonesia, Papua New Guinea, Malaysia), in both JJA and DJF (Fig. S11a and b for DJF), GPP simulated including thermal acclimation (AcKK.Med) compares best to both FluxCom and MOD17 across large spatially consistent areas. Outside of these areas, Fq.Jac also improves the simulation of GPP in the tropics, as does the Medlyn g_s model (Fq.Med) in JJA in South China and Indo-China. Also, in the high northern latitudes, dominated by evergreen needleleaf forests, inclusion of thermal acclimation more closely aligns simulated GPP with both FluxCom and MOD17 (Fig. 6a and b). Compared to FluxCom, ET in JJA is simulated best by thermal acclimation (AcKK.Med) in the northern temperate and boreal region, although this pattern is not consistent in comparison to GLEAM (Fig. 6c and d). In contrast to GPP, results are more mixed in the tropics for ET. In areas dominated by tropical tree cover, thermal acclimation (AcKK.Med) and Medlyn (Fq.Med) tend to give the lowest RMSE in JJA and DJF, and in tropical areas dominated by C₃ and C₄ grasses Farquhar (Fq.Jac) performs best (Fig. 6c and d), although in DJF the Medlyn model gives the lowest RMSE in these areas (Fig. S11c and d). In DJF for both GPP and ET, in northern temperate and boreal regions the Collatz with Jacobs (Clz.Jac) configuration performs the best (Fig. S11).

Figure 6 Colours indicate the JULES model configuration that gives the lowest RMSE compared to either the (a) FluxCom and (b) MOD17 global GPP (g C m² d⁻¹) products or (c) FluxCom and (d) GLEAM global ET (mm d⁻¹) products for JJA over the period 2002 to 2012. Actual RMSE values shown in Figs. S12 and S13.

4.3 Application under future climate

We run the new configurations forced by variables from a future climate scenario (HadGEM3-GC3.1 forcing under a high-end emission scenario of the SSPs) to investigate the response of simulated fluxes to long-term warming. Changing the photosynthesis scheme from Collatz to Farquhar results in lower GPP (up to 30 % decrease) by 2050 across the high-northern-latitude forests (Fig. 7a), with the impact on LE (decreased) and H (increase) less extensive (Fig. 7b and c). This area is dominated by NET, NDT and BDT PFTs in JULES. The different temperature sensitivity of photosynthesis parameterised with the Farquhar model compared to Collatz (Fig. S1g, j and m) means at lower leaf temperatures, photosynthesis is higher with Farquhar; however, as leaf temperature increases, photosynthesis falls in Farquhar relative to Collatz. The crossover point at which this occurs is relatively low for these PFTs, particularly NET. This impact of the change in temperature sensitivity was seen in the site-level simulations at FLUXNET NET and BDT sites. There, modelled GPP tended to be higher with Farquhar than Collatz in MAM but lower in the warmer conditions of JJA, and in this climate change scenario, the temperate and boreal region both experience large increases in mean annual air temperature (+5 ^∘C from 1980 to 2060; Fig. S4a and c).

Figure 7 The difference of difference approach (Eq. 12) to determine the impact on GPP (g C m² d⁻¹), LE and H (both W m²) of the individual changes to each JULES model configuration over the course of the future (HadGEMGC3.1) simulation (1980 to 2050) in June–July–August (JJA). The AcKK.Med acclimation effect is calculated from Fig. S16 AcKK.Med–Fq.Med, the effect of the Medlyn g_s scheme is calculated from Fig. S16 Fq.Med–Fq.Jac and the effect of the photosynthesis scheme is calculated from Fig. S16 Fq.Jac–Clz.Jac.

Replacing the Jacobs g_s scheme with Medlyn has the biggest impact on the surface energy fluxes, with increased LE of up to 30 % and a corresponding decrease in H by 2050 across the temperate region (Fig. 7e and f). This area is dominated by the C₃ grass PFT in JULES, which has a less conservative water-use strategy in the Medlyn scheme (high g₁) compared to Jacobs. This means in the Medlyn scheme, the C₃ grass PFT is less sensitive to increasing humidity deficit at the leaf surface; therefore as humidity deficit increases, Medlyn simulates higher c_i, leading to a higher rate of transpiration and LE compared to Jacobs (Fig. S7).

Thermal acclimation of photosynthesis leads to widespread increases in GPP by 2050 (Fig. 7g). This amounts to 10 % in the tropical forests, up to 30 % in northern temperate and boreal regions, and up to 40 % in Southeast Asia. In this long-term climate change scenario, with large increases in mean annual temperature (Fig. S4), the impact of thermal acclimation on GPP can clearly be seen. The flexibility in $T_{{\mathrm{opt}}_v}$, $T_{{\mathrm{opt}}_j}$ and the $J_{max}:V_{\mathrm{cmax}}$ ratio of photosynthesis that thermal acclimation allows through letting these parameters move with the prevailing T_growth allows for higher rates of photosynthesis and therefore GPP as temperatures increase. By contrast, in simulations where photosynthetic rates are controlled by fixed temperature sensitivities, vegetation may have moved past its thermal optimum. Time series of the area-weighted mean annual GPP show that in this simulation, across the tropical region, thermal acclimation enhances GPP by ∼ 7.5 PgC compared to no acclimation (Fig. 8a). In the temperate region and sub-tropics, thermal acclimation increases GPP by ∼ 1 PgC by 2050 (Fig. 8b and d), and in the boreal region, GPP is enhanced by ∼ 0.4 PgC (Fig. 8c). Thermal acclimation of photosynthesis also has a large impact on simulated energy fluxes, most notably in the northern temperate region, where LE is increased by up to 50 % to 60 % (decreased H up to 40 % to 50 %) (Fig. 7h and i).

Figure 8 Time series of the regional mean acclimation effect, i.e AcKK.Med–Fq.Med (black), and the effect of the Medlyn g_s model, i.e. Fq.Med–Fq.Jac (grey).

5 Discussion

Photosynthesis and g_s are central to the estimate of carbon and water fluxes in LSMs, and, when coupled in ESMs, these processes feedback onto the climate system to influence predictions of future climate change. Therefore improving the representation of these processes in LSMs is important, and previous studies have identified thermal acclimation of photosynthesis as a key missing process (Booth et al., 2012).

5.1 Performance of the new JULES plant physiology model configurations: thermal acclimation

Our results show that including thermal acclimation of photosynthesis in the JULES model improves simulated carbon and water fluxes in several key areas for the recent contemporary period. Firstly, the seasonal mean estimates of global GPP show that in most seasons (JJA, MAM and DJF), thermal acclimation of photosynthesis with Medlyn g_s (AcKK.Med) predicts GPP in closer agreement with estimates from FluxCom compared to the traditional “standard” JULES configuration of Collatz photosynthesis with Jacobs g_s (Clz.Jac). Secondly, thermal acclimation with Medlyn g_s improves the simulation of GPP (reduces GPP) in the tropical forests in JJA and DJF (i.e. the Amazon basin and central African rainforest region) and is in closest agreement with estimates of GPP from both FluxCom and MOD17 for these regions. Thirdly, in the high-northern-latitude forests dominated by evergreen needleleaf trees, thermal acclimation increases GPP in JJA and is again in closest agreement with the observational estimates. Finally, in JJA, AcKK.Med improves the simulation of ET across a large area of the temperate north and boreal regions.

Our evaluation therefore suggests that fixed, PFT-specific temperature dependencies for V_cmax (and J_max) do not accurately simulate GPP for the tropical tree and evergreen needleleaf tree PFTs for the present day in the JULES model. Thermal acclimation allows the temperature sensitivity of photosynthesis to adjust to the local temperature environment through flexibility in $T_{{\mathrm{opt}}_v}$, $T_{{\mathrm{opt}}_j}$ and the $J_{max}:V_{\mathrm{cmax}}$ ratio. In the tropical forests, for example, GPP is overestimated by both Clz.Jac and Fq.Jac. The configuration with thermal acclimation reduces GPP compared to both these model configurations. From the leaf-level plots in Fig. S1a, the fixed T_opt of photosynthesis in the Collatz scheme is ∼ 33 ^∘C and in Farquhar is ∼ 34 ^∘C. This is higher than observations from Fig. 1a of Kumarathunge et al. (2019b), where the T_opt for net leaf photosynthesis lies between ∼ 29 and 32 ^∘C, and other studies also show a lower T_opt for photosynthesis of around 30 ^∘C for mature tropical trees (Hernández et al., 2020; Mau et al., 2018). This supports our results and suggests the fixed temperature sensitivity of photosynthesis for tropical trees in the JULES model results in a T_opt of photosynthesis that is too high for the current day. Thermal acclimation results in a more realistic T_opt of photosynthesis for tropical trees because it is influenced by actual growth temperature and so can adjust to local environmental conditions.

Under the climate change scenario used in this study, thermal acclimation shows a sustained positive acclimation effect in all regions, increasing GPP in response to long-term warming (although this is less pronounced in the boreal region). By 2050 GPP was ∼ 10 % higher with thermal acclimation in the tropical forests, up to 30 % to 40 % higher across a large area of the Northern Hemisphere. Our findings broadly agree with Mercado et al. (2018), who implemented the Kattge and Knorr (2007) thermal acclimation scheme into JULES running as part of a coupled climate–carbon model and found that thermal acclimation increased land carbon storage in tropical and temperate regions. This is in contrast to Lombardozzi et al. (2015) and Smith et al. (2016), whose studies both found a negative impact of photosynthetic thermal acclimation in the tropics, again using the Kattge and Knorr (2007) thermal acclimation scheme. Mercado et al. (2018) attribute these differences to the method used to implement acclimation of the $J_{max}:V_{\mathrm{cmax}}$ ratio at 25 ^∘C, that is either reducing J_max alone as in the case of the latter two studies or by decreasing J_max and increasing V_cmax simultaneously whilst keeping the total amount of leaf nitrogen the same as used in the present study and in Mercado et al. (2018). The simulated response of thermal acclimation therefore appears to be sensitive to this subtlety in the parameterisation of the acclimation schemes and warrants further investigation. Yet a clear understanding of what drives the change in the $J_{max}:V_{\mathrm{cmax}}$ ratio in response to T_growth is still lacking. More recent results from the analysis by Kumarathunge et al. (2019b) highlight the difficulty in pinning down what drives this process. They found that the $J_{max}:V_{\mathrm{cmax}}$ ratio responded strongly and consistently to T_growth, but whether that was achieved by increasing V_cmax, decreasing J_max or both was highly variable.

The behaviour of the thermal acclimation scheme in JULES in response to long-term warming implies unlimited thermal resilience of vegetation, but how realistic is this? Observational studies suggest temperate tree species have sufficient capacity to acclimate to rising temperatures (e.g. Drake et al., 2015; Reich et al., 2018; Sendall et al., 2015), although large inter-specific variability in thermal tolerance is identified in co-occurring temperate tree species (Guha et al., 2018). Studies exploring thermal acclimation of photosynthesis for grasslands and C₃ herbaceous vegetation are more limited. For boreal tree species, experimental studies suggest high variability between species with respect to photosynthetic acclimation responses to increasing temperatures; for example, there is an increasing body of work suggesting that the evergreen boreal conifer species Picea might be particularly vulnerable to warming (Benomar et al., 2017; Dusenge et al., 2020; Kroner and Way, 2016; Kurepin et al., 2018; Way and Sage, 2008; Zhang et al., 2015). The 3-year open-air warming experiment of Reich et al. (2018) showed that for 11 temperate and boreal tree species studied, warming increased photosynthesis in most species on wet soils but not in drier conditions. Further, under moist soil conditions, all deciduous species showed an acclimation response to increased temperatures; however, the two boreal evergreen species, Abies and Picea, showed no thermal acclimation response at any soil moisture concentration. It is generally thought that evergreen species have a reduced capacity to acclimate growth and photosynthesis to warming compared to deciduous tree species (Dusenge et al., 2020; Way and Yamori, 2014). Therefore, the response of boreal forest ecosystems to warming will depend on species composition given the varied acclimation capacities shown and lower diversity of boreal forests and, as Reich et al. (2018) highlight, also on interaction with other climate changes such as precipitation. In contrast to temperate and boreal forests, tropical forests are thought to be more susceptible to climate change, having evolved under relatively narrow temperature regimes and experiencing less seasonal and day-to-day variation in temperature changes (Cunningham and Read, 2003). As a consequence, an increasing number of studies show that tropical trees have less capacity to physiologically acclimate photosynthesis to increasing temperatures (Carter et al., 2021; Dusenge et al., 2021; Mau et al., 2018; Miller et al., 2021; Vårhammar et al., 2015). Other studies have determined high temperature threshold responses of photosynthesis, indicating an ability of tropical trees to acclimate to moderate warming, but more severe warming decreases carbon gain (Doughty and Goulden, 2008; Pau et al., 2018; Slot and Winter, 2017; Sullivan et al., 2020). In two tropical understorey species, acclimation of the T_opt of photosynthesis was observed in the early successional species, whereas no acclimation capacity was shown by the mid-successional species (Carter et al., 2020). Our study demonstrates a large positive impact of thermal acclimation on GPP in tropical forests. However a notable uncertainty in the parameterisation is that the dataset used in the Kattge and Knorr (2007) scheme to construct the empirical relationships is heavily weighted towards temperate species, including only two boreal species and no tropical species (Kattge and Knorr, 2007). There is a significant gap in understanding tropical forest responses to increasing temperature. Observational studies are starting to address this gap, but this increasing knowledge is yet to be incorporated into models. Therefore, whilst results from this study demonstrate the importance of thermal acclimation of photosynthesis on simulation of the future global carbon cycle, they should be interpreted with some caution. The varied results from experimental studies highlight the research needed to further understand thermal acclimation responses in a variety of ecosystems, over different timescales and from leaf level through to the canopy and finally to translate that understanding so it is amenable to incorporation into ESMs.

5.2 Performance of the new JULES plant physiology model configurations: Medlyn g_s

In this study, the Medlyn g_s model had the biggest impact on surface energy fluxes simulated by the C₃ grass PFT and needleleaf evergreen tree PFT in JULES. This reflects a change to the water-use strategy of these PFTs as reported by Lin et al. (2015) that is not currently captured by parameterisations in the JULES Jacobs model. Global simulations with the Medlyn scheme for the recent contemporary period simulated a ∼ 10 % decrease in LE (increased H) across the high northern latitudes dominated by the NET PFT compared to the standard JULES Jacobs g_s scheme. The future climate change experiment showed a large response across the temperate region dominated by the C₃ PFT, where LE increased by ∼ 30 % (H decreased) with Medlyn. Our study for the current day is in agreement with De Kauwe et al. (2015) who found a large impact of the Medlyn model on transpiration fluxes in needle leaved evergreen trees (∼ 30 % reduction) in the CABLE LSM. Coupled simulations using CABLE within the Australian Community Climate and Earth Systems Simulator (ACCESSv1.3b) showed that the Medlyn g_s scheme reduced the LE flux from the land surface over the boreal forests during JJA by 0.5–1.0 mm d⁻¹, leading to warmer daily maximum and minimum temperatures by up to 1.0 ^∘C and warmer extreme maximum temperatures by up to 1.5 ^∘C (Kala et al., 2015). In future simulations, this new parameterisation of the stomatal scheme in ACCESS1.3 substantially increased the intensity of future heatwaves across northern Eurasia (Kala et al., 2016).

5.3 Implications for land–atmosphere feedbacks

Modifying the leaf-level stomatal behaviour in JULES impacts the simulated surface energy fluxes. In our study, a change in stomatal opening results from either a direct change in the parameterisation of g_s or through altered stomatal behaviour in response to temperature. In our offline climate change simulation, thermal acclimation increased stomatal opening in response to long-term warming, and in some regions this increased the rate of transpiration and evaporative cooling and decreased the sensible heat flux. When coupled to an atmospheric model, such behaviours have potential to feedback on the land surface via changes in temperature, cloud cover and precipitation, as, for example, modelled by de Arellano et al. (2012), Kala et al. (2015, 2016), Kooperman et al. (2018) and Zeng et al. (2017). The extent and amplitude of acclimation-induced perturbations to surface energy fluxes in our offline simulation suggest a potential impact on regional-scale circulations, for example, across the East Asian monsoon region. The impact of these changes to the plant physiology routines in JULES on land–atmosphere feedbacks will be investigated in future work through coupled simulations in the HadGEM global climate model.

5.4 Limitations of this study

Across all latitudes, the changes introduced to JULES by the new plant physiology routines did not degrade the performance of JULES. All model configurations compared reasonably well to the FluxCom and MOD17 GPP products and FluxCom and GLEAM ET products, given that there are also uncertainties inherent in estimates from these products. For example, the satellite-based products of GPP have recently been shown to incorrectly capture the response of photosynthesis to CO₂, which means they potentially underestimate the response of GPP to rising atmospheric CO₂ (Keenan et al., 2021). Nevertheless, some notable biases in the model were identified that were common to all JULES model configurations, for example, the over-prediction of GPP and ET in the temperate and boreal region in SON and the over-prediction of both fluxes in MAM in the southern tropics (0 to −20^∘ S). Potential sources of error to consider may be the use of a prescribed climatology of MODIS-based LAI, which some studies have reported to be inaccurate over forested areas (Shabanov et al., 2005). Other processes currently missing in the model may also contribute to these large biases, such as a lack of seasonality in photosynthetic capacity (i.e. V_cmax and J_max), which has been demonstrated for many different forest species (Croft et al., 2017; Wilson et al., 2001) and without which likely causes overestimation of forest carbon exchange. For example, in SON the high GPP and ET bias occur in the northern temperate and boreal region, which could be linked to a lack of photosynthetic phenology in the model. Towards the end of the growing season, leaves in this region have reduced nitrogen content and therefore lower photosynthetic capacity, but because JULES uses a fixed value for photosynthetic capacity, JULES maintains a high rate of carbon assimilation despite having seasonal LAI.

More generally, this study revealed limited data to inform the temperature sensitivity response functions of different PFTs for implementation into LSMs. We found only a few datasets for C₃ grass/herbaceous vegetation (e.g. Wohlfahrt et al., 1999; Joseph et al., 2014), which represents only limited geographical coverage. Consequently, we fitted the temperature response function for this PFT in the Farquhar scheme to that of the existing function in the JULES Collatz photosynthesis scheme. We also encountered an issue regarding uncertainty about the temperature response functions at low temperatures. The data-led functions we implemented for all PFTs (with the exception of the C₃ PFT) from Kumarathunge et al. (2019b) showed higher rates of leaf-level photosynthesis at low leaf temperatures compared to the existing functions in the JULES Collatz scheme, where photosynthesis was much lower and goes to zero at 0 ^∘C for most PFTs (see PFT leaf-level temperature sensitivity curves for gross photosynthesis in Fig. S1). In our simulations, this led to higher GPP in DJF when using the Farquhar scheme, which increased biases with respect to FluxCom and MOD17 global estimates of GPP. It is desirable to use the temperature response functions from Kumarathunge et al. (2019b) as these are entirely data-led. However for some PFTs, the resulting behaviour of photosynthesis at very low temperatures looks potentially unrealistic, and the question here is how well constrained the temperature sensitivity curves are by observations at low temperatures. For global modelling applications, understanding the response of photosynthesis to temperature over a wide temperature range is essential, including at low temperatures as well as around the T_opt of photosynthesis for different species and PFTs. Additionally, increasing the understanding and data availability of the temperature sensitivity of different species from different biomes will allow for greater representation within LSMs of the variation that exists across the globe.

The simulations presented in this work use a prescribed map of vegetation cover which means the extent and location of each PFT does not change over time. The model can alternatively be run with dynamic vegetation enabled, which means the model predicts the extent of each PFT, and therefore vegetation cover can change in space and time as PFTs compete with each other in response to changing climatic conditions. Yet to be explored as part of this work is how changes to the plant physiology routines, as implemented here, might affect the extent of different PFTs over time when vegetation dynamics is enabled. For example, changes to the temperature response of photosynthesis may lead to a competitive advantage of one PFT over another, and therefore the vegetation distribution may be very different as temperatures rise compared to simulations that either use the original Collatz temperature sensitivities or do not include thermal acclimation of photosynthesis. We hypothesise, for example, that allowing for thermal acclimation of the temperature sensitivity of photosynthesis would make the vegetation distribution more stable in a warmer climate as vegetation can adjust its photosynthetic capacity to function more efficiently as temperatures rise. Applied in a coupled ESM, a change in vegetation distribution would impact projections of future climate change.

The treatment of soil moisture stress in JULES is through a linear response function (the β function; Eq. 12 in Best et al., 2011), the use of which in JULES and other LSMs has been identified as a key source of uncertainty (Blyth et al., 2011; Verhoef and Egea, 2014; Vidale et al., 2021). Incorrect representation of soil moisture stress has large impacts for modelled carbon and water fluxes and is of particular importance as droughts are predicted to increase in frequency or intensity in the future. Work is ongoing to improve the representation of soil moisture stress in JULES. Harper et al. (2021) investigated alternative parameterisations for β and found that increasing modelled soil depth and therefore plant access to deep soil moisture improved the simulation of soil moisture stress at eddy covariance flux tower sites. In addition, using soil matric potential instead of volumetric water content in the β function allowed for PFT-specific parameterisation of soil moisture stress responses to further improve modelled fluxes. Vidale et al. (2021) explored combinations of non-linear β function responses applied at different points in the photosynthesis–g_s pathway (i.e. carbon assimilation, g_s, or mesophyll conductance). They found that treatments allowing β to act on vegetation fluxes via stomatal and mesophyll routes were able to better capture the spatiotemporal variability in water-use efficiency during the growing season. However, in addition to these alternative parameterisations of β, further developments to how the soil–plant hydraulic system is represented in JULES are being made, including an optimality-based plant hydraulic transport model recently implemented in JULES (Eller et al., 2020).

Whilst the development of multi-layer canopy radiation models in LSMs has improved the simulation of radiation and energy within vegetation canopies, the interception of light by plants in JULES, like most LSMs, is not well represented despite being critical to predicting the uptake of carbon by plants (Loew et al., 2014). LSMs generally make the simplifying assumption that leaves are randomly arranged in space, instead of being clustered into tree crowns or around branches, leaving gaps in and around the canopy. Shortwave radiation is used by plants to photosynthesise, and canopy structure has a direct impact on the fraction of this radiation absorbed. Therefore canopy architecture plays an important role in the partitioning of incident solar radiation, photosynthesis, transpiration and momentum fluxes (Braghiere et al., 2019). More recently, alternative approaches are being considered to represent the forest light environment in LSMs to account for the structural effects of vegetation on radiation partitioning, ranging from canopy clumping parameterisations (Braghiere et al., 2019, 2020, 2021) to 3-dimensional models of the canopy light environment (Hogan et al., 2018; Kobayashi et al., 2012), embedded in radiative transfer schemes, although the latter tend to be computationally expensive (Yang et al., 2001). Braghiere et al. (2019) incorporated canopy clumping from satellite data into JULES, which resulted in an increase in carbon uptake by photosynthesis. The greatest effect was in the tropics, where the canopy clumping parameterisation allowed more light to reach the lower layers of the canopy, where photosynthesis tends to be limited by light availability.

5.5 Conclusions

Here we introduce new representations of plant physiological processes into the JULES model, building enhanced capability and allowing for stronger links between model and field studies. This work (a) introduces updated understanding of plant physiological processes into JULES, (b) increases the flexibility of the modelling capacity within JULES by allowing for the use of two alternative photosynthesis and g_s schemes, in addition to thermal acclimation of photosynthesis, and (c) provides new parameters that are entirely based on large observational datasets. Testing and evaluation at site level and globally show some key improvements are made to the JULES model. Thermal acclimation of photosynthesis coupled with the optimality-based g_s scheme led to improved simulated carbon fluxes across much of the tropics for the present day. With about 40 % of the world's vegetation carbon residing in tropical forests, they play a crucial role in regulating both regional and global climate through water and carbon cycle dynamics (Erb et al., 2018; Pan et al., 2011). Therefore, accurate representation of tropical carbon fluxes within LSMs is important. Thermal acclimation and the optimality-based g_s scheme also improved simulated carbon fluxes in the high-northern-latitude forests in the Northern Hemisphere summer, and the same model configuration also improved simulated water fluxes across much of this region in the same season. The optimality-based Medlyn g_s scheme reduced the LE flux substantially across the northern boreal forests in JJA. This change reflects a more conservative water-use strategy for the needleleaf evergreen tree PFT that dominates in this region as suggested by the global synthesis of experimental data from Lin et al. (2015). The current JULES Jacobs scheme parameterisation does not accurately capture the water-use strategy of this PFT. Our future climate experiment highlights the impact of thermal acclimation on simulating carbon cycle dynamics and energy fluxes in response to long-term warming. The potential impact of this altered stomatal behaviour on land–atmosphere feedbacks via changes in surface energy fluxes will be examined in future coupled simulations.