The biophysics, ecology, and biogeochemistry of functionally diverse, vertically and horizontally heterogeneous ecosystems: the Ecosystem Demography model, version 2.2 – Part 1: Model description

Earth system models (ESMs) have been developed to represent the role of terrestrial ecosystems on the energy, water, and carbon cycles. However, many ESMs still lack representation of within-ecosystem heterogeneity and diversity. In this paper, we present the Ecosystem Demography model version 2.2 (ED-2.2). In ED-2.2, the biophysical and physiological processes account for the horizontal and vertical heterogeneity of the ecosystem: the energy, water, and carbon cycles are solved separately for a series of vegetation cohorts (groups of individual plants of similar size and plant functional type) distributed across a series of spatially implicit patches (representing collections of micro-environments that have a similar disturbance history). We define the equations that describe the energy, water, and carbon cycles in terms of total energy, water, and carbon, which simplifies the differential equations and guarantees excellent conservation of these quantities in long-term simulation (< 0.1 % error over 50 years). We also show examples of ED-2.2 simulation results at single sites and across tropical South America. These results demonstrate the model's ability to characterize the variability of ecosystem structure, composition, and functioning both at stand and continental scales. A detailed model evaluation was conducted and is presented in a companion paper (Longo et al., 2019a). Finally, we highlight some of the ongoing model developments designed to improve the model's accuracy and performance and to include processes hitherto not represented in the model.

1 Introduction

The dynamics of the terrestrial biosphere play an integral role in the Earth's carbon, water, and energy cycles (Betts and Silva Dias, 2010; Santanello Jr et al., 2018; Le Quéré et al., 2018), and consequently, how the Earth's climate system is expected to change over the coming decades due to the increasing levels of atmospheric carbon dioxide arising from anthropogenic activities (IPCC, 2014; Le Quéré et al., 2018). Models for the dynamics of the terrestrial biosphere and its bidirectional interaction with the atmosphere have evolved considerably over the past decades (Levis, 2010; Fisher et al., 2014, 2018). As described by Sellers et al. (1997), the first generation of land surface models (LSMs) was limited to provide boundary conditions to atmospheric models; they only solved a simplified energy and water budget, and accounted for the effects of surface on frictional effects on near-surface winds (e.g., Manabe et al., 1965; Somerville et al., 1974). These models, however, did not account for the active role of vegetation. The second generation of LSMs considered the active role of vegetation and represented the spectral properties of the canopy, the changes in roughness of vegetated surfaces, and the biophysical controls on evaporation and transpiration (Sellers et al., 1997); examples of these models include National Center for Atmospheric Research Biosphere-Atmosphere Transfer Scheme (NCAR/BATS) (Dickinson et al., 1986) and Simple Biosphere model (SiB) (Sellers et al., 1986). The increasing recognition of the role of vegetation in mediating the exchanges of carbon, water, and energy between the land and the atmosphere led to the third generation of LSMs, which incorporated explicit representations of plant photosynthesis and resulting dynamics of terrestrial carbon uptake, turnover, and release within terrestrial ecosystems (Sellers et al., 1997); examples of such models included LSM (Bonan, 1995) and SiB2 (Sellers et al., 1996). While the fluxes of carbon, water, and energy predicted by these models would change in response to changes in their climate forcing, the biophysical and biogeochemical properties of the ecosystem within each climatological grid cell were prescribed and thus did not change over time.

Subsequently, building upon previous work (Prentice et al., 1992; Neilson, 1995; Haxeltine and Prentice, 1996), Foley et al. (1996) adopted an approach to calculate the productivity of a series of plant functional types (PFTs), based on a leaf-level model of photosynthesis. The abundance of each PFT within each grid cell was dynamic, with the abundance changes being determined by the relative productivity of the PFTs. This allowed the fast-timescale exchanges of carbon, water, and energy within the plant canopy to be explicitly linked with the long-term dynamics of the ecosystem. This approach followed the concept of dynamic global vegetation model (DGVM), originally coined by Prentice et al. (1989) to describe this kind of terrestrial biosphere model in which changes in climate could drive changes in ecosystem composition, structure, and functioning. DGVMs, when run coupled to atmospheric models, would then feedback onto climate. The subsequent generation of terrestrial biosphere-based DGVMs (i.e., DGVMs incorporating couple carbon, water, and energy fluxes) such as the Lund–Potsdam–Jena (LPJ) model (Sitch et al., 2003), Community Land Model's Dynamic Global Vegetation Model (CLM-DGVM) (Levis et al., 2004), and Top-down Representation of Interactive Foliage and Flora Including Dynamics Joint UK Land Environment Simulator (TRIFFID/JULES) (Hughes et al., 2004; Clark et al., 2011; Mangeon et al., 2016) have included additional mechanisms such as disturbance through fires and multiple types of mortality.

Analyses have shown that most terrestrial biosphere models are capable of reproducing the current distribution of global biomes (e.g., Sitch et al., 2003; Blyth et al., 2011) and their carbon stocks and fluxes (Piao et al., 2013). However, they diverge markedly in their predictions of how terrestrial ecosystems will respond to future climate change (Friedlingstein et al., 2014). In fully coupled Earth system model simulations, some of these differing predictions arise from divergent predictions about the direction and magnitude of regional climate change. However, offline analyses, in which the models are forced with prescribed climatological forcing, have shown that there is also substantial disagreement between the models about how terrestrial ecosystems will respond to any shift in climate (e.g., Sitch et al., 2008; Zhang et al., 2015). In addition, the transitions between biome types, for example, the transition that occurs between closed-canopy tropical forests and grass- and shrub-dominated savannahs in South America, are generally far more abrupt in typical DGVM results than in observations (Good et al., 2011; Levine et al., 2016).

One important limitation of most DGVMs is that they do not represent within-ecosystem diversity and heterogeneity. The representation of plant functional diversity within terrestrial biosphere models is normally coarse, with broadly defined PFTs defined from a combination of morphological and leaf physiological attributes (Purves and Pacala, 2008). In addition, there is limited variation in the resource conditions (light, water, and nutrient levels) experienced by individual plants within the climatological grid cells of traditional DGVMs. Some models, such as CLM (Oleson et al., 2013), have options to represent a multi-layer plant canopy (e.g., two canopy layers allowing for Sun and shade leaves), and/or differences in rooting depth between PFTs; however, resource conditions are assumed to be horizontally homogeneous, meaning that there is no horizontal spatial variation in resource conditions experienced by individual plants. The lack of significant variability in resource conditions limits the range of environmental niches within the climatological grid cells of terrestrial biosphere and makes the coexistence between PFTs difficult. Consequently, models often predict ecosystems comprised of single homogeneous vegetation types (Moorcroft, 2003, 2006).

Field- and laboratory-based studies conducted over the past 30 years indicate that plant functional diversity significantly affects ecosystem functioning (Loreau and Hector, 2001; Tilman et al., 2014, and references therein), and variations in trait expression are strongly driven by disturbances and local heterogeneity of abiotic factors such as soil characteristics (Bruelheide et al., 2018; Both et al., 2019). In many cases, biodiversity increases ecosystem productivity and ecosystem stability (e.g., Tilman and Downing, 1994; Naeem and Li, 1997; Cardinale et al., 2007; García-Palacios et al., 2018), and biodiversity has also been shown to contribute to enhanced ecosystem functionality in highly stressed environments (e.g., Jucker and Coomes, 2012). Other studies have also established correlations between tropical forest diversity and carbon storage and primary productivity (Cavanaugh et al., 2014; Poorter et al., 2015; Liang et al., 2016; Huang et al., 2018).

In addition to the absence of within-ecosystem diversity in conventional terrestrial biosphere models, plants of each PFT are also assumed to be homogeneous in size, while, in contrast, most terrestrial ecosystems, particularly forests and woodlands, exhibit marked size structure of individuals within plant canopies (Hutchings, 1997). This size-related heterogeneity is important because plant size strongly affects the amount of light, water, and nutrients individual plants within the canopy can access, which, in turn, affects their performance, dynamics, and responses to climatological stress. It also allows representation of the dynamics of pervasive human-driven degradation of forest ecosystems (Lewis et al., 2015; Haddad et al., 2015), which affects carbon stocks and forest structure and composition, which cannot be easily represented in highly aggregated models (Longo and Keller, 2019).

An alternative approach to simulating the dynamics of terrestrial ecosystems has been individual-based vegetation models (Friend et al., 1997; Bugmann, 2001; Sato et al., 2007; Fischer et al., 2016; Maréchaux and Chave, 2017). Also known as forest gap models, due to the importance of canopy gaps for the dynamics of closed canopy forests, these models simulate the birth, growth, and death of individual plants, thereby incorporating diversity and heterogeneity of the plant canopy mechanistically. In forest gap models, the ecosystem properties such as total carbon stocks, and net ecosystem productivity are emergent properties resulting from competition of limiting resources and the differential ability of plants to survive and be productive under a variety of micro-environments (e.g., gaps or the understory of a densely populated patch of old-growth forest).

This approach has two main advantages. First, gap models represent the dynamic changes in the ecosystem structure caused by disturbances such as tree fall, selective logging, and fires. These disturbances create new micro-environments that are significantly different from old-growth vegetation areas and allow plants with different life strategies (for example, shade-intolerants) to coexist in the landscape. Second, because individual trees are represented in the model, the results can be directly compared with field measurements. Gap models have various degrees of complexity, with some models being able to represent the interactions between climate variability and gross primary productivity (Friend et al., 1997; Sato et al., 2007), as well as the impact of climate change in the ecosystem carbon balance (Fischer et al., 2016, and references therein). However, because the birth and death of individuals within a plant canopy are stochastic processes, multiple realizations of given model formulation are required to determine the long-term, large-scale dynamics of these models, which limits their applicability over large regions or global scales and has precluded their use in Earth system modeling studies.

The need to represent heterogeneity in vegetation structure and composition in terrestrial biosphere models, without the computational burden of simulating every tree at regional and global scales, led to the development of cohort-based models (Hurtt et al., 1998; Moorcroft et al., 2001; Fisher et al., 2018). In the cohort-based approach, individual trees are grouped according to their size (e.g., height or diameter at breast height); functional groups, which can be defined along trait axes (e.g., Reich et al., 1997; Wright et al., 2004; Fortunel et al., 2012); and micro-environment conditions (e.g., whether plants are living in a gap, recently burned fragment, or in a patch of old-growth forest). Over the past two decades, several cohort-based models have emerged, including the Ecosystem Demography model (ED; Moorcroft et al., 2001; Hurtt et al., 2002; Albani et al., 2006; Medvigy et al., 2009); the Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS; Smith et al., 2001; Ahlström et al., 2012; Lindeskog et al., 2013); the Land Model version 3 with perfect plasticity approximation (LM3-PPA; Weng et al., 2015); and the Functionally-Assembled Terrestrial Ecosystem Simulator (FATES; Fisher et al., 2015; Huang et al., 2019). Similar to gap models, these models represent functional diversity and heterogeneity of micro-environments, and consequently the ecosystem's structure, diversity, and functioning emerge from the interactions between plants with different life strategies under different resource availability, albeit at a lesser extent than individual-based models (Fisher et al., 2018).

The Ecosystem Demography model (ED; Moorcroft et al., 2001) is a cohort-based model. Through this approach, it addresses the need to incorporate heterogeneity into models of the long-term, large-scale response of terrestrial ecosystems to changes in climate and other environmental forcings within a deterministic modeling framework. The size- and age-structured partial differential equations that describe the plant community are derived from individual-level properties but are properly scaled to account for the spatially localized nature of interactions within plant canopies. The model was later extended by Hurtt et al. (2002) and Albani et al. (2006) to incorporate multiple forms of disturbance including land clearing, land abandonment, and forest harvesting. An important difference between ED and most DGVMs is that in ED, PFTs are defined not simply based on their biogeographic ranges but also represent diversity in plant life-history strategies within any given ecosystem. These different PFTs represent a suite of physiological, morphological, and life-history traits that mechanistically represent the ways different kinds of plants utilize resources (Fisher et al., 2010).

The original ED model formulation was an offline ecosystem model describing the coupled carbon and water fluxes of a heterogeneous tropical forest ecosystem (Moorcroft et al., 2001). Subsequently, Medvigy et al. (2009) applied a similar approach to develop the Ecosystem Demography model version 2 (ED-2) that describes coupled carbon, water, and energy fluxes of the land surface. Since then, the ED-2 model has been continuously developed to improve several aspects of the model (see Supplement Sect. S1 for further information): (1) the conservation and thermodynamic representation of energy, water, and carbon cycles of the ecosystems; (2) the representation of several components of the energy, water, and carbon cycles, including the canopy radiative transfer, aerodynamic conductances and eddy fluxes, and leaf physiology (photosynthesis); (3) the structure of the code, including efficient data storage, code parallelization, and version control and code availability. ED-2 has been used in many studies including offline simulations (e.g., Medvigy et al., 2009; Antonarakis et al., 2011; Kim et al., 2012; Zhang et al., 2015; Castanho et al., 2016; Levine et al., 2016) and simulations running interactively with a regional atmospheric model (e.g., Knox et al., 2015; Swann et al., 2015).

In this paper, we describe in detail the biophysical, physiological, ecological, and biogeochemical formulation of the most recent version of the ED-2 model (ED-2.2), focusing in particular on the model's formulation of the fast timescale dynamics of the heterogeneous plant canopy that occur at subdaily timescales. While many parameterizations and submodels in ED-2.2 are based on approaches that are also used in other DGVMs, their implementation in ED-2.2 has some critical differences from other ecosystem models and also previous versions of ED. (1) In ED-2.2, the fundamental budget equations use energy and total mass as the main prognostic variables; because we use equations that directly track the time changes of the properties we seek to conserve, we can assess the model conservation of such properties with fewer assumptions. (2) In ED-2.2, all thermodynamic properties are scalable with mass, and the model is constructed such that when individual biomass changes, due to growth and turnover, the thermodynamic properties are also updated to reflect changes in heat and water holding capacity. (3) The water and energy budget equations for vegetation are solved at the individual cohort level and the corresponding equations for environments shared by plants such as soils and canopy air space are solved for each micro-environment in the landscape, and thus ecosystem-scale fluxes are emergent properties of the plant community. This approach allows the model to represent both the horizontal and vertical heterogeneity of environments of the plant communities. It also links the individual's ability to access resources such as light and water and accumulate carbon under a variety of micro-environments, which ultimately drives the long-term dynamics of growth, reproduction, and survivorship.

2 Model overview

2.1 The representation of ecosystem heterogeneity in ED-2.2

In ED-2.2 the terrestrial ecosystem within a given region of interest is represented through a hierarchy of structures to capture the physical and biological heterogeneity in the ecosystem's properties (Fig. 1).

Figure 1 Schematic representation of the multiple hierarchical levels in ED-2.2, organized by increasing level of detail from top to bottom. Static levels (grid, polygons, and sites) are assigned during the model initialization and remain constant throughout the simulation. Dynamic levels (patches and cohorts) may change during the simulation according to the dynamics of the ecosystem.

Physical heterogeneity. The domain of interest (grid) is geographically divided into polygons. Within each polygon, the time-varying meteorological forcing above the plant canopy is assumed to be spatially uniform. For example, a single polygon may be used to simulate the dynamics of an ecosystem in the neighborhood of an eddy flux tower, or alternatively, a polygon may represent the lower boundary condition within one horizontal grid cell in an atmospheric model. Each polygon is subdivided into one or more sites that are designed to represent landscape-scale variation in other abiotic properties, such as soil texture, soil depth, elevation, slope, aspect, and topographic moisture index. Each site is defined as a fractional area within the polygon and represents all regions within the polygon that share similar time-invariant physical (abiotic) properties. Both polygons and sites are defined at the beginning of the simulation and are fixed in time, and no geographic information exists below the level of the polygon.

Biotic heterogeneity. Within each site, horizontal, disturbance-related heterogeneity in the ecosystem at any given time t is characterized through a series of patches that are defined by the time elapsed since last disturbance (i.e., age, a) and the type of disturbance that generated them. Like sites, patches are not physically contiguous: each patch represents the collection of canopy gap-sized (∼10 m) areas within the site that have a similar disturbance history, defined in terms of the type of disturbance experienced (represented by subscript q, $q\in \mathrm{1},\mathrm{2},\mathrm{\dots },N_Q$; a list of indices is available in Table 1) and time since the disturbance event occurred. The disturbance types accounted for in ED-2.2, and the possible transitions between different disturbance types, are shown in Fig. S1. The collection of gaps within each given site belonging to a polygon follows a probability distribution function α, which can be also thought of as the relative area within a site, that satisfies

$$\begin{array}{}\text{(1)}& \sum _{q=\mathrm{1}}^{N_Q}\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\alpha }}_q\left(a,t\right)\mathrm{d}a\right]=\mathrm{1}.\end{array}$$

Similarly, the plant community population is characterized by the number of plants per unit area (hereafter number density, n) and is further classified according to their PFT, represented by subscript f ($f\in \mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{F}}$; Table 1) and the type of gap (q). The number density distribution depends on the individuals' biomass characteristics (size, C), the age since last disturbance (a), and the time (t), and is expressed as $n_{fq}(\mathit{C},a,t)$. Size is defined as a vector $\mathit{C}=n_{fq}^{-\mathrm{1}}\left(C_{\mathrm{l}};C_{\mathrm{r}};C_{\mathit{\sigma }};C_{\mathrm{h}};C_{\mathrm{n}}\right)$ (units: kg_C plant⁻¹) corresponding to biomass of leaves, fine roots, sapwood, heartwood, and non-structural storage (starch and sugars), respectively.

Table 1 List of subscripts associated with ED-2.2's hierarchical levels for any variable X. N_T is the total number of cohorts, N_G is the total number of soil (ground) layers, N_S is the total number of temporary surface water/snowpack layers, and N_C is the total number of canopy air space layers, currently only used to obtain properties related to canopy conductance. The complete list of subscripts is available in Table S1.

Following Moorcroft et al. (2001), Albani et al. (2006), and Medvigy and Moorcroft (2012), the partial differential equations that describe the dynamics of plant density and probability distribution of patches within each site in the size-and-age structured model are defined as (dependencies omitted in the equations for clarity)

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\underset{\text{Change rate}}{\underbrace{{\displaystyle \frac{\partial n_{fq}}{\partial t}}}}=\underset{\text{Aging}}{\underbrace{-{\displaystyle \frac{\partial n_{fq}}{\partial a}}}}\underset{\text{Growth}}{\underbrace{-{\mathbf{\nabla }}_{\mathit{C}}\cdot \left({\mathit{g}}_f{n}_{fq}\right)}}\underset{\text{Mortality}}{\underbrace{-m_f{n}_{fq}}},\text{(3)}& {\displaystyle }& {\displaystyle }-\underset{\text{Change rate}}{\underbrace{{\displaystyle \frac{\partial {\mathit{\alpha }}_q}{\partial t}}}}=\underset{\text{Aging}}{\underbrace{{\displaystyle \frac{\partial {\mathit{\alpha }}_q}{\partial a}}}}\underset{\text{Disturbance}}{\underbrace{-\sum _{q^{\prime }=\mathrm{1}}^{N_Q}\left({\mathit{\lambda }}_{q^{\prime }q}{\mathit{\alpha }}_q\right)}},\end{array}$$

where m_f is mortality rate, which may depend on the PFT, size, and the individual carbon balance; g_f is the vector of the net growth rates for each carbon pool, which also may depend on the PFT, size, and carbon balance; ∇_C⋅ is the divergence operator for the size vector; and ${\mathit{\lambda }}_{q^{\prime }q}$ is the transition matrix from gaps generated by previous disturbance q^′ affected by new disturbance of type q, which may depend on environmental conditions. Boundary conditions are shown in Sect. S2.

Equations (2) and (3) cannot be solved analytically except for the most trivial cases; therefore, the age distribution is discretized into patches (subscript u, $u\in \mathrm{1},\mathrm{2},\mathrm{\dots },N_P$; Table 1) of similar age and same disturbance type, and the population size structure living in any given patch u is discretized into cohorts (subscript k, $k\in \mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{T}}$; Table 1) of similar size and same PFT (Fig. 1). Unlike polygons and sites, patches and cohorts are dynamic levels: changes in distribution (fractional area) of patches are driven by aging and disturbance rates, whereas changes in the distribution of cohorts in each patch are driven by growth, mortality, and recruitment (Fig. S2).

The environment perceived by each plant (e.g., incident light, temperature, vapor pressure deficit) varies across large scales as a consequence of changes in climate (macro-environment) but also varies at small scales (within the landscape; micro-environment) because of the horizontal and vertical position of each individual relative to other individuals in the plant community (e.g., Bazzaz, 1979) and the position of the local community in landscapes with complex terrains. Both macro- and micro-environmental conditions drive the net primary productivity of each individual, and ultimately determine growth, mortality, and recruitment rates for each individual. Likewise, they can also affect the disturbance rates: for example, during drought conditions (macro-environment), open canopy patches (micro-environment) may experience faster ground desiccation and consequently increase local fire risk. To account for the variability in micro-environments within the landscape and within local plant communities, in ED-2.2 the energy, water, and carbon dioxide cycles are solved separately for each patch, and within each patch, fluxes and storage associated with individual plants are solved for each cohort.

