LPJ-GUESS/LSMv1.0: a next-generation land surface model with high ecological realism

Land biosphere processes are of central importance to the climate system. Specifically, ecosystems interact with the atmosphere through a variety of feedback loops that modulate energy, water, and CO₂ fluxes between the land surface and the atmosphere across a wide range of temporal and spatial scales. Human land use and land cover modification add a further level of complexity to land–atmosphere interactions. Dynamic global vegetation models (DGVMs) attempt to capture land ecosystem processes and are increasingly incorporated into Earth system models (ESMs), which makes it possible to study the coupled dynamics of the land biosphere and the climate. In this work we describe a number of modifications to the LPJ-GUESS DGVM, aimed at enabling direct integration into an ESM. These include energy balance closure, the introduction of a sub-daily time step, a new radiative transfer scheme, and improved soil physics. The implemented modifications allow the model (LPJ-GUESS/LSM) to simulate the diurnal exchange of energy, water, and CO₂ between the land ecosystem and the atmosphere and thus provide surface boundary conditions to an atmospheric model over land. A site-based evaluation against FLUXNET2015 data shows reasonable agreement between observed and modelled sensible and latent heat fluxes. Differences in predicted ecosystem function between standard LPJ-GUESS and LPJ-GUESS/LSM vary across land cover types. We find that the emerging ecosystem composition and carbon fluxes are sensitive to both the choice of stomatal conductance model and the response of plant water uptake to soil moisture. The new implementation described in this work lays the foundation for using the well-established LPJ-GUESS DGVM as an alternative land surface model (LSM) in coupled land–biosphere–climate studies, where an accurate representation of ecosystem processes is essential.

1 Introduction

The land surface is of central importance in the climate system, as feedbacks between the land biosphere and the atmosphere impact climate across a wide range of temporal and spatial scales (Pitman, 2003). Biological processes affected by climate variations can feed back into the climate by modulating the fluxes of energy and water between vegetation and the atmosphere (Guo et al., 2006; Green et al., 2017). For example, the early, strong greening caused by the warming climate can enhance evapotranspiration, which may result in a seasonal cooling effect or in an amplification of heat waves, depending on regional characteristics and water availability (Peñuelas et al., 2009; Lorenz et al., 2013). On decadal timescales, decreased vegetation cover caused by reduced rainfall can further decrease local precipitation (Zeng et al., 1999). Large-scale shifts in vegetation cover in response to climate change can affect global and regional climate by altering the radiation and water budgets (O'ishi and Abe-Ouchi, 2009; Levis et al., 2000; Wramneby et al., 2010; Wu et al., 2021).

The climate and the biosphere are also coupled biogeochemically through the carbon cycle (Luo, 2007). Increased atmospheric carbon dioxide (CO₂) concentration promotes vegetation growth through CO₂ fertilization, which increases plant CO₂ absorption from the atmosphere. However, higher temperatures caused by a higher atmospheric CO₂ concentration enhance the release of CO₂ from respiration (Cramer et al., 2001; Piao et al., 2013). Other important effects relate to extreme events (Zscheischler et al., 2014), disturbances (Kurz et al., 2008; Metsaranta et al., 2010), or interaction with the nitrogen cycle (Arneth et al., 2010; Lamarque et al., 2013; Ciais et al., 2014).

Of particular importance is the added complexity arising from land use and land cover change. Conversion of forests into cropland or grassland increases surface albedo, which may promote surface cooling at temperate latitudes (e.g. de Noblet-Ducoudré et al., 2012) but is also a significant contributor to anthropogenic CO₂ emissions (Arneth et al., 2017; Le Quéré et al., 2018). Observations and model studies suggest that historical land cover changes over the industrial era have had a minor net impact on the climate system at the global scale, but regional effects are large (Brovkin et al., 2004; Pongratz et al., 2010; Christidis et al., 2013). Further complexity arises from the interaction between land use change and the water cycle (e.g. Narisma and Pitman, 2003; Kumar et al., 2013; Lawrence and Vandecar, 2015), atmospheric circulation (Swann et al., 2012; Wu et al., 2017), and atmospheric teleconnections (Werth and Avissar, 2002; Medvigy et al., 2013).

Incorporating DGVMs into ESMs allows the interactions between the biosphere and the rest of the climate system to be studied on the long timescales of vegetation dynamics and biogeochemical and biogeographical responses (Quillet et al., 2010; Fisher et al., 2018). There is considerable uncertainty regarding the carbon cycle response to future climate warming scenarios (Friedlingstein et al., 2006, 2014; Jones et al., 2013), which has been attributed to uncertainty in the representation of land surface processes (Huntingford et al., 2009; Booth et al., 2012; Friend et al., 2014) and differences between the global circulation models (GCMs) used to make such projections (Ahlström et al., 2013, 2017; Schurgers et al., 2018). Improved representations of land biosphere processes and land use change in ESMs are therefore essential to constrain climate change projections (Friend et al., 2014) and thus to support the assessment of mitigation and adaptation strategies.

DGVMs are frequently integrated into ESMs through an intermediary land surface model (LSM), which facilitates the sub-daily energy, water, and gas exchange calculations (e.g. Bonan et al., 2003; Krinner et al., 2005; Smith et al., 2011; Döscher et al., 2022). This is necessary because DGVMs normally run on a daily or longer time step, while atmospheric models may use time steps ranging from seconds to tens of minutes, depending on the required resolution. This indirect approach can, however, entail inconsistencies between the DGVM and the LSM, such as the use of different time steps and temperatures in photosynthetic calculations, duplicated or inconsistent soil water tracking, inconsistent carbon mass balance, or different characterization of vegetation types. In this work we modify the LPJ-GUESS DGVM (Smith et al., 2001, 2014) to enable coupling with an atmospheric model without the need for a mediating LSM. LPJ-GUESS simulates a wide range of land biosphere processes, including vegetation growth, establishment, and mortality, plant functional type (PFT) competition, disturbances, wildfires, and land use change. This model has been used in a broad range of applications, including coupled biosphere–atmosphere regional (Wramneby et al., 2010; Smith et al., 2011; Zhang et al., 2014, 2018; Wu et al., 2016, 2021) and global (Weiss et al., 2014; Alessandri et al., 2017; Forrest et al., 2020; Döscher et al., 2022) studies, although these suffer from the above-mentioned limitations of the indirect coupling approach. LPJ-GUESS is maintained by an international developer community and undergoes active development and evaluation, which makes it a suitable choice for studying climate–biosphere interactions.

Coupling LPJ-GUESS with an atmospheric model requires it to be able to calculate diurnal energy and water exchange rates between plant canopies and the atmosphere. To achieve this, we introduced several major modifications to LPJ-GUESS v4.0, namely (a) energy balance closure on a sub-daily time step and (b) a new radiative transfer scheme, capable of calculating upwelling short-wave radiation dynamically on a sub-daily time step as well as accounting for direct and diffuse solar radiation separately, and (c) an improved representation of heat and water transport in the soil. Section 2 describes these modifications in detail. A site-based comparison with the standard LPJ-GUESS model and an evaluation of the modelled fluxes against eddy covariance data are presented in Sect. 3. Finally, the work is discussed and summarized in Sect. 4.

2 Model description

2.1 LPJ-GUESS

LPJ-GUESS (Smith et al., 2001, 2014) is a process-based model of vegetation dynamics and ecosystem biogeochemistry and water cycling that incorporates tree demographic processes and competition for light, space, and soil resources among co-occurring PFTs. Capturing the establishment, growth, and death of individuals allows us to better represent the mechanisms underlying competition, population and community structural dynamics, carbon assimilation, and ecosystem carbon turnover (Smith et al., 2001; A. Wolf et al., 2011). In LPJ-GUESS, natural vegetation is represented as a co-occurring mixture of different PFTs (Table 1), divided into age classes or cohorts, in a modelled area or patch. New cohorts can establish in the patch when climatic conditions are within PFT-prescribed bioclimatic limits and compete with other cohorts for light, water, and soil nitrogen. Each cohort assimilates atmospheric CO₂ at a rate, updated daily in the standard model, that depends on the amount of photosynthetically active radiation (PAR) it absorbs, water availability, temperature, and the maximum rate of carboxylation, V_max. The maximum rate of carboxylation is estimated under the assumption that plants redistribute leaf nitrogen content across the canopy so as to maximize net assimilation at the canopy level (Haxeltine and Prentice, 1996) and is limited by nitrogen availability (Smith et al., 2014). The yearly assimilated carbon is distributed between roots, leaves, and, in the case of woody PFTs, sapwood, according to a set of PFT-specific allometric constraints. The phenological status of the cohorts (for summergreen and raingreen PFTs) is updated daily. Population dynamics (establishment and mortality) and non-fire-related disturbances are modelled as stochastic processes, influenced by environmental factors, vegetation structure, growth, and competition. Disturbances occur recurrently and destroy all vegetation in a patch, restarting the successional cycle. Wildfires are modelled explicitly (Thonicke et al., 2001). At any given geographical location (grid cell), a number of replicate patches with independent successional histories are simulated.

LPJ-GUESS can represent managed land (croplands, pastures/rangelands, and managed forest) and land use change (Lindeskog et al., 2013, 2021; Olin et al., 2015). Each grid cell contains different land cover types or stands, which are updated every simulation year (for example, to simulate conversion of forest to cropland). Croplands are represented as single PFT stands, distinguishing various rainfed and irrigated crop functional types. In pasture stands only grassy PFTs are allowed to establish. Simulated land management practices include crop sowing, irrigation, fertilization, harvest, rotation and abandonment, and pasture grazing.

2.2 Model modifications

Figure 1 shows a comparison of the daily loop in standard LPJ-GUESS and in the new LSM implementation. In both versions, phenology and soil organic matter dynamics are calculated daily, and carbon allocation (growth) and vegetation dynamics (establishment, mortality, and disturbance) are computed at the end of every simulation year.

Figure 1 Flowchart of the main daily simulation loop in standard LPJ-GUESS (red branch) and the modified version (LPJ-GUESS/LSM, blue branch). The shaded area indicates the sub-daily loop in the modified version. The dashed line encloses coupled iterative calculations.

Radiative transfer in standard LPJ-GUESS is based on Beer's law (Monsi and Saeki, 1953, 2005). The canopy is divided into vertical layers, each absorbing a fraction of the PAR let through by the layer above. The PAR absorbed by each layer is then split among cohorts according to their share of the leaf area index (LAI) in that layer. In this way, taller cohorts have access to more PAR and shade than the lower layers of the canopy. Daily unstressed values of V_max and canopy conductance g_pot are first computed for each cohort assuming well-watered conditions. The actual evapotranspiration rate in the patch is then calculated as the minimum of a potential rate, determined by atmospheric conditions and g_pot, and a supply rate, which depends on the amount of soil water available for uptake and the vegetation rooting profiles. For each cohort, the model calculates a daily assimilation rate that is consistent with its contribution to the total patch evapotranspiration. The soil column consists of a top layer of 0.5 m and a bottom layer of 1 m thickness. The fraction of root matter in each soil layer is PFT-specific. Soil water content is updated taking into account daily precipitation, interception, percolation between the two layers, evapotranspiration, and runoff. Daily soil temperature is calculated as a dampened, lagged oscillation around the annual mean of the forcing air temperature, as described in Sitch et al. (2003). More detailed descriptions of the radiative transfer, evapotranspiration, assimilation, and soil organic matter calculations can be found in the Supplement to Smith et al. (2001), Smith et al. (2014), and references therein. The hydrology scheme is described in Gerten et al. (2004).

In the LSM implementation, radiative transfer, energy balance, assimilation, and soil heat and water transport are all solved on a sub-daily basis. Elaborating on Dai et al. (2004), each cohort is conceptualized as two big leaves representing its sunlit and shaded parts. Sunlit leaves receive direct solar radiation and diffuse radiation, while shaded leaves receive only diffuse radiation. The total LAI for each cohort is calculated dynamically by LPJ-GUESS. A stem area index (SAI) was added to account for the impact of stems and branches in the energy balance and radiative transfer calculations. Whole-canopy leaf area and plant area (PAI) indices are obtained by aggregating over cohorts denoted by i:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\text{LAI}}_{\mathrm{c}}& {\displaystyle }=\sum _i{\text{LAI}}^{\left(i\right)};\text{(2)}& {\displaystyle }{\text{PAI}}_{\mathrm{c}}& {\displaystyle }=\sum _i[{\text{LAI}}^{\left(i\right)}+{\text{SAI}}^{\left(i\right)}].\end{array}$$

