2.28. Dynamic Global Vegetation¶
2.28.1. What has changed¶
Deprecation of the dynamic global vegetation model (DGVM): The CLM5.0 model contains the legacy ‘CNDV’ code, which runs the CLM biogeochemistry model in combination with the LPJderived dynamics vegetation model introduced in CLM3. While this capacity has not technically been removed from the model, the DGVM has not been tested in the development of CLM5 and is no longer scientifically supported.
Introduction of FATES: The Functionally Assembled Terrestrial Ecosystem Simulator (FATES) is the actively developed DGVM for the CLM5.
2.29. Technical Documentation for FATES¶
FATES is the “Functionally Assembled Terrestrial Ecosystem Simulator”. It is an external module which can run within a given “Host Land Model” (HLM). Currently (November 2017) implementations are supported in both the Community Land Model(CLM) and in the land model of the E3SM Dept. of Energy Earth System Model.
FATES was derived from the CLM Ecosystem Demography model (CLM(ED)), which was documented in:
Fisher, R. A., Muszala, S., Verteinstein, M., Lawrence, P., Xu, C., McDowell, N. G., Knox, R. G., Koven, C., Holm, J., Rogers, B. M., Spessa, A., Lawrence, D., and Bonan, G.: Taking off the training wheels: the properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED), Geosci. Model Dev., 8, 35933619, https://doi.org/10.5194/gmd835932015, 2015.
and this technical note was first published as an appendix to that paper.
https://pdfs.semanticscholar.org/396c/b9f172cb681421ed78325a2237bfb428eece.pdf
2.29.1. Introduction¶
The Ecosystem Demography (‘ED’), concept within FATES is derived from the work of Moorcroft et al. (2001)
and is a cohort model of vegetation competition and coexistence, allowing a representation of the biosphere which accounts for the division of the land surface into successional stages, and for competition for light between height structured cohorts of representative trees of various plant functional types.
The implementation of the Ecosystem Demography concept within FATES links the surface flux and canopy physiology concepts in the CLM/E3SM with numerous additional developments necessary to accommodate the new model also documented here. These include a version of the SPITFIRE (Spread and InTensity of Fire) model of Thonicke et al. (2010), and an adoption of the concept of Perfect Plasticity Approximation approach of Purves et al. 2008, Lichstein et al. 2011 and Weng et al. 2014, in accounting for the spatial arrangement of crowns. Novel algorithms accounting for the fragmentation of coarse woody debris into chemical litter streams, for the physiological optimisation of canopy thickness, for the accumulation of seeds in the seed bank, for multilayer multiPFT radiation transfer, for droughtdeciduous and colddeciduous phenology, for carbon storage allocation, and for tree mortality under carbon stress, are also included and presented here.
Numerous other implementations of the Ecosystem Demography concept exist (See Fisher et al. (2018) for a review of these) Therefore, to avoid confusion between the concept of ‘Ecosystem Demography’ and the implementation of this concept in different models, the CLM(ED) implementation described by Fisher et al. (2015) will hereafter be called ‘FATES’ (the Functionally Assembled Terrestrial Ecosystem Simulator).
2.29.2. The representation of ecosystem heterogeneity in FATES¶
The terrestrial surface of the Earth is heterogeneous for many reasons, driven by variations in climate, edaphic history, ecological variability, geological forcing and human interventions. Land surface models represent this variability first by introducing a grid structure to the land surface, allowing different atmospheric forcings to operate in each grid cell, and subsequently by representing ‘subgrid’ variability in the surface properties. In the CLM, the land surface is divided into numerous ‘landunits’ corresponding to the underlying condition of the surface (e.g. soils, ice, lakes, bare ground) and then ‘columns’ referring to elements of the surface that share below ground resources (water & nutrients). Within the soil landunit, for example, there are separate columns for crops, and for natural vegetation, as these are assumed to use separate resource pools. The FATES model at present only operates on the naturally vegetated column. The soil column is subdivided into numerous tiles, that correspond to statistical fractions of the potentially vegetated land area. In the CLM 4.5 (and all previous versions of the model), subgrid tiling operates on the basis of plant functional types (PFTs). That is, each piece of land is assumed to be occupied by only one plant functional type, with multiple PFTspecific tiles sharing a common soil water and nutrient pool. This PFTbased tiling structure is the standard method used by most land surface models deployed in climate prediction.
The introduction of the Ecosystem Demography concept introduces significant alterations to the representation of the land surface in the CLM. In FATES, the tiling structure represents the disturbance history of the ecosystem. Thus, some fraction of the land surface is characterized as ‘recently disturbed’, some fraction has escaped disturbance for a long time, and other areas will have intermediate disturbances. Thus the ED concept essentially discretizes the trajectory of succession from disturbed ground to ‘mature’ ecosystems. Within FATES, each “disturbance history class” is referred to as a ‘patch’. The word “patch” has many possible interpretations, so it is important to note that: there is no spatial location associated with the concept of a ‘patch’ . It refers to a fraction of the potential vegetated area consisting of all parts of the ecosystem with similar disturbance history.
The ‘patch’ organizational structure in CLM thus replaces the previous ‘PFT’ structure in the organization heirarchy. The original hierarchical land surface organizational structure of CLM as described in Oleson et al. 2013 may be depicted as:
and the new structure is altered to the following:
Thus, each gridcell becomes a matrix of ‘patches’ that are conceptualized by their ‘age since disturbance’ in years. This is the equivalent of grouping together all those areas of a gridcell that are ‘canopy gaps’, into a single entity, and all those areas that are ‘mature forest’ into a single entity.
2.29.2.1. Cohortized representation of tree populations¶
Each commondisturbancehistory patch is a notional ecosystem that might in reality contain numerous individual plants which vary in their physiological attributes, in height and in spatial position. One way of addressing this heterogeneity is to simulate a forest of specific individuals, and to monitor their behavior through time. This is the approach taken by “gap” and individualbased models (Smith et al. 2001, Sato et al. 2007, Uriarte et al. 2009, Fyllas et al. 2014). The depiction of individuals typically implies that the outcome of the model is stochastic. This is because we lack the necessary detailed knowledge to simulate the individual plant’s fates. Thus gap models imply both stochastic locations and mortality of plants. Thus, (with a genuinely random seed) each model outcome is different, and an ensemble of model runs is required to generate an average representative solution. Because the random death of large individual trees can cause significant deviations from the mean trajectory for a small plot (a typical simulated plot size is 30m x 30 m) the number of runs required to minimize these deviations is large and computationally expensive. For this reason, models that resolve individual trees typically use a physiological timestep of one day or longer (e.g. Smith et al. 2001, Xiaidong et al. 2005, Sato et al. 2007
The approach introduced by the Ecosystem Demography model Moorcroft et al. 2001 is to group the hypothetical population of plants into “cohorts”. In the notional ecosystem, after the landsurface is divided into commondisturbancehistory patches, the population in each patch is divided first into plant functional types (the standard approach to representing plant diversity in large scale vegetation models), and then each plant type is represented as numerous height classes. Importantly, for each PFT/height class bin, we model *one* representative individual plant, which tracks the average properties of this `cohort` of individual plants. Thus, each commondisturbancehistory patch is typically occupied by a set of cohorts of different plant functional types, and different height classes within those plant functional types. Each cohort is associated with a number of identical trees, (where denotes the identification or index number for a given cohort)..
The complete hierarchy of elements in FATES is therefore now described as follows:
2.29.2.2. Discretization of cohorts and patches¶
Newly disturbed land and newly recruited seedlings can in theory be generated at each new model timestep as the result of germination and disturbance processes. If the new patches and cohorts established at every timestep were tracked by the model structure, the computational load would of course be extremely high (and thus equivalent to an individualbased approach). A signature feature of the ED model is the system by which functionally equivalent patches and cohorts are fused into single model entities to save memory and computational time.
This functionality requires that criteria are established for the meaning of functional equivalence, which are by necessity slightly subjective, as they represent ways of abstracting reality into a more tractable mathematical representation. As an example of this, for heightstructured cohorts, we calculate the relativized differences in height (, m) between two cohorts of the same pft, and as
If is smaller than some threshold , and they are of the same plant functional type, the two cohorts are considered equivalent and merged to form a third cohort , with the properties of cohort and averaged such that they conserve mass. The model parameter can be adjusted to adjust the tradeoff between simulation accuracy and computational load. There is no theoretical optimal value for this threshold but it may be altered to have finer or coarser model resolutions as needed.
Similarly, for commondisturbancehistory patches, we again assign a threshold criteria, which is then compared to the difference between patches and , and if the difference is less than some threshold value () then patches are merged together, otherwise they are kept separate. However, in contrast with heightstructured cohorts, where the meaning of the difference criteria is relatively clear, how the landscape should be divided into commondisturbancehistory units is less clear. Several alternative criteria are possible, including Leaf Area Index, total biomass and total stem basal area.
In this implementation of FATES we assess the amount of aboveground biomass in each PFT/plant diameter bin. Biomass is first grouped into fixed diameter bins for each PFT () and a significant difference in any bin will cause patches to remain separated. This means that if two patches have similar total biomass, but differ in the distribution of that biomass between diameter classes or plant types, they remain as separate entities. Thus
is the binned aboveground biomass profile for patch , is the diameter class. and are the lower and upper boundaries for the diameter class. and depict the biomass (KgC m) and the number of individuals of each cohort respectively. A difference matrix between patches and is thus calculated as
If all the values of are smaller than the threshold, , then the patches and are fused together to form a new patch .
To increase computational efficiency and to simplify the coding structure of the model, the maximum number of patches is capped at . To force the fusion of patches down to this number, the simulation begins with a relatively sensitive discretization of patches ( = 0.2) but if the patch number exceeds the maximum, the fusion routine is repeated iteratively until the two most similar patches reach their fusion threshold. This approach maintains an even discretization along the biomass gradient, in contrast to, for example, simply fusing the oldest or youngest patches together.
The area of the new patch (, m) is the sum of the area of the two existing patches,
and the cohorts ‘belonging’ to patches and now cooccupy patch . The state properties of and (litter, seed pools, etc. ) are also averaged in accordance with mass conservation .
2.29.2.3. Linked Lists: the general code structure of FATES¶
The number of patches in each natural vegetation column and the number of cohorts in any given patch are variable through time because they are recalculated for each daily timestep of the model. The more complex an ecosystem, the larger the number of patches and cohorts. For a slowly growing ecosystem, where maximum cohort size achieved between disturbance intervals is low, the number of cohorts is also low. For fastgrowing ecosystems where many plant types are viable and maximum heights are large, more cohorts are required to represent the ecosystem with adequate complexity.
In terms of variable structure, the creation of an array whose size could accommodate every possible cohort would mean defining the maximum potential number of cohorts for every potential patch, which would result in very large amounts of wasted allocated memory, on account of the heterogeneity in the number of cohorts between complex and simple ecosystems (n.b. this does still happen for some variables at restart timesteps). To resolve this, the cohort structure in FATES model does not use an array system for internal calculations. Instead it uses a system of linked lists where each cohort structure is linked to the cohorts larger than and smaller than itself using a system of pointers. The shortest cohort in each patch has a ‘shorter’ pointer that points to the null value, and the tallest cohort has a ‘taller’ pointer that points to the null value.
Instead of iterating along a vector indexed by , the code structures typically begin at the tallest cohort in a given patch, and iterate until a null pointer is encountered.
Using this structure, it is therefore possible to have an unbounded upper limit on cohort number, and also to easily alter the ordering of cohorts if, for example, a cohort of one functional type begins to grow faster than a competitor of another functional type, and the cohort list can easily be reordered by altering the pointer structure. Each cohort has pointers indicating to which patch and gridcell it belongs. The patch system is analogous to the cohort system, except that patches are ordered in terms of their relative age, with pointers to older and younger patches where cp is the oldest:
2.29.2.4. Indices used in FATES¶
Some of the indices used in FATES are similar to those used in the standard CLM4.5 model; column (), land unit(), grid cell() and soil layer (). On account of the additional complexity of the new representation of plant function, several additional indices are introduced that describe the discritization of plant type, fuel type, litter type, plant height, canopy identity, leaf vertical structure and fuel moisture characteristics. To provide a reference with which to interpret the equations that follow, they are listed here.
Parameter Symbol 
Parameter Name 

ft 
Plant Functional Type 
fc 
Fuel Class 
lsc 
Litter Size Class 
coh 
Cohort Index 
patch 
Patch Index 
Cl 
Canopy Layer 
z 
Leaf Layer 
mc 
Moisture Class 
2.29.2.5. Cohort State Variables¶
The unit of allometry in the ED model is the cohort. Each cohort represents a group of plants with similar functional types and heights that occupy portions of column with similar disturbance histories. The state variables of each cohort therefore consist of several pieces of information that fully describe the growth status of the plant and its position in the ecosystem structure, and from which the model can be restarted. The state variables of a cohort are as follows:
Quantity 
Variable name 
Units 
Notes 

Plant Functional Type 
integer 

Number of Individuals 
n per 10000m:math:` ^{2}` 

Height 
m 

Diameter 
cm 

Structural Biomass 
KgC plant 
Stem wood (above and below ground) 

Alive Biomass 
KgC plant 
Leaf, fine root and sapwood 

Stored Biomass 
KgC plant 
Labile carbon reserve 

Leaf memory 
KgC plant 
Leaf mass when leaves are dropped 

Canopy Layer 
integer 
1 = top layer 

Phenological Status 
integer 
1=leaves off. 2=leaves on 

Canopy trimming 
fraction 
1.0=max leaf area 

Patch Index 
integer 
To which patch does this cohort belong? 
2.29.2.6. Patch State Variables¶
A patch, as discuss earlier, is a fraction of the landscape which contains ecosystems with similar structure and disturbance history. A patch has no spatial location. The state variables, which are ‘ecosystem’ rather than ‘tree’ scale properties, from which the model can be restarted, are as follows
Quantity 
Variable name 
Units 
Indexed By 


Area 
m 

Age 
years 

Seed 
KgC m 

Leaf Litter 
KgC m 

Root Litter 
KgC m 

AG Coarse Woody Debris 
:math:` {CWD}_{A G,patch}` 
KgC m 
Size Class (lsc) 

BG Coarse Woody Debris 
:math:` {CWD}_{B G,patch}` 
KgC m 
Size Class (lsc) 

Canopy Spread 
Canopy Layer 

Column Index 
integer 
2.29.2.7. Model Structure¶
Code concerned with the Ecosystem Demography model interfaces with the CLM model in four ways: i) During initialization, ii) During the calculation of surface processes (albedo, radiation absorption, canopy fluxes) each model time step (typically halfhourly), iii) During the main invokation of the ED model code at the end of each day. Daily cohortlevel NPP is used to grow plants and alter the cohort structures, disturbance processes (fire and mortality) operate to alter the patch structures, and all fragmenting carbon pool dynamics are calculated. iv) during restart reading and writing. The net assimilation (NPP) fluxes attributed to each cohort are accumulated throughout each daily cycle and passed into the ED code as the major driver of vegetation dynamics.
2.29.3. Initialization of vegetation from bare ground¶
If the model is restarted from a bare ground state (as opposed to a preexisting vegetation state), the state variables above are initialized as follows. First, the number of plants per PFT is allocated according to the initial seeding density (, individuals per m) and the area of the patch , which in the first timestep is the same as the area of the notional ecosystem . The model has no meaningful spatial dimension, but we assign a notional area such that the values of ‘’ can be attributed. The default value of is one hectare (10,000 m), but the model will behave identically irrepective of the value of this parameter.
Each cohort is initialized at the minimum canopy height , which is specified as a parameter for each plant functional type and denotes the smallest size of plant which is tracked by the model. Smaller plants are not considered, and their emergence from the recruitment processes is unresolved and therefore implicitly parameterized in the seedling establishment model.. The diameter of each cohort is then specified using the loglinear allometry between stem diameter and canopy height
where the slope of the loglog relationship, is 0.64 and the intercept is 0.37. The structural biomass associated with a plant of this diameter and height is given (as a function of wood density, , g cm)
taken from the original ED1.0 allometry Moorcroft et al. 2001 (values of the allometric constants in Table [table:allom]. The maximum amount of leaf biomass associated with this diameter of tree is calculated according to the following allometry
from this quantity, we calculate the active/fine root biomass as
where is the fraction of fine root biomass to leaf biomass, assigned per PFT
Parameter Symbol 
Parameter Name 
Units 
Default Value 

Minimum plant height 
m 
1.5 

Initial Planting density 
Individuals m 

Model area 
m 
10,000 
[table:init]
2.29.4. Allocation of biomass¶
Total live biomass is the state variable of the model that describes the sum of the three live biomass pools leaf , root and sapwood (all in kGC individual). The quantities are constrained by the following
Sapwood volume is a function of tree height and leaf biomass
where is the ratio of sapwood mass (kgC) to leaf mass per unit tree height (m). Also, root mass is a function of leaf mass
Thus
Rearranging gives the fraction of biomass in the leaf pool as
Thus, we can determine the leaf fraction from the height at the tissue ratios, and the phenological status of the cohort .
To divide the live biomass pool at restart, or whenever it is recalculated, into its consituent parts, we first
Because sometimes the leaves are dropped, using leaf biomass as a predictor of root and sapwood would produce zero live biomass in the winter. To account for this, we add the LAI memory variable to the live biomass pool to account for the need to maintain root biomass when leaf biomass is zero. Thus, to calculated the root biomass, we use
To calculated the sapwood biomass, we use
Parameter Symbol 
Parameter Name 
Units 
Default Value 

Allometry intercept 
0.37 

Allometry slope 
0.64 

Structural biomass multiplier 
0.06896 

Structural Biomass dbh exponent 
1.94 

Structural Biomass height exponent 
0.572 

Structural Biomass density exponent 
0.931 

Leaf biomass multiplier 
0.0419 

Leaf biomass dbh exponent 
1.56 

Leaf biomass density exponent 
0.55 

Ratio of sapwood mass to height 
m 

Ratio of fine root mass to leaf mass 
1.0 
[table:allom]
2.29.5. Canopy Structure and the Perfect Plasticity Approximation¶
During initialization and every subsequent daily ED timestep, the canopy structure model is called to determine how the leaf area of the different cohorts is arranged relative to the incoming radiation, which will then be used to drive the radiation and photosynthesis calculations. This task requires that some assumptions are made about 1) the shape and depth of the canopy within which the plant leaves are arranged and 2) how the leaves of different cohorts are arranged relative to each other. This set of assumptions are critical to model performance in EDlike cohort based models, since they determine how light resources are partitioned between competing plants of varying heights, which has a very significant impact on how vegetation distribution emerges from competition Fisher et al. 2010.
The standard ED1.0 model makes a simple ‘flat disk’ assumption, that the leaf area of each cohort is spread in an homogenous layer at one exact height across entire the ground area represented by each patch. FATES has diverged from this representation due to (at least) two problematic emergent properties that we identified as generating unrealistic behaviours espetially for largearea patches.
1. Overestimation of light competition . The vertical stacking of cohorts which have all their leaf area at the same nominal height means that when one cohort is only very slightly taller than it’s competitor, it is completely shaded by it. This means that any small advantage in terms of height growth translates into a large advantage in terms of light competition, even at the seedling stage. This property of the model artificially exaggerates the process of light competition. In reality, trees do not compete for light until their canopies begin to overlap and canopy closure is approached.
2. Unrealistic overcrowding. The ‘flatdisk’ assumption has no consideration of the spatial extent of tree crowns. Therefore it has no control on the packing density of plants in the model. Given a mismatch between production and mortality, entirely unrealistic tree densities are thus possible for some combinations of recruitment, growth and mortality rates.
To account for the filling of space in three dimensions using the onedimensional representation of the canopy employed by CLM, we implement a new scheme derived from that of Purves et al. 2008. Their argument follows the development of an individualbased variant of the SORTIE model, called SHELL, which allows the location of individual plant crowns to be highly flexible in space. Ultimately, the solutions of this model possess an emergent property whereby the crowns of the plants simply fill all of the available space in the canopy before forming a distinct understorey.
Purves et al. developed a model that uses this feature, called the ‘perfect plasticity approximation’, which assumes the plants are able to perfectly fill all of the available canopy space. That is, at canopy closure, all of the available horizontal space is filled, with negligible gaps, owing to lateral tree growth and the ability of tree canopies to grow into the available gaps (this is of course, an oversimplified but potential useful ecosystem property). The ‘perfect plasticity approximation’ (PPA) implies that the community of trees is subdivided into discrete canopy layers, and by extension, each cohort represented by FATES model is assigned a canopy layer status flag, . In this version, we set the maximum number of canopy layers at 2 for simplicity, although is possible to have a larger number of layers in theory. = 1 means that all the trees of cohort are in the upper canopy (overstory), and = 2 means that all the trees of cohort are in the understorey.
In this model, all the trees in the canopy experience full light on their uppermost leaf layer, and all trees in the understorey experience the same light (full sunlight attenuated by the average LAI of the upper canopy) on their uppermost leaves, as described in the radiation transfer section (more nuanced versions of this approach may be investigated in future model versions). The canopy is assumed to be cylindrical, the lower layers of which experience selfshading by the upper layers.
To determine whether a second canopy layer is required, the model needs to know the spatial extent of tree crowns. Crown area, , m, is defined as
where is the crown area of a single tree canopy (m:math:^{2}) and is the ‘canopy spread’ parameter (m cm^1) of this canopy layer, which is assigned as a function of canopy space filling, discussed below. In contrast to Purves et al. 2008 , we use an exponent, identical to that for leaf biomass, of 1.56, not 2.0, such that tree leaf area index does not change as a function of diameter.
To determine whether the canopy is closed, we calculate the total canopy area as:
where is the number of cohorts in a given patch. If the area of all crowns (m:math:^{2}) is larger than the total ground area of a patch (), then some fraction of each cohort is demoted to the understorey.
Under these circumstances, the extra crown area (i.e.,  ) is moved into the understorey. For each cohort already in the canopy, we determine a fraction of trees that are moved from the canopy () to the understorey. is calculated as Fisher et al. 2010
where is a weighting of each cohort determined by basal diameter (cm) and the competitive exclusion coefficient
The higher the value of the greater the impact of tree diameter on the probability of a given tree obtaining a position in the canopy layer. That is, for high values, competition is highly deterministic. The smaller the value of , the greater the influence of random factors on the competitive exclusion process, and the higher the probability that slower growing trees will get into the canopy. Appropriate values of are poorly constrained but alter the outcome of competitive processes.
The process by which trees are moved between canopy layers is complex because 1) the crown area predicted for a cohort to lose may be larger than the total crown area of the cohort, which requires iterative solutions, and 2) on some occasions (e.g. after fire), the canopy may open up and require ‘promotion’ of cohorts from the understorey, and 3) canopy area may change due to the variations of canopy spread values ( , see the section below for details) when fractions of cohorts are demoted or promoted. Further details can be found in the code references in the footnote.
2.29.5.1. Horizontal Canopy Spread¶
Purves et al. 2008 estimated the ratio between canopy and stem diameter as 0.1 m cm for canopy trees in North American forests, but this estimate was made on trees in closed canopies, whose shape is subject to space competition from other individuals. Sapling trees have no constraints in their horizontal spatial structure, and as such, are more likely to display their leaves to full sunlight. Also, prior to canopy closure, light interception by leaves on the sides of the canopy is also higher than it would be in a closed canopy forest. If the ‘canopy spread’ parameter is constant for all trees, then we simulate high levels of selfshading for plants in unclosed canopies, which is arguably unrealistic and can lower the productivity of trees in areas of unclosed canopy (e.g. low productivity areas of boreal or semiarid regions where LAI and canopy cover might naturally be low). We here interpret the degree of canopy spread, as a function of how much tree crowns interfere with each other in space, or the total canopy area . However itself is a function of , leading to a circularity. is thus solved iteratively through time.
Each daily model step, and the fraction of the gridcell occupied by tree canopies in each canopy layer ( = /) is calculated based on from the previous timestep. If is greater than a threshold value , is increased by a small increment . The threshold is, hypothetically, the canopy fraction at which light competition begins to impact on tree growth. This is less than 1.0 owing to the nonperfect spatial spacing of tree canopies. If is greater than , then is reduced by an increment , to reduce the spatial extent of the canopy, thus.
The values of are bounded to upper and lower limits. The lower limit corresponds to the observed canopy spread parameter for canopy trees and the upper limit corresponds to the largest canopy extent
This iterative scheme requires two additional parameters ( and ). affects the speed with which canopy spread (and hence leaf are index) increase as canopy closure is neared. However, the model is relatively insensitive to the choice of either or .
2.29.5.2. Definition of Leaf and Stem Area Profile¶
Within each patch, the model defines and tracks cohorts of multiple plant functional types that exist either in the canopy or understorey. Light on the top leaf surface of each cohort in the canopy is the same, and the rate of decay through the canopy is also the same for each PFT. Therefore, we accumulate all the cohorts of a given PFT together for the sake of the radiation and photosynthesis calculations (to avoid separate calculations for every cohort).
Therefore, the leaf area index for each patch is defined as a threedimensional array where is the canopy layer, is the functional type and is the leaf layer within each canopy. This threedimensional structure is the basis of the radiation and photosynthetic models. In addition to a leaf area profile matrix, we also define, for each patch, the area which is covered by leaves at each layer as .
Each plant cohort is already defined as a member of a single canopy layer and functional type. This means that to generate the matrix, it only remains to divide the leaf area of each cohort into leaf layers. First, we determine how many leaf layers are occupied by a single cohort, by calculating the ‘tree LAI’ as the total leaf area of each cohort divided by its crown area (both in m)
where is the specific leaf area in m KgC and is in kGC per plant.
Stem area index (SAI) is ratio of the total area of all woody stems on a plant to the area of ground covered by the plant. During winter in deciduous areas, the extra absorption by woody stems can have a significant impact on the surface energy budget. However, in previous big leaf versions of the CLM, computing the circumstances under which stem area was visible in the absence of leaves was difficult and the algorithm was largely heuristic as a result. Given the multilayer canopy introduced for FATES, we can determine the leaves in the higher canopy layers will likely shade stem area in the lower layers when leaves are on, and therefore stem area index can be calculated as a function of woody biomass directly.
Literature on stem area index is particularly poor, as it’s estimation is complex and not particularly amenable to the use of, for example, assumptions of random distribution in space that are typically used to calculate leaf area from light interception. Kucharik et al. 1998 estimated that SAI visible from an LAI2000 sensor was around 0.5 m^2 m^2. Low et al. 2001 estimate that the wood area index for Ponderosa Pine forest is 0.270.33. The existing CLM(CN) algorithm sets the minimum SAI at 0.25 to match MODIS observations, but then allows SAI to rise as a function of the LAI lost, meaning than in some places, predicted SAI can reach value of 8 or more. Clearly, greater scientific input on this quantity is badly needed. Here we determine that SAI is a linear function of woody biomass, to at very least provide a mechanistic link between the existence of wood and radiation absorbed by it. The nonlinearity between how much woody area exists and how much radiation is absorbed is provided by the radiation absorption algorithm. Specifically, the SAI of an individual cohort (, m m) is calculated as follows,
where is the coefficient linking structural biomass to SAI. The number of occupied leaf layers for cohort () is then equal to the rounded up integer value of the tree SAI () and LAI () divided by the layer thickness (i.e., the resolution of the canopy layer model, in units of vegetation index (+) with a default value of 1.0, ),
The fraction of each layer that is leaf (as opposed to stem) can then be calculated as
Finally, the leaf area in each leaf layer pertaining to this cohort is thus
and the stem area index is
where is the remainder of the canopy that is below the last full leaf layer
is the total canopy area occupied by plants in a given patch (m:math:^{2}) and is calculated as follows,
The canopy is conceived as a cylinder, although this assumption could be altered given sufficient evidence that canopy shape was an important determinant of competitive outcomes, and the area of ground covered by each leaf layer is the same through the cohort canopy. With the calculated SAI and LAI, we are able to calculate the complete canopy profile. Specifically, the relative canopy area for the cohort is calculated as
The total occupied canopy area for each canopy layer (), plant functional type () and leaf layer () bin is thus
where and
All of these quantities are summed across cohorts to give the complete leaf and stem area profiles,
2.29.5.3. Burial of leaf area by snow¶
The calculations above all pertain to the total leaf and stem area indices which charecterize the vegetation structure. In addition, the model must know when the vegetation is covered by snow, and by how much, so that the albedo and energy balance calculations can be adjusted accordingly. Therefore, we calculated a ‘total’ and ‘exposed’ and profile using a representation of the bottom and top canopy heights, and the depth of the average snow pack. For each leaf layer of each cohort, we calculate an ‘exposed fraction via consideration of the top and bottom heights of that layer and (m),
where is the plant functional type () specific fraction of the cohort height that is occupied by the crown. Specifically, the ‘exposed fraction is calculated as follows,
The resulting exposed () and total () leaf and stem area indicies are calculated as
and are used in the radiation interception and photosynthesis algorithms described later.
Parameter Symbol 
Parameter Name 
Units 
Notes 
Indexed by 