The ED-2.2 model represents processes that have inherently different timescales; therefore, the model also has a hierarchy of time steps, in order to attain maximum computational efficiency (Table 2). Processes associated with the short-term dynamics are presented in this paper. A summary of the phenological processes and those associated with longer-term dynamics is presented in Sects. S3 and S4 (see also Moorcroft et al., 2001; Albani et al., 2006; Medvigy et al., 2009).

Table 2 Time steps associated with processes resolved by ED-2.2. The thermodynamic substep is dynamic and it depends on the error evaluation of the integrator, but it cannot be longer than the biophysics step, which is defined by the user. Other steps are fixed as of ED-2.2. Processes marked with an asterisk are presented in the main paper. Other processes are described in Sects. S3 and S4.

2.2 Software requirements and model architecture

Software requirements. The ED-2.2 source code is mainly written in Fortran 90, with a few file management routines written in C. Most input and output files use the Hierarchical Data Format 5 (HDF5) format and libraries (The HDF Group, 2016). In addition, the Message Passing Interface (MPI) is highly recommended for regional simulations and is required for simulations coupled with the Brazilian Improvements on Regional Atmospheric Modeling System (BRAMS) atmospheric model (Knox et al., 2015; Swann et al., 2015; Freitas et al., 2017). The source code can be also compiled with shared memory processing (SMP) libraries, which enable parallel processing of thermodynamics and biophysics steps at the patch level and thus allow shorter simulation time.

Code design and parallel structure. ED-2.2 has been designed to be run in three different configurations: (1) as a stand-alone land surface model over a small list of specified locations (sites); (2) as a stand-alone land surface model distributed over a regional grid; (3) coupled with an atmospheric model distributed over a regional grid (e.g., ED-BRAMS; Knox et al., 2015; Swann et al., 2015). For regional stand-alone grids, the model partitions the grid into spatially contiguous tiles of polygons, which access the initial and boundary conditions and are integrated independently of each other but write the results to a unified output file using collective input/output functions from HDF5. In the case of simulations dynamically coupled with an atmospheric model such as BRAMS, polygons are defined to match each atmospheric grid cell.

Memory allocation. The code uses dynamic allocation of variables and extensive use of pointers to efficiently reduce the amount of data transferred between routines. To reduce the output file size, polygon-, site-, patch-, and cohort-level variables are always written as long vectors, and auxiliary index vectors are used to map variables from higher hierarchical levels to lower hierarchical levels (for example, to which patch a cohort-level variable belongs).

2.3 Model inputs

Every ED-2.2 simulation requires an initial state for forest structure and composition (initial state), a description of soil characteristics (edaphic conditions), and a time-varying list of meteorological drivers (atmospheric conditions).

Initial state. To initialize a plant community from inventory data, one must have either the diameter at breast height of every individual or the stem density of different diameter size classes, along with plant functional type identification and location; in addition, necromass from the litter layer, woody debris, and soil organic carbon are needed. Alternatively, initial conditions can be obtained from airborne lidar measurements (Antonarakis et al., 2011, 2014) or a prescribed near-bare-ground condition may be used for long-term spin up simulations. Previous simulations can be used as initial conditions as well.

Edaphic conditions. The user must also provide soil characteristics such as total soil depth, total number of soil layers, the thickness of each layer, as well as soil texture, color, and the bottom soil boundary condition (bedrock, reduced drainage, free drainage, or permanent water table). This flexibility allows the user to easily adjust the soil characteristics according to their regions of interest. Soil texture can be read from standard data sets (e.g., Global Soil Data Task, 2000; Hengl et al., 2017) or provided directly by the user. Soil layers, soil color, and bottom boundary condition must be provided directly by the user as of ED-2.2. In addition, simulations with multiple sites per polygon also need to provide the fractional areas of each site and the mean soil texture class, soil depth, slope, aspect, elevation, and topographic moisture index of each site.

Atmospheric conditions. Meteorological conditions needed to drive ED-2.2 include temperature, specific humidity, CO₂ molar fraction, pressure of the air above the canopy, precipitation rate, incoming solar (shortwave) irradiance (radiation flux), and incoming thermal (longwave) irradiance (Table 3) at a reference height that is at least a few meters above the canopy. Subdaily measurements (0.5–6 h) are highly recommended so the model can properly simulate the diurnal cycle and interdiurnal variability. Meteorological drivers can be either at a single location (e.g., eddy covariance towers) or gridded meteorological drivers such as reanalyses (e.g., Dee et al., 2011; Gelaro et al., 2017) or bias-corrected products based on reanalysis (e.g., Sheffield et al., 2006; Weedon et al., 2014). Whenever available, CO₂ must be provided at comparable temporal and spatial resolution to other meteorological drivers; otherwise, it is possible to provide spatially homogeneous, time-variant CO₂, or constant CO₂, although this may increase uncertainties in the model predictions (e.g., Wang et al., 2007). Alternatively, the meteorological forcing (including CO₂) may be provided directly by BRAMS (Knox et al., 2015; Swann et al., 2015).

Table 3 Atmospheric boundary conditions driving the ED-2.2 model. Variable names and subscripts follow a standard notation throughout the paper (Tables S1 and S2). Flux variables between two thermodynamic systems are defined by a dot and two indices separated by a comma, and they are positive when the net flux goes from the thermodynamic system represented by the first index to the one represented by the second index.

Plant functional types. The user must specify which PFTs are allowed to occur in any given simulation. ED-2.2 has a list of default PFTs, with parameters described in Tables S5–S6. Alternatively, the user can modify the parameters of existing PFTs or define new PFTs through an extensible markup language (XML) file, which is read during the model initialization.

3 Overview of enthalpy, water, and carbon dioxide cycles

Here, we present the fundamental equations that describe the biogeophysical and biogeochemical cycles. Because the environmental conditions are a function of the local plant community and resources are shared by the individuals, these cycles must be described at the patch level, and the response of the plant community can be aggregated to the polygon level once the cycles are resolved for each patch. In ED-2.2, patches do not exchange enthalpy, water, and carbon dioxide with other patches; thus, patches are treated as independent systems. Throughout this section, we will only refer to the patch and cohort levels, and indices associated with patches, sites, and polygons will be omitted for clarity.

3.1 Definition of the thermodynamic state

Each patch is defined by a thermodynamic envelope (Fig. 2), comprised of multiple thermodynamic systems: each soil layer (total number of layers N_G), each temporary surface water or snow layer (total number of layers N_S), leaves, and branch wood portion of each cohort (total number of cohorts N_T), and the canopy air space. For simplicity, roots are assumed to be in thermal equilibrium with the soil layers and have negligible heat capacity compared to the soil layers. Although patches do not exchange heat and mass with other patches, they are allowed to exchange heat and mass with the free air (i.e., the atmosphere above and outside of the air space control volume we deem as within canopy) and lose water and associated energy through surface and subsurface runoff. We also assume that intensive variables such as pressure and temperature are uniform within each thermodynamic system. Note that free air is not considered a thermodynamic system in ED-2 because the thermodynamic state is determined directly from the boundary conditions and thus external to the model.

Table 4 List of state variables solved in ED-2.2. Unless otherwise noted, the reference equation is the ordinary differential equation that defines the rate of change of the thermodynamic state. The list of fluxes that describe the thermodynamic state is presented in Table 5. For a complete list of subscripts and variables used in this paper, refer to Tables S1–S2.

^a Budget fluxes are in units of area, and the state variable is updated following the conversion described in Sect. 3.2.1. ^b Canopy air space pressure is not solved using ordinary differential equations but based on the atmospheric pressure from the meteorological forcing. ^c Canopy air space depth is determined from vegetation characteristics, not from an ordinary differential equation.

The fundamental equations that describe the system thermodynamics are the first law of thermodynamics in terms of enthalpy H (J m⁻²), and the mass continuity for incompressible fluids for total water mass W (kg_W m⁻²):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \underset{\text{Change in enthalpy}}{\underbrace{{\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}}}=\underset{\text{Net heat flux}}{\underbrace{\dot{Q}}}+\underset{\begin{array}{c}\text{Enthalpy flux due}\\ \text{ to mass flux}\end{array}}{\underbrace{\dot{H}}}\\ & +\underset{\text{Pressure change}}{\underbrace{\mathcal{V}{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(5)}& \underset{\text{Change in water mass}}{\underbrace{{\displaystyle \frac{\mathrm{d}W}{\mathrm{d}t}}}}=\underset{\text{Net water mass flux}}{\underbrace{\dot{W}}},\end{array}$$

where 𝒱 is the volume of the thermodynamic system and p is the ambient pressure. The components on the right-hand side of Eqs. (4) and (5) depend on the thermodynamic system and will be presented in detail in the following sections. Net heat fluxes ($\dot{Q}$) represent changes in enthalpy that are not associated with mass exchange (radiative and sensible heat fluxes), whereas the remaining enthalpy fluxes ($\dot{H}$) correspond to changes in heat capacity due to addition or removal of mass from each thermodynamic system.

The merit of solving the changes in enthalpy over internal energy is that changes in enthalpy are equivalent to the net energy flux when pressure is constant (Eq. 4). Pressure is commonly included in atmospheric measurements, making it easy to track changes in enthalpy not related to energy fluxes. In reality, the only thermodynamic system where the distinction between internal energy and enthalpy matters is the canopy air space. Work associated with thermal expansion of solids and liquids is several orders of magnitude smaller than heat (Dufour and van Mieghem, 1975), and changes in pressure contribute significantly less to enthalpy because the specific volumes of solids and liquids are comparatively small. Likewise, enthalpy fluxes that do not involve gas phase (e.g., canopy dripping and runoff) are nearly indistinguishable from internal energy flux, whereas differences between enthalpy and internal energy fluxes are significant when gas phase is involved (e.g., transpiration and eddy flux). For simplicity, from this point on, we will use the term “enthalpy” whenever internal energy is indistinguishable from enthalpy. The complete list of state variables in ED-2.2 is shown in Table 4.

Variations in enthalpy are more important than their actual values, but they must be consistently defined relative to a pre-determined and known thermodynamic state, at which we define enthalpy to be zero. For any material other than water, enthalpy is defined as zero when the material temperature is 0 K; for water, enthalpy is defined as zero when water is at 0 K and completely frozen. The general definitions of enthalpy and internal energy states used in all thermodynamic systems in ED-2.2 are described in Sect. S5. In ED-2.2, enthalpy is used as the prognostic variable because these are directly and linearly related to the governing ordinary differential equation (Eq. 4). Temperature is diagnostically obtained based on the heat capacity of each thermodynamic system, and the heat capacities of different thermodynamic systems are defined in Sect. S6.

3.2 Heat ($\dot{Q}$), water ($\dot{W}$), and enthalpy ($\dot{H}$) fluxes

The enthalpy and water cycles for each patch in ED-2.2 are summarized in Fig. 2, and these cycles are solved every thermodynamic substep (Δt_Thermo), using a fourth-order Runge–Kutta integrator with dynamic time steps to maintain the error within prescribed tolerance. For all fluxes and variables, we follow the subscript notation described in Table S1, and denote flux variables with a dot and two indices separated by a comma, denoting the systems impacted by the flux. For any variable X that has flux between a system m and a system n, we assume that ${\dot{X}}_{m,n}>\mathrm{0}$ when the net flux goes from system m to system n, and that ${\dot{X}}_{m,n}=-{\dot{X}}_{n,m}$. Arrows in Fig. 2 indicate the directions allowed in ED-2.2. The list of fluxes solved in ED-2.2 is provided in Table 5, and a complete list of variables is provided in Table S2. In addition, the values of global constants and global parameters are listed in Tables S3 and S4, respectively, and the default parameters specific for each tropical plant functional type are presented in Tables S5–S6.

Figure 2 Schematic of the fluxes that are solved in ED-2.2 for a single patch (thermodynamic envelope). In this example, the patch has N_T cohorts, N_G soil layers, and N_S=1 temporary surface water. Both N_G and the maximum N_S are specified by the user; N_T is dynamically defined by ED-2.2. Letters near the arrows are the subscripts associated with fluxes, although the flux variables have been omitted here for clarity. Solid red arrows represent heat flux with no exchange of mass, and dashed yellow arrows represent exchange of mass and associated enthalpy. Arrows that point to a single direction represent fluxes that can only go in one (non-negative) direction, and arrows pointing to both directions represent fluxes that can be positive, negative, or zero.

Table 5 List of energy, water, and carbon dioxide fluxes that define the thermodynamic state in ED-2.2, along with sections and equations that define them. Fluxes are denoted by a dotted letter, and two subscripts separated with a comma: ${\dot{X}}_{m,n}$. Positive fluxes go from thermodynamic system m to thermodynamic system n; negative fluxes go in the opposite direction. Acronyms in the description column: canopy air space (CAS); temporary surface water (TSW). The complete list of subscripts and variables used in this paper is available in Tables S1–S2, and the list of state variables is shown in Table 4.

^a Net flux between leaf respiration (positive) and gross primary productivity (negative). ^b When negative, this flux corresponds to dew or frost formation.

3.2.1 Soil

In ED-2.2, the soil characteristics (number of soil layers, thickness of each soil layer and total soil depth, soil texture, soil color) are defined by the user, and assumed constant throughout the simulation. Within each patch, each soil layer (comprised by soil matrix and soil water in each layer) is considered a separate thermodynamic system, with the main size dimension being the layer thickness $\mathrm{\Delta }{z}_{{\text{g}}_j}$, with j=1 being the deepest soil layer, and j=N_G being the topmost soil layer. Typically, the top layer thickness is set to $\mathrm{\Delta }{z}_{g_{N_{\mathrm{G}}}}=\mathrm{0.02}\mathrm{m}$, which is a compromise between computational efficiency and ability to represent the stronger gradients near the surface, and layers with increasing thickness ($\mathrm{\Delta }{z}_{{\text{g}}_j}$) are added for the entire rooting zone.

The thermodynamic state is defined in terms of the soil volume: the bulk specific enthalpy ${\mathcal{H}}_{{\text{g}}_j}$ (J m⁻³) and volumetric soil water content ${\mathit{\vartheta }}_{{\text{g}}_j}$ (${\mathrm{m}}_{\mathrm{W}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$), which can be related to Eqs. (4)–(5) by defining $H_{{\text{g}}_j}={\mathcal{H}}_{{\text{g}}_j}\mathrm{\Delta }{z}_{{\text{g}}_j}$ and $W_{{\text{g}}_j}={\mathit{\rho }}_{\mathrm{\ell }}\cdot {\mathit{\vartheta }}_{{\text{g}}_j}\cdot \mathrm{\Delta }{z}_{{\text{g}}_j}$, where ρ_ℓ is the density of liquid water (Table S3). Soil net fluxes for any layer j are defined as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \underset{\text{Net heat flux}}{\underbrace{{\dot{Q}}_{{\text{g}}_j}}}=\underset{\begin{array}{c}\text{Net sensible heat flux}\\ \text{between consecutive layers}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\dot{Q}}_{g_{j-\mathrm{1}},{\text{g}}_j}-{\dot{Q}}_{{\text{g}}_j,g_{j+\mathrm{1}}}}}+\underset{\begin{array}{c}\text{Absorbed irradiance}\\ (\text{Sect}.\mathrm{4.3.2})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{g}}_j{g}_{N_{\mathrm{G}}}}{\dot{Q}}_{a,g_{N_{\mathrm{G}}}}}}\\ & -\underset{\begin{array}{c}\text{Ground–CAS sensible heat}\\ (\text{Sect}.\mathrm{4.5.2}\text{and}\mathrm{4.5.3})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{g}}_j{g}_{N_{\mathrm{G}}}}{\dot{Q}}_{g_{N_{\mathrm{G}}},c}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \underset{\begin{array}{c}\text{Net enthalpy flux}\\ \text{due to water flux}\end{array}}{\underbrace{{\dot{H}}_{{\text{g}}_j}}}=\underset{\begin{array}{c}\text{Water percolation}\\ \text{between consecutive layers}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\dot{H}}_{g_{j-\mathrm{1}},{\text{g}}_j}-{\dot{H}}_{{\text{g}}_j,g_{j+\mathrm{1}}}}}-\underset{\begin{array}{c}\text{Ground evaporation}\\ (\text{Sect}.\mathrm{4.5.2}\text{and}\mathrm{4.5.3})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{g}}_j{g}_{N_{\mathrm{G}}}}{\dot{H}}_{g_{N_{\mathrm{G}}},c}}}\\ & -\underset{\begin{array}{c}\text{Water uptake}\\ \text{by cohorts}\\ (\text{Sect}.\mathrm{4.6})\end{array}}{\underbrace{\sum _{k=\mathrm{1}}^{N_T}{\dot{H}}_{{\text{g}}_j,l_k}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \underset{\text{Net water flux}}{\underbrace{{\dot{W}}_{{\text{g}}_j}}}=\underset{\begin{array}{c}\text{Water percolation}\\ \text{between consecutive layers}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\dot{W}}_{g_{j-\mathrm{1}},{\text{g}}_j}-{\dot{W}}_{{\text{g}}_j,g_{j+\mathrm{1}}}}}-\underset{\begin{array}{c}\text{Ground evaporation}\\ (\text{Sect}.\mathrm{4.5.2}\text{and}\mathrm{4.5.3})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{g}}_j{g}_{N_{\mathrm{G}}}}{\dot{W}}_{g_{N_G},c}}}\\ & -\underset{\begin{array}{c}\text{water uptake}\\ \text{by cohorts}\\ (\text{Sect}.\mathrm{4.6})\end{array}}{\underbrace{\sum _{k=\mathrm{1}}^{N_{\mathrm{T}}}{\dot{W}}_{{\text{g}}_j,l_k}}},\end{array}\end{array}$$

