Scaling from Individual Trees to Forests in an Earth System Modeling Framework Using a Mathematically Tractable Model of Height-structured Competition

The long-term and large-scale dynamics of ecosystems are in large part determined by the performances of individual plants in competition with one another for light, water, and nutrients. Woody biomass, a pool of carbon (C) larger than 50 % of atmospheric CO 2 , exists because of height-structured competition for light. However, most of the current Earth system models that predict climate change and C cycle feedbacks lack both a mechanistic formulation for height-structured competition for light and an explicit scaling from individual plants to the globe. In this study, we incorporate height-structured competition for light, competition for water, and explicit scaling from individuals to ecosystems into the land model version 3 (LM3) currently used in the Earth system models developed by the Geophysical Fluid Dynamics Laboratory (GFDL). The height-structured formulation is based on the perfect plasticity approximation (PPA), which has been shown to accurately scale from individual-level plant competition for light, water, and nutrients to the dynamics of whole communities. Because of the tractability of the PPA, the coupled LM3-PPA model is able to include a large number of phenomena across a range of spatial and temporal scales and still retain computational tractability, as well as close linkages to mathematically tractable forms of the model. We test a range of predictions against data from temperate broadleaved forests in the northern USA. The results show the model predictions agree with diurnal and annual C fluxes, growth rates of individual trees in the canopy and understory, tree size distributions, and species-level population dynamics during succession. We also show how the competitively optimal allocation strategy – the strategy that can competitively exclude all others – shifts as a function of the atmospheric CO 2 concentration. This strategy is referred to as an evolutionarily stable strategy (ESS) in the ecological literature and is typically not the same as a productivity-or growth-maximizing strategy. Model simulations predict that C sinks caused by CO 2 fertilization in forests limited by light and water will be down-regulated if allocation tracks changes in the competitive optimum. The implementation of the model in this paper is for temperate broadleaved forest trees, but the formulation of the model is general. It can be expanded to include other growth forms and physiologies simply by altering parameter values.


Introduction
Terrestrial ecosystems regulate biophysical exchanges of matter, energy, and momentum between the atmosphere and land surface and affect long-term climate dynamics by regulating the atmospheric CO 2 concentration (Chapin et al., 2008).Biogeochemical and biophysical interactions between terrestrial ecosystems and climate are now widely recognized as essential determinants of past and future climate change (Bonan, 2008).For this reason, global models of terrestrial Published by Copernicus Publications on behalf of the European Geosciences Union.

E. S. Weng et al.: LM3-PPA model
ecosystems are critical, but highly uncertain, components of Earth system models (ESMs) that predict climate and climate change (Friedlingstein et al., 2006).
In most ESMs, terrestrial vegetation is simulated by a dynamic global vegetation model (DGVM; e.g., Foley et al., 1996;Sitch et al., 2003) with global plant functional diversity represented by ∼ 10 plant functional types (PFTs; from Prentice et al., 1992).Vegetation in each model grid cell (e.g., 1 • latitude × 1 • longitude) is modeled as a set of pools describing different plant tissues (e.g., leaves, fine roots, sapwood, heart wood) belonging to one or more PFT (e.g., Sitch et al., 2008;Quillet et al., 2010).Mechanistic physiological and biophysical equations govern photosynthetic carbon gain, transpiration, respiration of all plant tissues, and uptake of water and (in some models) nutrients by fine roots.Modelspecific rules (often empirically derived) are used to allocate C to the different pools and to determine which PFT(s) dominates each grid cell or subgrid tile (Potter et al., 1993;Foley et al., 1996;Sitch et al., 2003).Dead plant tissues are sent to a decomposition submodel, which usually is a variant of the CENTURY model (Parton et al., 1987).Water availability is governed by a coupled hydrological model.Some DGVMs include dynamical models of important nutrients, such as nitrogen (Thornton et al., 2007;Zaehle and Friend, 2010;Gerber et al., 2010).In fully coupled implementations, plant canopies exchange carbon, water, energy, and momentum with the atmosphere through a boundary layer above the canopy airspace, and roots exchange matter and energy with one or more soil layers.
Although we have a sophisticated understanding of some important fine-scale processes, such as leaf-level photosynthesis, and a growing capacity to measure grid-scale fluxes and storage of carbon, most current ESMs lack a set of equations that explicitly scale physiological, population dynamic, and biogeochemical processes from individual plants to stands, communities, and grid cells.This may contribute to the high uncertainty about C sources and sinks predicted by the ESMs as revealed by model intercomparison studies (Shao et al., 2013;Todd-Brown et al., 2013;Friedlingstein et al., 2014).For example, some models predict that CO 2 fertilization and climate change will create a large terrestrial C sink, whereas others predict a large C source, with the spread between models large relative to global anthropogenic fossil fuel emissions (Friedlingstein et al., 2006).
Several DGVMs with explicit scaling have been developed from forest gap models (Friend et al., 1997;Sato et al., 2007;Haverd et al., 2014), which have been shown to scale from individual vital rates to stand dynamics with reasonable accuracy (Botkin et al., 1972;Shugart and West, 1977;Pacala et al., 1996), and are thus widely used to manage forests (e.g., Coates et al., 2003).Some gap models simulate heightstructured competition among individual seedlings, saplings, and adult trees for light, as well as competition for belowground resources.Because simulating every individual plant on Earth in this way is unfeasible, some models, such as HYBRID (Friend et al., 1993(Friend et al., , 1997)), LPJ-GUESS (Smith et al., 2001), and SEIB (Sato et al., 2007), simulate a sample of individuals in each grid cell that is small enough to allow reasonable run time but large enough to dampen random fluctuations in the underlying stochastic population dynamics.An alternative approach was developed by Moorcroft et al. (2001), who derived a set of integro-partial differential equations that approximately govern the dynamics of the first moment of the stochastic process (the mean population density of trees in the forest of each species and size) that is simulated in a gap model.Instead of averaging over the many individuals in a stochastic simulation, these equations directly predict the mean population densities of individuals of each species and size (height, diameter, or biomass) that would have been produced by a gap model of a large stand with the same functional forms and parameter values.Medvigy et al. (2009) and Fisher et al. (2010) coupled the ED model into full DGVMs, and several efforts are now underway to build models derived from ED into ESMs.
An important advantage of the DGVMs developed from gap models, such as HYBRID, LPJ-GUESS, SEIB, and ED, is that they include the mechanistic function of stem wood.Trees use stem wood to overtop their neighbors when in competition for light, and to avoid being overtopped by their neighbors.The wood of living trees is the largest vegetation carbon pool (363 ± 28 Pg C; Pan et al., 2011), equivalent to around half of the atmospheric carbon pool.Furthermore, a large fraction of soil organic matter (SOM) comes from wood litter.It is thus likely that predictions about the future of the terrestrial C sink will be improved in models that include the mechanistic function of wood.For example, to determine how the terrestrial C sink will change because of climate change and CO 2 fertilization, one needs to predict changes in plant C allocation patterns.Because of the large difference in residence time of wood, leaves, and fine roots in forests, changes in allocation can drastically change carbon sinks (Luo et al., 2003;Zhang et al., 2010).Theoretically, it has been shown that, under water limitation, competitively optimal shifts towards greater fine-root allocation can lead to greatly diminished vegetation C sinks despite significant increases in productivity (Farrior et al., 2013).Thus, mechanistic predictions of whether allocation to wood will increase, decrease, or stay the same under the altered environmental conditions are critical.However, competitively optimal plant allocation has not, to our knowledge, been rigorously studied in any of the previous gap-model-derived DGVMs.
Despite the advantages of gap-model DGVMs, it is difficult to understand the behavior of these models because they are analytically intractable even under idealized conditions, such as constant climate, and so can only be studied using numerical simulations.For example, competitively optimal plant C allocation could only be studied in these models by relying on computational experiments that may be difficult to interpret in the absence of any theoretical guidance.The price of added complexity in a DGVM is that it increases the number of ways in which model errors can interact and cause misleading predictions (e.g., model equifinality), which are especially difficult to diagnose and understand if one cannot study the model analytically.This problem is particularly acute when developing an ESM, which has many interacting components.For this reason, height-structured competition was not included in the GFDL land model version 3 (LM3) (Shevliakova et al., 2009;Milly et al., 2014).
In this paper, we present a new, biodiverse version of LM3 that includes height-structured competition among plants for light, as well as competition for water.Future versions will include competition for nitrogen and phosphorus.The new model, LM3-PPA, is based on the perfect plasticity approximation (PPA), a computationally simple and mathematically tractable model that scales from individuals to stand dynamics (Strigul et al., 2008).Like ED, the PPA allows one to derive integro-partial differential equations for the first moment of the stochastic process that defines an individual-based forest model (Strigul et al., 2008).But, unlike ED, these equations are analytically tractable under idealized conditions (e.g., constant climate).The PPA model closely matches the behavior of stochastic individual-based forest dynamics models (gap simulators; Strigul et al., 2008).More importantly, it has been shown to predict species-level succession across different soils in the US lake states (Purves et al., 2008) and to accurately predict canopy structure in temperate and tropical forests (Bohlman and Pacala, 2012;Purves et al., 2007;Zhang et al., 2014).Dybzinski et al. (2011Dybzinski et al. ( , 2013) ) and Farrior et al. (2013) have developed game-theoretic versions of the PPA that use analytical methods to identify the most competitive allocation strategy (investment in fine roots, wood, and leaves) of trees competing for light, water, and nitrogen.Although these game-theoretic models are physiologically simpler than most DGVMs, they yield quantitatively accurate predictions of net primary production (NPP) and plant allocation observed at Fluxnet sites (Luyssaert et al., 2007).These theoretical studies have guided the development of the new DGVM presented here, LM3-PPA.
Although the fast-timescale processes in LM3-PPA (e.g., exchanges of energy and matter between vegetation, atmosphere, and soil) render it analytically intractable, its close association with the stand-alone PPA model allows for a greater understanding of model behavior than is possible with other gap-model DGVMs, including how competition for multiple resources is expected to affect allocation of NPP among different plant tissues.Variation among individuals, species, or PFTs in how carbon is allocated to leaves, wood, fine roots, etc. is recognized as a key feature of nextgeneration DGVMs that aim to represent plant functional diversity (both within and between model grid cells) more accurately than the current suite of models (Scheiter et al., 2013;Wullschleger et al., 2014).LM3-PPA was specifically designed with allocational and other aspects of plant functional diversity in mind.
In particular, we developed LM3-PPA to 1. include the influence of height-structured competition for light on forest dynamics and dominant allocation strategies; 2. improve the representation of feedbacks that alter ecosystem-level allocation to wood; 3. include within-PFT biodiversity by allowing for multiple, competing variants or "species" that differ in their allocational strategy or other traits; 4. improve the scaling from individuals to landscapes using macroscopic equations from the literature on individual-based forest models; 5. provide a global land model that can be solved analytically in idealized cases (e.g., constant climate).
In what follows, we first present the equations that underpin the LM3-PPA model in their continuous (in time and plant size) form.The numerical machinery that is necessary to discretize and implement the model as a component of an ESM is described in the Appendix (A and B).The model structure allows for an arbitrary number of "species" (broadly defined to include different genotypes or PFTs) that may have fixed or plastic parameter values describing their physiological properties and how they allocate available carbon.We evaluate the model's behavior at a series of organizational scales in a temperate forest -physiological (photosynthetic carbon gain), individual (stem diameter and height growth rates), population (size structure and population densities), community (species-level successional dynamics), and ecosystem (C storage, NPP) -and at a series of temporal scales: diurnal, seasonal, interannual, and centennial.We also introduce a prototype algorithm for determining the most competitive allocation strategy (i.e., the evolutionarily stable strategy, ESS) within a functional type.We use this ESS algorithm to evaluate the expected shift in C allocation between fine roots and woody tissues caused by the leaf-level water use efficiency benefits of CO 2 fertilization and the impact of this shift on the predicted C sink.

E. S. Weng et al.: LM3-PPA model
model predicts the size-structured dynamics of each species (or PFT, for example) by predicting the fate of each and every individual.This spatial stochastic process is analogous to the dynamics of the atmosphere resulting from the stochastic movement of every gas molecule.In the same way that one can derive the Navier-Stokes equations from the stochastic process of molecular motion, it is possible to derive equations for the mean population densities of trees from a stochastic gap model.But because gap models are highly nonlinear, approximations must be used.One impediment to a tractable approximation has been the lack of a mechanistic and compact way of representing how the irregular spatial distribution of stems, which strongly affects the outcome of competition, results in a nearly continuous leaf canopy, which strongly affects gas exchange.ED, like the stochastic models from which it was derived (Shugart, 1984;Botkin et al., 1972), does so by simply partitioning space into adult-treesized cells and assuming that each individual's crown covers all the area in its cell (Moorcroft et al., 2001).As there can be many individuals per cell, there can be many overlapping canopies, and any tree that is not the tallest in the cell is not in full sun.The problem with this assumption is particularly evident in recently disturbed gaps, which in reality may contain multiple trees that are all in full sun.
In the forest gap simulator from which the PPA is derived (SORTIE with plastic crown shapes due to phototropism; Strigul et al., 2008), the crown of an open-grown tree of height Z is an envelope of leaves whose shape is defined by a function A(Z , Z) that gives the crown area above height Z .The potential crowns of trees in a closed-canopy forest overlap, so that, from above, the canopy looks like a patchwork of non-overlapping territories, with each territory being the portion of a canopy tree's crown that is in the sun (Mitchell, 1969).Strigul et al. (2008) studied the statistics of the places where the potential crowns of adjacent canopy trees join and showed that, if tree growth was realistically plastic because of mild phototropism of apical meristems, then the standard deviation of these canopy-crown join heights is an order of magnitude smaller than the mean join height.They derived approximate equations, taken in the limit of zero crown-joinheight standard deviation, for the time evolution of the first moments of the stochastic process in the gap simulator, i.e., the function N i (s, t) for each species i in the model, which gives the expectation of a species' population density for individuals of size s at time t.The derivation of these means used only the individual-level information in a gap simulator and thus scales from individual to stand.The approximation is called the perfect plasticity approximation (PPA) because it is derived from the limit of extreme flexibility of crown shape in the horizontal by trees in pursuit of light.
The PPA equations are a special case of a general sizestructured demographic model governing the time evolution of the population density of individuals of species i and size s, N i (s, t) (Strigul et al., 2008;von Foerster, 1959).One should think of N i (s, t) as the mean population density of individuals per unit ground area in a stochastic gap model.It is the limit of the expectation of n ist /(δ s δ x δ x ) as δ s and δ x approach zero, where n ist is the number of individuals of species i with size between s and s + δ s at time t in a randomly chosen quadrat with ground area (δ x ) 2 in runs of a stochastic individual-based forest model (Strigul et al., 2008).In reality, and in most of this paper, the size (s) of a tree is a vector describing its height, crown area, tissue pool sizes, etc.But for the moment, consider the simple case where there is only a single measure of size.The system of equations governing the time evolution of N i (s,t) is usually written as a system of nonlinear advection equations (advection in s) with a boundary condition governing the recruitment of new individuals at the smallest size (Moorcroft et al., 2001;Strigul et al., 2008).But we write them here in a mathematically equivalent form as implemented in the LM3-PPA code.

