Robust Ecosystem Demography (RED): a parsimonious approach to modelling vegetation dynamics in Earth System Models

. A signiﬁcant 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 cohort-based patch models are now becoming established in the latest DGVMs. These 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 5 impact of stochastic disturbance events on the forest structure and biomass, but at the cost of needing to update a probability density function in two-dimensions. Here we present the Robust Ecosystem Demography (RED) , in which the pdf is collapsed on to 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 soluble. The population of each Plant Functional Type (PFT) is partitioned into mass classes with a ﬁxed baseline mortality along with an assumed power-law scaling of growth-rate 10 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 15 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.

2 Description of the Model

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(m) of plants per unit area of mass m: ∂n ∂t + ∂ ∂m n g = −γ n.
(1) 5 Here g(m) is the growth-rate and γ(m) is the mortality rate of a plant of mass m. In general, g and γ could take any form, but we make simplifying assumptions for these functions consistent with observed n(m) from forest inventory data (Moore et al., 2018(Moore et al., , 2019. Simply, we assume that γ is independent of plant mass, and that g(m) is a power-law of plant mass: 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 10 ( Niklas and Spatz, 2004). We also follow (Niklas and Spatz, 2004) in assuming the scaling of plant canopy area a with plant mass: where φ a = 0.5 by default.
Solutions for n can be integrated over mass to derive the total plant number, N = n dm, the total growth-rate, G = 15 g n dm, the total biomass, M = m n dm, and the fractional area covered ν = a n dm.

Discrete Mass Classes
We wish to produce a model of vegetation demography that can be updated numerically and which explicitly conserves vegetation carbon. In order to do this we integrate Eq. (1) over finite mass ranges: 20 where i denotes the i th mass class; F i is the flux of plants growing out of the i th 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 i th mass class; and N i is the number of plants per unit area in the i th 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 Section 2.4 and 2.5. In order to conserve carbon (see below) the flux F i must take the form: where m i is the mean mass of a plant in the i th mass class, and g i is the growth-rate per plant of the i th mass class [kgC yr −1 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). The total vegetation carbon in each mass class is M i = m i N i . The update 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: 10 Now substituting Eq. (5) into Eq. (9) 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.

Seedling production and gap competition 15
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 being allocated to the growth of existing plants. In addition, we assume that only those seedlings growing in 'gaps' will survive. The net incoming flux of seedlings of mass m 0 is therefore: 20 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 Figure 1.
Equation (11) assumes a random overlap between seedlings and the existing vegetation. This lower boundary condition is the only place within RED where there is significant competition. Minimum overlap, which is broadly consistent with 'perfect plasticity' (Strigul et al., 2008), is assumed once seedlings have been injected into the cohort according to Eq. (11). The space available to the seedlings of the k th PFT is calculated from the area fractions of the PFTs to which it is subdominant: where ν l is the area fraction of the l th 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 domi-5 nant over PFT l, c kl = 0 ( Figure 1). This is similar to the competition regime in the TRIFFID model (Cox, 2001), and allows for the co-existence between inter-functional groups of PFTs. 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 1.  Figure 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 ): 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: These γ d values can in principle be both PFT-dependent and mass-dependent (e.g. to capture forestry practices).  Figure 2. Schematic of RED coupled to an ESM on 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 on to an assumed baseline mortality, γ b , as a function of both PFT and mass-class.

10
The input values of net assimilate for each PFT (P ), define the total structural growth-rate, G = (1 − α) P , and the seedling flux F 0 (via Eq. (11)), using PFT-specific values of the parameter α (see table 2). The definition of the total structural growthrate at a given time-step: 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 mass class (N i ) is updated using a discretised form of Eq. (4): where ∆t is the RED timestep (typically 1 month), and the superscript (j) denotes the j th timestep. The lower boundary seedling flux is calculated from Eq. (11) using Eq. (12). We do not 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 10 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): 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 15 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 20 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. (13)), to calculate the total litterfall flux into the soil/litter system.

