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

Abstract

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 . 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 , 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 . Models for the dynamics of the terrestrial biosphere and its bidirectional interaction with the atmosphere have evolved considerably over the past decades . As described by , 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.,. 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 ; examples of these models include National Center for Atmospheric Research Biosphere-Atmosphere Transfer Scheme (NCAR/BATS) and Simple Biosphere model (SiB) . 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 ; examples of such models included LSM and SiB2 . 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 , 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 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 , Community Land Model's Dynamic Global Vegetation Model (CLM-DGVM) , and Top-down Representation of Interactive Foliage and Flora Including Dynamics Joint UK Land Environment Simulator (TRIFFID/JULES) 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., and their carbon stocks and fluxes . However, they diverge markedly in their predictions of how terrestrial ecosystems will respond to future climate change . 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.,. 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 .

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

Field- and laboratory-based studies conducted over the past 30 years indicate that plant functional diversity significantly affects ecosystem functioning and references therein, and variations in trait expression are strongly driven by disturbances and local heterogeneity of abiotic factors such as soil characteristics . In many cases, biodiversity increases ecosystem productivity and ecosystem stability e.g.,, and biodiversity has also been shown to contribute to enhanced ecosystem functionality in highly stressed environments e.g.,. Other studies have also established correlations between tropical forest diversity and carbon storage and primary productivity .

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 . 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 , which affects carbon stocks and forest structure and composition, which cannot be easily represented in highly aggregated models .

An alternative approach to simulating the dynamics of terrestrial ecosystems has been individual-based vegetation models . 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 , as well as the impact of climate change in the ecosystem carbon balance 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 . 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.,; 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;; the Lund–Potsdam–Jena General Ecosystem Simulator LPJ-GUESS;; the Land Model version 3 with perfect plasticity approximation LM3-PPA;; and the Functionally-Assembled Terrestrial Ecosystem Simulator FATES;. 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 .

The Ecosystem Demography model ED; 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 and 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 .

The original ED model formulation was an offline ecosystem model describing the coupled carbon and water fluxes of a heterogeneous tropical forest ecosystem . Subsequently, 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., and simulations running interactively with a regional atmospheric model e.g.,.

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

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.

[Figure omitted. See PDF]

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 ($\sim \mathrm{10}\mathrm{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 ) 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 $\mathit{\alpha }$, which can be also thought of as the relative area within a site, that satisfies 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}.$$ 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 ) and the type of gap ($q$). The number density distribution depends on the individuals' biomass characteristics (size, $\mathit{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: ${\mathrm{kg}}_{\mathrm{C}}{\mathrm{plant}}^{-\mathrm{1}})$ 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_{\mathrm{T}}$ is the total number of cohorts, $N_{\mathrm{G}}$ is the total number of soil (ground) layers, $N_{\mathrm{S}}$ is the total number of temporary surface water/snowpack layers, and $N_{\mathrm{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.

Subscript	Description
$X_{\mathrm{c}}$	Canopy air space (single layer)
$X_{{\text{c}}_j}$	Canopy air space, layer $j$ ($j\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{C}}\right\}$)
$X_{e_j}$	Necromass pools: $e_{\mathrm{1}}$, metabolic litter (fast); $e_{\mathrm{2}}$, structural debris (intermediate); $e_{\mathrm{3}}$, humified/dissolved (slow)
$X_{\text{f}}$	Plant functional type
$X_{{\text{g}}_j}$	Soil (ground), layer $j$ ($j\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{G}}\right\}$)
$X_q$	Disturbance type
$X_{{\text{s}}_j}$	Temporary surface water/snowpack, layer $j$ ($j\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{S}}\right\}$)
$X_{t_k}$	Cohort $k$ ($k\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{T}}\right\}$)
$X_u$	Patch $u$ ($u\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{P}}\right\}$)
$X_{y_k}$	Property $y$ of cohort $k$ ($k\in \left\{\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{T}}\right\}$). Possible values of $y$: branch wood ($b$), structural tissue (heartwood) ($h$), leaves ($l$), non-structural carbon storage (starch, sugars) ($n$), roots ($r$), total living tissues ($\mathit{\alpha }$), branch boundary layer ($\mathit{\beta }$), carbon balance ($\mathrm{\Delta }$), leaf boundary layer ($\mathit{\lambda }$), reproductive tissues ($\mathit{\varrho }$), sapwood ($\mathit{\sigma }$)

Following , , and , 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; ${\mathit{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; ${\mathbf{\nabla }}_{\mathit{C}}\cdot$ 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^{\prime }$ affected by new disturbance of type $q$, which may depend on environmental conditions. Boundary conditions are shown in Sect. S2.

Equations () and () 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 ) 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 ) of similar size and same PFT (Fig. ). 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., 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 ). 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.

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.

Time step	Timescale	Processes
Thermodynamics		Energy and water fluxes${}^{*}$
($\mathrm{\Delta }{t}_{\text{Thermo}}$)	$\mathrm{1}\mathrm{s}$–$\mathrm{\Delta }{t}_{\text{Bio}}$	Eddy fluxes (including ${\mathrm{CO}}_{\mathrm{2}}$ flux)${}^{*}$
		Most thermodynamic state functions${}^{*}$
		Meteorological and ${\mathrm{CO}}_{\mathrm{2}}$ forcing${}^{*}$
Biophysics		Radiation model${}^{*}$
($\mathrm{\Delta }{t}_{\text{Bio}}$)	2–$\mathrm{15}\min$	Photosynthesis model${}^{*}$
		Respiration fluxes (autotrophic and heterotrophic)${}^{*}$
		Evaluation of energy, water, and ${\mathrm{CO}}_{\mathrm{2}}$ budgets${}^{*}$
		Maintenance of active tissues
		Update of the storage pool
