A stand-alone tree demography and landscape structure module for Earth system models : integration with inventory data from temperate and boreal forests

Poorly constrained rates of biomass turnover are a key limitation of Earth system models (ESMs). In light of this, we recently proposed a new approach encoded in a model called Populations-Order-Physiology (POP), for the simulation of woody ecosystem stand dynamics, demography and disturbance-mediated heterogeneity. POP is suitable for continental to global applications and designed for coupling to the terrestrial ecosystem component of any ESM. POP bridges the gap between first-generation dynamic vegetation models (DVMs) with simple large-area parameterisations of woody biomass (typically used in current ESMs) and complex second-generation DVMs that explicitly simulate demographic processes and landscape heterogeneity of forests. The key simplification in the POP approach, compared with second-generation DVMs, is to compute physiological processes such as assimilation at grid-scale (with CABLE (Community Atmosphere Biosphere Land Exchange) or a similar land surface model), but to partition the gridscale biomass increment among age classes defined at subgrid-scale, each subject to its own dynamics. POP was successfully demonstrated along a savanna transect in northern Australia, replicating the effects of strong rainfall and fire disturbance gradients on observed stand productivity and structure. Here, we extend the application of POP to wide-ranging temporal and boreal forests, employing paired observations of stem biomass and density from forest inventory data to calibrate model parameters governing stand demography and biomass evolution. The calibrated POP model is then coupled to the CABLE land surface model, and the combined model (CABLE-POP) is evaluated against leaf–stem allometry observations from forest stands ranging in age from 3 to 200 year. Results indicate that simulated biomass pools conform well with observed allometry. We conclude that POP represents an ecologically plausible and efficient alternative to large-area parameterisations of woody biomass turnover, typically used in current ESMs.


