Aboveground carbon in Quebec forests: stock quantification at the provincial scale and assessment of temperature, precipitation and edaphic properties effects on the potential stand-level stocking

Study area

The study area corresponds to the forest territory below the current northern limit of the managed forest in Quebec, Canada, which extends from approximately lat. 45° to 52°N and covers approximately 583,078 km², of which 434,667 km² is classified as productive forest. This territory is characterized by three different forest subzones from south to north: the hardwood forest, the mixed forest, and the continuous boreal forest (Fig. 1A, MRN, 2013). The temperate mixed forest (lat. 47°to 48°N) marks the transition between the hardwood forest to the south, which is dominated by sugar maple (Acer saccharum Marsh.), and the coniferous forest to the north, which is dominated by balsam fir (Abies balsamea [L.] Mill) and black spruce (Picea mariana (Mill) B.S.P.). This mixed forest corresponds to the yellow birch (Betula alleghaniensis Britt.)—balsam fir bioclimatic domain. The study area is mainly characterized (51.4%) by medium-textured mineral soil deposits thicker than 25 cm and a xeric to mesic soil moisture regime (soil class 2, Table 1). Data from 387 weather stations distributed throughout the studied territory indicated that normal (1971–2000) mean annual temperature ranges from approximately from −2.6 °C to 7.4 °C, and annual precipitation, from 770 mm to 1,600 mm (Figs. 1B and 1C).

The intensity of logging activities in Quebec has changed spatially through time. At the beginning of the 20th century, they were confined to the southern portion of the province (south of lat. 49°N), whereas they currently take place up to lat. 51°N. Besides forest management activities, fire and spruce budworm (Choristoneura fumiferana Clemens) outbreaks were the main disturbances regulating forest dynamics over the last century (Girard, Payette & Gagnon, 2008).

Area-weighted average current carbon stock

Within a conventional forest map, the productive forestland is stratified into polygons representing stands with different forest attributes (composition, density, age, height) and site characteristics (soil deposit, drainage, slope) that are interpreted from aerial photographs (MRNF, 2009a). Forest attributes (number of stems by species and 2 cm diameter class) for each polygon were estimated from the compiled information of the corresponding forest inventory stratum (see next section). Provincial forest managers in Quebec have elaborated a mapping system (Système d’information forestière par tesselle; SIFORT) in which a conventional forest map (vector or object-oriented images) is translated into a grid of tiles (mixed vector and raster images) separated by 15″ (∼375 m) (MRNF, 2007). The latest forest map divided the territory in approximately 7.7 million of polygons while the systematic sampling of the whole area results in approximately 4.1 million tiles, of which approximately 3.1 million were classified as productive forests. Land use, site, and forest attributes of each tile correspond to the information of the polygon of the conventional forest map at the center of the tile. This systematic sampling of the conventional forest map results in a relatively high-definition raster map of forest attributes at the provincial scale which has been used to portray the contemporary evolution of the managed forest in southern Quebec over the last three rounds of forest inventory (MRNF, 2009b).

To assess the contemporary area-weighted average and total CACS in live tree biomass for the studied territory, we estimated CACS for each tile characterizing the studied territory based on the third round of forest inventory (1990–2002). We estimated the biomass of standing live trees using the Canadian national species-specific allometric aboveground biomass equations developed by Lambert, Ung & Raulier (2005), then converted the biomass data to C by assigning a C content of 0.5 Mg C per Mg oven dry biomass (IPCC, 2006). We estimated the total C stock in aboveground live biomass by summing the C content of each living stem and reporting the cumulative value on a hectare basis.

Carbon stock estimates at the plot level