Population dynamics
LM3-PPA makes population dynamics predictions by simply simulating the birth, mortality, and growth of each age cohort of plants.The cohorts within the same place (tile within a grid cell; see Appendix A) interact with one another only indirectly by affecting resources levels -canopy trees shade understory trees, and all cohorts reduce available water.In addition, cohorts in the same place have indirect biophysical impacts on one another because they jointly affect the temperature and humidity of the sub-canopy airspace.These indirect effects are explained in later sections and a series of appendices.Here, we describe the population dynamics assuming that the resource levels and biophysical conditions affecting a cohort are known.For each species (or PFT) i, we describe its population density in the same age cohort N ib (t), individuals of species i at time t that were produced at time b, and the size of the individuals in this cohort (s ib (t)).If F i (s,t) is the rate of new seedling production for an individual of species i at time t and size s, then the population density of a newly produced cohort is (1) Eq. ( 1) simply sums the reproductive output of all cohorts of a given species to get the initial density of the new cohort produced at time t.We also need an equation for the loss of individuals in each cohort as it ages.After being produced, individuals die at rate µ i (s,t): Finally, we need an equation for the growth of individuals in each cohort.If g i (s,t) is the growth rate of individuals of species i and size s at time t, then ds ib (t) dt = g i s ib (t), t . (3) Equations (1-3) provide an efficient way to solve the model numerically, because one can simply discretize b and thus yield a set of ordinary differential equations that have much greater numerical stability than advection equations.The LM3-PPA model uses this numerical method and thus simulates a discrete number of cohorts.
To convert Eqs.(1-3) into the measures we need for a DGVM, we first divide each individual into five separate tissues or carbon (C) pools (leaf, fine root, sapwood, heartwood, and nonstructural carbohydrates; Fig. S1 in the Supplement) and introduce allometric relationships to calculate the amounts of C in these five pools, as well as other measures of size, from three quantities: stem diameter (D(t)), crown LAI (l(t)), and carbon in the nonstructural carbohydrate pool (NSC(t)).Stem diameter and crown LAI were chosen because these are easily observable, and NSC(t) because all plant carbon starts as nonstructural carbohydrates.In this paper, stem diameter is assumed to equal diameter at breast height (DBH) in any comparisons with DBH data.With the three measures of size, s, the right-hand sides of Eqs.(1-3) each become three separate equations -one for each measure of size.Also, because each cohort has a size vector, it is always possible to calculate the density of a species or PFT as a function of any measure of size, rather than as a function of birth date.In what follows, we switch to size-structured densities, N i (s,t), whenever convenient.

Vertical and horizontal spatial structure
Again, each cohort in LM3-PPA belongs to a species (or PFT, for example) and has three time-evolving measures of size: stem diameter, D(t); crown LAI, l(t); and amount of carbon in the nonstructural carbohydrate pool, NSC(t).We sometimes omit from the notation the time dependence from D(t), l(t), and NSC(t) to keep the formulae easy to read.These measures are related to other important measures of size by species-or PFT-specific allometric relationships.Height, Z(D); wood carbon mass, S(D) (including stem, branches, and coarse roots); and crown area, A CR (D), are functions of diameter: where α c , α Z , , and ρ W (wood carbon density; kg C m −3 ) are species-or PFT-specific constants; θ c and θ Z are constant across species/PFTs (1.5 and 0.5, respectively), though these could be made species/PFT-specific if necessary.A tree's total leaf mass, L(D, l), is its total leaf area (l × A CR ) times its species-or PFT-specific leaf mass per area, LMA, and -following the pipe model (Shinozaki et al., 1964) -fine-root carbon mass, FR(D,l), and sapwood cross-sectional area, A sw (D,L), are proportional to total leaf area: where ϕ RL , SRA, and α CSA are species/PFT-specific constants: ϕ RL is the ratio of total root surface area to the total leaf area, SRA is specific root area, and α CSA is an empirical ratio of target leaf area to sapwood cross-sectional area.Unless otherwise stated, units are mass in kilograms carbons, area in square meters, height in meters, and diameter in meters.All other size measures of structural pools can be calculated from these quantities.For example, heartwood carbon mass is S(D)-A sw (D,l)Z(D) ρ W.

Fine-root spatial structure
Because the area covered by a tree's root distribution is significantly larger than its crown area (Hruska et al., 1999), we assume that roots of competing individuals are uniformly distributed in the horizontal plane (Dybzinski et al., 2011 and refs therein).LM3 and LM3-PPA can be configured with an arbitrary number of vertical soil layers, with 20 layers in this study (see Appendix B for details).Each species or PFT has an empirical exponential depth distribution for its fine roots (Appendix B).

Canopy structure
A critical quantity in the PPA model is the crown join height that separates the upper canopy from the understory Z * 1 : where k =1 for the top canopy layer, η is the proportion of each canopy layer that remains open due to spacing between individual tree crowns, N i (Z, t) is the density (m −2 ) of trees of species i with height Z, and A i (Z * k , Z) is the area (m 2 ) of the portion of a tree's crown at a height greater than or equal to Z * k .If the right-hand side (RHS) of Eq. ( 6), the collective crown area of all trees per unit ground area, is less than the fraction of ground area that could potentially be filled (1-η) even for Z * 1 = 0, then plant density is too low to close the canopy.However, if the Z * 1 that solves Eq. ( 6) is greater than zero, then the trees close the canopy, by definition filling the canopy with the sun-exposed portion of the crowns of individuals taller than Z * 1 .Plants that are shorter than this value, Z * 1 , are in the understory.In many temperate and boreal forests, the potential crowns of all individuals add up to less than two (do not fill a second canopy), and so Eq. ( 6) has