Introduction
Changes in woody biomass storage in forest and savanna ecosystems, including woody ecosystems regenerating on abandoned agricultural lands, are the major driver of the terrestrial carbon sink, which currently amounts to around a quarter of anthropogenic emissions, mitigating climate change (Ahlström et al., 2012;Pan et al., 2011;Le Quéré et al., 2013).Such ecosystem dynamics and their feedbacks to atmospheric carbon content and radiative forcing are represented in Earth system models (ESMs) by incorporating dynamic vegetation models (DVMs).These attempt to describe changes in vegetation biomass components over time as the net effect of the allocation of net primary production (NPP), which increases or decreases biomass pools through phenological (seasonal) cycles of foliage and roots, mortality of plant individuals and disturbances such as wildfires and storms.The first-generation DVMs adopted by most current ESMs (Arora et al., 2013) employ large-area parameterisations designed for application on the scales of grid cells 10s to 100s of kilometres on a side.Typically these V. Haverd et al.: A stand-alone tree demography and landscape structure module for ESM parameterisations treat carbon flows associated with respiration and mortality as first-order decay processes, expressed as products of pool biomasses and bulk rate parameters independent of age structure (the "big wood" approximation; Wolf et al., 2011).These are computationally efficient -an important consideration for global-scale applications -but have the disadvantage of not resolving underlying population and community processes such as recruitment, mortality and competition between individuals and species for limiting resources (e.g.Sitch et al., 2003).This lack of mechanistic detail means that the models are unable to directly exploit the wealth of information on forest stand structure and dynamics available from forest inventories.These have been used to develop individual-based height-structured models that have been successfully used to simulate forest dynamics at the stand scale since the 1970s (e.g.Botkin et al., 1972;Bugmann, 2001;Smith et al., 2001).Different DVMs have also been shown to simulate widely different patterns and time evolution of biomass pools, especially under future climate projections (Cramer et al., 2001;Friedlingstein et al., 2006;Sitch et al., 2008;Friend et al., 2014) where models with a conservative response of biomass turnover to climatic forcing tend to retain a net biomass sink over the coming century, whereas others simulate a source or reduced sink by late 21st century (Ahlström et al., 2012).In ESM simulations with an active carbon cycle feedback to climate, such differences translate into divergence in the simulated global climate (Friedlingstein et al., 2006).It has been suggested that the representation of forest dynamics in ESMs may be one of the greatest sources of uncertainty in future climate projections (Purves and Pacala, 2008).
A handful of offline (not coupled to the atmosphere) second-generation DVMs exist that simulate demographic processes and landscape heterogeneity of forests using more explicit approaches that have been demonstrated to accurately replicate forest size structure and successional dynamics as predicted by community ecological theory.Examples include LPJ-GUESS (Smith et al., 2001) and ED (Moorcroft et al., 2001;Fisher et al., 2010).Such approaches are perceived as offering promise as an improved, second generation of DVMs (Purves and Pacala, 2008;Fisher et al., 2010).These models include stochastic representations of processes such as recruitment, mortality and large-scale disturbance, which requires replication of dynamic objects such as tree individuals and patches, and repeated computations of the same processes as applied to different objects, in order to obtain a representative average for the ecosystem as a whole.For studies of global and continental carbon balance, and for coupling to ESMs, however, this is a potential disadvantage because of their demand in terms of both of memory and processing power, and the additional complication that the results are not strictly deterministic, which complicates the analysis of results.In addition, the intricate internal representation of stand structure and its integration with plant physiological processes such as carbon assimilation, allocation and phenology implies that the enhancement of existing land surface models (LSMs) lacking or employing simpler parameterisations of vegetation dynamics may be a time-consuming, technically challenging task.Wolf et al. (2011) recently used global forest inventory data to assess forest biomass allometry in eight global land surface models, including two second-generation DVMs.Simulated relationships between stem and foliage biomass pools generally conformed poorly with observed allometry, indicative of model failure to consistently reproduce both structural and functional characteristics of vegetation.Best overall performance was noted for the models ED and ORCHIDEE-FM (Bellassen et al., 2010), which include an explicit parameterisation of self-thinning, which strongly controls biomass turnover rates in closed-forest ecosystems (Westoby, 1984).The study recommended the use of biomass allometry data from forest inventories as a simple approach to improving the characteristic behaviour of global land surface models with respect to structural dynamics.
To simultaneously overcome the limited ecological realism of simulated wood turnover in many first-generation DVMs, and the technical limitations of current secondgeneration DVMs, Haverd et al. (2013) proposed a new approach for the simulation of woody ecosystem stand dynamics, demography and disturbance-mediated heterogeneity.The approach, encoded in a model called Populations-Order-Physiology (POP), is designed to be modular, deterministic, computationally efficient and based on sufficient ecological realism for application at the grid scales typically employed by DVMs and ESMs for continental to global applications.Coupled to the CABLE (Community Atmosphere Biosphere Land Exchange; Wang et al., 2011) LSM, POP receives woody biomass increment from CABLE and returns an updated biomass state (an approach conceptually similar to that of ORCHIDEE-FM, but with key differences discussed in Sect.4.2).CABLE-POP was demonstrated along a savanna transect in northern Australia, successfully replicating the effects of strong rainfall and fire disturbance gradients on observed stand productivity and structure (Haverd et al., 2013).The key simplification in the POP approach, compared with second-generation DVMs, is to compute physiological processes such as carbon assimilation at grid-scale (with CABLE or a similar land surface model), but to partition the grid-scale biomass increment among age classes defined at sub-grid-scale, each subject to its own dynamics.POP is not a new DVM but a scheme for dynamically estimating size structure and turnover of woody vegetation, forced by productivity information from an external LSM.
In the present study, we extend the application of POP to globally distributed forests, heeding the recommendation of Wolf et al. (2011) to constrain and improve the performance of the model by using allometric scaling relationships from forest inventory data.Thus calibrated, the combined model (CABLE-POP) is evaluated against leaf-stem allometry and total biomass observations from forest stands.POP is described in Appendix 1 of Haverd et al. (2013), and the detailed description (Appendix A) and summary below are largely reproduced from that paper.For the purpose of the present application, which includes closed-forest ecosystems, we extended tree mortality in POP to include a crowding component as described below.A Fortran 90 version of the POP computer code is included in the supplementary material to this paper.
POP is designed to be modular, deterministic, computationally efficient, and based on defensible ecological principles.Parameterisations of tree growth and allometry, recruitment and mortality are broadly based on the approach of the LPJ-GUESS DVM (Smith et al., 2001).The time step ( t) is 1 year.
POP is designed to be coupled to a land surface model (LSM) or the land surface component of an ESM (Sect.2.1.2) which provides forcing in terms of the annual grid-scale stem biomass increment ( C (kg C m −2 )) for woody vegetation, as an average across a simulated tile or grid cell.In LSMs such as the CABLE model employed in the present study (see below), each tile represents the proportion of a grid cell dominated by one major vegetation type, such as evergreen or deciduous forest.For scaling to the landscape (tile or grid cell) scale, POP also requires mean return times of exogenous (large-scale) disturbances.For the present study, where we focus on the patch-scale size structural dynamics, we adopt a mean "catastrophic" disturbance return time of 100 years, which kills all individuals (cohorts) and removes all biomass in a given patch.POP can also account for "partial" disturbance, such as fire, which results in the loss of a size-dependent fraction of individuals and biomass, preferentially affecting smaller (younger) cohorts.However, this feature is not used here.Stem biomass increment is provided by the host LSM (here CABLE) or prescribed for stand-alone calibration.
State variables are the density of tree stems partitioned among age classes (cohorts) of trees and representative neighbourhoods (patches) of different age since last disturbance across a simulated landscape, representing a spatial unit (tile or grid cell) of an LSM.Hereinafter we use the term "grid cell" to refer to the spatial unit at which POP is coupled to the host LSM, in our study a vegetation tile comprising either evergreen needleleaved or deciduous broadleaved forest.A patch thus represents a stand of vegetation of sufficient extent to encompass a neighbourhood of individual woody plants, competing with one another in the uptake and utilisation of light, soil resources and space.Patches are not spatially referenced, but represent a statistical sample of local stand structure within the overall landscape of the grid cell.Trees are assumed to belong to the plant functional type (PFT) defined for the tile or grid cell in the host LSM.Individuals are not distinguished within a cohort, but each cohort has a diagnostically varying mean individual stem biomass (see below), from which other size metrics (height, stem diameter and crown area) can be derived (see Appendix A4).
POP thus simulates allometric growth of cohorts of trees that compete for light and soil resources within a patch.The annual stem biomass increment is partitioned among cohorts according to a power function of their current aggregate stem biomass (size), on the assumption that larger individuals preempt resources owing to a larger surface area and exploration volume of their resource-uptake surfaces (leaves and fine roots), and due to the advantage conferred on taller individuals by the shading of shorter ones in crowded stands (Westoby, 1984).A cohort's share of the total annual biomass increment is divided equally among individuals.Thus, between cohorts, there is shading implicit in the weighting of the biomass partitioning towards larger trees/cohorts.
The mortality parameterisation was specifically updated for this study and is therefore described here in detail.Key model parameters are listed in Table 1.
Population dynamics are governed by where N y is the stem density of the cohort established in year y, and m R,y and m C,y are cohort mortalities (yr −1 ) due to resource limitation and crowding respectively.N y is initialised as recruitment density, and is reset (according to disturbance intensity) when the patch experiences disturbance.
The mortality rate for a cohort depends on the growth efficiency (GE), closely related to the concept of the relative growth rate (RGR), given by where C y (kg C m −2 ) is the stem biomass and C y is the annual stem biomass increment of the yth cohort; s is set to 0.75, which is the power governing the mean proportionality between plant resource-uptake surfaces (leaves and roots) and stem biomass for a wide range of plant taxa and vegetation types according to Enquist and Niklas (2001).GE thus represents annual growth relative to the estimated area of resource-uptake surfaces for trees of a particular size.We characterise the response of resource-limitation mortality to GE by a logistic curve with the inflection point at GE min : where y is the index for a particular cohort and m R,max (yr −1 ) is the upper asymptote for mortality as GE declines, a proxy for the resilience of plants to extended periods of resource stress, and is set to the value of 0.3 adopted in the LPJ-GUESS DVM (Smith et al., 2001).The exponent p, assigned a default value of 5, governs the steepness of the response of m R to GE around GE = GE min .For this study GE min was set to its calibration value of 0.015, as determined previously by optimisation against northern Australian tree basal area data (Haverd et al., 2013).
The partitioning of the grid-cell-level annual biomass increment among cohorts in a patch is governed by where C is the grid-cell annual biomass increment (assumed to be partitioned equally among patches), t is the time step of 1 year and the summation is over all cohorts in the patch.Individuals within a patch are thus assumed to capture resources in proportion to the area of their resource-uptake surfaces, estimated as the s power of stem biomass following the allometric scaling theory of Enquist and Niklas (2001).
The additional crowding mortality component (m C,y ) was included to allow for self-thinning in forest canopies.Selfthinning is dependent on the assumption that some trees (within a cohort) have a slight advantage in pre-empting resources, creating a positive feedback to their growth and ultimately resulting in death of the most suppressed individuals.In contrast, in POP, the total stem biomass increment for a cohort is equally partitioned amongst all members.To compensate for this simplification, we use the following parameterisation, which emulates the contribution to self-thinning associated with within-cohort competition: Equation ( 5) implies that crowding mortality never exceeds growth.Here, c pc,y is crown projective cover (Appendix A5), A c,y crown projected area (m −2 m −2 ) of all crowns in the yth and taller cohorts, α C a coefficient which determines the onset of crowding mortality with respect to c pc and f c is a tunable scaling factor.α C was set to 10.0, corresponding to an onset of crowding mortality at c pc ∼ 0.8.This value implies that crowding mortality is insignificant in the Australian savanna simulations, thus retaining the validity of the parameters relating to m R,y in Eq. ( 5) as used in the earlier study of Haverd et al. (2013).Crown projected area is evaluated as where N y is stem density (m −2 ); D y is stem diameter at breast height (m) (Appendix A, Eq.A6-A9); and k allom and k rp are parameters set to respective values of 200 and 1.67 respectively, based on literature values compiled by Widlowski et al. (2003).
Additional mortality occurs as a result of disturbances.Replicate patches representing stands of differing age since last disturbance are simulated for each grid cell.It is assumed that each grid cell is large enough to accommodate a landscape in which the frequency of patches of different ages follows a negative exponential distribution with an expectation related to the current disturbance interval.This assumption is valid if grid cells are large relative to the average area affected by a single disturbance event and disturbances are a Poisson process, occurring randomly with the same expectation at any point across the landscape, independent of previous disturbance events.To account for disturbances and the resulting landscape structure, state variables of patches of different ages are averaged, and weighted by probability intervals from the negative exponential distribution.The resultant weighted average of, for example, total stem biomass or annual stem biomass turnover is taken to be representative for the grid cell as a whole.Strictly, the Poisson assumption demands that the mean disturbance interval is invariable over time, a difficult assumption to uphold in practice, as disturbance agents such as wildfires, windthrow and pest or pathogen attacks may increase or decrease depending on variations in climate and other drivers.A constant past disturbance regime was assumed in the present study.

