Equilibrium forest demography explains the distribution of tree sizes across North America

The distribution of tree sizes within a forest strongly influences how it will respond to disturbances and environmental changes such as future climate change and increases in atmospheric CO2. This means that global vegetation models must include variation in tree size to accurately represent carbon sinks, such as that seen in North America. Here we use an analytical model of large-scale forest demography which assumes tree growth varies as a power of tree diameter whilst tree mortality is independent of size. The equilibrium solutions of this model are able to accurately reproduce the tree-size distributions, for 61 species and four plant functional types, measured across North America, using just a single species-specific fitting parameter, μ, which determines the ratio of mortality to growth. The predictions of metabolic scaling theory for tree-size distributions are also tested and found to deviate significantly from observations and that maybe explained by the assumptions made about how individual trees fill the available space. We show that equilibrium forest demography implies a single curve that relates mean tree diameter to μ, and that this can be used to make reasonable estimates of the whole dataset mean trunk diameter by fitting only to the larger trees. Our analysis suggests that analytical solutions such as those in this paper may have a role in aiding the understanding and development of next-generation Dynamic Global Vegetation Models based on ecosystem demography.


Introduction
Understanding the size-dependent dynamics of forests is of crucial importance for being able to predict the future role of vegetation in the carbon cycle and hence climate change. The land biosphere currently performs an important role for humanity by absorbing about a quarter of anthropogenic CO 2 emissions (Ciais et al 2013). However, the fluxes of CO 2 between land and atmosphere are known to be sensitive to climate, leading to the possibility of significant climate-land carbon cycle feedbacks (Cox et al 2000). This motivated the inclusion of climate-carbon cycle feedbacks in many of the climate projections reported in the most recent Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR5) (Arora et al 2013).
Unfortunately, the projections from such Earth System Models (ESMs) differ markedly (Friedlingstein et al 2014). While the evolution of the global ocean carbon sink and its large-scale spatial pattern are similar amongst the models, the future land carbon sink has a huge-range of uncertainty (as much as 500 GtC by 2100 for 1% increase in CO 2 emissions per year). The divergences amongst model projections are even more significant when feasible changes in land use are included (Brovkin et al 2013). Such large uncertainties feed through into estimates of the emission reductions required to stabilise to a given level of global warming (Friedlingstein et al 2014). There is therefore an urgent need to improve the representation of vegetation and soil processes in ESMs.
To this end, dynamic global vegetation models (DGVMs) are increasingly used within climate models to capture biophysical and carbon cycle feedbacks associated with changes in vegetation distribution and structure (Sitch et al 2008). In the IPCC AR5 such Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. models also attempted to include the impacts of changes in land use on climate (Brovkin et al 2013). However, most first generation DGVMs currently fail to represent forest size distributions (Sitch et al 2015, Fisher et al 2018 and so are unable to realistically simulate sinks associated with forest regrowth, such as that seen in the US (Hurtt et al 2002), and recovery from disturbance (Zhu et al 2018). Fisher et al (2018) discusses the development of a new generation of DGVMs, which are now representing size either via individual based models (Shugart et al 2018) or using cohort-based ecosystem demography models (Moorcroft et al 2001).
DGVMs are by their nature increasingly complex numerical models, and it is often time-consuming to study the effects of changing parameters. They are also dependent of initial conditions and have to trade off numerical accuracy (rounding and truncation errors) against speed and simplicity.
Conversely analytical solutions often have to be simpler and only capture some essential features but apply to all parameter values and therefore provide a more complete view of the behaviour of those key features. Equilibrium solutions are especially parameter sparse, as they are typically independent of initial conditions (such as disturbance history). So analytical equilibrium solutions are especially generic. They can also be used to aid model initialisation, by allowing the model to start at an equilibrium state rather than having to 'spin-up' and they can be used to validate numerical schemes or to aid in checking a model against data.
Here we present an analytical model of equilibrium forest demography called demographic equilibrium theory (DET) and demonstrate that it fits the distributions of tree-sizes measured across North America. We also compare our model against metabolic scaling theory (MST) (West et al 2009), which is a theory that derives predictions about forest structure and dynamics based on principles of metabolism and allometry of the individual and also of how individuals fill the available space. A key prediction of the spacefilling assumption of this theory is that the size distribution must scale with the trunk diameter D as D −2 . We test this by choosing to use the allometry from MST in DET but not the space-filling assumptions.

Demographic equilibrium theory
The distribution of tree sizes in a forest is determined by demographic processes such as growth, mortality and recruitment. Growth is simply the extra tissue acquired via photosynthesis but limited by nutrient, water and light availability. Mortality is a single demographic process arising from a number of biotic and abiotic processes, such as competition, senescence and disturbance (including land use change). Recruitment is the process of new individuals joining the population due the reproductive process.
These size-structured forest dynamics can be thought of as analogous to the concept of continuity in fluid dynamics (Kohyama et al 2003) but with a loss term to account for mortality. In forests the number of trees in a given size range is determined by the rate at which individual trees join it via growth and recruitment and leave it through mortality and growth (van Sickle 1977, Coomes et al 2003, Kohyama et al 2003. The demographic process can be expressed dynamically for size-structured populations by a onedimensional drift equation or 'continuity equation' of fluid dynamics van Sickle (1977), Kohyama et al (2003). This equation is also variously called the Fokker-Planck equation or the Kolmogorov forward equation: where D is the trunk diameter at breast height (DBH) in cm, n is the number of trees in a given size class per unit area in trees cm −1 ha −1 , g the trunk diameter growth rate in cm yr −1 , γ the mortality rate in yr −1 and t is time in yr. In this equation the first term represents the change in tree density per size class at a specific size D, the second represents the 'flow' of trees in size, through the point D, due to growth and the right hand term is the loss due to mortality. We solve this governing equation, in a similar way to Muller-Landau et al (2006), to derive analytical descriptions of equilibrium size distributions. This assumes demographic equilibrium, i.e. unchanging size distribution. While this is not strictly true for US forests, which are expected to have emerging dynamical responses to climate and land use (Zhu et al 2018), nevertheless, equilibrium solutions are a useful first approximation for a model whose goal is to help with DGVM development and understanding. Furthermore, Zhu et al (2018) used extensive ground measurements to show that US forests are well on the way to recovery from past disturbance and the US carbon sink is approaching saturation.
The equilibrium assumption means the growth rate g(D) and the mortality rate γ(D) are now simply time-independent functions of only the trunk diameter D and the first term of equation (1) is set as zero.
For the growth rate we assume this scales as a power law of trunk diameter g(D) ∝ D f , which was previously suggested by Niklas and Spatz (2004) and West et al (2009). This is also consistent with the recent observation that growth rate increases continuously with tree size ( The growth function is then defined as: where g 0 is the growth rate at fixed trunk diameter D 0 , which is the smallest sampled trunk diameter (only trees with trunk diameter of 12.7 cm or greater were used in our analysis, so we set D 0 =12.7 cm). For simplicity, we assume that the mortality rate is independent of tree-size. As mortality tends to be high for smaller trees and as these are excluded due to the sampling of the US dataset (only trees greater than 12.7 cm are sampled), the constant mortality assumption is a reasonable first approximation. Similar assumptions have been used previously (Muller-Landau et al 2006) and shown to fit measured size distributions of 150 tropical species in Panama with reasonable accuracy (Lima et al 2016).
By substituting (2) into the equilibrium form of (1) the equation can be solved for the tree density per size class n(D). So the solution for the case of power law growth and constant mortality is found to be (see supplementary material available online at stacks.iop.org/ ERL/13/084019/mmedia for derivation): where f is the growth scaling power, μ = γ D 0 /g 0 is the ratio of mortality rate to relative trunk diameter growth rate for a tree of diameter D 0 and n 0 is the tree density per size class at D 0 . During this study we chose to make the simplifying assumption of a fixed allometry where the trunk diameter growth rate scaling power f=1/3. This corresponds to the growth scaling used in theories based on hydraulic principles (Niklas and Spatz 2004) and metabolism and biomechanics (West et al 2009). A recent study (Duncanson et al 2015) has shown that for the US, MST, which is the basis of this allometry, applies well for trees that have a steady state, particularly old growth forest plots without recent disturbance. Although many of the US plots in the study were found not to be in a steady state, there is considerable variation for these forests and the allometry did approach the theoretical values asymptotically with tree height. Hence even though steady state does not apply consistently over the whole US this choice of allometry is the most consistent simplifying choice for this particular study, and is very applicable to the larger trees.
We extend the analysis by deriving (see supplementary material) a new equation which relates the mean trunk diameterD, in cm, to the ratio of mortality rate to growth rate μ, for the assumption of f = 1/3: where Γ represents an incomplete Gamma function. Equation (4) therefore suggests that when the growth scaling power f is fixed at a common value for all species,D can be expressed as a function of just one parameter μ.

Methods
We test the DET (3)   Weibull distribution (Wingo 1989, Zhang and Xie 2011, Lima et al 2016, is fitted by using maximum likelihood estimation (MLE), to find the value of the single parameter μ that maximises the likelihood (see section 2.2). The resulting best estimate of μ, together with the mean diameter obtained directly from the data of each species, can be used to compare to the theory relating these two properties of the data. The analysis is also repeated for the grouped PFT data and also for all 61 species grouped together as one large dataset.

Forest inventory data
The USDA Forest Service's FIA programme is the primary source for information about the extent, condition, status, and trends of forest resources in the United States (Oswalt 2014). FIA applies a nationally consistent sampling protocol using a quasi-systematic design across the United States, resulting in a national sample intensity of one plot per 2428 ha (Bechtold and Patterson 2005). Classified satellite imagery is used to identify forested land, which is defined as areas with at least 10% forest cover, at least 0.4 ha in size, and at least 36.6 m wide. In forest land, FIA inventory plots consist of four 7.2 m fixed-radius sub-plots spaced 36.6 m apart in a triangular arrangement with one sub-plot in the centre (Bechtold and Patterson 2005). All standing live trees with a DBH of at least 12.7 cm are inventoried on these forested sub-plots. Trees smaller than 12.7 cm diameter are sampled separately on four smaller 2.07 m radius micro-plots.
In this analysis, forest inventory data was extracted from 70 703 fully forested natural plots (non-plantation and non-disturbance) in the contiguous United States consisting of the lower 48 states, from FIADB version 6 on 20 March 2015 (available online http:// fia.fs.fed.us/). We restricted analysis to 61 species with sufficient sample sizes (>1000 plots), leading to 1442 517 trees in total. The threshold of 12.7 cm was selected as a lower level for inclusion in the analysis, as different sampling plots and techniques were used for trees smaller than this value. We excluded both plantations and dead trees from the dataset.

Fitting methodology
To fit to the DET distribution, MLE was used, as this has been shown to be an effective method for parameter fitting of forest size distributions (White et al 2008, Taubert et al 2013, particularly when the effective binning due to measurement precision is small, in this case to the nearest cm (Taubert et al 2013). To achieve this the probability density function (pdf) f (D), assuming demographic equilibrium and power law growth and constant mortality, can be obtained from (3): where N is the total number of trees in the dataset being fitted and A the area of the plots containing that PFT or species. This is equivalent to the standard form of the left-truncated Weibull distribution . For our chosen allometry this is then where f i is data point i in the dataset for that species or PFT. L is then maximised using Brentʼs bounded root finding algorithm (Brent 1973), from the Python SciPy library, to find the value of μ that gives the maximum log-likelihood. Once the parameter μ is estimated, then this allows n 0 , the tree density per size class at D 0 , to be obtained from (9) as the total number of trees N and the plot area A are known from the data: To obtain a confidence interval for μ a non-parametric bootstrap was used where the DBH distribution was sampled with replacement and then fitted using MLE in each instance. This was repeated 1000 times and the resulting distribution of results for μ then allowed an estimate of the 95% confidence interval.

Comparison with MST
To compare to MST the following equation was used : This model has no free parameters to fit via MLE as the pdf is : n 0 , the tree density per size class at D 0 , is obtained from N/A=n 0 D 0 as the total number of trees N and the plot area A are known from the data :

Predicting whole dataset mean diameter from only larger trees
To test the robustness of the fit in a practical application, the data fitting was repeated on a subset of each species and PFTs data. The subset chosen was to truncate the data to a new value > D D T 0 , leaving a smaller subset of larger size trees where > D D T . Then the fit obtained is used to predict the mean trunk diameter of not just the truncated data but the original full dataset.
Once the data has been truncated by D T then the DET fit methodology, using (7) and (8) (but with D T instead of D 0 ), is used to obtain the fitted parameter m T . As the μ parameter is a function of the truncation point, (D 0 or D T ) then a conversion is needed from the value fitted to the truncated dataset m T to that compatible with the untruncated dataset : . Now the mean trunk diameter can be calculated for the whole dataset from (4), using the μ obtained from the fit to the truncated data (13). This analysis may be helpful for earth observation satellites in space, which only resolve the large trees.

Results
All 61 species, as well as the grouped datasets for the four PFTs and all the species together (total of 66 fits) were fitted to the DET by MLE to obtain estimates of the mortality to relative growth ratio μ, and hence n 0 . The full table of results can be found in supplementary material table 1.
The DET fits the PFT size distributions very well with just one fitting parameter for each PFT. This can be seen when plotted against the data (see figure 2). The only significant deviation seen is for the broadleaf deciduous PFT at large tree sizes. So for this PFT the simplifying choices made in this study may be not as applicable and a different growth scaling power f or relaxing the constant mortality assumption may improve the fit.
The fit to all species together as one dataset in figure 3 (of a total of 1442 517 trees across the whole continental US) shows a very robust fit, despite the large scale and being a mix of species. In fact this largescale fit appears better than many of the individual species or the PFT groupings. The MST distribution is also plotted for each PFT and is significantly worse fit to the data than the DET fit in every case.
The size distributions in figures 2 and 3 were plotted in the more intuitive tree density per size class n(D) (trees per cm per hectare) rather than as a probability distribution function f (D).
The quality of the fits to the 61 individual species varies somewhat (see figures 2-8 in the supplementary material), with the most common deviations being a peak in the centre of the distribution, and in some cases the fit diverges for larger trees. Divergences from the demographic equilibrium will occur in forests that are not at steady state, for example as a result of forest management or natural disturbances, where early successional forests could deviate from the equilibrium assumption. Where the fit diverges at larger tree sizes could also indicate the different choices for growth scaling power f or mortality assumptions may give a better fit. Again the MST is inferior to DET for each species. Table 2 in the supplementary material shows the statistical assessment of the relative quality of fit of the two models. The log likelihood from the best fits of both models to each species and PFT are presented along with both the Akaike information criterion (AIC) and Bayesian information criterion (BIC) for both the DET and MST models. Both the AIC and BIC are a way of determining from several models which represents the data better, with a lower value indicating a better fit. For every single PFT and species the DET was the better fit to the data for both AIC and BIC.
In figure 4(a) the fitted DET mortality:growth ratios μ of each species (circles) are plotted as a function of the species mean trunk diameter and similar plot for the four PFTs is shown separately in the supplementary material. The mean trunk diameterD was calculated directly from the trunk diameter data for that species. The DET theoretical curve was plotted by solving (4) for a range of μ values and obtaining the mean trunk diameter for each of those values. Error bars, representing the 95% confidence interval, were omitted from figure 4(a) as they were with the exception of one species (Ostrya virginiana) too small to be visible on the figure.
The close fit between theory and data in figure 4(a) is to be expected as the MLE algorithm essentially fits to mean trunk diameter. However, the quality of the fit to the entire tree-size distribution is sufficient to enable reasonable estimates of mean D even when only the larger trees are fitted-to, as shown in figure 4(b). This would not be the case for alternative theories for the tree-size distributions such as MST or pure selfthinning profiles (μ=0), both of which predict an infinite mean D when integrating the size distributions out to infinity.
In figure 5 the variation of the mean diameter prediction of the whole dataset when only the larger trees are fitted-to is shown in relation to the truncation point used. The prediction is fairly reasonable and within 20% of the actual value for the whole range of truncation values. The results are quite similar for all the PFTs, with broadleaf deciduous having slightly worst performance over the whole range, this is consistent with the deviation of the fit in figure 2 for this PFT at large tree size.

Discussion and conclusions
This study has shown that, even assuming a highly idealised allometry (Niklas andSpatz 2004, West et al 2009), power law growth and constant mortality, the resulting DET is a suitable model of size distributions of most species across the US. The theory fits even better for the data grouped together into PFTs except for broadleaf deciduous, where it fits well at small and medium tree sizes but deviates at large tree sizes. The fit for the whole dataset of all species grouped together is particularly good.
The DET model is compared to the MST model and in every case the DET distribution solution is a better fit than the MST distribution solution. This has implications for MST, as while it has been shown before in many studies that the MST size distribution does not fit well to the observed forest distributions, the theory has tended to be rejected as a whole (Coomes et al 2003, Muller-Landau et al 2006, Coomes and Allen 2009, Anfodillo et al 2013, Anderson-Teixeira et al 2015. Here the DET theory has used the allometry from the MST to specify the scaling of tree growth with size, and has generally produced good fits to the data. This suggests that the allometric part of MST may be stronger than the space- Another theory (Anfodillo et al 2013) has been developed to explain forest size distributions based on height rather than trunk diameter and using the concepts of finite size-scaling and crown shape. This theory does raise the interesting question of sizelimitation to growth, which this study does not address. This study also follows the space-filling Figure 5. The mean trunk diameter for the whole of each dataset predicted from just a subset of larger size trees (D>D T >D 0 ), as a function of the truncation point D T . The DET is fitted to the truncated (D>D T ) dataset and the μ obtained is used to predict mean diameter for the whole dataset D>12.7 cm. The y axis is the predicted mean diameter divided by the actual mean diameter obtained directly from the data. Figure 4. (a) Mortality:relative growth ratios from distribution fits, for each species, as a function of mean tree trunk diameter. Each of the 61 species are plotted (circles) and colour coded by plant functional type, the mean trunk diameter of each species is determined directly from the data. The red dashed line shows the expected distribution based on the theory presented in this paper. The individual species follow the trend of the theory line. The close fit shown is not by itself validation of the theory. (b) The mean trunk diameter for the whole of each dataset predicted from just a subset of larger size trees (diameter D>50 cm). The DET is fitted to the truncated (D>50 cm) dataset and the μ obtained is used to predict mean diameter for whole dataset D>12.7 cm. The y axis is the predicted mean diameter divided by the actual mean diameter obtained directly from the data. arguments of West et al (2009), so the results here appear to contradict the assumptions of Anfodillo et al (2013), leaving an interesting future question as to how this may be resolved. Both models make different assumptions so it would be interesting to test the assumptions of each model against the other.
The size distribution solution has been further integrated to obtain, for the first time, a relationship between the mean trunk diameter of a dataset and a single parameter μ, which represents the mortality to growth ratio at trunk diameter D 0 =12.7 cm. We show that this relationship can allow reasonable predictions of the mean trunk diameter of the whole dataset from just a subset of the data consisting of only larger trees, this may have an application in estimating biomass via remote sensing. The challenge for remote sensing is that it cannot always easily observe small individual trees (Hauglin andNaesset 2016, Kellner andHubbell 2017), yet small trees are important components of the total biomass. Hence by calculating the fit to just the upper part of the distribution, it may then be possible to infer the rest of the size distribution from just the larger trees, which in turn can be resolved by Earth observation satellites.
The simplifying assumptions used in this study have implications for large-scale modelling such as in DGVMs. DGVMs need to accurately predict largescale carbon cycle interactions (Le Quéré et al 2009, Sitch et al 2015 and so must capture large-scale forest dynamics while relying as little as possible on empirical relationships, which may only be valid under particular environmental conditions (Chave et al 2004). Instead, DGVMs need to be based, as much as possible, on processes and theoretically derived relationships, if they are to accurately predict vegetation dynamics under environmental change. However, many DGVMs have become overly complex and this especially raises issues of appropriate parameterisation which may be unknown.
This study shows that rather simple assumptions (demographic equilibrium, constant mortality, powerlaw dependence of growth-rates on tree-size) are able to reproduce the demographic profiles observed in a very wide range of forests, and with just one site-specific parameter of the mortality:growth ratio μ. This also gives us confidence that the underlying dynamical model will provide an excellent basis for representing forest demography in next-generation DGVMs and ESMs. The sparsity of free parameters, and the availability of analytical equilibrium states for initialisation, are very attractive features for such global-scale applications.