Phenology	$\mathrm{1}\mathrm{d}$	Leaf phenology
($\mathrm{\Delta }{t}_{\text{Phen}}$)		Plant carbon balance
		Integration of mortality rate due to cold
		Soil litter pools
		Growth of structural tissues
Cohort dynamics		Mortality rate
($\mathrm{\Delta }{t}_{\mathrm{CD}}$)	$\mathrm{1}\mathrm{month}$	Reproduction – cohort creation
		Integration of fire disturbance rate
		Cohort fusion, fission, and extinction
Patch dynamics	$\mathrm{1}\mathrm{year}$	Annual disturbance rates and patch creation
($\mathrm{\Delta }{t}_{\mathrm{PD}}$)		Patch fusion and termination

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 . 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 . 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;. 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 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., 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, ${\mathrm{CO}}_{\mathrm{2}}$ molar fraction, pressure of the air above the canopy, precipitation rate, incoming solar (shortwave) irradiance (radiation flux), and incoming thermal (longwave) irradiance (Table ) 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., or bias-corrected products based on reanalysis e.g.,. Whenever available, ${\mathrm{CO}}_{\mathrm{2}}$ must be provided at comparable temporal and spatial resolution to other meteorological drivers; otherwise, it is possible to provide spatially homogeneous, time-variant ${\mathrm{CO}}_{\mathrm{2}}$, or constant ${\mathrm{CO}}_{\mathrm{2}}$, although this may increase uncertainties in the model predictions e.g.,. Alternatively, the meteorological forcing (including ${\mathrm{CO}}_{\mathrm{2}}$) may be provided directly by BRAMS .

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.