CABLE-POP
CABLE is a global land surface model consisting of five components (Wang et al., 2011): (1) the radiation module describes direct and diffuse radiation transfer and absorption by sunlit and shaded leaves; (2) the canopy micrometeorology module describes the surface roughness length, zero-plane displacement height and aerodynamic conductance from the reference height to the air within canopy or to the soil surface; (3) the canopy module includes the coupled energy balance, transpiration, stomatal conductance and photosynthesis of sunlit and shaded leaves; (4) the soil module describes heat and water fluxes within soil and snow at their respective surfaces; and (5) the CASA-CNP biogeochemical model (Wang et al., 2010).In this study, we used CABLE-2.0 with the default the soil module replaced by the Soil-Litter-Iso (SLI) model (Haverd and Cuntz, 2010).
As illustrated in Fig. 1, coupling between CABLE and POP is achieved by exchange of two variables: CABLE supplies annual stem biomass increment to POP and POP returns an annual stem biomass loss to CABLE.To convert between stem biomass (POP) and tree biomass (CABLE), we assume a ratio of 0.7, a representative average for forest and woodland ecosystems globally (Poorter et al., 2012).The resulting tree biomass turnover is applied as an annual decrease in the CABLE tree biomass pool, and replaces the default fixed biomass turnover rate.
The model set-up in this study was designed to permit evaluation of CABLE-POP predictions of leaf-stem allometry.CABLE-POP was run offline at 1 • × 1 • spatial resolution for grid cells containing the locations of forests in the Cannell-Usoltsev (C-U) database (see Sect. 2.2 on data below).Simulations were forced using GSWP-2 (Globa Soil Wetness Project) 3-hourly meteorology for the 1986-1995 period (Dirmeyer et al., 2006).Leaf area index (LAI) was prescribed using a monthly climatology from the MODIS (Moderate Resolution Imaging Spectroradiometer) Collection 5 product (Ganguly et al., 2008).Vegetation cover was prescribed as one of three of the CABLE PFTs: evergreen needleleaf, evergreen broadleaf or deciduous broadleaf, each with its own set of physiological parameters.Note that CABLE does not distinguish between cold-and droughtdeciduous broadleaved vegetation.Needleleaf and broadleaf were distinguished based on the classification in the C-U database.All needleleaf forests were assumed evergreen, and broadleaf forests were classified as deciduous or evergreen according to the larger area fraction specified in the vegetation distribution data set by Lawrence et al. (2012).In cases with no information on either, a distinction was made by the location, with broadleaf forests north of 17 • N assumed deciduous.
The modelling protocol was as follows: (i) CABLE soil moisture and temperature were initialised by running CA-BLE (without CASA-CNP) once for 10 years (using the 10 year meteorological data record); (ii) CABLE (without CASA-CNP) was run a second time for 10 years from this initial state, this time with daily forcing inputs to CASA-CNP being saved, namely gross primary production, soil moisture and soil temperature; and (iii) CASA-CNP was run for 400 year (40 × 10 year of repeated forcing) at daily time step, with POP being called annually and initial biomass stores set to zero.

