the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Robust Ecosystem Demography (RED version 1.0): a parsimonious approach to modelling vegetation dynamics in Earth system models
Chris Huntingford
Andrew J. Wiltshire
Anna B. Harper
Chris D. Jones
A significant proportion of the uncertainty in climate projections arises from uncertainty in the representation of land carbon uptake. Dynamic global vegetation models (DGVMs) vary in their representations of regrowth and competition for resources, which results in differing responses to changes in atmospheric CO_{2} and climate. More advanced cohortbased patch models are now becoming established in the latest DGVMs. These models typically attempt to simulate the size distribution of trees as a function of both tree size (mass or trunk diameter) and age (time since disturbance). This approach can capture the overall impact of stochastic disturbance events on the forest structure and biomass – but at the cost of increasing the number of parameters and ambiguity when updating the probability density function (pdf) in two dimensions. Here we present the Robust Ecosystem Demography (RED), in which the pdf is collapsed onto the single dimension of tree mass. RED is designed to retain the ability of more complex cohort DGVMs to represent forest demography, while also being parameter sparse and analytically solvable for the steady state. The population of each plant functional type (PFT) is partitioned into mass classes with a fixed baseline mortality along with an assumed powerlaw scaling of growth rate with mass. The analytical equilibrium solutions of RED allow the model to be calibrated against observed forest cover using a single parameter – the ratio of mortality to growth for a tree of a reference mass (μ_{0}). We show that RED can thus be calibrated to the ESA LC_CCI (European Space Agency Land Cover Climate Change Initiative) coverage dataset for nine PFTs. Using net primary productivity and litter outputs from the UK Earth System Model (UKESM), we are able to diagnose the spatially varying disturbance rates consistent with this observed vegetation map. The analytical form for RED circumnavigates the need to spin up the numerical model, making it attractive for application in Earth system models (ESMs). This is especially so given that the model is also highly parameter sparse.
A key requirement of Earth system science is to estimate how much carbon the land surface will take up in the decades ahead (Ciais et al., 2014). This is an important component of the total carbon budget consistent with avoiding global warming thresholds, such as 2 ^{∘}C (Schleussner et al., 2016). Unfortunately, projections of future land carbon storage still span a wide range (Brovkin et al., 2013; Friedlingstein et al., 2014; Arora et al., 2019). Beyond the CO_{2} and nutrient fertilization effects and landuse change, significant uncertainty also arises from the representation of vegetation demographics such as recruitment, competition, and mortality (Brovkin et al., 2013; Ahlström et al., 2015). The representation of plant communities within Earth system models (ESMs) is achieved through the use of dynamic global vegetation models (DGVMs). DGVMs employ a variety of biophysical, biogeographical, and biochemical processes to simulate growth, competition, and recruitment of vegetation. The variety in the number and resolution of the processes contributes to the differences found at the Earth system level.
Within the context of modelling vegetation at a global level, there is a tradeoff between the complexity of ecological process representation and the necessity of parsimony at scale (Fisher et al., 2018). DGVMs range from the simplistic, older, topdown approaches to that of complex individualbased DGVMs. For example, in the first instance the TRIFFID model (Cox, 2001) simulates the fractional area of each plant functional type (PFT) using phenomenological Lotka–Volterra equations. The benefit of the TRIFFID approach is its simplicity and robustness. However, the model suffers from the lack of size representation and other processes, which results in the overestimation of regrowth time (Burton et al., 2019). In the second instance, individualbased models can explicitly represent a multitude of biological and ecosystem processes at an individual plant level (Smith, 2001; Sato et al., 2007). The benefit of this is that sizedependent physiology and spatial heterogeneity can be explicitly represented. However, multiple ensemble members are often needed to construct meaningful forest statistics, which makes such models computationally expensive to run at large scales. Compromises between the complexity of individualbased and topdown DGVMs exist as a class of tree cohort models. In the ED model (Moorcroft et al., 2001; Medvigy et al., 2009) the tree population is partitioned between patch disturbance and biomass classes allowing for the scaling of process to be represented in both age and size. ED2 can realistically model forests around the world (boreal, rainforest, and temperate) (Medvigy et al., 2009; Fisher et al., 2018). However, parameterization of competition within cohort DGVMs can result in a wide spread of outcomes when simulating climate change (Fisher et al., 2010; Scheiter et al., 2013).
In a similar vein other models have limited the number of cohort dimensions. The POP model (Haverd et al., 2014) uses standage cohorts as the dimension for population dynamics, every time step applying crowding and resource limited mortality rates. Another example is the ORCHIDEEMICT (Yue et al., 2018), which disaggregates the populations of a PFT into patch cohort functional types, with transitions between cohorts diagnosed when the average basal diameter passes a threshold.
This paper presents a simplified cohort model – Robust Ecosystem Demography (RED) – which updates the number of trees in each mass class but does not separately track tree age or patch age. RED assumes that the treesize distribution of a forest is determined by how the rates of tree growth and mortality vary with tree size (Kohyama et al., 2003; Coomes et al., 2003; MullerLandau et al., 2006; Lima et al., 2016). We follow many other studies in assuming that treegrowth rates vary with the threequarter power of tree mass (m^{3∕4}), as suggested by metabolic scaling theory (West et al., 1997). Where tree mortality rate can also be assumed to be approximately independent of tree mass, the demographic equation yields equilibrium treesize distributions which follow a Weibull distribution. This is sometimes termed demographic equilibrium theory (DET) (see Appendix B). These simplifications significantly reduce the number of free parameters in RED but still enable it to fit forest inventory data in North America (Moore et al., 2018) and South America (Moore et al., 2020).
A full list of variables, parameters, and units are given in Table 1.
2.1 Theory
The underlying theoretical model for RED is a continuity equation, for each PFT and spatial location, which describes the time evolution of the number density n of plants per unit area per unit mass m:
Here g is the growth rate, and γ is the mortality rate of a plant of mass m. In general, g and γ could be any reasonable function of tree size. For largescale applications we make simplifying assumptions for these functions consistent with observed n from forest inventory data (Moore et al., 2018, 2020). By default we assume that γ is independent of plant mass and that g follows a power law of plant mass as follows:
Here g_{0} is the growth rate of a plant with the reference mass, m_{0}. A value of ϕ_{g}=0.75 is assumed by default, consistent with the analysis of fieldbased measurements by Niklas and Spatz (2004). We also follow Niklas and Spatz (2004) in assuming the scaling of plant canopy area a with plant mass as follows:
where ϕ_{a}=0.5 by default. Solutions for n can be integrated over mass to derive the total plant number, $N={\int}_{\mathrm{0}}^{\mathrm{\infty}}n\text{d}m$, the total growth rate, $G={\int}_{\mathrm{0}}^{\mathrm{\infty}}gn\text{d}m$, the total biomass, $M={\int}_{\mathrm{0}}^{\mathrm{\infty}}mn\text{d}m$, and the fractional area covered $\mathit{\nu}={\int}_{\mathrm{0}}^{\mathrm{\infty}}an\text{d}m$.
2.2 Discrete mass classes
We wish to produce a model of vegetation demography that can be updated numerically and which explicitly conserves vegetation carbon, providing a constraint on the number of plants moving between mass classes in the discrete form. In order to do this we integrate Eq. (1) over finite mass ranges:
where i denotes the ith mass class; F_{i} is the flux of plants growing out of the ith mass class and into the (i+1)th mass class; F_{i−1} is the flux of plants growing out of the (i−1)th mass class and into the ith mass class; and N_{i} is the number of plants per unit area in the ith mass class. For clarity, Eq. (4) is deliberately presented as continuous in time at this stage, as the focus in this subsection is on discretization of the mass profile. The fully numerical version of RED, which includes discretization of time, is described in Sects. 2.4 and 2.5. In order to explicitly conserve carbon, the flux F_{i} must take the form (see Appendix A)
where m_{i} is the mean mass of a plant in the ith mass class, and g_{i} is the growth rate per plant of the ith mass class (kg C yr^{−1} plant^{−1}).
2.3 Seedling production and gap competition
To solve Eq. (4) we also require a lower boundary condition, which represents the rate at which seedlings of mass m_{0} are introduced into the cohort. Here we assume that a fixed fraction, α, of the total assimilate available to a PFT (P), is devoted to producing new seedlings, with the remainder $G=(\mathrm{1}\mathit{\alpha})P$ being allocated to the growth of existing plants. Spreading is homogeneous across the entirety of the grid box, but only seedlings established within “unoccupied” space will survive to join the plant cohort. The net incoming flux of seedlings of mass m_{0} is therefore
where s is the fractional gap area available for seedlings. The definition of s is assumed to differ by PFT to reflect an underlying tree–shrub–grass dominance hierarchy, as shown schematically in Fig. 1. Therefore, the rate of recruitment F_{0} is the ratio of a fraction of the carbon assimilate allocated to reproduction, αP, and m_{0}, multiplied by the gap area s.
The space available to the seedlings of the kth PFT is calculated from the area fractions of the PFTs to which it is subdominant as follows:
where ν_{l} is the area fraction of the lth PFT, and c_{kl} is the competition coefficient for the impact of PFT l on PFT k. If PFT l is within the same plant functional group (trees, shrubs or grasses) as PFT k or dominant over it, c_{kl}=1. If PFT k is dominant over PFT l, c_{kl}=0 (Fig. 1). This “gap” boundary condition results in there being no equilibrium solution where the amount of coverage exceeds 1. Doing so would halt the recruitment flux such that mortality processes would bring the fractional coverage back below unity. This is a similar competition regime to the Lotkainspired TRIFFID model (Cox, 2001) and allows for the coexistence between interfunctional groups (trees, shrubs, and grasses) of PFTs. For instance, a PFT such as broadleaf deciduous tree can coexist with a deciduous shrub and C3 grass. The hierarchy also enables the simulation of succession during regrowth. Fastergrowing species of grasses will not be able to expand into space occupied by trees and shrubs, unless there is space created by disturbance. A summary of the competition coefficients is given in Table 2.
2.4 Coupling to Earth system models
RED updates plant size distributions, biomass, and fractional areal coverage for an arbitrary number of PFTs at each spatial location and can be driven by variables provided by a land carbon cycle model, an Earth system model, or observations (see Fig. 2). For each PFT, the minimum required input is a time series of net carbon assimilate (P), defined as the difference between net primary productivity (Π_{N}) and local litter production due to turnover of leaves, stems, and roots (Λ_{l}):
We apply the m^{3∕4} scaling to P. We therefore implicitly assume the same scaling for both gross primary productivity and plant respiration. This is consistent with observations suggesting that plant production also scales approximately as m^{3∕4} (Enquist et al., 1998; Niklas and Enquist, 2001). Where available, additional mortality due to disturbance events such as droughts, fires, and anthropogenic deforestation (γ_{d}) can be added to the baseline mortality rates (γ_{b}), for each PFT as follows:
Disturbance rates γ_{d} can in principle be both PFT dependent and mass dependent (e.g. to capture forestry practices).
The input values of net assimilate for each PFT (P) define the total structural growth rate, $G=(\mathrm{1}\mathit{\alpha})P$, and the seedling flux F_{0} (via Eq. 6), using PFTspecific values of the parameter α (see Table 3). The definition of the total structural growth rate at a given time step is
which can be combined with the growth scaling given by Eq. (2) to derive the reference growth rate, g_{0}, from the net assimilate, P, which is a driving input:
This in turn enables the growth rate of each mass class to be calculated using Eq. (2). For each PFT, the number of plants in a mass class (N_{i}) is updated using a discretized form of Eq. (4):
where Δt is the RED time step (typically 1 month), and the superscript ^{(j)} denotes the jth time step. Our results are robust to changes in model time step so long as the time step remains small compared to the characteristic timescales associated with regrowth (${m}_{\mathrm{0}}/{g}_{\mathrm{0}}\sim $ 4 years) and plant mortality ($\mathrm{1}/\mathit{\gamma}\sim $ 20 years). The lowerboundary seedling flux is calculated from Eq. (6) using Eq. (7). We impose a zeroflux condition out of the upper mass class, under the assumption that there will be enough mass classes to ensure that this flux is negligible. However, to ensure carbon conservation on the land we add any plants that grow out of the upper mass class into a demographic litterfall term for each PFT, which is a RED output. This demographic litterfall term, Λ_{d}, keeps track of the carbon lost from the vegetation due to competition, mortality, and the carbon in any such plants that grow out of the largest resolved mass class (class I) as follows:
The first term on the righthand side of this equation represents carbon loss due to the shading of seedlings; the second term represents mortality of the resolved mass classes (which may include disturbance events); and the third term, which is normally very small, is the loss of vegetation carbon due to plants growing beyond the modelled mass classes. In order to initiate regrowth from bare soil, RED also assumes a minimum effective fractional area of each PFT. Where the net assimilate would be sufficiently negative to take the vegetation fraction below this minimum, the minimum value is maintained by subtraction from the demographic litter. The demographic litterfall term therefore represents the net addition litter production consistent with the prescribed net assimilate flux, the disturbance rate, and the change in vegetation carbon modelled by RED. When coupling to an ESM or land carbon model, the demographic litterfall term (Λ_{d}) should be added to the input local litterfall (Λ_{l}) (as used in Eq. 8) to calculate the total litterfall flux into the soil and/or litter system.
2.5 Steady state
The steady state of the continuum model defined by Eqs. (1) and (2) can be solved analytically for each PFT (Moore et al., 2018, 2020). The continuum analytical solutions for the equilibrium mass distribution n_{eq}(m), the total plant number (N_{eq}), biomass (M_{eq}), growth rate (G_{eq}), and fractional area (ν_{eq}) are summarized in Appendix B. The shape of the mass distribution and each of these parameters depend on the ratio of plant mortality to growth, which we choose to define for the reference mass class m_{0} as follows:
In order to initialize the numerical RED model in a driftfree initial state, we also derive the steady state of the discrete model (of Eq. 12), which will differ slightly from the continuum model for a finite number of mass classes. The equilibrium solution of Eq. (12) is derived in Appendix B2, based on the balance between seedling recruitment and total cohort mortality that defines the equilibrium state. The discretized version of RED thus yields formulae for the coverage (Eq. B28) and biomass densities (Eq. B30) which depend on the lowest mass class through the value of μ_{0}. Similarly, analytical expressions can be derived for total plant number and total growth rate of each PFT at equilibrium.