:math:` delta_ {vai}` 
Thickness of single canopy layer 
mm: math:^{2} 

Competitive Exclusion Parameter 
none 

Minimum canopy spread 
m cm:math:` ^{1}` 

Competitive Exclusion Parameter 
m cm:math:` ^{1}` 

Incremental change in 
m cm:math:` ^{1}` y 

Threshold canopy closure 
none 

Crown fraction 
none 

Stem area per unit woody biomass 
m^2 KgC^1 
2.29.6. Radiation Transfer¶
2.29.6.1. Fundamental Radiation Transfer Theory¶
The first interaction of the land surface with the properties of vegetation concerns the partitioning of energy into that which is absorbed by vegetation, reflected back into the atmosphere, and absorbed by the ground surface. Older versions of the CLM have utilized a “twostream” approximation Sellers 1985, Sellers et al. 1986 that provided an empirical solution for the radiation partitioning of a multilayer canopy for two streams, of diffuse and direct light. However, implementation of the Ecosystem Demography model requires a) the adoption of an explicit multiple layer canopy b) the implementation of a multiple plant type canopy and c) the distinction of canopy and understorey layers, inbetween which the radiation streams are fully mixed. The radiation mixing between canopy layers is necessary as the position of different plants in the understorey is not defined spatially or relative to the canopy trees above. In this new scheme, we thus implemented a onedimensional scheme that traces the absorption, transmittance and reflectance of each canopy layer and the soil, iterating the upwards and downwards passes of radiation through the canopy until a predefined accuracy tolerance is reached. This approach is based on the work of Norman 1979.
Here we describe the basic theory of the radiation transfer model for the case of a single homogenous canopy, and in the next section we discuss how this is applied to the multi layer multi PFT canopy in the FATES implementation. The code considers the fractions of a single unit of incoming direct and a single unit of incoming diffuse light, that are absorbed at each layer of the canopy for a given solar angle (, radians). Direct radiation is extinguished through the canopy according to the coefficient that is calculated from the incoming solar angle and the dimensionless leaf angle distribution parameter () as
where
and
The leaf angle distribution is a descriptor of how leaf surfaces are arranged in space. Values approaching 1.0 indicate that (on average) the majority of leaves are horizontally arranged with respect to the ground. Values approaching 1.0 indicate that leaves are mostly vertically arranged, and a value of 0.0 denotes a canopy where leaf angle is random (a ‘spherical’ distribution).
According to Beer’s Law, the fraction of light that is transferred through a single layer of vegetation (leaves or stems) of thickness , without being intercepted by any surface, is
and the incident direct radiation transmitted to each layer of the canopy () is thus calculated from the cumulative leaf area ( ) shading each layer ():
The fraction of the leaves that are exposed to direct light is also calculated from the decay coefficient .
where is the fraction of leaves that are shaded from direct radiation and only receive diffuse light.
Diffuse radiation, by definition, enters the canopy from a spectrum of potential incident directions, therefore the unintercepted transfer () through a leaf layer of thickness is calculated as the mean of the transfer rate from each of 9 different incident light directions () between 0 and 180 degrees to the horizontal.
The fraction (1) of the diffuse radiation is intercepted by leaves as it passes through each leaf layer. Of this, some fraction is reflected by the leaf surfaces and some is transmitted through. The fractions of diffuse radiation reflected from () and transmitted though () each layer of leaves are thus, respectively
where and are the fractions of incident light reflected and transmitted by individual leaf surfaces.
Once we know the fractions of light that are transmitted and reflected by each leaf layer, we begin the process of distributing light through the canopy. Starting with the first leaf layer (=1), where the incident downwards diffuse radiation () is 1.0, we work downwards for layers, calculating the radiation in the next layer down () as:
Here, calculates the fraction of incoming energy transmitted downwards onto layer . This flux is then increased by the additional radiation that is reflected upwards from further down in the canopy to layer , and then is reflected back downwards according to the reflected fraction . The more radiation in , the smaller the denominator and the larger the downwards flux. is also calculated sequentially, starting this time at the soil surface layer (where )
where is the soil albedo characteristic. The upwards reflected fraction for each leaf layer, moving upwards, is then Norman 1979
The corresponding upwards diffuse radiation flux is therefore the fraction of downwards radiation that is incident on a particular layer, multiplied by the fraction that is reflected from all the lower layers:
Now we have initial conditions for the upwards and downwards diffuse fluxes, these must be modified to account for the fact that, on interception with leaves, direct radiation is transformed into diffuse radiation. In addition, the initial solutions to the upwards and downwards radiation only allow a single ‘bounce’ of radiation through the canopy, so some radiation which might be intercepted by leaves higher up is potentially lost. Therefore, the solution to this model is iterative. The iterative solution has upwards and a downwards components that calculate the upwards and downwards fluxes of total radiation at each leaf layer ( and ) . The downwards component begins at the top canopy layer (). Here we define the incoming solar diffuse and direct radiation ( and respectively).
The first term of the righthand side deals with the diffuse radiation transmitted downwards, the second with the diffuse radiation travelling upwards, and the third with the direct radiation incoming at each layer () that is intercepted by leaves () and then transmitted through through the leaf matrix as diffuse radiation (). At the bottom of the canopy, the light reflected off the soil surface is calculated as
The upwards propagation of the reflected radiation is then
Here the first two terms deal with the diffuse downwards and upwards fluxes, as before, and the third deals direct beam light that is intercepted by leaves and reflected upwards. These upwards and downwards fluxes are computed for multiple iterations, and at each iteration, and are compared to their values in the previous iteration. The iteration scheme stops once the differences between iterations for all layers is below a predefined tolerance factor, (set here at ). Subsequently, the fractions of absorbed direct () and diffuse () radiation for each leaf layer then
and, the radiation energy absorbed by the soil for the diffuse and direct streams is is calculated as
Canopy level albedo is denoted as the upwards flux from the top leaf layer
and the division of absorbed energy into sunlit and shaded leaf fractions, (required by the photosynthesis calculations), is
2.29.6.2. Resolution of radiation transfer theory within the FATES canopy structure¶
The radiation transfer theory above, was described with reference to a single canopy of one plant functional type, for the sake of clarity of explanation. The FATES model, however, calculates radiative and photosynthetic fluxes for a more complex hierarchical structure within each patch/timesincedisturbance class, as described in the leaf area profile section. Firstly, we denote two or more canopy layers (denoted ). The concept of a ‘canopy layer’ refers to the idea that plants are organized into discrete over and understories, as predicted by the Perfect Plasticity Approximation (Purves et al. 2008, Fisher et al. 2010). Within each canopy layer there potentially exist multiple cohorts of different plant functional types and heights. Within each canopy layer, , and functional type, , the model resolves numerous leaf layers , and, for some processes, notably photosynthesis, each leaf layer is split into a fraction of sun and shade leaves, and , respectively.
The radiation scheme described in Section is solved explicitly for this structure, for both the visible and nearinfrared wavebands, according to the following assumptions.
A canopy layer () refers to either the over or understorey
A leaf layer () refers to the discretization of the LAI within the canopy of a given plant functional type.
All PFTs in the same canopy layer have the same solar radiation incident on the top layer of the canopy
Light is transmitted through the canopy of each plant functional type independently
Between canopy layers, the light streams from different plant functional types are mixed, such that the (undefined) spatial location of plants in lower canopy layers does not impact the amount of light received.
Where understorey layers fill less area than the overstorey layers, radiation is directly transferred to the soil surface.
All these calculations pertain to a single patch, so we omit the patch subscript for simplicity in the following discussion.
Within this framework, the majority of the terms in the radiative transfer scheme are calculated with indices of , and . In the following text, we revisit the simplified version of the radiation model described above, and explain how it is modified to account for the more complex canopy structure used by FATES.
Firstly, the light penetration functions, and are described as functions of , because the leaf angle distribution, , is a pftspecific parameter. Thus, the diffuse irradiance transfer rate, is also specific because , on which it depends, is a function of .
The amount of direct light reaching each leaf layer is a function of the leaves existing above the layer in question. If a leaf layer ‘’ is in the top canopy layer (the overstorey), it is only shaded by leaves of the same PFT so is unchanged from equation. If there is more than one canopy layer (), then the amount of direct light reaching the top leaf surfaces of the second/lower layer is the weighted average of the light attenuated by all the parallel tree canopies in the canopy layer above, thus.
where is the areal fraction of each canopy layer occupied by each functional type and is the index of the bottom canopy layer of each pft in each canopy layer (the subscripts and are implied but omitted from all references to avoid additional complications)
Similarly, the sunlit fraction for a leaf layer ‘’ in the second canopy layer (where ) is
where is the weighted average sunlit fraction in the bottom layer of a given canopy layer.
Following through the sequence of equations for the simple single pft and canopy layer approach above, the and fluxes are also indexed by , , and . The diffuse radiation reflectance ratio is also calculated in a manner that homogenizes fluxes between canopy layers. For the canopy layer nearest the soil ( = ). For the top canopy layer (=1), a weighted average reflectance from the lower layers is used as the baseline, in lieu of the soil albedo. Thus:
For the iterative flux resolution, the upwards and downwards fluxes are also averaged between canopy layers, thus where
and where =1, and
The remaining terms in the radiation calculations are all also indexed by , and so that the fraction of absorbed radiation outputs are termed and . The sunlit and shaded absorption rates are therefore
and
The albedo of the mixed pft canopy is calculated as the weighted average of the upwards radiation from the top leaf layer of each pft where =1:
The radiation absorbed by the soil after passing through through understorey vegetation is:
to which is added the diffuse flux coming directly from the upper canopy and hitting no understorey vegetation.
and the direct flux coming directly from the upper canopy and hitting no understorey vegetation.
These changes to the radiation code are designed to be structurally flexible, and the scheme may be collapsed down to only include on canopy layer, functional type and pft for testing if necessary.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Leaf angle distribution parameter 
none 
ft 