Data
Forest inventory data for total biomass, stem biomass, foliage biomass and stem density were sourced from the Biomass Compartments Database (Teobaldelli, 2008).This database contains data from around 5790 plots and represents a harmonised collection of existing data sets (Cannell, 1982;Usoltsev, 2001), covering the temperate and boreal forest region globally.The data include the following compartments: stem, bark, branches, foliage, roots, fruits, dead wood and understorey.Latitudes and longitudes were rounded to the nearest degree centred on the half-degree, and the data were separated into broadleaf and needleleaf groups, with "mixed forest" sites removed.Latitude/longitude duplicates were then removed separately for each of the needleleaf and broadleaf subsets, leaving all but one randomly selected occurrence in each 1 data required for our analysis.For comparison with model output, the data were further filtered, such that only plots with data for stem biomass, foliar biomass, stem density and age were retained, leaving 178 broadleaf plots and 304 needleleaf plots.Hereafter we refer to the data for these plots as the "C-U data".Their locations are denoted in Fig. 2. Average stem biomass M stem (kg tree −1 ) and foliar biomass M fol (kg tree −1 ) per tree were obtained by dividing the bulk quantities by stem density (N).Total biomass per tree (M) was estimated as the sum of woody, foliar and fine root biomass, assuming allometric ratios of stem biomass to total woody biomass (0.7) (Widlowski et al., 2003) and fine root to foliar biomass (1.0) (Luyssaert et al., 2007).
In this study, we use the C-U data in three ways: (i) to construct a biomass-density (log(M) vs. log(N)) plot for the purpose of calibrating the crowding mortality component of POP (Sect.2.3 below); (ii) to construct leaf-stem allometry plots (log(M fol ) vs. log(M stem )) for the purpose of evaluating the CABLE-POP scaling exponent (slope) relating M fol to M stem , and for tuning the CABLE allocation coefficients to leaves and stems, to which the intercept is sensitive (Sect.3 below); and (iii) to evaluate CABLE-POP predictions of stem biomass directly against data (Sect.3 below).

Calibration
The crowding mortality component of the POP model was calibrated using average biomass per tree (M) (kg dry matter per tree) and stem density (N) (trees ha −1 ) data from the combined broadleaf and needleleaf data sets.These variables can be plotted in the form of the self-thinning "law" (e.g.Westoby, 1984): which describes the ageing trajectory of forest stands after they exit the initial density-independent growth phase and before the stand is sufficiently self-thinned that mortality be- comes density-independent. (Hereafter all log functions refer to log 10 ).The self-thinning part of the trajectory forms the upper bound of a plot of log(M) vs. log(N ) (Fig. 3), with points below this upper bound resulting from young stands in the density-independent growth phase and additional disturbance-related mortality beyond that described by self-thinning.Thirty-three points along this upper bound were selected for POP calibration.This was done by binning the data in Fig. 3 into 39 evenly spaced bins spanning a log(N) range of 2.3-4.5, and selecting the observation corresponding to the maximum value in each non-empty bin.The coefficients in Eq. ( 7) were estimated from these points using reduced major axis (rma) regression (e.g.Sokal and Rohlf, 1995, Sect. 14.13) and treated as observations.The corresponding model observables were constructed from stand-alone POP simulations of stands with the same age and CABLE-estimated annual stem increment (hereafter Stem-NPP) as the observations, and with a high initial stem density (3 individuals m −2 ) to accelerate the progress of young stands towards self-thinning behaviour.The residuals between modelled and observed coefficients of Eq. ( 7) were minimised by optimising the f c parameter (Eq.A11) using the PEST parameter estimation software (Doherty, 2004) which implements the Levenberg-Marquardt down-gradient search algorithm.This returned a value of f c = 0.013 ± 0.007 (1σ ).All other POP parameters were held fixed at their prior values (Haverd et al., 2013) to ensure that the model parameter set is equally valid for the savanna landscape (to which the model was initially applied) as for simulation of forest stands in the present study.
Data points not lying on the upper bound were not selected for calibration because these observations are not expected to be influenced by crowding mortality.The key indicator of the calibrated model's skill is its representation of the leaf-stem allometry plots (log(Mfol) vs. log(Mstem)) (Sect. 3.1,Figs. 4 (ii) and (iv)), of which none of the observation points were used directly in calibration.
CABLE parameters were held fixed at their default PFTspecific values, except for allocation coefficients of evergreen needleleaf forests and deciduous broadleaf forests, which were manually tuned to match the intercept of the leafstem allometry plots.Resulting proportions of NPP allocated to leaf, wood and fine roots respectively are [0.210.29 0.50] (evergreen needleleaf) and [0.33 0.37 0.30] (deciduous broadleaf).Corresponding values for evergreen broadleaf forests were held fixed at their default values of [0.20 0.35 0.45] because this PFT is underrepresented in the data (Fig. 2).

