a next-generation land surface model with high ecological realism

. Land biosphere processes are of central importance to the climate system. Speciﬁcally, ecosystems inter-act with the atmosphere through a variety of feedback loops that modulate energy, water, and CO 2 ﬂuxes between the land surface and the atmosphere across a wide range of temporal and spatial scales. Human land use and land cover modiﬁ-cation add a further level of complexity to land–atmosphere interactions. Dynamic global vegetation models (DGVMs) attempt to capture land ecosystem processes and are increas-ingly 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 modiﬁcations 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 modiﬁcations allow the model (LPJ-GUESS/LSM) to simulate the diurnal exchange of energy, water, and CO 2 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 ﬂuxes. Differences in predicted ecosystem function between standard LPJ-GUESS and LPJ-GUESS/LSM vary across land cover types. We ﬁnd that the emerging ecosystem composition and carbon ﬂuxes 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 accu-rate representation of ecosystem processes is essential.


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 . 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 D. Martín Belda et al.: LPJ-GUESS/LSM 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 2 ) concentration promotes vegetation growth through CO 2 fertilization, which increases plant CO 2 absorption from the atmosphere. However, higher temperatures caused by a higher atmospheric CO 2 concentration enhance the release of CO 2 from respiration 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 2 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(Friedlingstein et al., , 2014Jones 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 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  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 biosphereatmosphere regional (Wramneby et al., 2010;Smith et al., 2011;Zhang et al., 2014Zhang et al., , 2018Wu et al., 2016Wu et al., , 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 sitebased 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 ) 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  A. . 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 2 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 . 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 . 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(Lindeskog et al., , 2021Olin 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. 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.

Model modifications
Radiative transfer in standard LPJ-GUESS is based on Beer's law 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. (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: 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 2 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.

Energy balance
The energy balance of the patch canopy is described by the following equations (e.g. Bonan, 2015): 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; λ = 2.44×10 6 J kg −1 ). 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 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 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): 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 LAI (i) sun is the leaf area index of the sunlit part of cohort i. The stomatal conductance of the sunlit leaves of cohort i, g (i) s,sun , is related to the leaflevel 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. 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 0 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 (i) h,sun [sha] is the conductance for sensible heat transport from the sunlit (shaded) part of cohort i to the canopy air. g (i) w,sun [sha] 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.
The energy balance equation for the ground surface is where G is heat conducted into the ground. The sensible heat from the ground surface to the canopy air space is 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 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 where κ (1) s , and z (1) 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: where H ↑ and λE ↑ are, respectively, the sensible and latent heat fluxes into the atmosphere, given by λE ↑ = −ρλg aa (q atm − q ca ).
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 0 , 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 (1994Raupach ( , 1995: 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).

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) Saeki, 1953, 2005): where I ↓ D0 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 (Î The direct-beam radiation absorbed in a canopy layer l between P and P + l P is calculated as the fraction of the decrease in direct-beam intensity in that layer that is not scattered: 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: where I ↓ d0 is incoming atmospheric diffuse radiation. The radiation absorbed by the sunlit and shaded parts of this layer is where the sunlit and shaded fractions of the layer are given by The total amount of short-wave radiation absorbed by the sunlit and shaded parts of the canopy is obtained by summing over layers: The short-wave radiation absorbed by the ground surface is calculated as the difference between the downward and upward beams at P = PAI c : The short-wave radiation reflected back at the atmosphere is obtained by evaluating the upward beams at P = 0: 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: α leaf, vis = 0.1; α stem, vis = 0.16, (29) τ leaf, vis = 0.05; τ stem, vis = 0.001, (30) α leaf, nir = 0.45; α stem, nir = 0.39, (31) τ leaf, nir = 0.25; τ stem, nir = 0.001, 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: The sunlit leaf and plant area indices of cohort i are obtained by aggregating over layers: The sunlit plant area index for the whole canopy is calculated by summing over cohorts: Equations analogous to Eqs. (33) to (36) apply to the shaded parts of the canopy.

Long-wave radiative transfer
The long-wave radiation emitted by the sunlit part of the canopy is (Dai et al., 2004) where σ is the Stefan-Boltzmann constant, L ↓ is the incoming atmospheric long-wave radiation, T sun , T g is expressed in Kelvin, and 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 The bulk long-wave radiation emitted by the land surface toward the atmosphere is

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. (1991Collatz et al. ( , 1992, the strongoptimality 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 (i) max , is recalculated at the end of every simulation day and depends linearly on the total amount of daily absorbed photosynthetic active radiation, PAR (i) day (Haxeltine and Prentice, 1996): In this equation, V (i) max,day 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 (i) max, day 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 (i) max,day 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 (i) leaf,dt , is weighted by the daily-averaged fractions of sunlit and shaded leaves for cohort i: 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 where PAR sun,day and PAR sha,day 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, The maximum carboxylation rate per unit leaf area is then calculated as sun,dt 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 −1 converts the rate from d −1 to s −1 . 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 PFTspecific vertical rooting profile: 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 av decreases linearly in each soil layer with volumetric water content θ (j ) down to the wilting point: 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 ψ wilt = −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): 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 where the parameter c 2 depends on PFT, and takes values between 4.36 and 6.37 (Xue et al., 1991). In this study, we set c 2 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 (j ) av 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 2 concentration inside the stomatal cavity. This concentration is related to the atmospheric CO 2 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 2 concentration at the leaf surface, c s . The stomatal conductance for sunlit leaves of cohort i is where A (i) n,sun is the net photosynthetic rate per unit leaf area, g min is a minimum stomatal conductance, and g 1 is a PFTspecific parameter. The Medlyn model (Medlyn et al., 2011) is derived from the assumption that stomata optimize CO 2 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 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 .

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

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 (J m −2 s −1 ), is calculated as 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 where t is the time step in seconds, and c (j ) h , z (j ) and T (j ) s 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 ↓ + L ↓ ). 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 %).
The water conservation error is computed as 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, ∼ 5 × 10 −12 % of the precipitation at any given time step (σ ∼ 9 × 10 −12 ), and never surpasses 3 × 10 −8 % of the water input. We therefore conclude that the magnitude of the errors in energy balance closure and water conservation are negligible.

