Royal Society Publishing

Predicting ecosystem dynamics at regional scales: an evaluation of a terrestrial biosphere model for the forests of northeastern North America

David Medvigy, Paul R. Moorcroft


Terrestrial biosphere models are important tools for diagnosing both the current state of the terrestrial carbon cycle and forecasting terrestrial ecosystem responses to global change. While there are a number of ongoing assessments of the short-term predictive capabilities of terrestrial biosphere models using flux-tower measurements, to date there have been relatively few assessments of their ability to predict longer term, decadal-scale biomass dynamics. Here, we present the results of a regional-scale evaluation of the Ecosystem Demography version 2 (ED2)-structured terrestrial biosphere model, evaluating the model's predictions against forest inventory measurements for the northeast USA and Quebec from 1985 to 1995. Simulations were conducted using a default parametrization, which used parameter values from the literature, and a constrained model parametrization, which had been developed by constraining the model's predictions against 2 years of measurements from a single site, Harvard Forest (42.5° N, 72.1° W). The analysis shows that the constrained model parametrization offered marked improvements over the default model formulation, capturing large-scale variation in patterns of biomass dynamics despite marked differences in climate forcing, land-use history and species-composition across the region. These results imply that data-constrained parametrizations of structured biosphere models such as ED2 can be successfully used for regional-scale ecosystem prediction and forecasting. We also assess the model's ability to capture sub-grid scale heterogeneity in the dynamics of biomass growth and mortality of different sizes and types of trees, and then discuss the implications of these analyses for further reducing the remaining biases in the model's predictions.

1. Introduction

Terrestrial biosphere models are critical tools for inferring the current state of terrestrial ecosystems, and predicting changes in their composition, structure and function over the coming century. Previous model inter-comparisons have found that terrestrial biosphere models are able to replicate spatial patterns of potential vegetation and seasonal patterns of changes in regional atmospheric CO2, but that the different models diverge in their predictions of ecosystem composition, structure and function under novel climates [13]. One recent model inter-comparison [2] showed particularly striking discrepancies in North America, where climate–carbon–vegetation feedbacks ranged from being minimal to causing large increases in forest cover. Because North America is still recovering from large-scale deforestation during the past century, there is large uncertainty surrounding the magnitude of the continent's current net carbon uptake [14], and its ability to sequester carbon in the future.

A powerful method of reducing model uncertainty is to use observational datasets to estimate model parameters [5,6]. This approach is particularly relevant for terrestrial biosphere models because a significant number of key parameters and aspects of model formulation, such as those determining patterns of carbon allocation, are difficult to measure directly. In North America, observations spanning many spatial and temporal scales and sensitive to many processes are available for constraining and evaluating the performance of terrestrial biosphere models. These include eddy-flux datasets associated with the Ameriflux network of eddy-flux towers [7], forest inventory datasets from the USDA Forest Inventory and Analysis (FIA) and the Canadian Forest Service [8], and phenology data from the USA National Phenology Network [9]. Canopy reflectance data [10] and measurements of soil carbon [11] may also be used to reduce model uncertainty. In the near future, many other datasets will be collected as part of the National Ecology Observatory Network (NEON) and will be broadly available [12].

Early model-data synthesis studies [1315] focused on evaluating terrestrial biosphere model predictions against flux-tower measurements of net ecosystem exchange (NEE). Current model evaluation initiatives such as the North American Carbon Programme (NACP) model-data intercomparison [16] have continued to focus on NEE measurements, despite the current availability of other diverse datasets. Although NEE is a relevant diagnostic of current ecosystem function, an exclusive focus on NEE is problematic because NEE measurements by themselves are unlikely to provide sufficient information to adequately constrain terrestrial biosphere model predictions [14].

As far as we are aware, the Ecosystem Demography version 2 (ED2) biosphere model was the first to be simultaneously constrained using eddy-flux measurements and forest inventory measurements [17]. Specifically, these measurements included eddy-flux tower measurements of CO2 and H2O and forest inventory measurements of tree growth and mortality spanning a 2 year period at Harvard Forest (42.5° N, 72.1° W). Subsequent comparison with independent datasets showed that the constrained model formulation produced realistic estimates of eddy-fluxes, tree growth and mortality dynamics on timescales ranging from hours to a decade. The generality of the model was illustrated by comparing model predictions against ecosystem measurements at Howland Forest (45.1° N, 68.8° W), where the model realistically predicted the observed patterns of carbon fluxes and tree growth without further parameter adjustment, despite the vegetation composition being markedly different from that of Harvard Forest.