E. S. Weng et al.: LM3-PPA model
no solution for k > 1.However, in some forests (e.g., tropical rainforests, and temperate forests with multiple understory layers), the sum of the crown areas of all individuals combined is typically 3 to 4 times the land area (Bohlman and Pacala, 2012), in which case Eq. ( 6) defines a Z * 2 separating the first full understory from the second understory beneath it, a Z * 3 separating the second from the third understory, and so on.
Mathematical and computational tractability is greatly facilitated in the PPA model by the assumption that trees have flat-topped crowns (Strigul et al., 2008), which allows for accurate predictions of observed succession and canopy structure in broad-leaved temperate forests (Purves et al., 2008;Zhang et al., 2014) and vital rates and canopy structure in a neotropical forest (Bohlman and Pacala, 2012).With a flattopped crown, all the leaves of a tree are assumed to be in one layer, either in the upper canopy or in a single understory layer (Figs.A1 and S1a).We assume flat tops in LM3-PPA and thus use A CR (D) as the sole measure of crown area, i.e., A i (Z',Z(D)) = A CR (D) for all Z ≤ Z(D).Each cohort in LM3-PPA (and all of its leaves) belongs to exactly one canopy layer.Again, the upper canopy layer includes the tallest cohorts of trees whose collective crown area sums to the fillable ground area (1-η, times the ground area; or less than this area if the canopy is not closed; Eq. 6).Trees within the same layer do not shade each other.The trees in each understory layer are shaded by the leaves of all taller canopy layers (Appendix B).In LM3-PPA, the assumption of flattopped crowns introduces a potential problem that does not occur in simpler versions of the PPA model that lack physiological mechanisms.Specifically, the NSC pool can, in some cases, be quickly consumed when a tree enters the upper canopy from the understory because of the sudden increase in target leaf and fine-root biomasses.This increase would be more gradual with other crown shapes (e.g., rounded).To address this problem (which we view as a model artifact), we introduced a parameter to limit the rate of increase of target leaf mass (and therefore fine-root mass, given the pipe-model constraint) for cohorts that recently entered the upper canopy (see Eq. A6 in Appendix A).

Fast-timescale exchanges of matter, energy, and momentum
Like other land models that are fully coupled to atmospheric models, LM3-PPA computes fluxes of matter, energy, and momentum between a plant's surface and the bottom of the boundary layer in the atmosphere on the fast timescale of the atmospheric model (e.g., every 30 min in most implementations of LM3 and LM3-PPA).This requires a network of interacting equations that are similar among many land models, including 1. energy and mass balance equations that govern leaf, canopy air, and soil temperatures; canopy vapor pres-sure deficit (VPD); wind speed in the canopy airspace; and long-and shortwave radiation transfer; 2. a photosynthesis model at the leaf level; 3. a model of respiration for all plant tissues; 4. a model of stomatal conductance and fine-root water uptake; 5. a model of soil water dynamics; 6. a model of the decomposition of soil organic matter.
The fast-timescale equations are described in Appendix B.
They are identical to those in the version of LM3 used in the ESMs of GFDL (Dunne et al., 2012(Dunne et al., , 2013)), except for a few key differences.First, whereas LM3 has only a single cohort in any one place, LM3-PPA has a multi-cohort canopy and fine-root distribution that (a) can be composed of more than one species/PFT, (b) may have one or more complete understory canopy, and (c) always has a partially full lowest understory layer if it has one or more full canopies.Second, the respiration parameterization for sapwood has been updated in LM3-PPA.Observations show that the respiration rate of sapwood per unit of biomass decreases with sapwood biomass (Ryan et al., 2004).Consistent with these observations, LM3-PPA assumes that respiration of sapwood is proportional to crown area, A CR (D).LM3-PPA handles radiation transfer through the crowns of each cohort in the same way that LM3 handles transfer through its single canopy.Radiation emanating from the bottoms of crowns in the same canopy or partial canopy layer is summed before hitting the next layer or the ground.All other calculations are made separately for each cohort, and summed where necessary.For example, sensible and latent heat fluxes from the leaves of each cohort into the sub-canopy airspace are summed in the energy balance for the airspace.Appendix B documents the details of the fast-timescale calculations in LM3-PPA.

Growth and reproduction
In this section, we briefly describe the fecundity function (F ) in Eq. ( 1) and the growth functions on the RHS of Eq. (3) for stem diameter, D(t); crown LAI, l(t); and amount of carbon in the nonstructural carbohydrate pool, NSC(t).The derivations and detailed discussion of these expressions are in Appendix A.
The carbon fluxes from the fast-timescale equations (Appendix B) are summed over the diurnal cycle to provide daily total carbon gain from photosynthesis, P s (t), and loss from respiration, R a (t), for each cohort.This carbon is added to or taken from the cohort's NSC pool once a day: where G L+FR (t) is the amount of carbon allocated to produce new leaves and fine roots minus the carbon retranslocated from senescing leaves and fine roots, and G W+F (t) is the carbon allocated to stem and seed production.Expressions for G L+FR (t) and G W+F (t) are derived in Appendix A. The derivations assume that a plant allocates carbon so that the LAI within its crown tracks a species-or PFT-specific target.This target crown LAI differs between understory and canopy individuals and seasonally because of a phenology function, p(t), which is unchanged from LM3, except that it is updated daily rather than once per month as in LM3 and LM3V (Shevliakova et al., 2009;Milly et al., 2014).Individuals also have a target root area per unit crown area, which is equal to the target crown LAI multiplied by ϕ RL (the ratio of total root surface area to total leaf area; see Eq. 5).Finally, there is a target ratio of wood to seed production, and a species-or PFT-specific NSC target, which scales with target leaf mass and tracks a plant's phenological state (Eq.A5 in Appendix A).
Our formulation for G L+FR (t) assumes that positive net production, P s (t) − R a (t), is allocated first to leaves and fine roots if these are beneath their target levels.Carbon is retranslocated back to NSC if leaves are above target (i.e., at the end of the growing season, or if a cohort falls into the understory from the overstory).Carbon is allocated to wood and seeds from NSC only if NSC is above its target level.The formulation also includes parameters that limit the maximum rate at which NSC can be converted into leaves and fine roots and wood.
Appendix A shows how the assumptions about allocation can be combined with the allometric equations (Eqs.4-5) to produce differential equations for the growth of stem diameter and crown LAI.All other measures of plant size (e.g., fine-root mass or leaf mass) can be calculated from NSC, diameter, and crown LAI using the allometric equations.

Mortality and disturbance
In this section, we specify the mortality functions on the RHS of Eq. (2).Mortality in the PPA reduces the population density of a cohort (i.e., by a fraction µδt in a time step δt if the individual mortality rate is µ).In LM3-PPA, mortality is assumed to occur due to carbon starvation if a cohort's NSC pool falls to zero.Because the target size of the NSC pool is assumed to be several times the size of the combined target leaf and fine-root masses (see Eq. A5 in Appendix A), trees rarely die of carbon starvation unless they experience prolonged drought (which was not simulated in the current study) or have chronic negative carbon balance due to shading.In addition to carbon starvation, each species/PFT has a canopy-layer-specific background mortality rate that is assigned from the literature (Runkle, 2000).These background rates are assumed to be size-independent for upper-canopy trees (µ C0 in Table 1) but size-dependent for understory trees  2).In panel (a), the solid line shows the mortality rates of understory trees (same for H0-H3); the dashed line shows the mortality rates of canopy trees in H0, H2, and H3; and the dotted line is for canopy trees in H1.In panel (b), the solid line is tree height (same for H0-H3); the dashed line shows crown area in H0, H1, and H3; and the dotted line is crown area in H2.
according to This functional form reduces mortality by a factor of 3.67 between germination and adulthood (Fig. 1a).It accounts for the additional sources of non-starvation mortality facing small individuals, including herbivory by large mammals and branch fall.For all canopy layers, the background mortality rate is assumed to be independent of the physiological state of the focal individual and the density of competing individuals, as these physiological and competitive effects are already accounted for by mortality due to carbon starvation.We also evaluated an alternative assumption for canopy trees in this paper, in which the mortality rate of large trees increases with size (see Sect. 2.2.2 below).
Stand-level disturbances (e.g., due to insect outbreaks, windstorms, fire, or land use) may be implemented in LM3-PPA using the land use tiling scheme described below and in Appendix A, but they were not implemented in the simulations presented in this paper.for trees with DBH < 0.8 and ≥ 0.8 m, respectively.Dµ C = 1.0 m The "Gap dynamics" algorithm labeled "Tallest" is the standard PPA assumption, in which the tallest understory trees fill the space vacated by the death of canopy trees (H0-H2; Strigul et al., 2008).The alternative assumption ("Randomly selected") selects understory trees at random (regardless of their height) to fill this vacated space (H3).
Grid structure, subgrid-scale heterogeneity, and relation to LM3 Like LM3, LM3-PPA is implemented on a flexible grid, whose cell size can be specified independently of the atmospheric model's grid.LM3-PPA also includes LM3's dynamic tiling scheme for land use, stand-level disturbance, and subgrid-scale heterogeneity (Shevliakova et al., 2009).
As explained in Appendix A, the tiling scheme can be used to implement the ED approximation for canopy gap dynamics (Moorcroft et al., 2001), but this feature was not used in the simulations presented in the current paper.The critical difference between LM3 and the LM3-PPA model described in this paper is that each tile in LM3-PPA can contain an arbitrary number of cohorts that compete with one another for light and water.Each cohort belongs to a single species or PFT, but different cohorts within the same tile can be from different species/PFTs.Thus, there is competition for light and water among cohorts belonging to the same species/PFT (intraspecific competition), as well as among cohorts belonging to different cooccurring species/PFTs (interspecific competition).Coexistence of multiple species/PFTs is not assumed but is rather a possible emergent outcome of the individual-level processes that determine the community dynamics.

Model evaluation and simulation tests
The model was evaluated in temperate deciduous forest in Wisconsin, USA.A variety of data are available in this region to evaluate the model's behavior, including forest inventory data from the U.S. Forest Inventory and Analy-sis (FIA) database (http://www.fia.fs.fed.us/),biometric data (Curtis et al., 2002), and eddy-covariance data (Desai et al., 2005).Furthermore, there are clear patterns of forest succession among some of the dominant tree species in the region (see below), which facilitates tests of predicted successional dynamics.Meteorological inputs were extracted from the Sheffield et al. (2006) 1 • latitude × 1 • longitude, 3-hourly,   climate reanalysis data set for the grid cell containing the Willow Creek AmeriFlux site (Desai et al., 2005).We forced the model with the Sheffield reanalysis data rather than the meteorological data from the AmeriFlux site because some model tests (e.g., forest size structure and successional chronosequences) were performed at a regional scale (see details below).
Models such as LM3-PPA are inevitably tuned during development so that they reproduce realistic behavior.We tuned physiological aspects of the model (photosynthesis, respiration, and NSC dynamics) to produce the observed magnitude of NPP and a single parameter affecting diameter growth rates (the taper constant, ).We also tuned the size dependence of background mortality (Fig. 1a) for small seedlings and saplings to reconcile large observed abundances of germinating seedlings with low observed abundances of saplings.We did not tune emergent behaviors such as differences among the growth rates of canopy and understory trees, differences among the growth rates of trees of different species, population densities of individuals above 0.1 m in diameter, successional turnover, and patterns of carbon storage.In what follows, comparisons of predicted and actual NPP should be viewed as demonstrations that the model is capable of exhibiting realistic behavior, because physiological aspects of the model were tuned.However, comparisons involving variation among individuals in whole-tree growth rates, population densities, and size structure for individuals above 0.1 m in diameter, and successional and ecosystem dynamics should be viewed as tests of emergent predictions of the model.

One-vs. three-species simulations
We implemented the model with three tree species -trembling aspen (Populus tremuloides Michx.), red maple (Acer rubrum L.), and sugar maple (Acer saccharum Marsh.)-to evaluate the model's capacity to capture successional dynamics and to quantify how successional diversity affects model behavior compared to one-species simulations.The three species are common in eastern North America and at the Willow Creek site in particular, and they differ in their successional status and shade tolerance (Burns and Barbara, 1990): trembling aspen is a pioneer species with a high growth and mortality rate and low shade tolerance, sugar maple is a latesuccessional species with a low growth and mortality rate and high shade tolerance, and red maple is an intermediate species.These three species are not intended to fully characterize the Willow Creek or other temperate tree commu-nities, and in this paper we do not attempt to determine the optimal number of species or functional types for ESM applications.In addition to the three-species simulations designed to evaluate successional dynamics and perform model-data comparisons at Willow Creek, we also performed a series of competition experiments with multiple functional variants defined by their allocational strategy (see Sect. 2.3,below) as an initial exploration of an axis of functional variation that could be incorporated into future global applications.We estimated model parameters for the three Willow Creek species using data from the literature (Table 1).Most of the other parameter values (Tables 1 and C1-C3 in Appendix C) were taken directly from LM3.
We compared carbon and population dynamics of runs with one species (sugar maple) and all three species.Simulations were initialized with a number of small seedlings for each species (Tables 1 and C4) and run for 1000 years.Runs simulating species succession were initialized with abundances and size distributions of each species from earlysuccessional FIA plots (plots less than 10 years of age, Table C4).We examined model predicted population densities; size distributions; annual gross primary production (GPP) and NPP; growth rates of diameter at breast height (DBH), foliage biomass, stems, and fine roots; and total C storage.
We compared model output both to published data of GPP, NPP, plant DBH growth rates, and forest composition at the Willow Creek AmeriFlux site and to FIA data on mesic soils from the Laurentian Mixed Forest Ecoprovince (Cleland et al., 2007), which spans northern Michigan, Wisconsin, and Minnesota, USA, and includes the Willow Creek site.Hereafter, we refer to this ecoprovince as the "northern lake states".Each FIA plot includes measurements on only 0.067 ha distributed over a 0.4 ha area; thus, data from many plots must be aggregated by stand age class to estimate successional patterns of biomass, density, and size distribution.

Sensitivity of LM3-PPA to alternative assumptions: mortality, allometry, and gap dynamics
Runs of LM3-PPA predict realistic size distributions for a few hundred years of succession but produce unrealistically large trees in old-growth forests (see results below).
Although there are only a few unrealistically large trees, they are so large that they store considerable carbon and skew predictions.We have encountered this problem before when working with forest gap simulators (e.g., SOR-TIE; Pacala et al., 1996), and we hypothesize two possible causes.First, although LM3-PPA assumes constant size-and density-independent death rates of canopy trees (aside from carbon starvation, which rarely occurred for canopy trees in the simulations presented here), many studies have documented increased mortality as trees become very large (Runkle, 2000).Xu et al. (2012)  in an old-growth temperate forest.We thus compared H0, the baseline LM3-PPA model with constant canopy tree mortality rates, with H1, the same model with upper-canopy mortality rates that increase with tree size as shown in Fig. 1a.Second, the allometry and respiration assumptions in LM3-PPA predict that a canopy tree's DBH growth rate increases monotonically to an asymptote as a tree becomes large.This prediction is supported by dendrochronological studies for the first one or two centuries, but actual growth rates subsequently decline in very old trees (Sillett et al., 2010).We compared output from H0 and H2, in which DBH growth rates decline for very old trees, as reported in dendrochronological studies.Rather than prescribing an arbitrary growth curve, the DBH growth rate decline results from a modified crown area allometry in H2, in which crown area becomes constant after a tree reaches 0.8 m in DBH (C.Canham, unpublished data), rather than continuing to increase with diameter according to the crown area allometry in H0 (see Eq. 4).The modified allometry in H2 results in declining DBH growth rate for DBH > 0.8 m because leaf area (and thus potential C gain) plateaus.All else equal, this causes sapwood volume growth to plateau, which causes decreasing diameter growth (because the volume is "stretched" around a growing circumference and along an increasing height).
Finally, the mathematical approximation behind the PPA leads to a sharp separation between canopy and understory, i.e., a single height at any one time separating all canopy trees from all understory trees in a given stand (or subgrid cell tile in LM3-PPA).The PPA thus predicts that old-growth recruitment into the canopy comes exclusively from saplings that have spent a long time in the understory (advance regeneration).While this is true for shade-tolerant species, it is not true for pioneers that exploit large gaps in old-growth forests.Section 5 of Appendix A describes how the subgridscale tiling scheme in LM3-PPA could be used to simulate gap dynamics (which were not implemented in the simulations presented in this paper).We suspect that this change will be necessary to maintain successional diversity indefinitely in old growth, but we do not expect that gap phase dynamics would substantially affect old-growth carbon storage because most trees in old growth belong to shade-tolerant species.To check this supposition, we compared runs of the baseline model with identical runs of H3 -a model in which understory cohorts were drawn at random (independent of size) to fill space in the canopy opened by canopy tree mortality.Comparisons between the three alternative models (H1-H3) and the baseline model (H0) were based on simulations with one species (sugar maple).

Comparison with a standard biogeochemical model
To explore how incorporating individual-level competition and successional diversity into land models affects carbon accumulation in vegetation and soil, we compared the LM3-PPA predictions to those of a CENTURY-like standard biogeochemical (BGC) model (Fig. S1b) as described in Parton et al. (1987) and Luo et al. (1999).Like most current DGVMs and land surface models, the standard BGC model that we implemented was formulated at the level of the grid cell without explicitly scaling from individual plants to ecosystemlevel dynamics.In such models, photosynthesis and respiration submodels simulate the net influx of C (NPP) at the level of the grid cell.NPP is then allocated to grid-cell-level plant C pools and, after senescence, plant carbon moves through litter and soil pools before returning to the atmosphere.Carbon allocation coefficients and residence times in the various pools determine total carbon storage (Weng and Luo, 2011).
We chose this BGC model because all of its C pools -leaves, fine roots, sapwood, heartwood, labile soil carbon, and recalcitrant soil carbon -can be precisely matched to quantities predicted by LM3-PPA.The BGC model simulations were forced with the NPP produced by the single-species runs of LM3-PPA, and so differed only in the patterns of allocation and residence times assumed in the standard BGC model and those that emerged by aggregating finer-scale patterns in LM3-PPA.

Competitive allocation strategies at different CO 2 concentrations
A competitively optimal allocation strategy is the one that can competitively exclude all others.This can be significantly different from the allocation strategy that most effectively uses available resources (i.e., the optimal monoculture strategy).The analytical model derived by Farrior et al., 2013) predicts that increased leaf-level water use efficiency from CO 2 fertilization should cause a shift in the competitively optimal allocation strategy among fine roots, leaves, and wood, which in turn causes the changes in carbon storage described in the Discussion section of this paper.We simulated competition among red maple variants with different target fine-root biomasses under each of two atmospheric CO 2 concentrations ([CO 2 ]) in LM3-PPA: 280 ppm for preindustrial and 560 ppm for doubled [CO 2 ].All runs shared the same meteorological forcing.All red maple variants shared all parameters except for the ratio of fine root to leaf surface area (ϕ RL ) for canopy individuals.Because the target crown LAI of a canopy tree (l * C ) was constant across red maple variants -and because the amount of carbon allocated to wood depends on the amount of NSC not taken by leaves and fine roots (see Appendix A) -variation in canopy tree ϕ RL among variants had little effect on leaf allocation but strong effects on fine-root and wood allocation.Across different monocultures that differ only in ϕ RL , fine-root allocation should increase and wood allocation should decrease with increasing ϕ RL , at least in the region of parameter space near the competitive optimum.Note that this fine-root vs. wood allocational tradeoff is not necessarily apparent when comparing allocational types in competition with each other.
For example, relatively high ϕ RL may offer a competitive advantage if trees are water-limited, which could increase carbon gain and fractional wood allocation compared to less competitive types with lower values of ϕ RL that have little NSC available for wood growth.
We performed three sets of experiments with different canopy tree variants with ϕ RL ranging from 0.5 to 1.0 (understory ϕ RL was 0.8 for all variants).Each experiment was performed at both preindustrial and doubled [CO 2 ] (Table 3): 1. Polyculture runs were initiated with five variants (ϕ RL = 0.5, 0.6, 0.7, 0.8, and 0.9) all having the same initial population density (250 seedlings ha −1 ).Polyculture runs simulated competition among the five variants for 500 years to identify the most competitive strategy.
2. Monoculture runs were performed for each of the five above variants (ϕ RL = 0.5, 0.6, 0.7, 0.8, and 0.9) to identify the most productive strategy in monoculture.Each run simulated the dynamics of a single variant for 500 years.
3. Invasion runs were performed for six pairwise combinations of four variants (ϕ RL = 0.6, 0.7, 0.9, and 1.0; see Table 3 for details of the combinations at the two [CO 2 ] levels) to confirm the identity of the most competitive strategy identified in the polyculture runs.Each invasion run included two different variants: a "resident" variant and an "invader" variant.We first ran the model with only the resident present for 400 years, which was long enough for it to come close to an equilibrium state.At the beginning of year 401, we converted 5 % of the population in each resident cohort into a new invader cohort by changing ϕ RL .We then ran the model for a further 240 years to get the DBH growth rates of invaders.To determine whether a ϕ RL = X was an evolutionarily stable strategy (ESS, the strategy when in monoculture that cannot be invaded) we examined runs in which the resident had ϕ RL = X and the invader had ϕ RL = X ± δ.
We also verified that the ESSs at the two CO2 concentrations are convergence stable (Geritz et al., 1998) by examining runs in which the resident had ϕ RL = X ± δ (with δ = 0.1 or 0.2) and the invader had ϕ RL = X.