Steady-State
The steady state of the continuum model defined by Eq. (1) and Eq.
(2) can be solved analytical for each PFT (Moore et al., 2018(Moore et al., , 2019. The continuum analytical solutions for the equilibrium mass distribution (n eq (m)), the total plant number (N eq ), 25 biomass (M eq ), growth-rate (G eq ) and fractional area (ν eq ) are summarised 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 : In order to initialise the numerical RED model in a drift-free initial state, we also derive the steady-state of the discrete model (Equation (17)), which will differ from the continuum model for a finite number of mass classes. The equilibrium solution of Eq. (17) defines a recursive relationship for the number of plants N i in each mass-class: where Thus for the discrete model the shape of the mass-number distribution also depends on the mortality-to-growth parameter, 10 Where X N : and I is the top mass-class. Where i and j are class indices over the sum product. Likewise, we can calculate: 1. the total structural growth at equilibrium:

15
where Eq. (2) implies: 2. the total biomass at equilibrium: where:  . We can use an observation to constrain (red dashed arrow) µ0 giving us an equilibrium coverage. If the addition of total growth is known (black double-headed arrow), we can rearrange µ0 to find an equilibrium mortality rate (or vice-versa).
3. the fractional area covered by the PFT at equilibrium: where Eq. (3) implies: The equations above therefore define the equilibrium state of the discrete system for given values of N 0 and µ 0 . The value 5 of µ 0 can be estimated from forest demographic data where this is available (Moore et al., 2018(Moore et al., , 2019. However, for global applications we rarely have more observations than the fractional coverage of each PFT. Under these circumstances, we use the condition that at equilibrium the rate of injection of seedlings (Equation (11)) must balance the rate of loss of plants due to mortality (γ N eq ): 10 Now substituting in Eq. (22), Eq. (24) and Eq. (12) yields a balance equation for the k th PFT: where X N and X G are given by Eq. (23) and Eq. (25) respectively. As the righthand-side of this equation depends only on prescribed constants and µ 0 , Eq. (31) can be inverted (by numerical iteration) to estimate µ 0 for observed values of the PFT fractions (ν l ) and an assumed value of α (see Table 2). Once the value of µ 0 has been derived in this manner, it can be used to 15 calculate X ν via Eq. (29), and therefore N 0 by inversion of Eq. (28): Equations (31) and (32) therefore allow us to define an initial equilibrium state (N i ) which is consistent with observed area fractions of each PFT ( Figure 3). 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 as part of the global tests of RED described in the next section.

5
For these tests, the numerical RED model is set-up to use the 9 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., in prep). RED is integrated forward using a one month time-step and successive mass-classes that differ by a multiplicative constant ξ, so that m i = ξ m i−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, 10 assuming µ 0 = 0.25 (see Appendix B2). Other PFT-specific parameters are assumed as summarised in Table 2. Table 2. The PFT list and their corresponding parameters (m0, a0, h0), seedling fraction (α) and competition coefficient (c pft,j ). The growth allometry of trees and shrubs across size is assumed to follow (Niklas and Spatz, 2004) The competition coefficients given describe which PFT functional group shades the current PFT, if c pft,j = 1, the PFT is shaded, otherwise it is not (Table 1).

Local: Simulating Succession
We begin by demonstrating the vegetation succession simulated by RED in an idealised spin-up from bare-soil (  The modelled evolution of the plant number versus mass distribution for each PFT is shown in panel c (after 6 years), panel d (after 13 years) and panel (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 5 to reach equilibrium.
The dashed-lines in Figure 4 represent a dynamical RED simulation from the diagnosed demographic equilibrium state.
This state is derived using the methodology described in Section 2.5, with one significant change. The competition rules given by Eq. (12) and Table 1 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 10 class, we choose to initialise the dominant PFT to have the total area fraction of all the PFTs in that vegetation class.

Global: Simulation of Current Vegetation Map
In this section we use a similar procedure to diagnose the map of PFT-specific mortality rates consistent with the current observed vegetation state, and rates of net assimilation (P ) provided by UKESM (Sellar et al., in prep Figure 6. Driving assimilate (Equation (13)) from UKESM, averaged between 2000-2010. The average is constructed by truncating any negative growth rates to zero.
We use the procedure outlined in Section 2.5 to estimate spatially-varying values of µ 0 for each PFT, using Eq. (31), and then Eq. (32) to estimate N 0 . This method successfully reproduces the ESA map of dominant PFT to good accuracy, as shown in Figure 7 and Table 3.
14 https://doi.org/10.5194/gmd-2019-300 Preprint.   There are a few reasons why the model equilibrium fractions are not an identical fit with the observed fractions. Firstly, the two datasets are not necessarily consistent -there are a few places (Central Asia, Sahel) were the average UKESM assimilates used is zero, not aligning with the positive coverage from the ESA dataset. Secondly, areas of mixed PFTs within the same vegetation class, as previously stated, will have an adjusted RED equilibrium fraction where the dominant PFT equilibrium will be the sum of the vegetation class. As stated within the previous section, Figure 8, demonstrates that the initialisation scales globally. There are two RED simulations, one being started from bare soil the other from the diagnosed equilibrium. Using a constant assimilate rate ( Figure   6) and the mortality distribution (Figure 10), we see convergence towards the two runs in global coverage (Panel a) and global biomass (Panel b). The equilibrium initialised run also demonstrates the practical advantage RED posses in avoiding spinning up. In instances that have a low spin-up growth rates, this saves significant time.

Global: Diagnosed Plant Mortality Rates
As 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. (19) : where g 0 is given by Eq. (16) Figure 9. Diagnosed 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 -mortality is assumed infinite within RED.
The mortality rate derived is very dependent on the overall coverage and the total assimilate. Having a high coverage with a low growth rate will result in RED compensating through having a low mortality rate (and vice-versa). This explains why some mortality rates of PFTs seem to posses large variations ( Figure 10). Furthermore, the choice of α (Equation (16)) and the m 0 is also influential when it comes to the value of γ. Under the assumption that high coverages are close to a 'healthy' 10 environment for a PFT. We can take a sub-sample of the grid-boxes that are within the top quartile of non-zero coverages (ν eq > 0.01) ( Table 4). The median µ 0 value diagnosed from the top quartile of BET-Tr of 0.232 +0.008 −0.007 (Table 3), is close to the values calculated by our previous paper (Moore et al., 2019) of approximately 0.235 for all of South America. The value within the paper 0.198 is converted to dry carbon mass through Eq. (B.9). Carbon dry mass constitutes approximate half of the total dry mass (Thomas and Martin, 2012). In addition to parametrisations (such as α), some of the differences between these µ 0 values can arise from the discretization of the model, as the discretised form will underestimate the diagnosed µ 0 values to meet the same observation -when compared with the continuous form ( Figure B1.a). A possible usage of RED might be to diagnose a µ 0 , fix a baseline mortality rate from surveys and diagnose the required carbon assimilate to match an additional observation (Figure 3). Potentially providing a future constraint on ESM growth rates for PFTs. There have been multiple site-level assessments of the rates of stand mortality within pan-tropical forests -typically the background rate is between 1 yr −1 to 4 yr −1 (Lugo and Scatena, 1996;Phillips, 1996;Phillips et al., 2004). (Phillips, 1996) estimates 60 mortality rates collected across 40 pan-tropical sites for tree sizes greater than 10 − 25 cm dbh, Later work by Phillips in 2004 used the demographic data from the RAINFOR dataset of trees ≥ 10cm 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 10 to see old growth forests. Within the top 25% of coverages, we assume represent areas of undisturbed forest. BET-Tr captures the baseline mortality rates seen post 2000 in the Amazon from site data ( Figure 10). ) (c) RED Figure 11. Comparision of observed tropical tree mortality with γ. Comparing datasets presented in (Phillips, 1996;Phillips et al., 2004) for an adjusted estimate of observed stand mortality with the equilibrium mortality rates for BET-Tr within the largest 25% of fractions. Next we investigate if the equilibrium mortality rates implicitly capture areas of disturbances. We compare the mean woody equilibrium mortality rate to fire and land-use surveys. The woody mortality is defined as the sum of γ weighted against the coverage of Trees and Shrubs to their collective coverage. Areas with large rates of disturbances area generally not expected to conform to the equilibrium assumptions, such as DET, used to initialise RED. Generally, biomes, that are humid and dry such as savannah or grassland, wildfires play a natural part in maintaining the balance of vegetation (Bond et al., 2005), therefore 5 for woody PFTs we expect to see a raised mortality/higher µ 0 . Using the ESA FIRE_CCI dataset (Chuvieco et al., 2019) we can estimate the rate of burnt fraction per year: where Burnt Area and Burnable Fraction are given from the dataset. Area is inferred from the longitude, latitude quadrant.
The resolution of FIRE_CCI is 0.25 • × 0.25 • is bi-linearly interpolated onto the simulated grid-boxes. Taking the averages of 10 the burnt fraction rate between the months between 2000 and 2010, and converting into an annual burn rate.
We also carry-out a comparison with agriculture, we expect that in area of land-use we will see raised mortality for woody PFTs as the industry will use the space for crops and pasture. There is difficulty in getting an explicit rate of clearance from land-use, however, by comparing with the fraction of cropland we achieve a non-direct geographic comparison to the to the dataset (Figure 12). For the crop fraction we use the 2000 ESA LC_CCI inferred PFT from (Li et al., 2019) of half a degree resolution -again interpolating onto the RED grids.  If the carbon assimilate is similar to that of forests, to compensate for the lower observed coverage of woody PFTs the trees will posses higher mortality rates. Such compensation will be a non-direct approach to estimating mortality of disturbance prone areas. We see that in comparison to the distribution of mortalities increases in areas of observed land-use and fire ( Figure 13), however it does underestimate the annual burn rate. 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 complexity. In principle, the number of age-defined patches grows indefinitely, and this can only be managed 5 by arbitrarily merging patches of different ages after a certain age.
In this paper, we represent an intermediate complexity second generation DGVM ('RED'), which is designed to capture important features of plant demography, and yet avoid unwieldy computation. Our guiding principles in the development of RED have been that the model should: (i) simulate forest tree-size distributions; (ii) be globally applicable for ESM applications ; (iii) be parameter sparse to minimise parameter uncertainties; (iv) be analytically soluble for steady-states to aid model initiali-10 sation. These design criteria, along with evaluation exercises against observed tree-size distributions in North America (Moore et al., 2018) and South America (Moore et al., 2019), has led us to make a number of simplifications.
Firstly, the number density for each PFT is treated as function of plant mass alone. This immediately eliminates the need to explicitly represent patches, and therefore removes age as an independent dimension. Secondly, we assume that plant growthrates vary as a power of plant mass. By default we assume a power of φ g = 3/4, which is consistent with Metabolic Scaling

15
Theory (Enquist et al., 1998) and the empirically determined allometric relationships of (Niklas and Spatz, 2004). Finally, we assume that light-competition is only significant for the lowest 'seedling' mass class. This enables us to capture the impacts of light competition on seedling emergence through a simple 'gap' boundary condition.
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 initialised 20 in drift-free pre-industrial 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.
Aside from the existence of analytical stead-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 25 plant growth-rate for each PFT, which is defined as net primary production minus the local litter-fall. 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 PFT-specific constant that is independent of mass, but which could be dependent on plant mass and time to represent individual disturbance events (e.g. forest fires, disease outbreaks). Despite its simplicity, the RED modelling is able to fit the global distribution of vegetation types (Figure 7), and simulates realistic suc-30 cession including changes in forest demography (Figure 4).
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 sub-dominant 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 Section 3.1. Such competitive exclusion is a common problem in DGVMs (Fisher et al., 2018). Adapting the 'gap' boundary condition (Equation (12)) appears to be a promising way to deal with this issue in RED, without unduly increasing model complexity.
We see this as a key priority for future research.
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, mean canopy height, and the 5 ratio of mortality-to-growth. This in turn allows RED to be calibrated using observations of any two of these quantities (Figure 3). The analytical solutions also allow optimality hypotheses to be explored (e.g. the hypothesis that the fraction of net assimilate allocated to seed production maximises stand-density and/or biomass)

Conclusions
In this  to estimate how plant growth-rate varies with plant mass, and perfect crown plasticity to minimise light competition) and also comparison to observed forest demography (Moore et al., 2018(Moore et al., , 2019. 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 15 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 publically available, and we hope that Earth System and land-surface modellers will make good use of this framework to further their own research.  Competition coefficient, the fraction a PFT, k, is shaded by the canopy of PFT l. − Equilibrium µ0 The boundary turnover parameter -the ratio of mass lost to growth gained in the boundary mass class. − λi The proportional population of the i th class to the i th − 1 class at equilibrium. − eq Subscript denoting a variable in equilibrium. The current time-step. − ξ The size scaling coefficient, mass classes are defined as mj = ξ mj−1, with ξ > 1. − with growth (height, basal diameter, etc.). The variable µ 0 can be thought of a parameter tied to the rate of biomass lost to biomass gained. The larger µ 0 is the greater the associated cost of replacing lost biomass -the smaller the total population density. Where n 0 is a boundary condition that describes the number density at the mass m 0 .
Integrating Eq. (B.1) from m 0 to ∞ gives us estimates for the total number density: We can also gain estimates of the total growth and biomass values by integrating with the allometric relationships: the total biomass: (B.4) and the total vegetation cover: Where Γ(a, b) is the incomplete upper gamma function. When we assume the constants presented in (Niklas and Spatz, 2004) φ g = 3 4 , φ a = 1 3 simplifies towards: G eq = g 0 N eq 1 + 3 4µ 0 + 3 8µ ν eq = a 0 N eq 1 + For converting a µ 0 found using total dry mass (µ 0,tdm ) of 1 kg to that of dry carbon mass: Using Eq. (B.2) and Eq. (B.6) with the competitive constraint, we find that the equilibrium fraction is given by: We can rearrange Eq. (B.10) into Eq. (B.5) allowing for the substitution of ν eq into the equilibrium solutions (Equation (B.6) and Eq.(B.7)). For instance, the exact solution for the total number density is given as: enables the calculation of total growth to Eq.(B.6)) as: There are a few differences between the numerical steady state and the continuous form of RED. Firstly, the truncation point of the top mass class results in a underestimation of the total coverage, biomass and number density for an identical µ 0 . The second source of difference arises from the binning of the mass classes. Discretising results in the continuous scaling of the growth and mass is not fully captured. In the current scheme a number density between the masses m j of m j+1 will have its physiological characteristics represented at the m j mass, this can lead to underestimating the total growth/biomass/coverage 20 within the class. This is demonstrated in figure B1.a, where the total coverage of the biomass is lower in the discrete model than the continuous solution. There is also convergence when the bins between the methods when the number of classes goes towards infinity (I → ∞) and class widths goes towards zero (ξ → 1).
From Figure B1.c, we see that there is a clear optimum amount of scaling for a given number of classes which minimise 25 the difference between continuous and discrete. This can be found by taking the difference of the continuous and discrete coverages and differentiating with respect to ξ to find the minima. It should be noted that the as the continuous form is not dependent on ξ, we get:

∂ ∂ξ
[ν eq,continuous − ν eq ] = − ∂ ∂ξ [ν eq ] . (B.14) Where ν eq corresponds with discrete equilibrium (Equation (31), with ν eq = (1 − s)). Setting Eq. (B.14) equal to zero we reduce the relationship to just a dependence of X N and X G : Finding the partial derivative of X N , we get: (B.16) and for X G : Finding λ i we get: − λ i i(φ g − 1)ξ φg−2 + µ 0 ξ (i−1)(1−φg) , (B.18) and for λ i : To numerically solve for the minimum, we must differentiate Eq. (B.15), with respect to ξ. Through the product rule we get: Eq. (B.16) differentiated simplifies towards: and Eq. (B.17): λ i is given as: For the double differential of λ i we get: We now possess the identities needed to perform a numerical root finding algorithm for the optimum bin scaling for a given class. Using a Newton root finding method for Eq. (B.15) with it's differential; Eq. (B.20), we find the optimum. On figure B1 the optimum line is shown as the bright dashed black line.
Author contributions. Originally the model framework 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 and PMC. Currently

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