One reason that this approach to model parameterization and evaluation has not been adopted previously is that conventional terrestrial biosphere models use a ‘canopy-as-big-leaf’ approximation that limits their ability to connect to tree-level measurements of growth and mortality. In contrast, the structured biosphere model used by Penner et al. [18] comprises a system of partial differential equations that approximate the behaviour of a spatially distributed ensemble of individual plants [19], enabling it to predict the growth and mortality rates of individual plants, and thus directly connect to forest inventory measurements of tree demography.

The goal of this paper is to evaluate ED2's ability to predict the regional decadal-scale biomass dynamics of the forests of northeastern North America. In §2, we describe the USA and Quebec (QC) forest inventory measurements that we used to evaluate the model's performance, summarize the ED2 model, and give details on our numerical simulation experiments. In §3, we compare the predicted regional patterns of growth and mortality to the observations. In §4, we explore how the optimization improved the predictions and investigate the reasons for the remaining discrepancies between the model and the observations. In §5, we present our conclusions.

2. Methods

(a) Datasets

Decadal scale forest inventory measurements covering the past two decades are available for the northeastern USA and Quebec. In the USA, these forest censuses have been undertaken under the auspices of the national FIA programme [8,20]. While the sampling designs vary, the methodology used in the northeastern USA is relatively uniform: a grid of photo points overlaid on aerial photographs is used to generate a series of plot locations stratified by land use and timber volume; then, individual measurements of tree diameter growth and tree status (alive or dead) are made over a census interval of approximately 10 years for trees having diameters greater than 5 inches (12.7 cm). Plots are assigned weighting factors specifying their proportion of the total landscape sample. For further details see

In a similar way, Canadian forest inventories have been undertaken on the provincial level and have been used to obtain a detailed picture of Canadian forests [18]. In each on-the-ground inventory plot, all trees with a diameter at breast height (DBH) greater than 10 cm in a one-twentieth hectare area were measured and identified to species. Unlike the FIA inventory data, the Quebec inventory contains information on a substantial fraction of the understory. Within each plot, trees between 1 and 10 cm DBH were censused in a 1/200th hectare sub-plot. The species and the size of each stem sub-plot were recorded; however, the trees were not tagged, and thus were not re-censused.

The regional distribution of the forest inventory plots is illustrated in figure 1a. Over 16 000 plots were sampled in the US northeast of 40° N, 81° W, and over 11 000 plots were sampled in Quebec south of 52° N. Since most plots were inventoried twice, once in the mid-1980s and once in mid-1990s, it was possible to compute rates of growth, natural mortality and harvesting across the region.

Figure 1.

Forest inventory plots for the USA and Quebec. (a) Locations of USA (green) and Quebec (blue) inventory plots. There are over 16 000 plots northeast of 40° N, 81° W in the USA and over 11 000 plots south of 52° N in Quebec. (b) Locations of states and provinces. Abbreviations are expanded in the text. (c) AGB from forest inventories aggregated at 0.25° resolution. (d) Decadal-mean harvesting rates of AGB for the forested plots. Data are missing for all Canadian grid cells outside of Quebec as well as for a few grid cells in northern Quebec; these grid cells are coloured black.

(b) Model description

ED2 is a terrestrial biosphere model providing a physically and biologically consistent framework for both short-term and long-term studies of terrestrial ecosystem dynamics. The model uses a system of size- and age-structured partial differential equations to closely approximate the first moment behaviour of a corresponding individual-based stochastic gap model to realistically represent fine-scale heterogeneity in canopy structure within each grid cell (figure 2). The biomass density of trees of type i, of size vector z existing in a patch of age a at time t is denoted Ci(z,a,t) and evolves according to:Embedded Image 2.1

Figure 2.