Evaluation set-up
We evaluated the revised model by comparing hourly and monthly simulated fluxes of sensible and latent heat and annual CO 2 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. 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 2 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.
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 4 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 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. the whole measurement period was repeated cyclically, with interannual trends in air temperature removed, and the atmospheric CO 2 concentration was kept at the level of the first year of observations at each site.

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: , 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 * < 0.2 m s −1 in order to avoid possibly biased eddy covariance measurements during periods of weak turbulence (Schroder et al., 2014).
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 (σ m /σ o ), which provides a measure of the agreement between observed and simulated variability.

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 4 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 3 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.
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 g C −2 yr −1 m −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 ∼ −62 g C −2 yr −1 m −2 . This result is inconsistent with measurements at the site, which indicate a carbon source of 143 g C −2 yr −1 m −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 4 grassland, savanna, and woody savanna sites. This results in an increased average NPP value in C 4 grasslands (∼ 18 % in the CLM/BB run and ∼ 31 % in the CLM/Med run) and a decreased average NPP at woody savanna sites (∼ −11 % and ∼ −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 ∼ −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 4 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 ∼ −129 % for the CLM/BB run and ∼ 122 % for the CLM/Med run. 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.  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).

Site
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 4 grasslands, savanna, woody sa-vanna, 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 4 sites simulations predict NEE values between −88 and −110 g C m −2 yr −1 , while observations indicate a less negative value of −33 g C m −2 yr −1 . For savanna, measured NEE is −221 g C m −2 yr −1 , while simulations predict an average between −48 and −56 g C m −2 yr −1 . For woody savanna, measured NEE averages to −238 g C m −2 yr −1 , while simulated fluxes range between −28 and 3 g C m −2 yr −1 . Simulated fluxes at evergreen broadleaf forests are, on average, between −84 and −130 g C m −2 yr −1 , while measurements indicate an average NEE of −396 g C m −2 yr −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 2 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 2 concentration. 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).