The total equilibrium stand density, N_{eq}:
$$\begin{array}{}\text{(15)}& {N}_{\mathrm{eq}}={N}_{\mathrm{0}}{X}_{N}.\end{array}$$ 
The total equilibrium structural growth, G_{eq}:
$$\begin{array}{}\text{(16)}& {G}_{\mathrm{eq}}=\sum _{i=\mathrm{0}}^{I}{N}_{i}{g}_{i}={N}_{\mathrm{0}}{g}_{\mathrm{0}}{X}_{G}.\end{array}$$ 
The total equilibrium coverage, ν_{eq}:
$$\begin{array}{}\text{(17)}& {\mathit{\nu}}_{\mathrm{eq}}=\sum _{i=\mathrm{0}}^{I}{N}_{i}{a}_{i}={N}_{\mathrm{0}}{a}_{\mathrm{0}}{X}_{\mathit{\nu}}.\end{array}$$ 
The total equilibrium carbon mass, M_{eq}:
$$\begin{array}{}\text{(18)}& {M}_{\mathrm{eq}}=\sum _{i=\mathrm{0}}^{I}{N}_{i}{m}_{i}={N}_{\mathrm{0}}{m}_{\mathrm{0}}{X}_{M}.\end{array}$$
Here X_{N}, X_{G}, X_{ν}, and X_{M} are functions of μ_{0} (see Appendix B2). This equilibrium state is derived by setting ${N}_{i}^{(j+\mathrm{1})}={N}_{i}^{\left(j\right)}$ in Eq. (B17), such that the flux entering into a mass class is equal to the flux leaving that class, due to growth out of the class and the loss of plants due to mortality.
The equations above therefore define the equilibrium state of the discrete system for given values of N_{0} and μ_{0}. The value of μ_{0} can be estimated from forest demographic data where this is available (Moore et al., 2018, 2020). However, for global applications we rarely have more observations than the fractional coverage of each PFT. Starting from the derived forms for N_{eq} (Eq. 15) and G_{eq} (Eq. 16) and requiring that the recruitment flux ($\mathit{\alpha}/(\mathrm{1}\mathit{\alpha}){G}_{\mathrm{eq}}s$) is equal to that of the total population dying (γN_{eq}), we can derive an equation for the total equilibrium coverage (full details in Appendix B2):
As the lefthand side of this equation depends only on prescribed constants and μ_{0}, Eq. (19) can be inverted (by numerical iteration) to estimate μ_{0} for observed values of the PFT fractions (ν_{k}, ν_{l}) and an assumed value of α (see Table 3). Once the value of μ_{0} has been derived in this manner, it can be used to calculate X_{ν} and, therefore, N_{0} by inversion of Eq. (B28):
Equations (19) and (20) therefore allow us to define an initial equilibrium state (N_{i}) which is consistent with observed area fractions of each PFT. Furthermore, when paired with an estimate of the net carbon assimilate (from a model or observations), the μ_{0} estimate can be converted into a map of the implied mortality (γ) by PFT. We demonstrate this capability globally in the next section.
For these runs, the numerical RED model is set up to use the nine PFTs which are currently used in JULES (Harper et al., 2018). This enables us to directly use driving data – time series of the rate of net assimilation (P) – from a previous UKESM model simulation that includes JULES (Sellar et al., 2019). RED is integrated forward using a 1month time step and successive mass classes that differ by a multiplicative constant ξ, so that ${m}_{i}=\mathit{\xi}{m}_{i\mathrm{1}}$. The value of ξ was chosen to optimally fit the analytical equilibrium solutions assuming 10 mass classes for trees, 8 mass classes for shrubs, and 1 mass class for grasses, assuming μ_{0}=0.25 (see Appendix B3). Other PFTspecific parameters are assumed as summarized in Table 3.
3.1 Global: diagnosed plant mortality rates
Here we use the analytical forms for the equilibrium state (Sect. 2.5) and observations of global vegetation cover to diagnose the corresponding map of PFTspecific mortality rates. These mortality rates are therefore consistent with the current observed vegetation state and rates of net assimilation (P) provided from UKESM (Sellar et al., 2019). The UKESM simulation provides net primary productivity (NPP) and local litterfall per unit area of each PFT. We multiply by PFT fraction to get the gridbox mean values required to drive RED (using ESA landcover data, as explained below). The observed maps of PFTs are provided by the ESA LC_CCI dataset for 2008–2012 (Poulter et al., 2015), projected onto the nine JULES PFTs (Fig. 3). Maps of the prescribed annual mean values of the rate of net assimilation (P) are shown in Fig. 4.
We use the procedure outlined in Sect. 2.5 to estimate spatially varying values of μ_{0} for each PFT, using Eq. (B32) and then Eq. (B34) to estimate N_{0}. This method successfully reproduces the ESA map of dominant PFT to good accuracy, as shown in Fig. 5 and Table 4.
The fit of the RED equilibrium vegetation coverage to the ESA observations is generally very good (Table 4). However, it is imperfect in some areas (e.g. Central Asia, Sahel) where the driving net assimilate from UKESM is zero or negative. Also, areas where the observational dataset indicates coexisting PFTs within the same vegetation class (e.g. broadleaf trees and needleleaf trees) are not well simulated by this first version of RED, which leads to competitive exclusion in the equilibrium state (see Discussion). Since we now have diagnosed values of μ_{0} and N_{0}, along with prescribed values of P, we can also diagnose the mean plant mortality rate γ, for each location and for each PFT, from Eq. (14) as follows:
where g_{0} is given by Eq. (11) combined with Eqs. (B18) and (B20). Maps of γ values, derived in this way, are shown in Fig. 6.
The mortality rate derived is dependent on the assumed areal coverage and the total assimilate. A high coverage with a low growth rate will result in a compensating low diagnosed mortality rate (and vice versa). Furthermore, the choice of α (Eq. 11) and m_{0} also influence the diagnosed value of γ. An analysis of the sensitivity of the inferred value of γ to these factors is presented in Appendix C. Assuming ±20 % uncertainty on assimilate, α, and m_{0} and ±5 % on the coverage gives an uncertainty bound of ±35 % on γ. Under the assumption that high coverages are indicative of the baseline mortality for a given PFT, we take a subsample of the grid boxes that are within the top quartile of nonzero coverages (ν_{eq}>0.01) (Table 5). The median μ_{0} value diagnosed from the top quartile of BETTr of ${\mathrm{0.232}}_{\mathrm{0.007}}^{+\mathrm{0.008}}$ (Table 5) is very close to the value calculated in our previous paper (Moore et al., 2020) of approximately 0.235 for all of South America using the RAINFOR sites.
Sitelevel assessments of the rates of stand mortality within pantropical forests conclude a range of background rates (Lugo and Scatena, 1996; Phillips, 1996; Phillips et al., 2004). Phillips (1996) estimates mortality rates collected across 40 pantropical sites for tree sizes greater than 10–25 cm dbh. Later work by Phillips et al. (2004) used the demographic data from the RAINFOR dataset of trees ≥ 10 cm dbh. Using these site assessments, we can make a comparison to BETTr equilibrium mortality rates by looking at the values of γ in areas where we would expect to see oldgrowth forests. We use the top 25 % of coverages of the BETTr PFT to represent plausible areas of undisturbed forest. Figure 7 shows that the diagnosed baseline mortality rates are in reasonable agreement with these observational estimates for Amazonia.
There is a need to better understand the influence of mortality arising from disturbance events such as droughts and fire in order to constrain model projections (Pugh et al., 2020). Here we investigate if the equilibrium mortality rates implicitly capture areas of disturbances, by comparing the mean tree mortality rate to fire and landuse surveys (the mean mortality is defined here by weighting gridbox γ values by gridbox fractional coverages). There are a number of surveys relating stand mortality in regions prone to wildfires (Swaine, 1992; Kinnaird and O’Brien, 1998; Peterson and Reich, 2001; Van Nieuwstadt and Sheil, 2005; Prior et al., 2009; Staver et al., 2009; Brando et al., 2014). In a broad sense, postfire mortality rates can range from 0.06 yr^{−1} to catastrophic rates around 0.8 yr^{−1} and can vary quite considerably depending on tree species, fire frequency, and drought severity. The drought–fire interaction is responsible for significantly increasing mortality post fire and can be a driving cause of regional dieback (Allen et al., 2010; Brando et al., 2014). Using the ESA FIRE_CCI dataset (Chuvieco et al., 2019) we can estimate the burnt vegetation fraction per year. Taking the average burnt vegetation fraction for the months between 2000 and 2010 and converting into annual burn rate, we gain an estimate of fire severity.
Another key issue is anthropogenic land use and landuse change (Nepstad et al., 2008; Haddad et al., 2015). Fragmentation of natural forests is understood to raise the mortality of the remaining forest and to decrease the overall resilience of the ecosystem (Esseen, 1994; Laurance et al., 1998; Jönsson et al., 2007). In order to maintain a nearconstant agricultural fraction, regular disruption such as grazing is needed to prevent recolonization and secondary succession (Dorrough and Moxham, 2005; Van Uytvanck et al., 2008; Chaturvedi et al., 2012). We carry out a comparison with land use using the 2000 ESA LC_CCI inferred crop coverages (Li et al., 2019).
In Fig. 9, we see the derived observations for burn area (panel a) and crop fraction (panel b), along with the derived mean γ for the tree PFTs (panel c). From Fig. 9d, we see that there are areas of large mortality ($\mathit{\gamma}>\mathrm{0.075}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{\mathrm{1}}$) that do correspond to areas where we see large fire activity (burn rate $>\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{\mathrm{1}}$) and increased crop fraction (> 0.25). However, large burn rates are seen to overlap in parts of central Brazil around the Cernado region, southern Africa, and north Western Australia where fires are understood to play a significant part within the ecosystem (Coutinho, 1990; Medeiros and Miranda, 2008; Prior et al., 2009; Staver et al., 2009). There are also some areas of agriculture which correspond to deforestation, such as in the Atlantic forests of Brazil and in Indonesia (Higuchi et al., 2008; Curran et al., 2004). Areas of increased disturbances result in grasses and shrubs dominating (Fig. 3).
Analysis of the RED equilibrium is an indirect approach to estimating tree mortality based on simple yet mechanistic principles of demography and relying on few inputs (vegetation cover and assimilate). It is however conditional on the assumed estimates of vegetation coverage and net rates of assimilation.
3.2 Dynamical simulations
3.2.1 Local: simulating succession
In this subsection we demonstrate the vegetation successional dynamics simulated by RED in an idealized spinup from bare soil, for a grid box at the edge of the Amazonian rainforest (Fig. 10). Under these circumstances, the diagnosed initial state is indeed the longterm equilibrium state, as evidenced by the horizontal dashed lines in panels (a) and (b) of Fig. 10.
Fastergrowing grass PFTs dominate the grid box within the first 12 years, before being replaced by evergreen shrubs which shade the grass seedlings. Eventually, broadleaf evergreen tropical trees replace much of the shrub and grass, on a timescale determined in large part by the parameter α and the reference mass class m_{0}. With the parameters used here, the vegetation fraction reaches close to its equilibrium value after about 20 years (panel a), but full spinup of the biomass takes around 150 years (panel b).
The modelled evolution of number density versus mass distribution for each PFT is shown in panel (c) (after 6 years), panel (d) (after 13 years), and panel (e) (after 100 years), with the eventual demographic equilibrium profiles shown by the dashed lines. It is clear that grass PFTs are close to their demographic equilibrium after only 6 years, but tree PFTs need more than 100 years to reach equilibrium.
The dashed lines in Fig. 10 represent a dynamical RED simulation from the diagnosed demographic equilibrium state. This state is derived using the methodology described in Sect. 2.5, with one significant change. The competition rules given by Eq. (7) and Table 2 result ultimately in equilibria which have a single dominate PFT in each class of cocompeting types (trees, shrubs, grasses). To avoid drifts associated with the competitive exclusion of the subdominant PFTs in each vegetation class, we choose to initialize the dominant PFT to have the total area fraction of all the PFTs in that vegetation class.
3.2.2 Global: spinup from bare soil
Transient simulations of global vegetation will be the subject of a future paper, but in the final subsection of this paper we wish to demonstrate the utility of the semianalytical equilibrium for initialization of global model runs. Figure 11 shows the time evolution of global mean PFT fractions and biomass from a global run driven by net assimilation rates from the UKESM model. Once again, two RED simulations are shown, one started from bare soil (solid lines) and the other from the semianalytical equilibrium state (dashed lines). Using a constant assimilate rate (Fig. 4) and the mortality distribution (Fig. 7), we see convergence of these two runs but only after more than 1000 years of simulated time. The ability to diagnose the equilibrium state therefore has the potential to reduce model spinup time hugely, especially for Earth system model (ESM) applications.
The response of the land surface to climate change is a key uncertainty in climate projections. Ambitious climate targets also rely on land management practices such as reforestation and afforestation to increase the storage of carbon on land. Firstgeneration dynamic global vegetation models (DGVMs) attempted to model the land surface in terms of bulk properties such as mean vegetation cover, vegetation carbon and leaf area index. These models lack information about the plant size distribution, which compromised their ability to represent recovery from disturbance and the impact of land management. Providing useful guidance on these issues requires improved DGVMs, which can represent changes in treesize distributions within forests (so called “demography”). A number of much more sophisticated secondgeneration DGVMs are now under development. These models often explicitly simulate the number of plants within different size or mass classes and on different patches of land, which are defined by the time since a disturbance event. Such secondgeneration models are therefore in principle able to simulate variations in plant number density as both a function of patch age and plant size. However, this completeness is at the expense of much computational and parameter complexity.
Our previous work in evaluating demographic equilibrium theory for regional forest inventory datasets in North America (Moore et al., 2018) and using RAINFOR sites for South America (Moore et al., 2020) has provided the theoretical basis for the development of RED. In those studies we found that the treesize distributions observed at a largescale in forests can be satisfactorily understood in terms of demographic equilibrium in the size dimension alone. This is a reduction in complexity compared to other cohort models which are based on patch age, and yet it is an improvement in ecological fidelity compared to older phenomenological DGVMs such as TRIFFID (Cox, 2001). The modular design of RED allows for easy coupling to landsurface schemes, merely requiring the per unit gridbox total carbon assimilate rate and any additional mortality disturbance rates as inputs for each grid box (Fig. 2). In principle, RED allows scope for more complex treesizedependent processes, although in this first study we chose to assume sizeindependent (but spatially varying) mortality rates for each PFT. Our previous work suggests that this is a good firstorder assumption (Moore et al., 2018, 2020).
Internally within the model we make a number of simplifications. Firstly, the number density for each PFT is treated as a function of plant mass alone. This immediately eliminates the need to explicitly represent patches and, therefore, removes age as an independent dimension. This is a distinct approach relative to cohort DGVMs which are based on patches defined by time since disturbance, such as the POP or ORCHIDEEMICT models (Haverd et al., 2014; Yue et al., 2018). Secondly, we assume that plant growth rates vary as a power of plant mass. By default we assume a power of ${\mathit{\varphi}}_{\text{g}}=\mathrm{3}/\mathrm{4}$, which is consistent with metabolic scaling theory (Enquist et al., 1998) and the empirically determined allometric relationships of Niklas and Spatz (2004).
Finally, we assume that competition is only significant for the lowest “seedling” mass class. This enables us to represent gap dynamics among plants and resultant stages in succession. This represents a significant simplification compared to other approaches involving the perfect plasticity assumption (PPA), as used within DGVMs such as LM3PPA or CLM(ED) (Fisher et al., 2015; Weng et al., 2015), where canopies are assumed to perfectly fill gaps through photomorphism (Strigul et al., 2008). In LM3PPA the radiative flux is limited by the available gap fraction in a given crown layer. PPA parallels our gap boundary condition at the lowest mass class (Eq. 6), but in RED the growth of a cohort is given by the disaggregation of total growth via metabolic scaling (Eq. 11).
These simplifications allow RED to be solved analytically for the steadystate vegetation cover given information on the mortality and growth rates per unit area for each PFT. Such analytical steadystate solutions mean that RED can be easily initialized in driftfree preindustrial states, which is vital to avoid spurious sources and sinks in climate–carbon cycle projections. The analytical solutions also enable RED to be calibrated to the observed vegetation cover, via a single parameter (μ_{0}) which represents the ratio of mortality to growth for a tree of an arbitrary reference mass. The existence of analytical steadystate solutions for RED also opens up other promising research avenues. For example, these solutions imply relationships between the fractional coverage of each PFT, total plant biomass, and the ratio of mortality to growth. This in turn allows RED to be calibrated using observations of any two of these quantities. The analytical solutions also allow optimality hypotheses to be explored (e.g. the hypothesis that the fraction of net assimilate allocated to seed production maximizes stand density and/or biomass).
Aside from the existence of analytical steadystate solutions, RED is attractive for largescale applications because it is both parameter sparse (“parsimonious”) and requires very few driving variables. The main driving variable is the timevarying net plant growth rate for each PFT, which is defined as net primary production minus the local litterfall. These driving data can be provided by a landsurface scheme, as we do in this study, or from observations. The only other driving variable for RED is the mortality rate, which we treat in this study as a geographically varying PFTspecific constant that is independent of mass. However, in principle RED could utilize mortality rates that depend on plant mass and time to represent individual disturbance events (e.g. forest fires, disease outbreaks). Despite its simplicity, the RED model is able to fit the global distribution of vegetation types (Fig. 5) and simulate successional dynamics, including changes in forest demography (Fig. 10).
There are inevitably weaknesses with any particular modelling approach. For RED, a current limitation is for competition to lead to a single PFT at each location within each cocompeting vegetation class (i.e. tree, shrub, grass). The PFT with the highest equilibrium fraction will end up excluding subdominant PFTs within the same vegetation class. It was necessary for us to account for this eventual competitive exclusion to derive zerodrift steady states for the global runs presented in Sect. 3.2.1. Such competitive exclusion is a common problem in DGVMs (Fisher et al., 2018). Currently, RED would therefore not be the most appropriate DGVM to answer important questions regarding the role of biodiversity in ecosystem function (Pavlick et al., 2013; Levine et al., 2016). More sophisticated DGVMs are required to simulate plant diversity, such as individualbased models (Fischer et al., 2016) and DGVMs specifically designed to capture subgridscale patch dynamics (Longo et al., 2019a, b). Adapting our gap boundary condition (Eq. 7) appears to be a promising way to allow greater PFT diversity in RED, without unduly increasing model complexity. We see this as a key priority for future research.
RED is currently being coupled to the JULES landsurface model, replacing TRIFFID as the default DGVM within that framework. In parallel, significant improvements are being made to the representation of physiological processes in JULES, most notably through the representation of nonstructural carbohydrate (“SUGAR”; Jones et al., 2019), and through the inclusion of a coupled model of stomatal conductance and hydraulic failure under drought stress (“SOX”; Eller et al., 2018, 2020). Plans are also being made to derive the mortality rates for RED from the INFERNO forestfire model (Burton et al., 2019). These developments will allow us to simulate the effects of sizedependent tree mortality rates within the near future.
In this paper we have presented a new intermediate complexity secondgeneration dynamic global vegetation model (DGVM), which captures important changes in forest demography. The Robust Ecosystem Demography (RED) model makes a number of important simplifications to achieve this. These simplifications are based on theoretical concepts (e.g. metabolic scaling theory to estimate how plant growth rate varies with plant mass and minimum crown overlap) and also comparison to observed forest demography (Moore et al., 2018, 2020). As a result, RED is parameter sparse and can be driven with time series of net plant growth rate (and optionally disturbance rates) for each plant functional type (PFT). We have demonstrated that RED can be calibrated effectively to observed global vegetation maps, using a single fitting parameter (representing the ratio of mortality to growth for a plant of an arbitrary reference mass). The next stage will be to use RED in coupled climate–carbon cycle projections so to assess how changes in vegetation demography impact future CO_{2} and climate. We have made the prototype RED code publicly available, and we hope that Earth system and landsurface modellers will make good use of this framework to further their own research.
For largescale application in ESMs, a primary concern is to ensure that the total vegetation carbon obeys carbon balance (i.e. only changes due to the net impact of total growth minus total mortality). Here we use that requirement to derive the functional form for F_{i} as given in Eq. (5).
The total vegetation carbon in each mass class is M_{i}=m_{i}N_{i}. The updated equation for M_{i} is therefore Eq. (4) multiplied by m_{i}:
The total carbon in the vegetation, M, is the sum of the carbon in each of the mass classes:
Thus the update equation for the total carbon is
which can be rewritten as
Now substituting Eq. (5) into Eq. (A4) gives
The first term on the righthand side of this equation is the total carbon uptake due to growth, and the second term represents the total carbon loss due to mortality, which is the required carbon conservation equation.
Equation (1) can be solved for the steady state if we assume metabolic scaling of growth using Eq. (2) and a sizeindependent mortality (Moore et al., 2018, 2020) as follows:
where n_{0} is a boundary condition that describes the number density at the mass m_{0}. The parameter μ_{0} is the ratio of the rate biomass loss due to mortality to the rate of biomass gain due to growth, for the reference mass class m_{0}. Similar analytical solutions can be derived for other measures of tree size, such as basal diameter or height (Moore et al., 2018, 2020).
Integrating Eq. (B1) from m_{0} to ∞ gives the total number density:
Other cohort integrals can be derived by integrating over the number density distribution, such as total growth rate (∫gndm):
total biomass (∫mndm):
and total vegetation cover (∫andm):
where Γ(a,b) is the incomplete upper gamma function. As we assume the allometric exponents presented in Niklas and Spatz (2004) (${\mathit{\varphi}}_{\text{g}}=\mathrm{3}/\mathrm{4}$, ${\mathit{\varphi}}_{\text{a}}=\mathrm{1}/\mathrm{3}$), these functional forms simplify to
Finally, to convert a μ_{0} found using biomass (μ_{0,tdm}) to one based on carbon mass, we use the formula
assuming that biomass is twice the carbon mass.
B1 Closed continuous form
The lowest population flux, n_{0}g_{0}, is equal to the seedling boundary condition, F_{0}, in Eq. (6):
Substituting the total number density, N_{eq} (Eq. B2), into the lefthand side and total growth, G_{eq} (Eq. B6), into the righthand side yields a solution for the equilibrium coverage, assuming $s=\mathrm{1}{\mathit{\nu}}_{\mathrm{eq}}$, as follows:
which simplifies to
Using Eq. (B8) we can write the total number density at equilibrium in terms of ν_{eq}:
This enables Eq. (B6) to be rewritten as
This equation in turn defines the total assimilate:
Finally the total biomass can be written in closed form as
B2 Discrete steady state
To solve for the discrete model equilibrium, we start from the flow equation from Eq. (4) with the term $\partial N/\partial t\to \mathrm{0}$:
Considering the population flux (Eq. 5), we find N_{i} in relation to the lower mass class, N_{i−1}:
Assuming no population grows out of the top class, λ_{I} is given as
λ_{i} can be simplified to depend only on μ_{0}, by using ${\mathit{\mu}}_{\mathrm{0}}=(\mathit{\gamma}{m}_{\mathrm{0}}/{g}_{\mathrm{0}})$ (Eq. 14) and applying the mass scaling of growth rates ${g}_{i}={g}_{\mathrm{0}}({m}_{i}/{m}_{\mathrm{0}}{)}_{\text{g}}^{\mathit{\varphi}}$. We can show that λ_{i} and λ_{I} are
An expression for the total stand density at equilibrium, N_{eq}, can be derived. Using Eq. (B18), we can represent any population of mass class i in terms of the lowest mass class N_{0} as follows:
Therefore, when finding the total number of stands relative to N_{0} we get
where X_{N} describes the sum of the all mass classes as a proportion of N_{0}. We can describe the total class growth rate in relation to N_{0} as
By using the allometric relationship (Eq. 2)
we describe the total class growth rate in relation to the lowest class growth rate, N_{0}g_{0}. Like N_{eq}, we can show the total growths across all classes is therefore
We can repeat the same process for coverage
and using allometric relationship (Eq. 3)
This gives the total coverage, ν_{eq}, as
Finally, for the total carbon mass within the class, we get
with the total carbon density equalling
In equilibrium, the rate of the recruitment of seedlings (Eq. 6) must balance the rate of loss of plants due to total mortality (γsN_{eq}):
Substituting in Eq. (B22), Eq. (B25) yields a balance equation for the kth PFT:
We can get the equilibrium fraction of a PFT, k, by rearranging the above equation, assuming c_{kk}=1 as follows:
Once the value of μ_{0} has been derived in this manner, we can find N_{0} by inversion of Eq. (B28) as follows:
Substituting Eq. (B33) into Eq. (B34) allows us to determine N_{0} and hence most other total densities in terms of purely μ_{0} and prescribed constants.
B3 Continuous–discrete convergence
Inevitably discretized models will not exactly reproduce exact continuum analytical solutions, as a result of numerical inaccuracies that arise from using a finite number of mass classes. However, where exact analytical solutions exist they can be used to benchmark numerical models and optimize discretization schemes, which is what we set out to do in this appendix. We compare the continuum analytical solution for the equilibrium coverage (Eq. B12) to results from RED with differing numbers of mass classes m_{i} and a geometric mass class scaling, ${m}_{i+\mathrm{1}}=\mathit{\xi}{m}_{i}$. Figure B2a shows how the relationship between ν_{eq} varies with μ_{0} for the exact continuum solution (black line) and variants of the numerical version of RED with different numbers of mass classes (coloured lines). As hoped, results from the discretized model converge on the exact solution as the number of mass classes increases.
The numerical versions of RED shown in Fig. B2a each use a value of ξ that is near optimum for the number of mass classes, as shown in panels (b) and (c) of Fig. B2. Optimum ξ values reduce from about 2.3 for 10 mass classes to 1.1 for 100 mass classes. This variation results from a tradeoff. For a given number of mass classes, small values of ξ give greater numerical accuracy but explicitly model less of the mass range, and the opposite is true of large ξ values. As a result, optimum values of ξ an be defined for each number of mass classes as outlined below.
For geometric scaling any mass can be expressed in terms of m_{0}, by writing m_{i}=m_{0}(ξ)^{i}. Therefore, by using ${m}_{i+\mathrm{1}}{m}_{i}={m}_{\mathrm{0}}\left(\mathit{\xi}{)}^{i}\right(\mathit{\xi}\mathrm{1})$, we find that our equilibrium form of λ_{i} is reduced to
From Fig. B2c, we see that there is an optimum value for ξ, the geometric scaling for a given number of classes, which minimizes the difference between the continuous and discrete forms. This can be found by taking the difference in the continuous and discrete coverages and differentiating with respect to ξ to find the minima. It should be noted that as the continuous form is not dependent on ξ, we get
where ν_{eq} corresponds with the discrete equilibrium, i.e. Eq. (B32), with ${\mathit{\nu}}_{\mathrm{eq}}=(\mathrm{1}s)$. Setting Eq. (B36) equal to zero we reduce the relationship to only a dependence on X_{N} and X_{G}:
Finding the partial derivative of X_{N}, using the geometric form of Eq. (B18), we get
and for X_{G} we get the following:
Finding ${\mathit{\lambda}}_{i}^{\prime}$ we get
and for the top class, ${\mathit{\lambda}}_{I}^{\prime}$ we get the following:
To numerically solve for the minimum, we must differentiate Eq. (B37), with respect to ξ. Through the product rule we get
Differentiating equation (Eq. B38) and simplifying gives
and doing the same for Eq. (B39) gives
${\mathit{\lambda}}_{i}^{\prime \prime}$ is given by
For the double differential of λ_{i} we get
We now possess the identities needed to numerically find the optimum bin scaling for a given number of classes. In Fig. B2c, the optimum scaling, ξ, is shown as the solid black line.
The diagnosed mortality rates in Fig. 6 are sensitive to variation in model inputs and parameters. The mortality rate, γ, can be found for the continuous solutions by rearranging the boundary condition equation (Eq. 6) and substituting in Eqs. (B2) and (B13) as follows:
The key external inputs to this equation are the observed PFT fraction ν_{eq} and the net assimilate P_{eq}. In addition, our estimates of γ are dependent on the internal model parameters, α and m_{0}.
The red lines in Fig. C1 demonstrate how the estimate of γ depends on these four inputs. The black dashed lines in Fig. C1 indicate how uncertainties in each input relate to uncertainties in γ, for “true” values typical of a tree PFT. We estimate uncertainties in the observed PFT fraction (e.g. from remote sensing) to be ±5 % and uncertainties in P (e.g. from JULES) to be ±20 %, leading to errors of ±17 % and ±20 % respectively. Likewise, ±20 % uncertainties in the internal parameters α and m_{0} lead to ±12 % and ± 20 % uncertainties in γ. Combining these sources of uncertainty leads to an overall uncertainty in our inferred estimate of γ of about ±35 %.
The RED model Python code is archived at https://doi.org/10.5281/zenodo.3548678 (Argles, 2019). Furthermore, RED is currently being coupled into JULES, where a basic integration currently exists as branch (vn5.8_red_jules_couple) – this requires registration for the JULES repository (https://code.metoffice.gov.uk/trac/home, last access: 2 September 2020).
Originally the model framework was in JRM's thesis (Moore, 2016) under the supervision of PMC and CH. The description of PFT competition, the numerical model, and the equilibrium solutions has been further developed by APKA, JRM, ABH, and PMC. Currently RED is being integrated into JULES with the supervision of AJW and CJ. AJW also provided and processed the UKESM growth rates needed to drive RED globally within this paper.
The authors declare that they have no conflict of interest.
We are grateful to the Met Office for considering implementation in JULES, via a ticket 1034 within a branch of the code repository.
This research has been supported by the Newton Fund (CSSPBrazil), the European Research Council (ECCLES), the University of Exeter/Met Office (CASE Studentship), Joint UK BEIS/Defra Met Office Hadley Centre Climate Programme (GA01101), and the European Commission Horizon 2020 research and innovation programme (grant agreement no. 641816).
This paper was edited by Carlos Sierra and reviewed by three anonymous referees.
Ahlström, A., Xia, J., Arneth, A., Luo, Y., and Smith, B.: Importance of vegetation dynamics for future terrestrial carbon cycling, Environ. Res. Lett., 10, 054019, https://doi.org/10.1088/17489326/10/5/054019, 2015. a
Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., McDowell, N., Vennetier, M., Kitzberger, T., Rigling, A., Breshears, D. D., Hogg, E. H., Gonzalez, P., Fensham, R., Zhang, Z., Castro, J., Demidova, N., Lim, J.H., Allard, G., Running, S. W., Semerci, A., and Cobb, N.: A global overview of drought and heatinduced tree mortality reveals emerging climate change risks for forests, Forest Ecol. Manage., 259, 660–684, 2010. a
Argles, A. P. K.: aargles/RED_DGVM: RED DGVM, version 1.0.1, Zenodo, https://doi.org/10.5281/zenodo.3548678, 2019. a
Arora, V. K., Katavouta, A., Williams, R. G., Jones, C. D., Brovkin, V., Friedlingstein, P., Schwinger, J., Bopp, L., Boucher, O., Cadule, P., Chamberlain, M. A., Christian, J. R., Delire, C., Fisher, R. A., Hajima, T., Ilyina, T., Joetzjer, E., Kawamiya, M., Koven, C., Krasting, J., Law, R. M., Lawrence, D. M., Lenton, A., Lindsay, K., Pongratz, J., Raddatz, T., Séférian, R., Tachiiri, K., Tjiputra, J. F., Wiltshire, A., Wu, T., and Ziehn, T.: Carbonconcentration and carbonclimate feedbacks in CMIP6 models, and their comparison to CMIP5 models, Biogeosciences Discuss., https://doi.org/10.5194/bg2019473, in review, 2019. a
Brando, P. M., Balch, J. K., Nepstad, D. C., Morton, D. C., Putz, F. E., Coe, M. T., Silvério, D., Macedo, M. N., Davidson, E. A., Nóbrega, C. C., Alencar, A., and SoaresFilho, B. S.: Abrupt increases in Amazonian tree mortality due to drought–fire interactions, P. Natl. Acad. Sciences, 111, 6347–6352, 2014. a, b
Brovkin, V., Boysen, L., Arora, V. K., Boisier, J. P., Cadule, P., Chini, L., Claussen, M., Friedlingstein, P., Gayler, V., Van den hurk, B. J., Hurtt, G. C., Jones, C. D., Kato, E., De noblet ducoudre, N., Pacifico, F., Pongratz, J., and Weiss, M.: Effect of anthropogenic landuse and landcover changes on climate and land carbon storage in CMIP5 projections for the twentyfirst century, J. Climate, 26, 6859–6881, https://doi.org/10.1175/JCLID1200623.1, 2013. a, b
Burton, C., Betts, R., Cardoso, M., Feldpausch, T. R., Harper, A., Jones, C. D., Kelley, D. I., Robertson, E., and Wiltshire, A.: Representation of fire, landuse change and vegetation dynamics in the Joint UK Land Environment Simulator vn4.9 (JULES), Geosci. Model Dev., 12, 179–193, https://doi.org/10.5194/gmd121792019, 2019. a, b
Chaturvedi, R., Raghubanshi, A., and Singh, J.: Effect of grazing and harvesting on diversity, recruitment and carbon accumulation of juvenile trees in tropical dry forests, Forest Ecol. Manage., 284, 152–162, 2012. a
Chuvieco, E., Pettinari, M., Lizundia Loiola, J., Storm, T., and Padilla Parellada, M.: ESA Fire Climate Change Initiative (Fire_cci): MODIS Fire_cci Burned Area Grid product, version 5.1, https://doi.org/10.5285/3628cb2fdba443588155e15dee8e5352, 2019. a
Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., Chhabra, A., DeFries, R., Galloway, J., Heimann, M., Jones, C., Le Quéré, C., Myneni, R. B., Piao, S., and Thornton, P.: Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 465–570, 2014. a
Coomes, D. A., Duncan, R. P., Allen, R. B., and Truscott, J.: Disturbances prevent stem sizedensity distributions in natural forests from following scaling relationships, Ecol. Lett., 6, 980–989, https://doi.org/10.1046/j.14610248.2003.00520.x., 2003. a
Coutinho, L. M.: Fire in the ecology of the Brazilian cerrado, in: Fire in the tropical biota, Springer, 82–105, 1990. a
Cox, P. M.: Description of the “TRIFFID” dynamic global vegetation model, Hadley Centre for Climate Prediction and Research, Hadley Centre, Met Office, Bracknell, UK, Technical note, no. 24, 2001. a, b, c
Curran, L. M., Trigg, S. N., McDonald, A. K., Astiani, D., Hardiono, Y., Siregar, P., Caniago, I., and Kasischke, E.: Lowland forest loss in protected areas of Indonesian Borneo, Science, 303, 1000–1003, 2004. a
Dorrough, J. and Moxham, C.: Eucalypt establishment in agricultural landscapes and implications for landscapescale restoration, Biol. Conserv., 123, 55–66, 2005. a
Eller, C. B., Rowland, L., Oliveira, R. S., Bittencourt, P. R. L., Barros, F. V., da Costa, A. C. L., Meir, P., Friend, A. D., Mencuccini, M., Sitch, S., and Cox, P.: Modelling tropical forest responses to drought and El Niño with a stomatal optimization model based on xylem hydraulics, Philos. T. Roy. Soc. B, 373, 20170315, https://doi.org/10.1098/rstb.2017.0315, 2018. a
Eller, C. B., Rowland, L., Mencuccini, M., Rosas, T., Williams, K., Harper, A., Medlyn, B. E., Wagner, Y., Klein, T., Teodoro, G. S., Oliveira, R. S., Matos, I. S., Rosado, B. H. P., Fuchs, K., Wohlfahrt, G., Montagnani, L., Meir, P., Sitch, S., and Cox, P. M.: Stomatal optimisation based on xylem hydraulics (SOX) improves land surface model simulation of vegetation responses to climate, New Phytologist, 226, 1622–1637, 2020. a
Enquist, B. J., Brown, J. H., and West, G. B.: Allometric scaling of plant energetics and population density, Nature, 395, 163–165, 1998. a, b
Esseen, P.A.: Tree mortality patterns after experimental fragmentation of an oldgrowth conifer forest, Biol. Conserv., 68, 19–28, 1994. a
Fischer, R., Bohn, F., de Paula, M. D., Dislich, C., Groeneveld, J., Gutiérrez, A. G., Kazmierczak, M., Knapp, N., Lehmann, S., Paulick, S., Pütz, S., Rödig, E., Taubert, F., Köhler, P., Huth, A.: Lessons learned from applying a forest gap model to understand ecosystem and carbon dynamics of complex tropical forests, Ecol. Model., 326, 124–133, 2016. a
Fisher, R., McDowell, N., Purves, D., Moorcroft, P., Sitch, S., Cox, P., Huntingford, C., Meir, P., and Woodward, F. I.: Assessing uncertainties in a secondgeneration dynamic vegetation model caused by ecological scale limitations, New Phytologist, 187, 666–681, 2010. a
Fisher, R. A., Muszala, S., Verteinstein, M., Lawrence, P., Xu, C., McDowell, N. G., Knox, R. G., Koven, C., Holm, J., Rogers, B. M., Spessa, A., Lawrence, D., and Bonan, G.: Taking off the training wheels: the properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED), Geosci. Model Dev., 8, 3593–3619, https://doi.org/10.5194/gmd835932015, 2015. a
Fisher, R. A., Koven, C. D., Anderegg, W. R., Christoffersen, B. O., Dietze, M. C., Farrior, C. E., Holm, J. A., Hurtt, G. C., Knox, R. G., Lawrence, P. J., Lichstein, J. W., Longo, M., Matheny, A. M., Medvigy, D., MullerLandau, H. C., Powell, T. L., Serbin, S. P., Sato, H., Shuman, J. K., Smith, B., Trugman, A. T., Viskari, T., Verbeeck, H., Weng, E., Xu, C., Xu, X., Zhang, T., and Moorcroft, P. R.: Vegetation demographics in Earth System Models: A review of progress and priorities, Glob. Change Biol., 24, 35–54, https://doi.org/10.1111/gcb.13910, 2018. a, b, c
Friedlingstein, P., Meinshausen, M., Arora, V. K., Jones, C. D., Anav, A., Liddicoat, S. K., and Knutti, R.: Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks, J. Climate, 27, 511–526, https://doi.org/10.1175/JCLID1200579.1, 2014. a
Haddad, N. M., Brudvig, L. A., Clobert, J., Davies, K. F., Gonzalez, A., Holt, R. D., Lovejoy, T. E., Sexton, J. O., Austin, M. P., Collins, C. D., Cook, W. M., Damschen, E. I., Ewers, R. M., Foster, B. L., Jenkins, C. N., King, A. J., Laurance, W. F., Levey, D. J., Margules, C. R., Melbourne, B. A., Nicholls, A. O., Orrock, J. L., Song, D.X., and Townshend, J. R.: Habitat fragmentation and its lasting impact on Earth's ecosystems, Sci. Adv., 1, e1500052, https://doi.org/10.1126/sciadv.1500052, 2015. a
Harper, A. B., Cox, P. M., Friedlingstein, P., Wiltshire, A. J., Jones, C. D., Sitch, S., Mercado, L. M., Groenendijk, M., Robertson, E., Kattge, J., Bönisch, G., Atkin, O. K., Bahn, M., Cornelissen, J., Niinemets, Ü., Onipchenko, V., Peñuelas, J., Poorter, L., Reich, P. B., Soudzilovskaia, N. A., and Bodegom, P. V.: Improved representation of plant functional types and physiology in the Joint UK Land Environment Simulator (JULES v4.2) using plant trait information, Geosci. Model Dev., 9, 2415–2440, https://doi.org/10.5194/gmd924152016, 2016. a
Harper, A. B., Wiltshire, A. J., Cox, P. M., Friedlingstein, P., Jones, C. D., Mercado, L. M., Sitch, S., Williams, K., and DuranRojas, C.: Vegetation distribution and terrestrial carbon cycle in a carbon cycle configuration of JULES4.6 with new plant functional types, Geosci. Model Dev., 11, 2857–2873, https://doi.org/10.5194/gmd1128572018, 2018. a
Haverd, V., Smith, B., Nieradzik, L. P., and Briggs, P. R.: A standalone tree demography and landscape structure module for Earth system models: integration with inventory data from temperate and boreal forests, Biogeosciences, 11, 4039–4055, https://doi.org/10.5194/bg1140392014, 2014. a, b
Higuchi, P., OliveiraFilho, A. T., Bebber, D. P., Brown, N. D., Silva, A. C., and Machado, E. L.: Spatiotemporal patterns of tree community dynamics in a tropical forest fragment in Southeast Brazil, Plant Ecol., 199, 125–135, 2008. a
Jones, S., Rowland, L., Cox, P., Hemming, D., Wiltshire, A., Williams, K., Parazoo, N. C., Liu, J., da Costa, A. C. L., Meir, P., Mencuccini, M., and Harper, A. B.: The impact of a simple representation of nonstructural carbohydrates on the simulated response of tropical forests to drought, Biogeosciences, 17, 3589–3612, https://doi.org/10.5194/bg1735892020, 2020. a
Jönsson, M. T., Fraver, S., Jonsson, B. G., Dynesius, M., Rydgård, M., and Esseen, P.A.: Eighteen years of tree mortality and structural change in an experimentally fragmented Norway spruce forest, Forest Ecol. Manage., 242, 306–313, 2007. a
Kinnaird, M. F. and O’Brien, T. G.: Ecological effects of wildfire on lowland rainforest in Sumatra, Conserv. Biol., 12, 954–956, 1998. a
Kohyama, T., Suzuki, E., Partomihardjo, T., Yamada, T., and Kubo, T.: Tree species differentiation in growth, recruitment and allometry in relation to maximum height in a Bornean mixed dipterocarp forest, J. Ecol., 91, 797–806, 2003. a
Laurance, W. F., Ferreira, L. V., Rankinde Merona, J. M., and Laurance, S. G.: Rain forest fragmentation and the dynamics of Amazonian tree communities, Ecology, 79, 2032–2040, 1998. a
Levine, N. M., Zhang, K., Longo, M., Baccini, A., Phillips, O. L., Lewis, S. L., AlvarezDávila, E., de Andrade, A. C. S., Brienen, R. J., Erwin, T. L., Feldpausch, T. R., Monteagudo Mendoza, A. L., Nuñez Vargas, P., Prieto, A., SilvaEspejo, J. E., Malhi, Y., and Moorcroft, P, R.: Ecosystem heterogeneity determines the ecological resilience of the Amazon to climate change, P. Natl. Acad. Sci. USA, 113, 793–797, 2016. a
Li, W., MacBean, N., Ciais, P., Defourny, P., Lamarche, C., Bontemps, S., and Peng, S.: Derivation of plant functional type (PFT) maps from the ESA CCI Land Cover product, Data set, Zenodo, https://doi.org/10.5281/zenodo.1048163, 2019. a
Lima, R. A., MullerLandau, H. C., Prado, P. I., and Condit, R.: How do size distributions relate to concurrently measured demographic rates? Evidence from over 150 tree species in Panama, J. Trop. Ecol., 32, 179–192, 2016. a
Longo, M., Knox, R. G., Medvigy, D. M., Levine, N. M., Dietze, M. C., Kim, Y., Swann, A. L. S., Zhang, K., Rollinson, C. R., Bras, R. L., Wofsy, S. C., and Moorcroft, P. R.: The biophysics, ecology, and biogeochemistry of functionally diverse, vertically and horizontally heterogeneous ecosystems: the Ecosystem Demography model, version 2.2 – Part 1: Model description, Geosci. Model Dev., 12, 4309–4346, https://doi.org/10.5194/gmd1243092019, 2019a. a
Longo, M., Knox, R. G., Levine, N. M., Swann, A. L. S., Medvigy, D. M., Dietze, M. C., Kim, Y., Zhang, K., Bonal, D., Burban, B., Camargo, P. B., Hayek, M. N., Saleska, S. R., da Silva, R., Bras, R. L., Wofsy, S. C., and Moorcroft, P. R.: The biophysics, ecology, and biogeochemistry of functionally diverse, vertically and horizontally heterogeneous ecosystems: the Ecosystem Demography model, version 2.2 – Part 2: Model evaluation for tropical South America, Geosci. Model Dev., 12, 4347–4374, https://doi.org/10.5194/gmd1243472019, 2019b. a
Lugo, A. E. and Scatena, F. N.: Background and catastrophic tree mortality in tropical moist, wet, and rain forests, Biotropica, 28, 585–599, 1996. a
Medeiros, M. and Miranda, H.: Postfire resprouting and mortality in cerrado woody plant species over a threeyear period, Edinburgh Journal of Botany, 65, 53–68, 2008. a
Medvigy, D., Wofsy, S. C., Munger, J. W., Hollinger, D. Y., and Moorcroft, P. R.: Mechanistic scaling of ecosystem function and dynamics in space and time: Ecosystem Demography model version 2, J. Geophys. Res.Biogeo., 114, G01002, https://doi.org/10.1029/2008JG000812, 2009. a, b
Moorcroft, P. R., Hurtt, G., and Pacala, S. W.: a Method for Scaling Vegetation Dynamics: the Ecosystem Demography Model (ED), Ecol. Monogr., 71, 557–586, https://doi.org/10.1890/00129615(2001)071[0557:AMFSVD]2.0.CO;2, 2001. a
Moore, J. R.: Aspects of Land Surface Modelling: Role of Biodiversity in Ecosystem Resilience to Environmental Change and a Robust Ecosystem Demography Model, PhD thesis, University of Exeter, 2016. a
Moore, J. R., Zhu, K., Huntingford, C., and Cox, P. M.: Equilibrium forest demography explains the distribution of tree sizes across North America, Environ. Res. Lett., 13, 084019, https://doi.org/10.1088/17489326/aad6d1, 2018. a, b, c, d, e, f, g, h, i
Moore, J. R., Argles, A. P. K., Zhu, K., Huntingford, C., and Cox, P. M.: Validation of demographic equilibrium theory against treesize distributions and biomass density in Amazonia, Biogeosciences, 17, 1013–1032, https://doi.org/10.5194/bg1710132020, 2020. a, b, c, d, e, f, g, h, i, j
MullerLandau, H. C., Condit, R. S., Harms, K. E., Marks, C. O., Thomas, S. C., Bunyavejchewin, S., Chuyong, G., Co, L., Davies, S., Foster, R., Gunatilleke, S., Gunatilleke, N., Hart, T., Hubbell, S. P., Itoh, A., Kassim, A. R., Kenfack, D., LaFrankie, J. V., Lagunzad, D., Lee, H. S., Losos, E., Makana, J. R., Ohkubo, T., Samper, C., Sukumar, R., Sun, I. F., Nur Supardi, M. N., Tan, S., Thomas, D., Thompson, J., Valencia, R., Vallejo, M. I., Muñoz, G. V., Yamakura, T., Zimmerman, J. K., Dattaraja, H. S., Esufali, S., Hall, P., He, F., Hernandez, C., Kiratiprayoon, S., Suresh, H. S., Wills, C., and Ashton, P.: Comparing tropical forest tree size distributions with the predictions of metabolic ecology and equilibrium models, Ecol. Lett., 9, 589–602, https://doi.org/10.1111/j.14610248.2006.00915.x, 2006. a
Nepstad, D. C., Stickler, C. M., Filho, B. S., and Merry, F.: Interactions among Amazon land use, forests and climate: prospects for a nearterm forest tipping point, Philos. T. Roy. Soc. B, 363, 1737–1746, 2008. a
Niklas, K. J. and Enquist, B. J.: Invariant scaling relationships for interspecific plant biomass production rates and body size, P. Natl. Acad. Sci. USA, 98, 2922–2927, 2001. a
Niklas, K. J. and Spatz, H.C.: From The Cover: Growth and hydraulic (not mechanical) constraints govern the scaling of tree height and mass, P. Natl. Acad. Sci. USA, 101, 15661–15663, https://doi.org/10.1073/pnas.0405857101, 2004. a, b, c, d, e
Pavlick, R., Drewry, D. T., Bohn, K., Reu, B., and Kleidon, A.: The Jena DiversityDynamic Global Vegetation Model (JeDiDGVM): a diverse approach to representing terrestrial biogeography and biogeochemistry based on plant functional tradeoffs, Biogeosciences, 10, 4137–4177, https://doi.org/10.5194/bg1041372013, 2013. a
Peterson, D. W. and Reich, P. B.: Prescribed fire in oak savanna: fire frequency effects on stand structure and dynamics, Ecol. Appl., 11, 914–927, 2001. a
Phillips, O. L.: Longterm environmental change in tropical forests: increasing tree turnover, Environ. Conserv., 23, 235–248, https://doi.org/10.1017/s0376892900038856, 1996. a, b, c
Phillips, O. L., Baker, T. R., Arroyo, L., Higuchi, N., Killeen, T. J., Laurance, W. F., Lewis, S. L., Lloyd, J., Malhi, Y., Monteagudo, A., Neill, D. A., Nüñez Vargas, P., Silva, J. N., Terborgh, J., Vásquez Martínez, R., Alexiades, M., Almeida, S., Brown, S., Chave, J., Comiskey, J. A., Czimczik, C. I., Di Fiore, A., Erwin, T., Kuebler, C., Laurance, S. G., Nascimento, H. E., Olivier, J., Palacios, W., Patiño, S., Pitman, N. C., Quesada, C. A., Saldias, M., Torres Lezama, A., and Vinceti, B.: Pattern and process in Amazon tree turnover, 1976–2001, Philos. T. Roy. Soc. B, 359, 381–407, https://doi.org/10.1098/rstb.2003.1438, 2004. a, b, c
Poulter, B., MacBean, N., Hartley, A., Khlystova, I., Arino, O., Betts, R., Bontemps, S., Boettcher, M., Brockmann, C., Defourny, P., Hagemann, S., Herold, M., Kirches, G., Lamarche, C., Lederer, D., Ottlé, C., Peters, M., and Peylin, P.: Plant functional type classification for earth system models: results from the European Space Agency's Land Cover Climate Change Initiative, Geosci. Model Dev., 8, 2315–2328, https://doi.org/10.5194/gmd823152015, 2015. a
Prior, L. D., Murphy, B. P., and RussellSmith, J.: Environmental and demographic correlates of tree recruitment and mortality in north Australian savannas, Forest Ecol. Manage., 257, 66–74, 2009. a, b
Pugh, T. A. M., Rademacher, T., Shafer, S. L., Steinkamp, J., Barichivich, J., Beckage, B., Haverd, V., Harper, A., Heinke, J., Nishina, K., Rammig, A., Sato, H., Arneth, A., Hantson, S., Hickler, T., Kautz, M., Quesada, B., Smith, B., and Thonicke, K.: Understanding the uncertainty in global forest carbon turnover, Biogeosciences, 17, 3961–3989, https://doi.org/10.5194/bg1739612020, 2020. a
Sato, H., Itoh, A., and Kohyama, T.: SEIBDGVM: A new Dynamic Global Vegetation Model using a spatially explicit individualbased approach, Ecol. Model., 200, 279–307, https://doi.org/10.1016/j.ecolmodel.2006.09.006, 2007. a
Scheiter, S., Langan, L., and Higgins, S. I.: Nextgeneration dynamic global vegetation models: Learning from community ecology, New Phytologist, 198, 957–969, https://doi.org/10.1111/nph.12210, 2013. a
Schleussner, C.F., Rogelj, J., Schaeffer, M., Lissner, T., Licker, R., Fischer, E. M., Knutti, R., Levermann, A., Frieler, K., and Hare, W.: Science and policy characteristics of the Paris Agreement temperature goal, Nat. Clim. Change, 6, 827–835, 2016. a
Sellar, A. A., Jones, C. G., Mulcahy, J., Tang, Y., Yool, A., Wiltshire, A., O'connor, F. M., Stringer, M., Hill, R., Palmieri, J., Woodward, S., de Mora, L., Kuhlbrodt, T., Rumbold, S. T., Kelley, D. I., Ellis, R., Johnson, C. E., Walton, J., Abraham, N. L., Andrews, M. B., Andrews, T., Archibald, A. T., Berthou, S., Burke, E., Blockley, E., Carslaw, K., Dalvi, M., Edwards, J., Folberth, G. A., Gedney, N., Griffiths, P. T., Harper, A. B., Hendry, M. A., Hewitt, A. J., Johnson, B., Jones, A., Jones, C. D., Keeble, J., Liddicoat, S., Morgenstern, O., Parker, R. J., Predoi, V., Robertson, E., Siahaan, A., Smith, R. S., Swaminathan, R., Woodhouse, M. T., Zeng, G., and Zerroukat, M.: UKESM1: Description and evaluation of the UK Earth System Model, J. Adv. Model. Earth Sy., 11, 45134558, 2019. a, b
Smith, B.: LPJGUESSan ecosystem modelling framework, Department of Physical Geography and Ecosystems Analysis, INES, Sölvegatan, 12, 22362, 2001. a
Staver, A. C., Bond, W. J., Stock, W. D., Van Rensburg, S. J., and Waldram, M. S.: Browsing and fire interact to suppress tree density in an African savanna, Ecol. Appl., 19, 1909–1919, 2009. a, b
Strigul, N., Pristinski, D., Purves, D., Dushoff, J., and Pacala, S.: Scaling from trees to forests: tractable macroscopic equations for forest dynamics, Ecol. Monogr., 78, 523–545, 2008. a
Swaine, M.: Characteristics of dry forest in West Africa and the influence of fire, J. Veg. Sci., 3, 365–374, 1992. a
Van Nieuwstadt, M. G. and Sheil, D.: Drought, fire and tree survival in a Borneo rain forest, East Kalimantan, Indonesia, J. Ecol., 93, 191–201, 2005. a
Van Uytvanck, J., Maes, D., Vandenhaute, D., and Hoffmann, M.: Restoration of woodpasture on former agricultural land: the importance of safe sites and time gaps before grazing for tree seedlings, Biol. Conserv., 141, 78–88, 2008. a
Weng, E. S., Malyshev, S., Lichstein, J. W., Farrior, C. E., Dybzinski, R., Zhang, T., Shevliakova, E., and Pacala, S. W.: Scaling from individual trees to forests in an Earth system modeling framework using a mathematically tractable model of heightstructured competition, Biogeosciences, 12, 2655–2694, https://doi.org/10.5194/bg1226552015, 2015. a
West, G. B., Brown, J. H., and Enquist, B. J.: A general model for the origin of allometric scaling laws in biology, Science, 276, 122–126, 1997. a
Yue, C., Ciais, P., Luyssaert, S., Li, W., McGrath, M. J., Chang, J., and Peng, S.: Representing anthropogenic gross land use change, wood harvest, and forest age dynamics in a global vegetation model ORCHIDEEMICT v8.4.2, Geosci. Model Dev., 11, 409–428, https://doi.org/10.5194/gmd114092018, 2018. a, b
 Abstract
 Introduction
 Description of the model
 Modelling results
 Discussion
 Conclusions
 Appendix A: Functional form of flux F_{i} in discretized RED
 Appendix B: Continuum solutions and demographic equilibrium theory
 Appendix C: Sensitivity of diagnosed mortality rates to model parameters
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Description of the model
 Modelling results
 Discussion
 Conclusions
 Appendix A: Functional form of flux F_{i} in discretized RED
 Appendix B: Continuum solutions and demographic equilibrium theory
 Appendix C: Sensitivity of diagnosed mortality rates to model parameters
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References