ED2 model structure and processes. (a) Each grid cell is subdivided into a series of tiles. The relative area of each tile is determined by the proportion of canopy-gap-sized areas within the grid cell having a similar canopy structure as a result of a common disturbance history. (b) ED2 computes the multi-layer canopy fluxes of water, internal energy and carbon within each sub-grid scale tile. (c) Summary of the long-term vegetation dynamics within each tile arising from the integration of short-term fluxes shown in (b).

Note that z has three components corresponding to the sizes of the plant's active, storage and dead biomass pools. Plant DBH and height are computed from these biomass pools via allometric relationships. The growth rate of each compartment of z is given by gi(z,a,t), while the mortality rate is μi(z,a,t).

The changing landscape age-structure and associated sub-grid scale heterogeneity arising from prior disturbance history is tracked in the model by:Embedded Image 2.2where pj(a,t) is a vector whose elements describe the distributions of times since disturbance for each of the m land-use states represented in the model, and Λjk is an m × m matrix whose elements describe the disturbance history forcing—the rate at which land transitions between the different land-use states as a function of time and age (time since last disturbance).

Under the assumption that a fraction ri of seeds is randomly dispersed between gaps within a grid cell and that all other seeds remain within their gap of origin, the recruitment of new seedlings fi(z,a,t) corresponds to a flux of individuals into the system at (z0,a), yielding the following boundary condition:Embedded Image 2.3

Equation (2.1) has a second boundary condition that describes the state of the ecosystem after the disturbance event:Embedded Image 2.4where the function si,k(z,a,t) describes the survivorship of individuals of size z and type i following a disturbance of type k. The rates of growth, mortality and recruitment within canopy are computed from a corresponding vertically and horizontally stratified set of equations that determine the land–atmosphere exchange of carbon, water and energy between the ecosystem and the atmosphere (figure 2). For a complete description of the model, see [17].

Following Albani et al. [21], five empirically defined plant functional types (PFTs) differing in their physiology and life-history characteristics were used to represent the range of vegetation found within the northeastern USA: three deciduous broad-leaved tree types (early-, mid- and late-successional), and two coniferous tree types (northern pines (NP) and late-successional conifers (LWC)). These PFTs exhibit correlated differences in their leaf physiological traits, such as specific leaf area, leaf longevity, foliar C : N ratio, photosynthetic rate per unit leaf area (Vm0), and in their growth, mortality and recruitment of life-history characteristics (see electronic supplementary material, table S1).

(c) Numerical experiments

We evaluated the ED2 model's regional-scale predictions of above-ground biomass (AGB) dynamics by conducting two simulations for the northeastern USA and southern Quebec at 0.25° × 0.25° resolution for 1983–1995, a time interval that corresponds to the approximate period of the measurements. The first simulation, which we refer to as ‘ORIG’, was a literature-based parametrization used previously to conduct simulations for eastern North America [21]. The second simulation, which we refer to as ‘HF-OPT’, used the constrained model parametrization that was developed by fitting the model's predictions to 2 years of eddy-flux and tree census datasets from Harvard Forest [17]. The first 2 years of simulation (1983–1984) were discarded to avoid the potential impact of transient dynamics arising from the soil moisture initial condition. The model was forced with solar radiation, long wave radiation, temperature, humidity, precipitation, wind speed and pressure data from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-40 re-analysis dataset [22], while soil textural class for each grid cell was specified from the 1° × 1° resolution USDA global soil database, since higher resolution data were unavailable for Quebec. As in Medvigy et al. [17], we drove the optimized model with average phenology obtained from moderate resolution imaging spectroradiometer (MODIS) between 2001 and 2004 because region-wide phenology data were not available for this period. The forest inventory measurements were also used to calculate the spatial pattern of forest harvesting (figure 1d), which were then applied as a disturbance forcing to the model.