GPP, NPP, tree growth rates, and abundances
Below, we focus on annual to successional timescales because diurnal and seasonal patterns are caused by the structure of the biophysical parameterizations in LM3-PPA (Appendix A and B), which are identical to those in LM3, have been under development for more than a decade, and are reviewed elsewhere (Shevliakova et al., 2009;Milly et al.,  S2. The model-simulated annual GPP and NPP for the Willow Creek AmeriFlux site are close to estimates from eddycovariance and biometric data collected at the same site (Fig. 2a; Desai et al., 2005;Curtis et al., 2002).NPP in the model was 48 % of GPP at the approximate steady state.The slight decline of GPP after forest closure was caused by self-thinning (Fig. S3a).Model predictions in Fig. 2 are taken from the monoculture sugar maple runs, but the threespecies runs predicted very similar values after the first 20 years (Fig. S4).
The allocation of NPP to leaves, fine roots, and woody biomass predicted by LM3-PPA is roughly similar to the measurements in Curtis et al. (2002), with the allocation to wood being too high and the allocation to leaves and roots too low (Fig. 2b).We did not tune the model to better predict the allocation data at Willow Creek, in part because the difference between the model and data could be caused by the fact that we simulated only 1 or 3 of the ∼ 10 species at Willow Creek.Because the allocation scheme assumes that NSC is allocated preferentially to the leaf and fine-root targets, interannual variation of sapwood and seed production is greater than that of leaves and fine roots (Fig. 2b).
DBH growth rates in the canopy layer are much higher than in the understory (Fig. 2c) because of shading (Fig. S2a).The predicted DBH growth rates of upper-canopy trees agree well with those derived from FIA data (Zhang et al., 2014) for all three species (Fig. 3).Predicted understory growth rates for sugar maple also agree well with estimates from FIA data, but predicted understory growth rates for red maple and trembling aspen are lower than estimates from FIA data (Fig. 3).
With initial population densities taken from earlysuccessional FIA plots (Table C4), the LM3-PPA model correctly predicts the subsequent successional turnover of trem- bling aspen, red maple, and sugar maple (compare Fig. 4a  and b).The transition from trembling aspen to sugar maple dominance is caused primarily by low survivorship of aspen in the understory, which was due to a combination of growth suppression from shading (which keeps cohorts in small size classes, where understory mortality rates are highest; Fig. 1a) and aspen's relatively high background rate of understory mortality (Table 1 and Eq. 8).Mortality due to carbon starvation rarely occurred in our simulations, although this may simply reflect our parameterization of mortality, which attributes high rates of mortality in small size classes to "back- Closed and open symbols are for uppercanopy ("Top") and understory ("Under") trees, respectively.The FIA data used to estimate observed growth rates are from the northern lake states (Michigan, Wisconsin, and Minnesota), USA.Canopy growth rates were estimated by combining trees with a reported crown class of "dominant" or "co-dominant", and understory growth rates were estimated from trees with a crown class of "overtopped" (Zhang et al., 2014).ground mortality" (Fig. 1a), with "starvation mortality" occurring in our model only if NSCs drop to zero.The timing of the transition from aspen to sugar maple is set primarily by the longevity of aspen canopy trees.
The model-predicted size distributions of both numbers and biomass for stands at 40-60 and 80-100 years are also qualitatively similar to FIA data (Fig. 5), despite significant quantitative differences in tree numbers.These differences are likely to be caused primarily by a combination of model error, the fact that our simulations included only a subset of species in the FIA plots, and differences between the initial conditions of early-successional plots today (which were used to initialize the simulations) and those 40-100 years ago (when succession began in the 40-100-year-old FIA plots).
The number of small trees in the baseline LM3-PPA model (H0; see Fig. 1 and Table 2) is significantly reduced near the late-successional equilibrium (Fig. 6a; mean model state from 600 to 1000 years).Moreover, the size distribution predicted for these old-growth forests has considerable biomass in trees larger than 1.2 m in diameter, which is unrealistic for these species (Fig. 6c).The alternative model H1 (high mortality rate for large trees) removes the unrealistically large trees.Like H1, cessation of crown area expansion at high DBH (H2) reduces the predicted number of very large trees.H2 also predicts a decline in DBH growth rate as trees become very large (Fig. S6), which is consistent with observations (Sillett et al., 2010;Lorimer et al., 1999).The random selection of understory trees to fill canopy layer gaps (H3) has little impact on size and biomass distributions (Fig. 6).GPP and NPP (Fig. S7a) and allocation of NPP to leaves, fine roots, and sapwood (Fig. S7b) simulated with the three alternative assumptions were close to those simulated by the default model (H0).The assumption of high mortality rates of very large trees (H1) led to reduced woody biomass since this assumption increased the mean turnover rate of wood, but it did not significantly affect equilibrium soil C. Assumptions H2 and H3 had little impact on C storage in wood or in the soil (Fig. S7c).

Effects of vegetation dynamics on vegetation and soil C storage
Comparisons of the predictions of LM3-PPA to those of the standard BGC model (Fig. S1b), forced with the same GPP and NPP from LM3-PPA, highlight the effects of successional diversity on carbon storage.The single-species runs of LM3-PPA include a dominant species for the region (sugar maple), which is dominant precisely because it is a long-lived late-successional species (Burns and Barbara, 1990).Parameters for the standard BGC model were chosen to be consistent with the one-species LM3-PPA model, and so, as expected, the BGC model and the single-species runs of LM3- PPA predict similar patterns of biomass and soil carbon storage (Fig. 7a and b).
In contrast, the three-species runs of LM3-PPA are dominated early in succession by a pioneer species (trembling aspen), which is short-lived, perhaps because its low wood density trades resistance to disease and windthrow for rapid height growth (Burns and Honkala, 1990).As a result, threespecies runs of LM3-PPA predict lower carbon storage in the woody biomass C pool (Fig. 7a) and higher soil carbon (Fig. 7b) early in succession than the standard BGC model or the single-species runs of LM3-PPA.The woody biomass C pool with one species needs ∼ 300 years to reach equilibrium, whereas the three-species runs need more than 500 years (Fig. 7a).
In the standard BGC model, the turnover rate of the woody biomass carbon pool was set as the mean mortality rate of sugar maple trees in the canopy layer (0.012 yr −1 ).In contrast, in the LM3-PPA simulation with one species, there was a peak in the biomass turnover rate because of the selfthinning of trees that had been pushed into the understory after canopy closure (red dashed line in Fig. 7c).In the LM3-PPA simulation with three species, the biomass turnover rates were much higher early in succession than in the singlespecies run because the mortality rates of aspen, and to a lesser extent red maple, are higher than that of sugar maple (green dashed line in Fig. 7c).The peak in the biomass turnover rate in the three-species run early in succession is caused by self-thinning following canopy closure, which occurs at a younger stand age than in the single-species run.As the models approached their equilibrium states, the carbon in biomass and soil pools converged because the inputs (NPP) and the residence times in biomass and soil C pools converged (Fig. 7).

Competitively optimal allocation strategy at different atmospheric CO 2 levels
After 500 years of competition among five allocation strategies of red maple (with the ratio of crown LAI to fineroot area, ϕ RL , ranging from 0.5 to 0.9 for upper-canopy trees) in the "polyculture runs", the variant with ϕ RL = 0.7 had the highest basal area at preindustrial [CO 2 ] (280 ppm), whereas ϕ RL = 0.9 had the highest basal area at doubled [CO 2 ] (560 ppm; Fig. 8).These results suggest that ϕ RL = 0.7 and ϕ RL = 0.9 are approximate competitive optima at 280 and 560 ppm, respectively.The precision of the approximations is limited by the resolution of the experiments (five discrete values of ϕ RL ).These approximate competitive optima were confirmed to be approximate ESSs by two-species "invasion runs" in  1 and Fig. 4a.The standard BGC model is summarized in Fig. S1b.which an equilibrium monoculture of one variant (a species with a given value of ϕ RL ) competed against an invading alternative variant (a species with a different value of ϕ RL ) that was initially rare.At [CO 2 ] = 280 ppm, ϕ RL = 0.7 was the competitively optimal strategy since it could not be invaded by any other variant and could invade all other variants (i.e., the convergence-stable ESS; Geritz et al., 1998), and at [CO 2 ] = 560 ppm ϕ RL = 0.9 was the competitively optimal strategy (Fig. 9).
Using the results in Farrior et al. (2013), it is possible to show mathematically that -for the case considered here, where understory traits are constant across species/PFTsthe competitive optimum (ESS) reduces to the strategy with the highest woody NPP when in the canopy and when in competition with the other strategies.Note also that species rankings of lifetime reproductive success, woody NPP, and DBH growth rate are equivalent here because all variants share the same other vital rates, wood density, and stem allometry.In the polyculture simulations, the strategy with the highest woody NPP or DBH growth rate in the canopy (over the last 60 simulation years) was ϕ RL = 0.7 at preindustrial [CO 2 ], and ϕ RL = 0.9 at doubled [CO 2 ] (Fig. 10), which further confirms the CO 2 -induced allocational shift implied by the results described above.The mechanisms causing this al- locational shift under elevated [CO 2 ] are explored in detail in the Discussion.Here, we simply note that these results imply that woody carbon sinks caused by elevated [CO 2 ] will be reduced by competitively optimal shifts in allocation away from long-lived woody tissues and toward short-lived fine roots, either because of an evolved plastic response or because a species or genotype with a larger ϕ RL will become competitively dominant under elevated [CO 2 ] (Farrior et al., 2013).
In contrast, among the "monoculture runs", the strategies with the highest canopy woody NPP and DBH growth rates were ϕ RL = 0.6 and ϕ RL = 0.7 for preindustrial and doubled [CO 2 ], respectively (Fig. 10).Both of these monoculture optima have higher allocation to wood and less allocation to fine roots than monocultures of the corresponding competitive optima (ϕ RL = 0.7 and ϕ RL = 0.9 at preindustrial and doubled [CO 2 ], respectively).Note that, in Fig. 10, competitively optimal growth rates are sometimes higher than those for the monoculture optima.This is because the competitively optimal growth rates in Fig. 10 are from polyculture runs, where individuals of the most competitive strategy have access to more water than in a monoculture of their own strategy; that is, in polyculture, individuals of the most com- In each experiment, the resident type was simulated for 400 years in monoculture, and then a small fraction of its density was converted to the invading type.The competitive optimum (ϕ RL = 0.7 and ϕ RL = 0.9 at 280 and 560 ppm, respectively) is the type (ϕ RL ) that cannot be invaded and can invade all other types (i.e., the convergence-stable evolutionarily stable strategy, ESS).
petitive strategy compete against individuals whose fine-root density is lower than that of the most competitive strategy.
To understand how differences between the monoculture and competitive optima arise, consider the following example.Under preindustrial [CO 2 ], ϕ RL = 0.7 had higher DBH growth rate than ϕ RL = 0.6 when invading a monoculture in which light and water availabilities were determined primarily by ϕ RL = 0.6.For this reason, the model predicts that ϕ RL = 0.7 will competitively exclude ϕ RL = 0.6, even though it will have a lower equilibrium growth rate once it has taken over the stand (because ϕ RL = 0.7 has a lower growth rate in conditions created by ϕ RL = 0.7 than ϕ RL = 0.6 has in conditions created by ϕ RL = 0.6).These differences between the competitive (polyculture) and non-competitive (monoculture) optima illustrate that plant strategies predicted by naïve (e.g., productivity-maximizing) optimization algorithms are often at odds with predictions from game-theoretic (ESS) competitive optimization (McNickle and Dybzinski, 2013;Farrior, 2014).
Figure 11 contains additional results that will be used in the Discussion to explain the predicted allocational shift caused by elevated [CO 2 ].It reports the percentage differ- ence between two runs of a monoculture of ϕ RL = 0.7 at [CO 2 ]= 560 ppm and at preindustrial [CO 2 ] for each of five quantities.A doubling of [CO 2 ] increased the fraction of each growing season in which canopy trees were watersaturated (defined as the fraction of days during the growing season in which water supply was greater than or equal to demand at 14:00 LT over the final 60 years of a 500year run) by 21 %.The water use efficiency (WUE; GPP per unit transpiration) of canopy trees during the water-limited period (days in which water supply was less than demand at 14:00 LT) increased by 79 %.The change in the length of the water-saturated period is relatively small (21 % increase, compared to a 79 % increase in WUE during the water-limited period) because of biophysical feedbacks in the model.Specifically, although a doubling of [CO 2 ] decreased transpiration by 4.55 % for the whole tile, this change was offset by a 1.78 % increase in the sum of evaporation and runoff.In absolute terms, the decrease in transpiration was 10.1 mm yr −1 , while the increase in evaporation plus runoff was 10.2 mm yr −1 , which canceled out the effect of increased [CO 2 ] on mean growing-season soil moisture (152.49mm at preindustrial [CO 2 ] and 152.91 mm at doubled [CO 2 ]).

Overview
In this paper, we describe the biophysical coupling between the height-structured PPA forest dynamics model and the GFDL LM3.The new model, LM3-PPA, was developed for future Earth system model (ESM) simulations in which vegetation dynamics are based on individual-level resource competition among size-structured cohorts of plants belonging to multiple species or PFTs.Our paper describes (1) the details of the biophysical coupling between LM3 and PPA, (2) preliminary model evaluation for a single site in the northeastern USA, (3) simulation experiments involving multiple allocational types at different atmospheric CO 2 concentrations, and (4) an interpretation of these competition experiments based on a mathematically tractable version of the PPA model.LM3-PPA is among the first land models to represent individual-level resource competition -including height-structured competition for light -and is the only land model to date that is closely tied to a mathematically tractable forest dynamics model, which affords a greater level of understanding of land model behavior than would be possible otherwise.Our paper is novel because we present novel land model predictions of how resource competition affects allocation to wood (a long-lived C pool) vs. fine roots (a shortlived C pool) at different CO 2 levels, and because we show how these land model predictions can be understood in the context of analytical predictions derived from a mathematically tractable version of the PPA model, as explained in Sect.4.5 below.

Model evaluation
The comparisons between the model's predictions and data at various scales (Figs.2-5, and S5) are intended as an initial evaluation and validation of LM3-PPA.The comparisons show that the model produces reasonable fast-timescale carbon and water dynamics (Supplement) as well as reasonable annual values for GPP and NPP (Fig. 2).The model also makes realistic predictions of individual growth rates, population structure (Fig. 5), and forest succession (Fig. 4).These comparisons must be evaluated in light of the tuning of the physiological model to produce observed NPP, the tuning of a single parameter affecting diameter growth, and the tuning of the elevated mortality of seedlings and small saplings.
The model formulation predicts tree-and ecosystem-level allocation patterns that are supported by a number of empirical studies.In LM3-PPA, the ratio of NPP to GPP and the fraction of NPP allocated to the three main plant structural C pools (foliate, fine roots, and wood) are not assumed to directly depend on tree size and stand age.Nonetheless, foliage and fine-root biomasses equilibrate in the model more than an order of magnitude more quickly than woody biomass.Experimental studies have indeed found that leaves and fine roots reach equilibrium quickly, long before total biomass reaches equilibrium (Goulden et al., 2011).Studies have also found that the ratio of autotrophic respiration to GPP is independent of age (Ryan et al., 2004), which is consistent with our model.Note that this is contrary to the expectation that maintenance respiration of stems should increase with tree size if it is proportional to sapwood biomass.Instead, LM3-PPA assumes that stem maintenance respiration is proportional to crown area, which -like fine-root surface area -is assumed to be proportional to DBH 1.5 (see Dybzinski et al., 2011;Farrior et al., 2013).This is consistent with the finding that bole respiration per unit of biomass decreases with age (Ryan et al., 2004).Also, it is possible to show that, if NPP and crown area are proportional to DBH 1.5 , and both DBH growth rate and fractional allocation of NPP to wood are size-independent, then wood biomass should be proportional to DBH 2.5 , as it is in the model and in empirical reports (e.g., Jenkins et al., 2003;Wang, 2006).
Because LM3-PPA is based on macroscopic equations from gap simulators (Strigul et al., 2008), forest inventory data can also be used to evaluate the model.LM3-PPA was tuned to reproduce canopy tree growth rates for three tree species near Willow Creek, but it was not tuned to fit un-derstory growth rates, which therefore provide useful tests of model performance.Observed understory growth rates for the two least shade-tolerant species were underpredicted (Fig. 3; note that uncertainties in mean growth rates are much smaller than the variances in the growth observations shown by the error bars in Fig. 3).One likely reason for this model-data discrepancy is that shade-intolerant species such as trembling aspen tend to experience darker understory conditions in our simulations (which assume homogeneous light conditions within each understory layer) than in real forests, where saplings of shade-intolerant species tend to occur in unusually bright understory locations (Clark and Clark, 1992;Davies, 2001;Poorter and Arets, 2003;Lichstein et al., 2010).
LM3-PPA also predicts the observed successional turnover of trembling aspen, red maple and sugar maple and size structure in the forests of the northern lake states, USA (Figs. 4 and 5; see also Woods, 2000;Purves et al., 2008).The model's ability to make detailed 100-year predictions that are consistent with data from successional chronosequences is not surprising, because forest simulators have been succeeding in this type of prediction for decades.However, it does reaffirm the value of constructing a DGVM from the scaling algorithms in forest gap simulators.
Although LM3-PPA successfully captures the main features of secondary forest succession in the northern lake states, USA (as does the PPA model; Purves et al., 2008), we would not expect LM3-PPA to maintain successional diversity indefinitely in old-growth forests.This is because LM3-PPA (like the PPA model) does not represent the gap-scale disturbances that shade-intolerant species require for persistence in old growth.Future implementations of LM3-PPA may include the gap-dynamics approximation from the ED model (Moorcroft et al., 2001), which should allow successional diversity to be maintained in old growth, and which may also capture other forms of spatial heterogeneity (e.g., the presence of emergent trees in some tropical forests).As explained in Sect. 5 of Appendix A, the ED gap-age approximation is already built into the LM3-PPA model code (but was not used in the simulations presented here).

Alternative assumptions about effects of size and age on growth and mortality
In the baseline LM3-PPA model (H0 in Table 2), canopy tree mortality rates are constant and independent of tree size and age, and canopy tree diameter growth rates remain roughly constant after approaching an asymptote when trees are still small (see text below Eq. 17).As a result, the model predicts unrealistically large trees in old forests (Fig. 6).Although this is a common problem of forest simulators, it is often not very important in regions of the world where little old growth remains (e.g., the temperate zone) or where stand-replacing natural disturbances are relatively common (e.g., fire-prone boreal forests).We explored alternative assumptions about growth and death rates of very large trees in this paper, primarily because LM3-PPA will ultimately need to perform in regions, such as the wet tropics, where old-growth forests are more common.Of the hypotheses examined (H0-H3), size-dependent decrease in the exponent relating crown area and diameter (H2) provides the best mix of empirical support and ability to produce realistic size distributions.Note, however, that none of the alternative assumptions about large trees has a large effect on predicted ecosystem-level carbon fluxes or storage in 600-1000-year-old forests that are at quasi-equilibrium (Fig. S7).

Effects of vegetation structure and successional diversity on C dynamics
For the tests that we have applied to date, the extra structure and diversity in LM3-PPA has relatively little effect on diurnal patterns of fluxes or annual NPP and GPP but does affect long-term carbon accumulation.The successional effects of size structure are best seen in the three-species run in Fig. 7c (green dashed line), where the biomass turnover rate first climbs by ∼ 30 % and then falls by more than a factor of 3 over the first 200 years of succession because of the successional transition from aspen, which has a high mortality rate, to sugar maple, which has a low mortality rate.As a result, carbon accumulation in the three-species run of LM3-PPA is significantly lower than that in the single-species run for more than 200 years (Fig. 7a).
The woody carbon accumulation rate after t years of succession in a simple biogeochemical box model is approximately α w NPP × e −µt (where α w is the fraction of NPP allocated to wood and µ is the annual tree mortality rate; Weng et al., 2012).Thus, the biomass growth rate in the standard BGC model exponentially decays over time to yield the asymptotic biomass accumulation curve in Fig. 7a (solid line).In contrast, in the PPA, an even-aged cohort of shadeintolerant saplings will self-thin so that the sum of their crown areas equals the area of the disturbance they are competing to fill.That is, the number of individuals in the cohort, n(t), tends to be proportional to the reciprocal of an individual's crown area, A CR (D(t)).Since total biomass is simply individual biomass, b(D(t)), multiplied by n(t), total stem biomass tends to be proportional to b(D(t))/A CR (D(t)), which -given the allometric constants for wood biomass, S(D(t)), and A CR (D(t)) -is simply proportional to diameter, D(t) (see Eq. 4).Finally, because diameter grows at an approximately constant rate after saplings reach ∼ 10 cm in diameter (around year 30 in Fig. 2c), LM3-PPA predicts linear biomass growth for an extended period when shadeintolerant species are present, like the green dashed line in Fig. 7a, and as observed in real chronosequences (Yang et al., 2011).

Competitive optimization and ecosystem C storage
When [CO 2 ] doubles from 280 to 560 ppm, the most competitive strategy in LM3-PPA shifts toward trees with greater allocation to fine roots and less allocation to wood .This is important because it would reduce the carbon sink caused by CO 2 fertilization.Thus, competitive optimization provides a way to discover carbon cycle feedbacks that involve changes in ecosystem-level allocation.
Elevated [CO 2 ] leads to greater leaf-level or intrinsic water use efficiency (WUE; carbon fixation per unit transpiration) in LM3-PPA, as observed in CO 2 enrichment experiments (Norby and Zak, 2011).Higher leaf-level WUE in LM3-PPA increases leaf productivity during the water-limited period of the growing season, while also decreasing the proportion of the growing season that plants spend in water limitation.These two responses to increased [CO 2 ] have opposing effects on the most competitive fine-root allocation strategy (i.e., the evolutionarily stable strategy, ESS; Farrior et al., 2013Farrior et al., , 2015)).ESS root allocation increases with increasing productivity (due to high water availability or high water use efficiency) during the water-limited period (up until the point where plants are water-saturated, and thus no longer waterlimited) for competitive reasons related to "the tragedy of the commons" for water use in plants (Gersani et al., 2001;Zea-Cabrera et al., 2006;Farrior et al., 2013).In contrast, ESS root allocation decreases as the length of the watersaturated period increases because roots represent a respiratory sink when plants are water-saturated.The net effect of an increase in [CO 2 ] on the ESS depends on the quantitative balance between these two opposing forces (Farrior et al., 2015), and thus depends on the full suite of biophysical feedbacks present in a model like LM3-PPA that must exchange matter, energy, and momentum with the atmosphere.In the case study presented here, increased evaporation and runoff largely compensate for reduced transpiration under elevated [CO 2 ], so that [CO 2 ] has little effect on mean soil moisture or the total number of hours each growing season during which plants are water-saturated (Fig. 11).In contrast, increased evaporation and runoff under elevated [CO 2 ] do not attenuate the expected increase in leaf productivity (due to increased WUE) during the period when water is limiting.The upshot, in our case study, is that of the two opposing forces on ESS fine-root allocation -(1) a decrease in root allocation due to an increased period of water saturation, vs. (2) an increase in root allocation due to increased leaf productivity during the water-limited period -the latter effect dominates, and the most competitive strategy shifts to one with greater allocation to fine roots .This result has now focused our attention on the strength of the biophysical feedbacks in LM3 and LM3-PPA, which might be too strong.The important point here is that we know what to focus on only because of the understanding afforded by the connection between LM3-PPA and the analytically tractable PPA model (Farrior et al., 2013(Farrior et al., , 2015)).We understand the pre-dicted feedback in LM3-PPA involving [CO 2 ], water, fineroot allocation, and carbon storage only because the model may be interrogated analytically.

