Articles | Volume 13, issue 9
Geosci. Model Dev., 13, 4067–4089, 2020

Special issue: Joint UK Land Environment Simulator (JULES) – configurations,...

Geosci. Model Dev., 13, 4067–4089, 2020

Model description paper 07 Sep 2020

Model description paper | 07 Sep 2020

Robust Ecosystem Demography (RED version 1.0): a parsimonious approach to modelling vegetation dynamics in Earth system models

Robust Ecosystem Demography (RED version 1.0): a parsimonious approach to modelling vegetation dynamics in Earth system models
Arthur P. K. Argles1, Jonathan R. Moore1, Chris Huntingford2, Andrew J. Wiltshire3, Anna B. Harper1, Chris D. Jones3, and Peter M. Cox1 Arthur P. K. Argles et al.
  • 1College of Engineering, Mathematics, and Physical Sciences, University of Exeter, Exeter EX4 4QF, UK
  • 2UK Centre for Ecology and Hydrology, Wallingford OX10 8BB, UK
  • 3Met Office Hadley Centre, Fitzroy Road, Exeter EX1 3PB, UK

Correspondence: Arthur P. K. Argles (, Peter M. Cox ( and Jonathan R. Moore (


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 CO2 and climate. More advanced cohort-based 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 power-law 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.

1 Introduction

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 2C (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 CO2 and nutrient fertilization effects and land-use 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 trade-off between the complexity of ecological process representation and the necessity of parsimony at scale (Fisher et al.2018). DGVMs range from the simplistic, older, top-down approaches to that of complex individual-based DGVMs. For example, in the first instance the TRIFFID model (Cox2001) 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, individual-based models can explicitly represent a multitude of biological and ecosystem processes at an individual plant level (Smith2001; Sato et al.2007). The benefit of this is that size-dependent 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 individual-based and top-down 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 stand-age cohorts as the dimension for population dynamics, every time step applying crowding and resource limited mortality rates. Another example is the ORCHIDEE-MICT (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 tree-size 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; Muller-Landau et al.2006; Lima et al.2016). We follow many other studies in assuming that tree-growth rates vary with the three-quarter power of tree mass (m3∕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 tree-size 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).

2 Description of the model

A full list of variables, parameters, and units are given in Table 1.

Table 1Model variables, parameters, and units.

Download Print Version | Download XLSX

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:

(1) n t + m n g = - γ n .

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 large-scale 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:

(2) g = g 0 m m 0 ϕ g .

Here g0 is the growth rate of a plant with the reference mass, m0. A value of ϕg=0.75 is assumed by default, consistent with the analysis of field-based 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:

(3) a = a 0 m m 0 ϕ a ,

where ϕa=0.5 by default. Solutions for n can be integrated over mass to derive the total plant number, N=0ndm, the total growth rate, G=0gndm, the total biomass, M=0mndm, and the fractional area covered ν=0andm.

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:

(4) N i t + F i - F i - 1 = - γ N i ,

where i denotes the ith mass class; Fi is the flux of plants growing out of the ith mass class and into the (i+1)th mass class; Fi−1 is the flux of plants growing out of the (i−1)th mass class and into the ith mass class; and Ni 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 Fi must take the form (see Appendix A)

(5) F i = N i g i ( m i + 1 - m i ) ,

where mi is the mean mass of a plant in the ith mass class, and gi 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 m0 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=(1-α)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 m0 is therefore

(6) F 0 = α P m 0 s = α ( 1 - α ) G m 0 s ,

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 F0 is the ratio of a fraction of the carbon assimilate allocated to reproduction, αP, and m0, multiplied by the gap area s.

Figure 1Schematic depicting the hierarchical PFT functional group regime within RED. Trees shade trees, shrubs, and grasses. Shrubs shade shrubs and grasses, while grasses only shade grasses.


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:

(7) s k = 1 - l c k l ν l ,

where νl is the area fraction of the lth PFT, and ckl 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, ckl=1. If PFT k is dominant over PFT l, ckl=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 Lotka-inspired TRIFFID model (Cox2001) and allows for the coexistence between inter-functional 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. Faster-growing 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.

Table 2Competition coefficients assumed for different plant functional groups. A more detailed example of this is given for specific PFTs in Table 3.

Download Print Version | Download XLSX

Figure 2Schematic of RED coupled to an ESM or land carbon cycle model. RED is driven by a time series of net carbon assimilate, P, which is then split between seedling production, αP, and the growth of existing plants, G=(1-α)P. The seedling flux is limited by the available free space, s. Additional mortality rates diagnosed from disturbance models, γd, can be added onto an assumed baseline mortality, γb, as a function of both PFT and mass class.


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

(8) P = Π N - Λ l .

We apply the m3∕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 m3∕4 (Enquist et al.1998; Niklas and Enquist2001). 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:

(9) γ = γ b + γ d .

Disturbance rates γd can in principle be both PFT dependent and mass dependent (e.g. to capture forestry practices).

Table 3List of PFT names and assumed allometric scaling parameters (m0,a0,h0), seedling fraction (α), and competition coefficient (cpft,j). The growth allometry of trees and shrubs across size is assumed to follow Niklas and Spatz (2004) (ϕg=0.75, ϕa=0.5, ϕh=0.25). The competition coefficients given describe which PFT functional group shades the current PFT; if cpft,j=1, the PFT is shaded; otherwise it is not (Table 2).

Download Print Version | Download XLSX

The input values of net assimilate for each PFT (P) define the total structural growth rate, G=(1-α)P, and the seedling flux F0 (via Eq. 6), using PFT-specific values of the parameter α (see Table 3). The definition of the total structural growth rate at a given time step is

(10) G = i N i g i ,

which can be combined with the growth scaling given by Eq. (2) to derive the reference growth rate, g0, from the net assimilate, P, which is a driving input:

(11) g 0 = ( 1 - α ) P i N i m i m 0 ϕ g .

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 (Ni) is updated using a discretized form of Eq. (4):

(12) N i ( j + 1 ) = N i ( j ) + Δ t F i - 1 ( j ) - F i ( j ) - γ ( j ) N i ( j ) ,

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 (m0/g0 4 years) and plant mortality (1/γ 20 years). The lower-boundary seedling flux is calculated from Eq. (6) using Eq. (7). We impose a zero-flux 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:

(13) Λ d = α P ( 1 - s ) + i γ i M i + g I N I .

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 neq(m), the total plant number (Neq), biomass (Meq), growth rate (Geq), 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 m0 as follows:

(14) μ 0 = γ m 0 g 0 .

In order to initialize the numerical RED model in a drift-free 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.

  1. The total equilibrium stand density, Neq:

    (15) N eq = N 0 X N .
  2. The total equilibrium structural growth, Geq:

    (16) G eq = i = 0 I N i g i = N 0 g 0 X G .
  3. The total equilibrium coverage, νeq:

    (17) ν eq = i = 0 I N i a i = N 0 a 0 X ν .
  4. The total equilibrium carbon mass, Meq:

    (18) M eq = i = 0 I N i m i = N 0 m 0 X M .

Here XN, XG, Xν, and XM are functions of μ0 (see Appendix B2). This equilibrium state is derived by setting Ni(j+1)=Ni(j) 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.

Figure 3Observation-based dataset of the PFT area fractions for the nine JULES PFTs (Harper et al.2016) as listed in Table 3.

Figure 4Mean net assimilate P assimilate (Eq. 8) from UKESM between 2000 and 2010. The mean is constructed by setting any negative growth rates to zero.

The equations above therefore define the equilibrium state of the discrete system for given values of N0 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 Neq (Eq. 15) and Geq (Eq. 16) and requiring that the recruitment flux (α/(1-α)Geqs) is equal to that of the total population dying (γNeq), we can derive an equation for the total equilibrium coverage (full details in Appendix B2):

(19) ν eq , k = 1 - 1 - α α μ 0 X N X G - l k c k l ν l .

As the left-hand 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, N0 by inversion of Eq. (B28):

(20) N 0 = ν eq a 0 X ν .

Equations (19) and (20) therefore allow us to define an initial equilibrium state (Ni) 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.

3 Modelling results

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 1-month time step and successive mass classes that differ by a multiplicative constant ξ, so that mi=ξmi-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 PFT-specific 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 PFT-specific 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 grid-box mean values required to drive RED (using ESA land-cover 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 N0. This method successfully reproduces the ESA map of dominant PFT to good accuracy, as shown in Fig. 5 and Table 4.

Figure 5Maps of dominant PFT for (a) ESA LC_CCI dataset and (b) RED model equilibrium fractions. Sparse area is defined as where the total vegetation coverage is less than 10 %.

Table 4Goodness of fits for the RED equilibrium coverages to the coverages from ESA LC_CCI dataset across PFTs. r represents the Pearson correlation coefficient, after weighting by the grid-box area to account for latitudinal variation in grid-box areas.

Download Print Version | Download XLSX

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 N0, 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:

(21) γ = μ 0 g 0 m 0 ,

where g0 is given by Eq. (11) combined with Eqs. (B18) and (B20). Maps of γ values, derived in this way, are shown in Fig. 6.

Figure 6Diagnosed maps of mortality rates γ for each PFT, as required for consistency with the ESA observations and the UKESM growth rates. White areas correspond with zero coverage and/or zero growth.

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 m0 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 m0 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 BET-Tr of 0.232-0.007+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.

Table 5The area-weighted median values of observed coverage and driving net assimilate against μ0 and γ for the upper quartile of grid boxes for each PFT.

Download Print Version | Download XLSX

Site-level assessments of the rates of stand mortality within pantropical forests conclude a range of background rates (Lugo and Scatena1996; Phillips1996; 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 BET-Tr equilibrium mortality rates by looking at the values of γ in areas where we would expect to see old-growth forests. We use the top 25 % of coverages of the BET-Tr 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 land-use surveys (the mean mortality is defined here by weighting grid-box γ values by grid-box fractional coverages). There are a number of surveys relating stand mortality in regions prone to wildfires (Swaine1992; Kinnaird and O’Brien1998; Peterson and Reich2001; Van Nieuwstadt and Sheil2005; Prior et al.2009; Staver et al.2009; Brando et al.2014). In a broad sense, post-fire mortality rates can range from 0.06yr−1 to catastrophic rates around 0.8yr−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.

Figure 7Diagnosed mortality rates for (a) trees, (b) grasses, and (c) shrubs in the top quartile of coverage. Notches within the box represent the confidence bounds of the median. The confidence bounds are estimated using a bootstrap method. Bracketed numbers represent the number of grid points.


Another key issue is anthropogenic land use and land-use 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 (Esseen1994; Laurance et al.1998; Jönsson et al.2007). In order to maintain a near-constant agricultural fraction, regular disruption such as grazing is needed to prevent recolonization and secondary succession (Dorrough and Moxham2005; 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 (γ>0.075yr-1) that do correspond to areas where we see large fire activity (burn rate >0.1yr-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 (Coutinho1990; Medeiros and Miranda2008; 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).

Figure 8Comparison of observation-based estimates of tropical tree mortality (Phillips1996; Phillips et al.2004) to γ values diagnosed from RED for the BET-Tr PFT (for the top 25 % of fractions for this PFT). (a) Location of observational sites (blue and green crosses) versus the chosen RED grid points (red circles); (b) distribution of mortality across grid boxes; (c) mortality distribution across the BET-Tr grid points. Bracketed numbers in panel (b) represent the number of measurements and in panel (c) the number of grid points.

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.

Figure 9Comparison of diagnosed mortality rates, with observation-based maps of fire and land use. (a) Annual burnt area fraction from the ESA FIRE_CCI dataset; (b) crop fraction from the ESA LC_CCI 2000 dataset; (c) diagnosed mortality rate γ for the tree PFTs (BET-Tr, BET-Te, BDT, NET, NDT); (d) overlap of areas of higher tree mortality rates (γ>0.075yr-1) with areas of fire (burnt area >0.yr-1) and agriculture (crop fraction ≥25 %).

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 spin-up 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 long-term equilibrium state, as evidenced by the horizontal dashed lines in panels (a) and (b) of Fig. 10.

Figure 10Dynamical runs of RED for a grid box at the edge of the Amazonian rainforest, starting from bare soil (solid lines) and the diagnosed equilibrium state (dashed lines). (a) PFT fractions versus time; (b) biomass versus time; (c–e) snapshots of the number density distribution of the PFTs across mass classes at different times. Lines marked as + are the equilibrium runs, while X indicates the spin-up run. The ultimate steady state is determined by the balance between recruitment and mortality (Eq. 6). Intra- and inter-PFT occurs here through the shading of seedlings, which implies that just a fraction of the grid box (s, gap fraction) is available to grow seedlings (Eq. 7).


Faster-growing 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 m0. With the parameters used here, the vegetation fraction reaches close to its equilibrium value after about 20 years (panel a), but full spin-up 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 co-competing 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: spin-up 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 semi-analytical 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 semi-analytical 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 spin-up time hugely, especially for Earth system model (ESM) applications.

Figure 11Global model spin-up from bare soil. As for Fig. 10, solid lines are spin-up from bare soil; dashed lines are the equilibrium instillation run. Panel (a) represents the fractional global coverage relative to the total land area; panel (b) represents the total biomass of the vegetation.


4 Discussion

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. First-generation 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 tree-size distributions within forests (so called “demography”). A number of much more sophisticated second-generation 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 second-generation 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 tree-size distributions observed at a large-scale 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 (Cox2001). The modular design of RED allows for easy coupling to land-surface schemes, merely requiring the per unit grid-box 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 tree-size-dependent processes, although in this first study we chose to assume size-independent (but spatially varying) mortality rates for each PFT. Our previous work suggests that this is a good first-order 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 ORCHIDEE-MICT 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 ϕg=3/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 LM3-PPA 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 LM3-PPA 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 steady-state vegetation cover given information on the mortality and growth rates per unit area for each PFT. Such analytical steady-state solutions mean that RED can be easily initialized in drift-free 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 steady-state 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 steady-state solutions, RED is attractive for large-scale applications because it is both parameter sparse (“parsimonious”) and requires very few driving variables. The main driving variable is the time-varying 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 land-surface 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 PFT-specific 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 co-competing 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 zero-drift 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 individual-based models (Fischer et al.2016) and DGVMs specifically designed to capture sub-grid-scale 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 land-surface 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 forest-fire model (Burton et al.2019). These developments will allow us to simulate the effects of size-dependent tree mortality rates within the near future.

5 Conclusions

In this paper we have presented a new intermediate complexity second-generation 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 CO2 and climate. We have made the prototype RED code publicly available, and we hope that Earth system and land-surface modellers will make good use of this framework to further their own research.

Appendix A: Functional form of flux Fi in discretized RED

For large-scale 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 Fi as given in Eq. (5).

The total vegetation carbon in each mass class is Mi=miNi. The updated equation for Mi is therefore Eq. (4) multiplied by mi:

(A1) M i t + m i F i - F i - 1 = - γ M i .

The total carbon in the vegetation, M, is the sum of the carbon in each of the mass classes:

(A2) M = i M i .

Thus the update equation for the total carbon is

(A3) M t + i m i F i - F i - 1 = - γ M ,

which can be rewritten as

(A4) M t + i F i m i - m i + 1 = - γ M .

Now substituting Eq. (5) into Eq. (A4) gives

(A5) M t = i N i g i - γ M .

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.

Appendix B: Continuum solutions and demographic equilibrium theory

Equation (1) can be solved for the steady state if we assume metabolic scaling of growth using Eq. (2) and a size-independent mortality (Moore et al.2018, 2020) as follows:

(B1) n = n 0 m m 0 - ϕ g exp μ 0 ( 1 - ϕ g ) 1 - m m 0 1 - ϕ g , μ 0 = γ m 0 g 0 .

where n0 is a boundary condition that describes the number density at the mass m0. 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 m0. Similar analytical solutions can be derived for other measures of tree size, such as basal diameter or height (Moore et al.2018, 2020).

Figure B1The quasi-Weibull number density solution to DET (Eq. B1), assuming the same initial n0 and growth scaling ϕg=0.75 but different μ0 values.


Integrating Eq. (B1) from m0 to gives the total number density:

(B2) N eq = n 0 g 0 γ = n 0 m 0 μ 0 .

Other cohort integrals can be derived by integrating over the number density distribution, such as total growth rate (gndm):

(B3) G eq = g 0 N eq μ 0 1 - ϕ g ϕ g ϕ g - 1 exp μ 0 1 - ϕ g Γ 1 1 - ϕ g , μ 0 1 - ϕ g ;

total biomass (mndm):

(B4) M eq = m 0 N eq μ 0 1 - ϕ g 1 ϕ g - 1 exp μ 0 1 - ϕ g Γ 1 1 - ϕ g + 1 , μ 0 1 - ϕ g ;

and total vegetation cover (andm):

(B5) ν eq = a 0 N eq μ 0 1 - ϕ g ϕ a ϕ g - 1 exp μ 0 1 - ϕ g Γ ϕ a 1 - ϕ g + 1 , μ 0 1 - ϕ g ,

where Γ(a,b) is the incomplete upper gamma function. As we assume the allometric exponents presented in Niklas and Spatz (2004) (ϕg=3/4, ϕa=1/3), these functional forms simplify to

(B6) G eq = g 0 N eq 1 + 3 4 μ 0 + 3 8 μ 0 2 + 3 32 μ 0 3

(B7) M eq = m 0 N eq 1 + 1 μ 0 + 3 4 μ 0 2 + 3 8 μ 0 3 + 3 32 μ 0 4

(B8) ν eq = a 0 N eq 1 + 1 2 μ 0 + 1 8 μ 0 2 .

Finally, to convert a μ0 found using biomass (μ0,tdm) to one based on carbon mass, we use the formula

(B9) μ 0 = 2 1 - ϕ g μ 0 , t d m ,

assuming that biomass is twice the carbon mass.

B1 Closed continuous form

The lowest population flux, n0g0, is equal to the seedling boundary condition, F0, in Eq. (6):

(B10) n 0 g 0 = α 1 - α G m 0 s .

Substituting the total number density, Neq (Eq. B2), into the left-hand side and total growth, Geq (Eq. B6), into the righthand side yields a solution for the equilibrium coverage, assuming s=1-νeq, as follows:

(B11) γ N eq = α 1 - α g 0 m 0 N eq 1 - ν eq 1 + 3 4 μ 0 + 3 8 μ 0 2 + 3 32 μ 0 3 ,

which simplifies to

(B12) ν eq = 1 - 1 - α α μ 0 1 + 3 4 μ 0 + 3 8 μ 0 2 + 3 32 μ 0 3 .

Using Eq. (B8) we can write the total number density at equilibrium in terms of νeq:

(B13) N eq = ν eq a 0 1 1 + 1 2 μ 0 + 1 8 μ 0 2 .

This enables Eq. (B6) to be rewritten as

(B14) G eq = ν eq g 0 a 0 1 + 3 4 μ 0 + 3 8 μ 0 2 + 3 32 μ 0 3 1 + 1 2 μ 0 + 1 8 μ 0 2 .

This equation in turn defines the total assimilate:

(B15) P eq = 1 1 - α G eq .

Finally the total biomass can be written in closed form as

(B16) M eq = ν eq m 0 a 0 1 + 1 μ 0 + 3 4 μ 0 2 + 3 8 μ 0 3 + 3 32 μ 0 4 1 + 1 2 μ 0 + 1 8 μ 0 2 .

B2 Discrete steady state

To solve for the discrete model equilibrium, we start from the flow equation from Eq. (4) with the term N/t0:

(B17) γ N i + F i = F i - 1 .

Considering the population flux (Eq. 5), we find Ni in relation to the lower mass class, Ni−1:

(B18) N i = N i - 1 g i - 1 / ( m i - m i - 1 ) g i / ( m i + 1 - m i ) + γ = N i - 1 λ i .

Assuming no population grows out of the top class, λI is given as

(B19) λ I = g i - 1 ( m i - m i - 1 ) γ .

λi can be simplified to depend only on μ0, by using μ0=(γm0/g0) (Eq. 14) and applying the mass scaling of growth rates gi=g0(mi/m0)gϕ. We can show that λi and λI are

(B20) λ i = m i - 1 / m 0 ϕ g m 0 / ( m i - m i - 1 ) m i / m 0 ϕ g m 0 / ( m i + 1 - m i ) + μ 0 , λ I = ( m i - 1 / m i ) ϕ g m 0 ( m i - m i - 1 ) μ 0 .

An expression for the total stand density at equilibrium, Neq, can be derived. Using Eq. (B18), we can represent any population of mass class i in terms of the lowest mass class N0 as follows:

(B21) N i = N 0 j = 1 i λ j .

Therefore, when finding the total number of stands relative to N0 we get

(B22) N eq = N 0 1 + i = 1 I j = 1 i λ j = N 0 X N ,

where XN describes the sum of the all mass classes as a proportion of N0. We can describe the total class growth rate in relation to N0 as

(B23) G i = N 0 g i j = 1 i λ i .

By using the allometric relationship (Eq. 2)

(B24) G i = N 0 g 0 m i m 0 ϕ g j = 1 i λ j ,

we describe the total class growth rate in relation to the lowest class growth rate, N0g0. Like Neq, we can show the total growths across all classes is therefore

(B25) G eq = N 0 g 0 1 + i = 1 I m i m 0 ϕ g j = 1 i λ j = N 0 g 0 X G .

We can repeat the same process for coverage

(B26) ν i = N 0 a i j = 1 i λ j

and using allometric relationship (Eq. 3)

(B27) ν i = N 0 a 0 m i m 0 ϕ a j = 1 i λ j .

This gives the total coverage, νeq, as

(B28) ν eq = N 0 a 0 1 + i = 1 I m i m 0 ϕ a j = 1 i λ j = N 0 a 0 X ν .

Finally, for the total carbon mass within the class, we get

(B29) M i = N 0 m i j = 1 i λ i ,

with the total carbon density equalling

(B30) M eq = N 0 m 0 1 + i = 1 I m i m 0 j = 1 i λ j = N 0 m 0 X M .

In equilibrium, the rate of the recruitment of seedlings (Eq. 6) must balance the rate of loss of plants due to total mortality (γsNeq):

(B31) γ N eq = α ( 1 - α ) G eq m 0 s .

Substituting in Eq. (B22), Eq. (B25) yields a balance equation for the kth PFT:

(B32) α 1 - α 1 - l c k l ν l = μ 0 X N X G .

We can get the equilibrium fraction of a PFT, k, by rearranging the above equation, assuming ckk=1 as follows:

(B33) ν eq , k = 1 - 1 - α α μ 0 X N X G - l k c k l ν l .

Once the value of μ0 has been derived in this manner, we can find N0 by inversion of Eq. (B28) as follows:

(B34) N 0 = ν eq a 0 X ν .

Substituting Eq. (B33) into Eq. (B34) allows us to determine N0 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 mi and a geometric mass class scaling, mi+1=ξmi. 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.

Figure B2Comparison of the discretized model to the continuum analytical solution, showing convergence for higher numbers of mass classes. This example uses parameters for broadleaf evergreen tropical trees (BET-Tr PFT) with α=0.1: (a) equilibrium coverage νeq versus μ0 for the exact continuum solution (black line) and discretizations of the mass dimension with varying numbers of mass classes and mass class width scaling (ξ); (b) absolute error in the modelled value of νeq against the number of mass classes using the optimum value of ξ for each case; (c) optimum ξ versus number of mass classes, with contours showing the absolute error in νeq. Panels (b) and (c) assume μ0=0.25. The white dots in (c) have the same number of classes and scaling as the discrete lines in (a).


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 trade-off. 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 m0, by writing mi=m0(ξ)i. Therefore, by using mi+1-mi=m0(ξ)i(ξ-1), we find that our equilibrium form of λi is reduced to

(B35) λ i = ξ ( ϕ g - 1 ) ( i - 1 ) ξ i ( ϕ g - 1 ) + μ 0 ( ξ - 1 ) , λ I = ξ ( ϕ g - 1 ) ( i - 1 ) μ 0 ( ξ - 1 ) .

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

(B36) ξ ν eq , continuous - ν eq = - ξ ν eq ,

where νeq corresponds with the discrete equilibrium, i.e. Eq. (B32), with νeq=(1-s). Setting Eq. (B36) equal to zero we reduce the relationship to only a dependence on XN and XG:

(B37) 0 = ξ X N X G = X G X N - X G X N .

Finding the partial derivative of XN, using the geometric form of Eq. (B18), we get

(B38) X N = j = 1 I i = 1 j λ i i = 1 j λ i λ i ,

and for XG we get the following:

(B39) X G = j = 1 I ξ j ϕ g i = 1 j λ i j ϕ g ξ - 1 + i = 1 j λ i λ i .

Finding λi we get

(B40) λ i = λ i ( 1 - i ) ( ϕ g - 1 ) ξ - 1 - λ i i ( ϕ g - 1 ) ξ ϕ g - 2 + μ 0 ξ ( i - 1 ) ( 1 - ϕ g ) ,

and for the top class, λI we get the following:

(B41) λ I = ( 1 - ξ - 1 ) ( I - 1 ) ( ϕ g - 1 ) - 1 ξ - 1 λ I .

To numerically solve for the minimum, we must differentiate Eq. (B37), with respect to ξ. Through the product rule we get

(B42) 2 ξ 2 X N X G = X G X N ′′ - X G ′′ X N .

Differentiating equation (Eq. B38) and simplifying gives

(B43) X N ′′ = j = 1 I i = 1 j λ i i = 1 j λ i ′′ λ i ,

and doing the same for Eq. (B39) gives

(B44) X G ′′ = j = 1 I ξ j ϕ g i = 1 j λ i j ϕ g ξ - 2 ( j ϕ g - 1 ) + i = 1 j 2 j ϕ g ξ - 1 λ i - λ i ′′ λ i .

λi′′ is given by

(B45) λ i ′′ = λ i - λ i λ i ( i - 1 ) ( ϕ g - 1 ) ξ - 1 - ( i - 1 ) ( ϕ g - 1 ) ξ - 2 - λ i ( ϕ g - 1 ) ξ - 1 ξ ( i - 1 ) ( 1 - ϕ g ) i ( ϕ g - 1 ) ξ ϕ g - 2 - μ 0 ( i - 1 ) ξ ( i - 1 ) ( 1 - ϕ g ) λ i λ i 2 .

For the double differential of λi we get

(B46) λ i ′′ = λ ′′ i 2 λ i + λ i ξ - 1 × ( I - 1 ) ( ϕ - 1 ) ξ 2 - λ i λ i .

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.

Appendix C: Sensitivity of diagnosed mortality rates to model parameters

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:

(C1) γ = α P eq a 0 m 0 1 - ν eq ν eq 1 + 1 2 μ 0 + 1 8 μ 0 2 .

The key external inputs to this equation are the observed PFT fraction νeq and the net assimilate Peq. In addition, our estimates of γ are dependent on the internal model parameters, α and m0.

Figure C1The sensitivity of the mortality rate to assumed input variables: (a) coverage, νeq; (b) Peq carbon assimilate rate; and model parameters (c) reseed fraction, α, and (d) boundary mass, m0. The solid black line indicates the fixed values with corresponding (b–d) ±20 % or (a) ±5 % variation (dotted black lines).


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 m0 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 %.

Code availability

The RED model Python code is archived at (Argles2019). 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 (, last access: 2 September 2020).

Author contributions

Originally the model framework was in JRM's thesis (Moore2016) 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.

Competing interests

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.

Financial support

This research has been supported by the Newton Fund (CSSP-Brazil), 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).

Review statement

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,, 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 heat-induced 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,, 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.: Carbon-concentration and carbon-climate feedbacks in CMIP6 models, and their comparison to CMIP5 models, Biogeosciences Discuss.,, 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 Soares-Filho, 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 land-use and land-cover changes on climate and land carbon storage in CMIP5 projections for the twenty-first century, J. Climate, 26, 6859–6881,, 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, land-use change and vegetation dynamics in the Joint UK Land Environment Simulator vn4.9 (JULES), Geosci. Model Dev., 12, 179–193,, 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,, 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 size-density distributions in natural forests from following scaling relationships, Ecol. Lett., 6, 980–989,, 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 landscape-scale 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,, 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 old-growth 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 second-generation 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,, 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., Muller-Landau, 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,, 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,, 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,, 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,, 2016. a

Harper, A. B., Wiltshire, A. J., Cox, P. M., Friedlingstein, P., Jones, C. D., Mercado, L. M., Sitch, S., Williams, K., and Duran-Rojas, 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,, 2018. a

Haverd, V., Smith, B., Nieradzik, L. P., and Briggs, P. R.: A stand-alone tree demography and landscape structure module for Earth system models: integration with inventory data from temperate and boreal forests, Biogeosciences, 11, 4039–4055,, 2014. a, b

Higuchi, P., Oliveira-Filho, A. T., Bebber, D. P., Brown, N. D., Silva, A. C., and Machado, E. L.: Spatio-temporal patterns of tree community dynamics in a tropical forest fragment in South-east 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 non-structural carbohydrates on the simulated response of tropical forests to drought, Biogeosciences, 17, 3589–3612,, 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., Rankin-de 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., Alvarez-Dá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., Silva-Espejo, 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,, 2019. a

Lima, R. A., Muller-Landau, 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,, 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,, 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.: Post-fire resprouting and mortality in cerrado woody plant species over a three-year 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,, 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,[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,, 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 tree-size distributions and biomass density in Amazonia, Biogeosciences, 17, 1013–1032,, 2020. a, b, c, d, e, f, g, h, i, j

Muller-Landau, 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,, 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 near-term 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,, 2004. a, b, c, d, e

Pavlick, R., Drewry, D. T., Bohn, K., Reu, B., and Kleidon, A.: The Jena Diversity-Dynamic Global Vegetation Model (JeDi-DGVM): a diverse approach to representing terrestrial biogeography and biogeochemistry based on plant functional trade-offs, Biogeosciences, 10, 4137–4177,, 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.: Long-term environmental change in tropical forests: increasing tree turnover, Environ. Conserv., 23, 235–248,, 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,, 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,, 2015. a

Prior, L. D., Murphy, B. P., and Russell-Smith, 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,, 2020. a

Sato, H., Itoh, A., and Kohyama, T.: SEIB-DGVM: A new Dynamic Global Vegetation Model using a spatially explicit individual-based approach, Ecol. Model., 200, 279–307,, 2007. a

Scheiter, S., Langan, L., and Higgins, S. I.: Next-generation dynamic global vegetation models: Learning from community ecology, New Phytologist, 198, 957–969,, 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, 4513-4558, 2019. a, b

Smith, B.: LPJ-GUESS-an 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 height-structured competition, Biogeosciences, 12, 2655–2694,, 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 ORCHIDEE-MICT v8.4.2, Geosci. Model Dev., 11, 409–428,, 2018. a, b

Short summary
The Robust Ecosystem Demography (RED) model simulates cohorts of vegetation through mass classes. RED establishes a framework for representing demographic changes through competition, growth, and mortality across the size distribution of a forest. The steady state of the model can be solved analytically, enabling initialization. When driven by mean growth rates from a land-surface model, RED is able to fit the observed global vegetation map, giving a map of implicit mortality rates.