Forest inventory programs conducted by the provincial forest authorities use temporary sample plots (TSPs) to portray existing forest resources comprehensively. The third round of forest inventory (1990–2002) provided data from more than 94,000 TSPs (Fig. 1D). Forest inventory of productive forest over the entire territory followed a stratified sampling design with proportional allocation. Forest stands interpreted from aerial photographs were first stratified based on stand characteristics (composition, density, height, age), edaphic properties (slope, drainage, deposit), and perturbations history (MRNF, 2009a). Circular plots (radius = 11.28 m, area = 400 m²) were then proportionally allocated in each strata according to their respective surface area (MRN, 2001). Within each plot, every tree with a diameter at breast height (DBH, measured 1.3 m above the highest root) larger than 90 mm was classified according to its species, and its DBH was measured to the nearest 2-cm class using a tree caliper. The number of saplings (11 mm ≤ DBH ≤ 90 mm) for each 2-cm class was also recorded by species in a subplot (radius = 3.57 m, area = 40 m²) (MRN, 2001). We estimated the biomass of standing live trees using the same method employed for C stock estimation at the scale of the forest map (see previous section).

Influence of climate and soil physical environment on MSAC

We simulated mean annual temperature and precipitation (1971–2000) for each plot with the stochastic weather generator of BioSIM software (Régnière, 1996). Its simulation models provide forecasts based on regional air temperature and precipitation, interpolated from nearby weather stations and adjusted for elevation and location differentials with regional gradients.

Climate appears to be the first factor determining forest biomass and forest distribution at the global scale. However, at the finer scales of landscapes and stands, other factors such as edaphic properties, species interactions and disturbances become important (Pan et al., 2013). Studied plots included widely variable stands in terms of age, structure, composition, and disturbance history. Many of the studied stands had been recently disturbed by forest management activities or natural disturbances and were not fully stocked, since they had not reached maturity. To document the influence of climate and soil physical environment on MSAC, we first stratified sampling plots into a total of 80 strata according to their soil physical environment (10 classes, Table 1), temperature (4 classes; <1 °C, 1–1.9 °C, 2–2.9 °C, ≥3 °C) and precipitation (2 classes; <1,000 mm yr⁻¹ and ≥1,000 mm yr⁻¹). Secondly, we selected the stands with C stock within the 95% confidence interval of the 90th percentile (Hahn & Meeker, 1991) of the C stock distribution for each of the 80 strata, assuming that they were mature, fully stocked stands in which CACS corresponds to the MSAC. We computed percentiles and related confidence intervals with the CAPABILITY procedure (SAS Institute Inc., 2011). Using this selection, we analyzed multiple regression models to investigate the relationships between mean annual temperature, annual precipitation and MSAC for each of the 10 soil classes. In order to determine the maximum amount of MSAC variance that could be explained as a function of temperature, precipitation, their squared value, and their interaction, we tested all possible subsets regression models with the REG procedure (SAS Institute Inc., 2011). We selected candidate models on the basis of two widely used criteria: Mallow’s Cp-statistic (Mallows, 1973) and Akaike’s information criteria (Akaike, 1973). Only subset models that had Mallow’s Cp values close to the number of parameters were considered. The final model selection was based on the lowest AIC. We tested for multicollinearity among the selected models’ climate variables using condition indexes and the variance inflection factor, to verify that dependencies among temperature and precipitation did not affect the regression estimates (Belsley, Kuh & Welsch, 1980). We also checked residuals graphically to assess the tenability of the regression assumptions.

The impact of climate change on MSAC