Variable	Description	Units
$u_x$	Zonal wind speed	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
$u_y$	Meridional wind speed	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$
$p_{\mathrm{a}}$	Free-air pressure	$\mathrm{Pa}$
$T_{\mathrm{a}}$	Free-air temperature	$\mathrm{K}$
$w_{\mathrm{a}}$	Free-air specific humidity	${\mathrm{kg}}_{\mathrm{W}}{\mathrm{kg}}^{-\mathrm{1}}$
$c_{\mathrm{a}}$	Free-air ${\mathrm{CO}}_{\mathrm{2}}$ mixing ratio	${\mathrm{mol}}_{\mathrm{C}}{\mathrm{mol}}^{-\mathrm{1}}$
$z_{\mathrm{a}}$	Height of the reference point above canopy	$\mathrm{m}$
${\dot{W}}_{\mathrm{\infty },\mathrm{a}}$	Precipitation mass rate	${\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
${\dot{Q}}_{\mathrm{TIR}(\mathrm{\infty },\mathrm{a})}^{\Downarrow }$	Downward thermal infrared irradiance	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\dot{Q}}_{\mathrm{PAR}(\mathrm{\infty },\mathrm{a})}^{\odot }$	Downward photosynthetically active irradiance, direct	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\dot{Q}}_{\mathrm{PAR}(\mathrm{\infty },\mathrm{a})}^{\Downarrow }$	Downward photosynthetically active irradiance, diffuse	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\dot{Q}}_{\mathrm{NIR}(\mathrm{\infty },\mathrm{a})}^{\odot }$	Downward near-infrared irradiance, direct	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\dot{Q}}_{\mathrm{NIR}(\mathrm{\infty },\mathrm{a})}^{\Downarrow }$	Downward near-infrared irradiance, diffuse	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$

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. ), comprised of multiple thermodynamic systems: each soil layer (total number of layers $N_{\mathrm{G}}$), each temporary surface water or snow layer (total number of layers $N_{\mathrm{S}}$), leaves, and branch wood portion of each cohort (total number of cohorts $N_{\mathrm{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 . For a complete list of subscripts and variables used in this paper, refer to Tables S1–S2.

State variable	Description	Units	Reference equation
$c_{\mathrm{c}}$	${\mathrm{CO}}_{\mathrm{2}}$ mixing ratio – canopy air space	${\mathrm{mol}}_{\mathrm{C}}{\mathrm{mol}}^{-\mathrm{1}}$	()
$h_{\mathrm{c}}$	Specific enthalpy – canopy air space	$\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$	()
${\mathcal{H}}_{{\text{g}}_j}$	Volumetric enthalpy – soil layer $j$	$\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}$	()${}^{\text{a}}$
$H_{{\text{s}}_j}$	Enthalpy – temporary surface water layer $j$	$\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$	()
$H_{t_k}$	Enthalpy – cohort $k$	$\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$	()
$p_{\mathrm{c}}$	Atmospheric pressure – canopy air space	$\mathrm{Pa}$	(S81)${}^{\text{b}}$
$w_{\mathrm{c}}$	Specific humidity – canopy air space	${\mathrm{kg}}_{\mathrm{W}}{\mathrm{kg}}^{-\mathrm{1}}$	()
$W_{{\text{s}}_j}$	Water mass – temporary surface water layer $j$	${\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}^{-\mathrm{2}}$	()
$W_{t_k}$	Intercepted/dew/frost water mass – cohort $k$	${\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}^{-\mathrm{2}}$	()
$z_{\mathrm{c}}$	Depth (specific volume) – canopy air space	$\mathrm{m}$	()${}^{\text{c}}$
${\mathit{\vartheta }}_{{\text{g}}_j}$	Volumetric soil moisture – soil layer $j$	${\mathrm{m}}_{\mathrm{W}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$	()${}^{\text{a}}$

${}^{\text{a}}$ Budget fluxes are in units of area, and the state variable is updated following the conversion described in Sect. . ${}^{\text{b}}$ Canopy air space pressure is not solved using ordinary differential equations but based on the atmospheric pressure from the meteorological forcing. ${}^{\text{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$ ($\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$), and the mass continuity for incompressible fluids for total water mass $W$ (${\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}^{-\mathrm{2}}$): 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}$$ 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}}},$$ where $\mathcal{V}$ is the volume of the thermodynamic system and $p$ is the ambient pressure. The components on the right-hand side of Eqs. () and () 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. ). 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 , 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 .

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 $\mathrm{0}\mathrm{K}$; for water, enthalpy is defined as zero when water is at $\mathrm{0}\mathrm{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. ). 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. , and these cycles are solved every thermodynamic substep ($\mathrm{\Delta }{t}_{\text{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.  indicate the directions allowed in ED-2.2. The list of fluxes solved in ED-2.2 is provided in Table , 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_{\mathrm{T}}$ cohorts, $N_{\mathrm{G}}$ soil layers, and $N_{\mathrm{S}}=\mathrm{1}$ temporary surface water. Both $N_{\mathrm{G}}$ and the maximum $N_{\mathrm{S}}$ are specified by the user; $N_{\mathrm{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.

[Figure omitted. See PDF]

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 .

Variable	Description	Section	Equation	Units
${\dot{C}}_{a,c}$	${\mathrm{CO}}_{\mathrm{2}}$ flux from turbulent mixing		()	$\mathrm{kgC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
${\dot{C}}_{e_j,c}$	Heterotrophic respiration flux (soil carbon pool $e_j$)		()
${\dot{C}}_{l_k,c}$	Net leaf (cohort $k$)–CAS ${\mathrm{CO}}_{\mathrm{2}}$ flux${}^{\text{a}}$		()
${\dot{C}}_{n_k,c}$	Storage turnover (cohort $k$) respiration flux		()
${\dot{C}}_{r_k,c}$	Fine-root (cohort $k$) metabolic respiration flux		()
${\dot{C}}_{{\mathrm{\Delta }}_k,c}$	Growth and maintenance (cohort $k$) respiration flux		()
${\dot{Q}}_{a,g_{\text{NG}}}$	Net absorbed irradiance (topmost soil layer)		()	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\dot{Q}}_{a,{\text{s}}_j}$	Net absorbed irradiance (TSW layer $j$)		()
${\dot{Q}}_{a,t_k}$	Net absorbed irradiance (cohort $k$)		()
${\dot{Q}}_{g_{\text{NG}},c}$	Ground–CAS net sensible heat flux		()
${\dot{Q}}_{g_{j-\mathrm{1}},{\text{g}}_j}$	Net sensible heat flux between two soil layers		()
${\dot{Q}}_{s_{\text{NS}},c}$	TSW–CAS net sensible heat flux		()
${\dot{Q}}_{s_{j-\mathrm{1}},{\text{s}}_j}$	Net sensible heat flux between two TSW layers		()
${\dot{Q}}_{t_k,c}$	Cohort $k$–CAS net sensible heat flux		()
${\dot{H}}_{a,c}$	Enthalpy flux from turbulent mixing at the top of CAS		()	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$
${\dot{H}}_{a,s_{\text{NS}}}$	Enthalpy flux to the top TSW layer associated with throughfall precipitation		()
${\dot{H}}_{a,t_k}$	Enthalpy flux associated with rainfall interception by cohort $k$		()
${\dot{H}}_{g_{\text{NG},c}}$	Enthalpy flux associated with ground–CAS evaporation${}^{\text{b}}$		()
${\dot{H}}_{g_{j-\mathrm{1}},{\text{g}}_j}$	Enthalpy flux associated with water percolation between two soil layers		()
${\dot{H}}_{{\text{g}}_j,l_k}$	Enthalpy flux associated with soil water extraction from soil layer $j$ by cohort $k$		()
${\dot{H}}_{g_{\mathrm{1}},g_{\mathrm{0}}}$	Enthalpy flux associated with subsurface runoff from the bottom soil layer		()
${\dot{H}}_{l_k,c}$	Enthalpy flux associated with transpiration by cohort $k$		()
${\dot{H}}_{s_{\text{NS},c}}$	Enthalpy flux associated with TSW–CAS evaporation${}^{\text{b}}$		()
${\dot{H}}_{s_{\text{NS},o}}$	Enthalpy flux associated with surface runoff from the top TSW layer		()
${\dot{H}}_{s_{j-\mathrm{1}},{\text{s}}_j}$	Enthalpy flux associated with water percolation between two TSW layers		()
${\dot{H}}_{t_k,c}$	Enthalpy flux associated with evaporation${}^{\text{b}}$ of intercepted water (cohort $k$)		()
${\dot{H}}_{t_k,s_{\text{NS}}}$	Enthalpy flux associated with canopy dripping from cohort $k$ to the top TSW layer		()
${\dot{W}}_{a,c}$	Water flux from turbulent mixing at the top of CAS		()	${\mathrm{kg}}_{\mathrm{W}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$
${\dot{W}}_{a,s_{\text{NS}}}$	Precipitation throughfall flux to the top TSW layer		()
${\dot{W}}_{a,t_k}$	Water flux from rainfall interception (cohort $k$)		()
${\dot{W}}_{g_{\text{NG},c}}$	Ground–CAS evaporation${}^{\text{b}}$ flux		()
${\dot{W}}_{g_{j-\mathrm{1}},{\text{g}}_j}$	Water percolation flux between two soil layers		()
${\dot{W}}_{{\text{g}}_j,l_k}$	Water flux associated with soil water extraction by plants		()
${\dot{W}}_{g_{\mathrm{1}},g_{\mathrm{0}}}$	Water flux associated with subsurface runoff from the bottom soil layer		()
${\dot{W}}_{l_k,c}$	Transpiration flux (cohort $k$)		()
${\dot{W}}_{s_{j-\mathrm{1}},{\text{s}}_j}$	Water percolation flux between two TSW layers		()–()
${\dot{W}}_{s_{\text{NS},o}}$	Surface runoff water flux from the top TSW layer		()
${\dot{W}}_{s_{\text{NS},c}}$	TSW–CAS evaporation${}^{\text{b}}$ flux		()
${\dot{W}}_{t_k,c}$	Evaporation${}^{\text{b}}$ flux from intercepted water (cohort $k$)		()
${\dot{W}}_{t_k,s_{\text{NS}}}$	Canopy dripping flux from cohort $k$ to the top TSW layer		()

${}^{\text{a}}$ Net flux between leaf respiration (positive) and gross primary productivity (negative). ${}^{\text{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=\mathrm{1}$ being the deepest soil layer, and $j=N_{\mathrm{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}$ ($\mathrm{J}{\mathrm{m}}^{-\mathrm{3}}$) 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. ()–() 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 ${\mathit{\rho }}_{\mathrm{\ell }}$ is the density of liquid water (Table S3). Soil net fluxes for any layer $j$ are defined as

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}$$ 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}$$ 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}$$ where ${\mathit{\delta }}_{{\text{g}}_j{g}_{j^{\prime }}}$ is the Kronecker delta for comparing two soil layers ${\text{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. ). 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 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_{\mathrm{S}}$ and the thickness of each layer depends on the total mass and the water phase, following . When the layer is in liquid phase, only one layer ($N_{\mathrm{S}}=\mathrm{1}$) is maintained. If a snowpack develops, the temporary surface water can be divided into several layers (subscript $j$, with $j=\mathrm{1}$ being the bottommost TSW layer, and $j=N_{\mathrm{S}}$ being the topmost TSW layer). Net TSW fluxes are defined as

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}$$ 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}$$ 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}$$ where ${\mathit{\delta }}_{{\text{s}}_j{s}_{j^{\prime }}}$ is the Kronecker delta for comparing two TSW layers ${\text{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. ), 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. ()–() 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_{\mathrm{S}}+\mathrm{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 ${\mathit{\rho }}_{\mathrm{\ell }}$ 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 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

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}}},$$ 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}$$ 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}$$ Each term is described in detail in the sections shown underneath each term on the right-hand side of Eqs. ()–().

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_{\mathrm{c}}$ and total water mass $W_{\mathrm{c}}$ of the canopy air space can be written in terms of air density ${\mathit{\rho }}_{\mathrm{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}$ (${\mathrm{cm}}^{\mathrm{2}}$) and $z_{t_k}$ ($\mathrm{m}$) are the basal area and the height of cohort $k$, respectively; and $N_{T\left(\text{canopy}\right)}$ 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\left(\text{canopy}\right)}$ is the total number of cohorts, and a minimum value of $\mathrm{5}\mathrm{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 ($\mathrm{\Delta }{t}_{\mathrm{CD}}$, Table ). If we substitute Eqs. () and () into Eqs. () and (), 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

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}$$ 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}$$ 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}$$

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_{\mathrm{c}}$. In ED-2.2, CAS pressure is not solved through a differential equation; instead, $p_{\mathrm{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. (). Likewise, CAS density (${\mathit{\rho }}_{\mathrm{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. . The canopy air space is the only thermodynamic system with ${\mathrm{CO}}_{\mathrm{2}}$ storage that is solved by ED-2.2; nonetheless, we assume that the contribution of ${\mathrm{CO}}_{\mathrm{2}}$ to density and heat capacity of the canopy air space is negligible; hence, only the molar ${\mathrm{CO}}_{\mathrm{2}}$ mixing ratio $c_{\mathrm{c}}$ (${\mathrm{mol}}_{\mathrm{C}}{\mathrm{mol}}^{-\mathrm{1}}$) is traced.

Figure 3

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

[Figure omitted. See PDF]

The change in ${\mathrm{CO}}_{\mathrm{2}}$ storage in the canopy air space is determined by the following differential equation: 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}},$$ 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}$$ where ${\mathcal{M}}_{\mathrm{d}}$ and ${\mathcal{M}}_{\mathrm{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. () 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 that combines the decomposition rates from three soil carbon pools, defined by their characteristic lifetime: fast (metabolic litter and microbial; $e_{\mathrm{1}}$), intermediate (structural debris; $e_{\mathrm{2}}$), and slow (humified and passive soil carbon; $e_{\mathrm{3}}$). Note that the soil carbon pools are not directly related to the soil layers used to describe the thermodynamic state (Sect. ).

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: 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}$$ 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 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-\mathrm{1}$ and $j$ is determined from thermal conductivity ${\mathrm{{\rm Y}}}_{\mathrm{Q}}$ and temperature gradient , with an additional term for temporary surface water to scale the flux when the temporary surface water covers only a fraction $f_{\text{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-\mathrm{1}$ and $j$ to the height at the interface. The bottom boundary condition of Eq. () 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. () 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-\mathrm{1}}$ and ${\text{g}}_j,j\in \left\{\mathrm{2},\mathrm{3},\mathrm{\dots },N_{\mathrm{G}}\right\}$ is determined from Darcy's law :

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},$$ where $\mathrm{\Psi }$ is the soil matric potential and ${\mathrm{{\rm Y}}}_{\mathrm{\Psi }}$ is the hydraulic conductivity, both defined after , 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. () 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. (). Let $\mathit{ð}$ 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 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{ð}.$$ Special cases of Eq. () are the zero-flow conditions ($\mathit{ð}=\mathrm{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 . Liquid water in excess of $\mathrm{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: 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}$$ 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}$$

Surface runoff of liquid water is simulated using a simple extinction function, applied only at the topmost temporary surface water layer: 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),$$ where $t_{\text{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_{\mathrm{\ell }}$ is the specific heat of liquid water (Table S3), and $T_{\mathrm{\ell }\mathrm{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: 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.,$$ where the subscript $x_j$ represents either soil (${\text{g}}_j$) or temporary surface water (${\text{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 ($\mathcal{O}$) 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_{\mathrm{a}}$), and we use $T_{\mathrm{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_{\mathrm{a}}$ is above the water triple point ($T_{\mathrm{3}}=\mathrm{273.16}\mathrm{K}$); in this case, the rain temperature is always assumed to be at $T_{\mathrm{a}}$. Pure snow occurs when the free-air temperature is below $T_{\mathrm{3}}$, and likewise snow temperature is assumed to be $T_{\mathrm{a}}$. When free-air temperature is only slightly above $T_{\mathrm{3}}$, a mix of rain and snow occurs, with the rain temperature assumed to be $T_{\mathrm{a}}$ and snow temperature assumed to be $T_{\mathrm{3}}$:

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],$$ where $(q_i;q_{\mathrm{\ell }})$ are the specific heats of ice and liquid, respectively, and $T_{\mathrm{\ell }\mathrm{0}}$ is temperature at which supercooled water would have enthalpy equal to zero (Eq. S53). The fraction of precipitation that falls as rain ${\mathrm{\ell }}_{\mathrm{a}}$ is based on the parameterization, slightly modified to make the function continuous: 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..$$

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. () and (): $$\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}}$ .

4.3 Radiation model

The radiation budget is solved using a multi-layer version of the two-stream model applied to three broad spectral bands: photosynthetically active radiation (PAR, wavelengths between 0.4 and 0.7 $\mathrm{\mu }$m), near-infrared radiation (NIR, wavelengths between 0.7 and 3.0 $\mathrm{\mu }$m) and thermal infrared radiation (TIR, wavelengths between $\mathrm{3.0}$ and 15 $\mathrm{\mu }$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 ). In the cases of PAR ($m=\mathrm{1}$) and NIR ($m=\mathrm{2}$), the downward irradiance is comprised of a beam (direct) and isotropic (diffuse) components, whereas TIR irradiance ($m=\mathrm{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 , 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

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 }}},$$ 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}$$ 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}$$ 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. ); interface $N_{\mathrm{T}}+\mathrm{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$; ${\mathit{\varsigma }}_{mk}$ is the scattering coefficient, and thus ($\mathrm{1}-{\mathit{\varsigma }}_{mk}$) is the absorptivity; ${\mathit{\beta }}_{mk}^{\odot }$ and ${\mathit{\beta }}_{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_{\mathrm{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 ); $Z^{\odot }$ 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. –); ${\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_k$ is the effective plant area index (Eq. S100); $\left({\mathit{\varsigma }}_{mk};{\mathit{\beta }}_{mk}\right)$ 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. ); and ${\dot{Q}}_{a,t_k}^{\odot }$ is the net absorbed irradiance (Eq. ).

[Figure omitted. See PDF]

Equations ()–() simplify for each spectral band. First, ${\dot{Q}}_{mk}^{\odot }\equiv \mathrm{0}$ for the TIR ($m=\mathrm{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=\mathrm{1}$) and NIR ($m=\mathrm{2}$) bands, because thermal emission is negligible at these wavelengths. The black-body emission for the TIR band is defined as 48$${\dot{Q}}_{m=\mathrm{3}k}^{\mathit{⧫}}={\mathit{\sigma }}_{\mathrm{SB}}{T}_{t_k}^{\mathrm{4}},$$ where ${\mathit{\sigma }}_{\mathrm{SB}}$ is the Stefan–Boltzmann constant (Table S3). Note that for emission of TIR radiation (Eqs. –), we assume that emissivity is the same as absorptivity Kirchhoff's law; and hence the $\left(\mathrm{1}-\mathit{\varsigma }\right)$ 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 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}$: 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}$$ This term is then used in the enthalpy budget of each cohort (Eqs.  and ).

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 (${\mathit{\varsigma }}_{m\mathrm{0}}$; 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:

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}$$ where ($\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}g}$) and ($\mathrm{1}-{\mathit{\varsigma }}_{\mathrm{3}s}$), are, respectively, the thermal-infrared emissivities of the top soil layer and the temporary surface water (Table S4), and $f_{\text{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 .

Once the irradiance profile for the canopy is determined from Eqs. () to (), 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: 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}$$ 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: 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}.$$

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 , which has been widely used by biosphere–atmosphere models representing a variety of biomes e.g.,, although this is often an extrapolation of the theory that was not originally developed for heterogeneous vegetation or tall vegetation .

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 ${\mathrm{CO}}_{\mathrm{2}}$. This assumption allows us to define a single canopy conductance $G_{\mathrm{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_{\mathrm{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_{\mathrm{0}}$ is the reference pressure, $\mathcal{R}$ is the universal gas constant, $q_{\mathrm{pd}}$ is the specific heat of dry air at constant pressure, and ${\mathcal{M}}_{\mathrm{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 ), although it is not directly applied to the energy balance in the canopy air space (${\dot{H}}_{\mathrm{a},\mathrm{c}}$ is used instead):

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}$$ 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}}.$$

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. , ). Throughout this section, we use subscripts ${\mathit{\lambda }}_k$ and ${\mathit{\beta }}_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. ). Let $G_{Q{\mathit{\lambda }}_k}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and $G_{W{\mathit{\lambda }}_k}$ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) 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 $\mathit{\pi }$ in Eq. () 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. (). 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_{\text{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 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_{\text{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 :

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.$$ 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\},$$ where $g$ is the gravity acceleration, ${\mathcal{M}}_{\mathrm{w}}$ is the water molar mass, and $\mathcal{R}$ 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 ${\mathit{\vartheta }}_{\mathrm{Fc}}$ and ${\mathit{\vartheta }}_{\mathrm{Re}}$ are the soil moisture at field capacity and the residual soil moisture, respectively. The exponential term in Eq. () corresponds to the soil pore relative humidity derived from the Kelvin equation , and $s_g$ is the soil wetness function, which takes a similar functional form as the relative humidity term from and the $\mathit{\beta }$ term from . 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_{lv\mathrm{3}}$ and $l_{iv\mathrm{3}}$ are the specific latent heats of vaporization and sublimation at the water triple point ($T_{\mathrm{3}}$), $q_{pv}$ is the specific heat of water vapor at constant pressure, and $q_i$ and $q_{\mathrm{\ell }}$ 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 ${\mathrm{\ell }}_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

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}$$ By using the definitions from Eq. (S54), Eq. () can be further simplified to 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)}},$$ which is consistent with the exchange of pure water vapor and enthalpy between the thermodynamic systems. Equation () 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 and for ${\mathrm{C}}_{\mathrm{3}}$ plants; for ${\mathrm{C}}_{\mathrm{4}}$ plants; and the model for stomatal conductance. This submodel ultimately determines the net leaf-level ${\mathrm{CO}}_{\mathrm{2}}$ 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^{-\mathrm{1}}$). Leaf-level sensible heat flux (${\dot{q}}_{{\mathit{\lambda }}_k,c}$; Eq. ) and leaf-level vapor flux between intercepted water and canopy air space (${\dot{w}}_{{\mathit{\lambda }}_k,c}$; Eq. ) are also shown for comparison. Cohort index $k$ is omitted from the figure for clarity.

[Figure omitted. See PDF]

The exchange of water and ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$) 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. ). When soil moisture is not a limiting factor, the fluxes of ${\mathrm{CO}}_{\mathrm{2}}$ and water can be written as 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}$$ 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}$$ where 78$$w_{l_k}=w_{\text{Sat}}\left(T_{t_k},p_{\mathrm{c}}\right)\text{(Sect. S15)},$$ and 79$${\widehat{G}}_{X{\mathit{\lambda }}_k}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{c}}{G}_{X{\mathit{\lambda }}_k}}{{\mathcal{M}}_d}},$$ and 80$${\widehat{G}}_{X{l}_k}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{c}}{G}_{X{l}_k}}{{\mathcal{M}}_d}},$$ where $G_{X{\mathit{\lambda }}_k}$ and $G_{X{l}_k}$ (units $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) 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 ${\mathrm{CO}}_{\mathrm{2}}$ mixing ratio and the specific humidity of the leaf boundary layer, respectively; $c_{l_k}$ and $w_{l_k}$ are the ${\mathrm{CO}}_{\mathrm{2}}$ and specific humidity of the leaf intercellular space, respectively; ${\mathcal{M}}_d$ is the molar mass of dry air; and ${\mathit{\rho }}_c$ is the air density. As stated in Eq. (), 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 ${\mathrm{CO}}_{\mathrm{2}}$ assimilation flux and stomatal conductances are described below.