Initial conditions for the vegetation structure and composition at the beginning of the simulation were prescribed using forest inventory measurements from approximately 1985. Each inventory plot was represented as a sub-grid-scale tile (figure 2a), and individual trees were assigned to the closest PFT [17]. However, our analysis excluded the relatively small number of plots that occurred on highly atypical soils and plots for which growth and mortality information was not available. The resulting initial condition for total AGB is shown in figure 1c. There is a general latitudinal gradient, with values of about 50 tC ha−1 typical in the south decreasing to about 20 tC ha−1 in the north. Lower AGB levels near the Great Lakes and in Maine (ME) were exceptions to the general trend. Hardwoods had more AGB than conifers in all US states, but hardwoods and conifers have roughly equal amounts of biomass in Quebec (figure 3a; see also figure 1b for a map of state and province locations). Early-successional hardwoods (ESHWs; figure 3b) comprise 10–20% of the AGB in all states and provinces (figure 3a), and are the dominant PFT in central Quebec. Mid-successional hardwoods (MSHWs; figure 3c) are dominant in Pennsylvania (PA), southern New York (NY), Connecticut (CT), Massachusetts (MA) and parts of southern Quebec. Late-successional hardwoods (LSHWs; figure 3d) are particularly common in northern NY, Vermont (VT), New Hampshire (NH) and ME. NP are generally found in eastern MA, southern NH and southern ME (figure 3e). They do not dominate any of the grid cells. LWC are found mainly in the northern half of the domain and are the dominant species north of about 48° N (figure 3f).

Figure 3.

AGB for different PFTs derived from forest inventories. (a) Shows the proportion of AGB contributed by each PFT, broken down by state and province. The remaining panels show the spatial distributions of PFT AGB: (b) early-successional hardwoods (ESHWs); (c) mid-successional hardwoods (MSHWs); (d) late-successional hardwoods (LSHWs); (e) northern pines (NP); (f) late-successional conifers (LWC).

3. Analysis and results

(a) Forest growth dynamics

The observed pattern of forest growth within the region calculated from the forest inventory measurements is shown in figure 4a. Rates of biomass growth are highest in PA and western NY, and generally decrease as one moves north and east across the region, with the lowest rates of accumulation being found in northern Quebec and ME. The corresponding predictions of the original (ORIG) and constrained ED2 (HF-OPT) model formulations are shown in figure 4b,c, respectively. ORIG systematically over-predicts biomass growth rates across the entire region (figure 4b), although it qualitatively reproduces the observed spatial gradient of higher biomass growth in NY and PA and lower biomass growth in ME and Quebec. Its area-averaged bias and root mean square error (r.m.s.e.) were 0.68 and 0.84 tC ha−1 y−1, respectively. The systematic over-prediction of growth rates was largely corrected in HF-OPT. (figure 4c; bias = −0.06 tC ha−1 y−1; r.m.s.e. = 0.41 tC ha−1 y−1). In particular, HF-OPT reproduces the observed variation in patterns of biomass growth across the New England states and Quebec (CT, RI, MA, VT, NH, ME and Quebec) both qualitatively and quantitatively, although it under-predicts biomass growth in NY and PA. Figure 4d summarizes the overall spatial pattern, showing the observed and predicted average growth rates for each state/province.

Figure 4.

Annual-average growth rates. (a) Forest inventory data (missing data are coloured black). (b) Results from ORIG simulation. (c) Results from HF-OPT simulation. (d) State- and province-averaged growth rates: black bars, ORIG (initial model); dark grey bars, observed; pale grey bars, HF-OPT (optimized model).

(b) Forest mortality dynamics

The observed pattern of forest biomass mortality across the region is shown in figure 5a. The observed rates of biomass mortality are higher in southern Quebec than in northern Quebec, the New England states, NY, or PA. The corresponding predictions of the original and constrained ED2 model formulations are shown in figure 5b,c. The area-averaged biomass mortality predicted by ORIG had a bias of 0.07 tC ha−1 y−1 and an r.m.s.e. of 0.35 tC ha−1 y−1. While the ORIG simulation captures the high-biomass mortality rates in southern Quebec, it under-predicts biomass mortality rates in northern Quebec, and its predictions of biomass mortality rates to the south (in the New England states, PA and NY) are unrealistically high, with discrepancies as high as 0.5–1 tC ha−1 y−1 in the southwest of the domain. The biomass mortality rates predicted by HF-OPT (bias = −0.11 tC ha−1 y−1; r.m.s.e. = 0.28 tC ha−1 y−1) are generally lower than those of ORIG. These lower biomass mortality rates are a closer match to the observations throughout the south, but the mismatch is increased in northern areas. Figure 5d shows the observed and predicted average biomass mortality rates for each state/province. As the figure indicates, the ORIG formulation has reasonable mortality rates in Quebec but systematically over-predicts biomass mortality in all the US states except ME. The predictions in the southernmost states (CT, NY and PA) are two to four times higher than the observed rates. In contrast, the HF-OPT model formulation has substantially improved predictions for all of the US states, but under-predicts biomass mortality rates in ME and Quebec.