To assess the area-weighted average MSAC response to anticipated climate change, we compared the MSAC predicted from climate data under present and future climate scenarios using previously calibrated models at the plot level. For the present climate scenario, we used the BioSIM software (Régnière, 1996) to compute mean annual temperature and precipitation (1971–2000) for each tile. For the future climate scenario, we obtained climate projections for each tile using the delta change method (Fowler, Blenkinsop & Tebaldi, 2007) calculated from the output of global climate models (GCM), then applied them to the observed values for each tile. Future values were constructed using a large ensemble of GCMs made available by the World Climate Research Programme’s (WCRP’s) Coupled Model Intercomparison Project phase 3 (CMIP3) multi-model dataset (Meehl et al., 2007). Simulated climate data was available for present-day (20th century) atmospheric conditions and for projected future climate in response to three projected future greenhouse gas emission scenarios (SRES family A2, A1B and B1 scenarios; Nakicenovic et al., 2000) which have been endorsed by the IPCC and form the basis of their 4th assessment report (AR4) published in 2007. More recently, the publication of the IPCC 5th assessment in 2013 has caused a general shift away from the use of CMIP3 data in impact studies. This shift was particularly marked by a replacement of SRES greenhouse gas emissions scenarios by Representative Concentration Pathways (Moss et al., 2010) or RCPs. AR4 results remain valid however and comparisons between RCP and the SRES family scenarios can been made. Approximate equivalents (in terms of average global temperature) can be seen from the results of Knutti & Sedláček (2012) where SRES B1 is roughly equivalent to RCP 4.5 and SRES A1B is most similar to RCP 6.0, while SRES A2 is seen to be somewhere between RCP 6.0 and 8.5, showing stronger warming than RCP 6.0 but not as strong as RCP 8.5. Gleckler, Taylor & Doutriaux (2008) demonstrated that using the median or average of a large ensemble of climate simulations was more robust than using individual projections to reproduce the observed climate. This approach also has the advantage of providing a more robust indication of uncertainty of the projected future conditions. In all, we selected 68 climate simulations for two future horizons: 2041–2070 and 2071–2100.

As previously mentioned for the actual MSAC, our assessment of the future MSAC does not take into account the potential changes in natural disturbance regimes and forest management practices as well as the potential effects of increasing CO₂ concentration and nitrogen deposition.

Uncertainty analysis

Raupach et al. (2005) and Keith et al. (2010) identified many potential sources of errors that may cause uncertainties in C accounting across heterogeneous landscapes. These include the extent to which the field sampling covers the diversity of forest ecosystems across the landscape, error associated to the spatial extrapolation of site data, accuracy of tree allometric equations for estimating tree biomass, and variability in the C to biomass ratio of trees. We apportioned the main sources of uncertainty related to area-weighted average CACS and MSAC estimates with Monte Carlo uncertainty analysis (Yanai et al., 2010). First, we documented uncertainty due to tree biomass equations based on national forest inventory data by randomly sampling allometric model parameters and residual variance estimates values from their specified distribution (100 iterations, Lambert, Ung & Raulier, 2005) to generate a distribution of tree biomass estimates. We limited our analysis to 100 iterations because of the very large data set involved (171.8 million records). From these results, we calculated the mean total biomass and its standard deviation. Then, we generated a distribution of CACS estimates by combining total biomass and wood C content variability of North American hardwood and softwood species (10,000 iterations, Lamlom & Savidge, 2003).

Secondly, we documented the uncertainty related to the MSAC modeled from temperature and precipitation data along with uncertainty associated with future climate scenarios (100 iterations, Logan et al., 2011) to assess these sources of uncertainty related to area-weighted average estimates of MSAC under present and future climate. Uncertainty was reported as the coefficient of variation (standard deviation divided by the mean) of the distributions of Monte Carlo iterations. Our uncertainty analysis focused on the main sources of uncertainty related to CACS and MSAC computation, but did not take into account uncertainty in forest inventory measurements, their up-scaling to the entire territory, or the uncertainty associated to the interpretation of stand and site characteristics from aerial photography.

Current aboveground carbon stock estimates

CACS in live tree biomass for the whole 43.5-M ha study area (Fig. 2A) was estimated to be 1,633 Tg of which, because of its size, the boreal forest alone represents 56% (Table 2). This represents an area-weighted average CACS of 37.6 Mg ha⁻¹, including only land classified as productive forest. On average, forest stands in the hardwood subzone stocked approximately 41% more C per ha than those in the boreal forest ecozone. The area-weighted average CACS decreases from the hardwood forest in the south to the mixed and continuous boreal forests in the north (Table 2).