Future challenges
In this paper, we do not provide parameter values needed to implement LM3-PPA at the global scale using PFTs or more flexible trait-based approaches (e.g., Scheiter et al., 2013;Wullschleger et al., 2014).The PPA has previously been applied to other temperate forest types that include conifers (e.g., Purves et al., 2008;Strigul et al., 2008), as well as tropical forests with more than two canopy layers (Bohlman and Pacala, 2012), and we are currently developing parameter values for non-tree vegetation types, such as shrubs and grasses (Weng et al., unpublished).The formalism we describe in this paper requires no structural changes to work in non-forested ecosystems, including those with open canopies or with no competition for light (i.e., because of severe water limitation).Furthermore, as explained in Appendix A, the current version of the LM3-PPA code can already accommodate land use change, secondary forest management, stand-replacing disturbance, and the ED approximation for canopy gap dynamics, which is required to maintain successional diversity in old-growth forests with low rates of standreplacing disturbance.In summary, LM3-PPA can, in principle, be extended to global-scale simulations in fully coupled ESM experiments with little modification to the processes already encoded in the model.
In addition to developing parameterizations for globalscale applications, another important area for future work is to better understand the transient dynamics of vegetation response to global change.Our results suggest potentially important effects of allocational shifts, driven by competition among plants for light and water under elevated CO 2 , on terrestrial carbon balance.However, our competition experiments were designed only to identify the eventual outcome of competition under a given set of conditions, and are therefore agnostic about the rate and pathway of response.In reality, allocational shifts could be potentially rapid (e.g., tracking environmental conditions on an annual timescale) if individual plasticity were sufficient (Franklin et al., 2012), would occur over intermediate timescales (e.g., decadal) if allocational shifts required shifts in relative abundances of species already present within a landscape, and would be even slower if allocational shifts required long-distance migration by dispersal-limited species (Lischke et al., 2006;Snell et al., 2014) and/or the evolution of novel types (Valladares et al., 2007).Empirical evidence suggests that intraspecific variation in allocation is often sufficient to accommodate the shift in competitively optimal allocation predicted by LM3-PPA under a doubling of atmospheric CO 2 (R. Dybzinski, unpublished analysis), and free-air CO 2 enrichment (FACE) experiments demonstrate considerable individual plasticity in allocation to leaves, wood, and fine roots (Jackson et al., 2009;McCarthy et al., 2010;Norby and Zak, 2011;Iversen et al., 2012).However, there are clearly limits to plasticity (Valladares et al., 2007), and it is unknown whether the plastic responses of individuals to environmental change (which evolved over the last ∼ 20 million years under relatively low atmospheric CO 2 concentrations; Zachos et al., 2001) would be the competitively optimal responses under future novel conditions.A key challenge, then, is to better understand the transient dynamics that ecosystems will undergo as they approach competitive equilibria from different initial conditions.

Conclusions
We present a model, LM3-PPA, which simulates vegetation dynamics and biogeochemical processes by explicitly scaling from individual plants to ecosystems using the perfect plasticity approximation (PPA).The model is formulated to be the land surface component of an Earth system model.It includes height-structured competition for light and root allocation-dependent competition for belowground resources (water in this study).The partitioning of space by plant crowns following the rules of the PPA to form canopy layers simplifies the simulation of light competition among trees and allows the LM3-PPA model to predict forest succession with an explicit description of the size distributions of individuals within each species or functional type, in addition to the predictions of carbon fluxes of an ecosystem (GPP, NPP, and R a ), the dynamics of soil organic matter and decomposition (heterotrophic respiration, R h ), evapotranspiration, and soil hydrology.Because of the tractability of the PPA, the coupled LM3-PPA model is computationally efficient (relative to existing alternatives to modeling height-structured, individual-level competition within ESMs) and retains close linkages to mathematically tractable special cases (e.g., constant climate).
Comparisons of model simulations with data show that the model makes reasonable predictions for diurnal and annual carbon and water fluxes, growth rates of individual trees, and population sizes and species turnover during succession.The model marginally underpredicts the growth rates of shadeintolerant species in the understory and seriously overpredicts of abundances of very large trees in old growth.The overestimate of large trees can be corrected by adding either size-specific mortality or size-specific crown area allometry, both of which are supported by some studies.The model also shows that within-functional-type successional diversity has significant ecosystem-level effects at timescales up to a century or more.Finally, simulation experiments show that the dominant competitor's root-leaf-stem allocation pattern shifts as a function of the atmospheric CO 2 concentration and predict that carbon sinks caused by CO 2 fertilization in forests limited by light and water will be down-regulated if allocation tracks changes in the competitive optimum.These results indicate that the ecological strategies functioning at the scales of individuals and communities, which are usually missing in ESMs, have strong impacts on biogeochemical processes and their responses to climate changes.
The implementation of the model in this paper is for temperate broadleaved forest trees, but the formulation of the model is general and can be expanded to include other growth forms and physiologies.The model can accommodate an arbitrary number of functional types, species, and/or genotypes in competition with one another across the terrestrial regions of the globe.