Figure 5.

Annual-average mortality rates. (a) Forest Inventory data. (b) Results from ORIG simulation. (c) Results from HF-OPT simulation. (d) State- and province-averaged mortality rates: black bars, ORIG; dark grey bars, observed; pale grey bars, HF-OPT.

(c) Sub-grid scale variability in growth and mortality dynamics

One of the distinguishing features of ED2's structured biosphere model equations (equations (2.1)–(2.4)) is that they predict AGB dynamics not only at the resolution of climatological grid cells, but also at the finer sub-grid scales that characterize biotic heterogeneity in canopy structure and dynamics. In this section, we evaluate the model's predictions of size- and PFT-related sub-grid scale variability in patterns of biomass growth and mortality dynamics. Our objective in doing so is to identify patterns of sub-grid scale variability in model bias that can assist in identifying mechanisms and processes responsible for mis-matches in predictions, For each grid cell, the relative bias in AGB growth (calculated as predicted growth/observed growth − 1) for each tree size class was computed, and then the bias values for each grid cell were stratified according to the observed mean growth rate. This grouping by size and by growth rate enables assessment of the accuracy of the model's predictions for plants with different levels of resource availability because the observed growth rate incorporates both the direct impacts of a plant's physical environment and its biological environment arising from local competition for light, water and nutrients.

The growth rate bias of ORIG is shown in figure 6a. As would be expected from figure 4, the bias is generally large and uniformly positive, with an average bias across all growth rate and size classes of 150 per cent. However, as the plot shows, the bias is larger for small trees, and is generally larger for trees growing in unfavourable environments (i.e. in sites with low growth rates). The corresponding growth rate bias for HF-OPT model formulation (figure 6b) shows that the bias is generally much smaller in magnitude than in ORIG: the average bias in HF-OPT is only 4.7 per cent. The bias of the trees in the largest size class has shifted from a positive bias to a negative bias of similar magnitude, but in all other size classes the bias is shifted towards zero. In terms of variation with respect to grid-cell productivity, HF-OPT outperformed ORIG in all productivity classes except those with the highest growth rates. As the figure indicates, the re-parametrization improved region-wide growth rates by reducing the growth for trees of all size classes in all environments, but particularly for trees in moderately productive sites. This remaining bias impacts northern Quebec, where growth is over-predicted, and also PA and NY, where growth is under-predicted. In the north, where observed growth rates are generally low, the predicted growth of large trees is reasonably accurate but the predicted growth rates of smaller trees remain too high. Conversely, in areas with high growth rates such as PA and NY, the predicted growth rate of small trees is reasonable, but the large trees are predicted to grow too slowly.

Figure 6.

Growth rate biases, 100% (simulated/observed − 1), for (a) ORIG and (b) HF-OPT. For each grid cell, the bias was computed for each diameter size class. Then, grid cells were grouped according the grid cell-level growth rate. The mean bias from ORIG is 150% and the mean bias from HF-OPT is 4.7%.

The dependence of biomass mortality on tree size and grid cell-averaged growth rate is shown in figure 7a,b for ORIG and HF-OPT. For small trees in sites with moderate to high rates of growth, ORIG over-predicts mortality. However, the magnitude of the over-prediction is reduced in larger size classes and in grid cells with lower growth rates. For trees with a DBH greater than 50 cm, ORIG under-predicts mortality regardless of grid-cell productivity class. HF-OPT significantly improves the biomass mortality rates of small trees in moderate- to high-productivity grid cells; however, there is very little difference between ORIG and HF-OPT for trees in grid cells with the highest productivity, where both model formulations over-predict biomass mortality. Mortality rates for trees with DBH greater than about 40 cm are consistently under-predicted. These results indicate that, in the generally high-productivity sites found in the southern part of the domain (figure 5), the improved biomass mortality rates in the constrained parametrization reflect a cancelling of errors: small trees experience excessive mortality, while larger trees experience too little. In contrast, in the generally low-productivity north, total mortality in the constrained model remained too low, primarily because of insufficient large-tree mortality.