Annual and diurnal cycles of turbulent heat fluxes
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 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 2 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.
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 −2 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 −2 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 −2 (observed: ∼ 500 W m −2 ), 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 −2 above measured values from March to May and ∼ 30-45 W m −2 below the measurements between August and October. The average simulated latent heat diurnal cycle peaks at ∼ 350 W m −2 in May (observed: ∼ 175 W m −2 ) and at ∼ 25 W m −2 in September (observed: ∼ 145 W m −2 ; 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 −2 from September to October, while measured values stay at ∼ 60 W m −2 throughout the dry season. The average sensible heat diurnal cycle peaks at ∼ 530 W m −2 in September, while the observed average diurnal peak is slightly under ∼ 300 W m −2 (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 −2 for most of the year and increases to ∼ 30 W m −2 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 −2 . Sensible and latent heat fluxes are systematically overestimated by the model by ∼ 10-20 W m −2 . 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 −2 (measured: ∼ 100-150 W m −2 ). Latent heat flux peaks daily between ∼ 280 and 340 W m −2 (measured: ∼ 220-320 W m −2 , Fig. 12a-f).
At the GF-Guy site, another tropical rainforest, monthly averages of sensible heat are overestimated by ∼ 20 W m −2 throughout the year, while latent heat flux is underestimated by about the same amount. The simulated sensible heat diur- 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. nal cycle peaks, on average, by ∼ 100 W m −2 above the measured values, while the peak of the simulated latent heat diurnal cycle is ∼ 130 W m −2 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 −2 (measured: ∼ 190 W m −2 ), while the average latent heat diurnal cycle peaks at ∼ 180 W m −2 (measured: ∼ 360 W m −2 ).  (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 2 = RMSE 2 − Bias 2 ), the correlation coefficient, and the standard deviation of observed and modelled data: E 2 = σ 2 o + σ 2 m − 2σ o σ 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 /σ 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 3 grasslands (C3G), C 4 grasslands (C4G), evergreen broadleaf forest (EBF), and deciduous broadleaf forest (DBF).

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.
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 σ m /σ o ∼ 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.
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 3 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.

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. 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.  (7)    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 −2 , while the average value for the Med runs is ∼ 35 W m −2 . Average errors are also smaller in the Med simulations. For hourly fluxes, the average RMSE is ∼ 90 W m −2 for BB runs and ∼ 81 W m −2 for Med runs. For monthly fluxes, RMSE averages are ∼ 31 and ∼ 28 W m −2 , 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 −2 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 ∼ −10 W m −2 (BB: ∼ −20 W m −2 ) for hourly fluxes and ∼ −5 W m −2 (BB: ∼ −11 W m −2 ) 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 ∼ −2 W m −2 , 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 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.
standard deviation of the sample of modelled monthly fluxes is, on average, ∼ 1.7 times larger than observed.

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 ecosystematmosphere carbon fluxes.

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 * < 0.2 m s −1 ) 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.
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 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 −2 . 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) 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;V max : normalized maximum carboxylation rate without nitrogen limitation;V 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.
of sensible heat in the simulations. A detailed discussion of these issues is provided in Haughton et al. (2016). 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.

Why does C 4 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 4 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 4 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 4 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.

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 4 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(Lindeskog et al., , 2021Olin et al., 2015), the range of applications includes exploring impacts of management on regional cli- 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: 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 where w c is the current canopy water content (kg m −2 ). The evaporation rates in the above equations can be expressed as follows: where the index i runs over cohorts. Inserting these expressions into Eq. (A1), dividing both sides by t, simplifying, and rearranging terms yields g w,sun = f wet η sun PAI c,sun g b where η sun = t wet,sun / t. Identical equations apply to the shaded part of the canopy.  Plant matter area index of patch canopy layer l m 2 m −2 LAI (i) Leaf area index of cohort i in the patch m 2 m −2 SAI (i) Stem area index of cohort i in the patch m 2 m −2 LAI (i,l) Leaf area index of cohort i in patch canopy layer l m 2 m −2 SAI (i,l) Stem area index of cohort i in patch canopy layer l m 2 m −2 h c Patch canopy height m Net long-wave radiation emitted by the (sunlit, shaded) part of the canopy W m −2 H [sun,sha] Sensible heat flux from the (sunlit, shaded) canopy to the canopy air space W m −2 λE [sun,sha] Sensible heat flux from the (sunlit, shaded) canopy to the canopy air space W m −2 S g Short-wave radiation absorbed by the ground surface W m −2 L g Net long-wave radiation emitted by the ground surface W m −2 H g Sensible heat flux from the ground surface to the canopy air space W m −2 λE g Latent heat flux from the ground surface to the canopy air space W m −2 G Heat flux conducted into the soil W m −2 H ↑ Sensible heat flux from the canopy air space to the atmosphere W m −2 λE ↑ Sensible heat flux from the canopy air space to the atmosphere W m −2 T [sun,sha] Temperature of the (sunlit, shaded) part of the canopy • C, K T g Temperature of the ground surface • C, K T ca Temperature of the canopy air • C T atm Temperature of the atmosphere at the reference level • C q [sun,sha] Specific humidity of the stomatal cavity air for (sunlit, shaded) leaves kg kg −1 q * (T g ) Saturated specific humidity at the ground surface temperature kg kg −1 α Ground surface specific humidity as a fraction of q * (T g ) q ca Specific humidity of the canopy air kg kg −1 q atm Specific humidity of the atmosphere at the reference level kg kg −1 g b Leaf boundary layer conductance m s −1 g w, [sun,sha] Conductance to water vapour between the (sunlit, shaded) canopy and the canopy air m s −1 g surf Conductance to water vapour from near-surface soil pores to the ground surface m s −1 g ab Aerodynamic conductance from the ground surface to the canopy air m s −1 g aa Aerodynamic conductance from the canopy air to the atmosphere reference level m s −1 f wet Wet fraction of the canopy η [sun,sha] Factor limiting evaporation from the (sunlit, shaded) canopy z 0 Canopy roughness length m z d Zero plane displacement height m z (1) Top soil layer thickness m T  Photosynthetically active radiation absorbed by the (sunlit, shaded) part of cohort i W m −2 α [sun,sha], [vis,nir] Reflectivity of (sunlit, shaded) leaves in the (visible, infrared) waveband τ [sun,sha], [vis,nir] Transmissivity of (sunlit, shaded) leaves in the (visible, infrared) waveband -L ↓ Incoming (atmospheric) long-wave radiation W m 2 L ↑ Outgoing long-wave radiation W m 2 γ [sun,sha] Effective thermal emissivity of the (sunlit, shaded) part of the canopy - Table B4. Assimilation and stomatal conductance.