Based on Kucharik et al. (1998), we set the stem area index of woody PFTs to 10 % of their leaf area index at full leaf coverage. Grasses do not have a stem area index. The sunlit and shaded fractions of leaf and plant area indices are updated in the radiative transfer routine on a sub-daily basis (Sect. 2.2.2).

We replaced the original two-layer soil column with a new profile consisting of nine layers. The top four layers have thicknesses of 7, 10, 13, and 20 cm, in order of increasing depth, and correspond to the top soil layer in the original soil column. The next three layers have thicknesses of 30, 30, and 40 cm and correspond to the original bottom layer. These seven layers constitute the rooting zone. The new water transport scheme assumes, for simplicity, free gravitational drainage at the bottom of the soil column, which can lead to excessive soil dryness during dry periods. Additionally, no heat flux is allowed through the bottom boundary, an approximation better met at higher soil depths. In order to mitigate spurious effects derived from this choice of boundary conditions, we extended the soil column with two additional layers of 50 and 100 cm, reaching a total depth of 3 m.

The sunlit and shaded leaves of each cohort have different assimilation rates and stomatal conductances. The temperatures of sunlit and shaded leaves are different but common to all the cohorts in the patch. The vertical layering of the canopy is kept in the radiation calculations, but the new scheme distinguishes direct and diffuse radiation and two separate wavebands (visible and near infrared). Infrared radiation does not contribute to photosynthetic assimilation but needs to be accounted for in the energy balance calculations. A separate treatment of diffuse and direct radiation allows us to resolve sunlit and shaded leaves. This approach has been shown to lead to predictions of fluxes of energy, water, and CO₂ that are comparable in accuracy with those made by more complex, and considerably more computationally expensive, multi-layered canopy models (Wang and Leuning, 1998).

Each cohort exchanges sensible and latent heat with a common canopy air space, which in turn exchanges sensible and latent heat with the atmosphere (Fig. 2). Assimilation and evapotranspiration are calculated consistently in the energy balance routine. Daily averages of absorbed PAR are used to update V_max for each cohort. The new energy balance, radiative transfer, and soil physics calculations are detailed in Sects. 2.2.1 to 2.2.5.

Figure 2 Networks of sensible (red) and latent (blue) heat exchange between the ground surface, the canopy, and the atmosphere in the patch. Light green indicates the sunlit fraction of the cohorts, dark green the shaded fraction. g_aa is the aerodynamic conductance from the canopy air to the atmospheric reference level (z_atm). z₀ and z_d are, respectively, roughness length and zero plane displacement. g_ab is the aerodynamic conductance for heat/moisture flux from the ground surface to the canopy air. g_surf is the surface conductance for moisture flux. $g_{\text{h,sun[sha]}}^{\left(i\right)}$ is the conductance for sensible heat transport from the sunlit (shaded) part of cohort i to the canopy air. $g_{\text{w,sun[sha]}}^{\left(i\right)}$ is the conductance for moisture transport from the sunlit (shaded) part of cohort i to the canopy air, with the contributions from stomatal conductance and leaf boundary layer conductance represented explicitly. A dry canopy (f_wet=0) is represented for clarity.

2.2.1 Energy balance

The energy balance of the patch canopy is described by the following equations (e.g. Bonan, 2015):

$$\begin{array}{}\text{(3)}& {\displaystyle }{S}_{\text{sun}}& {\displaystyle }=L_{\text{sun}}+H_{\text{sun}}+\mathit{\lambda }{E}_{\text{sun}},\text{(4)}& {\displaystyle }{S}_{\text{sha}}& {\displaystyle }=L_{\text{sha}}+H_{\text{sha}}+\mathit{\lambda }{E}_{\text{sha}},\end{array}$$

where the S terms are absorbed short-wave radiation, L is net emitted long-wave radiation, H is sensible heat flux towards the canopy air space, E is water vapour flux towards the canopy air space, and λ is latent heat of vaporization (here taken to be constant; $\mathit{\lambda }=\mathrm{2.44}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹). The sub-indices “sun” and “sha” refer to the sunlit and shaded parts of the canopy. The calculation of the short-wave and long-wave radiation terms is detailed in Sects. 2.2.2 and 2.2.3.

The sensible heat flux from the sunlit part of the canopy to the canopy air space is formulated as

$$\begin{array}{}\text{(5)}& H_{\text{sun}}=-\mathrm{2}{\text{PAI}}_{\text{c,sun}}\mathit{\rho }{c}_P{g}_{\mathrm{b}}(T_{\text{ca}}-T_{\text{sun}}),\end{array}$$

where PAI_c,sun is the plant area index of the sunlit canopy, ρ is air density, c_P is the specific heat of air at constant pressure, g_b is average leaf boundary layer conductance (e.g. Bonan, 2015), T_ca is the temperature of the canopy air, and T_sun is the temperature of the sunlit canopy. The factor 2 expresses heat loss from both sides of the leaf and stem elements.

The latent heat flux from the sunlit part of the canopy to the canopy air is

$$\begin{array}{}\text{(6)}& \mathit{\lambda }{E}_{\text{sun}}=-\mathit{\rho }\mathit{\lambda }{g}_{\text{w,sun}}[q_{\text{ca}}-q^{\ast }(T_{\text{sun}}\left)\right],\end{array}$$

where q_ca is the specific humidity of the canopy air, q^∗(T_sun) is the specific humidity inside the stomatal cavity, taken to be the saturated humidity at the leaf temperature, and g_w,sun is the conductance for water vapour flux from the sunlit part of the canopy to the canopy air space. The latter is calculated as a weighted average of the contributions from evaporation of intercepted water and transpiration through the stomata (Appendix A):

$$\begin{array}{ll}{\displaystyle }{g}_{\text{w,sun}}& {\displaystyle }=f_{\text{wet}}{\mathit{\eta }}_{\text{sun}}{\text{PAI}}_{\text{c,sun}}{g}_{\mathrm{b}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}-f_{\text{wet}}{\mathit{\eta }}_{\text{sun}})\sum _i{\text{LAI}}_{\text{sun}}^{\left(i\right)}{\displaystyle \frac{g_{\text{s,sun}}^{\left(i\right)}{g}_{\mathrm{b}}}{g_{\text{s,sun}}^{\left(i\right)}+g_{\mathrm{b}}}}.\end{array}$$

In this equation f_wet is the wet fraction of the canopy, the factor η_sun limits evaporation to the amount of intercepted water present in the canopy, and ${\text{LAI}}_{\text{sun}}^{\left(i\right)}$ is the leaf area index of the sunlit part of cohort i. The stomatal conductance of the sunlit leaves of cohort i, $g_{\text{s,sun}}^{\left(i\right)}$, is related to the leaf-level net photosynthetic rate through a semi-empirical model. We implemented two selectable stomatal conductance models: the Ball–Berry model (Ball et al., 1987) and the Medlyn model (Medlyn et al., 2011). Equations analogous to Eqs. (5) through (7) apply to the shaded part of the canopy.

The energy balance equation for the ground surface is

$$\begin{array}{}\text{(8)}& S_{\mathrm{g}}=L_{\mathrm{g}}+H_{\mathrm{g}}+\mathit{\lambda }{E}_{\mathrm{g}}+G,\end{array}$$

where G is heat conducted into the ground. The sensible heat from the ground surface to the canopy air space is

$$\begin{array}{}\text{(9)}& H_{\mathrm{g}}=-\mathit{\rho }{c}_P{g}_{\text{ab}}(T_{\text{ca}}-T_{\mathrm{g}}),\end{array}$$

where g_ab is the aerodynamic conductance from the ground surface to the canopy air space, which is calculated following Sakaguchi and Zeng (2009). The latent heat from the ground surface to the canopy air is given by

$$\begin{array}{}\text{(10)}& \mathit{\lambda }{E}_{\mathrm{g}}=-\mathit{\rho }\mathit{\lambda }{\displaystyle \frac{g_{\text{surf}}{g}_{\text{ab}}}{g_{\text{surf}}+g_{\text{ab}}}}[q_{\text{ca}}-\mathit{\alpha }{q}^{\ast }(T_{\mathrm{g}}\left)\right],\end{array}$$

where we used the model of Sakaguchi and Zeng (2009) for the surface conductance g_surf, and αq^∗(T_g) is the specific humidity of the air at the ground surface (Philip, 1957).

The heat conducted into the ground is calculated as

$$\begin{array}{}\text{(11)}& G=-{\mathit{\kappa }}_{\mathrm{s}}^{\text{(1)}}{\displaystyle \frac{T_{\mathrm{s}}^{\text{(1)}}-T_{\mathrm{g}}}{\mathrm{\Delta }{z}^{\text{(1)}}/\mathrm{2}}},\end{array}$$

where ${\mathit{\kappa }}_{\mathrm{s}}^{\text{(1)}}$, $T_{\mathrm{s}}^{\text{(1)}}$, and Δz⁽¹⁾ are, respectively, the thermal conductivity, the temperature, and the thickness of the top soil layer.

The following two equations express conservation of latent and sensible heat:

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{H}^{\uparrow }=H_{\text{sun}}+H_{\text{sha}}+H_{\mathrm{g}},\text{(13)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }{E}^{\uparrow }=\mathit{\lambda }{E}_{\text{sun}}+\mathit{\lambda }{E}_{\text{sha}}+\mathit{\lambda }{E}_{\mathrm{g}},\end{array}$$

where H^↑ and λE^↑ are, respectively, the sensible and latent heat fluxes into the atmosphere, given by

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{H}^{\uparrow }=-\mathit{\rho }{c}_P{g}_{\text{aa}}(T_{\text{atm}}-T_{\text{ca}}),\text{(15)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }{E}^{\uparrow }=-\mathit{\rho }\mathit{\lambda }{g}_{\text{aa}}(q_{\text{atm}}-q_{\text{ca}}).\end{array}$$

Here, T_atm and q_atm are the temperature and specific humidity of the air at the atmospheric reference level, and g_aa is the aerodynamic conductance above the canopy. The latter is calculated by applying the Monin–Obukov similarity theory (see e.g. Bonan, 2015), which requires knowledge of the surface roughness length, z₀, and the zero plane displacement, z_d. These are calculated as a function of the canopy plant area index, PAI_c, and the canopy height, h_c, according to the model of Raupach (1994, 1995):

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{z_{\mathrm{d}}}{h_{\mathrm{c}}}}& {\displaystyle }=\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\exp (-\sqrt{\mathrm{7.5}{\text{PAI}}_{\mathrm{c}}})}{\sqrt{\mathrm{7.5}{\text{PAI}}_{\mathrm{c}}}}},\text{(17)}& {\displaystyle \frac{z_{\mathrm{0}}}{h_{\mathrm{c}}}}& {\displaystyle }=\left(\mathrm{1}-{\displaystyle \frac{z_{\mathrm{d}}}{h_{\mathrm{c}}}}\right)\exp \left(-{\displaystyle \frac{k}{\mathit{\beta }}}+\mathrm{0.193}\right),\end{array}$$

where k=0.4 is the von Karman constant and $\mathit{\beta }=min(\sqrt{\mathrm{0.003}+\mathrm{0.15}{\text{PAI}}_{\mathrm{c}}},\mathrm{0.3})$. Canopy height is calculated, following Forrest et al. (2020), as an average of cohort heights weighted by their foliar projective cover (FPC).