Fraction of light reflected by leaf surface 
none 
ft 

Fraction of light transmitted by leaf surface 
none 
ft 

Fraction of light reflected by soil 
none 
direct vs diffuse 
2.29.7. Photosynthesis¶
2.29.7.1. Fundamental photosynthetic physiology theory¶
In this section we describe the physiological basis of the photosynthesis model before describing its application to the FATES canopy structure. This description in this section is largely repeated from the Oleson et al. CLM4.5 technical note but included here for comparison with its implementation in FATES. Photosynthesis in C3 plants is based on the model of Farquhar 1980 as modified by Collatz et al. (1991). Photosynthetic assimilation in C4 plants is based on the model of Collatz et al. (1991). In both models, leaf photosynthesis, (mol CO m s) is calculated as the minimum of three potentially limiting fluxes, described below:
The RuBP carboxylase (Rubisco) limited rate of carboxylation (mol CO m s) is determined as
where is the internal leaf CO partial pressure (Pa) and ) is the O partial pressure (Pa). and are the MichaelisMenten constants (Pa) for CO and O. These vary with vegetation temperature (C) according to an Arrhenious function described in Oleson et al. 2013. is the leaf layer photosynthetic capacity ( mol CO m s).
The maximum rate of carboxylation allowed by the capacity to regenerate RuBP (i.e., the lightlimited rate) (mol CO m s) is
To find , the electron transport rate ( mol CO m s), we solve the following quadratic term and take its smaller root,
where is the maximum potential rate of electron transport (mol m s), is the is the light utilized in electron transport by photosystem II (mol m s) and is is curvature parameter. is determined as
where is the absorbed photosynthetically active radiation (Wm:math:^{2}) for either sunlit or shaded leaves ( and ). is converted to photosynthetic photon flux assuming 4.6 mol photons per joule. Parameter values are = 0.7 for C3 and = 0.85 for C4 plants.
The export limited rate of carboxylation for C3 plants and the PEP carboxylase limited rate of carboxylation for C4 plants (also in mol CO m s) is
is the triosephosphate limited rate of photosynthesis, which is equal to . is the initial slope of C4 CO response curve. The MichaelisMenten constants and are modeled as follows,
where = 30.0 and = 30000.0 are values (Pa) at 25 C, and = 2.1 and =1.2 are the relative changes in and respectively, for a 10C change in temperature. The CO compensation point (Pa) is
where the term 0.21 represents the ratio of maximum rates of oxygenation to carboxylation, which is virtually constant with temperature Farquhar, 1980.
2.29.7.2. Resolution of the photosynthesis theory within the FATES canopy structure.¶
The photosynthesis scheme is modified from the CLM4.5 model to give estimates of photosynthesis, respiration and stomatal conductance for a three dimenstional matrix indexed by canopy level (), plant functional type () and leaf layer (). We conduct the photosynthesis calculations at each layer for both sunlit and shaded leaves. Thus, the model also generates estimates of and indexed in the same three dimensional matrix. In this implementation, some properties (stomatal conductance parameters, topofcanopy photosynthetic capacity) vary with plant functional type, and some vary with both functional type and canopy depth (absorbed photosynthetically active radiation, nitrogenbased variation in photosynthetic properties). The remaining drivers of photosynthesis (, , , , temperature, atmospheric CO) remain the same throughout the canopy. The rate of gross photosynthesis ()is the smoothed minimum of the three potentially limiting processes (carboxylation, electron transport, export limitation), but calculated independently for each leaf layer:
For , we use
where now varies with PFT, canopy depth and layer (see below). Internal leaf ( is tracked seperately for each leaf layer. For the light limited rate , we use
where is calculated as above but based on the absorbed photosynthetically active radiation( ) for either sunlit or shaded leaves in Wm. Specifically,
The export limited rate of carboxylation for C3 plants and the PEP carboxylase limited rate of carboxylation for C4 plants (also in mol CO m s) is calculated in a similar fashion,
2.29.7.3. Variation in plant physiology with canopy depth¶
Both and vary with vertical depth in the canopy on account of the welldocumented reduction in canopy nitrogen through the leaf profile, see Bonan et al. 2012 for details). Thus, both and are indexed by by , and according to the nitrogen decay coefficient and the amount of vegetation area shading each leaf layer ,
where and are the topofcanopy photosynthetic rates. is the sum of exposed leaf area index () and the exposed stem area index ()( m m ). Namely,
The vegetation index shading a particular leaf layer in the top canopy layer is equal to
For lower canopy layers, the weighted average vegetation index of the canopy layer above () is added to this withincanopy shading. Thus,
where is calculated as
is the coefficient of nitrogen decay with canopy depth. The value of this parameter is taken from the work of Lloyd et al. 2010 who determined, from 204 vertical profiles of leaf traits, that the decay rate of N through canopies of tropical rainforests was a function of the at the top of the canopy. They obtain the following term to predict ,
where is again in mol CO m s.
2.29.7.4. Water Stress on gas exchange¶
The top of canopy leaf photosynthetic capacity, , is also adjusted for the availability of water to plants as
where the adjusting factor ranges from one when the soil is wet to zero when the soil is dry. It depends on the soil water potential of each soil layer, the root distribution of the plant functional type, and a plantdependent response to soil water stress,
where is a plant wilting factor for layer and is the fraction of roots in layer .The plant wilting factor is
where is the soil water matric potential (mm) and and are the soil water potential (mm) when stomata are fully closed or fully open, respectively. The term in brackets scales the ratio of the effective porosity (after accounting for the ice fraction) relative to the total porosity. = 0 when the temperature of the soil layer ( ) is below some threshold (2:math:^{o}C) or when there is no liquid water in the soil layer (). For more details on the calculation of soil matric potential, see the CLM4.5 technical note.
2.29.7.4.1. Variation of water stress and water uptake within tiles¶
The remaining drivers of the photosynthesis model remain constant (atmospheric CO and O and canopy temperature) throughout the canopy, except for the water stress index . must be indexed by , because plants of differing functional types have the capacity to have varying root depth, and thus access different soil moisture profile and experience differing stress functions. Thus, the water stress function applied to gas exchange calculation is now calculated as
where is the water stress at each soil layer and is the root fraction of each PFT’s root mass in layer . Note that this alteration of the parameter also necessitates recalculation of the vertical water extraction profiles. In the original model, the fraction of extraction from each layer () is the product of a single root distribution, because each patch only has one plant functional type. In FATES, we need to calculate a new weighted patch effective rooting depth profile as the weighted average of the functionaltype level stress functions and their relative contributions to canopy conductance. Thus for each layer , the extraction fraction is summed over all PFTs as
where is the number of soil layers, is the total canopy (see section 9 for details) and is the canopy conductance for plant functional type ,
2.29.7.5. Aggregation of assimilated carbon into cohorts¶
The derivation of photosynthetic rates per leaf layer, as above, give us the estimated rate of assimilation for a unit area of leaf at a given point in the canopy in mol CO m s. To allow the integration of these rates into fluxes per individual tree, or cohort of trees (gCO:math:_2 tree s), they must be multiplied by the amount of leaf area placed in each layer by each cohort. Each cohort is described by a single functional type, and canopy layer flag, so the problem is constrained to integrating these fluxes through the vertical profile ().
We fist make a weighted average of photosynthesis rates from sun (, mol CO m s) and shade leaves ( , mol CO m s) as
The assimilation per leaf layer is then accumulated across all the leaf layers in a given cohort (coh) to give the cohortspecific gross primary productivity (),
The is the exposed leaf area which is present in each leaf layer in m m. (For all the leaf layers that are completely occupied by a cohort, this is the same as the leaf fraction of ). The fluxes are converted from mol into mol and then multiplied by 12 (the molecular weight of carbon) to give units for GPP of KgC cohort s. These are integrated for each timestep to give KgC cohort day
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Maximum carboxylation capacity 
mol CO m s 
ft 