Comparison with observations
Figure 4 shows results of CABLE-POP simulations.Each simulation point represents a single patch with age matched to the age of the corresponding C-U data point.POP parameters were kept the same as for the POP calibration run, with no distinction between PFTs.The initial stem density was set to 2 stems m −2 , to approximately match the upper limit of N in the observations.
The CABLE-POP simulations in the density-biomass plots (i and ii) lie along the upper bound of the observations.It is not expected that the model should capture the distribution of the scattered observations below this bound: these observations are likely to correspond to (i) young stands in the density-independent growth phase (with density highly dependent on initial stem density, the variability of which is not captured in POP); (ii) very old stands of declining biomass in which N is decreasing while M is approximately constant (the total biomass of a very old stand ultimately declines, in part due to reduced productivity arising from physiological decline and nutrient limitations (Dewar, 1993;Gower et al., 1996).In a global context, this effect will be important mainly in the few global regions in which natural disturbance regimes and human management scarcely limit mean stand age.For this reason, and in the interest of model parsimony, we do not attempt to represent a declining biomass trend in very old stands here); and (iii) stands which have undergone managed thinning, particularly prevalent amongst needleleaf stands.Hence the discrepancy between linear fits to the predictions and observations (Table 2) is expected.
In contrast, the linear fits to the CABLE-POP predictions and observations in the leaf-stem allometry plots (iii and iv) (see also Table 1) agree very well, and generally better than the corresponding fits derived for other LSMs by Wolf et al. (2011).Note here that M fol and M stem are average foliage and stem biomass per tree (kg DM tree −1 ).
As noted by Wolf et al. (2011), a major impediment to validating models directly against measured biomass is the need to consider the many idiosyncrasies of each forest stand (e.g.species mix, climate, water/nutrient limitations, timing of disturbances, management).Nonetheless CABLE-POP simulations of biomass in broadleaf forest stands (Fig. 3v) are largely unbiased (slope = 0.94 ± 0.04 when intercept set to zero) and capture a high proportion of the variance (r 2 = 0.57).Total biomass in needleleaf stands (Fig. 3vi) is less well predicted (slope = 0.94 ± 0.3 when intercept set to zero, r 2 = 0.23), a likely consequence of intensive management, particularly deliberate thinning.Thinning is performed for economic reasons (e.g.Aruga et al., 2013) or to promote stand health (e.g.Ronnberg et al., 2013) and would reduce tree density while leaving the average stem biomass initially unaffected, resulting in a shift to the right for affected stands in Fig. 3iv (consistent with the high density of outliers).This would also explain the overestimation of low-biomass stands by CABLE-POP in Fig. 3vi because biomass would be removed in the early stages of a stand.Furthermore, as stated by Law et al. (2013), multi-stage thinning can also lead to an enhanced storage of long-term biomass, which explains why CABLE-POP overestimates the younger stands' biomass but does not reach the maximum values of the needleleaf stands in the C-U data.Uncertainties in woody increment prediction may also contribute to poor predictions of total biomass in needleleaf stands.
Figure 5 shows biomass component fractions extracted from CABLE-POP patch-scale simulations, and compared with estimates derived from the Cannell and Usoltsev databases by Wolf et al. (2011).The CABLE-POP simulations reproduce the major features of the data, particularly the sharp decline in the fraction of foliage biomass, and the relatively large foliage biomass fraction associated with needleleaf stands compared to broadleaf stands.In contrast, default CABLE simulations fail to capture the sharp decline in the fraction of foliage biomass, because they assume a fixed turnover time for biomass turnover (40 year for broadleaf and 70 year for needleleaf).The root component is dominated by the coarse-root fraction, which in our simulations was a constant proportion of woody biomass.Therefore we do not expect CABLE-POP or default CABLE to reproduce the observed decline in root : shoot ratio with total biomass.(i) and (ii) biomass-density plot; (iii) and (iv) leaf-stem allometry plot (M fol and M stem are average foliage and stem biomass per tree (kg DM tree −1 )), including results derived for other LSMs by Wolf et al. (2011); and (v) and (vi) total biomass: predictions vs. observations.In (i-iv), lines denote linear fits to the observations and predictions.Solid lines denote linear regression fits to the data points (see Table 1 for regression coefficients).

POP mortality dynamics
Figure 6 illustrates the dynamic behaviour of POP mortality via stand-alone POP simulations of two undisturbed patches with low and high extremes of annual stem biomass increment: StemNPP = 0.05 kg C m −2 yr −1 (low production) and StemNPP = 0.20 kg C m −2 yr −1 (high production).For reference, CABLE simulations of the C-U stands give average annual stem biomass increments of (0.17±0.06 1σ ) (broadleaf) and (0.16 ± 0.05 1σ ) (needleleaf) kg C m −2 year −1 .Figure 6i shows the ageing trajectory of each patch in biomassdensity space, with points representing every fifth year.The low-production patch exhibits an initial increase in density as recruitment augments the population during initial years (Fig. 6vii), before rapidly transitioning to a regime of declining stem density, characterised by a slope of −1, as resourcemediated stress induces mortality, cancelling any net increase in stand biomass (Fig. 6v).The ageing trajectory of the lowproduction patch never reaches the upper bound of the C-U data (representing self-thinning due to crowding mortality) because the stand is relatively sparse (Fig. 6ii) and resourcestress mortality prevents crowding (Fig. 6iii).In contrast, the high-production patch (Fig. 6i) experiences only a brief increase in density (Fig. 6viii), before transitioning to a regime of declining stem density following a trajectory analogous to the upper bound of the C-U data (slope: −1.45), corresponding to domination of crowding mortality (Fig. 6iv) due to high crown projective cover (Fig. 6ii) and initially low resource limitation mortality (Fig. 6iv).Resource-limitation mortality governs the level and rate of approach to equilibrium biomass at which mortality (population level) cancels the aggregate effects of individual tree growth, i.e. a slope of −1 on the ageing trajectory (Fig. 6i).Resource-limitation mortality rises more slowly initially than crowding mortality (Fig. 6iii, iv), responding to declining growth efficiency with size (trees investing more in maintenance costs of growing tall), while crowding occurs relatively rapidly as individuals spread horizontally to maximise their light uptake and growth.Turnover rate coefficients (Fig. 6iii, iv) simulated by POP increase with age, in contrast to the timeinvariant turnover rate coefficients assumed by default CA-BLE, and other LSMs employing the big wood assumption.
The big wood assumption leads to relative higher mortality in younger stands, particularly in the low-productivity example (Fig. 6iii).The big wood assumption also leads to turnover rate coefficients being invariant with productivity.One consequence of this is illustrated in Figs. 5 and 6, which show that for the fixed rate coefficient chosen here (0.02 year −1 ) the biomass of the low-productivity patch is adequately simulated, but that of the higher productivity patch is significantly underestimated compared with the POP simulations.Figure 6vii and viii show that there are fewer size (age) classes at higher NPP, reflecting stronger dominance by the tallest cohort in high-production stands, where deep shading beneath the upper canopy promotes high resource-limitation mortality for all but the dominant cohort, as expected from Eqs. (3-4).