Equations (3), (4), and (8), subject to constraints (12) and (13), are solved simultaneously every time step with a multidimensional Newton–Rhapson method (e.g. Press, 2003).

2.2.2 Short-wave radiative transfer

We adapted the two-big-leaf model of Dai et al. (2004), based on the two-stream model of Dickinson (1983) and Sellers (1985), to LPJ-GUESS's multiple-cohort, vertically layered canopy. This approach considers direct solar radiation and diffuse atmospheric radiation separately. The intensity of the direct solar radiation beam in the canopy decreases exponentially with cumulative plant area index P (measured from the top of the canopy, increasing downwards) (Monsi and Saeki, 1953, 2005):

$$\begin{array}{}\text{(18)}& I_{\mathrm{D}}^{\downarrow }\left(P\right)=I_{\text{D0}}^{\downarrow }{e}^{-kP},\end{array}$$

where $I_{\text{D0}}^{\downarrow }$ is incoming direct solar radiation and k is the direct-beam extinction coefficient. The profile of diffuse radiation in the canopy results from the multiple scattering and backscattering of incoming radiation by leaves and stems. Corrected profiles (normalized by incoming radiation) of scattered direct-beam (${\widehat{I}}_{\mathrm{b}}^{\uparrow }$ and ${\widehat{I}}_{\mathrm{b}}^{\downarrow }$) and scattered atmospheric diffuse radiation (${\widehat{I}}_{\mathrm{a}}^{\uparrow }$ and ${\widehat{I}}_{\mathrm{a}}^{\downarrow }$) are given in analytic form in Dai et al. (2004) (the arrows indicate the direction of propagation).

The direct-beam radiation absorbed in a canopy layer l between P and P+Δ_lP is calculated as the fraction of the decrease in direct-beam intensity in that layer that is not scattered:

$$\begin{array}{}\text{(19)}& S_{\mathrm{D}}^{\left(l\right)}=-(\mathrm{1}-\mathit{\omega }){\mathrm{\Delta }}_l{I}_{\mathrm{D}}^{\downarrow },\end{array}$$

where ω is the direct-beam scattering coefficient, and Δ_l denotes change across layer l. The diffuse radiation absorbed in the layer is the sum of the radiation from the direct beam that is scattered and reabsorbed in the layer and the contribution from the diffuse beams:

$$\begin{array}{ll}{\displaystyle }{S}_{\mathrm{d}}^{\left(l\right)}& {\displaystyle }=-\mathit{\omega }{\mathrm{\Delta }}_l{I}_{\mathrm{D}}^{\downarrow }+I_{\text{D0}}^{\downarrow }({\mathrm{\Delta }}_l{\widehat{I}}_{\mathrm{b}}^{\uparrow }-{\mathrm{\Delta }}_l{\widehat{I}}_{\mathrm{b}}^{\downarrow })\\ \text{(20)}& {\displaystyle }& {\displaystyle }+I_{\text{d0}}^{\downarrow }({\mathrm{\Delta }}_l{\widehat{I}}_{\mathrm{a}}^{\uparrow }-{\mathrm{\Delta }}_l{\widehat{I}}_{\mathrm{a}}^{\downarrow }),\end{array}$$

where $I_{\text{d0}}^{\downarrow }$ is incoming atmospheric diffuse radiation. The radiation absorbed by the sunlit and shaded parts of this layer is

$$\begin{array}{}\text{(21)}& {\displaystyle }{S}_{\text{sun}}^{\left(l\right)}& {\displaystyle }=S_{\mathrm{D}}^{\left(l\right)}+f_{\text{sun}}^{\left(l\right)}{S}_{\mathrm{d}}^{\left(l\right)},\text{(22)}& {\displaystyle }{S}_{\text{sha}}^{\left(l\right)}& {\displaystyle }=f_{\text{sha}}^{\left(l\right)}{S}_{\mathrm{d}}^{\left(l\right)},\end{array}$$

where the sunlit and shaded fractions of the layer are given by

$$\begin{array}{}\text{(23)}& {\displaystyle }{f}_{\text{sun}}^{\left(l\right)}& {\displaystyle }=-{\displaystyle \frac{e^{-k(P+{\mathrm{\Delta }}_{l}P)}-e^{-kP}}{k{\mathrm{\Delta }}_{l}P}},\text{(24)}& {\displaystyle }{f}_{\text{sha}}^{\left(l\right)}& {\displaystyle }=\mathrm{1}-f_{\text{sun}}^{\left(l\right)}.\end{array}$$

The total amount of short-wave radiation absorbed by the sunlit and shaded parts of the canopy is obtained by summing over layers:

$$\begin{array}{}\text{(25)}& {\displaystyle }{S}_{\text{sun}}& {\displaystyle }=\sum _l{S}_{\text{sun}}^{\left(l\right)},\text{(26)}& {\displaystyle }{S}_{\text{sha}}& {\displaystyle }=\sum _l{S}_{\text{sha}}^{\left(l\right)}.\end{array}$$

The short-wave radiation absorbed by the ground surface is calculated as the difference between the downward and upward beams at P=PAI_c:

$$\begin{array}{ll}{\displaystyle }{S}_{\mathrm{g}}=I_{\mathrm{D}}^{\downarrow }\left({\text{PAI}}_{\mathrm{c}}\right)& {\displaystyle }+I_{\text{D0}}^{\downarrow }\left[{\widehat{I}}_{\mathrm{b}}^{\downarrow }\right({\text{PAI}}_{\mathrm{c}})-{\widehat{I}}_{\mathrm{b}}^{\uparrow }({\text{PAI}}_{\mathrm{c}}\left)\right]\\ \text{(27)}& {\displaystyle }& {\displaystyle }+I_{\text{d0}}^{\downarrow }\left[{\widehat{I}}_{\mathrm{a}}^{\downarrow }\right({\text{PAI}}_{\mathrm{c}})-{\widehat{I}}_{\mathrm{a}}^{\uparrow }({\text{PAI}}_{\mathrm{c}}\left)\right].\end{array}$$

The short-wave radiation reflected back at the atmosphere is obtained by evaluating the upward beams at P=0:

$$\begin{array}{}\text{(28)}& I^{\uparrow }=I_{\mathrm{b}}^{\uparrow }\left(\mathrm{0}\right)+I_{\mathrm{a}}^{\uparrow }\left(\mathrm{0}\right).\end{array}$$

The optical elements in the canopy have different properties in the visible and near-infrared wavebands, so the equations above are applied separately to these two parts of the spectrum, and the contributions are summed to calculate total absorption. In order to keep the model development process tractable, we set the optical properties of the canopy to the following values, regardless of PFT:

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\text{leaf, vis}}=\mathrm{0.1};{\mathit{\alpha }}_{\text{stem, vis}}=\mathrm{0.16},\text{(30)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\text{leaf, vis}}=\mathrm{0.05};{\mathit{\tau }}_{\text{stem, vis}}=\mathrm{0.001},\text{(31)}& {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\text{leaf, nir}}=\mathrm{0.45};{\mathit{\alpha }}_{\text{stem, nir}}=\mathrm{0.39},\text{(32)}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\text{leaf, nir}}=\mathrm{0.25};{\mathit{\tau }}_{\text{stem, nir}}=\mathrm{0.001},\end{array}$$

where α is absorptivity, τ is transmissivity, “vis” refers to visible radiation, and “nir” refers to near-infrared. These values were taken from the ones assigned to tropical trees by Oleson et al. (2010). Soil albedo is calculated from the soil dry- and moisture-saturated reflectances and the water content of the top soil layer following Oleson et al. (2010). Soil colour classes are from Lawrence and Chase (2007) and were obtained from the dataset included in the CLM4.0 code (Lawrence et al., 2011).

The PAR absorbed by the sunlit leaves of a cohort i is obtained as the sum over layers of the absorbed visible radiation weighted by the cohort's fractional leaf area index in each layer:

$$\begin{array}{}\text{(33)}& {\displaystyle }{\text{PAR}}_{\text{sun}}^{\left(i\right)}=\sum _l{S}_{\text{sun,vis}}^{\left(l\right)}{\displaystyle \frac{{\text{LAI}}^{(i,l)}}{{\text{PAI}}^{\left(l\right)}}}.\end{array}$$

The sunlit leaf and plant area indices of cohort i are obtained by aggregating over layers:

$$\begin{array}{}\text{(34)}& {\displaystyle }{\text{LAI}}_{\text{sun}}^{\left(i\right)}& {\displaystyle }=\sum _l{f}_{\text{sun}}^{\left(l\right)}{\text{LAI}}^{(i,l)},\text{(35)}& {\displaystyle }{\text{PAI}}_{\text{sun}}^{\left(i\right)}& {\displaystyle }=\sum _l{f}_{\text{sun}}^{\left(l\right)}[{\text{LAI}}^{(i,l)}+{\text{SAI}}^{(i,l)}].\end{array}$$

The sunlit plant area index for the whole canopy is calculated by summing over cohorts:

$$\begin{array}{}\text{(36)}& {\displaystyle }{\text{PAI}}_{\text{sun, c}}=\sum _i{\text{PAI}}_{\text{sun}}^{\left(i\right)}.\end{array}$$

Equations analogous to Eqs. (33) to (36) apply to the shaded parts of the canopy.

2.2.3 Long-wave radiative transfer

The long-wave radiation emitted by the sunlit part of the canopy is (Dai et al., 2004)

$$\begin{array}{}\text{(37)}& L_{\text{sun}}={\mathit{\gamma }}_{\text{sun}}(\mathrm{2}\mathit{\sigma }{T}_{\text{sun}}^{\mathrm{4}}-L^{\downarrow }-\mathit{\sigma }{T}_{\mathrm{g}}^{\mathrm{4}}),\end{array}$$

where σ is the Stefan–Boltzmann constant, L^↓ is the incoming atmospheric long-wave radiation, T_sun, T_g is expressed in Kelvin, and

$$\begin{array}{}\text{(38)}& {\mathit{\gamma }}_{\text{sun}}=(\mathrm{1}-e^{-{\text{PAI}}_{\mathrm{c}}}){\displaystyle \frac{{\text{PAI}}_{\text{sun,c}}}{{\text{PAI}}_{\mathrm{c}}}}.\end{array}$$

The thermal emissivity of plants and soil is assumed to be 1. The net emission of long-wave radiation by the shaded part of the canopy is described by analogous equations.

The long-wave radiation emitted by the ground surface is

$$\begin{array}{}\text{(39)}& L_{\mathrm{g}}=\mathit{\sigma }{T}_{\mathrm{g}}^{\mathrm{4}}-{\mathit{\gamma }}_{\text{sun}}\mathit{\sigma }{T}_{\text{sun}}^{\mathrm{4}}-{\mathit{\gamma }}_{\text{sha}}\mathit{\sigma }{T}_{\text{sha}}^{\mathrm{4}}+(\mathrm{1}-{\mathit{\gamma }}_{\text{sun}}-{\mathit{\gamma }}_{\text{sha}})L^{\downarrow }.\end{array}$$

The bulk long-wave radiation emitted by the land surface toward the atmosphere is

$$\begin{array}{}\text{(40)}& L^{\uparrow }={\mathit{\gamma }}_{\text{sun}}\mathit{\sigma }{T}_{\text{sun}}^{\mathrm{4}}+{\mathit{\gamma }}_{\text{sha}}\mathit{\sigma }{T}_{\text{sha}}^{\mathrm{4}}+(\mathrm{1}-{\mathit{\gamma }}_{\text{sun}}-{\mathit{\gamma }}_{\text{sha}})\mathit{\sigma }{T}_{\mathrm{g}}^{\mathrm{4}}.\end{array}$$

2.2.4 Assimilation and stomatal conductance