Base Rate of Respiration 
gC gN) 

Temp. Response of stem and root respiration 

CN ratio of leaf matter 
gC/gN 
ft 

CN ratio of root matter 
gC/gN 
ft 

Growth Respiration Fraction 
none 

Water content when stomata close 
Pa 
ft 

Water content above which stomata are open 
Pa 
ft 
2.29.8. Plant respiration¶
Plant respiration per individual (KgC individual s) is the sum of two terms, growth and maintenance respiration and
Maintenance respiration is the sum of the respiration terms from four different plant tissues, leaf, , fine root , coarse root and stem , all also in (KgC individual s) .
To calculate canopy leaf respiration, which varies through we canopy, we first determine the topofcanopy leaf respiration rate (, gC s m) is calculated from a base rate of respiration per unit leaf nitrogen derived from Ryan et al. 1991. The base rate for leaf respiration () is 2.525 gC/gN s,
where is the base rate of metabolism (2.525 x 10 gC/gN s. This base rate is adjusted assuming a Q of 1.5 to scale from the baseline of 20C to the CLM default base rate temperature of 25C. For use in the calculations of net photosynthesis and stomatal conductance, leaf respiration is converted from gC s m, into mol CO m s ().
This topofcanopy flux is scaled to account for variation in through the vertical canopy, in the same manner as the values are scaled using .
Leaf respiration is also adjusted such that it is reduced by drought stress, , and canopy temperature, . For details of the temperature functions affecting leaf respiration see the CLM4 technical note, Section 8, Equations 8.13 and 8.14. The adjusted leaf level fluxes are scaled to individuallevel (gC individual s) in the same fashion as the calculations
The stem and the coarseroot respiration terms are derived using the same base rate of respiration per unit of tissue Nitrogen.
Here, is a temperature relationship based on a value of 1.5, where is the vegetation temperature. We use a base rate of 20 here as, again, this is the baseline temperature used by Ryan et al. 1991. The 10 converts from gC invididual s to KgC invididual s
The tissue N contents for live sapwood are derived from the leaf CN ratios, and for fine roots from the root CN ratio as:
and
where and are the biomass pools of sapwood and live root biomass respectively (KgC individual) and is the fraction of coarse root tissue in the root pool (0.5 for woody plants, 0.0 for grasses and crops). We assume here that stem CN ratio is the same as the leaf C:N ratio, for simplicity. The final maintenance respiration term is derived from the fine root respiration, which accounts for gradients of temperature in the soil profile and thus calculated for each soil layer as follows:
is a function of soil temperature in layer that has the same form as that for stem respiration, but uses vertically resolved soil temperature instead of canopy temperature. In the CLM4.5, only coarse and not fine root respriation varies as a function of soil depth, and we maintain this assumption here, although it may be altered in later versions. The growth respiration, is a fixed fraction of the carbon remaining after maintenance respiration has occurred.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Rate of reduction of N through the canopy 
none 

Base Rate of Respiration 
gC gN) 