From , the net ${\mathrm{CO}}_{\mathrm{2}}$ assimilation flux is defined as 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}}.$$ Oxygenation releases $\mathrm{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 ${\mathrm{CO}}_{\mathrm{2}}$ compensation point ${\mathrm{\Gamma }}_k$ : 82$${\dot{V}}_{O_k}={\displaystyle \frac{\mathrm{2}{\mathrm{\Gamma }}_k}{c_{l_k}}}{\dot{V}}_{C_k},$$ where $c_{l_k}$ is the ${\mathrm{CO}}_{\mathrm{2}}$ mixing ratio in the leaf intercellular space. The ${\mathrm{CO}}_{\mathrm{2}}$ compensation point is determined after : 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.,$$ where $o_{\oplus }$ is the reference ${\mathrm{O}}_{\mathrm{2}}$ mixing ratio (Table S3), and $\mathit{\varpi }$ represents the ratio between the rates of carboxylase to oxygenase and is a function of temperature. The general form of the function $\mathcal{T}$ describing the metabolic dependence upon temperature for any variable $x$ (including $\mathit{\varpi }$) is 84$$\mathcal{T}\left(T,x\right)=x_{\mathrm{15}}\times {\mathcal{Q}}_{{\mathrm{10}}_x}^{\frac{T-T_{\mathrm{15}}}{\mathrm{10}}},$$ where $x_{\mathrm{15}}$ is the value of variable $x$ at temperature $T_{\mathrm{15}}=\mathrm{288.15}$ $\mathrm{K}$ (15 ${}^{\circ }$C), and ${\mathcal{Q}}_{{\mathrm{10}}_x}$ is the parameter which describes temperature dependence (Table S4).

Because ${\mathrm{C}}_{\mathrm{4}}$ plants have a mechanism to concentrate ${\mathrm{CO}}_{\mathrm{2}}$ near the ${\mathrm{CO}}_{\mathrm{2}}$-fixing enzyme RuBisCO (ribulose-1,5-biphosphate carboxylase oxygenase), photorespiration is nearly nonexistent in ${\mathrm{C}}_{\mathrm{4}}$ plants and hence the assumption that ${\mathrm{\Gamma }}_k$ is zero. For ${\mathrm{C}}_{\mathrm{4}}$ 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 ${\mathrm{C}}_{\mathrm{3}}$, this rate is unattainable even under RuBP-saturated conditions because carboxylation and oxygenation are mutually inhibitive reactions . Therefore, the maximum attainable carboxylation (${\dot{V}}_{C_k}={\dot{V}}_{C_k}^{\text{RuBP}}$) is expressed by a modified Michaelis–Menten kinetics equation: 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.,$$ 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. () (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 : 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}$$ where $f_{\text{Cold}}$, $f_{\text{Hot}}$, $T_{\text{Cold}}$, and $T_{\text{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 ${\mathrm{CO}}_{\mathrm{2}}$ (${\dot{V}}_{C_k}^{\mathrm{InSl}}$) for ${\mathrm{C}}_{\mathrm{4}}$ plants by has been modified later e.g., to explicitly include ${\dot{V}}_{C_k}^{\text{max}}$; this is the same expression used in ED-2.2: 87$${\dot{V}}_{C_k}^{\mathrm{InSl}}=k_{\mathrm{PEP}}{\dot{V}}_{C_k}^{\text{max}}{c}_{l_k},$$ where $k_{\mathrm{PEP}}$ represents the initial slope of the response curve to increasing ${\mathrm{CO}}_{\mathrm{2}}$; the default value in ED-2.2 (Table S4) is the same value used by .

From the total photosynthetically active irradiance absorbed by the cohort ${\dot{Q}}_{\mathrm{PAR}:a,t_k}$ (Eq. ), 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}}$): 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},$$ where $\mathrm{Ein}$ is the average photon-specific energy in the PAR band (0.4–0.7 $\mathrm{\mu }\mathrm{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 , 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 (${\mathit{ϵ}}_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 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..$$ Carboxylation may also be limited by the export rate of starch and sucrose that is synthesized by triose phosphate, especially when ${\mathrm{CO}}_{\mathrm{2}}$ concentration and irradiance are simultaneously high, at low temperatures, or ${\mathrm{O}}_{\mathrm{2}}$ concentration is low . 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 : 90$${\dot{R}}_k=f_R{\dot{V}}_{C_k}^{max},$$ 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. proposed a semi-empirical stomatal conductance expression for water based on these trade-offs: 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.,$$ 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 $\mathrm{\Delta }{w}_k$ is an empirical coefficient controlling conductance under severe leaf-level water deficit; ${\widehat{G}}_{W{l}_k}^{\mathit{\varnothing }}$, $M_k$, and $\mathrm{\Delta }{w}_k$ are PFT-dependent parameters (Tables S5–S6). From , stomatal conductance of ${\mathrm{CO}}_{\mathrm{2}}$ is estimated by the ratio $f_{Gl}$ between the diffusivities of water and ${\mathrm{CO}}_{\mathrm{2}}$ in the air (Table S4): 92$${\widehat{G}}_{W{l}_k}=f_{Gl}{\widehat{G}}_{C{l}_k}.$$

Variables $w_{l_k}$, ${\dot{V}}_{C_k}^{\text{max}}$, ${\dot{R}}_k$, ${\mathit{\varpi }}_k$, ${\mathcal{K}}_{O_k}$, ${\mathcal{K}}_{C_k}$, ${\mathrm{\Gamma }}_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${}_{\mathrm{2}}$-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 (Eq. ) is regulated by leaf vapor pressure deficit. However, Eqs. () and () 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 ${\mathrm{CO}}_{\mathrm{2}}$ (${\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 ${\mathit{þ}}_k$ is 1 if the PFT is hypostomatous or 2 if the PFT is amphistomatous or needleleaf (Tables S5–S6). Alternatively, 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. (), 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. (): $$\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 . 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:

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}}},$$ where $r_{r_k}$ (${\mathrm{s}}^{-\mathrm{1}}$) is the PFT-dependent factor that describes the relative metabolic activity of fine roots at the reference temperature (15 ${}^{\circ }$C) (Tables S5–S6), and ${\mathcal{T}}^{\prime }$ is the same temperature-dependent function from Eq. (); 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 , 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}$; . 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}$ (${\mathrm{kg}}_{\mathrm{C}}{\mathrm{m}}^{-\mathrm{2}}$) 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_{\mathrm{1}}$), structural litter (intermediate, $e_{\mathrm{2}}$), and humified and passive soil carbon (slow, $e_{\mathrm{3}}$), determined from :

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),$$ 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 (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 $\mathrm{0.2}\mathrm{m}$ of soil; the relative soil moisture for each layer is defined as 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}}}};$$ 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: 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}$$ 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}$$ 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. ). 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 $\mathrm{0.1}$ % of the total enthalpy storage – sum of enthalpy stored at the canopy air space, cohorts, temporary surface water and soil layers, (Fig. a) and $\mathrm{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 $\mathrm{0.04}$ % of the total water stored in the ED-2.2 thermodynamic systems, or $\mathrm{0.0006}$ % of the total water input by precipitation (Fig. b), and the accumulated residual of carbon was $\mathrm{0.008}$ % of the total carbon storage or $\mathrm{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.

[Figure omitted. See PDF]

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 $\mathrm{15.2}$ % (ED-2.0.12) or $\mathrm{5.7}$ % (ED-2.1) to 6.1 % $\times$ ${\mathrm{10}}^{-\mathrm{5}}$ % (ED-2.2) (Fig. S3a–c). Similarly, the cumulative violation of perfect water budget closure, relative to total precipitation input, decreased from $\mathrm{3.4}$ % (ED-2.0.12) or $\mathrm{1.1}$ % (ED-2.1) to 1.2 % $\times$ ${\mathrm{10}}^{-\mathrm{4}}$ % (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 1969–2008, and with active fires (Sect. S3.4). For the last 100 years, we also prescribed land use changes derived from and . The distribution of basal area binned by diameter at breast height ($\mathrm{DBH}$) classes shows high variability across the domain and even within biome boundaries (Fig. ). For example, larger trees ($\mathrm{DBH}\ge \mathrm{50}\mathrm{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. d,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. b, 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. ). 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 $\mathrm{cm}$; (b) 10–30 $\mathrm{cm}$; (c) 30–$\mathrm{50}\mathrm{cm}$; (d) 50–80 $\mathrm{cm}$; (e) $\ge \mathrm{80}\mathrm{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. ): ($\mathit{\square }$) 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.

[Figure omitted. See PDF]

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. , 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 ; 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. a) and open areas with shorter, mostly early successional trees dominating the landscape at BSB (Fig. b). 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. c, d) with late-successional (GYF) or mid-successional trees (BSB) increasing in size only at the older patches ($>$ 30 years, Fig. c, d), and the variation of basal area as a function of age since last disturbance show great similarity at both sites (Fig. e). However, the disturbance regimes are markedly different: at GYF, fires never occurred and disturbance was driven exclusively by tree fall (prescribed at $\mathrm{1.11}\mathrm{\%}{\mathrm{year}}^{-\mathrm{1}}$), 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. f). 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.  for the location of both example sites.

[Figure omitted. See PDF]

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. , 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. a). The polygon-to-patch variability ratio was similar for both friction velocity (19 %–39 %) and water fluxes (17 %–44 %) (Fig. b, 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 $\mathrm{3.7}$ % during the dry season to $\mathrm{17}$ % during the wet season (Fig. d). 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.  for the location of the example site.

[Figure omitted. See PDF]

6 Discussion

6.1 Conservation of biophysical and biogeochemical properties

As demonstrated in Sect. , 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 $\mathrm{0.1}$ % of the total energy, water, and carbon stored in the pools resolved by the model (Fig. ), which is significantly less than the error accepted in each time step ($\mathrm{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.,. 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. ), 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. –; ). 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 ${\mathrm{CO}}_{\mathrm{2}}$ mixing ratio.

Unlike most existing terrestrial biosphere models but see SiB2, e.g.,, 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. 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 .

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 . In ED-2.2, we aggregate individuals and forest communities according to similar characteristics (Fig. ). 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. , 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., and in the extratropics e.g.,.

Previous analysis by 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, 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.,. 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 , 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 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 .

Importantly, since its inception, the ED model accounts for the disturbance-driven horizontal heterogeneity of ecosystems . As demonstrated in , 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. a, 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. e); however, because fire disturbances frequently affect large areas in the savannas, fragments of old-growth vegetation are nearly absent in these regions (Fig. f), 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. –). 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. ).

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 . We are currently implementing the carboxylation limitation by the maximum electron transport rate and by the triose phosphate utilization , and constrained by observations , and nitrogen and phosphorus limitation have been recently incorporated . 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 . 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 .

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 . Moreover, an integrated approach of hydraulic routing based on TOPMODEL , 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 and 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 , and it is currently being extended to incorporate plant functional types from different biogeographic regions, such as temperate semi-arid shrublands , as well as boreal ecosystems, building on previous works using ED-1 .