In what follows, variables that are updated daily are denoted with the subscript “day”. Daytime averages are denoted with the subscript “dt”. All the other variables are computed on a sub-daily basis. Photosynthetic assimilation is now calculated within the sub-daily energy balance routine (Fig. 1). A net photosynthetic rate is computed for the sunlit and shaded leaves of each cohort separately by calling the photosynthesis routine built in LPJ-GUESS. This calculation is based on the biochemical model of Collatz et al. (1991, 1992), the strong-optimality model of light use efficiency at the canopy level of Haxeltine and Prentice (1996), and the nitrogen limitation of the maximum carboxylation rate introduced in Smith et al. (2014). The net photosynthetic assimilation is accumulated over the diurnal cycle and subtracted from heterotrophic respiration (R_h, computed daily) to calculate daily net ecosystem exchange (NEE).

For a given cohort i, the maximum carboxylation rate, $V_{\text{max}}^{\left(i\right)}$, is recalculated at the end of every simulation day and depends linearly on the total amount of daily absorbed photosynthetic active radiation, ${\text{PAR}}_{\text{day}}^{\left(i\right)}$ (Haxeltine and Prentice, 1996):

$$\begin{array}{}\text{(41)}& V_{\text{max,day}}^{\left(i\right)}=f_{\mathrm{v}}(T_{\text{leaf,dt}}^{\left(i\right)},\mathrm{\cdots })\times {\text{PAR}}_{\text{day}}^{\left(i\right)}.\end{array}$$

In this equation, $V_{\text{max,day}}^{\left(i\right)}$ is expressed per unit patch area. The slope of the relationship, f_v, encodes the influence of temperature and nitrogen limitation. The calculation of $V_{\text{max, day}}^{\left(i\right)}$ uses the original LPJ-GUESS PAR absorption estimation as described in Sect. 2.2. However, the influence of temperature is calculated using the newly simulated canopy temperature rather than the daily average air temperature. Updating $V_{\text{max,day}}^{\left(i\right)}$ on sub-daily timescales is not necessary because readjustment of leaf nitrogen content and photosynthetic traits occurs on timescales of days to weeks (e.g. Reich et al., 1991; Irving and Robinson, 2006) and therefore cannot follow diurnal environmental variations.

The daytime averaged leaf temperature, $T_{\text{leaf,dt}}^{\left(i\right)}$, is weighted by the daily-averaged fractions of sunlit and shaded leaves for cohort i:

$$\begin{array}{}\text{(42)}& T_{\text{leaf,dt}}^{\left(i\right)}={\displaystyle \frac{\mathrm{1}}{n_{\text{dt}}}}\sum _{\text{dt}}{\displaystyle \frac{{\text{PAI}}_{\text{sun}}^{\left(i\right)}{T}_{\text{sun}}+{\text{PAI}}_{\text{sha}}^{\left(i\right)}{T}_{\text{sha}}}{{\text{PAI}}^{\left(i\right)}}},\end{array}$$

where n_dt is the number of daytime sub-daily periods.

Separating the contributions to daily absorbed PAR from sunlit and shaded leaves, maximum carboxylation rates for the sunlit and shaded parts of the cohort are estimated as

$$\begin{array}{ll}{\displaystyle }{V}_{\text{max,sun,day}}^{\left(i\right)}& {\displaystyle }=f_{\mathrm{v}}(T_{\text{leaf,dt}}^{\left(i\right)},\mathrm{\cdots })\times {\text{PAR}}_{\text{sun,day}}^{\left(i\right)},\\ \text{(43)}& {\displaystyle }{V}_{\text{max,sha,day}}^{\left(i\right)}& {\displaystyle }=f_{\mathrm{v}}(T_{\text{leaf,dt}}^{\left(i\right)},\mathrm{\cdots })\times {\text{PAR}}_{\text{sha,day}}^{\left(i\right)},\end{array}$$

where ${\text{PAR}}_{\text{sun,day}}^{\left(i\right)}$ and ${\text{PAR}}_{\text{sha,day}}^{\left(i\right)}$ are the total daily PAR absorbed by the sunlit and shaded leaves of cohort i, respectively. Combining Eqs. (41) and (43) yields, for sunlit leaves,

$$\begin{array}{}\text{(44)}& V_{\text{max,sun,day}}^{\left(i\right)}=V_{\text{max,day}}^{\left(i\right)}{\displaystyle \frac{{\text{PAR}}_{\text{sun,day}}^{\left(i\right)}}{{\text{PAR}}_{\text{day}}^{\left(i\right)}}}.\end{array}$$

The maximum carboxylation rate per unit leaf area is then calculated as

$$\begin{array}{}\text{(45)}& V_{\text{max,sun,day,leaf}}^{\left(i\right)}=\mathrm{86}{\mathrm{400}}^{-\mathrm{1}}\mathit{\beta }{\displaystyle \frac{V_{\text{max,sun,day}}^{\left(i\right)}}{{\text{LAI}}_{\text{sun,dt}}^{\left(i\right)}}},\end{array}$$

where ${\text{LAI}}_{\text{sun,dt}}^{\left(i\right)}$ is the daily-averaged sunlit LAI of cohort i, and we have introduced a factor β to limit the photosynthetic rate under conditions of water stress. The prefactor 86 400⁻¹ converts the rate from d⁻¹ to s⁻¹. Analogous equations apply to shaded leaves.

The water stress factor β is formulated as a sum over soil layers of a water uptake response function weighed by a PFT-specific vertical rooting profile:

$$\begin{array}{}\text{(46)}& \mathit{\beta }=\sum _j{r}^{\left(j\right)}{W}_{\text{av}}^{\left(j\right)},\end{array}$$

where r^(j) is the fraction of roots in soil layer j. In order to study the impact of the β factor on the model predictions, we implemented four different options for the water uptake response function $W_{\text{av}}^{\left(j\right)}$. In the Noah type (Niu et al., 2011), $W_{\text{av}}^{\left(j\right)}$ decreases linearly in each soil layer with volumetric water content θ^(j) down to the wilting point:

$$\begin{array}{}\text{(47)}& W_{\text{av}}^{\left(j\right)}={\displaystyle \frac{{\mathit{\theta }}^{\left(j\right)}-{\mathit{\theta }}_{\text{wilt}}}{{\mathit{\theta }}_{\text{fc}}-{\mathit{\theta }}_{\text{wilt}}}},\end{array}$$

where θ_wilt and θ_fc are volumetric water content at wilting point and field capacity, respectively. In LPJ-GUESS, the wilting point is assumed to be at a matric potential of ${\mathit{\psi }}_{\text{wilt}}=-\mathrm{45}$ m, and the corresponding soil water content is calculated following Prentice et al. (1992).

The CLM-type water uptake response function is formulated in terms of matric potential (Oleson et al., 2010):

$$\begin{array}{}\text{(48)}& W_{\text{av}}^{\left(j\right)}={\displaystyle \frac{{\mathit{\psi }}_{\text{wilt,CLM}}-{\mathit{\psi }}^{\left(j\right)}}{{\mathit{\psi }}_{\text{wilt,CLM}}-{\mathit{\psi }}_{\text{sat}}}},\end{array}$$

where ψ^(j) is the matric potential of layer j, ψ_sat is the matric potential at saturation, and ψ_wilt,CLM is the matric potential at wilting point, set to −150m. In this case, the water uptake response is flatter than in the Noah-type case when the soil is wet, and decreases more steeply when the soil gets drier. We also implemented a modified version of the CLM-type uptake function, with the same functional form but using LPJ-GUESS's −45m wilting matric potential instead of CLM's −150m.

The SSiB type water uptake response function is

$$\begin{array}{}\text{(49)}& W_{\text{av}}^{\left(j\right)}=\mathrm{1}-e^{-c_{\mathrm{2}}\ln [{\mathit{\psi }}_{\text{wilt}}/{\mathit{\psi }}^{\left(j\right)}]},\end{array}$$

where the parameter c₂ depends on PFT, and takes values between 4.36 and 6.37 (Xue et al., 1991). In this study, we set c₂ to a fixed value of 5.8 for all PFTs, which results in high β values in most of the water availability range, and a steep decrease when approaching the wilting point.

Figure 4 shows the behavior of the different formulations of $W_{\text{av}}^{\left(j\right)}$ as a function of volumetric water content. This type of formulations, which are widely used in LSMs (see Damour et al., 2010, for an overview), are phenomenological relationships that attempt to capture the response of plants to water stress in a rather simplified way (Egea et al., 2011; De Kauwe et al., 2013). Transpiration of soil water by plants is primarily driven by the water potential gradient along the soil–plant–atmosphere continuum. Plants regulate this gradient by opening and closing their stomata in response to environmental factors, including vapour pressure deficit, leaf water potential, and soil water availability, in a way that depends on their hydraulic strategy (a detailed discussion can be found in Lambers et al., 2008). Including a more explicit representation of soil–plant–air hydraulics as well as physiological constraints in an LSM leads to better agreement with observations than the above formulations under soil water stress conditions (Bonan et al., 2014). However, implementing these more complex formulations in ESMs remains a challenge due to a lack of data for broader applicability and computational efficiency tradeoffs (Clark et al., 2015).

Stomatal conductance and photosynthetic rate are related through a semi-empirical model. The photosynthesis rate depends on the CO₂ concentration inside the stomatal cavity. This concentration is related to the atmospheric CO₂ concentration through a diffusion process across the stomatal opening and the leaf boundary layer and therefore depends upon stomatal conductance, which in turn depends on the photosynthetic rate. Hence, photosynthetic rate and stomatal conductance are calculated simultaneously by iteration. A detailed description of the algorithm can be found in Bonan (2019).

As noted above, we implemented two selectable stomatal conductance models. In the Ball–Berry model (Ball et al., 1987), stomatal conductance depends linearly on net assimilation and the fractional humidity at the leaf surface h_s, and inversely on CO₂ concentration at the leaf surface, c_s. The stomatal conductance for sunlit leaves of cohort i is

$$\begin{array}{}\text{(50)}& g_{\text{s,sun}}^{\left(i\right)}=g_{\text{min}}+g_{\text{1,BB}}{\displaystyle \frac{A_{\text{n,sun}}^{\left(i\right)}{h}_{\text{s,sun}}}{c_{\text{s,sun}}}},\end{array}$$

where $A_{\text{n,sun}}^{\left(i\right)}$ is the net photosynthetic rate per unit leaf area, g_min is a minimum stomatal conductance, and g₁ is a PFT-specific parameter. The Medlyn model (Medlyn et al., 2011) is derived from the assumption that stomata optimize CO₂ uptake while minimizing water loss. In this model, stomatal conductance depends inversely on the square root of the vapour pressure deficit at the leaf surface, D_s. The stomatal conductance for sunlit leaves of cohort i is

$$\begin{array}{}\text{(51)}& g_{\text{s},\text{sun}}^{\left(i\right)}=g_{\text{min}}+\mathrm{1.6}\left(\mathrm{1}+{\displaystyle \frac{g_{\text{1,Med}}}{\sqrt{D_{\text{s,sun}}}}}\right){\displaystyle \frac{A_{\text{n,sun}}^{\left(i\right)}}{c_{\text{s,sun}}}}.\end{array}$$

Values of the parameters g_1,BB and g_1,Med for specific PFTs were obtained following Sellers et al. (1996) for the Ball–Berry model and De Kauwe et al. (2015) for the Medlyn model. Figure 3 shows the different behaviour of the stomatal conductance models as a function of D_s.

Figure 3 Stomatal conductance as a function of water vapour deficit at the leaf surface.

Figure 4 Factor limiting plant water uptake as a function of volumetric soil water content. The dashed vertical lines represent, from left to right, the volumetric soil water content at a wilting matric potential of −150 m, at a wilting matric potential of −45 m, and at field capacity.

2.2.5 Soil physics

In standard LPJ-GUESS, soil temperature is used in calculations related to ecosystem respiration and nitrogen cycling, while soil water content influences plant water uptake and evapotranspiration. Both quantities affect soil organic matter decomposition rates.

Soil temperature T_s is now calculated by solving the heat transport equation:

$$\begin{array}{}\text{(52)}& {\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{h}}}}{\displaystyle \frac{\partial }{\partial z}}\left({\mathit{\kappa }}_{\mathrm{s}}{\displaystyle \frac{\partial T_{\mathrm{s}}}{\partial z}}\right),\end{array}$$

where c_h(z) and κ_s(z) are soil heat capacity and thermal conductivity, respectively. The top boundary condition is given by the heat flux into the ground, G, calculated in the energy balance routine (Eq. 11). Heat flow through the bottom boundary of the soil column is neglected. Thermal conductivity is calculated following the method of Johansen (1977). Soil heat capacity is computed as a weighted sum of the heat capacities of the dry soil, which depends on soil texture, and water (de Vries, 1963). Soil organic matter does not contribute to soil heat capacity in the current version of the model.

Vertical water transport in the soil column is described by the Richards equation (Richards, 1931), which can be expressed in the following form:

$$\begin{array}{}\text{(53)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left[{\mathit{\lambda }}_{\mathrm{w}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}-{\mathit{\gamma }}_{\mathrm{w}}\right]+S_{\mathit{\theta }}\left(z\right).\end{array}$$

Here, θ is volumetric water content, λ_w(θ) is hydraulic diffusivity, γ_w(θ) is hydraulic conductivity, and S_θ(z) is a volumetric sink term that accounts for plant water uptake (S_θ≤0). Hydraulic diffusivity and conductivity are calculated as a function of soil texture and soil water content by using the expressions derived by Clapp and Hornberger (1978) and Cosby et al. (1984). Rainwater that is not intercepted by the canopy infiltrates into the soil at a rate limited by the soil's infiltration capacity as given by the Green–Ampt equation (Green and Ampt, 1911). Free gravitational drainage is assumed at the bottom of the soil column.

Soil temperature, water content, ecosystem respiration, plant water uptake and evapotranspiration are calculated in the sub-daily loop. Equations (52) and (53) are solved with a Crank–Nicolson scheme (e.g. Press, 2003). Daily averages of water content and temperature over the layers corresponding to standard LPJ-GUESS top and bottom layers are then used as inputs to the original soil organic matter and nitrogen cycling routines.

3 Model verification and evaluation

3.1 Model verification

The revised model was verified by performing energy and water conservation tests. At any given time step, the energy conservation error per unit time and per unit patch area, Δu_err ($\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), is calculated as

$$\begin{array}{}\text{(54)}& \mathrm{\Delta }{u}_{\text{err}}=S^{\downarrow }+L^{\downarrow }-\langle L^{\uparrow }+H^{\uparrow }+\mathit{\lambda }{E}^{\uparrow }+\mathrm{\Delta }{u}_{\text{soil}}\rangle ,\end{array}$$

where 〈⋅〉 indicates an average over patches, and Δu_soil is the rate of change of energy stored in the soil column per unit patch area. The latter is calculated as

$$\begin{array}{}\text{(55)}& \mathrm{\Delta }{u}_{\text{soil}}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\sum _j{c}_{\mathrm{h}}^{\left(j\right)}\mathrm{\Delta }{z}^{\left(j\right)}{T}_{\mathrm{s}}^{\left(j\right)},\end{array}$$

where Δt is the time step in seconds, and $c_{\mathrm{h}}^{\left(j\right)}$, Δz^(j) and $T_{\mathrm{s}}^{\left(j\right)}$ are, respectively, the heat capacity, thickness and temperature of soil layer j. Figure 5 shows the frequency of the energy conservation error relative to the energy input to the system (i.e. the total incoming irradiance, $S^{\downarrow }+L^{\downarrow }$). The vast majority of the time steps (∼99.8 %) the error is smaller than 0.2 % of the incoming radiation, and the error is never larger than 0.85 % of the energy input. The mean energy balance closure error is ∼0.013 % of the energy input (σ∼0.035 %).

Figure 5 The histogram shows the energy conservation error, as a percentage of the energy input, incurred at every time step. The symbols indicate the mean absolute error corresponding to each bin. The error bars indicate ±1σ around the mean. The plots are derived from data from the historical period of all LSM simulations in this study.

The water conservation error is computed as

$$\begin{array}{}\text{(56)}& \mathrm{\Delta }{w}_{\text{err}}=P-\langle R+E^{\uparrow }+\mathrm{\Delta }{w}_{\text{soil}}+\mathrm{\Delta }{w}_{\mathrm{c}}\rangle ,\end{array}$$

where P is precipitation, R is runoff (including surface runoff and base flow), E^↑ is evapotranspiration, Δw_soil is the change in soil water content per unit patch area, per unit time, and Δw_c is the change in canopy water content. This error is, on average, $\sim \mathrm{5}\times {\mathrm{10}}^{-\mathrm{12}}$ % of the precipitation at any given time step ($\mathit{\sigma }\sim \mathrm{9}\times {\mathrm{10}}^{-\mathrm{12}}$), and never surpasses $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{8}}$ % of the water input.

We therefore conclude that the magnitude of the errors in energy balance closure and water conservation are negligible.

3.2 Evaluation set-up

We evaluated the revised model by comparing hourly and monthly simulated fluxes of sensible and latent heat and annual CO₂ fluxes, with flux tower measurements from 21 FLUXNET2015 (Pastorello et al., 2020) sites. The current version of the model does not simulate snow or frozen soil water, so we restricted our study to sites where the air temperature remained above 0 ^∘C throughout the measuring period. We additionally discarded wetland sites, which require a more detailed representation of soil water and groundwater hydrology (Wania et al., 2009). A list of the selected sites is presented in Table 2. The locations of the sites are represented on the world map in Fig. 6.

Table 1 List of plant functional types in the standard configuration of LPJ-GUESS (only PFTs predicted by the simulations in this study are listed).

Schroder et al. (2014)López-Ballesteros et al. (2017)Hutley et al. (2011)Beringer et al. (2011)Merbold et al. (2009)S. Wolf et al. (2011)Hutley et al. (2011)Beringer et al. (2011)Ardö et al. (2008)Beringer et al. (2011)Beringer et al. (2016)Beringer et al. (2011)Bristow et al. (2016)Beringer et al. (2016)Saleska et al. (2003)Saleska et al. (2003)Bonal et al. (2008)Stefani et al. (2009)Kosugi et al. (2008)S. Wolf et al. (2011)Merbold et al. (2009)

Table 2 Brief description of selected sites. The land cover classification was taken from the FLUXNET site description web pages. The reference level height is taken as the height of the measuring sensors above the canopy. A dash indicates that we were not able to find an observed LAI value for the site.

Figure 6 FLUXNET sites selected for model evaluation. Different symbols indicate different land cover types. The sites are labelled according to their site code (Table 2).

For each site, we ran eight simulations, covering all possible configurations of the water uptake response functions and stomatal conductance schemes described in Sect. 2.2.4 (Table 3). We used the climate data collected at the tower sites to force the model. Half-hourly forcing data were converted to hourly averages to use a fixed time step of 1 h in all simulations. We set a lower boundary of 10 % of the dataset median on the air humidity to correct for physically invalid negative values. Nitrogen deposition data are from Lamarque et al. (2013). The soil texture data used to calculate soil hydraulic and thermal properties (as described in Sect. 2.2.5) at each site were as in Sitch et al. (2003), based on the Digitized Soil Map of the World (Zobler, 1986; FAO, 1991). Atmospheric CO₂ concentration data are from McGuire et al. (2001). Additionally, we ran a standard (non-LSM) LPJ-GUESS simulation to compare both model versions' predictions of monthly evapotranspiration and a number of ecosystem composition and function variables. The number of replicate patches was set to 100 in all the simulations to avoid spurious effects of the stochastic ecosystem processes on the modelled fluxes. Simulation of wildfires was switched off in all simulations.

Table 3 Summary of the LPJ-GUESS/LSM simulations carried out. Simulations with different stomatal conductance schemes are arranged in columns: Ball–Berry (BB) and Medlyn (Med). Simulations with different water uptake response function types are arranged in rows: NOAH, CLM, modified CLM, and SSiB.

All natural PFTs were allowed to establish in forest and savanna sites. Since the focus of the model evaluation was placed on the turbulent fluxes, we restricted the simulated PFTs to grassy types at sites classified as grasslands, which limits modelled surface roughness. This was also done for Spain-Amoladeras and Congo-Tchizalamou. Amoladeras is classified as an open shrubland on the FLUXNET reference, but the vegetation is short and the most abundant species is Machrocloa Tenacissima, a type of grass (López-Ballesteros et al., 2017). Tchizalamou, which is classified as savanna, is actually a C₄ grassland (Merbold et al., 2009).

The simulations were spun up from a bare ground state following a standard procedure that combines 500 simulation years with a semi-analytic calculation of the equilibrium size of the soil organic matter pools (see Supplement), to bring C and N soil and vegetation pools to near equilibrium with the climate. During the spin-up phase, the site climate spanning the whole measurement period was repeated cyclically, with interannual trends in air temperature removed, and the atmospheric CO₂ concentration was kept at the level of the first year of observations at each site.

3.3 Analysis

Half-hourly measured fluxes were converted to hourly averages for direct comparison with model outputs. Sub-daily FLUXNET data are classified into four quality categories: 0 (measured), 1 (good-quality gap fill), 2 (poor-quality gap fill), and 3 (downscale from ERA reanalysis data). In our analysis, we only used sub-daily fluxes with a quality flag of 0 or 1. For monthly and annual fluxes, the quality flag varies between 0 and 1 and indicates the fraction of the sub-daily values in that month/year whose quality is either 0 or 1. We only used monthly and annual fluxes with a quality flag equal to or greater than 0.75. Following Stöckli et al. (2008), we further discarded fluxes with friction velocity $u^{\ast }<\mathrm{0.2}$ m s⁻¹ in order to avoid possibly biased eddy covariance measurements during periods of weak turbulence (Schroder et al., 2014).

Figure 7 LAI values for the spin-up period at three selected sites: PA-SPs (panels a–c), BR-Sa1 (panels d–f), and AU-Dry (panels g–i). The columns correspond to standard LPJ-GUESS (a, d, g), CLM/BB (b, e, h), and CLM/Med (c, f, i) simulations. The time series were smoothed for better visualization by applying a 15-year running average.

To evaluate the agreement between measured and simulated turbulent heat fluxes at each site for all different model configurations we used standard statistical metrics: correlation coefficient (r), mean bias, and root mean square error (RMSE). We also considered the standard deviation of the modelled fluxes normalized by the standard deviation of the observed fluxes (${\mathit{\sigma }}_{\mathrm{m}}/{\mathit{\sigma }}_{\mathrm{o}}$), which provides a measure of the agreement between observed and simulated variability.

Figure 8 Panels (a, b): percent change in average gross primary production (blue), autotrophic respiration (orange), and net primary production (green), simulated by the LSM version with respect to standard LPJ-GUESS. Panels (c, d): percent change in predicted average net primary production (green), heterotrophic respiration (brown), and net ecosystem exchange (pink).

Table 4 Foliar projective cover, averaged over the whole simulated period, of the plant functional types predicted for each site, given as a percentage. The LSM simulations use the CLM-type water uptake response factor. The dominant PFT for each site is highlighted in bold font.

3.4 Results

3.4.1 Ecosystem composition and function

The emerging ecosystem composition in both LSM runs is similar to the standard LPJ-GUESS prediction over forests and grasslands, but it is sensitive to the choice of stomatal conductance scheme at some savanna and woody savanna sites and at ZM-Mon (Table 4). Figure 7 shows the LAI evolution of the established PFTs over the spin-up period for the CLM/BB, CLM/Med, and standard LPJ-GUESS simulations at three selected sites. All three simulations predict a C₄ grassland at PA-SPs, but LAI values are much higher in the LSM simulations (∼11) than the LPJ-GUESS prediction (∼6.5). At BR-Sa1 (a tropical rainforest), the species composition is similar for the three simulations, but LAI values are lower in the LSM runs (∼5.5 vs. ∼6.2). At AU-Dry, the use of different stomatal conductance schemes causes a shift in PFT composition. The BB simulation favours evergreen trees, while the PFT mix is dominated by raingreen trees in the Med simulation, a prediction closer to standard LPJ-GUESS. We found this behaviour to be representative of how the soil water uptake response factor and the stomatal conductance scheme influence the PFT composition at most savanna and woody savanna sites in the LSM simulations. A stronger limitation on transpiration (e.g. the NOAH-type water uptake response factor or the Ball–Berry stomatal conductance model) results in higher soil water content throughout the year, which promotes stronger growth of evergreen trees.