where ${\mathit{\delta }}_{{\text{g}}_j{g}_{j^{\prime }}}$ is the Kronecker delta for comparing two soil layers g_j and $g_{j^{\prime }}$ (1 if ${\text{g}}_j=g_{j^{\prime }}$; 0 otherwise), CAS is the canopy air space, and subscript o denotes the loss through runoff. References in parentheses underneath the terms correspond to the sections in which each term is presented in detail. In the equations above, we assume ${\dot{Q}}_{g_{\mathrm{0}},g_{\mathrm{1}}}$ to be zero (bottom boundary condition in thermal equilibrium) and (${\dot{H}}_{g_{\mathrm{1}},g_{\mathrm{0}}}=-{\dot{H}}_{g_{\mathrm{0}},g_{\mathrm{1}}}$; ${\dot{W}}_{g_{\mathrm{1}},g_{\mathrm{0}}}=-{\dot{W}}_{g_{\mathrm{0}},g_{\mathrm{1}}}$) to be subsurface runoff fluxes (see Sect. 4.1). In addition, (${\dot{Q}}_{g_{N_{\mathrm{G}}},g_{N_{\mathrm{G}}+\mathrm{1}}}$; ${\dot{H}}_{g_{N_{\mathrm{G}}},g_{N_{\mathrm{G}}+\mathrm{1}}}$; ${\dot{W}}_{g_{N_{\mathrm{G}}},g_{N_{\mathrm{G}}+\mathrm{1}}})$ are equivalent to (${\dot{Q}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$; ${\dot{H}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$; ${\dot{W}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$), which are the fluxes between the topmost soil layer and the bottommost temporary surface water layer (see also Sect. 3.2.2).

3.2.2 Temporary surface water (TSW)

Temporary surface water (TSW) exists whenever water falls to the ground, or dew or frost develops on the ground. The layer will be maintained only when the amount of water that reaches the ground exceeds the water holding capacity of the top soil layer (a function of the soil porosity), or when precipitation falls as snow. The maximum number of temporary surface water layers $N_{\mathrm{S}}^{\text{max}}$ is defined by the user, but the actual number of layers N_S and the thickness of each layer depends on the total mass and the water phase, following Walko et al. (2000). When the layer is in liquid phase, only one layer (N_S=1) is maintained. If a snowpack develops, the temporary surface water can be divided into several layers (subscript j, with j=1 being the bottommost TSW layer, and j=N_S being the topmost TSW layer). Net TSW fluxes are defined as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \underset{\text{Net heat flux}}{\underbrace{{\dot{Q}}_{{\text{s}}_j}}}=\underset{\begin{array}{c}\text{Net sensible heat flux}\\ \text{between consecutive layers}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\dot{Q}}_{s_{j-\mathrm{1}},{\text{s}}_j}-{\dot{Q}}_{{\text{s}}_j,s_{j+\mathrm{1}}}}}+\underset{\begin{array}{c}\text{Absorbed irradiance}\\ \left(\text{Sect.}\mathrm{4.3.2}\right)\end{array}}{\underbrace{{\dot{Q}}_{a,{\text{s}}_j}}}\\ & -\underset{\begin{array}{c}\text{Ground–CAS sensible heat}\\ (\text{Sect}.\mathrm{4.5.2}\text{and}\mathrm{4.5.3})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{Q}}_{s_{N_{\mathrm{S}}},c}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \underset{\begin{array}{c}\text{Net enthalpy flux}\\ \text{due to water flux}\end{array}}{\underbrace{{\dot{H}}_{{\text{s}}_j}}}=\underset{\begin{array}{c}\text{Water percolation}\\ \text{between consecutive layers}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\dot{H}}_{s_{j-\mathrm{1}},{\text{s}}_j}-{\dot{H}}_{{\text{s}}_j,s_{j+\mathrm{1}}}}}+\underset{\begin{array}{c}\text{Throughfall}\\ \text{precipitation}\\ (\text{Sect}.\mathrm{4.2})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{H}}_{a,s_{N_{\mathrm{S}}}}}}\\ & +\underset{\begin{array}{c}\text{Canopy dripping}\\ \text{from cohorts}\\ \text{Sect}.\left(\mathrm{4.2}\right)\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}\left(\sum _{k=\mathrm{1}}^{N_{\mathrm{T}}}{\dot{H}}_{t_k,s_{N_{\mathrm{S}}}}\right)}}-\underset{\begin{array}{c}\text{Surface runoff}\\ \text{Sect}.\left(\mathrm{4.1}\right)\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{H}}_{s_{N_{\mathrm{S}}},o}}}-\underset{\begin{array}{c}\text{Surface water}\\ \text{evaporation}\\ (\text{Sect}.\mathrm{4.5.2}\text{and}\mathrm{4.5.3})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{H}}_{s_{N_{\mathrm{S}}},c}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \underset{\text{Net water flux}}{\underbrace{{\dot{W}}_{{\text{s}}_j}}}=\underset{\begin{array}{c}\text{Water percolation}\\ \text{between consecutive layers}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\dot{W}}_{s_{j-\mathrm{1}},{\text{s}}_j}-{\dot{W}}_{{\text{s}}_j,s_{j+\mathrm{1}}}}}+\underset{\begin{array}{c}\text{Throughfall}\\ \text{precipitation}\\ (\text{Sect}.\mathrm{4.2})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{W}}_{a,s_{N_{\mathrm{S}}}}}}\\ & +\underset{\begin{array}{c}\text{Canopy dripping}\\ \text{from cohorts}\\ (\text{Sect}.\mathrm{4.2})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}\left(\sum _{k=\mathrm{1}}^{N_{\mathrm{T}}}{\dot{W}}_{t_k,s_{N_{\mathrm{S}}}}\right)}}-\underset{\begin{array}{c}\text{Surface runoff}\\ (\text{Sect}.\mathrm{4.1})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{W}}_{s_{N_{\mathrm{S}}},o}}}\\ & -\underset{\begin{array}{c}\text{Surface water}\\ \text{evaporation}\\ (\text{Sect}.\mathrm{4.5.2}\text{and}\mathrm{4.5.3})\end{array}}{\underbrace{{\mathit{\delta }}_{{\text{s}}_j{s}_{N_{\mathrm{S}}}}{\dot{W}}_{s_{N_{\mathrm{S}}},c}}},\end{array}\end{array}$$

where ${\mathit{\delta }}_{{\text{s}}_j{s}_{j^{\prime }}}$ is the Kronecker delta for comparing two TSW layers s_j and $s_{j^{\prime }}$ (1 if ${\text{s}}_j=s_{j^{\prime }}$; 0 otherwise), CAS is the canopy air space, and subscript o denotes loss from the thermodynamic envelope through runoff. Terms are described in detail in the sections shown underneath each term. Similarly to the soil fluxes (Sect. 3.2.1), we assume that (${\dot{Q}}_{s_{\mathrm{0}},s_{\mathrm{1}}}$; ${\dot{H}}_{s_{\mathrm{0}},s_{\mathrm{1}}}$; ${\dot{W}}_{s_{\mathrm{0}},s_{\mathrm{1}}}$) is equivalent to (${\dot{Q}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$; ${\dot{H}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$; ${\dot{W}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$), the fluxes between the topmost soil layer and the bottommost TSW layer. When solving Eqs. (9)–(11) for layer $s_{N_{\mathrm{S}}}$, we assume the terms ${\dot{Q}}_{{\text{s}}_j,s_{j+\mathrm{1}}}$, ${\dot{H}}_{{\text{s}}_j,s_{j+\mathrm{1}}}$ and ${\dot{W}}_{{\text{s}}_j,s_{j+\mathrm{1}}}$ to be all zero, as layer N_S+1 does not exist.

In the case of liquid TSW, the layer thickness of the single layer is defined as $\mathrm{\Delta }{z}_{s_{\mathrm{1}}}={\mathit{\rho }}_{\mathrm{\ell }}^{-\mathrm{1}}{W}_{s_{\mathrm{1}}}$, where ρ_ℓ is the density of liquid water (Table S3). In the case of snowpack development, the snow density and the layer thickness of the TSW are solved as described in Sect. S7. The thickness of each layer of snow ($\mathrm{\Delta }{z}_{{\text{s}}_j}$) is defined using the same algorithm as LEAF-2 (Walko et al., 2000) and is described in Sect. S7.

3.2.3 Vegetation

In ED-2.2, vegetation is solved as an independent thermodynamic system only if the cohort is sufficiently large. The minimum size is an adjustable parameter and the typical minimum heat capacity solved by ED-2.2 is on the order of $\mathrm{10}\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$ and total area index of $\mathrm{0.005}{\mathrm{m}}_{\mathrm{leaf}+\mathrm{wood}}^{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}$. Cohorts smaller than this are excluded from all energy and water cycle calculations and assumed to be in thermal equilibrium with canopy air space. The net fluxes of heat, enthalpy, and water for each cohort k that can be resolved are

$$\begin{array}{}\text{(12)}& \underset{\text{Net heat flux}}{\underbrace{{\dot{Q}}_{t_k}}}=\underset{\begin{array}{c}\text{Cohort's net}\\ \text{absorbed irradiance}\\ (\text{Sect}.\mathrm{4.3.1})\end{array}}{\underbrace{{\dot{Q}}_{a,t_k}}}-\underset{\begin{array}{c}\text{Cohort–CAS}\\ \text{sensible heat}\\ \left(\text{Sect.}\mathrm{4.5.1}\right)\end{array}}{\underbrace{{\dot{Q}}_{t_k,c}}},\end{array}$$

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \underset{\begin{array}{c}\text{Net enthalpy flux}\\ \text{due to water flux}\end{array}}{\underbrace{{\dot{H}}_{t_k}}}=\underset{\begin{array}{c}\text{Rainfall}\\ \text{interception}\\ (\text{Sect}.\mathrm{4.2})\end{array}}{\underbrace{{\dot{H}}_{a,t_k}}}-\underset{\begin{array}{c}\text{Canopy dripping}\\ \left(\text{Sect.}\mathrm{4.2}\right)\end{array}}{\underbrace{{\dot{H}}_{t_k,s_{N_{\mathrm{S}}}}}}+\underset{\begin{array}{c}\text{Soil moisture}\\ \text{uptake (transpiration)}\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{\left(\sum _{j=\mathrm{1}}^{N_{\mathrm{G}}}{\dot{H}}_{{\text{g}}_j,l_k}\right)}}\\ & -\underset{\begin{array}{c}\text{Transpiration}\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{{\dot{H}}_{l_k,c}}}-\underset{\begin{array}{c}\text{Evaporation of}\\ \text{intercepted water}\\ \left(\text{Sect.}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{H}}_{t_k,c}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \underset{\text{Water flux}}{\underbrace{{\dot{W}}_{t_k}}}=\underset{\begin{array}{c}\text{Rainfall}\\ \text{interception}\\ \left(\text{Sect.}\mathrm{4.2}\right)\end{array}}{\underbrace{{\dot{W}}_{a,t_k}}}-\underset{\begin{array}{c}\text{Canopy dripping}\\ \left(\text{Sect.}\mathrm{4.2}\right)\end{array}}{\underbrace{{\dot{W}}_{t_k,s_{N_{\mathrm{S}}}}}}+\underset{\begin{array}{c}\text{Soil moisture}\\ \text{uptake (transpiration)}\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{\left(\sum _{j=\mathrm{1}}^{N_{\mathrm{G}}}{\dot{W}}_{{\text{g}}_j,l_k}\right)}}\\ & -\underset{\begin{array}{c}\text{Transpiration}\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{{\dot{W}}_{l_k,c}}}-\underset{\begin{array}{c}\text{Evaporation of}\\ \text{intercepted water}\\ \left(\text{Sect.}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{W}}_{t_k,c}}}.\end{array}\end{array}$$

Each term is described in detail in the sections shown underneath each term on the right-hand side of Eqs. (12)–(14).

3.2.4 Canopy air space (CAS)

The canopy air space is a gas; therefore, extensive properties akin to the other thermodynamic systems are not intuitive because total mass and total volume cannot be directly compared to observations. Therefore, all prognostic and diagnostic variables are solved in the intensive form. Total enthalpy H_c and total water mass W_c of the canopy air space can be written in terms of air density ρ_c and the equivalent depth of the canopy air space ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}$ as

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{H}_{\mathrm{c}}={\mathit{\rho }}_{\mathrm{c}}{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}{h}_{\mathrm{c}},\text{(16)}& {\displaystyle }& {\displaystyle }{W}_{\mathrm{c}}={\mathit{\rho }}_{\mathrm{c}}{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}{w}_{\mathrm{c}},\text{(17)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}=max\left(\mathrm{5.0},{\displaystyle \frac{{\sum }_{k=\mathrm{1}}^{N_{T\left(\text{canopy}\right)}}{n}_{t_k}{\mathrm{BA}}_{t_k}{z}_{t_k}}{{\sum }_{k=\mathrm{1}}^{N_{T\left(\text{canopy}\right)}}{n}_{t_k}{\mathrm{BA}}_{t_k}}}\right),\end{array}$$

where ${\mathrm{BA}}_{t_k}$ (cm²) and $z_{t_k}$ (m) are the basal area and the height of cohort k, respectively; and N_T(canopy) is the number of cohorts that are in the canopy, and we assume that cohorts are ordered from tallest to shortest. In the event that the canopy is open, N_T(canopy) is the total number of cohorts, and a minimum value of 5 m is imposed when vegetation is absent or too short to prevent numerical instabilities. Because the equivalent canopy depth depends only on the cohort size, ${\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}$ is updated at the cohort dynamics step (Δt_CD, Table 2). If we substitute Eqs. (15) and (16) into Eqs. (4) and (5), respectively, and assume that changes in density over short time steps are much smaller than changes in enthalpy or humidity, and then we obtain the following equations for the canopy air space budget:

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{c}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{c}}{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}}}\left({\dot{Q}}_{\mathrm{c}}+{\dot{H}}_{\mathrm{c}}+{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}{\displaystyle \frac{\mathrm{d}{p}_{\mathrm{c}}}{\mathrm{d}t}}\right),\text{(19)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{w}_{\mathrm{c}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{c}}{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}}}{\dot{W}}_{\mathrm{c}},\end{array}$$

where

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& \underset{\text{Net heat flux}}{\underbrace{{\dot{Q}}_{\mathrm{c}}}}=\underset{\begin{array}{c}\text{Cohort–CAS}\\ \text{sensible heat}\\ \left(\text{Sect.}\mathrm{4.5.1}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{\left(\sum _{k=\mathrm{1}}^{N_T}{\dot{Q}}_{t_k,c}\right)}}+\underset{\begin{array}{c}\text{Surface water–CAS}\\ \text{sensible heat}\\ \left(\text{Sect.}\mathrm{4.5.2}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{Q}}_{s_{N_{\mathrm{S}}},c}}}\\ & +\underset{\begin{array}{c}\text{Ground–CAS}\\ \text{sensible heat}\\ \left(\text{Sect.}\mathrm{4.5.2}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{Q}}_{g_{N_{\mathrm{G}}},c}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& \underset{\text{Net enthalpy flux}}{\underbrace{{\dot{H}}_{\mathrm{c}}}}=\underset{\begin{array}{c}\text{Enthalpy flux from}\\ \text{turbulent mixing}\left(\text{Sect.}\mathrm{4.4}\right)\end{array}}{\underbrace{{\dot{H}}_{a,c}}}+\underset{\begin{array}{c}\text{Evaporation of}\\ \text{intercepted water}\\ \left(\text{Sect.}\mathrm{4.5.1}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{\left(\sum _{k=\mathrm{1}}^{N_T}{\dot{H}}_{t_k,c}\right)}}\\ & +\underset{\begin{array}{c}\text{Transpiration}\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{\left(\sum _{k=\mathrm{1}}^{N_T}{\dot{H}}_{l_k,c}\right)}}+\underset{\begin{array}{c}\text{Surface water}\\ \text{evaporation}\\ \left(\text{Sect.}\mathrm{4.5.2}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{H}}_{s_{N_{\mathrm{S}}},c}}}+\underset{\begin{array}{c}\text{Ground evaporation}\\ \left(\text{Sect.}\mathrm{4.5.2}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{H}}_{g_{N_{\mathrm{G}}},c}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& \underset{\text{Net water flux}}{\underbrace{{\dot{W}}_{\mathrm{c}}}}=\underset{\begin{array}{c}\text{Water flux from}\\ \text{turbulent mixing}\left(\text{Sect.}\mathrm{4.4}\right)\end{array}}{\underbrace{{\dot{W}}_{a,c}}}+\underset{\begin{array}{c}\text{Evaporation of}\\ \text{intercepted water}\\ \left(\text{Sect.}\mathrm{4.5.1}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{\left(\sum _{k=\mathrm{1}}^{N_T}{\dot{W}}_{t_k,c}\right)}}\\ & +\underset{\begin{array}{c}\text{Transpiration}\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{\left(\sum _{k=\mathrm{1}}^{N_T}{\dot{W}}_{l_k,c}\right)}}+\underset{\begin{array}{c}\text{Surface water}\\ \text{evaporation}\\ \left(\text{Sect.}\mathrm{4.5.2}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{W}}_{s_{N_{\mathrm{S}}},c}}}+\underset{\begin{array}{c}\text{Ground evaporation}\\ \left(\text{Sect.}\mathrm{4.5.2}\text{and}\mathrm{4.5.3}\right)\end{array}}{\underbrace{{\dot{W}}_{g_{N_{\mathrm{G}}},c}}}.\end{array}\end{array}$$

Unlike in the other thermodynamic systems (soil, temporary surface water, and vegetation), the net enthalpy flux of the canopy air space is not exclusively due to associated water flux: the eddy flux between the free air and the canopy air space (${\dot{H}}_{\mathrm{a},\mathrm{c}}$) includes both water transport and flux associated with mixing of air with different temperatures, and thus enthalpy, between canopy air space and free air.

In addition, we must also track the CAS pressure p_c. In ED-2.2, CAS pressure is not solved through a differential equation; instead, p_c is updated whenever the meteorological forcing is updated, using the ideal gas law and hydrostatic equilibrium following the method described in Sect. S8. The rate of change of canopy air pressure is then applied in Eq. (18). Likewise, CAS density (ρ_c) is updated at the end of each thermodynamic step to ensure that the CAS conforms to the ideal gas law.

3.3 Carbon dioxide cycle

In ED-2.2, the carbon dioxide cycle is a subset of the full carbon cycle, which is shown in Fig. 3. The canopy air space is the only thermodynamic system with CO₂ storage that is solved by ED-2.2; nonetheless, we assume that the contribution of CO₂ to density and heat capacity of the canopy air space is negligible; hence, only the molar CO₂ mixing ratio c_c (mol_C mol⁻¹) is traced.

Figure 3 Schematic of the patch-level carbon cycle solved in ED-2.2 for a patch containing N_T cohorts. Like Fig. 2, letters near the arrows are the subscripts associated with fluxes. Fluxes shown in solid yellow lines are part of the CO₂ cycle discussed in this paper, and dashed red lines are part of the carbon cycle but do not directly affect the CO₂ flux; these fluxes are summarized in Sects. S3 and S4.

The change in CO₂ storage in the canopy air space is determined by the following differential equation:

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\mathrm{d}{c}_{\mathrm{c}}}{\mathrm{d}t}}={\displaystyle \frac{{\mathcal{M}}_{\mathrm{d}}}{{\mathcal{M}}_{\mathrm{C}}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{c}}{\stackrel{\mathrm{‾}}{z}}_{\mathrm{c}}}}{\dot{C}}_{\mathrm{c}},\end{array}$$

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& \underset{\text{Net carbon flux}}{\underbrace{{\dot{C}}_{\mathrm{c}}}}=\underset{\begin{array}{c}\text{Carbon flux from}\\ \text{turbulent mixing}\\ \left(\text{Sect.}\mathrm{4.4}\right)\end{array}}{\underbrace{{\dot{C}}_{\mathrm{a},\mathrm{c}}}}+\underset{\begin{array}{c}\text{Net leaf–CAS flux}\\ \left(\text{respiration–GPP}\right)\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{\sum _{k=\mathrm{1}}^{N_T}{\dot{C}}_{l_k,c}}}+\underset{\begin{array}{c}\text{Fine-root}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.7}\right)\end{array}}{\underbrace{\sum _{k=\mathrm{1}}^{N_{\mathrm{T}}}{\dot{C}}_{r_k,c}}}\\ & +\underset{\begin{array}{c}\text{Storage turnover}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.7}\right)\end{array}}{\underbrace{\sum _{k=\mathrm{1}}^{N_{\mathrm{T}}}{\dot{C}}_{n_k,c}}}+\underset{\begin{array}{c}\text{Growth and maintenance}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.7}\right)\end{array}}{\underbrace{\sum _{k=\mathrm{1}}^{N_{\mathrm{T}}}{\dot{C}}_{{\mathrm{\Delta }}_k,c}}}+\underset{\begin{array}{c}\text{Heterotrophic}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.8}\right)\end{array}}{\underbrace{\sum _{j=\mathrm{1}}^{\mathrm{3}}{\dot{C}}_{e_j,c}}},\end{array}\end{array}$$

where ℳ_d and ℳ_C are the molar masses of dry air and carbon, respectively, used to convert mass to molar fraction ($\mathrm{1}{\mathrm{mol}}_{\mathrm{C}}=\mathrm{1}{\mathrm{mol}}_{{\mathrm{CO}}_{\mathrm{2}}}$). The terms on the right-hand side of Eq. (24) are described in detail in the sections displayed underneath each term. The net leaf–CAS flux (${\dot{C}}_{l_k,c}$) for any cohort k is positive when leaf respiration exceeds photosynthetic assimilation. The heterotrophic respiration is based on a simplified implementation of the CENTURY model (Bolker et al., 1998) that combines the decomposition rates from three soil carbon pools, defined by their characteristic lifetime: fast (metabolic litter and microbial; e₁), intermediate (structural debris; e₂), and slow (humified and passive soil carbon; e₃). Note that the soil carbon pools are not directly related to the soil layers used to describe the thermodynamic state (Sect. 3.2.1).

In addition to canopy air space, we also define a virtual cohort pool of carbon corresponding to the accumulated carbon balance ($C_{{\mathrm{\Delta }}_k}$). The accumulated carbon balance links short-term carbon cycle components such as photosynthesis and respiration with long-term dynamics that depend on carbon balance such as such as carbon allocation to growth and reproduction, and mortality (long-term dynamics described in Sect. S3). The accumulated carbon balance is defined by the following equation:

$$\begin{array}{}\text{(25)}& \begin{array}{rl}& \underset{\begin{array}{c}\text{Change in}\\ \text{carbon balance}\end{array}}{\underbrace{{\displaystyle \frac{\mathrm{d}{C}_{{\mathrm{\Delta }}_k}}{\mathrm{d}t}}}}=-\underset{\begin{array}{c}\text{Net leaf–CAS flux}\\ \left(\text{respiration–GPP}\right)\\ \left(\text{Sect.}\mathrm{4.6}\right)\end{array}}{\underbrace{{\dot{C}}_{l_k,c}}}-\underset{\begin{array}{c}\text{Fine-root}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.7}\right)\end{array}}{\underbrace{{\dot{C}}_{r_k,c}}}-\underset{\begin{array}{c}\text{Storage turnover}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.7}\right)\end{array}}{\underbrace{{\dot{C}}_{n_k,c}}}\\ & -\underset{\begin{array}{c}\text{Growth and maintenance}\\ \text{respiration}\\ \left(\text{Sect.}\mathrm{4.7}\right)\end{array}}{\underbrace{{\dot{C}}_{{\mathrm{\Delta }}_k,c}}}-\underset{\begin{array}{c}\text{Turnover of}\\ \text{non-lignified litter}\\ (\text{Sects.}\mathrm{4},\text{S4})\end{array}}{\underbrace{{\dot{C}}_{t_k,e_{\mathrm{1}}}}}-\underset{\begin{array}{c}\text{Turnover of}\\ \text{lignified litter}\\ (\text{Sects.}\mathrm{4},\text{S4}),\end{array}}{\underbrace{{\dot{C}}_{t_k,e_{\mathrm{2}}}}}\end{array}\end{array}$$