In sampling plots, CACS averaged 47.8 Mg ha⁻¹ and ranged from 0.03 to 100 Mg ha⁻¹ in 97.3% of the 94,268 forested sampling plots, although some reached exceptional values of up to 283.1 Mg ha⁻¹ (Fig. 3). The 18 most stocked stands (>200 Mg ha⁻¹) included old Norway spruce (Picea abies [L.] H. Karst.) plantations (one of the first non-native species introduced to North America) (39%), uneven-aged mature sugar maple stands (39%), and mature quaking aspen (Populus tremuloides Michaux) stands (22%). Although CACS was highly variable spatially, it was larger in warmer sites in the south of the province than in more northerly sites (Table 3). On average, CACS increased from 41.3 Mg ha⁻¹, for sites within the colder temperature class (<1 °C), to 58.0 Mg ha⁻¹ for sites within the warmer class (≥3 °C). On average, CACS was also slightly larger (49.6 Mg ha⁻¹) on wetter (≥1,000 mm) than on drier sites (<1,000mm, 45.3 Mg ha⁻¹). On average, the highest CACS values were observed on xeric to mesic sites with fine to medium-textured soils (soil class 2 and 3), while the lowest values were observed on hydric mineral soils and organic soils (soil class 7, 8, and 9).

MSAC and their relationships with climate

Table 4 presents the average CACS of the selected fully stocked stands (within the 95% confidence interval of the 90th percentile of the C stock distribution, and assumed to correspond to the MSAC) according to each stratum. The selected stands represent 2.9% of all plots. They vary widely in terms of composition (28% hardwood, 27% mixed, 45% coniferous), are mostly (90%) higher than 12 m, and generally (95%) have more than 40% forest cover. Depending on the stratum, average CACS of the fully stocked stands were 51.4–89.6% larger and generally less variable than CACS estimated from all plots (Tables 3 and 4). As previously mentioned, we assumed that these stands had reached their MSAC with respect to the present climate and their soil physical environment.

The MSAC was significantly related to both mean annual temperature and precipitation, or to the interaction between these variables (Table 5), for all soil classes except hydric mineral soils (soil class 7) and hygric coarse-textured soils (soil class 4), which were not related significantly to precipitation. Temperature and precipitation generally explained over 60% of the variance in MSAC, except for hydric mineral soils (40%, soil class 7) and bogs (48%, soil class 9). Overall, most of the variance explanation is due to the temperature gradient (Fig. 4). The explained variance was also higher (74–81%) for the soil classes 0–3 and 5, which cover the major proportion of the studied area as compared to the more marginal sites (40–64% of variance explanation for soil classes 4, 6, 7, 8, and 9) (Table 5). Overall, observed and predicted MSAC values for all the observations (n = 2, 775) were tightly correlated and well distributed along the 1:1 line (correlation coefficient = 0.91, Fig. 5).

MSAC response to climate change

Based on the previously defined MSAC models presented in Table 5, area-weighted average MSAC of individual stands over the whole territory (under the present climate: 1971–2000) was estimated to 72.5 Mg ha⁻¹ (Table 2, Fig. 2B). On average, future climate projections for the study area propose a median increase of air temperature of 2.4 °C in 2041–2070 and of 3.4 °C in 2071–2100 (Table 6). Concurrently, median annual precipitation is expected to increase by 7.5% in 2041–2070 and by 10.5% in 2071–2100. Based on the previously-defined models that relate MSAC and climate, under future climate scenarios the MSAC is expected to increase, on average by 20% in 2041–2070 and by 31% in 2071–2100 (Table 2, Fig. 2C). Regardless of the study period, average MSAC values on a per ha basis decreased from south to north and from warmer to colder climates, with the greatest values found in the hardwood forest subzone, followed by the mixed forest and continuous boreal forest (Table 2).

Uncertainty analysis