Temp. Response of stem and root respiration 

CN ratio of leaf matter 
gC/gN 
ft 

CN ratio of root matter 
gC/gN 
ft 

Growth Respiration Fraction 
none 
2.29.9. Stomatal Conductance¶
2.29.9.1. Fundamental stomatal conductance theory¶
Stomatal conductance is unchanged in concept from the CLM4.5 approach. Leaf stomatal resistance is calculated from the BallBerry conductance model as described by Collatz et al. (1991) and implemented in a global climate model by Sellers et al. 1996. The model relates stomatal conductance (i.e., the inverse of resistance) to net leaf photosynthesis, scaled by the relative humidity at the leaf surface and the CO concentration at the leaf surface. The primary difference between the CLM implementation and that used by Collatz et al. (1991) and Sellers et al. (1996) is that they used net photosynthesis (i.e., leaf photosynthesis minus leaf respiration) instead of gross photosynthesis. As implemented here, stomatal conductance equals the minimum conductance () when gross photosynthesis () is zero. Leaf stomatal resistance is
where is leaf stomatal resistance (s m mol), is a plant functional type dependent parameter equivalent to in the BallBerry model literature. This parameter is also scaled by the water stress index . Similarly, is the slope of the relationship between the assimilation, and humidty dependant term and the stomatal conductance, and so is equivalent to the term in the stomatal literature. is leaf photosynthesis (mol CO m s), is the CO partial pressure at the leaf surface (Pa), is the vapor pressure at the leaf surface (Pa), is the saturation vapor pressure (Pa) inside the leaf at the vegetation temperature conductance (mol m s) when = 0 . Typical values are = 9 for C plants and = 4 for C plants ( Collatz et al. 1991, Collatz, 1992, Sellers et al 1996). Sellers et al. 1996 used = 10000 for C plants and = 40000 for C plants. Here, was chosen to give a maximum stomatal resistance of 20000 s m. These terms are nevertheless plant strategy dependent, and have been found to vary widely with plant type Medlyn et al. 2011.
Resistance is converted from units of s m mol to s m as: 1 s m = R (mol m s), where R is the universal gas constant (J K kmol) and is the atmospheric potential temperature (K).
2.29.9.2. Resolution of stomatal conductance theory in the FATES canopy structure¶
The stomatal conductance is calculated, as with photosynthesis, for each canopy, PFT and leaf layer. The CLM code requires a single canopy conductance estimate to be generated from the multilayer multiPFT array. In previous iterations of the CLM, sun and shadeleaf specific values have been reported and then averaged by their respective leaf areas. In this version, the total canopy condutance , is calculated as the sum of the cohortlevel conductance values.
Cohort conductance is the sum of the inverse of the leaf resistances at each canopy layer ( ) multipled by the area of each cohort.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Slope of BallBerry term 
none 
ft 

