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 reﬁning the temperature sensitivity of photosynthesis, impacting modelled carbon, energy and water ﬂuxes. We test (1) an implementation of the Farquhar et al. (1980) photosynthesis scheme and associated plant functional type-dependent photosynthetic

Abstract.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 3 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 highnorthern-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

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 2 ) 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 2 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 2 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 2 emissions to atmospheric CO 2 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 3 and C 4 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 3 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 3 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(Franks et al., , 2018)), using evidence that stomata may behave to maximise CO 2 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 Ja-cobs (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., 2000Gunderson et al., , 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 multilayer 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) (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. (1991Collatz et al. ( , 1992) ) (for C 3 and C 4 plants) and Jacobs (1994) respectively.Photosynthesis and g s are modelled to respond to changes in environmental drivers of temperature, humidity deficit, light, CO 2 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(Atkin et al., , 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.

Farquhar photosynthesis for C 3 plants and parameterisation
We implement the Farquhar photosynthesis scheme (Farquhar et al., 1980) to describe the leaf-level biochemistry of photosynthesis for C 3 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) lightlimited 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 TPUlimited (triose phosphate utilisation) rate.Further, recent empirical studies suggest that TPU limitation rarely limits photosynthesis under present-day CO 2 concentrations and is also unlikely to limit photosynthesis at elevated CO 2 (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, where A p is the net potential (i.e.unstressed) leaf photosynthetic carbon uptake (mol m 2 s −1 ), R d is the rate of leaf respiration in the dark (mol m 2 s −1 ) and A n is the net photosynthetic rate (mol m 2 s −1 ), 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 2 s −1 ) is calculated as in Eq. ( 2).The maximum rate of carboxylation of Rubisco is determined by V cmax (mol m 2 s −1 ), c i and o a are the intercellular concentrations of CO 2 and O 2 (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 2 compensation point in the absence of mitochondrial respiration.
The light-limited rate of photosynthesis (A j , mol m 2 s −1 ) (Eq. 3) is a function of the rate of electron transport J (mol m 2 s −1 ) which is represented in Eq. ( 4).J depends on the incident photosynthetically active photon flux density Q (mol quanta m 2 s −1 ), the potential rate of electron transport J max (mol m 2 s −1 ), the apparent quantum yield of electron transport α (mol electrons mol −1 photon) fixed at 0.3 (mol electrons mol −1 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.
JULES currently uses Q 10 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: Here ), H a and H d (J mol −1 ) are the activation and deactivation energies respectively, and S (J mol −1 K −1 ) 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): To find new estimates for V cmax and the J max : V 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 2 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 2 response curves, which entrained data from 141 C 3 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 3 grasses; therefore to parameterise the temperature dependency of V cmax and J max for this PFT, we fitted both to the existing V cmax temperature response function in the Collatz scheme for C 3 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.

Medlyn model of g s and parameterisation
In JULES, g s (m s −1 ) is represented in Eq. (7).
where the factor 1.6 accounts for g s being the conductance for water vapour rather than CO 2 , R is the universal gas constant (J mol −1 K −1 ), T l is the leaf surface temperature (K), c a and c i (both Pa) are the leaf surface and internal CO 2 partial pressures respectively, and A n is the net photosynthetic rate.
Collatz Farquhar Collatz Farquhar Jacobs scheme as in Eq. ( 8) and relates the ratio of ambient (c a ) to leaf intercellular (c i ) partial pressure of CO 2 (c i /c a ), to leaf humidity deficit: where (Pa) is the CO 2 photorespiration compensation point, d q is the specific humidity deficit at the leaf surface (kg kg −1 ), and d q crit (kg kg −1 ) and f 0 are PFT-specific parameters representing the critical humidity deficit at the leaf surface and the leaf internal-to-atmospheric CO 2 ratio (c i /c 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 1 (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: PFT-specific values of the g 1 parameter were derived for the nine JULES PFTs from the global database of Lin et al. (2015) (Table 3).The g 1 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 1 , for each PFT, where individual species were assumed to be the random effect to account for the differences in the g 1 slope among species within the same group, following Lin et al. (2015).

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 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 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 opt v and T opt j ) increase by 0.44 and 0.33 • C per degree increase in T growth respectively, and the J max : V cmax ratio at 25 • C decreases by 0.035 • C per degree increase in T growth .
The x is either S v , S j or the J max : V 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 cmax ratio.T growth is the growth temperature (calculated online as the mean temperature of the previous 30 d).
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 2 ) 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:

Global-scale simulations
Four JULES simulations were performed globally for the period 1979-2013 as outlined in Table 5.These global presentday 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) https://doi.org/10.5194/gmd-15-5567-2022 Geosci.Model Dev., 15, 5567-5592, 2022  Kattge and Knorr (2007).Parameters describing the temperature sensitivity of photosynthesis ( S for V cmax and J max and the J max : V cmax ) allowed to acclimate to the temperature of the growth environment (T growth ).T opt adjusts to changes in T growth so varies spatially and temporally.
Yes n/a -not applicable.
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 2 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 −1 .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 physi-ology 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 2 concentration.Atmospheric CO 2 concentrations were prescribed based on observations up to 2014 as described in historical CMIP6 simulations (Eyring et al., 2016).From 2015 onwards, atmospheric CO 2 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): 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 = Fq.Jac and Y = Clz.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 = Fq.Med and Y = Fq.Jac.The impact of thermal acclimation is assessed where X = AcKK.Med and Y = Fq.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.

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 opt v , T opt j and the J max : V 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).At the C 3 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 3 , 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 3 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 3 grass PFT in the JULES model.In JJA, the adjustment of the temperature sensitivity of photosynthesis https://doi.org/10.5194/gmd-15-5567-2022 Geosci.Model Dev., 15, 5567-5592, 2022 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 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 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).

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 2 d −1 (up to 10 % reduction), whilst in the high northern latitudes, GPP is increased by up to 1.5 g C m 2 d −1 (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 2 , 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 2 d −1 (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 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 3 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).

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

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  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).
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 3 and C 4 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).

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 https://doi.org/10.5194/gmd-15-5567-2022Geosci.Model Dev., 15, 5567-5592, 2022 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 sitelevel 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).
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 3 grass PFT in JULES, which has a less conservative water-use strategy in the Medlyn scheme (high g 1 ) compared to Jacobs.This means in the Medlyn scheme, the C 3 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 opt v , T opt j and the J max : V 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).

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

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 opt v , T opt j and the J max : V 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 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 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 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 3 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 ac-climation 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.(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).

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. (2015Kala et al. ( , 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.

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 https://doi.org/10.5194/gmd-15-5567-2022Geosci.Model Dev., 15, 5567-5592, 2022 example, the satellite-based products of GPP have recently been shown to incorrectly capture the response of photosynthesis to CO 2 , which means they potentially underestimate the response of GPP to rising atmospheric CO 2 (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 overprediction 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 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 3 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 3 PFT) from Kumarathunge et al. (2019b) showed higher rates of leaflevel 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(Braghiere et al., , 2020(Braghiere et al., , 2021) ) to 3dimensional 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.

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.
Author contributions.RJO performed simulations and analysis and wrote the first version of the manuscript.DBC, LMM and RJO developed the model.PLV, PCM and MT provided data for the future climate runs, help with developing the JULES suites and general expertise.CH assisted with analysis.SF and VSS provided ancillary data for forcing the model.LMM, CMT, CH, PLV, BEM, PCM and MT contributed to editing the manuscript.All authors contributed to discussions throughout to develop the work.
Competing interests.The contact author has declared that none of the authors has any competing interests.
Disclaimer.Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figure 1 .
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 3 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).

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

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

Figure 4 .
Figure 4. Mean (2002 to 2012) GPP (g C m 2 d −1 ) 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 .
Figure 5. Mean (2002 to 2012) evapotranspiration (ET mm d −1 ) 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 TableS6(FluxCom) and TableS7(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 toJung et al. (2019).

Figure 6 .
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 2 d −1 ) products or (c) FluxCom and (d) GLEAM global ET (mm d −1 ) products for JJA over the period 2002 to 2012.Actual RMSE values shown in Figs.S12 and S13.

Figure 7 .
Figure7.The difference of difference approach (Eq.12) to determine the impact on GPP (g C m 2 d −1 ), LE and H (both W m 2 ) 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.S16AcKK.Med-Fq.Med, the effect of the Medlyn g s scheme is calculated from Fig.S16Fq.Med-Fq.Jac and the effect of the photosynthesis scheme is calculated from Fig.S16Fq.Jac-Clz.Jac.

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

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

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

Table 5 .
Description of the four model experiments performed both at site level and globally, with the JULES land surface model.
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 3 grass PFT and needleleaf evergreen tree PFT in JULES.This reflects a change to the water-use strategy of these PFTs as reported byLin et al. (2015)that is not currently captured by param-eterisations 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 3 PFT, where LE increased by ∼ 30 % (H decreased) with Medlyn.Our study for the current day is in agreement with DeKauwe 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.0mm d −1 , 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