Uncertainty related to total and area-weighted average estimates of CACS live tree biomass due to allometric equations was 0.4%, while uncertainty related to average wood C-biomass fraction was estimated to 0.9%. When both sources of uncertainty were combined, uncertainty related to total (1,633 Tg C, Table 2) and area-weighted averages CACS estimates (37.6 Mg ha⁻¹, Table 2) was 1.0%. These results are coherent with reported uncertainty of 1% for total wood volume estimates over the same territory (Commission d’étude sur la Gestion de la Forêt publique Québécoise, 2004; Bureau du Forestier en Chef, 2013). The uncertainty related to the MSAC modeled from current temperature and precipitation data was 0.4% while the uncertainty associated with future MSAC estimates was 5.4% and 9.9% for the 2041–2070 and the 2071–2100 periods, respectively. Clearly, the variability among the climate scenarios was the main source of uncertainty in future MSAC estimates.

Current aboveground carbon stock estimates

In the last forest resources assessment report of the Food and Agriculture Organization of the United Nations (FAO, 2010) the average CACS in live biomass of the Canadian managed forest was estimated to 45 Mg ha⁻¹, while the Intergovernmental Panel on Climate Change reported an average aboveground C stock of 48 Mg ha⁻¹ for the American temperate and coniferous forest (excluding the forest tundra and very young forests (<20 years old), and assuming half of the biomass is C; IPCC, 2003). According to Fang et al. (2006) who conducted a literature review, inventory-based forest C stocks in middle and high northern latitudes ranged from 36 to 56 Mg ha⁻¹, and averaged 43.6 Mg ha⁻¹. These estimates are very close to the average C stock of 47.8 Mg ha⁻¹ we estimated from the sampling plot network. Our area-weighted average CACS estimate (37.6 Mg ha⁻¹), which accounted for spatial heterogeneity in forests and site attributes, was however 21% lower than the average C stock of 47.8 Mg ha⁻¹ obtained when not considering spatial heterogeneity and simply averaging the values of the 94,268 sampling plots (these plots followed a stratified, not a random, sampling). Such an evaluation (i.e., one that takes into account the biomass stock for each tile of a static grid covering the entire territory) can be used as a reference to estimate changes in aboveground tree C storage in the future, and also to document the role of forest and management practices on regional C stock and budget.

MSAC and their relationships with climate

Sampling plots represent a wide range of variability in terms of stand age, structure, composition, productivity, edaphic qualities, and disturbance history. Because of the large number of observations, we were able to categorize the stands within 80 climate–soil physical environment strata. We assumed that selecting only the most stocked stands, within each of these 80 strata, provided a good approximation of the MSAC as constrained by soil physical properties and climatic environment. Analyzing so many strata has the advantage of constraining the potential spatial variability in MSAC that could be associated with site edaphic qualities, productivity, vegetation types and disturbance histories. Nevertheless, some bias can exist for certain specific strata. For instance, on some marginal sites, the MSAC can be limited by non-climatic factors (e.g., edaphic qualities such as soil chemical properties), and thus be overestimated. Conversely, the MSAC can be underestimated for exceptional sites with very large MSAC that might stock more C than the confidence interval formed by the 90th percentile values of the C distribution, for a given climate–soil physical environment stratum. Despite these considerations, we are confident that our selection of the most stocked stands within a specific climatic zone for a given soil class represents a robust approach for estimating the MSAC of forest ecosystems.

The area-weighted average MSAC under the present climate was evaluated at 72.5 Mg ha⁻¹, being twice higher than the CACS (37.6 Mg ha⁻¹). Indeed, this is clearly a potential difference given that the MSAC definition implies a forest landscape composed of fully stocked mature stands with no natural disturbance regime and human intervention. Each year over the 2007–2012 period, large-scale fires and insect outbreaks respectively affected 0.3% (1,242 km² yr⁻¹) and 0.4% (1,831.2 km² yr⁻¹ severely defoliated) of Quebec’s managed forestland (Boulay, 2013). At the landscape scale, such a disturbance regime maintains a portion of the forest area in younger age classes. A modelling exercise based on available data of past stand-replacing disturbance histories (fire, spruce budworm outbreaks, windthrow) revealed that under the preindustrial natural disturbance regime, old-growth forests would compose, on average, 65% of the forest area in Quebec (49–86%, depending on the ecoregion considered) (Boucher et al., 2011).