Slope of BallBerry term 
none 
ft 
2.29.10. Allocation and Growth¶
Total assimilation carbon enters the ED model each day as a cohortspecific Net Primary Productivity , which is calculated as
This flux of carbon is allocated between the demands of tissue turnover, of carbohydrate storage and of growth (increase in size of one or many plant organs). Priority is explicitly given to maintenance respiration, followed by tissue maintenance and storage, then allocation to live biomass and then to the expansion of structural and live biomass pools. All fluxes here are first converted into in KgC individual year and ultimately integrated using a timesteps of 1/365 years for each day.
2.29.10.1. Tissue maintenance demand¶
We calculate a ‘tissue maintenance’ flux. The magnitude of this flux is such that the quantity of biomass in each pool will remain constant, given background turnover rates. For roots, this maintenenace demand is simply
Where is the root turnover rate in y^1. Given that, for deciduous trees, loss of leaves is assumed to happen only one per growing season, the algorithm is dependent on phenological habit (whether or not this PFT is evergreen), thus
Leaf litter resulting from deciduous senescence is handled in the phenology section. The total quantity of maintenance demand (. KgC individual y) is therefore
2.29.10.2. Allocation to storage and turnover¶
The model must now determine whether the NPP input is sufficient to meet the maintenance demand and keep tissue levels constant. To determine this, we introduce the idea of ‘carbon balance’ (KgC individual) where
where is the minimum fraction of the maintenance demand that the plant must meet each timestep, which is indexed by ft and represents a lifehistorystrategy decision concerning whether leaves should remain on in the case of low carbon uptake (a risky strategy) or not be replaced (a conservative strategy). Subsequently, we determine a flux to the storage pool, where the flux into the pool, as a fraction of , is proportional to the discrepancy between the target pool size and the actual pool size where
The allocation to storage is a fourth power function of to mimic the qualitative behaviour found for carbon allocation in arabidopsis by Smith et al. 2007.
If the carbon remaining after the storage and minimum turnover fluxes have been met, the next priority is the remaining flux to leaves . If the quantity of carbon left is insufficient to supply this amount of carbon, then the store of alive carbon is depleted (to represent those leaves that have fallen off and not been replaced)
correspondingly, the carbon left over for growth (: (KgC individual year) is therefore
to allocate the remaining carbon (if there is any), we first ascertain whether the live biomass pool is at its target, or whether is has been depleted by previous low carbon timesteps. Thus
where the target leaf biomass ((Kg C individual)) is the allometric relationship between dbh and leaf biomass, ameliorated by the leaf trimming fraction (see ‘control of leaf area’ below)
is the wood density, in g cm.
2.29.10.3. Allocation to Seeds¶
The fraction remaining for growth (expansion of live and structural tissues) is 1 minus that allocated to seeds.
Allocation to seeds only occurs if the alive biomass is not below its target, and then is a predefined fixed fraction of the carbon remaining for growth. Allocation to clonal reproduction (primarily for grasses) occurs when is achieved.
the total amount allocated to seed production ( in KgC individual y) is thus
2.29.10.4. Allocation to growing pools¶
The carbon is then partitioned into carbon available to grow the and pools. A fraction is available to live biomass pools, and a fraction is available to structural pools.
If the alive biomass is lower than its ideal target, all of the available carbon is directed into that pool. Thus:
In this case, the division of carbon between the live and structural pools is derived as the inverse of the sum of the rates of change in live biomass with respect to structural:
To calculate all these differentials, we first start with , where
The rates of change of dbh with respect to leaf and structural biomass are the differentials of the allometric equations linking these terms to each other. Hence,
and where
If then we must also account for allocation for growing taller as:
where
Once we have the , we calculate as
and the sapwood differential as
where
In all of the above terms, height in in m, is in cm, and all biomass pools are in KgCm. The allometric terms for the growth trajectory are all taken from the ED1.0 model, but could in theory be altered to accomodate alternative allometric relationships. Critically, the nonlinear relationships between live and structural biomass pools are maintained in this algorithm, which diverges from the methodology currently deployed in the CLM4.5.
2.29.10.5. Integration of allocated fluxes¶
All of the flux calculations generate differential of the biomass state variables against time (in years). To integrate these differential rates into changes in the state variables, we use a simple simple forward Euler integration. Other methods exist (e.g. ODEINT solvers, Runge Kutta methods etc.), but they are more prone to errors that become difficult to diagnose, and the typically slow rates of change of carbon pools mean that these are less important than they might be in strongly nonlinear systems (soil drainage, energy balance, etc.)
In this case, is set to be one day ( years).
Parameter Symbol 
Parameter Name 
Units 
indexed by 

S 
Target stored biomass as fraction of 
none 
ft 
f 
Minimum fraction of turnover that must be met 
none 
ft 
R 
Fraction allocated to seeds 
none 
ft 
C 
Fraction allocated to clonal reproduction 
none 
ft 
Diameter at which maximum height is achieved 
m 
ft 

P 
Does this cohort have an evergreen phenological habit? 
1=yes, 0=no 
ft 
2.29.11. Control of Leaf Area Index¶
The leaf area (m:math:^{2}) of each cohort is calculated from leaf biomass (kgC individual) and specific leaf area (SLA, m kg C)
For a given tree allometry, leaf biomass is determined from basal area using the function used by Moorcroft et al. 2001 where is wood density in g cm.
However, using this model, where leaf area and crown area are both functions of diameter, the leaf area index of each tree in a closed canopy forest is always the same (where = , irrespective of the growth conditions. To allow greater plasticity in tree canopy structure, and for tree leaf area index to adapt to prevailing conditions, we implemented a methodology for removing those leaves in the canopy that exist in negative carbon balance. That is, their total annual assimilation rate is insufficient to pay for the turnover and maintenance costs associated with their supportive root and stem tissue, plus the costs of growing the leaf. The tissue turnover maintenance cost (KgC m of leaf is the total maintenance demand divided by the leaf area:
The net uptake for each leaf layer in (KgC m year) is
where is the GPP of each layer of leaves in each tree (KgC m year), is the rate of leaf dark respiration (also KgC m year). We use an iterative scheme to define the cohort specific canopy trimming fraction , on an annual timestep, where
If the annual maintenance cost of the bottom layer of leaves (KgC m2 year1) is less than then the canopy is trimmed by an increment (0.01), which is applied until the end of next calander year. Because this is an optimality model, there is an issue of the timescale over which net assimilation is evaluated, the timescale of response, and the plasticity of plants to respond to these pressures. These properties should be investigated further in future efforts.
We impose an arbitrary minimum value on the scope of canopy trimming of (0.5). If plants are able simply to drop all of their canopy in times of stress, with no consequences, then tree mortality from carbon starvation is much less likely to occur because of the greatly reduced maintenance and turnover requirements.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Fraction by which leaf mass is reduced next year 
none 

Minimum fraction to which leaf mass can be reduced 
2.29.12. Phenology¶
2.29.12.1. Cold Deciduous Phenology¶
2.29.12.1.1. Cold Leafout timing¶
The phenology model of Botta et al. 2000 is used in FATES to determine the leafon timing. The Botta et al. model was verified against satellite data and is one of the only globally verified and published models of leafout phenology. This model differs from the phenology model in the CLM4.5. The model simulates leafon date as a function of the number of growing degree days (GDD), defined by the sum of mean daily temperatures ( C) above a given threshold (0 C).
Budburst occurs when exceeds a threshold (). The threshold is modulated by the number of chilling days experienced (NCD) where the mean daily temperature falls below a threshold determined by Botta et al. 2000<botta2000> as 5C. A greater number of chilling days means that fewer growing degree days are required before budburst:
where a = 68, b= 638 and c=0.01 Botta et al. 2000<botta2000>. In the Northern Hemisphere, counting of degree days begins on 1st January, and of chilling days on 1st November. The calendar opposite of these dates is used for points in the Southern Hemisphere.
If the growing degree days exceed the critical threshold, leafon is triggered by a change in the gridcell phenology status flag where ‘2’ indicates that leaves should come on and ‘1’ indicates that they should fall.
2.29.12.1.2. Cold Leafoff timing¶
The leafoff model is taken from the Sheffield Dynamic Vegetation Model (SDGVM) and is similar to that for LPJ Sitch et al. 2003 and IBIS Foley et al. 1996 models. The average daily temperatures of the previous 10 day period are stored. Senescence is triggered when the number of days with an average temperature below 7.5 () rises above a threshold values , set at 5 days.
2.29.12.1.3. Global implementation modifications¶
Because of the global implementation of the colddeciduous phenology scheme, adjustments must be made to account for the possibility of colddeciduous plants experiencing situations where no chilling period triggering leafoff ever happens. If left unaccounted for, these leaves will last indefinitely, resulting in highly unrealistic behaviour. Therefore, we implement two additional rules. Firstly, if the number of days since the last senescence event was triggered is larger than 364, then leafoff is triggered on that day. Secondly, if no chilling days have occured during the winter accumulation period, then leafon is not triggered. This means that in effect, where there are no cold periods, leaves will fall off and not come back on, meaning that colddeciduous plants can only grow in places where there is a cold season.
Further to this rule, we introduce a ‘buffer’ time periods after leafon of 30 days, so that coldsnap periods in the spring cannot trigger a leaf senescence. The 30 day limit is an arbitrary limit. In addition, we constrain growing degree day accumulation to the second half of the year (Jult onwards in the Northern hemisphere, or JanJune in the Southern) and only allow GDD accumulation while the leaves are off.
2.29.12.2. Droughtdeciduous Phenology: TBD¶
In the current version of the model, a drought deciduous algorithm exists, but is not yet operational, due to issue detected in the existing CN and soil moisture modules, which also affect the behaviour of the native ED drought deciduous model. This is a priority to address before the science tag is released.
2.29.12.3. Carbon Dynamics of deciduous plants¶
In the present version, leaf expansion and senescence happen over the course of a single day. This is clearly not an empirically robust representation of leaf behaviour, whereby leaf expansion occurs over a period of 1014 days, and senescence over a similar period. This will be incorporated in later versions. When the cold or drought phenological status of the gridcell status changes () from ‘2’ to ‘1’, and the leaves are still on ( =2 ), the leaf biomass at this timestep is ’remembered’ by the model state variable . This provides a ‘target’ biomass for leaf onset at the beginning of the next growing season (it is a target, since depletion of stored carbon in the off season may render achieving the target impossible).
Leaf carbon is then added to the leaf litter flux (KgC individual)
The alive biomass is depleted by the quantity of leaf mass lost, and the leaf biomass is set to zero
Finally, the status is set to 1, indicating that the leaves have fallen off.
For bud burst, or leafon, the same occurs in reverse. If the leaves are off (=1) and the phenological status triggers budburst (=2) then the leaf mass is set the maximum of the leaf memory and the available store
this amount of carbon is removed from the store
and the new leaf biomass is added to the alive pool
Lastly, the leaf memory variable is set to zero and the phenological status of the cohort back to ‘2’. No parameters are currently required for this carbon accounting scheme.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Threshold of cold days for senescence 
none 

Threshold for counting growing degree days 
C 
2.29.13. Seed Dynamics and Recruitment¶
The production of seeds and their subsequent germination is a process that must be captured explicitly or implicitly in vegetation models. FATES contains a seed bank model designed to allow the dynamics of seed production and germination to be simulated independently. In the ED1.0 model, seed recruitment occurs in the same timestep as allocation to seeds, which prohibits the survival of a viable seed bank through a period of disturbance or low productivity (winter, drought). In FATES, a plant functional type specific seed bank is tracked in each patch ( KgC m), whose rate of change (KgC m y) is the balance of inputs, germination and decay:
where , and are the production, germination and decay (or onset of inviability) of seeds, all in KgC m year.
Seeds are assumed to be distributed evenly across the site (in this version of the model), so the total input to the seed pool is therefore the sum of all of the reproductive output of all the cohorts in each patch of the correct PFT type.
Seed decay is the sum of all the processes that reduce the number of seeds, taken from Lischke et al. 2006. Firstly, the rate at which seeds become inviable is described as a constant rate (y:math:^{1}) which is set to 0.51, the mean of the parameters used by Lischke et al. 2006.
The seed germination flux is also prescribed as a fraction of the existing pool (), but with a cap on maximum germination rate , to prevent excessive dominance of one plant functional type over the seed pool.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Maximum seed mass 
kgC m 

Proportional germination rate 

Maximum germination rate 
KgC m y 

Decay rate of viable seeds 
none 
FT 

Fraction of devoted to reproduction 
none 
FT 
2.29.14. Litter Production and Fragmentation¶
and are indexed by plant functional type (). The rational for indexing leaf and fine root by PFT is that leaf and fine root matter typically vary in their carbon:nitrogen ratio, whereas woody pools typically do not.
Rates of change of litter, all in kGC m year, are calculated as
2.29.14.1. Litter Inputs¶
Inputs into the litter pools come from tissue turnover, mortality of canopy trees, mortality of understorey trees, mortality of seeds, and leaf senescence of deciduous plants.
where is the leaf turnover rate for evergreen trees and is the leaf loss from phenology in that timestep (KgC . is the total mortality flux in that timestep (in individuals). For fine root input. is the number of cohorts of functional type ‘’ in the current patch.
where is the root turnover rate. For coarse woody debris input ( , we first calculate the sum of the mortality and turnover ) fluxes, then separate these into size classes and above/below ground fractions using the fixed fractions assigned to each ( and )
2.29.14.2. Litter Outputs¶
The fragmenting litter pool is available for burning but not for respiration or decomposition. Fragmentation rates are calculated according to a maximum fragmentation rate ( or ) which is ameliorated by a temperature and water dependent scalar . The form of the temperature scalar is taken from the existing CLM4.5BGC decomposition cascade calculations). The water scaler is equal to the water limitation on photosynthesis (since the CLM4.5BGC water scaler pertains to the water potential of individual soil layers, which it is difficult to meaningfully average, given the nonlinearities in the impact of soil moisture). The scaler code is modular, and new functions may be implemented trivially. Rate constants for the decay of the litter pools are extremely uncertain in literature, as few studies either separate litter into size classes, nor examine its decomposition under nonlimiting moisture and temperature conditions. Thus, these parameters should be considered as part of sensitivity analyses of the model outputs.
2.29.14.3. Flux into decompsition cascade¶
Upon fragmentation and release from the litter pool, carbon is transferred into the labile, lignin and cellulose decomposition pools. These pools are vertically resolved in the biogeochemistry model. The movement of carbon into each vertical layer is obviously different for above and belowground fragmenting pools. For each layer and chemical litter type , we derive a flux from ED into the decomposition cascade as (kGC m s)
where is the time conversion factor from years to seconds, , and are the fractions of labile, cellulose and lignin in leaf litter, and , and are their counterparts for root matter. Similarly, , and are the fractions of leaf, coarse root and fine root matter that are passed into each vertical soil layer , derived from the CLM(BGC) model.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Maximum fragmentation rate of CWD 
y 

Maximum fragmentation rate of leaf litter 
y 

Maximum fragmentation rate of fine root litter 
y 

Fraction of leaf mass in labile carbon pool 
none 

Fraction of leaf mass in cellulose carbon pool 
none 

Fraction of leaf mass in lignin carbon pool 
none 

Fraction of root mass in labile carbon pool 
none 

Fraction of root mass in cellulose carbon pool 
none 

Fraction of root mass in lignin carbon pool 
none 

Fraction of leaf matter directed to soil layer z 
none 
soil layer 

Fraction of coarse root matter directed to soil layer z 
none 
soil layer 

Fraction of fine root matter directed to soil layer z 
none 
soil layer 
2.29.15. Plant Mortality¶
Total plant mortality per cohort , (fraction year) is simulated as the sum of four additive terms,
where is the background mortality that is unaccounted by any of the other mortality rates and is fixed at 0.014. is the carbon starvation derived mortality, which is a function of the nonstructural carbon storage term and the PFTspecific ‘target’ carbon storage, , as follows:
where is the stress mortality parameter, or the fraction of trees in a landscape that die when the mean condition of a given cohort triggers mortality. This parameter is needed to scale from individuallevel mortality simulation to gridcell average conditions.
Mechanistic simulation of hydraulic failure is not undertaken on account of it’s mechanistic complexity (see McDowell et al. 2013 for details). Instead, we use a proxy for hydraulic failure induced mortality () that uses a water potential threshold beyond mortality is triggered, such that the tolerance of low water potentials is a function of plant functional type (as expressed via the parameter). For each day that the aggregate water potential falls below a threshold value, a set fraction of the trees are killed. The aggregation of soil moisture potential across the root zone is expressed using the function. We thus determine plant mortality caused by extremely low water potentials as
The threshold value of 10 represents a state where the average soil moisture potential is within 10 of the wilting point (a PFT specific parameter ).
is the fireinduced mortality, as described in the fire modelling section.
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Stress Mortality Scaler 
none 

Target carbon storage fraction 
none 
ft 
2.29.16. Fire (SPITFIRE)¶
The influence of fire on vegetation is estimated using the SPITFIRE model, which has been modified for use in ED following it’s original implementation in the LPJSPITFIRE model (Thonicke et al. 2010, Pfeiffer et al. 2013). This model as described is substantially different from the existing CLM4.5 fire model Li et al. 2012, however, further developments are intended to increase the merging of SPITFIRE’s natural vegetation fire scheme with the fire suppression, forestclearing and peat fire estimations in the existing model. The coupling to the ED model allows fires to interact with vegetation in a sizestructured manner, so small fires can burn only understorey vegetation. Also, the patch structure and representation of succession in the ED model allows the model to track the impacts of fire on different forest stands, therefore removing the problem of areaaveraging implicit in areabased DGVMs. The SPITFIRE approach has also been coupled to the LPJGUESS individualbased model (Forrest et al. in prep) and so this is not the only implementation of this type of scheme in existence.
The SPITFIRE model operates at a daily timestep and at the patch level, meaning that different litter pools and vegetation charecteristics of open and closed forests can be represented effectively (we omit the patch subscript throughout for simplicity).
2.29.16.1. Properties of fuel load¶
Many fire processes are impacted by the properties of the litter pool in the SPITFIRE model. There are one live (live grasses) and five dead fuel categories (dead leaf litter and four pools of coarse woody debris). Coarse woody debris is classified into 1h, 10h, 100h, and 1000h fuels, defined by the order of magnitude of time required for fuel to lose (or gain) 63% of the difference between its current moisture content and the equilibrium moisture content under defined atmospheric conditions. Thonicke et al. 2010. For the purposes of describing the behaviour of fire, we introduce a new index ‘fuel class’ fc, the values of which correspond to each of the six possible fuel categories as follows.
fc index 
Fuel type 
Drying Time 

1 
dead grass 
n/a 
2 
twigs 
1h fuels 
3 
small branches 
10h fuel 
4 
large branches 
100h fuel 
5 
stems and trunks 
1000h fuel 
6 
live grasses 
n/a 
2.29.16.2. Nesterov Index¶
Dead fuel moisture (), and several other properties of fire behaviour, are a function of the ‘Nesterov Index’ () which is an accumulation over time of a function of temperature and humidity (Eqn 5, Thonicke et al. 2010),
where is the daily mean temperature in C and is the dew point calculated as .
where is the relative humidity (%).
On days when the total precipitation exceeds 3.0mm, the Nesterov index accumulator is reset back to zero.
2.29.16.3. Fuel properties¶
Total fuel load for a given patch is the sum of the above ground coarse woody debris and the leaf litter, plus the alive grass leaf biomass multiplied by the nonmineral fraction (1).
Many of the model behaviours are affected by the patchlevel weighted average properties of the fuel load. Typically, these are calculated in the absence of 1000h fuels because these do not contribute greatly to fire spread properties.
2.29.16.3.1. Dead Fuel Moisture Content¶
Dead fuel moisture is calculated as
where is a parameter defining how fuel moisture content varies between the first four dead fuel classes.
2.29.16.3.2. Live grass moisture Content¶
The live grass fractional moisture content() is a function of the soil moisture content. (Equation B2 in Thonicke et al. 2010)
where is the fractional moisture content of the top 30cm of soil.
2.29.16.3.3. Patch Fuel Moisture¶
The total patch fuel moisture is based on the weighted average of the different moisture contents associated with each of the different live grass and dead fuel types available (except 1000h fuels).
2.29.16.3.4. Effective Fuel Moisture Content¶
Effective Fuel Moisture Content is used for calculations of fuel consumed, and is a function of the ratio of dead fuel moisture content and the moisture of extinction factor,
where the is a function of surfacearea to volume ratio.
2.29.16.3.5. Patch Fuel Moisture of Extinction¶
The patch ‘moisture of extinction’ factor () is the weighted average of the of the different fuel classes
2.29.16.3.6. Patch Fuel Bulk Density¶
The patch fuel bulk density is the weighted average of the bulk density of the different fuel classes (except 1000h fuels).
where is the bulk density of each fuel size class (kG m)
2.29.16.3.7. Patch Fuel Surface Area to Volume¶
The patch surface area to volume ratio () is the weighted average of the surface area to volume ratios () of the different fuel classes (except 1000h fuels).
2.29.16.4. Forward rate of spread¶
For each patch and each day, we calculate the rate of forward spread of the fire ros (nominally in the direction of the wind).
is the effective heating number (). is the heat of preignition (). is the propagating flux calculated as (see Thonicke et al. 2010 Appendix A).
is the influence of windspeed on rate of spread.
Where , and are all functions of surfaceareavolume ratio : , , . where is the the windspeed in ms, and where is the particle density (513).
is the reaction intensity, calculated using the following set of expressions (from Thonicke et al. 2010 Appendix A).:
is the residence time of the fire, and is the mineral damping coefficient (=0.174:math:S_e^{0.19} , where is 0.01 and so = 0.41739).
2.29.16.5. Fuel Consumption¶
The fuel consumption (fraction of biomass pools) of each dead biomass pool in the area affected by fire on a given day () is a function of effective fuel moisture and size class fc (Eqn B1, B4 and B5, Thonicke et al. 2010). The fraction of each fuel class that is consumed decreases as its moisture content relative to its moisture of extinction () increases.
and are parameters, the value of which is modulated by both size class and by the effective fuel moisture class , defined by . and are defined for low, mid, and highmoisture conditions, the boundaries of which are also functions of the litter size class following Peterson and Ryan 1986 (page 802). The fuel burned, (Kg m day) iscalculated from for each fuel class:
Where 0.45 converts from carbon to biomass. The total fuel consumption, (Kg m), used to calculate fire intensity, is then given by
There is no contribution from the 1000 hour fuels to the patchlevel used in the fire intensity calculation.
2.29.16.6. Fire Intensity¶
Fire intensity at the front of the burning area (, kW m) is a function of the total fuel consumed () and the rate of spread at the front of the fire, (m min) (Eqn 15 Thonicke et al. 2010)
where is the energy content of fuel (Kj/Kg  the same, 18000 Kj/Kg for all fuel classes). Fire intensity is used to define whether an ignition is successful. If the fire intensity is greater than 50Kw/m then the ignition is successful.
2.29.16.7. Fire Duration¶
Fire duration is a function of the fire danger index with a maximum length of (240 minutes in Thonicke et al. 2010 Eqn 14, derived from Canadian Forest Fire Behaviour Predictions Systems)
2.29.16.8. Fire Danger Index¶
Fire danger index (fdi) is a representation of the effect of meteorological conditions on the likelihood of a fire. It is calculated for each gridcell as a function of the Nesterov Index . is calculated from as
where = 0.00037 following Venevsky et al. 2002.
2.29.16.9. Area Burned¶
Total area burnt is assumed to be in the shape of an ellipse, whose major axis (m) is determined by the forward and backward rates of spread ( and respectively).
is a function of and windspeed (Eqn 10 Thonicke et al. 2010)
The minor axis to major axis ratio of the ellipse is determined by the windspeed. If the windspeed () is less than 16.67 ms then . Otherwise (Eqn 12 and 13, Thonicke et al. 2010)
and are the fractions of the patch surface covered by grass and trees respectively.
The total area burned ( in m) is therefore (Eqn 11, Thonicke et al. 2010)
where is the number of fires.
2.29.16.10. Crown Damage¶
is the fraction of the crown which is consumed by the fire. This is calculated from scorch height , tree height and the crown fraction parameter (Eqn 17 Thonicke et al. 2010):
The scorch height (m) is a function of the fire intensity, following Byram, 1959, and is proportional to a plant functional type specific parameter (Eqn 16 Thonicke et al. 2010):
where is the fraction of the aboveground biomass in each plant functional type.
2.29.16.11. Cambial Damage and Kill¶
The cambial kill is a function of the fuel consumed , the bark thickness , and , the duration of cambial heating (minutes) (Eqn 8, Peterson and Ryan 1986):
Bark thickness is a linear function of tree diameter , defined by PFTspecific parameters and (Eqn 21 Thonicke et al. 2010):
The critical time for cambial kill, (minutes) is given as (Eqn 20 Thonicke et al. 2010):
The mortality rate caused by cambial heating of trees within the area affected by fire is a function of the ratio between and (Eqn 19, Thonicke et al. 2010):
Parameter Symbol 
Parameter Name 
Units 
indexed by 

Intercept of bark thickness function 
mm 
FT 

Slope of bark thickness function 
mm cm 
FT 

Ratio of crown height to total height 
none 
FT 

Fuel moisture parameter 
C 
fc 

Fuel Bulk Density 
kG m 
fc 

Surface area to volume ratio 
cm 
fc 

Intercept of fuel burned 
none 
, moisture class 

Slope of fuel burned 
none 
, moisture class 

Fuel Mineral Fraction 

Maximum Duration of Fire 
Minutes 

Energy content of fuel 
kJ/kG 

Flame height parameter 
FT 