Figure 7.

Panels (a) and (b) show the mortality rate bias, 100% (simulated/observed − 1), for ORIG and HF-OPT, respectively. For each site, the bias was computed for each diameter size class; then sites were grouped according the site-level growth rate.

(d) Biomass dynamics of different plant functional types

Growth: The ESHW PFT was parametrized using birch and aspen as canonical species [21]. The patterns of biomass growth largely reflected the pattern of occurrence (cf. figures 3b and 8a). The ORIG model formulation fared very poorly in predicting early hardwood growth, over-predicting by a factor of 2–3 (figure 8b and table 1). In contrast, the HF-OPT formulation had far more realistic biomass growth rates (figure 8c and table 1) throughout the region, though in a few areas such as northern and western NY, the HF-OPT growth rates were slightly lower than observed.

View this table:
Table 1.

Area-averaged biases and r.m.s.e.s for biomass growth. Results are listed for all trees, and then for individual PFTs. Grid cells not containing a particular PFT were excluded from the calculations corresponding to that PFT. Units are tC ha−1 y−1.

Figure 8.

(ac) Growth rates (tC ha−1 y−1) for ESHWs; (df) MSHWs; (gi) LSHWs; (jl) NP and (mo) LSC. The left column shows the forest inventory observations, the centre column shows the results from the ORIG simulation, and the right column shows results from the HF-OPT simulation.

MSHWs were parametrized on the basis of red oak and red maple. Observed growth rates decrease with increasing latitude (figure 8d). As with the ESHWs, the ORIG formulation substantially over-predicts MSHW biomass growth (figure 8e and table 1). HF-OPT's predictions are dramatically improved (figure 8f and table 1), especially in New England and Quebec. However, in NY and PA, the growth rates predicted by HF-OPT are lower than observed, with the magnitude of under-prediction comparable with the under-prediction of ESHWs. It is the under-prediction of these two PFTs that is responsible for the under-prediction of the biomass growth in the southwest of the domain seen in figure 4. LSHWs were parametrized on the basis of sugar maple and beech [21]. As with all other PFTs, ORIG over-predicted the biomass growth of the late-successional PFT (compare figure 8g,h and table 1), and this over-prediction was corrected in HF-OPT (figure 8i and table 1), which closely matched the observed pattern of late-successional growth (compare figure 8g,i).

Evergreen needleleaf PFTs were classed either as NP or as LWC (principally hemlock, spruce and fir). The observed biomass growth rates of the northern pine PFT is highest in areas of southern New England and Quebec, where it is relatively abundant (figure 8j). The ORIG formulation systematically over-predicts the rates of northern pine biomass growth (figure 8k and table 1), while HF-OPT successfully corrects the over-prediction of northern pine growth, yielding a spatial pattern that is quantitatively and qualitatively accurate (figure 8l and table 1). The observed growth rate of LWC increases with latitude, reaching a maximum in the northeast corner of the domain (figure 8m). The ORIG formulation systematically over-predicts the magnitude of late-successional conifer growth (figure 8n and table 1). The growth rates predicted by the HF-OPT formulation are more realistic (figure 8o and table 1), but are still over-predicted, especially in the north of the domain. It is this discrepancy in the biomass growth of the late-successional conifer PFT that is primarily responsible for HF-OPT's over-prediction of total biomass growth in northern regions (figure 4).

(e) Mortality

The observed biomass mortality of the early hardwood PFT was generally higher in the north (Quebec) than in the south (figure 9a). In contrast, the observed biomass mortality of MSHWs was highest in the southernmost part of the domain (figure 9d), and the LSHWs had their greatest biomass mortality rates in southern Quebec (figure 9g). In both ORIG (figure 9b) and HF-OPT (figure 9c), the ESHW biomass mortality rates are too low in Quebec, but their mortality rates in the US states match the observations more closely. In the ORIG model formulation, the mid- and LSHW biomass mortality rates are over-predicted across the US states (figure 9e,h, respectively). HF-OPT reduces these rates for both the MSHWs and LSHWs (figure 9f,i, respectively) to the level of the observations throughout much of the domain. The main exception is in western NY, where HF-OPT continues to over-predict mid- and late-successional biomass mortality; however, this is relatively small compared to the under-prediction of ESHW mortality in Quebec.