Model predictions for the rest of the selected variables are shown in Table 5. The two C₃ grassland sites show different behaviour with respect to ecosystem productivity and respiration. At AU-Emr, LSM simulations predict substantially lower gross primary production (GPP) and autotrophic respiration (R_a) than standard LPJ-GUESS, which results in lower estimates of net primary production (NPP). This site is a net carbon source (positive NEE) in all three simulations, which agrees with observations. At ES-Amo, the NPP increase in the LSM runs is larger than the decrease in heterotrophic respiration (R_h), resulting in an enhanced carbon sink compared with standard LPJ-GUESS.

Table 5 Comparison of selected variables related to simulated ecosystem structure and function between standard LPJ-GUESS and the LSM version at the selected sites. The LSM values are from the CLM/BB and CLM/Med simulations. Gross primary production (GPP), autotrophic respiration (Ra), net primary production (NPP), heterotrophic respiration (Rh), and net ecosystem exchange (NEE) are given in $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Bold fonts in the LAI and NEE columns indicate the closest match to the observed value. Bold fonts in the rest of the columns indicate the LSM prediction closest to standard LPJ-GUESS.

At PA-SPn, both NPP and R_h decrease in the LSM simulations, but the former decreases less than the latter, resulting in slightly weaker carbon sinks in the LSM simulations. The three simulations predict carbon fluxes much smaller than the measured value of −458 $\mathrm{g}{\mathrm{C}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$. Predictions for ZM-Mon show some differences between runs, but NPP and R_h are similar in all three simulations, resulting in carbon sinks of $\sim -\mathrm{62}$ $\mathrm{g}{\mathrm{C}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$. This result is inconsistent with measurements at the site, which indicate a carbon source of 143 $\mathrm{g}{\mathrm{C}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$.

Differences in simulated carbon fluxes between standard LPJ-GUESS and the CLM/BB and CLM/Med runs for the remaining land cover types are summarized in Fig. 8. Both LSM runs predict, on average, higher GPP and R_a values than the non-LSM simulation over C₄ grassland, savanna, and woody savanna sites. This results in an increased average NPP value in C₄ grasslands (∼18 % in the CLM/BB run and ∼31 % in the CLM/Med run) and a decreased average NPP at woody savanna sites ($\sim -\mathrm{11}$ % and $\sim -\mathrm{7}$ % in the CLM/BB and CLM/Med runs, respectively). At savanna sites, the increase in GPP in both LSM simulations is similar (∼10 %), but the increase in R_a is much higher for CLM/BB, which leads to changes in NPP of $\sim -\mathrm{10}$ % in the CLM/BB run and ∼6 % in the CLM/Med run. At forest sites, the balance between decreased values of GPP and R_a results in lower NPP values in the LSM simulations. Average values of R_h in the CLM/BB simulation increase over C₄ grasslands and decrease over woody savannas and evergreen forests. The CLM/Med simulation shows the same pattern except over savanna sites, where R_h increases by ∼5 % with respect to standard LPJ-GUESS. Over woody savannas, the average NEE change is $\sim -\mathrm{129}$ % for the CLM/BB run and ∼122 % for the CLM/Med run.

The above-described discrepancies between standard LPJ-GUESS and the LSM versions stem from the different physical environments simulated in the models. Calculating assimilation at the newly simulated canopy temperature, rather than the air temperature, can lead to either higher or lower productivity, depending on the optimal photosynthetic temperature ranges of each PFT and the impact of temperature on nitrogen limitation (Sect. 4.2). Canopy temperature also affects autotrophic respiration, while differences in the simulated soil humidity and temperature impact organic matter decomposition rates and heterotrophic respiration. The combination of these effects results in differences in simulated carbon and nitrogen pools and NEE (we have included a comparison between soil carbon and nitrogen pools simulated by standard LPJ-GUESS and LPJ-GUESS/LSM in the Supplement).

The large relative changes in NEE between simulations result from small discrepancies in magnitude. Figure 9 shows a comparison between land cover averages of measured and modelled NEE for C₄ grasslands, savanna, woody savanna, and evergreen forests. Average measured NEE is negative for all land cover types and substantially more negative than in the simulations for savanna, woody savanna, and evergreen broadleaf forests, implying an average underestimation of the C sink by the models at these sites. At C₄ sites simulations predict NEE values between −88 and −110 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, while observations indicate a less negative value of −33 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. For savanna, measured NEE is −221 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, while simulations predict an average between −48 and −56 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. For woody savanna, measured NEE averages to −238 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, while simulated fluxes range between −28 and 3 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. Simulated fluxes at evergreen broadleaf forests are, on average, between −84 and −130 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, while measurements indicate an average NEE of −396 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$. However, this is the result of very large negative values measured at AU-Rob and MY-PSO (Table 5). In general, differences in simulated fluxes between standard LPJ-GUESS and the LPJ-GUESS/LSM simulations are small compared with the magnitude of observed fluxes, and the interannual and cross-site variabilities of the measured fluxes are much greater than in the simulations. The discrepancies between observed and simulated NEE magnitude and variability reflect the fact that, in the simulations, the carbon pools are all close to equilibrium with the climate and atmospheric CO₂ concentration as a result of the spin-up procedure described in Sect. 3.2. Differences between observed and simulated NEE values are to be expected because we did not attempt to reproduce site history, including age, disturbance, and legacies arising from historical trends in CO₂ concentration.

Figure 9 Comparison between observed and modelled annual NEE. The symbols indicate averages over sites with the same land cover type. Red triangles correspond to flux tower CO₂ measurements. Blue dots, green squares, and purple crosses correspond, respectively, to the CLM/BB, CLM/Med, and standard LPJ-GUESS simulations. The bars represent 1 standard deviation above and below the average.

3.4.2 Annual and diurnal cycles of turbulent heat fluxes

Figure 10 shows examples of simulated and observed monthly averages of turbulent and latent heat fluxes over the course of a year at four sites: Gingin (AU-Gin), Daly River Savanna (AU-DaS), Santarem-Km67 (BR-Sa1), and Guyaflux (GF-Guy). Examples of the monthly-averaged diurnal cycle for the same sites are shown in Figs. 11 and 12. We chose these sites and years to illustrate situations with varying degrees of agreement between simulations and measurements. The simulated fluxes are from the run using the CLM-type water uptake response function and the Medlyn model of stomatal conductance (CLM/Med).

Figure 10 Observed and simulated annual cycles of sensible (H) and latent (λE) heat flux at four selected sites: Gingin (AU-Gin, panel a), Daly River Savanna (AU-DaS, panel b), Santarem-Km67 (BR-Sa1, panel c), and Guyana (GF-Guy, panel d). Mean monthly precipitation and observed and modelled net radiation (R_n) are plotted for reference.

Figure 11 Monthly-averaged diurnal cycle of sensible and latent heat flux at the AU-DaP (panels a–f) and AU-DaS (panels g–l) sites in selected months. The red and blue lines represent simulated sensible and latent heat fluxes, respectively. The shaded areas around each curve delimit 1 standard deviation above and below it. The symbols represent monthly-averaged fluxes. The error bars indicate a ±1σ deviation from the observed mean.

Figure 12 Monthly-averaged diurnal cycle of sensible and latent heat flux at the BR-Sa1 (panels a–f) and GF-Guy (panels g–l) sites in selected months. The red and blue lines represent simulated sensible and latent heat fluxes, respectively. The shaded areas around each curve delimit 1 standard deviation above and below it. The symbols represent monthly-averaged fluxes. The error bars indicate a ±1σ deviation from the observed mean.

At the AU-Gin site, the shapes of the annual cycles of latent and sensible heat are well-reproduced in the simulations (Fig. 10a). Sensible heat is largest at the beginning of the year, decreases steeply to its minimum around June–July, and starts increasing again around August. The simulation agrees well with measurements most of the year but overestimates sensible heat by ∼45 W m⁻² in the first 2 months. Observed latent heat increases at the start of the wet season and dominates the turbulent exchange from May to September. Simulated latent heat is overestimated by up to ∼25 W m⁻² in the wet season and underestimated in the dry season. The shift from larger sensible heat to larger latent heat in May is well captured in the simulation, but, due to the overestimation of latent heat, the shift back to larger sensible heat flux is delayed by about 2 months with respect to the observations. The average simulated diurnal cycle of sensible heat is overestimated in January, peaking at ∼700 W m⁻² (observed: ∼500 W m⁻²), while it agrees very well with observations in May and September, both in terms of magnitude and day-to-day variability (Fig. 11a–c).

At the AU-DaS site (Fig. 10b), observed and simulated heat fluxes diverge substantially during the dry season (July–November). Simulated monthly averages of latent heat are ∼20–30 W m⁻² above measured values from March to May and ∼30–45 W m⁻² below the measurements between August and October. The average simulated latent heat diurnal cycle peaks at ∼350 W m⁻² in May (observed: ∼ 175 W m⁻²) and at ∼ 25 W m⁻² in September (observed: ∼ 145 W m⁻²; Fig. 11j–l). This marked divergence from measured values happens in very dry periods, when the simulated soil moisture in the rooting zone drops close to the wilting point and there is not enough precipitation to replenish it until the start of the wet season. As a consequence, sensible heat is greatly overestimated. Simulated monthly averages rise sharply and peak at ∼ 120–140 W m⁻² from September to October, while measured values stay at ∼60 W m⁻² throughout the dry season. The average sensible heat diurnal cycle peaks at ∼530 W m⁻² in September, while the observed average diurnal peak is slightly under ∼300 W m⁻² (Fig. 11g–i).

Monthly averages of sensible and latent heat at the BR-Sa1 tropical rainforest site show little variability throughout the year (Fig. 10c). Measured sensible heat flux stays at ∼20 W m⁻² for most of the year and increases to ∼30 W m⁻² around August and September, when measured precipitation reaches its minimum. During this period, the soil retains enough moisture in the rooting zone to maintain average latent heat levels at ∼80–90 W m⁻². Sensible and latent heat fluxes are systematically overestimated by the model by ∼10–20 W m⁻². This overestimation takes place even when simulated net radiation is very close to observations (June to November), so assuming the measurements do not underestimate the fluxes, it must be compensated by an underestimation of ground heat. One possible contributing factor is an underestimation of the heat conductivity, which could be caused by an underestimation of soil moisture in the top soil layer. Unfortunately, soil moisture measurements are not available for this site, so we were not able to test this hypothesis. Average sensible heat flux peaks daily between ∼160 and 275 W m⁻² (measured: ∼100–150 W m⁻²). Latent heat flux peaks daily between ∼280 and 340 W m⁻² (measured: ∼220–320 W m⁻², Fig. 12a–f).

At the GF-Guy site, another tropical rainforest, monthly averages of sensible heat are overestimated by ∼20 W m⁻² throughout the year, while latent heat flux is underestimated by about the same amount. The simulated sensible heat diurnal cycle peaks, on average, by ∼100 W m⁻² above the measured values, while the peak of the simulated latent heat diurnal cycle is ∼130 W m⁻² below the measured values. There is a marked decrease in simulated latent heat in October and a corresponding sharp increase in sensible heat, due to excessively low soil moisture in the rooting zone in the model. The simulated October average diurnal sensible heat cycle peaks at ∼350 W m⁻² (measured: ∼190 W m⁻²), while the average latent heat diurnal cycle peaks at ∼180 W m⁻² (measured: ∼360 W m⁻²).

Figure 13 Performance of the CLM/BB (a, c) and CLM/Med (b, d) runs for sensible heat flux. The Taylor diagrams (Taylor, 2001) in panels (a) and (b) summarize the degree of agreement between observed and simulated hourly fluxes by relying on the geometrical relationship between the “centred pattern” root mean square difference (defined as $E^{\prime \mathrm{2}}={\text{RMSE}}^{\mathrm{2}}-{\text{Bias}}^{\mathrm{2}}$), the correlation coefficient, and the standard deviation of observed and modelled data: $E^{\prime \mathrm{2}}={\mathit{\sigma }}_{\mathrm{o}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{m}}^{\mathrm{2}}-\mathrm{2}{\mathit{\sigma }}_{\mathrm{o}}{\mathit{\sigma }}_{\mathrm{m}}r$. Each point in the polar diagram represents a simulation. The radial coordinate indicates the ratio between modelled and observed standard deviations. The correlation between observed and modelled values is encoded by the polar angle; it decreases anticlockwise from r=1 (perfect correlation) for points situated on the x axis to r=0 for points situated on the y axis. The distance between a point in the diagram and the reference value 1 on the x axis equals the centred pattern RMSE normalized by the standard deviation of the observed values, $E^{\prime }/{\mathit{\sigma }}_{\mathrm{o}}$, and is therefore a measure of the agreement between observed and simulated data. The scatter plots (panels c, d) show a direct comparison of observed and modelled monthly-averaged fluxes. The different symbols refer to different land cover types: savanna (SAV), woody savanna (WSA), C₃ grasslands (C3G), C₄ grasslands (C4G), evergreen broadleaf forest (EBF), and deciduous broadleaf forest (DBF).

3.4.3 Influence of different stomatal conductance schemes on the simulated heat fluxes

Table 6 and Fig. 13 show model performance statistics for sensible heat fluxes for the CLM/BB and CLM/Med simulations. Correlations between modelled and observed sensible heat fluxes are very high and similar for both runs. For hourly fluxes, r is between ∼0.84 and 96. Correlations between monthly-averaged fluxes are weaker but still high at most sites (r>0.75), but they are very low for SD-Dem, AU-RDF, and ZM-Mon. At these three sites, the CLM/Med simulation shows better correlations and smaller errors than the one using the BB scheme.

Table 6 Model performance statistics for simulated hourly (left) and monthly (right) sensible heat fluxes for the CLM/BB and CLM/Med simulations. Bold fonts indicate the model configuration that performed better. The mean and standard deviation of the observed fluxes (${\stackrel{\mathrm{‾}}{H}}_{\mathrm{o}}$ and σ_o, respectively), shown for reference, are given in W m⁻². The RMSE and bias have been normalized by the mean of the observed fluxes for easier cross-site comparison.

The model tends to overestimate average sensible heat. The hourly and monthly mean biases are non-negative at all sites (except at CG-Tch for monthly fluxes, where it is slightly negative), but normalized RMSE and mean bias are smaller for the CLM/Med run at most sites. The simulations seem to perform comparatively better in grasslands; for the Med simulation, RMSE is between 0.4 and 0.9 of the sample average (hourly fluxes), whereas it is in the 0.6–1.6 range at savanna sites and in the 0.5–2.5 range at forest sites.

The variability of sensible heat flux is also overestimated by the model. In this case, the CLM/Med run performs better than the CLM/BB run for hourly fluxes, but the situation is the reversed for monthly average fluxes. Again, the simulations show better performance at grassland sites; for hourly fluxes, the Med simulation predicts ${\mathit{\sigma }}_{\mathrm{m}}/{\mathit{\sigma }}_{\mathrm{o}}\sim \mathrm{1.0}$–1.5 in grasslands, ∼1.3–1.8 in savanna and woody savanna sites, and ∼1.1–2.2 in forest sites.

Model performance statistics for latent heat fluxes are presented in Table 7 and Fig. 14. Correlations for hourly fluxes are between 0.7 and 0.9 for most sites. For monthly fluxes, correlations are poorer at forest sites but errors are comparatively small; normalized RMSE is in the 0.1–0.5 range. Hourly correlations are similar for both model configurations at most sites.

Table 7 Model performance statistics for simulated hourly (left) and monthly (right) latent heat fluxes for the CLM/BB and CLM/Med simulations. Bold fonts indicate the model configuration that performed better. The mean and standard deviation of the observed fluxes ($\mathit{\lambda }{\stackrel{\mathrm{‾}}{E}}_{\mathrm{o}}$ and σ_o, respectively), shown for reference, are given in W m⁻². The RMSE and bias have been normalized by the mean of the observed fluxes for easier cross-site comparison.

Figure 14 Performance of the CLM/BB (a, c) and CLM/Med (b, d) runs for latent heat flux. The Taylor diagram shows statistical metrics calculated from hourly observed and simulated fluxes. The scatter plots show a direct comparison of observed and modelled monthly-averaged fluxes. The different symbols refer to different land cover types: SAV, WSA, C3G, C4G, EBF, and DBF.

Latent heat fluxes tend to be underestimated in forest and savanna sites and overestimated over grasslands. The CLM/BB configuration seems to perform better at grassland sites, while the CLM/Med configuration performs slightly better at forest sites. Results for savanna sites are mixed in terms of RMSE, but the CLM/Med scheme yields somewhat smaller biases.

The variability of simulated latent heat fluxes is larger in the CLM/Med run than in the CLM/BB run. Over C₃ grasslands, both LSM runs predict a much larger variability than the observed one. At savanna and woody savanna sites, the CLM/Med run predicts a larger variability than the observed one, both in the hourly and monthly cases, whilst the CLM/BB simulation tends to produce a lower variability. For forest sites, both runs yield σ_m ≲ σ_o, with the exception of BR-Sa3 in the CLM/Med run, where the variability of the monthly fluxes is ∼1.6 times larger than the observed one. On average, the CLM/BB simulation shows better agreement with measured variability over grasslands, while CLM/Med performs somewhat better at forest sites.

3.4.4 Alternative model configurations

To evaluate the overall performance of the different model configurations, we considered the cross-site averaged statistics of each simulation (Table 8). Since the covered period differs across sites, this method ensures all sites contribute equally to the result.

Table 8 Cross-site averaged model performance statistics for simulated hourly and monthly sensible and latent heat fluxes. RMSE and bias are given in W m⁻². Bold fonts indicate the best-performing simulations in each metric.

Figure 15 shows the cross-site averaged metrics on a Taylor diagram. The clumping and clear separation of simulations using different stomatal conductance schemes suggest that this component of the model has a greater influence than the soil water uptake response function on the behaviour of the model, with the possible exception of the linear response function (Noah type), which is much more restrictive than the other three in terms of water uptake. In this case, the variability of modelled latent heat fluxes in the Med simulation is closer to the observed average, lying somewhere between the Med and BB clumps.

Figure 15 Average performance statistics for each model configuration, obtained from modelled and measured hourly (a, b) and monthly (c, d) fluxes of sensible and latent heat flux. The different symbol shapes represent different water uptake response functions (squares: NOAH type; circles: CLM type; crosses: modified CLM type; triangles: SSiB type), and the different colours represent different stomatal conductance schemes (green: Ball–Berry type; black: Medlyn). The red star represents average performance statistics for monthly latent heat fluxes derived from the standard LPJ-GUESS simulation.

Simulated sensible heat fluxes display similar correlation with observations in all runs. The correlation coefficient is very high (r∼0.9) for hourly fluxes and moderately high (r∼0.75) for monthly averages.

Sensible heat is overestimated in all model configurations; the average bias is always positive, but the Med simulations perform better in this respect. In the case of hourly averages, BB runs show an average bias of ∼46 W m⁻², while the average value for the Med runs is ∼35 W m⁻². Average errors are also smaller in the Med simulations. For hourly fluxes, the average RMSE is ∼90 W m⁻² for BB runs and ∼81 W m⁻² for Med runs. For monthly fluxes, RMSE averages are ∼31 and ∼28 W m⁻², respectively.

The model also generally overestimates the variability of sensible heat. For hourly fluxes, the standard deviation of the sample is, on average, ∼1.6 times greater than the measurements for the BB runs and ∼1.5 for the Med runs. In the case of monthly variability, BB runs perform slightly better; the average standard deviations of modelled fluxes are ∼1.4 and ∼1.6 for the BB and Med runs, respectively.

Correlations between modelled and measured latent heat fluxes are lower than for sensible heat: r∼0.8 for hourly fluxes and ∼0.7 for monthly fluxes. All runs show a similar RMSE: ∼75 and ∼24 W m⁻² for hourly and monthly fluxes, respectively. Latent heat is underestimated on average in all configurations. However, the Med runs perform significantly better than the BB runs with this metric. The average bias is $\sim -\mathrm{10}$ W m⁻² (BB: $\sim -\mathrm{20}$ W m⁻²) for hourly fluxes and $\sim -\mathrm{5}$ W m⁻² (BB: $\sim -\mathrm{11}$ W m⁻²) for monthly fluxes. The variability of hourly latent heat fluxes is underestimated in the BB runs by about the same amount that it is overestimated in the Med runs, but in the case of monthly fluxes, BB simulations seem to reproduce the measured variability better (σ_m∼σ_o, while σ_m∼1.3σ_o for Med runs).

Monthly averages of latent heat simulated by the non-LSM version of LPJ-GUESS show a slightly worse correlation with measurements than the LSM version of the model. The average bias is $\sim -\mathrm{2}$ W m⁻², in line with the CLM/Med simulation and significantly lower than the BB simulations, and the RMSE is slightly higher but close to the LSM runs. However, the predicted variability is significantly exaggerated; the standard deviation of the sample of modelled monthly fluxes is, on average, ∼1.7 times larger than observed.

4 Discussion and summary

In this work we described a number of modifications to the LPJ-GUESS DGVM aimed at making the model suitable for direct coupling with an atmospheric model. The newly incorporated energy balance module resolves the diurnal cycle of energy and water fluxes between the canopy and the atmosphere, as opposed to LPJ-GUESS's daily calculations. Calculating these fluxes on a sub-daily basis is necessary to match the shorter time step at which atmospheric models operate (typically 1 h or shorter, depending on resolution). The original daily PAR absorption calculations were replaced with a more sophisticated radiative transfer scheme by adapting the models of Sellers (1985) and Dai et al. (2004) to LPJ-GUESS's multi-cohort, multi-layer canopy (some differences in PAR absorption calculated by both schemes are shown in the Supplement). This enables the model to simulate the upwelling short-wave radiation flux on sub-daily timescales. Direct and diffuse radiations are treated separately, which allows us to resolve sunlit and shaded leaves in the canopy. This approach offers a reasonable compromise between accuracy of the modelled fluxes and computational efficiency (Wang and Leuning, 1998). The representation of soil physical processes was modified in two ways. Firstly, the original 1.5 m-deep, two-layer soil column was replaced with a 3 m-deep, nine-layer column. Secondly, the soil heat and water transport schemes were replaced with less parametrized formulations. Soil heat transport is now calculated by solving the heat diffusion equation, while soil water transport is solved by applying Richards' equation. These formulations are better fit to resolve near-surface heat and water fluxes on the sub-daily timescales introduced in the model.

The new physical schemes introduced in this work lead to discrepancies between the LSM and standard LPJ-GUESS, which stem from the different physical environments simulated in the models. Calculating assimilation at the newly simulated canopy temperature, rather than the air temperature, can lead to either higher or lower productivity, depending on the optimal photosynthetic temperature ranges of each PFT and the effect of temperature on N limitation (see Sect. 4.2). Canopy temperature also affects autotrophic respiration, while differences in the simulated soil humidity and temperature impact organic matter decomposition rates and heterotrophic respiration. The combination of these effects results in differences in the simulated equilibrium carbon and nitrogen pools (see the Supplement) and ecosystem–atmosphere carbon fluxes.

4.1 Evaluation of the simulated heat fluxes

The new model was evaluated by comparing simulated fluxes of sensible and latent heat with flux tower measurements at 21 FLUXNET sites. Stöckli et al. (2008) used a similar analysis (including the filtering of fluxes with $u^{\ast }<\mathrm{0.2}$ m s⁻¹) to evaluate the improvement in the performance of CLM 3.5 after introducing nitrogen limitation of photosynthesis, a groundwater model, and an updated formulation of surface resistance. Owing to the different site selection, a rigorous comparison of that study with the results presented in Sect. 3.4.3 is not possible, but average statistics of model performance can provide an overview of how the models compare. Table 9 shows averaged values of the correlation coefficient and RMSE for CLM 3.0, CLM 3.5, and LPJ-GUESS/LSM. Our model seems to yield stronger correlations between measured and observed sensible heat fluxes and RMSE values similar to CLM 3.5, while CLM 3.5 appears to perform better in terms of RMSE for hourly latent heat fluxes, and the correlation between measured and observed monthly latent heat fluxes is stronger. In order to ascertain the significance of these findings, a comparison using the same site-measured fluxes and forcing climate would be needed. Nevertheless, the values presented in Table 9 suggest that the performance of our model is closer to CLM 3.5 than to CLM 3.0.

Table 9 Average correlation and RMSE between observed and simulated sensible and latent heat fluxes for LPJ-GUESS/LSM, CLM 3.0, and CLM 3.5. RMSE values (between brackets) are given in W m⁻². Hourly and monthly values are labelled (h) and (m), respectively. CLM values correspond to averages over temperate, tropical, and grassland sites as reported in Stöckli et al. (2008). LPJ-GUESS/LSM values correspond to the CLM/Med simulation.

Figure 16 Effect of temperature (T) on modelled V_max; N_leaf,opt: leaf nitrogen content necessary to attain the maximum carboxylation rate; N_leaf: representative leaf nitrogen concentration; ${\widehat{V}}_{\text{max}}$: normalized maximum carboxylation rate without nitrogen limitation; ${\widehat{V}}_{\text{max,lim}}$: normalized maximum carboxylation rate with nitrogen limitation. The histograms show the frequency of temperatures in the AU-DaP simulation; T_atm,day: daily average of air temperature; T_leaf,dt: daytime average of leaf temperature. The shaded area indicates the temperature range where V_max is nitrogen-limited.

Table 10 Daily average climate measured at the four C₄ grassland sites in this study. T_atm,day: daily average temperature (^∘C); P_day: daily average precipitation (mm d⁻¹); $I_{\text{0,day}}^{\downarrow }$: daily average incoming solar irradiance (W m⁻²).

Sensible heat is generally overestimated by the LSM. Poor performance in sensible heat flux estimation is a common issue of many land surface models (Best et al., 2015). The reason for this is not well-understood. It has been suggested that the models, the majority of which use similar methods to calculate the turbulent fluxes, do not extract all the information available in the climate forcing data. However, eddy covariance measurements often fail to close the energy balance and might systematically underestimate sensible heat much more than latent heat, which would appear as an overestimation of sensible heat in the simulations. A detailed discussion of these issues is provided in Haughton et al. (2016).

Table 11 Observed and simulated LAI at the four C₄ grassland sites when all natural PFTs are allowed to grow in the patch. The BB and Med values correspond to LSM simulations using the CLM-type water uptake response function.

One issue in our simulations is the marked underestimation of latent heat flux during extremely dry periods, when the rooting zone is nearly depleted of water available for plant uptake. This, in turn, causes a strong spike in sensible heat (Fig. 10b). All eight model configurations show this behaviour. One possible reason for this is the choice of free drainage boundary conditions at the bottom of the soil column. Simulating groundwater in the model may promote the retention of some soil moisture during dry periods and thus help alleviate this problem (Stöckli et al., 2008). Deeper root profiles and lateral access to soil water may also be important for supporting evapotranspiration in dry periods (Schenk and Jackson, 2002).

We implemented two different stomatal conductance schemes: the Ball–Berry model (Ball et al., 1987) and the Medlyn model (Medlyn et al., 2011). One notable difference between these two models concerns the behaviour of the stomatal conductance when the vapour pressure deficit at the leaf surface (VPD) is small. In the Ball–Berry model, stomatal conductance increases linearly with decreasing VPD, while in the Medlyn model stomatal conductance increases much more rapidly as VPD approaches zero (Fig. 3). Larger stomatal conductance leads to generally higher evapotranspiration values (less negative bias values, Table 8) and enhanced GPP (Fig. 8) in simulations using the Medlyn model. A statistical evaluation of the impact of these differences on the model output was not carried out, but the clumping of symbols representing the two different stomatal conductance models seen in Fig. 15 suggests that the stomatal conductance scheme has a greater impact on the model's behaviour than the choice of soil water uptake response function.

4.2 Why does C₄ productivity increase so much in LSM simulations?

The results presented in Sect. 3.4.1 show that predictions of PFT composition, ecosystem productivity, respiration, and carbon dioxide exchange vary between LSM simulations using different water uptake response functions and stomatal conductance schemes and with respect to standard LPJ-GUESS. Very notably, both the gross and net productivities of C₄ grasses are substantially enhanced in the LSM simulations compared with the non-LSM model. This results in unrealistically high simulated LAI values at three out of the four sites where the grasses are allowed to grow without competition.

We found that the main reason for this behaviour is the occurrence of higher photosynthetic rates in the LSM simulations due to the mitigation of biochemical N limitation at higher leaf temperatures. Standard LPJ-GUESS uses a daily average of the forcing air temperature as a proxy for leaf temperature in the V_max calculation. By contrast, LPJ-GUESS/LSM simulates leaf temperature explicitly and uses a daytime average (Eq. 42) to estimate V_max. This average leaf temperature can be several degrees above the forcing air temperature, which makes it possible to reach the optimal maximum carboxylation rate at lower leaf nitrogen concentrations (see Haxeltine and Prentice, 1996). This makes it easier for the plants to attain the optimal V_max with the available nitrogen, which enhances productivity. Exceedingly high leaf temperatures can have a negative impact on V_max due to the thermal breakdown of the biochemical reactions. However, the simulated leaf temperatures are still within the optimal temperature range for C₄ grasses (20 to 45 ^∘C in LPJ-GUESS). The temperature dependence of V_max, including the effect on nitrogen limitation, is illustrated in Fig. 16.

At AU-Stp, all three simulations predict much lower productivities. In this case, water availability is the limiting factor. This site receives, on average, considerably less rainwater than the other three C₄ grasslands (Table 10), which leads to lower values of the β factor (Eq. 46) and brings photosynthetic rates down.

As pointed out in Sect. 3.2, the simulated PFTs were restricted to grassy types at these sites. Table 11 shows a summary of LAI values predicted by standard LPJ-GUESS and two representative LSM simulations when establishment is not restricted to grassy PFTs. All three experiments predict a mixture of trees and grasses, but total LAI in the LSM runs is much lower than in the simulations where only grasses were allowed to establish. In these runs, competition with trees limits the resources available to grasses, and shading from the taller trees helps lower the average leaf temperature of the grassy understory, all of which helps counteract the effect described above.

4.3 Conclusion and outlook

The developments presented in this paper will enable us to study feedbacks between the climate and the biosphere using the state-of-the-art DGVM LPJ-GUESS directly coupled to an atmospheric model. Work is in progress regarding the development of a flexible interface to enable such coupling and extending the model's ability to simulate cold-climate ecosystems. More work is also needed to characterize and fully understand the model's response to the switch from using air temperature as a proxy for leaf temperature to simulating leaf temperature explicitly, particularly as these concern the productivity of C₄ plants in well-watered, no-competition situations (e.g. monoculture crops or managed pastures) and the way the new schemes affect the simulation of the carbon cycle on regional and global scales. These developments will allow us to use LPJ-GUESS/LSM in regional as well as global studies. Given the capacity of LPJ-GUESS to represent land use change and management (Lindeskog et al., 2013, 2021; Olin et al., 2015), the range of applications includes exploring impacts of management on regional climate, which can be an important tool to help devise and assess climate change mitigation policies.

Appendix A: Derivation of the canopy conductance for water vapour flux

During a given time step Δt, the total amount of water evapotranspired from the sunlit part of the canopy can be expressed as the sum of the contributions from the dry and wet parts:

$$\begin{array}{ll}{\displaystyle }{E}_{\text{sun}}\mathrm{\Delta }t& {\displaystyle }=(\mathrm{1}-f_{\text{wet}})E_{\text{sun,tr}}\mathrm{\Delta }t+f_{\text{wet}}[E_{\text{sun,ev}}\mathrm{\Delta }{t}_{\text{wet,sun}}\\ \text{(A1)}& {\displaystyle }& {\displaystyle }+E_{\text{sun,tr}}(\mathrm{\Delta }t-\mathrm{\Delta }{t}_{\text{wet,sun}})],\end{array}$$

where E_sun is the actual evapotranspiration rate, E_sun,tr is the potential rate of transpiration, E_sun,ev is the potential rate of evaporation, and Δt_wet,sun is the time that the wet part of the sunlit canopy remains wet at the potential evaporation rate. The latter is calculated as

$$\begin{array}{}\text{(A2)}& \mathrm{\Delta }{t}_{\text{wet,sun}}=min\left({\displaystyle \frac{w_{\mathrm{c}}{\text{PAI}}_{\text{c,sun}}/{\text{PAI}}_{\mathrm{c}}}{f_{\text{wet}}{E}_{\text{sun,ev}}}},\mathrm{\Delta }t\right),\end{array}$$

where w_c is the current canopy water content (kg m⁻²).

The evaporation rates in the above equations can be expressed as follows:

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }{E}_{\text{sun}}=-\mathit{\rho }{g}_{\text{w,sun}}[q_{\text{ca}}-q^{\ast }(T_{\text{sun}}\left)\right],\text{(A4)}& {\displaystyle }& {\displaystyle }{E}_{\text{tr,sun}}=-\mathit{\rho }\sum _i{\text{LAI}}_{\text{sun}}^{\left(i\right)}{\displaystyle \frac{g_{\text{s,sun}}^{\left(i\right)}{g}_{\mathrm{b}}}{g_{\text{s,sun}}^{\left(i\right)}+g_{\mathrm{b}}}}[q_{\text{ca}}-q^{\ast }(T_{\text{sun}}\left)\right],\text{(A5)}& {\displaystyle }& {\displaystyle }{E}_{\text{ev,sun}}=-\mathit{\rho }{\text{PAI}}_{\text{c,sun}}{g}_{\mathrm{b}}[q_{\text{ca}}-q^{\ast }(T_{\text{sun}}\left)\right],\end{array}$$

where the index i runs over cohorts. Inserting these expressions into Eq. (A1), dividing both sides by Δt, simplifying, and rearranging terms yields

$$\begin{array}{ll}{\displaystyle }{g}_{\text{w,sun}}& {\displaystyle }=f_{\text{wet}}{\mathit{\eta }}_{\text{sun}}{\text{PAI}}_{\text{c,sun}}{g}_{\mathrm{b}}\\ \text{(A6)}& {\displaystyle }& {\displaystyle }+(\mathrm{1}-f_{\text{wet}}{\mathit{\eta }}_{\text{sun}})\sum _i{\text{LAI}}_{\text{sun}}^{\left(i\right)}{\displaystyle \frac{g_{\text{s,sun}}^{\left(i\right)}{g}_{\mathrm{b}}}{g_{\text{s,sun}}^{\left(i\right)}+g_{\mathrm{b}}}},\end{array}$$

where ${\mathit{\eta }}_{\text{sun}}=\mathrm{\Delta }{t}_{\text{wet,sun}}/\mathrm{\Delta }t$. Identical equations apply to the shaded part of the canopy.

Appendix B: List of symbols, parameters, and variables used in the model description

Table B1 Ecosystem structure.

Table B2 Energy balance.

Table B3 Radiative transfer.

Table B4 Assimilation and stomatal conductance.

Table B5 Soil physics.