Parameter
Description Units g (i) s, [sun,sha] Stomatal conductance of the (sunlit, shaded) leaves of cohort i m s −1 g min Minimum stomatal conductance m s −1 g 1,BB Stomatal conductance parameter for the Ball-Berry model g 1,Med Stomatal conductance parameter for the Medlyn model kPa 0.5 A (i) n, [sun,sha] Net photosynthetic assimilation rate of (sunlit, shaded) leaves of cohort i µmol C m −2 s −1 h s, [sun,sha] Fractional humidity at the surface of (sunlit, shaded) leaves kPa kPa −1 D s, [sun,sha] Water vapour deficit at the surface of (sunlit, shaded) leaves kPa c s, [sun,sha] Carbon dioxide concentration at the surface of (sunlit, shaded) leaves µmol mol −1 PAR Water stress factor limiting the assimilation rate r (j ) Fraction of roots in soil layer j -W (j ) av Soil water uptake response function (layer j ) - Volumetric water content of soil layer l m 3 m −3 θ wilt Volumetric soil water content at wilting point m 3 m −3 θ fc Volumetric soil water content at field capacity m 3 m −3 ψ (l) Matric potential of soil layer l m ψ wilt Matric potential of soil water at wilting point m ψ fc Matric potential of soil water at saturation point m γ w Hydraulic conductivity m s −1 λ w Hydraulic diffusivity m 2 s −1 S θ Volumetric soil water uptake sink term m 3 m −3 Code and data availability. LPJ-GUESS is a worldwide developed and refined DGVM. The model code is managed and maintained by the Department of Physical Geography and Ecosystem Science, Lund University, Sweden. The source code can be made available with a collaboration agreement under the acceptance of certain conditions. For this reason, a DOI for the model code cannot be provided. The code with the augmentations developed for this paper is available to the editor and reviewers via a restricted link, on the condition that the code is used only for review purposes and is deleted after the review process. Additional details and information can be found at the LPJ-GUESS website (http://web.nateko.lu.se/lpj-guess; Smith et al., 2014). The forcing data, evaluation data, model output, and analysis scripts used in this study have been uploaded to a public repository with DOI https://doi.org/10.5281/zenodo.6856036 (Martin-Belda, 2021).
Author contributions. AA had the original idea and motivated the development. DMB designed the model augmentations described in this work with input from all the coauthors. DMB implemented the code, ran the experiments, and performed the model evaluation analysis. All the coauthors provided input regarding the analysis of the results and the discussion and helped shape the final form of the manuscript.
Competing interests. The contact author has declared that none of the authors has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Financial support. This research has been supported by the Deutsche Forschungsgemeinschaft (grant no. SCHM 2376/2-1) and the Helmholtz Gemeinschaft (grant no. W3 -Impuls-und Vernetzungsfond). David Wårlind, Stefan Olin, Jing Tang, and Benjamin Smith were financially supported by the strategic research area "Modeling the Regional and Global Earth System" (MERGE). David Wårlind, Stefan Olin, and Benjamin Smith were financially supported by the Swedish national strategic e-science research programme eSSENCE.
The article processing charges for this open-access publication were covered by the Karlsruhe Institute of Technology (KIT).
Review statement. This paper was edited by David Lawrence and reviewed by two anonymous referees.