Figure 9.

Mortality rates (tC ha−1 y−1) for (ac) ESHWs; (df) MSHWs; (gi) LSHWs; (jl) NP and (mo) LWC. The left column shows the forest inventory observations, the centre column shows the results from the ORIG simulation and the right column shows results from the HF-OPT simulation.

The observed northern pine mortality rates (figure 9j) are largest in north central Quebec but are generally a small component of the total mortality. The late-successional conifer biomass mortality rates are highest in the northeast of the domain (figure 9m), where they are a significant component of the total AGB (figure 3e). The ORIG model formulation over-predicts northern pine and late-successional conifer biomass mortality in most of the USA, but under-predicts their biomass mortality in Quebec (figure 9k,n, respectively). In contrast, HF-OPT predicts more realistic mortality rates across the US states; however, the under-prediction of northern pine and late-successional conifer biomass mortality in Quebec increases (figure 9l,o, respectively). These results are summarized in table 2, which presents the area-averaged biomass mortality biases and r.m.s.e.s from ORIG and HF-OPT.

View this table:
Table 2.

Area-averaged biases and r.m.s.e.s for biomass mortality. Results are listed for all trees, and then for individual PFTs. Grid cells not containing a particular PFT were excluded from the calculations corresponding to that PFT. Units are tC ha−1 y−1.

4. Discussion

We found that the HF-OPT formulation of the ED2 biosphere model offered marked improvements over the ORIG model parametrization in predicting regional decadal-scale biomass dynamics in the northeastern North America. HF-OPT predicted AGB growth much more realistically than ORIG (figure 4), decreasing the region-wide bias from 150 per cent to less than 5 per cent, and also resulted in significant improvements in rates of biomass mortality across the region (figure 5). As discussed in more detail in Medvigy et al. [17], there were a number of important differences between the ORIG and HF-OPT model formulations that contribute to this improved predictive ability. In particular, the improvements in AGB growth in HF-OPT relative to ORIG reflected higher rates of fine root turnover, a shorter growing season for hardwoods and a lower maximum photosynthetic rate for conifers, which all reduce rates of AGB growth.

The predictions of sub-grid scale heterogeneity highlight how further improvements in the model's predictions of AGB growth will require simultaneously decreasing the growth rates of small trees, particularly in low-productivity sites, while increasing growth of high-productivity sites, especially for large trees (figure 6b). Several candidate mechanisms may explain this sub-grid scale signature of the remaining inaccuracies in HF-OPT AGB growth predictions. One source of model error may be the current model formulation's relatively simple representation of plant diversity [2325]. ED2 differs from other biosphere models in having several temperate deciduous and coniferous PFTs represented; however, the particular species that comprise each PFT class varies by location. For example, the predominant late-successional conifer species at Harvard Forest is eastern hemlock (Tsuga canadensis), while the dominant late-successional tree in northern Quebec is black spruce (Picea mariana), a species that is not found at Harvest Forest. Thus, the over-prediction of late-successional conifer growth and under-prediction of its biomass mortality, seen in figures 8 and 9, respectively, implies that ED2 may benefit from the introduction of a separate boreal conifer PFT. Similarly, in NY and PA, where growth rates are under-predicted, Prunus spp. comprise over 90 per cent of the ESHW PFT, while in the New England states and Quebec, the ESHWs are primarily birches (Betula spp.) and poplar (Populus spp.).

An alternative explanation for the above discrepancies is the absence of physiological, morphological or biochemical acclimation of PFT traits in different regions [2629]. Distinguishing between these two potential explanations is important because the time scales of ecosystem responses arising from acclimation responses within individuals to changes in climate forcing are likely to be considerably shorter than ecosystem responses arising from climate-induced shifts in the species-composition. Furthermore, the simulations did not account for the effects of nitrogen deposition and increasing CO2 concentrations. In particular, the absence of nitrogen deposition could potentially explain the under-prediction of growth rates in PA and western NY [30], where nitrogen deposition rates are high. We also expect that an improved radiative transfer scheme that accounted for the structure of neighbouring patches would likely reduce the amount of light captured by small trees, reducing their growth rates, and thus reducing the bias with the observations (figure 6b). The under-estimation of AGB growth of large trees may be corrected by allowing leaf nitrogen to vary as a function of an individual's position in the canopy [31,32], resulting in increased photosynthetic rates for the largest trees.