Limitations of big wood models
The term "big wood approximation" was coined by Wolf et al. (2011) to describe the representation of woody vegetation biomass dynamics in the majority of LSMs considered in their review.These effectively treat woody biomass (M wood ) as a single carbon pool obeying first-order kinetics, with a rate coefficient (k) linearly scaling a single woody biomass pool: dM wood /dt = −kM wood + α wood NPP, where α wood is the fraction of NPP allocated to wood.The big wood approach is unrealistic as it compounds the differential responses to environmental drivers and system state of tree growth (generally positive) and population growth (positive or negative, depending on the balance between recruitment and mortality).Cohorts of different age will face different mortality rates depending on the microenvironment imposed by realised stand structure, and this in turn will vary among patches in a landscape, depending on the disturbance history.This suggests that it is important to specifically simulate the size distribution of trees in forest vegetation and account for how changes in this distribution may alter the response of ecosystem functions like biomass carbon turnover to forcing variables.
As plotted in Fig. 6iii and iv modelled biomass turnover rate -an emergent property in POP -increases with stand age (as commonly observed in real forest stands; Franklin et al., 1987) before reaching an equilibrium value, which differs between stands of differing productivity.In accordance with the discussion of Wolf et al. (2011), this suggests that big wood models cannot be expected to perform well for globally important young forest stands, instead being more applicable to relatively rare older stands at equilibrium biomass.
A second limitation of big wood models, also emphasised by Wolf et al. (2011), is that they are not readily amenable to validation against forest inventory data, such as those used in this study.This is because they do not carry information about tree density.Wolf et al. (2011)  this problem by applying a fit to biomass-density (log(M) vs. log(N)) as a post hoc estimate of density, which was then used to evaluate component biomass per tree, as required e.g. for the leaf-stem allometry plots (Fig. 4iii, iv).(Results are denoted "other LSMs" in these plots.)N was thus estimated using log where M × N is total biomass, and β and α are respectively the slope and intercept of a reduced major axis regression fit to log(M) vs. log(N) observations for the whole (global) data set (Fig. 4).As an approach for estimating N at a gridcell level, this assumes that forests throughout the world are following the same log density-log biomass trajectory.We suggest this is unlikely, as individual stands may be expected to follow different trajectories depending on productivity and age, as illustrated in Fig. 6 and discussed above.We applied Eq. ( 8) to CABLE-POP grid cell estimates of total biomass to derive a post hoc prediction of stem density, and hence average foliage and stem biomass per tree, analogous to the approach of Wolf et al. (2011).Results are shown in Fig. 7, and indicate that the grid-cell results (deduced via Eq.( 8) with α = 6.22 and β = −1.32(Wolf et al., 2011, Table 5)) lie on a significantly different line to the patch-level CABLE-POP simulations (individual points shown in Fig. 4iii, iv), for which the internal model stem density was used to deduce average foliage and stem biomass per tree.The same was true when values of α = 6.9 and β = −1.44 were used, corresponding to a fit to log(M) vs. log(N ) used in the present study (Fig. 3).Allometric data from global forest inventories are thus of limited value for evaluating/constraining big wood models, which do not carry number-density information.
Predicted global carbon fluxes diverge markedly between Earth system models, both in the magnitude and shape (sign, timing of change) of the subsequent trajectory (Ahlström et al., 2012;Anav et al., 2013).One cause of such divergence is that models, many employing a big wood simplification for biomass dynamics, differ markedly in terms of woody biomass turnover and its response to future climate and CO 2 forcing.In projections with ESMs that account for carbon cycle feedbacks, this translates into divergence in the simulated climate (Friedlingstein et al., 2006(Friedlingstein et al., , 2013)).We concur with Wolf et al. (2011) and argue that big wood models lack ecological realism and cannot be expected to simulate woody biomass turnover in a realistic manner under changing climate forcing.They should be phased out from use in carbon cycle studies.

Limitations of the POP approach
The approach presented and demonstrated in this paper offers a potential alternative, suitable as a replacement for the big wood approximation, for the representation of biomass structural dynamics for woody vegetation in large-scale models.POP is not a replacement for a full-featured DVM.It does not represent biogeochemical processes, nor in its current form competitive interactions among PFTs, nor biophysical feedbacks associated with age-structured vegetation.Further, age effects on LAI and NPP are not accounted for because of the simplifying assumption that NPP (and therefore LAI) is uniform among patches, while structure is assumed to vary.
POP is designed to be readily coupled to a biogeochemical LSM as implemented in many current climate and Earth system models.Such LSMs generally do not feature competing plant types, but may have fixed tiles representing gridcell fractions dominated by different types of vegetation.POP simulates size structure dynamics separately for each (woody) vegetation tile, based on the principle of asymmetric (i.e.size-dependent) competition between co-occurring individuals, but with no competition among PFTs (tiles).Competition between trees and grasses, deciduous and evergreen vegetation, and C 3 and C 4 plants provides an important explanation for global biome distributions and may modulate the responses of vegetation to future climate and [CO 2 ] forc-ing (Smith et al., 2014).We plan to introduce PFTs and to distinguish canopy and understorey strata in a later development of the approach.
In its current form, vegetation structure represented in POP only feeds back on woody biomass turnover represented in  3 (iii) and (iv).CABLE-POP grid-cell results are from deducing stem density from total grid-cell biomass (average over multiple patches, weighted by probability of time since last disturbance with a mean disturbance interval of 100 year) via Eq.( 8), and using this to compute M fol and M stem from grid-cell foliage and stem biomass components.the host LSM.Other feedbacks, such as the acclimation of allocation with age/size (Lloyd, 1999) and age-related decline in forest growth (e.g.Zaehle et al., 2006), could also be implemented using the modular approach presented here.Indeed, such feedbacks are represented in the ORCHIDEE-FM model (Bellassen et al., 2010), which adopts a similarly modular approach: the forest management module (FM) in ORCHIDEE-FM is supplied with stem biomass increment from the host LSM (ORCHIDEE) and returns biomass turnover.In other respects FM and POP are different.In particular, POP is distinguished by its representation of subgrid-scale disturbance-mediated heterogeneity.Further, the stand dynamics of POP are applicable from sparse woody savannas to closed forest.In contrast, FM was developed specifically for forests: it represents management effects, and natural mortality is based exclusively on stand self-thinning.
The parameterisation of crowding mortality Eq. ( 5) is independent of growth rates and assumes only that the proportionality between log biomass and log density revealed by the observed data continues to hold true in a future simulation.This assumption is reasonable if mortality depends more on biomass than on growth rate.The robustness of this proportionality across wide climatic gradients is clear from the observed data as plotted in Fig. 3.
The present paper is concerned with the patch-scale dynamics -not the landscape scaling component of POP, which was presented and applied in Haverd et al. (2013).However, we expect that the equilibrium assumption implicit in the Poisson-based method of weighting patches of different time since last disturbance will be robust to a gradually changing disturbance regime (typical for most climate change studies), though the approach will break down in the case of a step change in mean disturbance interval.