Anthropogenic forest degradation is pervasive throughout the tropics . 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 . In addition, the original fire model has been recently improved to account for size- and bark-thickness-dependent survivorship , 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 .

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 , 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.,, 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 LM3-PPA; and the Functionally Assembled Terrestrial Ecosystem Simulation FATES;, which incorporated the patch and cohort structure of ED-2.2 into the Community Land Model CLM; 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 but also solves for the associated energy, water, and ${\mathrm{CO}}_{\mathrm{2}}$ fluxes of plants living in horizontally and vertically stratified micro-environments within the plant canopy, which was initially implemented by (ED-2) by adapting the big-leaf land surface model LEAF-3 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. ). Importantly, the current formulation of the model allows representation of functional and structural diversity both at local and regional scales (Figs. –; S4–S5), and the effect of the heterogeneity on energy, water, carbon, and momentum fluxes (Fig. ). 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 .

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.

Code availability

The ED-2.2 software and further developments are publicly available. The most up-to-date source code, post-processing R scripts, and an open discussion forum are available on the GitHub repository (The ED-2 Model Development Team, 2014). The code described in this paper, along with a wiki-based technical manual, is stored as a permanent release at https://github.com/mpaiao/ED2/releases/tag/rev-86 (last access: 25 September 2019) and permanently stored at .

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-12-4309-2019-supplement.