where ${\dot{C}}_{t_k,e_{\mathrm{1}}}$ and ${\dot{C}}_{t_k,e_{\mathrm{2}}}$ are the individual carbon losses caused by leaf shedding and turnover of living tissues that become part of the litter (${\dot{C}}_{t_k,e_{\mathrm{1}}}$) and structural debris (${\dot{C}}_{t_k,e_{\mathrm{2}}}$). The transfer of carbon from plants to the soil carbon pools and between the soil carbon pools does not directly impact the carbon dioxide budget but contributes to the long-term ecosystem carbon stock distribution and carbon balance. These components have been discussed in previous ED and ED-2 publications (Moorcroft et al., 2001; Albani et al., 2006; Medvigy, 2006; Medvigy et al., 2009) and are summarized in Sect. S4.

4 Submodels and parameterizations of terms of the general equations

4.1 Hydrology submodel and ground energy exchange

The ground model encompasses heat, enthalpy, and water fluxes between adjacent layers of soil and temporary surface water, as well as losses of water and enthalpy due to surface runoff and drainage. Fluxes between adjacent layers are positive when they are upwards, and runoff and drainage fluxes are positive or zero.

Sensible heat flux between two adjacent soil or temporary surface water layers j−1 and j is determined from thermal conductivity Υ_Q and temperature gradient (Bonan, 2008), with an additional term for temporary surface water to scale the flux when the temporary surface water covers only a fraction f_TSW of the ground:

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{\dot{Q}}_{g_{j-\mathrm{1}},{\text{g}}_j}=-{〈{\mathrm{{\rm Y}}}_{\mathrm{Q}}〉}_{g_{j-\mathrm{1}},{\text{g}}_j}{\left({\displaystyle \frac{\partial T_g}{\partial z}}\right)}_{g_{j-\mathrm{1}},{\text{g}}_j},\text{(27)}& {\displaystyle }& {\displaystyle }{\dot{Q}}_{s_{j-\mathrm{1}},{\text{s}}_j}=-f_{\text{TSW}}{〈{\mathrm{{\rm Y}}}_{\mathrm{Q}}〉}_{s_{j-\mathrm{1}},{\text{s}}_j}{\left({\displaystyle \frac{\partial T_s}{\partial z}}\right)}_{s_{j-\mathrm{1}},{\text{s}}_j},\end{array}$$

where the operator 〈〉 is the log-linear interpolation from the midpoint height of layers j−1 and j to the height at the interface. The bottom boundary condition of Eq. (26) is ${\left(\frac{\partial T}{\partial z}\right)}_{g_{\mathrm{0}},g_{\mathrm{1}}}\equiv \mathrm{0}$. The interface between the top soil layer and the first temporary surface water (${\dot{Q}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}$) is found by applying Eq. (27) with $\left(T_{s_{\mathrm{0}}};{\mathrm{{\rm Y}}}_{{\mathrm{Q}}_{s_{\mathrm{0}}}};\mathrm{\Delta }{z}_{s_{\mathrm{0}}}\right)=\left(T_{g_{N_{\mathrm{G}}}};{\mathrm{{\rm Y}}}_{{\mathrm{Q}}_{g_{N_{\mathrm{G}}}}};\mathrm{\Delta }{z}_{g_{N_{\mathrm{G}}}}\right)$. Soil thermal conductivity depends on soil moisture and texture properties, and the parameterization is described in Sect. S9. Both the fraction of ground covered by the temporary surface water and the thermal conductivity of the temporary surface water are described in Sect. S10.

Soil moisture exchange between layers occurs only if water is in liquid phase. The water flux between soil layers g_j−1 and ${\text{g}}_j,j\in \left\{\mathrm{2},\mathrm{3},\mathrm{\dots },N_{\mathrm{G}}\right\}$ is determined from Darcy's law (Bonan, 2008):

$$\begin{array}{}\text{(28)}& {\dot{W}}_{g_{j-\mathrm{1}},{\text{g}}_j}=-{\mathit{\rho }}_{\mathrm{\ell }}{〈{\mathrm{{\rm Y}}}_{\mathrm{\Psi }}〉}_{g_{j-\mathrm{1}},{\text{g}}_j}{\left[{\displaystyle \frac{\partial \mathrm{\Psi }}{\partial z}}+{\displaystyle \frac{\mathrm{d}{z}_g}{\mathrm{d}z}}\right]}_{g_{j-\mathrm{1}},{\text{g}}_j},\end{array}$$

where Ψ is the soil matric potential and Υ_Ψ is the hydraulic conductivity, both defined after Brooks and Corey (1964), with an additional correction term applied to hydraulic conductivity to reduce conductivity in the event that the soil is partially or completely frozen (Sect. S9). The bottom boundary condition for soil matric potential gradient is ${\left(\frac{\partial \mathrm{\Psi }}{\partial z}\right)}_{g_{\mathrm{0}},g_{\mathrm{1}}}\equiv \mathrm{0}$.

The term $\frac{\mathrm{d}{z}_g}{\mathrm{d}z}$ in Eq. (28) is the flux due to gravity, and it is 1 for all layers except the bottom boundary condition, which depends on the subsurface drainage. Subsurface drainage at the bottom boundary depends on the type of drainage and is determined using a slight modification of Eq. (28). Let ð be an angle-like parameter that controls the drainage beneath the deepest soil layer. Because we assume zero gradient in soil matric potential between the deepest soil layer and the boundary condition, the subsurface drainage flux (${\dot{W}}_{g_{\mathrm{1}},g_{\mathrm{0}}}$) becomes

$$\begin{array}{}\text{(29)}& {\dot{W}}_{g_{\mathrm{1}},g_{\mathrm{0}}}=-{\dot{W}}_{g_{\mathrm{0}},g_{\mathrm{1}}}={\mathit{\rho }}_{\mathrm{\ell }}{\mathrm{{\rm Y}}}_{{\mathrm{\Psi }}_{g_{\mathrm{1}}}}\sin \mathit{ð}.\end{array}$$

Special cases of Eq. (29) are the zero-flow conditions (ð=0) and free drainage ($\mathit{ð}=\frac{\mathit{\pi }}{\mathrm{2}}$).

For the temporary surface water, water flux between layers through percolation is calculated similarly to LEAF-2 (Walko et al., 2000). Liquid water in excess of 10  % is in principle free to percolate to the layer below, although the maximum percolation of the first surface water layer is limited by the amount of pore space available at the top ground layer:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& {\dot{W}}_{s_{\mathrm{1}},g_{N_{\mathrm{G}}}}=-{\dot{W}}_{g_{N_{\mathrm{G}}},s_{\mathrm{1}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{t}_{\text{Thermo}}}}\\ & max\left\{\mathrm{0},min\left[W_{s_{\mathrm{1}}}\left({\displaystyle \frac{{\mathrm{\ell }}_{s_j}-\mathrm{0.1}}{\mathrm{0.9}}}\right),{\mathit{\rho }}_{\mathrm{\ell }}\left({\mathit{\vartheta }}_{\mathrm{Po}}-{\mathit{\vartheta }}_{g_{N_{\mathrm{G}}}}\right)\mathrm{\Delta }{z}_{g_{N_G}}\right]\right\},\end{array}\end{array}$$

$$\begin{array}{}\text{(31)}& \begin{array}{rl}& {\dot{W}}_{{\text{s}}_j,s_{j-\mathrm{1}}}=-{\dot{W}}_{s_{j-\mathrm{1}},{\text{s}}_j}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{t}_{\text{Thermo}}}}\\ & max\left(\mathrm{0},W_{{\text{s}}_j}{\displaystyle \frac{{\mathrm{\ell }}_{{\text{s}}_j}-\mathrm{0.1}}{\mathrm{0.9}}}\right),\text{ for }j>\mathrm{1}.\end{array}\end{array}$$

Surface runoff of liquid water is simulated using a simple extinction function, applied only at the topmost temporary surface water layer:

$$\begin{array}{}\text{(32)}& {\dot{W}}_{s_{N_{\mathrm{S}}},o}={\mathrm{\ell }}_{s_{N_{\mathrm{S}}}}{W}_{s_{N_{\mathrm{S}}}}\exp \left(-{\displaystyle \frac{\mathrm{\Delta }{t}_{\text{Thermo}}}{t_{\text{Runoff}}}}\right),\end{array}$$

where t_Runoff is a user-defined e-folding decay time, usually on the order of a few minutes to a few hours (Table S4).

In addition to the water fluxes due to subsurface drainage, surface runoff, and the transport of water between layers, we must account for the associated enthalpy fluxes. Enthalpy fluxes due to subsurface drainage and surface runoff are defined based on the water flux and the temperature of the layers where water is lost by applying the definition of enthalpy (Sect. S5):

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{g_{\mathrm{1}},g_{\mathrm{0}}}={\dot{W}}_{g_{\mathrm{1}},g_{\mathrm{0}}}{q}_{\mathrm{\ell }}\left(T_{g_{\mathrm{1}}}-T_{\mathrm{\ell }\mathrm{0}}\right),\text{(34)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{s_{N_{\mathrm{S}}},o}={\dot{W}}_{s_{N_{\mathrm{S}}},o}{q}_{\mathrm{\ell }}\left(T_{s_{N_{\mathrm{S}}}}-T_{\mathrm{\ell }\mathrm{0}}\right),\end{array}$$

where q_ℓ is the specific heat of liquid water (Table S3), and T_ℓ0 is defined in Eq. (S53). The enthalpy flux between two adjacent layers is solved similarly, but it must account for the sign of the flux in order to determine the water temperature of the donor layer:

$$\begin{array}{}\text{(35)}& {\dot{H}}_{x_{j-\mathrm{1}},x_j}=\left\{\begin{array}{ll}{\dot{W}}_{x_{j-\mathrm{1}},x_j}{q}_{\mathrm{\ell }}\left(T_{x_j}-T_{\mathrm{\ell }\mathrm{0}}\right),& \text{ if}{\dot{W}}_{x_{j-\mathrm{1}},x_j}<\mathrm{0}\\ {\dot{W}}_{x_{j-\mathrm{1}},x_j}{q}_{\mathrm{\ell }}\left(T_{x_{j-\mathrm{1}}}-T_{\mathrm{\ell }\mathrm{0}}\right),& \text{ if}{\dot{W}}_{x_{j-\mathrm{1}},x_j}\ge \mathrm{0}\end{array}\right.,\end{array}$$

where the subscript x_j represents either soil (g_j) or temporary surface water (s_j).

4.2 Precipitation and vegetation dripping

In ED-2.2, precipitating water from rain and snow increases the water storage of the thermodynamic systems, as rainfall can be intercepted by the canopy or reach the ground. This influx of water also affects the enthalpy storage due to the enthalpy associated with precipitation, although no heat exchange is directly associated with precipitation.

To determine the partitioning of total incoming precipitation (${\dot{W}}_{\mathrm{\infty },\mathrm{a}}$) into interception by each cohort (${\dot{W}}_{a,t_k}$) and direct interception by the ground (throughfall, ${\dot{W}}_{a,s_{N_{\mathrm{S}}}}$), we use the fraction of open canopy (𝒪) and the total plant area index of each cohort (${\mathrm{\Phi }}_{t_k}$):

$$\begin{array}{}\text{(36)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{a,t_k}=\left(\mathrm{1}-\mathcal{O}\right){\dot{W}}_{\mathrm{\infty },\mathrm{a}}{\displaystyle \frac{{\mathrm{\Phi }}_{t_k}}{{\sum }_{k^{\prime }=\mathrm{1}}^{N_T}{\mathrm{\Phi }}_{t_{k^{\prime }}}}},\text{(37)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{a,s_{N_{\mathrm{S}}}}=\mathcal{O}{\dot{W}}_{\mathrm{\infty },\mathrm{a}},\text{(38)}& {\displaystyle }& {\displaystyle }\mathcal{O}=\prod _{k=\mathrm{1}}^{N_T}\left(\mathrm{1}-X_{t_k}\right),\end{array}$$

where ${\mathrm{\Phi }}_{t_k}={\mathrm{\Lambda }}_{t_k}+{\mathrm{\Omega }}_{t_k}$ is the total plant area index, ${\mathrm{\Lambda }}_{t_k}$ and ${\mathrm{\Omega }}_{t_k}$ being the leaf and wood area indices, both defined from PFT-dependent allometric relations (Sect. S18); $X_{t_k}$ is the crown area index of each cohort, also defined in Sect. S18. Throughfall precipitation is always placed on the topmost temporary surface water layer. In the event that no temporary surface water layer exists, a new layer is created, although it may be extinct in the event that all water is able to percolate down to the top soil layer.

Precipitation is a mass flux, but it also has an associated enthalpy flux (${\dot{H}}_{\mathrm{\infty },\mathrm{a}}$) that must be partitioned and incorporated into the cohorts and temporary surface water. Similar to the water exchange between soil layers, the enthalpy flux associated with rainfall uses the definition of enthalpy (Sect. S5). Because precipitation temperature is seldom available in meteorological drivers (towers or gridded meteorological forcing data sets), we assume that precipitation temperature is closely associated with the free-air temperature (T_a), and we use T_a to determine whether the precipitation falls as rain, snow, or a mix of both. Importantly, the use of free-air temperature partly accounts for the thermal difference between precipitation temperature and the temperature of intercepted surfaces. Rain is only allowed when T_a is above the water triple point (T₃=273.16 K); in this case, the rain temperature is always assumed to be at T_a. Pure snow occurs when the free-air temperature is below T₃, and likewise snow temperature is assumed to be T_a. When free-air temperature is only slightly above T₃, a mix of rain and snow occurs, with the rain temperature assumed to be T_a and snow temperature assumed to be T₃:

$$\begin{array}{}\text{(39)}& {\dot{H}}_{\mathrm{\infty },\mathrm{a}}={\dot{W}}_{\mathrm{\infty },\mathrm{a}}\left[\left(\mathrm{1}-{\mathrm{\ell }}_a\right)q_{i}min\left(T_{\mathrm{3}},T_{\mathrm{a}}\right)+{\mathrm{\ell }}_a{q}_{\mathrm{\ell }}\left(T_{\mathrm{a}}-T_{\mathrm{\ell }\mathrm{0}}\right)\right],\end{array}$$

where (q_i;q_ℓ) are the specific heats of ice and liquid, respectively, and T_ℓ0 is temperature at which supercooled water would have enthalpy equal to zero (Eq. S53). The fraction of precipitation that falls as rain ℓ_a is based on the Jin et al. (1999) parameterization, slightly modified to make the function continuous:

$$\begin{array}{}\text{(40)}& {\mathrm{\ell }}_{\mathrm{a}}=\left\{\begin{array}{ll}\mathrm{1.0},& \text{ if}{T}_{\mathrm{a}}>\mathrm{275.66}\mathrm{K}\\ \mathrm{0.4}+\mathrm{1.2}\left(T_{\mathrm{a}}-T_{\mathrm{3}}-\mathrm{2.0}\right),& \text{ if}\mathrm{275.16}\mathrm{K}<T_{\mathrm{a}}\le \mathrm{275.66}\mathrm{K}\\ \mathrm{0.2}\left(T_{\mathrm{a}}-T_{\mathrm{3}}\right),& \text{ if}{T}_{\mathrm{3}}<T_{\mathrm{a}}\le \mathrm{275.16}\mathrm{K}\\ \mathrm{0.0},& \text{ if}{T}_{\mathrm{a}}\le T_{\mathrm{3}}\end{array}\right..\end{array}$$

The enthalpy flux associated with precipitation is then partitioned into canopy interception (${\dot{H}}_{a,t_k}$) and throughfall (${\dot{H}}_{a,s_{N_{\mathrm{S}}}}$) using the same scaling factor as in Eqs. (36) and (37):

$$\begin{array}{}\text{(41)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{a,t_k}=\left(\mathrm{1}-\mathcal{O}\right){\dot{H}}_{\mathrm{\infty },\mathrm{a}}{\displaystyle \frac{{\mathrm{\Phi }}_{t_k}}{{\sum }_{k^{\prime }=\mathrm{1}}^{N_T}{\mathrm{\Phi }}_{t_{k^{\prime }}}}},\text{(42)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{a,s_{N_{\mathrm{S}}}}=\mathcal{O}{\dot{H}}_{\mathrm{\infty },\mathrm{a}}.\end{array}$$

Leaves and branches can accumulate only a finite amount of water on their surfaces, proportional to their total area. When incoming precipitation rates are too high (or more rarely when dew or frost formation is excessive), any water amount that exceeds the holding capacity is lost to the ground as canopy dripping. Similarly to incoming precipitation, the excess water lost through dripping also has an associated enthalpy that must be taken into account, although dripping has no associated heat flux. The canopy dripping fluxes of water (${\dot{W}}_{t_k,s_{N_{\mathrm{S}}}}$) and the associated enthalpy (${\dot{H}}_{t_k,s_{N_{\mathrm{S}}}}$) are defined such that the leaves and branches lose the excess water within one time step:

$$\begin{array}{}\text{(43)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{t_k,s_{N_{\mathrm{S}}}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{t}_{\text{Thermo}}}}max\left(\mathrm{0},W_{t_k}-{\widehat{w}}_{max}{\mathrm{\Phi }}_{t_k}\right),\text{(44)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{t_k,s_{N_{\mathrm{S}}}}={\dot{W}}_{t_k,s_{N_{\mathrm{S}}}}\left[\left(\mathrm{1}-{\mathrm{\ell }}_{t_k}\right)q_i{T}_{t_k}+{\mathrm{\ell }}_{t_k}\left(T_{t_k}-T_{\mathrm{\ell }\mathrm{0}}\right)\right],\end{array}$$

where ${\mathrm{\ell }}_{t_k}$ is the liquid fraction of surface water on top of cohort k and ${\widehat{w}}_{max}$ is the cohort holding capacity, which is an adjustable parameter (Table S4) but is typically of the order of $\mathrm{0.05}-\mathrm{0.40}{\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}_{\mathrm{Leaf}+\mathrm{Wood}}^{-\mathrm{2}}$ (Wohlfahrt et al., 2006).

4.3 Radiation model

The radiation budget is solved using a multi-layer version of the two-stream model (Sellers, 1985; Liou, 2002; Medvigy, 2006) applied to three broad spectral bands: photosynthetically active radiation (PAR, wavelengths between 0.4 and 0.7 µm), near-infrared radiation (NIR, wavelengths between 0.7 and 3.0 µm) and thermal infrared radiation (TIR, wavelengths between 3.0 and 15 µm).

4.3.1 Canopy radiation profile

For each spectral band m, the canopy radiation scheme assumes that each cohort corresponds to one layer of vegetation within the canopy, and within each layer the optical and thermal properties are assumed constant. For all bands, the top boundary condition for each band is provided by the meteorological forcing (Table 3). In the cases of PAR (m=1) and NIR (m=2), the downward irradiance is comprised of a beam (direct) and isotropic (diffuse) components, whereas TIR irradiance (m=3) is assumed to be all diffuse. Direct irradiance that is intercepted by the cohorts can be either backscattered or forward-scattered as diffuse radiation, and direct radiation reflected by the ground is assumed to be entirely diffuse.

Following Sellers (1985), the extinction of downward direct irradiance and the two-stream model for hemispheric diffuse irradiance for each of the spectral bands ($m=\mathrm{1},\mathrm{2},\mathrm{3}$) is given by

$$\begin{array}{}\text{(45)}& \underset{\begin{array}{c}\text{Downward direct}\\ \text{profile}\end{array}}{\underbrace{{\mathit{\mu }}_k^{\odot }{\displaystyle \frac{\mathrm{d}{\dot{Q}}_{mk}^{\odot }}{\mathrm{d}\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}}}}=\underset{\text{Interception}}{\underbrace{-{\dot{Q}}_{mk}^{\odot }}},\end{array}$$

$$\begin{array}{}\text{(46)}& \begin{array}{rl}& \underset{\begin{array}{c}\text{Downward diffuse}\\ \text{profile}\end{array}}{\underbrace{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_k{\displaystyle \frac{\mathrm{d}{\dot{Q}}_{mk}^{\Downarrow }}{\mathrm{d}\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}}}}=\underset{\text{Interception}}{\underbrace{-{\dot{Q}}_{mk}^{\Downarrow }}}+\underset{\begin{array}{c}\text{Forward scattering}\\ \text{(downward diffuse)}\end{array}}{\underbrace{\left(\mathrm{1}-{\mathit{\beta }}_{mk}\right){\mathit{\varsigma }}_{mk}{\dot{Q}}_{mk}^{\Downarrow }}}\\ & +\underset{\begin{array}{c}\text{Backscattering}\\ \text{(upward diffuse)}\end{array}}{\underbrace{{\mathit{\beta }}_{ik}{\mathit{\varsigma }}_{mk}{\dot{Q}}_{mk}^{\Uparrow }}}+\underset{\begin{array}{c}\text{Forward scattering}\\ \text{(downward direct)}\end{array}}{\underbrace{{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_k}{{\mathit{\mu }}_k^{\odot }}}{\mathit{\varsigma }}_{mk}\left(\mathrm{1}-{\mathit{\beta }}_{mk}^{\odot }\right){\dot{Q}}_{mk}^{\odot }}}+\underset{\text{Emission}}{\underbrace{\left(\mathrm{1}-{\mathit{\varsigma }}_{mk}\right){\dot{Q}}_{mk}^{\mathit{⧫}}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(47)}& \begin{array}{rl}& \underset{\begin{array}{c}\text{Upward diffuse}\\ \text{profile}\end{array}}{\underbrace{-{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_k{\displaystyle \frac{\mathrm{d}{\dot{Q}}_{mk}^{\Uparrow }}{\mathrm{d}\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}}}}=-\underset{\text{Interception}}{\underbrace{{\dot{Q}}_{mk}^{\Uparrow }}}+\underset{\begin{array}{c}\text{Forward scattering}\\ \text{(upward diffuse)}\end{array}}{\underbrace{\left(\mathrm{1}-{\mathit{\beta }}_{mk}\right){\mathit{\varsigma }}_{mk}{\dot{Q}}_{mk}^{\Uparrow }}}\\ & +\underset{\begin{array}{c}\text{Backscattering}\\ \text{(downward diffuse)}\end{array}}{\underbrace{{\mathit{\beta }}_{mk}{\mathit{\varsigma }}_{mk}{\dot{Q}}_{mk}^{\Downarrow }}}+\underset{\begin{array}{c}\text{Backscattering}\\ \text{(downward direct)}\end{array}}{\underbrace{{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_k}{{\mathit{\mu }}_k^{\odot }}}{\mathit{\varsigma }}_{mk}{\mathit{\beta }}_{mk}^{\odot }{\dot{Q}}_{mk}^{\odot }}}+\underset{\text{Emission}}{\underbrace{\left(\mathrm{1}-{\mathit{\varsigma }}_{mk}\right){\dot{Q}}_{mk}^{\mathit{⧫}}}},\end{array}\end{array}$$

where index $k\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{T}}\right\}$ corresponds to each cohort k or its lower interface (Fig. 4); interface N_T+1 is immediately above the tallest cohort; ${\dot{Q}}_{mk}^{\odot }$ is the downward direct irradiance incident at interface k; (${\dot{Q}}_{mk}^{\Downarrow }$ and ${\dot{Q}}_{mk}^{\Uparrow }$) are the downward and upward (hemispheric) diffuse irradiance incident at interface k; ς_mk is the scattering coefficient, and thus (1−ς_mk) is the absorptivity; ${\mathit{\beta }}_{mk}^{\odot }$ and β_mk are the backscattered fraction of scattered direct and diffuse irradiance, respectively; $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ is the effective cumulative plant area index, assumed zero at the top of each layer, and increasing downwards (${\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_k$ is the total for layer k); ${\mathit{\mu }}_k^{\odot }$ and ${\stackrel{\mathrm{‾}}{\mathit{\mu }}}_k$ are the inverse of the optical depth per unit of effective plant area index for direct and diffuse radiation, respectively; and ${\dot{Q}}_{mk}^{\mathit{⧫}}$ is the irradiance emitted by a black body at the same temperature as the cohort ($T_{t_k}$).

Figure 4 Schematic of the radiation module for a patch with N_T cohorts, showing the grid arrangement of the irradiance profiles relative to the cohort positions. Index k corresponds to each cohort, or, in the case of $\left({\dot{Q}}_{mk}^{\odot },{\dot{Q}}_{mk}^{\Downarrow },{\dot{Q}}_{mk}^{\Uparrow }\right)$, the interface beneath each cohort; m corresponds to each spectral band; $\left({\dot{Q}}_{m(\mathrm{\infty },a)}^{\odot },{\dot{Q}}_{m(\mathrm{\infty },a)}^{\Downarrow }\right)$ are the incoming direct and diffuse irradiance (Table 3); Z^⊙ is the Sun's zenith angle; $\left({\dot{Q}}_{mk}^{\odot },{\dot{Q}}_{mk}^{\Downarrow },{\dot{Q}}_{mk}^{\Uparrow }\right)$ are the downward direct, downward diffuse, and upward diffuse irradiance (Eqs. 45–47); ${\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_k$ is the effective plant area index (Eq. S100); (ς_mk;β_mk) are the scattering coefficients and the backscattering fraction for diffuse irradiance (Eqs. S101, S102); $\left({\mathit{\varsigma }}_{mk}^{\odot };{\mathit{\beta }}_{mk}^{\odot }\right)$ are the scattering coefficients and the backscattering fraction for direct irradiance (Eqs. S104–S105); ${\dot{Q}}_{mk}^{\mathit{⧫}}$ is the black-body irradiance (Eq. 48); and ${\dot{Q}}_{a,t_k}^{\odot }$ is the net absorbed irradiance (Eq. 49).

Equations (45)–(47) simplify for each spectral band. First, ${\dot{Q}}_{mk}^{\odot }\equiv \mathrm{0}$ for the TIR (m=3) band, because we assume that all incoming TIR irradiance is diffuse. Likewise, the black-body emission ${\dot{Q}}_{mk}^{\mathit{⧫}}=\mathrm{0}$ for the PAR (m=1) and NIR (m=2) bands, because thermal emission is negligible at these wavelengths. The black-body emission for the TIR band is defined as

$$\begin{array}{}\text{(48)}& {\dot{Q}}_{m=\mathrm{3}k}^{\mathit{⧫}}={\mathit{\sigma }}_{\mathrm{SB}}{T}_{t_k}^{\mathrm{4}},\end{array}$$

where σ_SB is the Stefan–Boltzmann constant (Table S3). Note that for emission of TIR radiation (Eqs. 45–47), we assume that emissivity is the same as absorptivity (Kirchhoff's law; Liou, 2002) and hence the (1−ς) term.

The effective plant area index ${\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_k$ is the total area (leaves and branches) that is corrected to account for that leaves are not uniformly distributed in the layer. It is defined as ${\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_k={\mathrm{\Omega }}_k+f_{{\text{Clump}}_k}{\mathrm{\Lambda }}_k$, where $f_{{\text{Clump}}_k}$ is the PFT-dependent clumping index (Chen and Black, 1992, default values in Tables S5–S6), ${\mathrm{\Lambda }}_{t_k}$ is the leaf area index, and ${\mathrm{\Omega }}_{t_k}$ is the wood area index. $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}$ is assumed zero at the top of each layer, increasing downwards.

The optical properties of the leaf layers – optical depth and scattering parameters for direct and diffuse radiation for each of the three spectral bands – are assumed constant within each layer. These properties are determined from PFT-dependent characteristics such as mean orientation factor, spectral-band-dependent reflectivity, transmissivity, and emissivity (Sect. S11). Because the properties are constant within each layer, it is possible to analytically solve the full profile of both direct and diffuse radiation using the solver described in Sect. S12.

Once the profiles of ${\dot{Q}}_{mk}^{\odot }$, ${\dot{Q}}_{mk}^{\Downarrow }$ and ${\dot{Q}}_{mk}^{\Uparrow }$ are determined, we obtain the irradiance that is absorbed by each cohort ${\dot{Q}}_{a,t_k}$:

$$\begin{array}{}\text{(49)}& \begin{array}{rl}& {\dot{Q}}_{a,t_k}=\sum _{m=\mathrm{1}}^{\mathrm{3}}\left[\left({\dot{Q}}_{m(k+\mathrm{1})}^{\odot }-{\dot{Q}}_{mk}^{\odot }\right)+\left({\dot{Q}}_{m(k+\mathrm{1})}^{\Downarrow }-{\dot{Q}}_{mk}^{\Downarrow }\right)\right.\\ & \left.+\left({\dot{Q}}_{mk}^{\Uparrow }-{\dot{Q}}_{m(k+\mathrm{1})}^{\Uparrow }\right)\right].\end{array}\end{array}$$

This term is then used in the enthalpy budget of each cohort (Eqs. 4 and 12).

4.3.2 Ground radiation

The ground radiation submodel determines the irradiance emitted by the ground surface and the profile of irradiance through the temporary surface water layers and top soil layer. Note that the ground radiation and the canopy radiation model are interdependent: the incoming radiation at the top ground layer is determined from the canopy radiation model, and the ground scattering coefficient (ς_m0; see Sect. S11) is needed for the canopy radiation bottom boundary condition (Sect. S12). However, since the scattering coefficient does not depend on the total incoming radiation, the irradiance profile can be solved for a standardized amount of incoming radiation, and once the downward radiation at the bottom of the canopy has been calculated, the absorbed irradiance for each layer can be scaled appropriately.

Black-body emission from the ground (${\dot{Q}}_{m\mathrm{0}}^{\mathit{⧫}}$) is calculated as an area-weighted average of the emissivities of exposed soil and temporary surface water:

$$\begin{array}{}\text{(50)}& \begin{array}{rl}& {\dot{Q}}_{m\mathrm{0}}^{\mathit{⧫}}=\mathrm{0}\text{if }m\in \left(\mathrm{1},\mathrm{2}\right),\\ & \left[\left(\mathrm{1}-f_{\text{TSW}}\right)\left(\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}g}\right)\left({\mathit{\sigma }}_{\text{SB}}{T}_{g_{N_G}}^{\mathrm{4}}\right)+f_{\text{TSW}}\left(\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}s}\right)\left({\mathit{\sigma }}_{\text{SB}}{T}_{s_{N_S}}^{\mathrm{4}}\right)\right]\\ & {\left[\left(\mathrm{1}-f_{\text{TSW}}\right)\left(\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}g}\right)+f_{\text{TSW}}\left(\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}s}\right)\right]}^{-\mathrm{1}}\text{ if }m=\mathrm{3}\end{array}\end{array}$$

where (1−ς_3g) and (1−ς_3s), are, respectively, the thermal-infrared emissivities of the top soil layer and the temporary surface water (Table S4), and f_TSW is the fraction of ground covered by temporary surface water. In ED-2.2, the soil and snow scattering coefficients for the TIR band are assumed constant (Table S4), following Walko et al. (2000).

Once the irradiance profile for the canopy is determined from Eqs. (45) to (47), the irradiance absorbed by each temporary surface water layer ($j\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{S}}\right\}$) is calculated by integrating the transmissivity profile for each layer, starting from the top layer:

$$\begin{array}{}\text{(51)}& \begin{array}{rl}& {\dot{Q}}_{a,{\text{s}}_j}=\\ & \left\{\begin{array}{l}{\sum }_{m=\mathrm{1}}^{\mathrm{2}}\left\{f_{\mathrm{TSW}}\left({\dot{Q}}_{m\mathrm{1}}^{\Downarrow }+{\dot{Q}}_{m\mathrm{1}}^{\odot }\right)\left[\mathrm{1}-\exp \left(-\frac{\mathrm{\Delta }{\stackrel{\mathrm{‾}}{z}}_{s_{N_{\mathrm{S}}}}}{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{s}}}\right)\right]\right\}\\ +f_{\mathrm{TSW}}\left(\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}s}\right)\left({\dot{Q}}_{m=\mathrm{3}k=\mathrm{1}}^{\Downarrow }-{\mathit{\sigma }}_{\mathrm{SB}}{T}_{s_{N_{\mathrm{S}}}}^{\mathrm{4}}\right),\text{if }j=N_{\mathrm{S}}\\ {\sum }_{m=\mathrm{1}}^{\mathrm{2}}\left\{f_{\mathrm{TSW}}\left({\dot{Q}}_{m\mathrm{1}}^{\Downarrow }+{\dot{Q}}_{m\mathrm{1}}^{\odot }\right)\left[\exp \left(-\frac{{\sum }_{j^{\prime }=j+\mathrm{1}}^{N_{\mathrm{S}}}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{z}}_{{\text{s}}_j}}{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{s}}}\right)\right.\right.\\ \left.\left.-\exp \left(-\frac{{\sum }_{j^{\prime }=j}^{N_{\mathrm{S}}}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{z}}_{s_{j^{\prime }}}}{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{s}}}\right)\right]\right\}\text{, otherwise}\end{array}\right.,\end{array}\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{s}}$ is the inverse of the optical depth of temporary surface water.

The irradiance absorbed by the ground is a combination of irradiance of exposed soil and irradiance that is transmitted through all temporary surface water layers and the net absorption of longwave radiation:

$$\begin{array}{}\text{(52)}& \begin{array}{rl}& {\dot{Q}}_{a,g_{N_{\mathrm{G}}}}=\\ & \sum _{m=\mathrm{1}}^{\mathrm{2}}\left\{\left[\mathrm{1}-f_{\mathrm{TSW}}+f_{\mathrm{TSW}}\exp \left(-{\displaystyle \frac{{\sum }_{j^{\prime }=\mathrm{1}}^{N_{\mathrm{S}}}\mathrm{\Delta }{\stackrel{\mathrm{‾}}{z}}_{s_{j^{\prime }}}}{{\stackrel{\mathrm{‾}}{\mathit{\mu }}}_{\mathrm{s}}}}\right)\right]\right.\\ & \left.\left({\dot{Q}}_{m\mathrm{1}}^{\Downarrow }+{\dot{Q}}_{m\mathrm{1}}^{\odot }\right)\right\}+\left(\mathrm{1}-f_{\mathrm{TSW}}\right)\left(\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}g}\right)\\ & \left({\dot{Q}}_{m=\mathrm{3}k=\mathrm{1}}^{\Downarrow }-{\mathit{\sigma }}_{\mathrm{SB}}{T}_{g_{N_{\mathrm{G}}}}^{\mathrm{4}}\right)\end{array}.\end{array}$$

4.4 Surface layer model

The surface layer model determines the fluxes of enthalpy, water, and carbon dioxide between the canopy air space and the free air above. It is based on the Monin–Obukhov similarity theory (Monin and Obukhov, 1954; Foken, 2006), which has been widely used by biosphere–atmosphere models representing a variety of biomes (e.g., Walko et al., 2000; Best et al., 2011; Oleson et al., 2013), although this is often an extrapolation of the theory that was not originally developed for heterogeneous vegetation or tall vegetation (Foken, 2006).

In order to obtain the fluxes, we assume that the eddy diffusivity of buoyancy is the same as the diffusivity of enthalpy, water vapor, and CO₂. This assumption allows us to define a single canopy conductance G_c for the three variables, following the algorithm described in Sects. S13 and S14.1. We then obtain the following equations for fluxes between canopy air space and the free atmosphere:

$$\begin{array}{}\text{(53)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{\mathrm{a},\mathrm{c}}={\mathit{\rho }}_{\mathrm{c}}{G}_{\mathrm{c}}\left(w_{\mathrm{a}}-w_{\mathrm{c}}\right),\text{(54)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{\mathrm{a},\mathrm{c}}={\mathit{\rho }}_{\mathrm{c}}{G}_{\mathrm{c}}\left({\stackrel{\mathrm{̃}}{h}}_{\mathrm{a}}-h_{\mathrm{c}}\right),\text{(55)}& {\displaystyle }& {\displaystyle }{\dot{C}}_{\mathrm{a},\mathrm{c}}={\displaystyle \frac{{\mathcal{M}}_{\mathrm{C}}}{{\mathcal{M}}_{\mathrm{d}}}}{\mathit{\rho }}_{\mathrm{c}}{G}_{\mathrm{c}}\left(c_{\mathrm{a}}-c_{\mathrm{c}}\right),\end{array}$$

where ${\stackrel{\mathrm{̃}}{h}}_{\mathrm{a}}$ is the equivalent enthalpy of air at reference height z_a when the air is adiabatically moved to the top of the canopy air space, using the definition of potential temperature:

$$\begin{array}{}\text{(56)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{h}}_{\mathrm{a}}=h\left({\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}},w_{\mathrm{a}}\right),\text{ from Eq. (S50)},\text{(57)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}}={\mathit{\theta }}_{\mathrm{a}}{\left({\displaystyle \frac{p_{\mathrm{c}}}{p_{\mathrm{0}}}}\right)}^{\frac{\mathcal{R}}{{\mathcal{M}}_{\mathrm{d}}{q}_{\mathrm{pd}}}},\end{array}$$

where p₀ is the reference pressure, ℛ is the universal gas constant, q_pd is the specific heat of dry air at constant pressure, and ℳ_d is the molar mass of dry air (Table S3).

Sensible heat flux between the free atmosphere and canopy air space (${\dot{Q}}_{\mathrm{a},\mathrm{c}}$) can be derived from the definition of enthalpy and enthalpy flux (Eqs. S50 and 54), although it is not directly applied to the energy balance in the canopy air space (${\dot{H}}_{\mathrm{a},\mathrm{c}}$ is used instead):

$$\begin{array}{}\text{(58)}& \begin{array}{rl}& {\dot{H}}_{\mathrm{a},\mathrm{c}}={\mathit{\rho }}_{\mathrm{c}}{G}_{\mathrm{c}}\left[\left(\mathrm{1}-w_{\mathrm{a}}\right)q_{\mathrm{pd}}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}}+w_{\mathrm{a}}{q}_{\mathrm{pv}}\left({\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}}-T_{v\mathrm{0}}\right)\right.\\ & \left.-\left(\mathrm{1}-w_{\mathrm{c}}\right)q_{\mathrm{pd}}{T}_{\mathrm{c}}-w_{\mathrm{c}}{q}_{\mathrm{pv}}\left(T_{\mathrm{c}}-T_{v\mathrm{0}}\right)\right]\\ & =\underset{{\dot{Q}}_{\mathrm{a},\mathrm{c}}}{\underbrace{{\mathit{\rho }}_{\mathrm{c}}{G}_{\mathrm{c}}\left(q_{p_{\mathrm{a}}}{\stackrel{\mathrm{̃}}{T}}_{\mathrm{a}}-q_{p_{\mathrm{c}}}{T}_{\mathrm{c}}\right)}}-{\mathit{\rho }}_{\mathrm{c}}{G}_{\mathrm{c}}\left(w_{\mathrm{a}}-w_{\mathrm{c}}\right)q_{\mathrm{pv}}{T}_{v\mathrm{0}},\end{array}\end{array}$$

$$\begin{array}{}\text{(59)}& {\displaystyle }{\dot{Q}}_{\mathrm{a},\mathrm{c}}={\dot{H}}_{\mathrm{a},\mathrm{c}}+{\dot{W}}_{\mathrm{a},\mathrm{c}}{q}_{\mathrm{pv}}{T}_{v\mathrm{0}}.\end{array}$$

4.5 Heat and water exchange between surfaces and canopy air space

4.5.1 Leaves and branches

Fluxes of sensible heat (${\dot{Q}}_{t_k,c}$) and water vapor (${\dot{W}}_{t_k,c}$) between the leaf surface and wood surface and the canopy air space follow the same principle of conductance and gradient that define the eddy fluxes between the free atmosphere and canopy air space (Eqs. 53, 54). Throughout this section, we use subscripts λ_k and β_k to denote leaf and wood boundary layers of cohort k, respectively; the different subscripts are needed to differentiate fluxes coming from the leaves' intercellular space (e.g., transpiration; see also Sect. 4.6). Let $G_{Q{\mathit{\lambda }}_k}$ (m s⁻¹) and $G_{W{\mathit{\lambda }}_k}$ (m s⁻¹) be the conductances of heat and water between the leaf boundary layer of cohort k and the canopy air space, and $G_{Q{\mathit{\beta }}_k}$ and $G_{W{\mathit{\beta }}_k}$ be the wood boundary layer counterparts. The surface sensible heat and surface water vapor fluxes are