Appendix A: Vegetation dynamics and subgrid-scale heterogeneity A1 Vegetation structure: cohorts and canopy layers
In the coupled LM3-PPA model, the vegetation is represented as a set of cohorts arranged in different vertical canopy layers according to the perfect plasticity approximation (PPA) model.Each cohort is a collection of identical individual trees of the same size and of the same PFT or species.In describing the model, we sometimes refer to "individuals" or "trees" to provide biological context and most of the equations below are at the individual level, but our implementation of the PPA model is in reality at the level of cohorts, each of which is defined by its PFT or species identity, the size of the identical individuals in the cohort, and the spatial density of these individuals (number per unit ground area).Any reference to "individuals" or "trees" should be understood to refer to the properties of the identical individuals in a given cohort.The PPA model allows for flexibility in the shapes of individual tree crowns (Strigul et al., 2008;Purves et al., 2007), but, for simplicity, here we assume that trees have flat-topped crowns, which allows for accurate predictions of observed succession and canopy structure in broadleaved temperate forests (Purves et al., 2008;Zhang et al., 2014) and canopy structure in a neotropical forest (Bohlman and Pacala, 2012).Individual tree height is defined as the height at the top of the crown, and all foliage of a given cohort is assumed to belong to a single canopy layer, which simplifies the energy balance equations for multi-layered canopies (Appendix B).The height of canopy closure for layer k (k = 1 is the top layer, k = 2 is the second layer, etc.) is referred to as Z * k , the height of the shortest tree in the layer, and is defined implicitly by the following equation: where N i (Z, t) is the density of PFT i trees of height Z per unit ground area, A CR,i (Z * k , Z) is the crown area of an individual PFT i tree of height Z, and η is the proportion of each canopy layer that remains open on average due to spacing between individual tree crowns.The top layer includes the tallest cohorts of trees whose collective crown area sums to 1 − η times the ground area, and lower layers are similarly defined.Trees within the same layer do not shade each other, but there is self-shading among the leaves within individual crowns.Cohorts in a sub-canopy layer are shaded by the leaves of all taller canopy layers using a mean field approximation; that is, in a given canopy layer, all cohorts are assumed to have the same incident radiation on the top of their crowns (Fig. A1; see also Appendix B for details on radiation transfer within and between canopy layers).The gap fraction η increases light penetration through each canopy layer and allows for the persistence of understory trees in monoculture forests in which the upper canopy builds a physiologically optimal number of leaf layers, i.e., one in which its lowest leaves are at zero carbon balance (Dybzinski et al., 2011;Farrior et al., 2013).

A2 Allometry, allocation, growth, and respiration
In this section, we describe the detail of the growths of leaves, stems, and roots of the trees in the model.Each individual is composed of five tissues: leaf, fine root, sapwood, heartwood, and labile carbon stores (nonstructural carbohydrates, NSC).Empirical allometric equations relate woody biomass (including coarse roots, bole, and branches), crown area, and stem diameter.Photosynthate and retranslocated carbon enter the NSC pool, and carbon for respiration and growth are removed from it.The carbon allocation rules track PFT-specific targets for leaf area per unit crown area (l * ), fine-root area per unit leaf area, and the NSC pool size, where PFTs may be defined at any level (e.g., ecotype, species, conifer/angiosperm, N fixer/non-fixer).The diameter growth rates and sizes of individual trees are calculated from individual-level allometry and allocation of assimilated carbon.To scale up to a cohort from individual trees, individual pools and fluxes are multiplied by the spatial density of individuals in a cohort.
As described in the main text, the individual-level dimensions of a tree, i.e., height (Z), biomass (S), and crown area (A CR ), are given by empirical allometries (Dybzinski et al., 2011;Farrior et al., 2013): where Z is tree height, S is total woody biomass carbon (including bole, coarse roots, and branches) of a tree, α c and α Z are PFT-specific constants, θ c = 1.5 and θ Z = 0.5 (Farrior et al., 2013; although they could be made PFT-specific if necessary), π is the circular constant, is a PFT-specific constant, and ρ W is PFT-specific wood density (kg C m −3 ).

Allocation of assimilated carbon (growth functions)
The carbon fluxes from the fast-timescale equations (hourly or half-hourly) are summed over the diurnal cycle to provide daily total carbon gain from photosynthesis (P s (t)) and loss from respiration (R a (t)) for each cohort.This carbon is added to or taken from the NSC pool once a day: where G L+FR (t) is the amount of carbon allocated to produce new leaves and fine roots minus the carbon retranslocated from senescing leaves and fine roots, and G W+F (t) is the carbon allocated to stem and seed production.The computations of G L+FR (t) and G W+F (t) are based on the "target" amount of leaf mass, L * (D,p); fine-root mass, FR * (D); and nonstructural carbohydrates, NSC * (D, p) for each cohort.These quantities change with the trunk diameter (D) and its phenological state (p).Following the pipe model (Shinozaki et al., 1964), the target leaf and fine-root biomass and sapwood cross-sectional area are related by the following equations, which are identical in form to the equations that relate the actual values of these variables (main text Eq. 5) with the inclusion of the time-dependent phenological control on target leaf biomass, p(t): where L * k (D, p) is the target leaf mass of canopy level k at given stem diameter (D), l * k is the target leaf area per unit crown area of a given PFT at canopy level k, A CR (D) is the crown area of a tree with diameter D, LMA is PFTspecific leaf mass per unit area, and p(t) is a PFT-specific function ranging from zero to one that governs leaf phenology.The phenology function p(t) is unchanged from LM3, except that it is updated daily rather than once per month as in LM3 and LM3V (Shevliakova et al., 2009;Milly et al., 2014).The onset of a growing season is controlled by two variables -growing degree days (GDD) and a weighted mean daily temperature (T pheno ) -while the end of a growing season is controlled by T pheno (see Sect. 3,below).FR * k (D) is the target fine-root biomass at diameter D and canopy level k, ϕ RL is the ratio of total root surface area to the total leaf area, SRA is specific root area, A * SW,k (D) is the target crosssectional area of sapwood at canopy level k, and α CSA is an empirical constant (the ratio of sapwood cross-sectional area to target leaf area).All plant tissues are assumed to be 50 % carbon by mass.Unless otherwise stated, units are mass in kilograms carbons, area in square meters, height in meters, and diameter in meters.
The target nonstructural carbohydrate pool, NSC*(D, p), is a multiple (q) of leaf target during the growing season plus the retranslocated carbon from leaves when a growing season ends: where l * is the target crown LAI in a given canopy layer (see Eq. A4) and q is a species-or PFT-specific parameter, with q chosen so that NSC * is several times the size of the combined target leaf and fine-root masses to reflect the reality that a healthy plant often stores enough NSC to refoliate after defoliation (Hoch et al., 2003;Richardson et al., 2013).The term 0.25(1 − p(t)) represents the retranslocated carbon stored by a tree in the non-growing season, when p(t) < 1, assuming that a quarter of leaf carbon is retranslocated after senescence (Vergutz et al., 2012).
The leaf mass target for an upper-canopy cohort depends on how recently it has moved from the understory to the canopy.After it has been in the canopy continuously for at least τ C years, the species-or PFT-specific target for a cohort's crown LAI during the growing season becomes l * C , which is greater than l * U .To avoid an unrealistically rapid depletion of the NSC pool when a tree moves from the understory to the canopy and its crown LAI target (l * ) increases, we define the growing-season l * for a canopy cohort that was last in the understory t * years ago, where t * < τ C , as LM3-PPA assumes that plants keep their leaves and fine roots tracking their targets if they have enough carbon in NSC: where the subscript k denotes the canopy layer (below, we use the values k = C and k = U for canopy and understory layers, respectively, but this can be generalized to an arbitrary number of layers, k).The first term in the Min function causes leaf mass (L) and fine-root mass (FR) to track their targets, L * k (D, p) and F R * k (D), respectively.Its Max function allocates from NSC to new leaves and fine roots only if the targets exceed the current masses.The bottom term in the Max function (0.25 ) is the carbon retranslocated from leaves to NSC when targets are smaller than the current masses at the end of a growing season.The second term in the Min function (f NSC NSC, with f NSC = 0.2 day −1 ) similarly reflects the maximum rate at which plants can convert NSC to structural tissue and prevents NSC from being suddenly depleted if an individual moves from the understory to the upper canopy, in which case leaf and fine-root targets increase suddenly from L * U (D, p) and FR * U (D) to L * C (D, p) and FR * C (D).It also define the rate of leaf flush at the beginning of a growing season (Wesolowski and Rowinski, 2006;Polgar and Primack, 2011).
The dynamics of leaf biomass of an individual is where γ L is the PFT-specific rate of leaf senescence triggered by the ending of a growing season.The new leaf biomass is converted into the change in leaf area by dividing by LMA.The total leaf area of a tree is converted to the tree's crown LAI, l k , by dividing by crown area.
Similarly, fine-root biomass is governed by The final term in Eq. (A3) (G W+F ) gives withdrawals from NSC to grow new wood and to produce seeds: The Max function causes stem and seed production to cease if NSC falls beneath its target.This typically happens in the model only when trees have negative carbon balance, in which case they stop growing in size and devote all labile carbon withdrawals to meet respiratory demand and replace leaves and fine roots that senesce.When the carbon balance is positive, plants spend their surplus on stems and seeds.
The parameter f SF is set at a value large enough (0.2 day −1 ) to keep NSC close to its target as the target increases because of plant growth.Seed production rate is assumed to be zero for understory trees and for canopy trees to equal: where v is the fraction of wood plus seed production that is devoted to seed (v = 0.1 for individuals in the canopy and v = 0.0 for individual in the understory).The cumulative biomass of seeds produced by a canopy cohort over a growing season of length T is converted to seedlings by dividing by the initial plant biomass (S 0 ) and multiplying by germination and establishment probabilities (p g and p e , respectively): where N (S 0 , t) is the spatial density of newly generated seedlings, and N (τ ) is the spatial density of the parent cohort at time τ .Finally, biomass of new wood growth is By differentiating the stem biomass allometry in Eq. ( A3) with respect to time, using the fact that dS/dt equals new sapwood biomass, and rearranging, we have which is the differential equation for diameter growth of the individuals in a cohort.The RHS of Eq. ( A15) is approximately independent of D (Farrior et al., 2013) because the numerator and denominator on the RHS are both usually roughly proportional to D 1.5 .The numerator tends to be proportional to D 1.5 because carbon gain is proportional to crown area, NSC surpluses tend to be a fraction of carbon gain, and crown area is usually roughly proportional to D 1.5 (Zhang et al., 2014).The denominator tends to be proportional to D 1.5 because θ Z tends to be about 0.5 (Zhang et al., 2014).The approximate diameter independence of Eq. (A15) allows many aspects of the behavior of the LM3-PPA model to be understood by referring to analytically tractable versions of the PPA model where dS/dt is assumed to be independent of D (e.g., Strigul et al., 2008;Dybzinski et al., 2011;Farrior et al., 2013).
Equations for height and crown area growth rates are obtained by differentiating the allometries for Z and A CR in Eq. (A4) using the chain rule: Finally, after summing the leaf and fine-root allometries in main-text Eq. ( 5), we have The time derivative of the RHS of Eq. (A19) must equal withdrawals from the NSC pool (from the RHS of Eq. (A7) without retranslocation): This equation can be solved for the dynamics of crown LAI (dl/dt):

Autotrophic respiration
The total autotrophic respiration rate of an individual is the sum of maintenance respiration of living tissues and the growth respiration for building new tissues: where R L , R SW , and R FR are the maintenance respirations of leaves, sapwood, and fine roots, and r g is a growth respiration constant (r g = 0.33 g C g −1 C).Maintenance respiration terms are calculated as where γ Leaf is a respiration coefficient of leaves; ε is a factor converting the unit of carboxylation rate V cmax (mol m −2 s −1 ) to kg C m −2 yr −1 ; β SW and β FR are respiration coefficients of sapwood and fine roots, respectively (kg C m −2 yr −1 for sapwood and kg C kg −1 C yr −1 for fine roots); ε is a factor converting the unit of kg C m −2 yr −1 to kg C m −2 day −1 ); A CB is cambium surface area (m 2 ), which we assume scales with diameter with an exponent 1.5 (A CB ∝ D 1.5 ), consistent with the height allometry exponent θ Z = 0.5; and f T is a temperature-dependent function adapted from Collatz et al. (1991Collatz et al. ( , 1992) ) that scales respiration rate with temperature: where T is • C.

Conversion from sapwood to heartwood
As trees grow, sapwood (SW) is transformed to heartwood (HW).This unidirectional process does not affect the size of the woody biomass C pool.We assume that, if the actual sapwood cross-sectional area A SW is larger than its target value, A * SW (D) (Eq.A2), the excess portion of sapwood biomass is converted to heartwood.Thus, to determine the amount of sapwood converted to heartwood in a given time step (HW), we simply calculate the difference between SW and the target sapwood C (SW * ) needed to balance L * and FR * : dHW = max(0, SW − SW * ). (A25) Using the equation for total tree biomass (main text Eq. 4), the target biomass of sapwood is where D is the diameter of the trunk and D HW is the heartwood diameter, which is given by where A HW is the cross-sectional area of heartwood.Assuming A SW is at its target value, The cross-sectional area of a trunk (A t ) is In addition, according to Eqs. ( A2) and (A4), the target crosssectional area of sapwood is

A3 Phenology
Here, the phenology for cold-deciduous plants used in the examples presented in this paper is described.The onset of a growing season is controlled by two variables, growing degree days (GDD), and a weighted mean daily temperature (T pheno ).The two variables are computed by the following equations: where t is the number of days from the end date of the last growing season and T d (t) is the daily mean temperature at day t.There are two thresholds for these two variables, GDD crit (320 day • C) and T crit (10 • C), respectively.When the criteria GDD(t) > GDD crit and T pheno (t) > T crit are met, a growing season is initiated by setting p(t) = 1.The ending of a growing season is controlled by T pheno .When T pheno (t) falls below T crit , the growing season is turned off (p(t) = 0), and leaves begin to senesce at an assumed rate (γ L , Eq. A8).A fraction of carbon (0.25) of senesced leaves is retranslocated to the NSC pool.

A4 Decomposition of soil organic matter
The soil carbon model is the same as in LM3V (Shevliakova et al., 2009), which describes soil carbon dynamics with a simplified variant of the CENTURY model (Parton et al., 1987).Dead plant tissues, including senesced leaves and fine roots, woody necromass due to mortality, and failed seeds, enter into two soil carbon pools (a fast turnover pool and a slow turnover pool) by different fractions according to their chemical composition, and then decompose at the rates regulated by soil temperature and moisture (see Sect. 2.7 of Shevliakova et al., 2009, for details).