Conclusions
By discriminating individual and population growth and explicitly representing asymmetric competition among age/size classes of trees co-occurring within forest stands, POP overcomes a key limitation of the big wood approach and proved able to reproduce allometric relationships reflecting linkages between productivity, biomass and density in widely distributed temperate and boreal forests.Coupled to a biogeochemical land surface model, able to prognose woody biomass productivity at stand (or grid cell) level, POP may be calibrated and evaluated against forest inventory data, as demonstrated in this study.This is achieved without a marked increase in model complexity or computational demand, thanks to a modular design that separates the role of the parent land surface model (prognosing whole-ecosystem production) and the population dynamics model (partitioning the production among cohorts, computing mortality for each and returning the stand-level integral as whole-ecosystem biomass turnover to the parent model) (Fig. 1).
The present paper focuses on stand-level demographics and its influence on the accumulation and turnover of stem carbon biomass.At the landscape level, the incidence and intensity of disturbance events -such as wildfires, storms or anthropogenic interventions, such as forest harvest or land use conversions -provide an additional, regionally important control on biomass accumulation and turnover (Shevliakova et al., 2009).Such landscape-level effects are accounted for by our model, for example along a rainfall-mediated wildfire and biomass gradient in savannah vegetation of northern Australia (Haverd et al., 2013).

V. Haverd et al.:
A stand-alone tree demography and landscape structure module for ESM of the LPJ-GUESS DVM (Smith et al., 2001) we aggregate multiple causes of mortality in real tree stands within a generalised resource-limitation mortality that increases with size (and therefore age) and under conditions of declining productivity, whether caused by interference in resource uptake among neighbours (Westoby, 1984), abiotic factors such as drought or secondary biotic factors (e.g.pathogen attacks) (Franklin et al., 1987).Resource limitation leads to a decline in growth efficiency (GE, i.e. growth rate as a function of size), characterised by the ratio of stem biomass increment to current stem biomass for a given cohort.Below a certain threshold GE min , mortality increases markedly (Pacala et al., 1993).We characterise the response of resource-limitation mortality to growth efficiency by a logistic curve with the inflection point at this threshold: where GE y = C y /C s y , s takes the same value as in Eq. (A5) and m R,max is the upper asymptote for mortality as GE declines, set to 0.3 year −1 following Smith et al. (2001).The exponent p, assigned a default value of 5, governs the steepness of response around GE = GE min , which may be regarded as a calibration parameter.
An additional crowding mortality component (m C,y ) is included to allow for self-thinning in forest canopies.Selfthinning is dependent on the assumption that some trees (within a cohort) have a slight advantage in pre-empting resources, creating a positive feedback to growth and ultimately resulting in death of the most suppressed individuals.In contrast, in POP, the total stem biomass increment for a cohort is equally partitioned amongst all members.Therefore we require the following new parameterisation which emulates the contribution to self-thinning associated with withincohort competition.
Crowding mortality is expressed as such that it never exceeds growth.Here c pc,y = 1 − exp −A c,y is crown projective cover, A c,y crown projected area (m −2 m −2 ) of all crowns in the yth and taller cohorts, α C a coefficient which determines the onset of crowding mortality with respect to c pc and f c is a tunable scaling factor.α C was set to 10.0, corresponding to an onset of crowding mortality at c pc ∼ 0.8.This value was chosen to be sufficiently high such that crowding mortality is insignificant in the Australian savanna simulations (Haverd et al., 2013), thus retaining the validity of the parameters relating to m R in Eq. (A10) as used in this earlier study.Crown projected area is evaluated as where N y is stem density (m −2 ); D y is stem diameter at breast height (m), and k allom and k rp are parameters set to respective values of 200 and 1.67, based on literature values compiled by Widlowski et al. (2003).
Additional mortality occurs as a result of disturbances; see below.

