The global distribution of vegetation is broadly determined by climate, and where bioclimatic parameters are favorable for several plant functional types (PFTs), by the competition between them. Most current dynamic global vegetation models (DGVMs) do not, however, explicitly simulate inter-PFT competition and instead determine the existence and fractional coverage of PFTs based on quasi-equilibrium climate–vegetation relationships. When competition is explicitly simulated, versions of Lotka–Volterra (LV) equations developed in the context of interaction between animal species are almost always used. These equations may, however, exhibit unrealistic behavior in some cases and do not, for example, allow the coexistence of different PFTs in equilibrium situations. Coexistence may, however, be obtained by introducing features and mechanisms such as temporal environmental variation and disturbance, among others.
A generalized version of the competition equations is proposed that includes the LV equations as a special case, which successfully models competition for a range of climate and vegetation regimes and for which coexistence is a permissible equilibrium solution in the absence of additional mechanisms. The approach is tested for boreal forest, tropical forest, savanna, and temperate forest locations within the framework of the Canadian Terrestrial Ecosystem Model (CTEM) and successfully simulates the observed successional behavior and the observed near-equilibrium distribution of coexisting PFTs.
The terrestrial biosphere plays an important role in determining the climate through a range of biophysical and biogeochemical processes. Terrestrial vegetation affects the way in which energy and water are exchanged between the land and the atmosphere on time scales of minutes to months, and the climatic implications of such biophysical interactions have been studied by various authors (Charney 1975; Lean and Rowntree 1993; Xue 1997; Douville et al. 2000; Heck et al. 2001). Our understanding of biophysical interactions between land and atmosphere is incorporated in complex soil–vegetation–atmosphere transfer (SVAT) parameterizations that now form an integral part of climate models (Henderson-Sellers et al. 1996).
The ways in which the terrestrial biosphere affects the climate through the biogeochemical pathway are less well known and modeled. The biogeochemical pathway incorporates the land–atmosphere exchange of CO2, other radiatively active trace gases, and aerosols and their precursors. The associated radiative processes affect the global climate at time scales from decades to centuries. Natural (e.g., fires and the competition between vegetation types) and anthropogenic (e.g., land use) processes that change the geographical distribution of vegetation affect climate through both biophysical and biogeochemical pathways. Attempts to better incorporate these ecosystem processes into climate models are currently underway (e.g., Jones et al. 2001; Friedlingstein et al. 2001; Bonan et al. 2003).
The biosphere responds to changes in climate by changing structural attributes [e.g., leaf area index (LAI) and vegetation height] and the distribution of vegetation. The distribution of terrestrial vegetation is affected by changes in climate in a manner determined by the competitive ability of plant species and the nature of the climatic disturbances, and this is the focus of the current study. The dynamic vegetation subcomponents of climate models have evolved from earlier approaches in which vegetation was assumed to be in equilibrium with and determined by climate (e.g., Henderson-Sellers 1993; Bonan et al. 2003; Delire et al. 2003). By their nature, equilibrium vegetation models are unsuitable for use in transient climate change simulations, although aspects of equilibrium vegetation models may be found in some current dynamic global vegetation models (DGVMs).
The vegetation in DGVMs is characterized in terms of plant functional types (PFTs) in which species with similar attributes are treated together. For instance, species that are characterized by similar leaf form (e.g., needleleaf or broadleaf) and physiological behavior (e.g., deciduous or evergreen) are grouped together as PFTs (e.g., needleleaf evergreen, broadleaf evergreen, broadleaf deciduous, etc.). The geographical distribution of plant species and hence PFTs is determined by competition between them, mediated by climate, soil conditions, and so on. Competition is not explicitly simulated in most current DGVMs, however, and among models that do simulate competition the most common approach is based on the Lotka–Volterra (LV) equations that originated in the context of predator–prey interaction (Cox 2001; Brentall et al. 2005). In some cases at least, the LV equations lead to unrealistic outcomes when used for competition between PFTs since, for instance, their equilibrium solutions prohibit the coexistence of PFTs in the absence of appreciable disturbance.
Generalized competition–colonization (cc) equations that overcome some of the limitations of the existing approaches, and which include the LV equations as a special case, are proposed and are incorporated and tested in the Canadian Terrestrial Ecosystem Model (CTEM). CTEM is a dynamic global vegetation model designed for inclusion in the Canadian Centre for Climate Modeling and Analysis (CCCma) coupled general circulation model (CGCM). CTEM is able to grow vegetation from bare ground and incorporates the processes of photosynthesis, autotrophic and heterotrophic respiration, allocation, phenology, turnover, mortality, fire, and land-use change (Arora 2003; Arora and Boer 2003; Arora and Boer 2005b). Section 2 provides a brief overview of approaches that are currently used to model competition in DGVMs. The modified Lotka–Volterra version of the competition–colonization equations is developed in section 3. The structure of the CTEM and the climate and soil data used to drive the proposed competition approach and CTEM are discussed in sections 4 and 5, respectively. The proposed competition–colonization model is tested for different climates and biomes within the framework of CTEM in section 6 and compared with the special case of the LV model. The successional behavior and equilibrium values of simulated fractional coverages of PFTs are tested against available observations for boreal forest, tropical forest, savanna, and temperate forest biomes to illustrate the overall performance of CTEM and of the competition–colonization model. Summary and discussion are presented in section 7.
2. Competition in DGVMs
Both explicit and implicit forms of competition amongst plants are seen in an array of models and approaches, but the treatment differs with the space scale involved. At the small forest patch scale (100–1000 m2), for instance, gap models explicitly simulate competition by tracking tree height and light competition (Bugmann 2001). However, the use of such models at the comparatively large spatial scales (>10 000 km2) characteristic of climate models is precluded because of the large computational expense involved. Ecosystem models operating on small scales but which do not model forest patch dynamics have used the LV equations (Lotka 1925; Volterra 1926) to simulate competition (e.g., Fernandez-Illescas and Rodriguez-Iturbe 2003; Fernandez-Illescas and Rodriguez-Iturbe 2004). Variants of these equations have also been used in numerous studies addressing biodiversity at the landscape scale (e.g., Tilman 1994; Pacala and Rees 1998; Bampfylde et al. 2005). At large spatial scales (>10 000 km2), DGVMs have generally used quasi-equilibrium climate–vegetation relationships in conjunction with simple functional relationships that determine PFT existence and their fractional coverage (e.g., Kucharik et al. 2000; Potter 2004). More recently, some DGVMs have applied the LV equations to simulate competition at climate model scales (Cox 2001; Brentall et al. 2005).
Current DGVMs have evolved from predecessors that assume climate–vegetation equilibrium. In equilibrium vegetation models, such as those of the BIOME family (Prentice et al. 1992; Haxeltine et al. 1996; Haxeltine and Prentice 1996), the geographical distribution of PFTs is a function of bioclimatic parameters that reflect the long-term mean climate. The work of Woodward (Woodward 1987) and Box (Box 1981), who relate vegetation distribution to environmental constraints, forms the basis of this approach. The BIOME1 model (Prentice et al. 1992), for instance, generates the distribution of seven tree and seven nontree PFTs as functions of the mean temperature of the coldest and warmest months, growing degree-days above 0° and 5°C, and a moisture availability index.
Equilibrium vegetation models have been used to predict vegetation distribution and the resulting changes in carbon storage under future warmer climate conditions (Prentice and Fung 1990; Smith et al. 1992; Solomon et al. 1993), but this approach has been criticized (e.g., Solomon and Kirilenko 1997; Woodward and Beerling 1997) because new PFTs replace the existing PFTs and occupy and dominate regions of favorable climate in step with climate change. Brubaker (Brubaker 1986) points out that shifts in the composition of tree stands depend on the disappearance of existing PFTs so the change in PFTs will lag many years behind the environmental change that ultimately caused it. Even if climatic conditions are unfavorable and the existing PFTs have died out, the invading PFTs require time to mature and become fully established at the newly available sites (Shugart 1984). This time frame is of the order of 50–150 yr, not including the time required for the dispersion of seeds. Woodward and Lomas (Woodward and Lomas 2004a) illustrate the importance of competition in determining the current vegetation distribution by pointing out that in absence of competition the boundaries of species would correspond to the boundaries of their climatic tolerances. The continued long-term productivity and survival of species beyond their observed distributions, such as in botanic gardens, is a testament to this reality. Thus where the climatic tolerance boundaries of PFTs overlap, competition determines the observed distribution.
Some DGVMs introduce time dependence into equilibrium vegetation models by using moving average values of the bioclimatic constraints that determine the existence of PFTs. For example, the Carnegie–Ames–Stanford Approach (CASA) model (Potter and Klooster 1999; Potter 2004) uses 3–5-yr moving average values of its bioclimatic variables, and the Integrated Biosphere Simulator (IBIS) (Kucharik et al. 2000) updates the values of its bioclimatic parameters in an e-folding sense using a time scale of 30 yr. The CASA model determines fractional coverage of trees on the basis of a 5-yr moving average value of soil moisture in the deep layers on the basis of the Walter hypothesis (Walter 1971; Walker and Noy-Meir 1982) that proposes preferential soil water use from shallow and deep soil layers by grasses and tree PFTs, respectively. The IBIS model assumes that the fractional coverage of woody PFTs in a grid cell is a function of woody biomass while LAI is used to specify the coverage of herbaceous PFTs (model code version 2.5). The Sheffield DGVM estimates fractional coverage of C3 and C4 grasses in proportion to their net primary productivities (Woodward and Lomas 2004a; Woodward and Lomas 2004b). The Lund–Potsdam–Jena (LPJ) model (Sitch et al. 2003) uses 20-yr moving average values of its bioclimatic parameters to determine PFT establishment and mortality. Competition is simulated in a simple manner and primarily realized through the capture of the bare fraction in a grid cell (J. Kaplan 2003, personal communication; model code available at http://www.pik-potsdam.de/lpj/). The relationships for determining the fractional coverage of PFTs are implicit rather than explicit in this class of models.
Competition is explicitly simulated on the basis of LV equations in the large-scale Top-down Representation of Interactive Foliage and Flora Including Dynamics (TRIFFID) model (Cox 2001) that serves as the dynamic vegetation component of the Hadley Centre's climate model.
3. Generalized competition–colonization equations
A general approach to competition–colonization would express the fractional coverages of PFTs in terms of the number of individuals plants n(a) in each crown area class (a), where a denotes the area of the plant crown projected onto the surface. The number of such plants in a grid square of area A is n̂ = ∫A0 n(a) da, the area covered is â = ∫A0 n(a) a da, and the fractional coverage is f = (1/A)∫n(a) a da. The evolution of the fractional coverage is
where dn(a)/dt has contributions from all of the different processes that act to change the number density of the PFT. These include the evolution of the number density distribution through self-interaction (e.g., self-thinning), the colonization of bare ground, the invasion of areas occupied by subdominant PFTs, loss through invasion by dominant PFTs, loss through mortality, and evolution to a higher crown area class due to growth. Of these, self-interaction processes, such as self-thinning (due to density-dependent mortality) in which fewer individuals with larger crown areas replace a larger number of individuals with smaller crown areas, will change the number density distribution but not the fractional coverage and so do not contribute to (1). Processes that are proportional to the number density, as is assumed to be the case for mortality, for instance, are easily treated in terms of (1). Thus, for dn(a)/dt|mort ∝ n,
where m is a PFT-dependent mortality rate. The remaining components of (1) that involve competition and colonization involve interactions with other PFTs and are correspondingly more complex.
We consider the fractional coverage equations for N PFTs ( fα) plus the bare ground ( fN+1 = fB) where ΣN+1α=1fα = 1. We reexpress (1) as
where competition–colonization (cc) is symbolically expressed by the function g involving interaction between PFT α and other PFTs β and . The increase in fractional coverage of a dominant PFT due to the establishment of new plants will reduce the fraction of bare ground and/or the coverage of other subdominant PFTs through invasion. The conceptual form of (3) is similar to that of Amarasekare et al. (Amarasekare et al. 2004).
In the absence of a full treatment of (1) we adopt a phenomenological approach based on (3). The basic assumption is that the PFTs may be ordered in terms of “dominance” such that PFT α invades PFT β and the bare ground, where α < β, and such that PFT mortality generates bare ground. For example, woody tree PFTs are dominant over grass PFTs since trees are capable of invading grasses because of their ability to shade them (Siemann and Rogers 2003). Within these two broad functional groups the dominance is calculated on the basis of colonization rates, such that the PFT with the highest colonization rate is considered to be the most dominant. The form adopted for (3) is
and for the bare fraction
where b is an exponent. Here PFT α, with fractional coverage fα, invades subdominant PFTs with fractional coverages fβ at a rate given by
The competition–colonization process depends on a PFT-dependent “invasion rate” cα,β of PFT α into the subdominant PFTs β and the product of their fractional coverages in the form fbαfβ (discussed further below).
This process is alternatively expressed as the product of the unimpeded colonization rate cα of PFT α into bare ground and a measure of the relative efficiency ɛαβ of colonization into other PFTs, ranging between 0 and 1. Thus all PFTs invade bare ground with ɛαB = 1, and woody tree PFTs invade grasses with little impediment (ɛ ≈ 1). However, woody tree PFTs invade other tree PFTs at a lower rate than they invade bare ground or grasses, and this is reflected by values of ɛ < 1. Putting ɛαβ = 0 for β ≤ α (i.e., subdominant PFTs do not invade themselves or dominant PFTs), and ɛαB = 1 for the invasion of the bare ground, allows (4) to be written somewhat more compactly as
The choice of the functional form for competition–colonization is motivated by a few simple considerations. In Equations (4)–(7) b is considered to be an empirical parameter. The usual form of the LV equations for competition among plants (Tilman 1994; Cox 2001; Fernandez-Illescas and Rodriguez-Iturbe 2003; Fernandez-Illescas and Rodriguez-Iturbe 2004; Bampfylde et al. 2005) is recovered from (4) and (7) when b = ɛ = 1. In that case, dominant PFTs invade all subdominant PFTs at the same rate as they colonize bare ground, cαβ = cα, and in proportion to the product of their coverages so that (6) becomes cαfαfβ. If the fractional coverages fα, fβ are taken to indicate the probability of finding a particular PFT in some region of the grid square, then the probability of independently finding both is the product fαfβ. Interaction occurs in these common regions in what might be termed the “random interaction” case.
The choice of b = 0, by contrast, implies that the dominant PFT has full access to all subdominant PFTs and invades them in proportion to their coverage in what might be termed the “full interaction” case. This case may be thought of as corresponding to the general availability of the seeds of the dominant PFT that may germinate and invade the coverage of the subdominant PFT provided the climate is favorable. Other values of b give varying levels of “access” of the dominant PFT to the coverage of a subdominant PFT; b may be regarded as a empirical parameter chosen to improve results in (4) and (5) or to satisfy some other constraint. In particular, coexistence of PFTs is not possible in an equilibrium situation for b = ɛ = 1, and this is an argument for choosing b < 1.
3.1. Equilibrium solution for the LV equations
For the usual form of the LV equations with b = ɛ = 1 and two PFTs with fractional coverages f1 and f2 together with the bare fraction fB the competition–colonization equations are
where dominant PFT 1 invades PFT 2 and the bare fraction, PFT 2 invades the bare fraction, and c1, c2, and m1, m2 are the colonization and mortality rates, respectively. The equilibrium solutions for the fractional coverages of these two PFTs are f1 = max[(c1 − m1)/c1, 0] and f2 = max[(c2 − c2f1 − c1f1 − m2)/c2, 0], respectively. In this version of the equations and in the absence of large frequent disturbances or strong interannual variability in climate, the dominant PFT is always successful in excluding the subdominant PFTs as is seen by rearranging the expression for the equilibrium value of f2 as
As long as (c1 − m1) > (c2 − m2) the equilibrium solution for f2 will always be zero. This is also the case for numerical simulations, as seen in section 6. This lack of coexistence differs from the application of the LV equations in the predator–prey case where stable coexisting nonzero equilibrium populations may be achieved. In that case, the predator depends on the prey for its nourishment so that the predator population suffers as the prey population declines. This is not the case for plants, since the dominant PFT does not depend on subdominant PFTs for its survival and is eventually successful in excluding subdominant PFTs.
Is the eventual exclusion of subdominant PFTs by the dominant PFT a realistic state under equilibrium conditions? The humid Tropics experience year-round favorable growth conditions and are therefore in near equilibrium. They are also well known for their biodiversity. If dominant species were generally successful in driving other plant species out of competition, tropical ecosystems would be expected to be less biologically diverse than they are. Species are apparently well able to coexist in this competitive environment and coexistence is widely observed in nature. We distinguish equilibrium coexistence in a homogenous environment from that which arises in nonequilibrium or heterogeneous environments. Under nonequilibrium conditions transient coexistence may be achieved due to a competition–colonization trade-off associated with climate variability and/or disturbance (Kisdi and Geritz 2003). Temporal environmental variability aids coexistence by providing favorable temporal niches to different plants types at different times. Disturbance aids coexistence by allowing subdominant plant types to temporarily invade gaps that are otherwise occupied by dominants. Coexistence is also realized when dominance rankings vary spatially due to heterogeneity in biotic (e.g., soil characteristics and nutrient availability) and abiotic (or climatic) environments. For instance, in the Canadian boreal region both broadleaf deciduous (aspen) and needleleaf evergreen trees (black spruce and jack pine) exist together, but they preferentially occupy different soil niches. Aspen occurs on fine textured upland soils, black spruce commonly occurs on poorly drained organic lowland soils, and jack pine occurs on excessively drained sandy upland soils (Gower et al. 1997). Here, climate plays the dominant role in the existence of these species, but the second-order control on their spatial distribution is exerted by local soil and topographic conditions. Coexistence under these conditions is aided by temporal and spatial heterogeneity of abiotic and biotic environments, respectively, and may not hold under more homogenous conditions.
Since the LV form of competition–colonization equations (with b = ɛ = 1) does not yield coexistence, under equilibrium conditions, additional mechanisms including temporal environmental variability (Fernandez-Illescas and Rodriguez-Iturbe 2003), random fruiting events (Bampfylde et al. 2005), disturbance (Bond et al. 2003), successional niche dynamics (Pacala and Rees 1998), and spatial heterogeneity in biotic and abiotic environments are invoked to explain coexisting species/PFTs. The application of the LV model for plants without these additional mechanisms, for instance, leads to eventual exclusion or submission of subdominant PFTs as is seen in the TRIFFID model (Cox 2001) used in coupled carbon–climate simulations (Cox et al. 2004; Crucifix et al. 2005). In these simulations, excess simulated tree fractions yield forests in locations that are currently characterized as savannas presumably due to the eventual exclusion or submission of grasses by broadleaf trees.
The problem of excess coverage by trees (or generally by the dominant PFT) also arises in the competition module of the LPJ model (Sitch et al. 2003) in which PFTs primarily compete for the bare fraction and invasion is not represented explicitly. Krinner et al. (Krinner et al. 2005) and Lunt et al. (Lunt et al. 2005) use the LPJ competition submodule in the Organizing Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) vegetation model, Levis et al. (Levis et al. 2004) use the module in the National Center for Atmospheric Research's (NCAR) Community Climate System Model 2, and Beer et al. (Beer et al. 2005) use the module for understanding vegetation dynamics in the Eurasian region. All express concern that in the LPJ model the simulated coverage of the dominant PFT is higher than observed.
3.3. Equilibrium solution for the generalized equations
For b = 0 and ɛ = 1, the competition–colonization equations for two PFTs with fractional coverages f1 and f2 are
and the equilibrium fractions are f1 = c1/(c1 + m1) and f2 = c2(1 − f1)/(c1 + c2 + m2), respectively. Equations (11) and (12) imply that as long as m1 > 0, then PFT 2 will always exist and equilibrium coexistence is possible. Equations (4) through (12) imply that relaxing the assumptions of random interaction in favor of full interaction, predicated on the availability of seeds, leads to coexistence of PFTs without invoking any additional mechanisms. The full interaction assumption is justified on the basis of availability of seeds owing to their long life spans and dispersion. Clark (Clark 1998), for example, illustrates that the fast migration rates at the end of the Pleistocene could not have been achieved without long-distance dispersal. He suggests that tree migrations were likely constrained by factors other than migration potential, implying that dispersion and thus seed availability were not a constraint. Pake and Venable (Pake and Venable 1995; Pake and Venable 1996) suggest that coexistence is promoted via the seed storage effect because long-lived seeds allow species to exploit favorable environmental conditions whenever they occur. Indeed in competition–colonization Equations (4) and (7), if b ≠ 0 (implying restricted seed availability) and a “parent” does not exist, then PFTs cannot expand despite favorable environmental conditions.
Analytical expressions for equilibrium fractional coverages for PFT 1 and 2 are not easily derived for noninteger values of b or ɛ ≠ 1 although, of course, any values of b and ɛ may be used for the numerical simulations discussed subsequently. Finally, the generalized competition–colonization equations are modified slightly to take into account the presence of crops in a grid cell. Noncrop PFTs are not allowed to invade agricultural areas. When the fraction of a grid cell covered by crops ( fc) is prescribed then the competition between noncrop PFTs for available space is restricted to the remaining fraction (1 − fc) of the grid cell.
3.4. The “start-up” problem
A difficulty with versions of the competition–colonization Equations (4)−(7) with b ≠ 0, including the usual LV version with b = 1, is the “start-up” problem. For b ≠ 0 the initial absence of a PFT in a grid cell means that dfα/dt = 0 and the PFT cannot expand. Cox (Cox 2001) overcomes the start-up problem by prescribing a minimum fractional coverage of 1% for each PFT termed the seeding fraction. For b ≠ 0 in Equations (4)–(7) a parent must exist before a PFT can colonize part of the grid cell. With b = 0, however, a PFT has the potential to expand into a grid cell if climatic conditions are favorable even if no parent PFT exists. Since seeds have long life spans and experience dispersion, they may exist at a location even in the absence of a parent and this is an argument for choosing b = 0.
3.5. Colonization and mortality rates
The PFT-dependent colonization rate for bare ground (day–1) is estimated as a function of daily net primary productivity (NPP; kg C m−2 day –1) as
where ξ is the inverse seedling density (m2 kg C−1), and Λ is the fraction of NPP allocated for expansion/reproduction. The term Λ NPP determines how much carbon is available for expansion, the inverse seedling density ξ determines the area that can be invaded with this available carbon, and their product ξ Λ NPP yields the colonization rate. Seedling density is assumed to be the maximum of a PFT-dependent fraction (μ) of the live biomass density (kg C m−2) and 0.5 kg C m−2. In a manner similar to that of Cox (Cox 2001), the fraction of NPP allocated for expansion/reproduction is given by
where Lα is the LAI of PFT α, Lmin is the LAI below which no expansion occurs, Lmax is the LAI above which the maximum fraction Λmax of NPP is used for expansion, and all these parameters are PFT dependent.
The PFT-dependent mortality rate m = ma + mg + md + mb represents the net effect of four processes: 1) intrinsic or age-related mortality ma (only a given fraction of trees exceed a characteristic specified age); 2) growth or stress mortality mg (negative or near-zero NPP values due to adverse climatic conditions lead to mortality); 3) mortality due to exogenous processes of disturbance md; and 4) mortality associated with the climate tolerances of PFTs mb.
Intrinsic or age-related mortality ma (Keane et al. 2001) is calculated as the annual death probability (yr−1) such that only 1% of the trees exceed their maximum age (amax)
For example, for only 1% of trees to attain or exceed a maximum age of 220 yr the annual death probability is 2.0%. Intrinsic mortality is taken to implicitly account for mortality due to hailstorms, ice damage, wind throw, and frost that are not explicitly modeled in CTEM. Growth- or stress-related mortality is modeled following Sitch et al. (Sitch et al. 2003) as
where ge is the growth efficiency (g C m−2), ΔC is the change in biomass over a period of a year (g C m−2), LAImax is the maximum leaf area index over the same year, and 0.01 yr–1 is the maximum growth mortality rate. The fire submodule of CTEM is used to calculate exogenous mortality due to disturbance (Arora and Boer 2005a) where area burned is a function of aboveground vegetation biomass, root zone soil moisture, and lightning frequency.
Mortality mb associated with the climatic tolerances of PFTs is prescribed at 0.1 yr–1 and is initiated only when long-term climatic conditions (averaged in an e-folding sense using a time scale of 25 yr) become unfavorable for a given PFT according to physiological bioclimatic constraints (Table 1) that are obtained following the LPJ model (Sitch et al. 2003). More complete dynamic vegetation models would not initiate mortality (as in this study) or exclude PFTs (as in LPJ, CASA, and IBIS models) on the basis of bioclimatic constraints expressed in terms of long-term climatic conditions. They would instead include explicit representation of biological processes that induce mortality on the basis of current conditions and plant “health.” However, our current understanding of plant physiological processes is inadequate to allow this in current DGVMs.
4. The Canadian Terrestrial Ecosystem Model
The competition–colonization equations described in section 3 are incorporated in the CTEM, which provides daily values of NPP, LAI, and biomass density for each of the nine PFTs in Table 1. CTEM contains three live vegetation pools for leaves, stem, and root and two dead carbon pools for litter and soil organic carbon as shown in Figure 1. CTEM is coupled to the Canadian Land Surface Scheme (CLASS 2.7) (Verseghy 1991; Verseghy et al. 1993) (Figure 2), and the two together produce fluxes of energy, water, and CO2 at the land–atmosphere boundary. Values of soil moisture, canopy temperature, and aerodynamic resistance provided by CLASS are used by the photosynthesis submodule of CTEM to provide estimates of canopy conductance and CO2 uptake. Fast biophysical land surface processes that affect the energy and water balance, including photosynthesis and evapotranspiration, are modeled on short time scales. The photosynthesis submodule is based on the biochemical model of Farquhar et al. (Farquhar et al. 1980), Collatz et al. (Collatz et al. 1991), and Collatz et al. (Collatz et al. 1992) and operates at the same time step as CLASS (∼30 min). The current version of the model uses a single-leaf photosynthesis approach, and coupling between photosynthesis and canopy conductance is based on vapor pressure deficit (Leuning 1995). Soil moisture and canopy and soil temperatures from the land surface scheme enter the calculation of the maintenance respiration of the three vegetation components and the heterotrophic respiration from litter and soil carbon pools. The photosynthesis, autotrophic, and heterotrophic respiration submodules of CTEM, as described in Arora (Arora 2003), are used to calculate net primary and net ecosystem productivities.
Positive NPP is allocated between leaves, stem, and root components based on light, root water, and leaf phenological status. The phenology submodule of CTEM uses a carbon-gain approach in which leaf onset is initiated when it is beneficial for the plant, in carbon terms, to produce new leaves. Leaf offset is initiated by unfavorable environmental conditions that incur carbon losses, and these include shorter day length, cooler temperatures, and dry soil moisture conditions. Leaf litter is generated due to the normal turnover of leaves, and cold and drought stress. Stem and root litter is generated based on prescribed PFT-dependent turnover times for the stem and the root components. The allocation, phenology, and turnover submodules of CTEM are described in Arora and Boer (Arora and Boer 2005b). Other than photosynthesis, all biogeochemical submodules of CTEM operate at a daily time step (Figure 2).
Allocation to, and the litter and respiratory losses from, the three vegetation components (leaves, stem, and root) results in time-varying biomasses. The change in biomass of the components is reflected in the structural vegetation attributes used in the energy and water balance calculations of the land surface scheme (Arora and Boer 2005b). In brief, the leaf biomass is converted to LAI via specific leaf area (SLA; m2 kg C–1). The aboveground biomass is related to vegetation height that affects the roughness length and therefore the turbulent transfer of heat and moisture. Finally, time-varying root biomass is related to an evolving root distribution profile (Arora and Boer 2003) used to calculate the fraction of roots present in the each soil layer. CTEM has been validated against observations of LAI, stem and root biomasses, and leaf onset and offset times for temperate and boreal deciduous, tropical evergreen, and tropical deciduous plant functional types in Arora and Boer (Arora and Boer 2005b).
5. Application and testing
CTEM and the competition–colonization component are tested at five locations (shown in Figure 3) that are characterized by different climate regimes and PFT distributions (or biomes), namely, Siberia (boreal forest); Amazon (tropical evergreen forest); Vancouver Island, Canada (temperate forest); Manitoba, Canada (boreal forest); and Botswana (savanna). We have confined our initial evaluation to these sites in order to be able to use information available from local studies regarding successional behavior of PFTs as their fractional coverages evolve in time, time to equilibrium after disturbance, and fire return intervals. The model is also tested by comparing equilibrium values of simulated fractional coverage of its PFTs with observation-based data. The observation-based data are obtained from two sources: the land-cover data of Wilson and Henderson-Sellers (Wilson and Henderson-Sellers 1985) currently used in the CCCma climate model, and the land-cover data based on the Global Land Cover (GLC 2000) initiative (Mayaux et al. 2004; Eva et al. 2004). Comparing equilibrium values of vegetation distributions with comparatively short-term averages of observed land cover has its difficulties, of course. In the real world, ecosystems are seldom fully in equilibrium and the fractional coverages of PFTs are continually adjusting to competition, natural disturbances, and the interannual variability in climate. Land clearing for establishing croplands, setting fire to grass and shrublands for fresh pasture, and establishment of managed forests for timber and paper also affect land cover. The observed land cover in the areas chosen nevertheless indicates whether the correct PFTs exist in a grid cell in reasonable proportion and provides insight into model behavior.
5.1. Climate and soil data
The data used to drive the coupled CLASS and CTEM models are obtained from Dirmeyer and Tan (Dirmeyer and Tan 2001) whose global land surface data (GOLD) set contains 6-hourly values of the required meteorological variables from 1979 to 1999 gridded at 3.75° resolution. These data are based on the National Centers for Environmental Prediction (NCEP) reanalysis but have been corrected for several known biases. The CCCma climate model is run operationally at horizontal resolutions of 3.75° and 2.8° so the resolution of the GOLD data is suitable for assessing CTEM's competition–colonization component where the atmosphere and soil data represent the average over a comparatively large grid cell. Since CLASS operates at a time step of half an hour, the 6-hourly climate data are disaggregated. Temperature, longwave radiation, wind speed, specific humidity, and pressure are linearly interpolated. Following Chen et al. (Chen et al. 1999) the shortwave radiation flux changes with solar zenith angle with a maximum value at solar noon. Precipitation data are disaggregated following Arora (Arora 1997) where the 6-hourly precipitation amount is used to estimate the number of wet half hours, and the total precipitation amount is randomly distributed. CLASS requires the fraction of sand and clay for each of its three soil layers together with total soil depth, and these are derived from the standard dataset used in the CCCma climate model based on Zobler (Zobler 1986) and Wilson and Henderson-Sellers (Wilson and Henderson-Sellers 1985).
The model is initialized with zero biomass and zero fractional coverage for all PFTs at all five sites, and the 21 yr of GOLD data are used repeatedly until the vegetation biomass and fractional coverages of all PFTs come into equilibrium. This spinup requires 50–150 yr depending on the vegetation type and the climate. At equilibrium the average fractional coverages of PFTs from the last 21 yr are compared with observation-based estimates. The parameters for the colonization–competition model used in these simulations are summarized in Table 2. The values of Lmin and Lmax are similar to those used by Cox (Cox 2001); Λmax is set to 0.10 following Sitch et al. (Sitch et al. 2003) such that a maximum of 10% of NPP is used for reproduction/expansion. The average maximum tree age generally increases from about 200 to 600 yr as elevation increases from 0 to 3000 m (Luo et al. 2004, their Figure 9); amax is set to 250 yr (corresponding to 100–500-m elevation range) for all PFTs following Luo et al. (Luo et al. 2004). Three simulations are performed for each site using values of parameter b equal to 1.0, 0.5, and 0.0. For simulations with b ≠ 0 a seeding fraction of 0.5% is used. In this initial approach a single value of ɛ = 0.5 is chosen for “intertree” competition while ɛ = 1 is used for tree–grass competition.
6.1. Siberia (boreal forest)
The evolution of the simulated fractional coverages of PFTs, with the three values of parameter b, for the Siberia grid cell, is shown in Figure 4. The LV version of the equations corresponds to b = 1. As the value of b decreases from 1 to 0, coexistence at equilibrium becomes increasingly possible. At equilibrium, for b = 1 and 0.5, only needleleaf deciduous trees and C3 grasses are simulated to exist while for b = 0 needleleaf evergreen and broadleaf cold deciduous trees also occupy portions of the grid cell. Although covering a small fraction, C3 grasses are able to coexist in the LV simulation for b = 1 because of the interannual variability in climate.
In Figure 4, the equilibrium values of the fractional coverages of PFTs can be compared with the two observation-based datasets, both of which show little crop area in this grid cell. Human influence is therefore likely to be minimal. When compared to observations, the simulation with b = 0 performs best in that it correctly simulates the existence of the four PFTs: needleleaf evergreen trees, needleleaf deciduous trees, broadleaf cold deciduous trees, and C3 grasses. Simulated fractional coverages of these four PFTs also compare reasonably well with observation-based estimates as shown. The LV model, on the other hand, leads to a higher-than-observed fractional coverage of needleleaf deciduous trees (∼78%) and a lower coverage of C3 grasses (∼6%). The results for the intermediate value of b = 0.5 fall in between. The case of C3 grass is somewhat difficult to assess since observation-based datasets exhibit considerable disagreement. In particular, the Wilson and Henderson-Sellers (Wilson and Henderson-Sellers 1985) dataset generally gives a higher fractional coverage of grasses.
The evolution of coverage, for all the values of b, shows an initial rapid expansion of C3 grasses that are eventually invaded by tree/woody PFTs. This behavior is similar to that of the LPJ model (Sitch et al. 2003, their Figure 2). The time needed for the grasses to reach their maximum coverage is ∼8 yr after the start of the simulation for b = 0, which is also comparable to Sitch et al. (Sitch et al. 2003). The LV equations, on the other hand, produce an unrealistic dominance of C3 grasses (covering more than 90% of the grid cell) for up to 50 yr after the start of the simulation. This is because of the dependence of colonization on the existing fractional coverage (i.e., colonization is restricted by the absence of parents) that results in an initial slow rate of increase of coverage for the needleleaf deciduous tree PFT.
Overall, the simulated fractional coverages of the four PFTs for b = 0, including the small fractional coverages of needleleaf evergreen trees and broadleaf deciduous trees, are in reasonable agreement with observations. The b = 0.5 case does somewhat less well while the b = 1 case does not allow for coexisting PFTs and the simulated fractional coverages do not compare well with observations.
6.2. Manitoba (boreal forest)
Figure 5 shows the evolution of fractional coverages of PFTs for the Manitoba grid cell using the three values of b and compares equilibrium fractional coverages with observation-based estimates. Comparison is difficult here because of the presence of the crops. The Wilson and Henderson-Sellers (Wilson and Henderson-Sellers 1985) dataset suggests a crop fraction of ∼17%, while the newer GLC 2000 dataset suggests a crop fraction of ∼33%. An additional simulation is therefore performed with a specified C3 crop fraction of 25% (an average of the two observed estimates) with b = 0.
Coexistence is increasingly possible as the value of b decreases from 1 to 0. At equilibrium, simulations with b = 0.5 and b = 0 correctly simulate the existence of the three PFTs: needleleaf evergreen trees, broadleaf cold deciduous trees, and C3 grasses, while in the LV model, the broadleaf cold deciduous PFT drives other PFTs out of competition. For the simulation with b = 0.5, the equilibrium fractional coverage of needleleaf evergreen trees is higher than both observation-based estimates, while fractional coverage of broadleaf deciduous trees compares well. The simulated equilibrium fractional coverage of needleleaf evergreen trees (∼44% in the simulation with b = 0) compares best with the estimate based on the GLC 2000 dataset (∼45%) but is higher than the Wilson and Henderson-Sellers (Wilson and Henderson-Sellers 1985) estimate (∼30%). Equilibrium fractional coverage of broadleaf cold deciduous tree (∼25%) in this simulation is also higher than observation-based estimates (∼2.5% to 9%). The equilibrium fractional coverages of PFTs are reduced, and their values compare better with observation-based estimates for the simulation with prescribed fractional coverage of C3 crop of 25% (Figure 5d). However, even this simulation has difficulties, since we do not model deforestation in a time-varying manner and do not have information about when and which PFTs were deforested in order to increase crop area. As in the case for the Siberia grid cell, the simulated equilibrium value of fractional coverage of C3 grasses (∼14%–20%) lies in between the two observation-based estimates for simulations with b = 0 and b = 0.5, while the LV version (b = 1) leads to unrealistic dominance of C3 grasses for more than 60 yr after the start of the simulation.
The model again simulates the rapid expansion of C3 grasses that are subsequently dominated by the tree/woody PFTs. Amongst the tree PFTs, in simulations with b = 0, the broadleaf cold deciduous trees colonize first and are later dominated by the needleleaf deciduous trees (∼50–60 yr after the start of the simulation), eventually leading to their coexistence. This behavior of the model is supported by observations from boreal forests after fire where early-successional deciduous trees give way to needleleaf leaf evergreen trees (black spruce) leading to their canopy closure at about 50 yr (Bond-Lamberty et al. 2004). The time when needleleaf deciduous trees dominate broadleaf deciduous trees is, however, longer in the simulation with b = 0.5 (∼70 yr). In the absence of explicit inter-PFT competition this successional behavior is not captured in the LPJ model (Sitch et al. 2003, their Figure 2 for the boreal grid cell), which simulates a gradual and simultaneous increase in fractional coverages of both needleleaf evergreen and broadleaf cold deciduous PFTs until they reach their equilibrium values after about 50 yr.
Overall, the model with b = 0 performs reasonably well in terms of successional behavior when compared to qualitative information from Bond-Lamberty et al. (Bond-Lamberty et al. 2004). It correctly predicts the occurrence of PFTs, and despite the caveat associated with presence of crops, the equilibrium values of noncrop PFTs lie within the range of observation-based estimates.
6.3. Amazon (tropical evergreen)
In Figure 6 for the Amazon grid cell, all simulations behave in a similar fashion, leading to dominance of the broadleaf evergreen tree PFT at equilibrium in agreement with observations. The transient behavior between initialization and equilibrium is, of course, different in the simulations. The simulated time to equilibrium for broadleaf evergreen trees in the simulation with b = 0 (∼25 yr) is smaller than that for tree PFTs in the Siberia and the Manitoba grid cell (∼40–50 yr) (consistent with the higher NPP of tropical forests) and is comparable to the 27-yr time scale achieved by the LPJ model for a tropical forest site (Sitch et al. 2003, their Figure 2). Kartawinata et al. (Kartawinata et al. 1981) reported that after 30 yr of natural regeneration in an abandoned pepper plantation in Indonesia, the resulting forest was physiognomically indistinguishable from the primary forest. Equilibrium times for simulations with b = 1 and b = 0.5, between 12 and 15 yr, are unrealistically lower than estimates from Kartawinata et al. (Kartawinata et al. 1981) and those based on the LPJ model (Sitch et al. 2003). In the simulations with b = 1 and b = 0.5, the full domination of broadleaf evergreen trees reduces the fraction of C4 grasses to its minimum prescribed seeding fraction, while with b = 0, C4 grasses are able to occupy a fraction of the grid cell (∼3%), which is slightly higher than observed. The simulated equilibrium fractional coverage of broadleaf evergreen trees compares reasonably well with observation-based estimates (∼93% to 96%) in all cases. In the model, C4 grasses capture almost the entire grid cell in about 1–2 yr after the start of the simulation, compared to ∼8 yr for C3 grasses in the Siberia and Manitoba grid cells. This is consistent with the highly productive nature of C4 grasses in warm and wet regions (Santos et al. 2004; Priante-Filho et al. 2004) that are able to establish themselves completely in about 2 years' time (Cerri et al. 2004).
6.4. Vancouver Island (temperate forest)
Simulated successional dynamics for the Vancouver Island grid cell are shown in Figure 7. The LV model (b = 1) leads to exclusion of all subdominant PFTs, and the simulation with b = 0.5 allows C3 grasses to coexist with needleleaf evergreen trees, while the simulation with b = 0 leads to coexistence of all the three PFTs: needleleaf evergreen trees, broadleaf cold deciduous trees, and C3 grasses. A further simulation is also performed with a prescribed fractional crop cover of 4% (an average of the two observation-based estimates). A reduction in value of b from 1 to 0, and the inclusion of crop fraction, leads to better agreement with observations. Model simulations (with b = 0) and observation-based estimates suggest that broadleaf cold deciduous trees occupy a small fraction (∼2%–5%) of the grid cell. This is because needleleaf evergreen trees perform better than broadleaf cold deciduous trees in regions where maximum precipitation and temperature are out of phase with each other (Stephenson 1998). Broadleaf cold deciduous trees do not perform well in regions where warm/hot summers are associated with little or no rainfall as is the case on the western coast of Canada. The dry summer growing season means that evergreen conifers outcompete deciduous trees because of their low transpiration rates in summer as well as their ability to take advantage of the moisture available during the wet winter months (Stephenson 1998). In CTEM, transpiration rates for conifers are lower than for deciduous trees (even for the same amount of photosynthesis) because of the lower value of coefficient m in the Leuning (Leuning 1995) stomatal conductance–photosynthesis formulation. As with other grid cells, the simulated fractional coverage of C3 grasses lies in between the two observation-based estimates but is closer to the estimate based on the GLC 2000 dataset.
6.5. Botswana (savanna)
In tropical and subtropical savanna regions the competition between trees and grasses is mediated by fire (Bond and Midgley 2000). Grasses promote fire in the dry season, causing tree mortality, and are able to capture the area available after fire due to their ability to colonize rapidly. Shorter fire-return intervals favor the expansion of grasses. We did not include fire in our simulations for Siberia, Manitoba, Vancouver Island, and Amazon grid cells because fire-return intervals at these locations are >400, ∼140, >500, and >500–1000 yr, respectively (Stocks et al. 2003; Thonicke et al. 2001), all of which are longer than the time to equilibrium. In semiarid tropical and subtropical regions, however, fire-return intervals are of the order of 1–12 yr (Thonicke et al. 2001) and fire plays an important role in the evolution of fractional coverage of PFTs. The fire module of CTEM (Arora and Boer 2005a) calculates the area burned that is distributed amongst existing PFTs in proportion to their fractional coverages and is subsequently available for colonization by all PFTs.
Figures 8a–c show the simulated successional dynamics for the Botswana grid cell without fire for the three values of the parameter b. Figure 8d shows results from a simulation with b = 0 that includes fire. Coexistence is increasingly possible as the value of b decreases from 1 to 0 and, in the simulation with b = 0, broadleaf dry deciduous trees, broadleaf evergreen trees, and C4 grasses are correctly simulated to exist although the simulated coverage is poor. Since precipitation shows distinct dry and wet seasons in Botswana, tropical broadleaf dry deciduous trees are able to perform substantially better than broadleaf evergreen trees. In Figure 8, model results for b = 0 and observations-based estimates show that evergreen trees occupy only a small fraction of the grid cell.
In Figure 8d, as expected, inclusion of fire increases the fractional coverage of C4 grasses and bare ground and decreases the fractional coverage of broadleaf dry deciduous trees, resulting in better agreement with observation-based estimates. Model values (for b = 0 with fire) and Wilson and Henderson-Sellers (Wilson and Henderson-Sellers 1985) estimates of bare fraction of ∼17% and ∼15%, respectively, are about half that of the GLC 2000 estimate (∼30%). The fractional coverage of C4 grass in the GLC 2000 estimate is correspondingly lower. This is consistent with the anomalously large number of fires in southern Africa (including northern Botswana) during the dry season of year 2000 (August–September) due to above-average high rainfalls the preceding growing season (December 1999–March 2000) that lead to prolific vegetation growth (Swap et al. 2003). Since our model is repeatedly driven with climate data from 1979 to 1999 we do not capture this event, although the simulated bare fraction can reach 25% in some years (Figure 8d).
7. Summary and discussion
The current generation of dynamic vegetation models has evolved from predecessors based on climate–vegetation equilibrium that relate the existence and coverage of PFTs to long-term average values of bioclimatic parameters. Many current dynamic global vegetation models (DGVMs) do not explicitly simulate competition but use these bioclimatic constraints in a moving-average sense together with other functional relationships to determine the existence and fractional coverage of PFTs. DGVMs that do simulate competition between PFTs typically use a version of the Lotka–Volterra (LV) equations that originated in the context of predator–prey interaction. The equilibrium solution of these equations yields a single dominant PFT, and subdominant PFTs are excluded. Although equilibrium may not be fully realized, coexistence of PFTs/species is widely observed in the real world. For example, the high biodiversity found in humid tropical ecosystems, which are in near equilibrium because of year-round favorable growth conditions, supports coexistence at near equilibrium.
Generalized competition–colonization equations are proposed that overcome some of the limitations of existing approaches and that include the LV equations as a special case. The assumption of random interaction and restricted seed availability (i.e., for the LV model with b = 1) generally leads to an unrealistic lack of coexistence, and relaxing this assumption (i.e., for b < 1) yields better agreement with observation-based estimates of sucessional behavior and evolution. Results with b = 0.5 yield coexistence and agree reasonably with available observations. The case with b = 0 avoids the requirement of a specified minimum seeding fraction that is necessary for b≠0 cases to overcome the “start-up” problem and is able to generate the observed mix of PFTs from initially bare ground. This version of the equations is also able to reasonably simulate successional behavior, time to equilibrium, and equilibrium values of fractional coverages of coexisting PFTs that compare reasonably well with observation-based estimates. The comparison of equilibrium values of fractional coverages of model PFTs with observed land-cover data can only be illustrative, however, since real-world ecosystems are seldom in equilibrium and humans may also modify the land cover. In this version of the equations, coexistence is obtained without invoking disturbance, random fruiting events, spatial heterogeneity in biotic and abiotic environments, or other special mechanisms. Thus, on the basis of the tests performed here and the consideration of seeding fraction, time to equilibrium and coexistence, the b = 0 version of the competition–colonization Equations (4)−(7) is chosen for inclusion in CTEM. The choice of b = 0 assumes that seed availability does not restrict colonization, an assumption also implicit in approaches that use bioclimatic constraints to determine PFT existence.
Of course, limitations remain with the proposed dynamic vegetation model. As is the case for other current dynamic vegetation models, bioclimatic constraints are used to induce mortality in PFTs that exist outside their climatic tolerance envelope. Ideally, the biological processes that induce mortality would be explicitly modeled on the basis of current conditions and plant state. The proposed dynamic vegetation model also does not take into account that certain areas of a grid cell may be unsuitable for colonization by plants (including rivers, lakes, roads, and rock and stony outcrops). The subgrid-scale heterogeneity of biotic and abiotic variables is not explicitly taken into account in our approach but could nominally be partially accounted for by empirically adjusting b and ɛ, although we do not attempt that here. CTEM is designed for use in a global coupled climate model where values apply to the average over large grid cells. This may affect the NPP of the PFTs and their successional dynamics. The use of grid average input also implies that the approach will miss the existence of PFTs that occupy small spatial climate niches, for example, the high elevation evergreen needleleaf forest that comprise a small fraction (∼3%) of Arizona's forests (O'Brien 2002).
Since simulation of the current vegetation distribution by DGVMs does not necessarily imply that these models will correctly predict future vegetation distribution (Woodward and Beerling 1997), they must be validated by other means. The successional behavior of PFTs and time to equilibrium after fire or other disturbances offers useful qualitative information in this respect.
In summary, the dynamic vegetation model developed here is based on a generalized version of the competition–colonization equations and is able to successfully simulate the competition between PFTs at large spatial scales. Colonization rates and tree/grass distinction determine the order of dominance, and colonization and mortality rates are functions of basic climate-dependent vegetation attributes (NPP, LAI, and vegetation biomass). Coexistence of PFTs is obtained in the “full interaction” version of the competition–colonization equations without invoking additional mechanisms. The explicit simulation of competition between PFTs allows their successional behavior to be modeled as fractional coverages evolve in time with changing climate. The proposed approach lies in between the inherently small-scale forest gap models and the large-scale quasi-equilibrium vegetation–climate functional relationship approach and offers a reasonable alternative for explicitly simulating inter-PFT competition and estimating time-varying fractional coverage of PFTs in dynamic vegetation models at large spatial scales.
We gratefully acknowledge Greg Flato, Charles Curry, Oleg Saenko, and Jiangnan Li for providing useful comments on an earlier version of this manuscript. This work was done as part of the Canadian Global Coupled Climate Carbon Modelling Research Network (CGC3M) supported by the Canadian Foundation for Climate and Atmospheric Sciences and Environment Canada. Comments from an anonymous reviewer are greatly appreciated and the help from the Chief Editor (Prof. Jon Foley) is also acknowledged.
* Corresponding author address: Vivek K. Arora, Canadian Centre for Climate Modelling and Analysis, Meteorological Service of Canada, University of Victoria, Victoria BC V8W 2Y2, Canada. Vivek.Arora@ec.gc.ca