$$\begin{array}{}\text{(60)}& {\displaystyle }& {\displaystyle }{\dot{Q}}_{t_k,c}={\dot{Q}}_{{\mathit{\lambda }}_k,c}+{\dot{Q}}_{{\mathit{\beta }}_k,c}=\mathrm{2}{\mathrm{\Lambda }}_k{\dot{q}}_{{\mathit{\lambda }}_k,c}+\mathit{\pi }{\mathrm{\Omega }}_k{\dot{q}}_{{\mathit{\beta }}_k,c},\text{(61)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{t_k,c}={\dot{W}}_{{\mathit{\lambda }}_k,c}+{\dot{W}}_{{\mathit{\beta }}_k,c}={\mathrm{\Lambda }}_k{\dot{w}}_{{\mathit{\lambda }}_k,c}+{\mathrm{\Omega }}_k{\dot{w}}_{{\mathit{\beta }}_k,c},\text{(62)}& {\displaystyle }& {\displaystyle }{\dot{q}}_{{\mathit{\lambda }}_k,c}=G_{Q{\mathit{\lambda }}_k}{\mathit{\rho }}_{\mathrm{c}}{q}_{p_{\mathrm{c}}}\left(T_{l_k}^{\mathrm{Sfc}}-T_{\mathrm{c}}\right),\text{(63)}& {\displaystyle }& {\displaystyle }{\dot{q}}_{{\mathit{\beta }}_k,c}=G_{Q{\mathit{\beta }}_k}{\mathit{\rho }}_{\mathrm{c}}{q}_{p_{\mathrm{c}}}\left(T_{b_k}^{\mathrm{Sfc}}-T_{\mathrm{c}}\right),\text{(64)}& {\displaystyle }& {\displaystyle }{\dot{w}}_{{\mathit{\lambda }}_k,c}=G_{W{\mathit{\lambda }}_k}{\mathit{\rho }}_{\mathrm{c}}\left(w_{l_k}^{\mathrm{Sfc}}-w_{\mathrm{c}}\right),\text{(65)}& {\displaystyle }& {\displaystyle }{\dot{w}}_{{\mathit{\beta }}_k,c}=G_{W{\mathit{\beta }}_k}{\mathit{\rho }}_{\mathrm{c}}\left(w_{b_k}^{\mathrm{Sfc}}-w_{\mathrm{c}}\right),\end{array}$$

where (${\dot{q}}_{{\mathit{\lambda }}_k,c};{\dot{q}}_{{\mathit{\beta }}_k,c};{\dot{w}}_{{\mathit{\lambda }}_k,c};{\dot{w}}_{{\mathit{\beta }}_k,c}$) are the leaf surface and branch surface heat and water fluxes by unit of leaf and branch area, respectively; the factors 2 and π in Eq. (60) mean that sensible heat is exchanged on both sides of the leaves and on the longitudinal area of the branches, which are assumed cylindrical. Intercepted water and dew and frost formation is allowed only on one side of the leaves, and an area equivalent to a one-sided flat plate for branches, and therefore only the leaf and wood area indices are used in Eq. (61). Canopy air space temperature, specific humidity, density, and specific heat, leaf temperature, and wood temperature are determined diagnostically. We also assume surface temperature of leaves and branches to be the same as their internal temperatures (i.e., $T_{l_k}^{\mathrm{Sfc}}\equiv T_{l_k}$ and $T_{b_k}^{\mathrm{Sfc}}\equiv T_{b_k}$). Specific humidity at the leaf surface $w_{l_k}^{\mathrm{Sfc}}=w_{\text{Sat}}\left(T_{l_k}^{\mathrm{Sfc}},p_{\mathrm{c}}\right)$ and branch surface $w_{{\mathit{\beta }}_k}^{\mathrm{Sfc}}=w_{\text{Sat}}\left(T_{b_k}^{\mathrm{Sfc}},p_{\mathrm{c}}\right)$ are assumed to be the saturation specific humidity w_Sat (Sect. S15).

Heat conductance for leaves and branches is based on the convective heat transfer, as described in Sect. S14.2. Further description of the theory can be found in Monteith and Unsworth (2008, Sect. 10.1).

4.5.2 Temporary surface water and soil

Sensible heat and water fluxes between the temporary surface water and soil and the canopy air space are calculated similarly to leaves and branches. Surface conductance G_Sfc is assumed to be the same for both heat and water, and also the same for soil and temporary surface water:

$$\begin{array}{}\text{(66)}& {\displaystyle }& {\displaystyle }{\dot{Q}}_{s_{N_{\mathrm{S}}},c}=f_{\text{TSW}}{G}_{\text{Sfc}}{\mathit{\rho }}_{\mathrm{c}}{q}_{p_{\mathrm{c}}}\left(T_{s_{N_{\mathrm{S}}}}-T_{\mathrm{c}}\right),\text{(67)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{s_{N_{\mathrm{S}}},c}=f_{\text{TSW}}{G}_{\text{Sfc}}{\mathit{\rho }}_{\mathrm{c}}\left(w_{s_{N_{\mathrm{S}}}}-w_{\mathrm{c}}\right),\text{(68)}& {\displaystyle }& {\displaystyle }{\dot{Q}}_{g_{N_{\mathrm{G}}},c}=\left(\mathrm{1}-f_{\text{TSW}}\right)G_{\text{Sfc}}{\mathit{\rho }}_{\mathrm{c}}{q}_{p_{\mathrm{c}}}\left(T_{g_{N_{\mathrm{G}}}}-T_{\mathrm{c}}\right),\text{(69)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{g_{N_{\mathrm{G}}},c}=\left(\mathrm{1}-f_{\text{TSW}}\right)G_{\text{Sfc}}{\mathit{\rho }}_{\mathrm{c}}\left(w_{g_{N_{\mathrm{G}}}}-w_{\mathrm{c}}\right).\end{array}$$

Specific humidity for temporary surface water is computed exactly as leaves and branches, $w_{s_{N_{\mathrm{S}}}}=w_{\text{Sat}}\left(T_{s_{N_{\mathrm{S}}}},p_{\mathrm{c}}\right)$ (Sect. S15). For soils, the specific humidity also accounts for the soil moisture and the sign of the flux, using a method similar to Avissar and Mahrer (1988):

$$\begin{array}{}\text{(70)}& w_{g_{N_{\mathrm{G}}}}=\left\{\begin{array}{l}{s}_g\exp \left(\frac{{\mathcal{M}}_{w}g{\mathrm{\Psi }}_{g_{N_{\mathrm{G}}}}}{\mathcal{R}{T}_{g_{N_{\mathrm{G}}}}}\right)w_{\text{Sat}}(T_{g_{N_{\mathrm{G}}}},p_{\mathrm{c}})+\left(\mathrm{1}-s_g\right)w_{\mathrm{c}},\\ \mathrm{\text{ if }}{w}_{\text{Sat}}(T_{g_{N_{\mathrm{G}}}},p_{\mathrm{c}})>w_{\mathrm{c}}\\ w_{\text{Sat}}(T_{g_{N_{\mathrm{G}}}},p_{\mathrm{c}}),\\ \mathrm{\text{ if }}{w}_{\text{Sat}}(T_{g_{N_{\mathrm{G}}}},p_{\mathrm{c}})\le w_{\mathrm{c}}\end{array},\right.\end{array}$$

$$\begin{array}{}\text{(71)}& {\displaystyle }{\displaystyle }{s}_g={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\mathrm{1.0}-\cos \left[\mathit{\pi }{\displaystyle \frac{min\left({\mathit{\vartheta }}_{g_{N_{\mathrm{G}}}},{\mathit{\vartheta }}_{\text{Fc}}\right)-{\mathit{\vartheta }}_{\text{Re}}}{{\mathit{\vartheta }}_{\text{Fc}}-{\mathit{\vartheta }}_{\text{Re}}}}\right]\right\},\end{array}$$

where g is the gravity acceleration, ℳ_w is the water molar mass, and ℛ is the universal gas constant (Table S3); $T_{g_{N_{\mathrm{G}}}}$, ${\mathit{\vartheta }}_{g_{N_{\mathrm{G}}}}$, and ${\mathrm{\Psi }}_{g_{N_{\mathrm{G}}}}$ are the temperature, soil moisture, and soil matric potential of the topmost soil layer, respectively; and ϑ_Fc and ϑ_Re are the soil moisture at field capacity and the residual soil moisture, respectively. The exponential term in Eq. (70) corresponds to the soil pore relative humidity derived from the Kelvin equation (Philip, 1957), and s_g is the soil wetness function, which takes a similar functional form as the relative humidity term from Noilhan and Planton (1989) and the β term from Lee and Pielke (1992). The total resistance between the surface and the canopy air space is a combination of the resistance if the surface was bare, plus the resistance due to the vegetation, as described in Sect. S14.3.

4.5.3 Enthalpy flux due to evaporation and condensation

Dew and frost are formed when water in the canopy air space condenses or freezes on any surface (leaves, branches, or ground); likewise, water that evaporates and ice that sublimates from these surfaces immediately become part of the canopy air space. In terms of energy transfer, two processes occur, the phase change and the mass exchange, and both must be accounted for the enthalpy flux. Phase change depends on the specific latent heat of vaporization (l_lv) and sublimation (l_iv), which are linear functions of temperature, based on Eqs. (S48) and (S49):

$$\begin{array}{}\text{(72)}& {\displaystyle }& {\displaystyle }{l}_{lv}\left(T\right)=l_{lv\mathrm{3}}+\left(q_{pv}-q_{\mathrm{\ell }}\right)\left(T-T_{\mathrm{3}}\right),\text{(73)}& {\displaystyle }& {\displaystyle }{l}_{iv}\left(T\right)=l_{iv\mathrm{3}}+\left(q_{pv}-q_i\right)\left(T-T_{\mathrm{3}}\right),\end{array}$$

where l_lv3 and l_iv3 are the specific latent heats of vaporization and sublimation at the water triple point (T₃), q_pv is the specific heat of water vapor at constant pressure, and q_i and q_ℓ are the specific heats of ice and liquid water, respectively (Table S3). The temperature for phase change must be the surface temperature because this is where the phase change occurs. In the most generic case, if a surface x at temperature T_x has a liquid water fraction ℓ_x, the total enthalpy flux between the surface and canopy air space ${\dot{H}}_{x,c}$ associated with the water flux W_x,c is

$$\begin{array}{}\text{(74)}& \begin{array}{rl}{\dot{H}}_{x,c}& ={\dot{W}}_{x,c}\left\{\underset{\text{Enthalpy flux due to mass exchange}}{\underbrace{\left[\left(\mathrm{1}-{\mathrm{\ell }}_x\right)q_i{T}_x+{\mathrm{\ell }}_x{q}_{\mathrm{\ell }}\left(T_x-T_{l\mathrm{0}}\right)\right]}}\right.\\ & \left.+\underset{\text{Enthalpy flux due to phase change}}{\underbrace{\left[\left(\mathrm{1}-{\mathrm{\ell }}_x\right)l_{iv}\left(T_x\right)+{\mathrm{\ell }}_x{l}_{lv}\left(T_x\right)\right]}}\right\}.\end{array}\end{array}$$

By using the definitions from Eq. (S54), Eq. (74) can be further simplified to

$$\begin{array}{}\text{(75)}& {\dot{H}}_{x,c}={\dot{W}}_{x,c}\left[q_{pv}\left(T_x-T_{v\mathrm{0}}\right)\right]={\dot{W}}_{x,c}\underset{\text{Eq. (S50)}}{\underbrace{h\left(T_x,w_x=\mathrm{1}\right)}},\end{array}$$

which is consistent with the exchange of pure water vapor and enthalpy between the thermodynamic systems. Equation (75) is used to determine ${\dot{H}}_{g_{N_{\mathrm{G}}},c}$, ${\dot{H}}_{s_{N_{\mathrm{S}}},c}$, and ${\dot{H}}_{t_k,c},k\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_T\right\}$.

4.6 Leaf physiology

In ED-2.2, leaf physiology is modeled following Farquhar et al. (1980) and Collatz et al. (1991) for C₃ plants; Collatz et al. (1992) for C₄ plants; and the Leuning (1995) model for stomatal conductance. This submodel ultimately determines the net leaf-level CO₂ uptake rate of each cohort k (${\dot{A}}_k$, ${\mathrm{mol}}_{\mathrm{C}}{\mathrm{m}}_{\mathrm{Leaf}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), controlled exclusively by the leaf environment, and the corresponding water loss through transpiration (${\dot{E}}_k$, ${\mathrm{mol}}_{\mathrm{W}}{\mathrm{m}}_{\mathrm{Leaf}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

Figure 5 Schematic of fluxes between a leaf and the surrounding canopy air space for a hypostomatous plant during the photo period, as represented in ED-2.2. Conductances are represented by the resistances between the different environments (G⁻¹). Leaf-level sensible heat flux (${\dot{q}}_{{\mathit{\lambda }}_k,c}$; Eq. 60) and leaf-level vapor flux between intercepted water and canopy air space (${\dot{w}}_{{\mathit{\lambda }}_k,c}$; Eq. 61) are also shown for comparison. Cohort index k is omitted from the figure for clarity.

The exchange of water and CO₂ between the leaf intercellular space and the canopy air space is mediated by the stomata and the leaf boundary layer, which imposes an additional resistance to fluxes of these substances. For simplicity, we assume that the leaf boundary layer air has low storage capacity, and thus the fluxes of any substance (water or CO₂) entering and exiting the boundary layer must be the same. Fluxes of water and carbon between the leaf intercellular space and the canopy air space must overcome both the stomatal resistance and the boundary layer resistance, whereas sensible heat flux and water flux from leaf surface water must overcome the boundary layer resistance only (Fig. 5). When soil moisture is not a limiting factor, the fluxes of CO₂ and water can be written as

$$\begin{array}{}\text{(76)}& \begin{array}{rl}{\dot{A}}_k& ={\widehat{G}}_{C{\mathit{\lambda }}_k}\left(c_{\mathrm{c}}-c_{{\mathit{\lambda }}_k}\right)={\widehat{G}}_{C{l}_k}\left(c_{{\mathit{\lambda }}_k}-c_{l_k}\right)\\ & ={\displaystyle \frac{{\widehat{G}}_{C{\mathit{\lambda }}_k}{\widehat{G}}_{C{l}_k}}{{\widehat{G}}_{C{\mathit{\lambda }}_k}+{\widehat{G}}_{C{l}_k}}}\left(c_{\mathrm{c}}-c_{l_k}\right),\end{array}\end{array}$$

$$\begin{array}{}\text{(77)}& \begin{array}{rl}{\dot{E}}_k& ={\widehat{G}}_{W{\mathit{\lambda }}_k}\left(w_{\mathrm{c}}-w_{{\mathit{\lambda }}_k}\right)={\widehat{G}}_{W{l}_k}\left(w_{{\mathit{\lambda }}_k}-w_{l_k}\right)\\ & ={\displaystyle \frac{{\widehat{G}}_{W{\mathit{\lambda }}_k}{\widehat{G}}_{W{l}_k}}{{\widehat{G}}_{W{\mathit{\lambda }}_k}+{\widehat{G}}_{W{l}_k}}}\left(w_{\mathrm{c}}-w_{l_k}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(78)}& w_{l_k}=w_{\text{Sat}}\left(T_{t_k},p_{\mathrm{c}}\right)\text{(Sect. S15)},\end{array}$$

and

$$\begin{array}{}\text{(79)}& {\widehat{G}}_{X{\mathit{\lambda }}_k}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{c}}{G}_{X{\mathit{\lambda }}_k}}{{\mathcal{M}}_d}},\end{array}$$

and

$$\begin{array}{}\text{(80)}& {\widehat{G}}_{X{l}_k}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{c}}{G}_{X{l}_k}}{{\mathcal{M}}_d}},\end{array}$$

where $G_{X{\mathit{\lambda }}_k}$ and $G_{X{l}_k}$ (units m s⁻¹) are the leaf boundary layer and stomatal conductances for element X (either water W or carbon C), respectively; $c_{{\mathit{\lambda }}_k}$ and $w_{{\mathit{\lambda }}_k}$ are the CO₂ mixing ratio and the specific humidity of the leaf boundary layer, respectively; $c_{l_k}$ and $w_{l_k}$ are the CO₂ and specific humidity of the leaf intercellular space, respectively; ℳ_d is the molar mass of dry air; and ρ_c is the air density. As stated in Eq. (78), we assume the leaf intercellular space to be at water vapor saturation. The leaf boundary layer conductances are obtained following the algorithm shown in Sect. S14.2. The net CO₂ assimilation flux and stomatal conductances are described below.

From Farquhar et al. (1980), the net CO₂ assimilation flux is defined as

$$\begin{array}{}\text{(81)}& {\dot{A}}_k=\underset{\text{Carboxylation}}{\underbrace{{\dot{V}}_{C_k}}}-\underset{\begin{array}{c}\text{Oxygenation}\\ \left(\text{photorespiration}\right)\end{array}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\dot{V}}_{O_k}}}-\underset{\text{Leaf respiration}}{\underbrace{{\dot{R}}_k}}.\end{array}$$

Oxygenation releases 0.5 ${\mathrm{mol}}_{{\mathrm{CO}}_{\mathrm{2}}}$ for every ${\mathrm{mol}}_{{\mathrm{O}}_{\mathrm{2}}}$, hence the half multiplier, and it is related to carboxylation by means of the CO₂ compensation point Γ_k (Lambers et al., 2008):

$$\begin{array}{}\text{(82)}& {\dot{V}}_{O_k}={\displaystyle \frac{\mathrm{2}{\mathrm{\Gamma }}_k}{c_{l_k}}}{\dot{V}}_{C_k},\end{array}$$

where $c_{l_k}$ is the CO₂ mixing ratio in the leaf intercellular space. The CO₂ compensation point is determined after Collatz et al. (1991, 1992):