The relationship between climate and MSAC

In good agreement with previous studies (see ‘Introduction’ section), we also found strong relationships between climate and MSAC for the various soil classes studied (Table 5). Temperature was responsible for most of the variation in MSAC, presumably because precipitation was relatively abundant within most of the studied territory. Overall, our data show the crucial importance of climate for controlling the MSAC of temperate and boreal forests, for which few results were available to date. Our results also suggest that along with climate, soil properties such as moisture regime and texture can influence the MSAC of forest ecosystems and weaken the MSAC–climate relationship. Indeed, MSAC was lower on marginal, poorly drained sites (hydric moisture regime, classes 7, 8 and 9) than on drier soil types. At the other end of the spectrum, fine to medium-textured mineral soil with xeric to mesic moisture regimes (soil classes 2 and 3) exhibited the greatest MSAC.

Projection of MSAC response to climate change

The projections show that forest MSAC has the potential to increase considerably in response to climate change. Area-weighted average MSAC over the studied territory is projected to increase from 72.5 Mg ha⁻¹ in 1971–2000 to 86.7 Mg ha⁻¹ in 2041–2070 and to 95.0 Mg ha⁻¹ in 2071–2100 (Table 2). Increased temperature is clearly the main driver of the projected increase in forest MSAC under future climate. In relative terms, the projected MSAC increase in 2071–2100 is greater for stands in the hardwood and mixed forests subzones (34–36%) than for stands in the boreal subzone (25%). Although they may appear surprisingly high, MSAC projections for stands in the temperate and boreal subzone over the 2071–2100 period compare well with contemporary MSAC estimates of mature forests located in areas characterized by climatic conditions similar to those expected in 2071–2100 for the studied territory (Hoover, Leak & Keel, 2012; Keith, Mackey & Lindenmayer, 2009; Liu et al., 2014). By 2071–2100, mean annual temperature is expected to reach 6.9 °C in the hardwood subzone with a mean annual precipitation value of 1,124 mm (Table 6). The projected MSAC for the hardwood subzone in 2071–2100 (119.7 Mg ha⁻¹) compares very well with the average aboveground C stock of 120.5 Mg ha⁻¹ found in live tree biomass of 12 old-growth forest stands in New England, USA (Hoover, Leak & Keel, 2012). Despite regional discrepancies in terms of nitrogen deposition, solar radiation and other physiological constraints, the studied sites in New England were characterized by climatic conditions similar to those expected in 2071–2100 for the hardwood subzone in Québec with mean annual temperature and annual precipitation ranging respectively from 6.1 to 7.6 °C and from 993 to 1,270 mm (Hoover, Leak & Keel, 2012).

The projected increases in MSAC, though quantitatively important, is a potential value that will not be reached. For instance, a part of the projected MSAC increase may require changes in forest composition and soil properties that are unlikely to happen within the next few decades (Lafleur et al., 2010). Also, the projected MSAC values do not take into account potential C losses that could be due to changes in disturbance regimes that may lead to an increase in burnt or defoliated areas, or to increase moisture stress associated to warming (Gauthier et al., 2015). Nevertheless, the large increase in MSAC values projected under future climate point to an important potential for Quebec forests to store more C in biomass in the future. Along with changes in climate and disturbance regime, forest management may have an important impact on future carbon stock at the landscape scale. This should encourage forest management practices that may optimize at the same time the carbon stocking of future forests, the production of long lived forest products, and forest resilience to climate change and natural perturbations.