A6 Disturbance and landscape heterogeneity
Landscape heterogeneity is accounted for by simulating a number of patches (within a grid cell) which differ by age (time since last disturbance).Each grid-cell state variable (e.g.stem biomass) is computed as a weighted mean of the patch state variables, with patch weight corresponding to the probability (p) of the patch occurring in the landscape.Disturbance is treated as a Poisson process, with time (x) since last disturbance distributed exponentially amongst the patches, according to the grid-cell mean disturbance interval (λ): p(x) = λ exp(−λx). (A13) In this work we held λ fixed in time, but it could also be prescribed a time-dependent variable.Disturbances (periodic events that recur randomly at the local scale, destroying all (catastrophic disturbance) or a fraction of biomass (partial disturbance) in a patch) significantly affect biomass residence time, vegetation structure and thereby resource use and productivity at a large spatial scale.At the level of a grid cell or large landscape, patches of different ages (years since last disturbance) should occur at frequencies corresponding to the expected likelihood of a local disturbance sometime in the corresponding period, given a known expected disturbance return time.The latter is prescribed in this work, but could be computed prognostically, e.g. by a wildfire module incorporated within CABLE.
A partial disturbance event in an affected patch is simulated by removing a size-dependent fraction of the biomass (and hence stems) in each cohort.In this work this fraction is specified using an observation-based relationship between tree size and the probability of tree survival following a fire of specified intensity (see Appendix A3).Following either type of disturbance event (catastrophic or partial), recruitment occurs and the patch persists in the landscape, and its age since the particular type of disturbance is set to zero.
Importantly, patch weightings by age are evaluated after growth, recruitment, non-disturbance-related mortality and disturbance.It is particularly important that the weightings be calculated after disturbance: otherwise patches of age zero would not be represented.
Vegetation structure changes more from year to year early in vegetation development following a disturbance than later.After a period corresponding to a few average tree generation times, a steady state is reached with no further net changes in vegetation structure (provided there are no trends in forcing data).To strike a balance between computational efficiency and "accuracy" in characterising landscape structure, we simulate patches with a sequence of (n age_max = 5 − 7) maximum ages with equally spaced cumulative probabilities (assuming an exponential distribution of patch ages, with expected value equal to the mean disturbance interval).To ensure a wide spread of ages in any given year, we additionally simulate (n patch_reps = 4 − 7) replicate patches for each maximum age, with the first disturbance occurring in year 1a, 2a, 3a...n patch_reps a, and thereafter every "n patch_reps a" years (where "n patch_reps a" is the maximum age of the patch).

A7 Patch weightings
Each patch is characterised by time since last disturbance, and its weighting in the landscape is given by the probability of this age occurring in an exponential distribution.When two disturbance types (catastrophic and partial) are considered, each patch is characterised by two ages corresponding to times since each disturbance type and two weights, one for each age.The patch weighting is then the product of the weights for each age, divided by the number of patches with the same combination of two ages, and normalised such that patch weights sum to one.
We evaluate weights for each unique time since each disturbance as follows: first we construct an ordered list of n age unique ages since last disturbance.Each unique age a i is associated with lower and upper integer bounds (b l,i and b u,i ), spanning the range of ages to be represented by a i .These bounds are set as follows: a i = 0 a i , (i = 1, a i > 0) or (i > 1, a i−1 = a i − 1) int[(a i + a i+1 )/2], 1 < i < n age a i , i = n age .
The weighting for each age a i is evaluated as the sum of the exponential frequencies, w = λ exp (−λx) (where λ is the expectation value of the time since disturbance and x the integral age), of integral ages from b l,i to b u,I inclusive.

A8 Tree foliage projective cover
Tree foliage projective cover is calculated following (Haverd et al., 2012) as 1−P gap , where P gap is defined here the probability of radiation penetrating the entire canopy from directly above (i.e.zero zenith angle):

Figure 1 .
Figure 1.Coupling of CABLE and POP, along with key inputs and outputs.
Figure 3. POP calibration.Biomass-density plot, showing all the points in the C-U data; a linear fit to all the C-U data; a linear fit to 30 data points lying along the upper bound; POP simulations of patches with the same age and StemNPP as the 30 data points lying along the upper bound.

Figure 4 .
Figure 4. Evaluation of CABLE-POP predictions against Cannell-Usoltsev data, separated into broadleaf and needleleaf classes respectively:(i) and (ii) biomass-density plot; (iii) and (iv) leaf-stem allometry plot (M fol and M stem are average foliage and stem biomass per tree (kg DM tree −1 )), including results derived for other LSMs byWolf et al. (2011); and (v) and (vi) total biomass: predictions vs. observations.In (i-iv), lines denote linear fits to the observations and predictions.Solid lines denote linear regression fits to the data points (see Table1for regression coefficients).

Figure 5 .
Figure 5. Biomass component fractions as a function of total biomass: top: CABLE-POP patches; second row: default CABLE; third and bottom rows: data-derived estimates reproduced from Wolf et al. (2011).

Figure 6 .
Figure 6.POP dynamics, as illustrated by undisturbed patch simulations at low and high extremes of annual stem increment (StemNPP = 0.05 kg Cm −2 year −1 and StemNPP = 0.20 kg Cm −2 year −1 ).(i) Biomass-density plot, showing the ageing trajectories of each patch (every 5th year plotted).Linear fits to the C-U data and their upper bound are shown for reference: (ii) time course of crown projective cover for each patch; and (iii) and (iv) time course of mortality, its components and turnover rate coefficient for patches with low (iii) and high (iv) annual stem increments.Dotted lines (blue) represent mortality based on a fixed turnover rate (dotted black line) of 0.02 year −1 , and (v) and (vi) time course of stem biomass for patches with low (v) and high (vi) annual stem increments.Solid lines indicate contributions form different cohorts.Dotted lines represent biomass accumulation assuming a fixed turnover rate of 0.02 year −1 , and (vii) and (viii) time course of stem density for patches with low (vii) and high (viii) annual stem increments.Solid lines indicate contributions from different cohorts.

Figure 7 .
Figure 7. Leaf-Stem allometry plots for (i) broadleaf and (ii) needleleaf PFTs.C-U data and CABLE-POP patch results are the same as in Fig.3(iii) and (iv).CABLE-POP grid-cell results are from deducing stem density from total grid-cell biomass (average over multiple patches, weighted by probability of time since last disturbance with a mean disturbance interval of 100 year) via Eq.(8), and using this to compute M fol and M stem from grid-cell foliage and stem biomass components.

Table 2 .
Reduced major axis regression coefficients associated with biomass-density, leaf-stem allometry and stem biomass plots of Fig.4.
* Standard linear regression, with intercept forced through zero.