A5 Subgrid-scale heterogeneity, land use change, and gap dynamics
The LM3-PPA model is implemented on a flexible grid (gridcell size is flexible and can be specified independently of the atmospheric model's grid), and also includes the dynamic tiling scheme for land use and subgrid-scale heterogeneity from LM3 (Shevliakova et al., 2009).Each grid cell is divided into tiles that differ in their history of land use and disturbance.Tiles may be cropland, pasture, primary vegetation, or secondary vegetation.Secondary vegetation is agestructured: each secondary vegetation tile has a different time since it was last harvested timber or since agricultural abandonment.The areas of the tiles in each grid cell change dynamically through time to simulate land use change.Each tile has its own set of cohorts, canopy airspace, soil texture, soil moisture, and undecomposed organic matter.Matter and energy are exchanged among grid cells and tiles within a grid cell because of coupling with the atmosphere, and because of water runoff, which collects to feed a river network.Tiles within a grid cell may be coupled by land use change.For example, when part of a forested tile is cleared for cropland, its soil organic matter is removed and deposited in the grid cell's cropland tile and the area of cleared land is added to the area of the cropland tile.Similarly, harvesting a portion of a forested tile causes a new tile to be produced, with successional age zero, and the material remaining in the harvested area is transferred to the new tile.To speed computation, two forested tiles with different successional ages may be merged if they become sufficiently similar.When two tiles are merged, all cohorts are retained.Cohorts within a tile are also merged if they become sufficiently similar.Finally, although the case study described in this paper involves only a single tile in a single grid cell, LM3-PPA is designed to allow seed dispersal among tiles within a grid cell and among different grid cells.This would require a modification of maintext Eq. (1).All transfers, merges, and divisions involving tiles and/or grid cells conserve matter and energy.
The model of land use change is unchanged from LM3 to LM3-PPA (details in Shevliakova et al., 2009).The land area that is converted from one land use type to another and/or harvested is specified by land use transition probabilities for each year and grid cell, which are generated as part of an external land use scenario.Because LM3 already tracks the age structure of secondary vegetation tiles, it is possible to configure its land use scheme to add the gap-disturbance approximation in the ED model (although this was not implemented in the current paper).ED divides the land surface into tree-sized cells, and combines forested cells into cohorts that differ in "gap phase age" -the time since the mortality of a tree larger than a threshold size (Moorcroft et al., 2001).This simulates the creation of light gaps in a fully spatial forest model.LM3-PPA would include the ED approximation if one were to create a new "secondary" tile containing only sub-canopy cohorts each time canopy trees were killed in a model grid cell.LM3-PPA, like its parent model LM3 (Shevliakova et al., 2009;Milly et al., 2014), represents subgrid-scale heterogeneity of the land surface by splitting each land grid cell (of flexible size) into multiple tiles.Each tile has distinct physical and biological properties as well as its own exchanges of energy, water, and CO 2 with the atmosphere.For example, different tiles within the same grid cell may represent natural vegetation, cropland, and multiple tiles for secondary vegetation last disturbed at different times in the past.The energy and mass exchanges are calculated separately for each tile, and the fluxes are aggregated to atmospheric grid cells, which may have different spatial resolution than land grid cells (Shevliakova et al., 2009).In LM3-PPA, the vegetation is represented by a set of cohorts arranged in different canopy layers according to the perfect plasticity approximation (PPA) within each tile (Appendix A).In the present study, we consider vegetation dynamics in a single grid cell with just one tile (natural vegetation).However, the model description presented here also applies to the general case with multiple grid cells, each with multiple tiles.
We refer to the canopy of cohort i as "canopy i."Each canopy i has its own temperature T v as well as amounts of intercepted water w l and snow w s .All cohorts exchange water, energy, and carbon dioxide with the common canopy airspace of mass m c , temperature T c , and specific humidity q c .Each cohort i is composed of individuals with density n i (individuals per m 2 of tile).With these assumptions, the energy balance of canopy i can be expressed as where Ĉi is the total heat capacity of canopy i, R Sv,i and R Lv,i are the net shortwave and longwave radiative balances of canopy k, Ĥv,i is the total sensible heat balance of canopy k, Lv,i is the total latent heat loss by canopy i, L f is the latent heat of water fusion, M i is the rate of melt of intercepted snow (or freezing of intercepted water, if negative), and the product L f M i is the heat associated with the phase transitions of the intercepted water.All the terms of Eq. (B1) are calculated per unit canopy area.These units are convenient for the energy balance calculations, especially for radiative transport in the multi-cohort canopy.The total heat capacity of the canopy i ( Ĉi ) is the sum of heat capacities of leaves (C v,k ), intercepted water (w l,k ), and intercepted snow (w s,k ): where c l and c s are the specific heat capacities of water and ice, respectively.The heat capacity of dry leaves (C v,i ) is assumed to be zero in the simulations presented in this manuscript.
The sensible heat term Ĥv,i in Eq. ( 1) is where H v,i is the sensible heat flux from the canopy airspace to canopy i due to turbulent exchange; H pl,k and H ps,k are the fluxes of heat carried by liquid and solid precipitation onto layer k; γ l,i and γ s,i are the fractions of liquid and solid precipitation on canopy i that are intercepted, equal to 1 − exp(LAI i ); and D l,i and D s,i are the rates of water and snow drip from canopy i, calculated as D l,i = w l,i /τ l and D s,i = w s,i /τ s , with timescales τ l = 6 h and τ s = 24 h for liquid and snow, respectively.The heat fluxes H pl,k and H ps,k carried by liquid and solid precipitation onto layer k are calculated simply as the heat content of precipitating water or snow multiplied by the intensity of rain P l,k or snow P s,k on top of this canopy layer (see Eqs. B7 and B8 below).We assume that the precipitation rates P l,k , P s,k and associated heat fluxes H pl,k and H ps,k are the same for all cohorts in layer k.
The latent heat term in Eq. ( B1) is where E t,i is transpiration, E l,i is evaporation of liquid intercepted water, and E s,i is sublimation of intercepted snow.L e (T ) and L s (T ) are the temperature-dependent specific heats of evaporation and sublimation, respectively.T u,i is the temperature of the water that canopy i transpires.Since we do not explicitly consider the heat exchange between water and the environment that occurs while the water is being transported to the leaves, the energy conservation law dictates that this temperature must be equal to the average temperature of water the plant uptakes from the soil.Liquid and solid water balance, respectively, for canopy i are dw l,i dt where P l,k and P s,k are liquid and solid precipitation on top of canopy layer k, γ l,i and γ s,i are the fractions of liquid and solid precipitation on canopy i that are intercepted, and M i is the rate of snow melt in canopy i. P l,1 is the rainfall on the upper canopy layer (k = 1), and P s,1 is the snowfall rate.For layers k > 1, we can write where f i is the area fraction of layer k − 1 occupied by cohort i.We assume that the drips from the canopy are not in-tercepted by the layers below, and so contribute directly to the water and energy balance of the underlying ground surface.With these assumptions, the canopy air water (specific humidity) balance equation, for all layers combined, is where E g is water vapor flux from the ground surface, E a is the water vapor flux to the atmosphere, m c is the mass of canopy air, and q c is the specific humidity of canopy air.The total water vapor flux E v,i from canopy i to the canopy air is a sum of three components: The canopy air energy balance is where c p and c v are specific heats of dry air and water vapor, respectively, H v,i is the flux of sensible heat from canopy k to the canopy air, H g is the sensible heat flux from the ground surface to the canopy air, and H a is the sensible heat flux from canopy air to the atmosphere.The energy balance of the ground surface is where R Sg is the net shortwave radiation absorbed by the ground, R Lg is the net longwave radiation absorbed by the ground, H g is the sensible heat flux between the ground and canopy air, L g is the latent heat of vaporization at the ground temperature, E g is the water vapor flux from the ground, G is heat flux from the ground surface to the underlying layers, L f is the latent heat of fusion, and M g is the rate of surface snow melt (or surface water freezing, if negative).Note that the "ground" can refer to either the soil surface or the top surface of the snow layer covering the soil.

B2 Propagation of solar radiation within cohorts and across canopy layers
The propagation of shortwave (R Sv,i ) and longwave radiation (R Lv,i ) in vegetation layers is calculated using a twostream approximation (Meador and Weaver, 1980), with the assumption of spherical leaf angular distribution (Pinty et al., 2006; Fig. B1).
To calculate the shortwave radiative balance R Sv,i for each canopy i (i.e., the canopy of cohort i), the equations of radiation transfer inside a single cohort's canopy can be written as follows: where L is the canopy depth from the top of the cohort canopy expressed in terms of LAI (leaf area per crown area), with l = 0 at the top of canopy i and l = LAI i at the bottom of canopy i; I ↑ and I ↓ are the upward and downward fluxes of diffuse radiation at depth L; G(µ 0 ) is the Ross function (Ross, 1975) for the solar zenith angle µ 0 ; γ i are the two-stream approximation coefficients (Pinty et al., 2006); F is the flux of direct solar radiation; and ω l is the singlescattering albedo of leaves.Eq. (B13) are solved following Liou (2002) to obtain the vertical profile of I ↓ and I ↑ within canopy i (which is assumed to reside in a single canopy layer; see Appendix A1).Given the boundary conditions of canopy i, we obtain its integral radiative properties: α i is the reflectance for diffuse radiation, τ k is the transmittance for diffuse radiation, τ i is the transmittance for the downward direct (solar beam) radiation, and s ↑ i and s ↓ i are coefficients of direct solar beam scattering upward and downward, respectively (these scattering coefficients, s, should not be confused with fluxes, S, in Fig. B1).
With the radiative properties of each canopy i, the equations for the radiation propagation through canopy layer k consisting of one or more cohorts (each with its own canopy k) can now be written as where f i is the fraction of area in layer k occupied by canopy i.Note that the subscripts k and k−1 refer to downward and upward radiation from layer k, respectively (Fig. B1).The summation is over all cohorts that belong to layer k, and η k is the total fraction of gaps in layer k.Solving the system of Eqs.(B14) yields the upward radiation from the ground surface: where α g and α g are the ground reflectance for diffuse and direct light, and index k = N refers to the fluxes underneath the entire vegetation canopy.The shortwave radiative balance for each cohort is Equations (B13-B16) are solved for two spectral bands (visible and near infrared) separately, with respective leaf and ground radiative properties, and the total shortwave radiative balance is calculated as the sum of the two.The longwave radiative balance terms, R Lv,i , are calculated similarly as the shortwave, except there is no contribution of direct solar light, and each of the cohort canopies emits longwave radiation according to the Stefan-Boltzmann law.

B3 Photosynthesis and stomatal conductance
We first calculate the net assimilation rate and stomatal conductance of leaves, integrated through the leaf area within a cohort's canopy, in the absence of soil water limitation.These values of assimilation and stomatal conductance imply a certain water demand.We then calculate available water supply and reduce the demand-based assimilation and stomatal conductance accordingly if water supply is less than water demand.The water-demand-based photosynthesis and stomatal conductance equations for a well-watered plant are modified from Farquhar et al. (1980), Collatz et al. (1991Collatz et al. ( , 1992)), and Leuning et al. (1995).We present equations for both C 3 and C 4 plants, although only the former are included in the examples presented in this paper.Consistent with the energy balance equations (Sect.2), we assume that the entire canopy of a given cohort is isothermal with temperature T v , and the air in the intercellular spaces is water-saturated with specific humidity equal to saturated specific humidity q*(T v ).The link between stomatal conductance (g s , mol m −2 s −1 ), the rate of net photosynthesis (A n , mol CO 2 m −2 s −1 ), intercellular concentration of CO 2 (C i , mol CO 2 mol −1 air), and the difference in specific humidity between the intercellular spaces and the canopy air (q a , kg H 2 O kg −1 air) can be expressed as a simplification of Leuning's (1995) empirical relationship assuming negligible cuticular conductance: where m is the slope of the stomatal conductance relationship, d 0 is a reference value of canopy air water vapor deficit (kg H 2 O kg −1 air), and * (mol CO 2 mol −1 air) is the CO 2 compensation point: where Net photosynthesis A n can be expressed as a CO 2 diffusive flux between canopy air and the stomata (Leuning et al., 1995): where C a is the concentration of CO 2 in the canopy air, and the factor 1.6 is the ratio of diffusivities for water vapor and CO 2 .We assume that the diffusion of CO 2 is mostly limited by stomatal conductance and not by the leaf boundary layer conductance, which we ignore for simplicity, following the formulation of the ED model (Moorcroft et al., 2001;Medvigy et al., 2009).Combining Eqs.(B17 and B21) yields the intercellular concentration of CO 2 : Following the mechanistic photosynthesis model of Farquhar et al. (1980), with extensions introduced by Collatz et al. (1991Collatz et al. ( , 1992)), we can also express net photosynthesis (A n ) as the difference between gross photosynthesis and leaf respiration, and assume gross photosynthesis is the minimum of several physiological process rates: where f T (T v ) is a thermal inhibition factor (see below); J E , J C , and J j are light-limited, Rubisco (CO 2 ) -limited, and export-limited rates of carboxylation, respectively; V m (T v ) is the maximum carboxylation velocity (mol CO 2 m −2 s −1 ); and γ is a constant relating leaf respiration to V m .The thermal inhibition factor, assumed to affect carbon acquisition and respiration equally, is . (B23) The maximum carboxylation velocity, V m , depends on the temperature of the leaf: where V max (the reference value of V m , at 25 • C) is a speciesspecific constant, f A (T v ) is given by Eq. (B19), and E v is the activation energy (see Appendix C).
For C 3 plants, Collatz et al. (1991): where a is the leaf absorptance of photosynthetically active radiation (PAR), Q is incident PAR per unit leaf area (E m −2 s −1 ), α LUE is the intrinsic quantum efficiency of photosynthesis (mol CO 2 E −1 ), p is atmospheric pressure, and P ref is the reference atmospheric pressure (1.01 × 10 5 Pa).
For C 4 plants, A n is calculated using a similar equation to Eq. (B22) according to Collatz et al. (1992).The rate of carboxylation is calculated by the minimum of the rates limited by light, maximum carboxylation velocity, and CO 2 as shown in the following: The solution of the Eqs.(B22-B26) yields net (A n ) and gross photosynthesis rates for a thin canopy layer with incident PAR flux Q per unit leaf area.We now solve for the photosynthesis integrated through the depth of a cohort's canopy, given incident PAR flux Q calculated according to the twostream approximation described in Sect. 2. Q is assumed to decrease exponentially, according to Beer's law, through the depth of a cohort's canopy: Q(l) = Q 0 exp(−κl), where Q 0 is incident PAR at the top of the cohort's canopy (obtained from the two-stream approximation in Sect.2) and l is the overlying leaf area per crown area at a given depth within the cohort's canopy, with L = 0 at the top of the cohort's crown, and l = LAI at the bottom (here, "LAI" is the total leaf area per crown area of a cohort's canopy).The Beer's law extinction coefficient κ is calculated as a function of the zenith angle of solar radiation (which varies by latitude, time of day, and day of year) and leaf angle distribution in the canopy (assumed spherical) to approximate the attenuation of photosynthetically active radiation within a single cohort's canopy according to the two-stream approximation described in Sect.B2 of this appendix.We can define a depth l eq where the light-limited rate J E is equal to the minimum of other limiting rates.Gross photosynthesis below depth l eq (the integral in Eq.B27 below) is a function of light availability, while above this depth it is equal to the minimum of other limiting rates.The net photosynthesis averaged over the entire canopy depth can be expressed as where If incident light Q 0 is so low that no part of canopy is lightsaturated, then l eq = 0.
Using the Beer's law approximation of the light profile within a cohort's canopy, we can obtain the following expressions for the integral in Eq. (B27): and for l eq , where α = α C i − * C i +2 * for C 3 plants and α = α for C 4 plants.Average stomatal conductance is calculated from Eqs. (B17) and (B27): where g s, min = 0.01 mol H 2 O m −2 s −1 is the minimum stomatal conductance allowed in the model.The model applies some further corrections to the net photosynthesis and stomatal conductance calculations above for a well-watered plant in order to take into account limitations imposed by water availability and other factors: where f l and f s are the fractions of canopy covered by liquid water and snow, respectively; α wet is the down-regulation coefficient, assumed to be 0.3; that is, photosynthesis of leaves fully covered by water or snow is reduced by 30 % compared to dry leaves.The model also imposes an upper limit on stomatal conductance.If the calculated g s is higher than the limit g max s = 0.25 mol m −2 s −1 , then stomatal conductance and net photosynthesis are adjusted: Finally, stomatal conductance and photosynthesis are adjusted down if available water supply is greater than water demand.Given mean stomatal conductance g s (Eq.B30), the water demand per individual (kg s −1 ) is where M air is the mass of air per mole (g mol −1 ), used to convert stomatal conductance to mass units, and A leaf is the total area of leaves in the individual's canopy.Given the water supply (i.e., the maximum plant water uptake rate, U max ; see Sect. 4 below), which is defined as the uptake rate when root water potential is at the plant permanent wilting point, net photosynthesis and stomatal conductance are adjusted for water limitation according to Eqs. (B31) and (B32) using the factor φ w = min(U max /U d , 1). (B36)

B4 Root water uptake and soil water dynamics
Calculations of water uptake by roots closely follow the model described in Milly et al. (2014), except in this study we consider multiple cohorts.Consequently, the Milly et al. (2014) formulation was expanded to allow roots from different cohorts to compete for water.We define maximum water uptake rate (or "water supply", U max ) as the amount of water an individual plant can potentially uptake from soil.Water demand (U d , Eq. B35) is the amount of water needed for non-water-limited photosynthesis, and uptake is the amount of water the plant actually gets.If supply (U max ) is greater than demand (U d ), then the plant is not waterlimited, and uptake will equal demand.If supply is less than demand, then the plant is water-limited, and uptake will be equal to supply.U max is calculated following Darcy's law, with a two-dimensional radial flow formulation in the approximation of quasi-steady flow in a small vicinity of fine roots.We use this model to derive an expression for water uptake as a function of xylem water potential.Setting xylem water potential equal to the plant permanent wilting point yields the value of U max needed for Eq.(B36).In the following, u is the water uptake rate per unit length of fine root (kg m −1 s −1 ) at a given soil depth, R the characteristic radial half-distance between fine roots (m), r r the root radius (m), and r the distance from the root axis (m).
For steady flow toward the root, where K(ψ) is unsaturated hydraulic conductivity (kg m −2 s −1 ), where K s is the conductivity of saturated soil, b is an empirical coefficient, ψis the soil water potential (m), and ψ * is the air entry water potential.Note that, since the flow is assumed to be in steady state, u does not depend on r; that is, ψ and thus K(ψ) are functions of r such that u(r) (Eq.B37) is constant (Gardner, 1960).Integrating from the root-soil interface (i.e., the root surface) to the half-distance R between fine roots (with potential ψ s at that distance from the root axis, and ψ r at the root surface): where ψ r is water potential at the root surface and ψ s that at distance R from the root.The macroscale water movement in the soil, and, consequently, water potential ψ s is calculated as in the LM3.0 model (see Milly et al., 2014, for details).
Eq. (B39) can be rewritten in a more convenient form: This relationship is assumed to hold at a macroscopic point, i.e., a soil layer at a given vertical depth in our case.The integral on the RHS of Eq. (B40) is sometimes called matric flux potential (Raats, 2007).The water flux through the root surface (i.e., membrane of surface cells) per unit length of root can also be expressed as where K r is permeability of root membrane per unit membrane area (kg m −2 area m −1 water potential gradient s −1 = kg m −3 s −1 ), and ψ x is root xylem water potential (m).The Supplement related to this article is available online at doi:10.5194/bg-12-2655-2015-supplement.

Figure 1 .
Figure 1.Scaling of mortality rates (a) and tree height and crown area (b) with DBH in four alternative versions of the LM3-PPA model (H0-H3; see Table2).In panel (a), the solid line shows the mortality rates of understory trees (same for H0-H3); the dashed line shows the mortality rates of canopy trees in H0, H2, and H3; and the dotted line is for canopy trees in H1.In panel (b), the solid line is tree height (same for H0-H3); the dashed line shows crown area in H0, H1, and H3; and the dotted line is crown area in H2.

Figure 2 .
Figure 2. GPP, NPP, allocation, and DBH growth rate.Panel (a) shows GPP (closed circles) and NPP (open circles) simulated by LM3-PPA for one species (sugar maple) in the 1 • × 1 • grid cell containing the Willow Creek AmeriFlux site in Wisconsin, USA.The red open circles with error bars are GPP estimates from the Willow Creek eddy flux data (Desai et al., 2005).The red open diamond is NPP estimated from biometric data at Willow Creek (Curtis et al., 2002).Panel (b) shows the simulated allocation of NPP to leaves, fine roots, woody tissues (including stems, branches, and coarse roots), and seeds.The green open circle, red open triangle, and black open circle are NPP of wood, fine roots, and leaves, respectively, estimated from biometric data(Curtis et al., 2002).Panel (c) shows the DBH growth rates of canopy trees (closed circles) and understory trees (open circles) simulated for sugar maple.The red circle and diamond show growth rates of canopy and understory trees for sugar maple in the northern lake states, USA, estimated from FIA forest inventory data(Zhang et al., 2014).The error bars represent 1 standard deviation.

Figure 3 .
Figure 3. Simulated vs. observed DBH growth rates of three tree species in the upper canopy and the understory.Circles, triangles, and diamonds are for Populus tremuloides, Acer saccharum, and A. rubrum, respectively.Closed and open symbols are for uppercanopy ("Top") and understory ("Under") trees, respectively.The FIA data used to estimate observed growth rates are from the northern lake states (Michigan, Wisconsin, and Minnesota), USA.Canopy growth rates were estimated by combining trees with a reported crown class of "dominant" or "co-dominant", and understory growth rates were estimated from trees with a crown class of "overtopped"(Zhang et al., 2014).

Figure 4 .
Figure 4. Forest succession.Panel (a) shows simulated forest succession for three species (Populus tremuloides, Acer saccharum, and A. rubrum), with parameters and initial densities in Table 1.Panel (b) shows the successional dynamics estimated from FIA inventory data in the northern lake states, USA.The basal areas of the three species are normalized relative to the maximum of their summed basal areas because the three species in the model runs account for only approximately one-half of the total basal area in the data.This normalization only changes the y-axis scale.The nonnormalized predictions and data are in Fig. S5.

Figure 5 .
Figure 5. Distributions of tree size (a) and biomass (b) in different stand age classes.Black symbols with dashed lines are from the FIA data of the northern lake states, USA, and blue symbols with solid lines are from the three-species LM3-PPA simulations in Fig. 4a.

Figure 6 .
Figure 6.Simulated distributions of tree size and biomass at quasiequilibrium in one-species (Acer saccharum) LM3-PPA simulations under alternative models assumptions (H0-H3).Size and biomass distributions are averaged over the last 400 years of 1000-year simulations.(a) Tree density of trees in 10 cm DBH bins.(b) Total tree density and basal area, summed over the size distribution in panel (a).The error bars represent one standard deviation.(c) Woody biomass in 10 cm DBH bins.Different colors in the figure refer to differ alternative model assumptions (see Table 2 and Fig. 1 for details): H0 is the baseline LM3-PPA model; H1 assumes that mortality rate increases with size for large trees; H2 assumes a maximum individual crown area, which causes a decline in DBH growth rate for large trees; and H3 assumes that open canopy space is filled by randomly chosen understory trees rather than the tallest understory trees as in the PPA model.

Figure 7 .
Figure 7. Simulated dynamics of biomass (a), soil carbon (b), and biomass turnover rate (c) in LM3-PPA and in a standard biogeochemical cycle (BGC) model that does not represent individuallevel processes.LM3-PPA was simulated with either one species (Acer saccharum) or all three species in Table1and Fig.4a.The standard BGC model is summarized in Fig.S1b.

Figure 8 .
Figure 8. Competition among PFTs that differ only their allocation to fine roots.Competition experiments were performed at two atmospheric CO 2 concentrations: 280 (a) and 560 ppm (b).In each experiment, a simulation was initialized with equal seedling densities of five PFTs that differed only in their ratio of fine-root area to leaf area (ϕ RL ).Because all PFTs shared the same target crown LAI, ϕ RL primarily determines allocation to fine roots and wood.

Figure 9 .
Figure 9. DBH growth rates of residents and invaders in pairwise invasion simulations.This figure shows DBH growth rates in pairwise competition experiments at (a) preindustrial [CO 2 ] (280 ppm) and (b) doubled [CO 2 ] (560 ppm) for residents (black bars) and invaders (gray bars) that differed only in their fine-root allocation (ϕ RL ; see Fig. 8 legend for explanation).In each experiment, the resident type was simulated for 400 years in monoculture, and then a small fraction of its density was converted to the invading type.The competitive optimum (ϕ RL = 0.7 and ϕ RL = 0.9 at 280 and 560 ppm, respectively) is the type (ϕ RL ) that cannot be invaded and can invade all other types (i.e., the convergence-stable evolutionarily stable strategy, ESS).

Figure 10 .
Figure 10.Woody NPP and DBH growth rates in monoculture and polyculture models runs at (a-b) [CO 2 ] = 280 ppm, and (cd) [CO 2 ] = 560 ppm.PFTs differed only their allocation to fine roots (ϕ RL ; see Fig.8legend for explanation).The optimal monoculture is defined as the type with the highest woody NPP (which, given the allometries in LM3-PPA, is also the type with the highest DBH growth rate) when grown in monoculture.In this figure, the competitive optimum is identified as the type with the highest woody NPP (or highest DBH growth rate) in polyculture model runs.Figures 8-10 present multiple ways to identify the competitive optimum (i.e., the convergence-stable ESS), and all yield consistent results: ϕ RL = 0.7 and ϕ RL = 0.9 at 280 and 560 ppm, respectively.

Figure 11 .
Figure 11.Changes in wet period length, water use efficiency (WUE), hydrological fluxes, and soil moisture due to a doubling of atmospheric CO 2 concentration.The bars show the percentage differences between LM3-PPA run with a single PFT (ϕ RL = 0.7) at [CO 2 ] = 560 ppm and at preindustrial [CO 2 ] (280 ppm).The "Wet season" bar shows the effect of a doubling of preindustrial [CO 2 ] on the fraction of each growing season in which canopy trees in the monoculture simulation are water-saturated (defined as the fraction of days during the growing season in which water supply was greater than or equal to demand at 14:00 over the final 60 years of a 500-year run).The "WUE" bar shows the change in the water use efficiency of canopy trees during the water-limited period (days in which water supply was less than demand at 14:00).The "Transp" and "Evap + Runoff" bars show the changes in water transpired by plants and lost via evaporation and runoff over the last 60 years of the model runs; when expressed in absolute amounts (mm yr −1 ), the decrease in transpiration and the increase in evaporation plus runoff almost exactly cancel each other out (see Sect. 3.3 in the Results).The "Soil moisture" bar shows change in growing-season mean soil moisture at doubled [CO 2 ].

Figure B1 .
Figure B1.Propagation of solar radiation through canopy layers.Each layer k is comprised of one or more cohorts according to the PPA model.I ↑ and I ↓ are the upward and downward fluxes of diffuse radiation, F is direct solar beam radiation flux, and S ↑ and S ↓ are upward and downward diffuse radiation fluxes due to canopy scattering of direct beam radiation.

Table 1 .
Parameter values for the three tree species in the LM3-PPA simulations presented in Figs.2-7.Initial densities in Table1are approximate and are summed across size classes.See TableC4in Appendix C for details of the initial size distributions used in the simulations. *

Table 3 .
Experimental design for model runs used to identify fineroot allocation strategies that are competitively optimal (evolutionarily stable strategies, ESSs) and that maximize wood production in monoculture.In these experiments, the plant functional types (PFTs) varied only in the ratio of fine-root surface area to leaf area (ϕ RL ).Because all PFTs shared the same target crown LAI, the parameter ϕ RL primarily controls allocation to fine roots and wood.Predicted diurnal and seasonal patterns of GPP, NPP, and evapotranspiration by the model are shown in Fig.
B18)where α c = 0.21 is the maximum ratio of oxygenation to carboxylation, [O 2 ] is the concentration of oxygen in canopy air (0.209 mol O 2 mol −1 air), and K C (mol CO 2 mol −1 air) and K O (mol O 2 mol −1 air) are the Michaelis-Menten constants for CO 2 and O 2 , respectively.K C and K O depend on temperature according to an Arrhenius function:

2694, 2015 2684 E. S. Weng et al.: LM3-PPA model leaves
, and φ m is the imposed maximum conductance limitation.If there is water or snow on the canopy, the photosynthesis is reduced proportionally to the covered fraction of leaves: B31)g s = φ w φ i φ m g s , (B32)where φ w is the reduction due to water limitations, φ i is reduction due to presence of intercepted water and snow on www.biogeosciences.net/12/2655/2015/Biogeosciences, 12, 2655-

Weng et al.: LM3-PPA model Appendix C: Tables for model variables, parameters, and initial plant sizes and density distributionsTable C1 .
Plant growth, respiration, and vegetation structure.

Table C1 .
Continued.G L + FR carbon used for the kg C tree −1 day −1 growth of leaves and fine roots G L carbon used for the growth of leaves kg C tree −1 day −1 G FR carbon used for the growth of fine roots kg C tree −1 day −1 G W + F carbon used for the growth kg C tree −1 day −1 of woody tissues and seed production G W carbon used for the growth of woody tissues kg C tree −1 day −1 www.biogeosciences.net/12/2655/2015/Biogeosciences, 12, 2655-2694, 2015

Table C3 .
Parameters of the three species (succession tests).

Table C4 .
Size distribution of initial trees in the succession test.