HF-OPT offered better mortality rate predictions than the default model, especially for small trees in the southern part of the region (figures 5 and 7). The total mortality rate in ED2 (equation 2.1) is the sum of a density-dependent rate and a density-independent rate. The density-dependent rate becomes large as trees fall into negative carbon balance. Because ORIG and HF-OPT have identical density-independent rates, the differences between their mortality predictions must be due to differences in the density-dependent mortality rates. The sites in NY and PA have relatively large AGB (figure 1c), and may also be expected to have large leaf area indices. The small trees in such sites are relatively shaded and, in the model, would have relatively large density-dependent mortality rates. Re-parametrization in HF-OPT decreased the impact of shading on mortality, though figure 7b suggests that it has not yet been reduced enough in the highest productivity sites. One way to address this problem would be to reduce leaf nitrogen for understorey trees, which would reduce leaf respiration and thus improve the carbon balance of these trees [33]. More realistic representations of crown shape and accounting for phototropism would also be likely to reduce mortality of shaded individuals [34].

Trees in larger size classes are not shaded, and so experience mainly density-independent mortality. In both the ORIG and HF-OPT, density-independent mortality was treated as a constant rate for each PFT, independent of size and was calibrated on the basis of the mortality of trees solely from the FIA. The Quebec inventory was not used to calibrate these mortality rates. In the highly productive sites of the northeastern USA, over-prediction of small tree mortality balanced the under-prediction of large tree mortality, leading to a realistic average. However, in Canada, biomass densities are lower (figure 1c), and so smaller trees are less shaded and do not experience excessive density-dependent mortality. Since the large trees continue to experience insufficient density-independent mortality, the overall result is an under-prediction of mortality. Several avenues will likely have to be explored in order to address this issue. First, model error may be reflecting the absence of climatic variation in density-independent mortality [35,36]. Allowing mortality rates to increase in areas experiencing cold winter, or springtime temperatures would increase mortality rates in Quebec relative to states like NY or PA. Second, ED2 did not explicitly include the impact of other sources of density-dependent mortality, such as insect outbreaks. In particular, omitting the impacts of the spruce budworm (Choristoneura fumiferana Clem.), for example, may have caused the model to under-predict mortality of LWC in Quebec [37].

5. Conclusions

We found that the HF-OPT model parametrization of the ED2 biosphere model developed by Medvigy et al. [17] offers a marked improvement over the ORIG literature-based model parametrization of Albani et al. [21] in predicting regional decadal-scale biomass dynamics in the northeastern USA and in Quebec. The biomass growth rates of all trees were reduced, and the mean model bias was nearly eliminated. Mortality improved mainly for small trees in moderate to high-productivity grid cells. These results of the analysis are encouraging as they demonstrate how short-term flux tower measurements, when combined with forest inventory measurements of forest demography can be successfully used to develop terrestrial biosphere model parametrizations that are general, not site-specific, and that these parametrizations can yield dramatic increases in the accuracy of a terrestrial biosphere model's predictions of long-term, large-scale, terrestrial ecosystem dynamics.

An important next milestone in assessing the predictive capabilities of terrestrial biosphere models is assessing their ability to capture inter-annual variability in biomass dynamics. The infrequent, approximately once per decade, re-measurement period in historical forest inventory data precludes an examination of inter-annual variability in AGB dynamics. However, the FIA has also recently been changed, and now the re-censusing of plots is staggered in time. This will help forest inventories to be sensitive to interannual variability, and, as a result, the regional impacts of climate and disturbances such as ice storms and pathogen attacks on vegetation growth and mortality rates should be accessible to observation. Evaluations of terrestrial biosphere model predictions of interannual variability in AGB dynamics against these datasets offer the promise of providing a critical assessment of whether current terrestrial biosphere model formulations have the accurate levels of climate sensitivity that are essential for accurate prediction of terrestrial ecosystem responses to future climate change.


We thank the USDA FIA programme and the Quebec Forest Inventory for providing the forest inventory measurements that were used in this analysis. These datasets were an invaluable source of information for assessing the regional-scale predictions of the ED2 model.



View Abstract