$$\begin{array}{}\text{(83)}& {\mathrm{\Gamma }}_k=\left\{\begin{array}{ll}{\displaystyle \frac{o_{\oplus }}{\mathrm{2}\mathit{\varpi }}},& \text{in case cohort}k\text{is a}{\mathrm{C}}_{\mathrm{3}}\text{plant}\\ \mathrm{0},& \text{in case cohort}k\text{is a}{\mathrm{C}}_{\mathrm{4}}\text{plant}\end{array}\right.,\end{array}$$

where o_⊕ is the reference O₂ mixing ratio (Table S3), and ϖ represents the ratio between the rates of carboxylase to oxygenase and is a function of temperature. The general form of the function 𝒯 describing the metabolic dependence upon temperature for any variable x (including ϖ) is

$$\begin{array}{}\text{(84)}& \mathcal{T}\left(T,x\right)=x_{\mathrm{15}}\times {\mathcal{Q}}_{{\mathrm{10}}_x}^{\frac{T-T_{\mathrm{15}}}{\mathrm{10}}},\end{array}$$

where x₁₅ is the value of variable x at temperature T₁₅=288.15 K (15 ^∘C), and ${\mathcal{Q}}_{{\mathrm{10}}_x}$ is the parameter which describes temperature dependence (Table S4).

Because C₄ plants have a mechanism to concentrate CO₂ near the CO₂-fixing enzyme RuBisCO (ribulose-1,5-biphosphate carboxylase oxygenase), photorespiration is nearly nonexistent in C₄ plants (Lambers et al., 2008) and hence the assumption that Γ_k is zero. For C₄ plants, the carboxylation rate under ribulose-1,5-biphosphate (RuBP) saturated conditions becomes the maximum capacity of RuBisCO to perform the carboxylase function (${\dot{V}}_{C_k}={\dot{V}}_{C_k}^{\text{max}}$). For C₃, this rate is unattainable even under RuBP-saturated conditions because carboxylation and oxygenation are mutually inhibitive reactions (Lambers et al., 2008). Therefore, the maximum attainable carboxylation (${\dot{V}}_{C_k}={\dot{V}}_{C_k}^{\text{RuBP}}$) is expressed by a modified Michaelis–Menten kinetics equation:

$$\begin{array}{}\text{(85)}& {\dot{V}}_{C_k}^{\text{RuBP}}=\left\{\begin{array}{ll}{\dot{V}}_{C_k}^{\text{max}}\frac{c_{l_k}}{c_{l_k}+{\mathcal{K}}_{{\mathrm{ME}}_k}},& \text{if cohort}k\text{is a}{\mathrm{C}}_{\mathrm{3}}\text{plant}\\ {\dot{V}}_{C_k}^{\text{max}},& \text{if cohort}k\text{is a}{\mathrm{C}}_{\mathrm{4}}\text{plant}\end{array}\right.,\end{array}$$

where ${\mathcal{K}}_{{\mathrm{ME}}_k}={\mathcal{K}}_{C_k}\left(\mathrm{1}+o_{\oplus }/{\mathcal{K}}_{O_k}\right)$ is the effective Michaelis constant, and ${\mathcal{K}}_{C_k}$ and ${\mathcal{K}}_{O_k}$ are the Michaelis constants for carboxylation and oxygenation, respectively. Both ${\mathcal{K}}_{C_k}$ and ${\mathcal{K}}_{O_k}$ are dependent on temperature, following Eq. (84) (default parameters in Table S4), whereas ${\dot{V}}_{C_k}^{\text{max}}$ follows a modified temperature-dependent function to account for the fast decline of both productivity and respiration at low and high temperatures (Sellers et al., 1996; Moorcroft et al., 2001):

$$\begin{array}{}\text{(86)}& \begin{array}{rl}& {\mathcal{T}}^{\prime }\left(T,x\right)\\ & ={\displaystyle \frac{\mathcal{T}\left(T,x\right)}{\left\{\mathrm{1}+\exp \left[-f_{\text{Cold}}\left(T-T_{\text{Cold}}\right)\right]\right\}\left\{\mathrm{1}+\exp \left[+f_{\text{Hot}}\left(T-T_{\text{Hot}}\right)\right]\right\}}},\end{array}\end{array}$$

where f_Cold, f_Hot, T_Cold, and T_Hot are PFT-dependent phenomenological parameters to reduce the function value at low and high temperatures (Tables S5–S6).

The original expression for the initial slope of the carboxylation rate under near-zero CO₂ (${\dot{V}}_{C_k}^{\mathrm{InSl}}$) for C₄ plants by Collatz et al. (1992) has been modified later (e.g., Foley et al., 1996) to explicitly include ${\dot{V}}_{C_k}^{\text{max}}$; this is the same expression used in ED-2.2:

$$\begin{array}{}\text{(87)}& {\dot{V}}_{C_k}^{\mathrm{InSl}}=k_{\mathrm{PEP}}{\dot{V}}_{C_k}^{\text{max}}{c}_{l_k},\end{array}$$

where k_PEP represents the initial slope of the response curve to increasing CO₂; the default value in ED-2.2 (Table S4) is the same value used by Foley et al. (1996).

From the total photosynthetically active irradiance absorbed by the cohort ${\dot{Q}}_{\mathrm{PAR}:a,t_k}$ (Eq. 49), we define the photon flux that is absorbed by the leaf (${\dot{q}}_k^{\mathrm{PAR}},\mathrm{mol}{\mathrm{m}}_{\mathrm{Leaf}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$):

$$\begin{array}{}\text{(88)}& {\dot{q}}_k^{\mathrm{PAR}}={\displaystyle \frac{\mathrm{1}}{\mathrm{Ein}}}{\displaystyle \frac{f_{{\text{Clump}}_k}}{{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_k}}{\dot{Q}}_{\mathrm{PAR}:a,t_k},\end{array}$$

where Ein is the average photon-specific energy in the PAR band (0.4–0.7 µm; Table S3). Even though a high fraction ${\mathit{ϵ}}_k^{\star }$ of the absorbed irradiance is used to transport electrons needed by the light reactions of photosynthesis (Lambers et al., 2008), only a fraction of the irradiance absorbed by the leaf is absorbed by the chlorophyll; in addition, the number electrons needed by each carboxylation and oxygenation reaction poses an additional restriction to the total carboxylation rate. The product of these three factors is combined into a single scaling factor for total absorbed PAR, the quantum yield (ϵ_k), which is a PFT-dependent property in ED-2.2 (Tables S5–S6). The maximum carboxylation rate under light limitation ${\dot{V}}_{C_k}^{\text{PAR}}$ is

$$\begin{array}{}\text{(89)}& {\dot{V}}_{C_k}^{\text{PAR}}=\left\{\begin{array}{l}{\mathit{ϵ}}_k{\dot{q}}_k^{\mathrm{PAR}}\frac{c_{l_k}}{c_{l_k}+\mathrm{2}{\mathrm{\Gamma }}_k},\\ \mathrm{\text{ if cohort}}k\mathrm{is}\mathrm{a}{\mathrm{C}}_{\mathrm{3}}\mathrm{plant}\\ {\mathit{ϵ}}_k{\dot{q}}_k^{\mathrm{PAR}},\\ \mathrm{\text{ if cohort}}k\mathrm{is}\mathrm{a}{\mathrm{C}}_{\mathrm{4}}\mathrm{plant}\end{array}\right..\end{array}$$

Carboxylation may also be limited by the export rate of starch and sucrose that is synthesized by triose phosphate, especially when CO₂ concentration and irradiance are simultaneously high, at low temperatures, or O₂ concentration is low (von Caemmerer, 2000; Lombardozzi et al., 2018). This limitation is not included in ED-2.2.

The leaf respiration term ${\dot{R}}_k$ comprises all leaf respiration terms that are not dependent on photosynthesis, and is primarily constituted of mitochondrial respiration; it is currently represented as a function of the maximum carboxylation rate, following Foley et al. (1996):

$$\begin{array}{}\text{(90)}& {\dot{R}}_k=f_R{\dot{V}}_{C_k}^{max},\end{array}$$

where f_R is a PFT-dependent parameter (Tables S5–S6).

A plant's stomatal conductance is a result of a trade-off between the amount of carbon that leaves uptake and the amount of water that they lose. Leuning (1995) proposed a semi-empirical stomatal conductance expression for water based on these trade-offs:

$$\begin{array}{}\text{(91)}& {\widehat{G}}_{W{l}_k}=\left\{\begin{array}{ll}{\widehat{G}}_{W{l}_k}^{\mathit{\varnothing }}+\frac{M_k{\dot{A}}_k}{\left(c_{{\mathit{\lambda }}_k}-{\mathrm{\Gamma }}_k\right)\left(\mathrm{1}+\frac{w_{l_k}-w_{{\mathit{\lambda }}_k}}{\mathrm{\Delta }{w}_k}\right)},& \text{ if }{\dot{A}}_k>\mathrm{0}\\ {\widehat{G}}_{W{l}_k}^{\mathit{\varnothing }},& \text{ if }{\dot{A}}_k\le \mathrm{0}\end{array}\right.,\end{array}$$

where ${\widehat{G}}_{W{l}_k}^{\mathit{\varnothing }}$ is the residual conductance when stomata are closed, M_k is the slope of the stomatal conductance function, and Δw_k is an empirical coefficient controlling conductance under severe leaf-level water deficit; ${\widehat{G}}_{W{l}_k}^{\mathit{\varnothing }}$, M_k, and Δw_k are PFT-dependent parameters (Tables S5–S6). From Cowan and Troughton (1971), stomatal conductance of CO₂ is estimated by the ratio f_Gl between the diffusivities of water and CO₂ in the air (Table S4):

$$\begin{array}{}\text{(92)}& {\widehat{G}}_{W{l}_k}=f_{Gl}{\widehat{G}}_{C{l}_k}.\end{array}$$

Variables $w_{l_k}$, ${\dot{V}}_{C_k}^{\text{max}}$, ${\dot{R}}_k$, ϖ_k, ${\mathcal{K}}_{O_k}$, ${\mathcal{K}}_{C_k}$, Γ_k, and ${\mathcal{K}}_{{\mathrm{ME}}_k}$ are functions of leaf temperature and canopy air space pressure and thus can be determined directly. In contrast, nine variables are unknown for each limitation case, namely the RuBP-saturated (Eq. 85), CO₂-limited (Eq. 87), and light-limited (Eq. 89) variables, as well as for the case when the stomata are closed (Eq. 91 for when ${\dot{A}}_k\le \mathrm{0}$): ${\dot{E}}_k$, ${\dot{A}}_k$, ${\dot{V}}_{C_k}$, ${\dot{V}}_{O_k}$, $c_{l_k}$, $c_{{\mathit{\lambda }}_k}$, $w_{{\mathit{\lambda }}_k}$, ${\widehat{G}}_{W{l}_k}$, and ${\widehat{G}}_{C{l}_k}$. These are determined numerically for each limiting case, and the value of ${\dot{A}}_k$ is taken to be the limiting case that yields the lowest ${\dot{A}}_k$, following the algorithm described in Sect. S16.

The stomatal conductance model by Leuning (1995) (Eq. 91) is regulated by leaf vapor pressure deficit. However, Eqs. (76) and (77) do not account for soil moisture limitation of photosynthesis. To represent this effect, we define a soil-moisture-dependent scaling factor ($f_{W{l}_k}$; Sect. S17) to reduce productivity and transpiration as soil available water decreases. Because stomatal conductance cannot be zero, the scaling factor $f_{W{l}_k}$ interpolates between the fully closed case (${\dot{A}}_k^{\mathit{\varnothing }};{\dot{E}}_k^{\mathit{\varnothing }}$) and the solution without soil moisture limitation (${\dot{A}}_k;{\dot{E}}_k$), yielding to the actual fluxes of CO₂ (${\dot{C}}_{l_k,c}$, ${\mathrm{kg}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and water (${\dot{W}}_{l_k,c}$, ${\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$):

$$\begin{array}{}\text{(93)}& {\displaystyle }& {\displaystyle }{\dot{C}}_{l_k,c}=-{\mathit{þ}}_k{\mathcal{M}}_{\mathrm{C}}{\mathrm{\Lambda }}_k\left[\left(\mathrm{1}-f_{W{l}_k}\right){\dot{A}}_k^{\mathit{\varnothing }}+f_{W{l}_k}{\dot{A}}_k\right],\text{(94)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{l_k,c}={\mathit{þ}}_k{\mathcal{M}}_w{\mathrm{\Lambda }}_k\left[\left(\mathrm{1}-f_{W{l}_k}\right){\dot{E}}_k^{\mathit{\varnothing }}+f_{W{l}_k}{\dot{E}}_k\right],\end{array}$$

where þ_k is 1 if the PFT is hypostomatous or 2 if the PFT is amphistomatous or needleleaf (Tables S5–S6). Alternatively, Xu et al. (2016) implemented a process-based plant hydraulics scheme that solves the soil–stem–leaf water flow in ED-2.2; details of this implementation are available in the above-mentioned paper.

For simplicity, we assume that the water content in the leaf intercellular space and the plant vascular system is constant; therefore, the amount of water lost by the intercellular space through transpiration always matches the amount of water absorbed by roots. Plants may extract water from all layers to which they have access, and the amount of water extracted from each layer is proportional to the available water in the layer relative to the total available water ($W_{{\text{g}}_j}^{\star }$):

$$\begin{array}{}\text{(95)}& {\displaystyle }& {\displaystyle }\sum _{j=j_{\mathrm{0}k}}^{N_{\mathrm{G}}}{\dot{W}}_{{\text{g}}_j,l_k}={\dot{W}}_{l_k,c},\text{(96)}& {\displaystyle }& {\displaystyle }{\dot{W}}_{{\text{g}}_j,l_k}={\dot{W}}_{l_k,c}{\displaystyle \frac{W_{{\text{g}}_j}^{\star }-W_{g_{j+\mathrm{1}}}^{\star }}{W_{g_{j\mathrm{0}}}^{\star }}},\end{array}$$

where $W_{{\text{g}}_j}^{\star }$ is defined following Sect. S17 and $W_{g(N_{\mathrm{G}}+\mathrm{1})}^{\star }\equiv \mathrm{0}$. The net water flux in the leaf intercellular space due to transpiration is assumed to be zero; however, the associated net energy flux is not: water enters the leaf intercellular space as liquid water at the soil temperature, reaches thermal equilibrium with leaves, and is lost to the canopy air space as water vapor at the leaf temperature. Therefore, the enthalpy flux between the soil layers and the cohort (${\dot{H}}_{g_j,l_k}$) is calculated in a similar manner to Eq. (35), whereas the enthalpy flux between the leaf intercellular space and the canopy air space (${\dot{H}}_{l_k,c}$) is solved in a similar manner to Eq. (75):

$$\begin{array}{}\text{(97)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{{\text{g}}_j,l_k}={\dot{W}}_{{\text{g}}_j,l_k}{q}_{\mathrm{\ell }}\left(T_{{\text{g}}_j}-T_{\mathrm{\ell }\mathrm{0}}\right),\text{(98)}& {\displaystyle }& {\displaystyle }{\dot{H}}_{l_k,c}={\dot{W}}_{l_k,c}{q}_{pv}\left(T_{t_k}-T_{v\mathrm{0}}\right).\end{array}$$

4.7 Non-leaf autotrophic respiration

Respiration from fine roots is defined using a phenomenological function of temperature that has the same functional form as leaf respiration (Moorcroft et al., 2001). Because roots are allowed in multiple layers, and in ED-2.2 roots have a uniform distribution of mass throughout the profile, the total respiration (${\dot{C}}_{r_k,c}$: ${\mathrm{kg}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is the integral of the contribution from each soil layer, weighted by the layer thickness:

$$\begin{array}{}\text{(99)}& {\dot{C}}_{r_k,c}=C_{r_k}{\displaystyle \frac{{\sum }_{j=j_{\mathrm{0}k}}^{N_{\mathrm{G}}}\left[{\mathcal{T}}^{\prime }\left(T_{{\text{g}}_j},r_{r_k}\right)\mathrm{\Delta }{z}_{{\text{g}}_j}\right]}{{\sum }_{j=j_{\mathrm{0}k}}^{N_{\mathrm{G}}}\mathrm{\Delta }{z}_{{\text{g}}_j}}},\end{array}$$

where $r_{r_k}$ (s⁻¹) is the PFT-dependent factor that describes the relative metabolic activity of fine roots at the reference temperature (15 ^∘C) (Tables S5–S6), and 𝒯^′ is the same temperature-dependent function from Eq. (86); default parameters are listed in Tables S5–S6.

Plants have two additional respiration terms: a phenomenological term that represents the long-term turnover rate of the accumulated storage pool (${\dot{C}}_{n_k,c}$), assumed constant (Medvigy et al., 2009), and a term related to the losses associated with the assimilated carbon for growth and maintenance of the living tissues (${\dot{C}}_{{\mathrm{\Delta }}_k,c}$;  Amthor, 1984). The latter term is a function of the plant metabolic rate, which has strong daily variability and hence is a function of the daily carbon balance:

$$\begin{array}{}\text{(100)}& {\displaystyle }& {\displaystyle }{\dot{C}}_{n_k,c}={\mathit{\tau }}_{n_k}{C}_{n_k},\text{(101)}& {\displaystyle }& {\displaystyle }{\dot{C}}_{{\mathrm{\Delta }}_k,c}={\mathit{\tau }}_{{\mathrm{\Delta }}_k}{C}_{{\mathrm{\Delta }}_k},\end{array}$$

where $({\mathit{\tau }}_{n_k},{\mathit{\tau }}_{{\mathrm{\Delta }}_k})$ are the PFT-dependent decay rates associated with storage turnover and consumption for growth, respectively (Tables S5–S6); and $C_{{\mathrm{\Delta }}_k}$ (kg_C m⁻²) is the total accumulated carbon from the previous day as defined in Eq. (25). The transport from non-structural storage and the accumulated carbon for maintenance, growth, and storage is summarized in Sect. S3.

4.8 Heterotrophic respiration

Heterotrophic respiration comes from the decomposition of carbon in the three soil/litter carbon pools. For each carbon pool $e_j;j\in \left(\mathrm{1},\mathrm{2},\mathrm{3}\right)$, we determine the maximum carbon loss based on the characteristic decay rate, which corresponds to the typical half-life for metabolic and microbial litter (fast, e₁), structural litter (intermediate, e₂), and humified and passive soil carbon (slow, e₃), determined from Bolker et al. (1998):

$$\begin{array}{}\text{(102)}& {\dot{C}}_{e_j,c}=C_{e_j}{f}_{h{e}_j}{B}_{e_j}{\mathcal{E}}_T\left({\stackrel{\mathrm{‾}}{T}}_{g_{\mathrm{20}}}\right){\mathcal{E}}_{{\mathit{\vartheta }}^{\prime }}\left({\stackrel{\mathrm{‾}}{\mathit{\vartheta }}}_{\mathrm{20}}^{\prime }\right),\end{array}$$

where f_he is the fraction of decay that is lost through respiration (Table S4), and by definition $f_{h{e}_{\mathrm{3}}}$ must always be 1 (slow soil carbon can only be lost through heterotrophic respiration); $B_{e_j}$ is the decay rates at optimal conditions of soil carbon e_j, based on Bolker et al. (1998) (Table S4); ${\stackrel{\mathrm{‾}}{T}}_{g_{\mathrm{20}}}$ and ${\stackrel{\mathrm{‾}}{\mathit{\vartheta }}}_{\mathrm{20}}^{\prime }$ are the average temperature and relative soil moisture of the top 0.2 m of soil; the relative soil moisture for each layer is defined as

$$\begin{array}{}\text{(103)}& \mathit{\vartheta }{{}^{\prime }}_{{\text{g}}_j}={\displaystyle \frac{{\mathit{\vartheta }}_{{\text{g}}_j}-{\mathit{\vartheta }}_{\text{Re}}}{{\mathit{\vartheta }}_{\text{Po}}-{\mathit{\vartheta }}_{\text{Re}}}};\end{array}$$

and ${\mathcal{E}}_T\left({\stackrel{\mathrm{‾}}{T}}_{g_{\mathrm{20}}}\right)$ and ${\mathcal{E}}_{{\mathit{\vartheta }}^{\prime }}\left({\stackrel{\mathrm{‾}}{\mathit{\vartheta }}}_{\mathrm{20}}^{\prime }\right)$ are functions that reduces the decomposition rate due to temperature or soil moisture under non-optimal conditions:

$$\begin{array}{}\text{(104)}& \begin{array}{rl}& {\mathcal{E}}_T\left({\stackrel{\mathrm{‾}}{T}}_{g_{\mathrm{20}}}\right)=\\ & {\displaystyle \frac{\mathrm{1}}{\left\{\mathrm{1}+\exp \left[-{\widehat{f}}_{\text{Cold}}\left({\stackrel{\mathrm{‾}}{T}}_{g_{\mathrm{20}}}-T_{g_{\text{Cold}}}\right)\right]\right\}\left\{\mathrm{1}+\exp \left[+{\widehat{f}}_{\text{Hot}}\left({\stackrel{\mathrm{‾}}{T}}_{g_{\mathrm{20}}}-T_{g_{\text{Hot}}}\right)\right]\right\}}},\end{array}\end{array}$$

$$\begin{array}{}\text{(105)}& \begin{array}{rl}& {\mathcal{E}}_{{\mathit{\vartheta }}^{\prime }}\left({\stackrel{\mathrm{‾}}{\mathit{\vartheta }}}_{\mathrm{20}}^{\prime }\right)=\\ & {\displaystyle \frac{\mathrm{1}}{\left\{\mathrm{1}+\exp \left[-{\widehat{f}}_{\text{Dry}}\left({\stackrel{\mathrm{‾}}{\mathit{\vartheta }}}_{\mathrm{20}}^{\prime }-{\mathit{\vartheta }}_{\text{Dry}}^{\prime }\right)\right]\right\}\left\{\mathrm{1}+\exp \left[+{\widehat{f}}_{\text{Wet}}\left({\stackrel{\mathrm{‾}}{\mathit{\vartheta }}}_{\mathrm{20}}^{\prime }-{\mathit{\vartheta }}_{\text{Wet}}^{\prime }\right)\right]\right\}}},\end{array}\end{array}$$

where (${\widehat{f}}_{\text{Cold}}$;$T_{g_{\text{Cold}}}$), (${\widehat{f}}_{\text{Hot}}$;$T_{g_{\text{Hot}}}$), (${\widehat{f}}_{\text{Dry}}$;${\mathit{\vartheta }}_{\text{Dry}}^{\prime }$), and (${\widehat{f}}_{\text{Wet}}$;${\mathit{\vartheta }}_{\text{Wet}}^{\prime }$) are phenomenological parameters to decrease decomposition rates at low and high temperatures, and dry and saturated soils, respectively (Table S4). The decay fraction from fast and structural soil carbon that is not lost through heterotrophic respiration is transported to the slow soil carbon (Sect. S4).

5 Results

5.1 Conservation of energy, water, and carbon dioxide

The ED-2.2 simulations show a high degree of conservation of the total energy, water, and carbon (Fig. 6). In the example simulation for one patch at Paracou, French Guiana (GYF), a tropical forest site, the accumulated deviation from perfect closure (residual) of the energy budget over 50 years (2 629 800 time steps) was 0.1 % of the total enthalpy storage – sum of enthalpy stored at the canopy air space, cohorts, temporary surface water and soil layers, (Fig. 6a) and 0.002 % of the accumulated losses through eddy flux, the largest cumulative flux of enthalpy. Results for the water budget were even better, with maximum accumulated residuals of 0.04 % of the total water stored in the ED-2.2 thermodynamic systems, or 0.0006 % of the total water input by precipitation (Fig. 6b), and the accumulated residual of carbon was 0.008 % of the total carbon storage or 0.017 % the total accumulated loss through eddy flux. The average absolute residual errors by time step, relative to the total storage, ranged from $\mathrm{3.6}\times {\mathrm{10}}^{-\mathrm{11}}$ (carbon) to $\mathrm{3.8}\times {\mathrm{10}}^{-\mathrm{10}}$ (energy) and were thus orders of magnitude less than the truncation error of single-precision numbers ($\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{7}}$) and the model tolerance for each time step ($\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{5}}$).

Figure 6 Examples of (a) enthalpy, (b) water, and (c) carbon conservation assessment in ED-2.2, for a single-patch simulation at GYF for 50 years. Terms are presented as the cumulative contribution to the change storage. Total storage is the combination of canopy air space, cohorts, and temporary surface water and soil layers in the cases of enthalpy and water, and canopy air space, cohorts, seed bank, and soil carbon pools in the case of carbon. Positive (negative) values mean accumulation (loss) by the combined storage pool over the time. Pressure change accounts for changes in enthalpy when pressure from the meteorological forcing is updated, and density change accounts for changes in mass to ensure the ideal gas law. CAS change and vegetation heat capacity (Veg Hcap) change reflect the addition/subtraction of carbon, water, and enthalpy due to the vegetation dynamics modifying the canopy air space depth and the total heat capacity of the vegetation due to biomass accumulation or loss. Storage change is the net gain or loss of total storage, and the residual corresponds to the deviation from the perfect closure. Note that we present the y axis in cube root scale to improve visualization of the smallest terms.

The conservation of energy and water of ED-2.2 also represents a substantial improvement from previous versions of the model. We carried out additional decade-long simulations with ED-2.2 and two former versions of the model (ED-2.0.12 and ED-2.1) and the most similar configuration possible among versions, and found that cumulative residual of enthalpy relative to eddy flux loss decreased from 15.2 % (ED-2.0.12) or 5.7 % (ED-2.1) to 6.1 % × 10⁻⁵ % (ED-2.2) (Fig. S3a–c). Similarly, the cumulative violation of perfect water budget closure, relative to total precipitation input, decreased from 3.4 % (ED-2.0.12) or 1.1 % (ED-2.1) to 1.2 % × 10⁻⁴ % (ED-2.2) (Fig. S3d–f).

5.2 Simulated ecosystem heterogeneity

Because ED-2.2 accounts for the vertical distribution of the plant community and the local heterogeneity of ecosystems, it is possible to describe the structural variability of ecosystems using continuous metrics. To illustrate this, we show the results of a 600-year simulation (1400–2002) carried out for tropical South America, starting from near-bare-ground conditions and driven by the Princeton Global Meteorological Forcing (Sheffield et al., 2006, 1969–2008), and with active fires (Sect. S3.4). For the last 100 years, we also prescribed land use changes derived from Hurtt et al. (2006) and Soares-Filho et al. (2006). The distribution of basal area binned by diameter at breast height (DBH) classes shows high variability across the domain and even within biome boundaries (Fig. 7). For example, larger trees (DBH≥50 cm) are nearly absent outside the Amazon biome, with the exception of more humid regions such as the Atlantic forest along the Brazilian coast, western Colombia, and Panama (Fig. 7d,e). In contrast, in seasonally dry areas such as the Brazilian Cerrado, intermediate-sized trees ($\mathrm{10}\le \mathrm{DBH}<\mathrm{50}\mathrm{cm}$) contribute the most to the basal area (e.g., areas near the Brasília (BSB) site; Fig. 7b, c). Even within the Amazon ecoregion, basal area shows variability in the contribution of trees with different sizes, including the areas outside the arc of deforestation along the southern and eastern edges of the biome (Fig. 7). Similarly, the abundance of different plant functional groups shows great variability across the region, with dominance of grasses and early successional tropical trees in deforested regions and in drier areas in the Brazilian Cerrado, whereas late-successional tropical trees dominate the tropical forests, albeit with lower dominance in parts of central Amazonia (Fig. S4).

Figure 7 Simulated distribution of size-dependent basal area across tropical South America, aggregated for the following diameter at breast height (DBH) bins: (a) 0–10 cm; (b) 10–30 cm; (c) 30–50 cm; (d) 50–80 cm; (e) ≥80 cm. Maps were obtained from the final state of a 6000-year simulation (1400–2002), initialized with near-bare-ground conditions, active fires, and prescribed land use changes between 1900 and 2002. Points indicate the location of the example sites (Fig. 8): (□) Paracou (GYF), a tropical forest site; (◦) Brasília (BSB), a woody savanna site. White contour is the domain of the Amazon biome, and grey contours are the political borders.

The variability of forest structural and functional composition observed in regional simulations emerges from both the competition among cohorts in the local micro-environment and the environmental controls on the disturbance regime. In Fig. 8, we present the impact of different disturbance regimes modulating the predicted ecosystem structure and composition for two sites: Paracou (GYF), a tropical forest region in French Guiana, and Brasília (BSB), a woody savanna site in central Brazil. Both sites were simulated for 500 years using a 40-year meteorological forcing developed from local meteorological observations, following the methodology described in Longo et al. (2018); we allowed fires to occur but for simplicity we did not prescribe land use change. After 500 years of simulation, the structure at the two sites is completely different, with large, late-successional trees dominating the canopy at GYF (Fig. 8a) and open areas with shorter, mostly early successional trees dominating the landscape at BSB (Fig. 8b). For GYF, the structural and functional composition is achieved only after 200 years of simulation, whereas in BSB a dynamic steady state caused by the strong fire regime is achieved in about 100 years (Fig. S5). At both sites, early successional trees dominate the canopy at recently disturbed areas (Fig. 8c, d) with late-successional (GYF) or mid-successional trees (BSB) increasing in size only at the older patches (> 30 years, Fig. 8c, d), and the variation of basal area as a function of age since last disturbance show great similarity at both sites (Fig. 8e). However, the disturbance regimes are markedly different: at GYF, fires never occurred and disturbance was driven exclusively by tree fall (prescribed at 1.11 % year⁻¹), whereas fires substantially increase the disturbance rates at BSB (average fire return interval of 19.3 years). Consequently, old-growth patches (older than 100 years) are nonexistent at BSB and abundant at GYF (Fig. 8f). In addition, the high disturbance regime at BSB meant that large trees and late-successional trees (slow growers) failed to establish but succeeded and maintained a stable population at GYF (Fig. S5).

Figure 8 Examples of size, age, and functional structure simulated by ED-2.2, after 500 years of simulation using local meteorological forcing and active fires. (a, b) Individual realization of simulated stands for sites (a) Paracou (GYF, tropical forest); (b) Brasília (BSB, woody savanna). The number of individuals shown is proportional to the simulated stem density, the distribution in local communities is proportional to the patch area, the crown size and stem height are proportional to the cohort size, and the crown color indicates the functional group. (c, d) Distribution of cohorts as a function of size (DBH and height) and age since last disturbance (patch age) for sites (c) GYF and (d) BSB. Crown sizes are proportional to the logarithm of the stem density within each patch. (e, f) Patch-specific properties as a function of age since last disturbance (patch age) for sites GYF and BSB after 500 years of simulation: (e) basal area and (f) probability density function of patch age (fractional patch area). See Fig. 7 for the location of both example sites.

The impacts of simulating structurally and functionally diverse ecosystems are also observed in the fluxes of energy, water, carbon, and momentum. For example, in Fig. 9, we show the monthly average fluxes from the last 40 years of simulation at GYF, along with the interannual variability of the fluxes aggregated to the polygon level (hereafter polygon variability; error bars) and the interannual variability of the fluxes accounting for the patch probability (hereafter patch variability; colors in the background). The polygon-level variability can be thought as the variability attributable exclusively to climate variability, whereas the patch variability also incorporates the impact of the structural heterogeneity in the variability. Most highly aggregated (“big-leaf”) models characterize the polygon-level variability but not the patch variability. However, in all cases, the patch variability greatly exceeded the polygon variability, indicating that structural variability is as important as the interannual variability in complex ecosystems. In the case of sensible heat, the polygon variable was between 39 % and 64 % of the patch variability (Fig. 9a). The polygon-to-patch variability ratio was similar for both friction velocity (19 %–39 %) and water fluxes (17 %–44 %) (Fig. 9b, c). In the case of gross primary productivity, the relevance of patch variability was even higher, with the polygon-to-patch variability ratio ranging from 3.7 % during the dry season to 17 % during the wet season (Fig. 9d). Importantly, the broader range of fluxes across patches in the site can be entirely attributed to structural and functional diversity, because all patches were driven by the same meteorological forcing.

Figure 9 Monthly averages and variability of fluxes attributable to meteorological conditions and plant community heterogeneity combined with interannual variability. Results are shown for GYF, a tropical forest site, for (a) sensible heat flux; (b) friction velocity (momentum flux); (c) water vapor flux; and (d) gross primary productivity. The variability was calculated for the last 40 years of a 500-year simulation starting from near-bare ground. Points correspond to the 40-year monthly averages for the entire polygon, line bars correspond to the 2.5 %–97.5 % quantile of monthly averages aggregated at the polygon level (polygon interannual variability), and background colors represent the 40-year probability density function of monthly means for each simulated patch, and scaled by the area of each patch (patch interannual variability). Density function colors outside the 2.5 %–97.5 % quantile interval are not shown. Note that the density function scale is logarithmic. See Fig. 7 for the location of the example site.

6 Discussion

6.1 Conservation of biophysical and biogeochemical properties

As demonstrated in Sect. 5.1, it is possible to represent the long-term, large-scale dynamics of heterogeneous and functionally diverse plant canopy, while still accurately conserving the fluxes of carbon, water, and energy fluxes that occur in the ecosystem. ED-2.2 exhibits excellent conservation of energy, water, and carbon dioxide even in multi-decadal scales. After 50 years of simulation, the accumulated residuals from perfect closure never exceeded 0.1 % of the total energy, water, and carbon stored in the pools resolved by the model (Fig. 6), which is significantly less than the error accepted in each time step (1 %).

The model's excellent conservation of these three key properties is possible because the ordinary differential equations are written directly in terms of the variables that we sought to conserve, thus reducing the effects of non-linearities. A key feature that facilitates the model's high level of energy conservation is the use of enthalpy as the primary energy-related state variable within the model. This contrasts with most terrestrial biosphere models, which use temperature as their energy state variable (e.g., Best et al., 2011; Oleson et al., 2013). By using enthalpy, the model can seamlessly incorporate energy storage changes caused by rapid changes in water content and consequently heat capacity. It also reduces errors near phase changes (freezing or melting), when changes in energy may not correspond to changes in temperature. Nonetheless, the residual errors in ED-2.2 are larger than the error of each time step after integrating the model over multiple decades (Fig. 6), which suggests that the errors may have a systematic component that deserves further investigation. The main contribution to the remaining residual errors in carbon, water, and energy fluxes comes from the linearization of the prognostic equations due to changes in density in the canopy air space (Eqs. 18–19; 23). The magnitude of these residuals would likely be further reduced by using the bulk enthalpy, water content, and carbon dioxide content in the canopy air space as the state variables instead of the specific enthalpy, specific humidity, and CO₂ mixing ratio.

Unlike most existing terrestrial biosphere models (but see SiB2, e.g., Baker et al., 2003; Vidale and Stöckli, 2005), in ED-2.2, we explicitly include the dynamic storage of energy, water, and carbon dioxide in the canopy air space. Canopy air space storage is particularly important in tall, dense tropical forests; accounting for this storage term, as well as the energy storage of vegetation, allows a more realistic representation of the fluxes between the ecosystem and the air above (see also Haverd et al., 2007). In addition, the separation of the ecosystem fluxes in the model into eddy fluxes and change in canopy air space storage allows a thorough evaluation of the model's ability to represent both the total exchange and the ventilation of water, energy, and carbon in and out of the ecosystem with eddy covariance towers, as shown in the companion paper (Longo et al., 2019a).

6.2 Heterogeneity of ecosystems

It has been long advocated that terrestrial biosphere models must incorporate demographic processes and ecosystem heterogeneity to improve their predictive ability in a changing world (Moorcroft, 2006; Purves and Pacala, 2008; Evans, 2012; Fisher et al., 2018). In ED-2.2, we aggregate individuals and forest communities according to similar characteristics (Fig. 1). For example, individuals are only aggregated into cohorts if they are of similar size, same functional group, and live in comparable micro-environments. Likewise, local plant communities are aggregated only if their disturbance history and their vertical structure are similar. The level of aggregation of ED-2.2 still allows mechanistic representation of ecological processes such as how individuals' access to and competition for resources vary depending on their size, adaptation, and presence of other individuals. This approach allows representing a broad range of structure and composition of ecosystems (Figs. 7, S4), as opposed to simplified biome classification. In this paper, we presented the functional diversity using only the default tropical PFTs, which describe the functional diversity along a single functional trait axis of broadleaf tropical trees. However, the ED-2.2 framework allows users to easily modify the traits and trade-offs of existing PFTs, or include new functional groups; previous studies using ED-2.2 have leveraged this feature of the code to define PFTs according to the research question both in the tropics (e.g., Xu et al., 2016; Trugman et al., 2018; Feng et al., 2018) and in the extratropics (e.g., Raczka et al., 2018; Bogan et al., 2019).

Previous analysis by Levine et al. (2016) has shown that the dynamic, fine-scale heterogeneity and functional diversity of the plant canopy in ED-2.2 is essential for capturing macro-scale patterns in tropical forest properties. Specifically, Levine et al. (2016) found that ED-2.1 was able to characterize the smoother observed transition in tropical forest biomass across a dry-season length gradient in the Amazon, whereas a highly aggregated (big-leaf-like) version of ED-2.1 predicted abrupt shifts in biomass, which is commonly observed in many dynamic global vegetation models (e.g., Good et al., 2011). Results from two related studies have shown that the incorporation of subgrid-scale heterogeneity and diversity within ED-2 also improves its ability to correctly capture the responses of terrestrial ecosystems to environmental perturbation. First, in an assessment of the ability of four terrestrial biosphere models to capture the impact of rainfall changes on biomass in Amazon forests (Powell et al., 2013), ED-2.1 was the only model that captured the timing and average magnitude of aboveground biomass loss that was observed in two experimental drought treatments, while all three big-leaf model formulations predicted minimal impacts of the drought experiment. Second, a recent analysis by Longo et al. (2018) on the impact of recurrent droughts in the Amazon found that drought-induced carbon losses in ED-2.2 arose mostly from the death of canopy trees, a characteristic that is consistent with field and remote-sensing observations of drought impacts in the region (Phillips et al., 2010; Yang et al., 2018).

Importantly, since its inception, the ED model accounts for the disturbance-driven horizontal heterogeneity of ecosystems (Moorcroft et al., 2001). As demonstrated in Moorcroft et al. (2001), the continuous development of tree-fall gaps is fundamental to explaining the long-term trajectory of biomass accumulation in tropical forests; for example, by representing both recently disturbed and old-growth fragments of forests, it is possible to simulate micro-environments where either shade-intolerant plants thrive or slow-growing, shade-tolerant individuals dominate the canopy (Fig. 8a, c). Moreover, ED-2.2 can also represent dynamic and diverse disturbance regimes, which ultimately mediate the regional variation of ecosystem properties. For example, tropical forests and woody savannas may share similarities in local communities with similar age since disturbance (Fig. 8e); however, because fire disturbances frequently affect large areas in the savannas, fragments of old-growth vegetation are nearly absent in these regions (Fig. 8f), which creates an environment dominated mostly by smaller trees (Fig. S5c).

Furthermore, the heterogeneity of ecosystems in ED-2.2 is integrated across all timescales, because we solve the biophysical and biogeochemical cycles for each cohort and each patch separately (Figs. 2–3). While solving the cycles at subgrid scale adds complexity, it also improves the characterization of heterogeneity of available water and energy for plants of different sizes, even within the same polygon: for example, the light profile and soil water availability are not only determined by meteorological conditions but also by the number of individuals, their height and rooting depth, and their traits and trade-offs that determine their ability to extract soil moisture or assimilate carbon. As a result, the variability in ecosystem functioning represented by ED-2.2 is significantly increased relative to the variability that a highly aggregated model based on the average ecosystem structure would be able to capture (Fig. 9).

6.3 Current and future developments

In this paper, we focused on describing the biophysical and biogeochemical core of the ED-2.2 model, and appraising its ability to represent both short-term (intra-annual and interannual) and long-term (decades to century) processes. However, the ED-2.2 community is continuously developing and improving the model. In this section, we summarize some of the recent and ongoing model developments being built on top of the ED-2.2 dynamic core.

Terrestrial biosphere models still show significant uncertainties in representing photosynthesis due to missing processes and inconsistencies in parameter estimations (Rogers et al., 2017). We are currently implementing the carboxylation limitation by the maximum electron transport rate and by the triose phosphate utilization (von Caemmerer, 2000; Lombardozzi et al., 2018), and constrained by observations (Norby et al., 2017), and nitrogen and phosphorus limitation have been recently incorporated (Medvigy et al., 2019). In addition, the model has also been recently updated to mechanistically represent plant hydraulics, and initial results indicate a significant improvement of the model's predictions of water use efficiency and water stress in tropical forests in Central America (Xu et al., 2016). Also, to better represent the dynamics of soil carbon in ED-2.2, we are implementing and optimizing a more detailed version of the CENTURY decomposition model (Bolker et al., 1998).

To improve the representation of surface and soil water dynamics, the model has been coupled with a hydrological routine model that accounts for lateral flux of water as a function of terrain characteristics and simulates river discharge (Pereira et al., 2017; Arias et al., 2018). Moreover, an integrated approach of hydraulic routing based on TOPMODEL (Walko et al., 2000; Beven and Freer, 2001), which allows exchange of water and internal energy exchange between different sites as a function of topographic characteristics, is being implemented in ED-2.2.

The ED-2.2 model framework is designed to simulate functionally diverse ecosystems, but trait values within each functional group are fixed. To account for the observed plasticity in many leaf traits, a new parameterization of leaf trait variation as function of the light level, based on the parameterization by Lloyd et al. (2010) and Xu et al. (2017) is being implemented. In addition, the ED-2.2 model has also been recently updated to represent the light competition and parasite–host relationships between lianas and trees (di Porcia e Brugnera et al., 2019), and it is currently being extended to incorporate plant functional types from different biogeographic regions, such as temperate semi-arid shrublands (Pandit et al., 2018), as well as boreal ecosystems, building on previous works using ED-1 (Ise et al., 2008).

Anthropogenic forest degradation is pervasive throughout the tropics (Lewis et al., 2015). To improve the model's ability to represent damage and recovery from degradation, we are implementing a selective logging module that represents the direct impact of felling of marketable tree stems, and accounts the damage associated with skid trails, roads, and decks, which are modulated by logging intensity and logging techniques (Pereira Jr. et al., 2002; Feldpausch et al., 2005). In addition, the original fire model has been recently improved to account for size- and bark-thickness-dependent survivorship (Trugman et al., 2018), and is being developed to account for natural and anthropogenic drivers of ignition, fire intensity, fire spread, and fire duration, building on existing process-based fire models (Thonicke et al., 2010; Le Page et al., 2015).

The complexity and sophistication of ED-2.2 also creates important scientific challenges. For example, the multiple processes for functionally diverse ecosystems represented by the model also require a large number of parameters, with some of them being highly uncertain given the scarcity of data. To explore the effect of parameter uncertainty on model results and leverage the growing number of observations, the ED-2.2 model has been fully integrated with the Predictive Ecosystem Analyzer (LeBauer et al., 2013; Dietze et al., 2014), a hierarchical Bayesian-based framework that constrains model parameters based on available data and quantifies the uncertainties on model predictions due to parameter uncertainty.

Importantly, the need to incorporate terrestrial ecosystem heterogeneity in Earth system models has been long advocated (e.g., Moorcroft, 2006; Purves et al., 2008; Evans, 2012), but only recently have global models been incorporating ecological mechanisms that allow representing functionally diverse and heterogeneous biomes at global scale without relying on artificial climate envelopes. Two examples are the Geophysical Fluid Dynamics Laboratory (GFDL) Land Model version 3 with perfect plasticity approximation (Weng et al., 2015, LM3-PPA;) and the Functionally Assembled Terrestrial Ecosystem Simulation (FATES; Fisher et al., 2015), which incorporated the patch and cohort structure of ED-2.2 into the Community Land Model (CLM; Oleson et al., 2013) framework.

7 Conclusions

ED-2.2 represents a significant advance in how to integrate a variety of processes ranging across multiple timescales in heterogeneous landscapes: it retains all the detailed representation of the long-term dynamics of functionally diverse, spatially heterogeneous landscapes and long-term dynamics from the original ED ecosystem model (Moorcroft et al., 2001; Hurtt et al., 2002; Albani et al., 2006) but also solves for the associated energy, water, and CO₂ fluxes of plants living in horizontally and vertically stratified micro-environments within the plant canopy, which was initially implemented by Medvigy et al. (2009) (ED-2) by adapting the big-leaf land surface model LEAF-3 (Walko et al., 2000) to the cohort-based structure of ED-2.

The results presented in the model description demonstrated that ED-2.2 has a high degree of conservation of carbon, energy, and water, even over multi-decadal scales (Fig. 6). Importantly, the current formulation of the model allows representation of functional and structural diversity both at local and regional scales (Figs. 7–8; S4–S5), and the effect of the heterogeneity on energy, water, carbon, and momentum fluxes (Fig. 9). In the companion paper, we use data from eddy covariance towers, forest inventory, bottom-up estimates of carbon cycles, and remote-sensing products to assess the strengths and limitations of the current model implementation (Longo et al., 2019a).

This paper focused on the major updates to the energy, water, and carbon cycle within the ED-2.2 framework; the model continues to be actively developed. Some of the further developments include implementing more mechanisms that influence photosynthesis and water cycle, such as plant hydraulics; additional nutrient cycles; expanding the representation of plant functional diversity, including trait plasticity and lianas; and expanding the types of natural and anthropogenic disturbances. ED-2.2 is a collaborative, open-source model that is readily available from its repository, and the scientific community is encouraged to use the model and contribute with new model developments.