Author contributions

ML, RGK, DMM, MCD, YK, RLB, SCW, and PRM designed the ED-2.2 model. ML, RGK, DMM, NML, MCD, YK, ALSS, KZ, CR, and PRM developed the model. ML, RGK, NML, and ALSS carried out the ED-2.2 simulations. ML, RGK, DMM, NML, MCD, YK, ALSS, and PRM wrote the paper.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

The research was partially carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We thank the reviewers Ian Baker and Stefan Olin, as well as Miriam Johnston, Luciana Alves, John Kim, and Shawn Serbin for suggestions that improved the manuscript, and Alexander Antonarakis, Fabio Berzaghi, Carl Davidson, Istem Fer, Miriam Johnston, Geraldine Klarenberg, Robert Kooper, Félicien Meunier, Manfredo di Porcia e Brugnera, Afshin Pourmokhtarian, Thomas Powell, Daniel Scott, Shawn Serbin, Alexey Shiklomanov, Anna Trugman, Toni Viskari, and Xiangtao Xu for contributing to the code development. The model simulations were carried out at the Odyssey cluster, supported by the FAS Division of Science, Research Computing Group at Harvard University. Marcos Longo was supported the NASA Postdoctoral Program, administered by Universities Space Research Association under contract with NASA. Abigail L. S. Swann was supported as a Giorgio Ruffolo Fellow in the Sustainability Science Program at Harvard University, for which support from Italy's Ministry for Environment, Land and Sea is gratefully acknowledged.

Financial support

This research has been supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (grant no. 200686/2005-4), the NASA Earth and Space Science Fellowship (grant no. NNX08AU95H), the National Science Foundation, Office of International Science and Engineering (grant no. OISE-0730305), the National Science Foundation (grant no. ATM-0449793), the National Aeronautics and Space Administration (grant no. NNG06GD63G), and the Fundação de Amparo à Pesquisa do Estado de São Paulo (grant no. 2015/07227-6).

Review statement

This paper was edited by Christoph Müller and reviewed by Ian Baker and Stefan Olin.