Comparison of carbon-stock changes, eddy-covariance carbon fluxes and model estimates in coastal Douglas-fir stands in British Columbia

Abstract

Background

The global network of eddy-covariance (EC) flux-towers has improved the understanding of the terrestrial carbon (C) cycle, however, the network has a relatively limited spatial extent compared to forest inventory data and plots. Developing methods to use inventory-based and EC flux measurements together with modeling approaches is necessary evaluate forest C dynamics across broad spatial extents.

Methods

Changes in C stock change (ΔC) were computed based on repeated measurements of forest inventory plots and compared with separate measurements of cumulative net ecosystem productivity (ΣNEP) over four years (2003 – 2006) for Douglas-fir (Pseudotsuga menziesii var menziesii) dominated regeneration (HDF00), juvenile (HDF88 and HDF90) and near-rotation (DF49) aged stands (6, 18, 20, 57 years old in 2006, respectively) in coastal British Columbia. ΔC was determined from forest inventory plot data alone, and in a hybrid approach using inventory data along with litter fall data and published decay equations to determine the change in detrital pools. These ΔC-based estimates were then compared with ΣNEP measured at an eddy-covariance flux-tower (EC-flux) and modelled by the Carbon Budget Model - Canadian Forest Sector (CBM-CFS3) using historic forest inventory and forest disturbance data. Footprint analysis was used with remote sensing, soils and topography data to evaluate how well the inventory plots represented the range of stand conditions within the area of the flux-tower footprint and to spatially scale the plot data to the area of the EC-flux and model based estimates.

Results

The closest convergence among methods was for the juvenile stands while the largest divergences were for the regenerating clearcut, followed by the near-rotation stand. At the regenerating clearcut, footprint weighting of CBM-CFS3 ΣNEP increased convergence with EC flux ΣNEP, but not for ΔC. While spatial scaling and footprint weighting did not increase convergence for ΔC, they did provide confidence that the sample plots represented site conditions as measured by the EC tower.

Conclusions

Methods to use inventory and EC flux measurements together with modeling approaches are necessary to understand forest C dynamics across broad spatial extents. Each approach has advantages and limitations that need to be considered for investigations at varying spatial and temporal scales.

Background

Forests are a large component of the global carbon (C) stocks, containing an estimated 1146 Pg C (Dixon et al. 1994). Forest processes, which may be influenced by forest management, can therefore have a large impact on the global C budget, either by storing or releasing C. It is therefore critical to understand forest C dynamics, including forest C stock components and transfer mechanisms in order to develop accurate forest C models such as the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) (Kurz et al. 2009), 3PG (Landsberg and Waring 1997), and Ecosys (Grant et al. 2007), amongst others, as well as to inform forest management policy and for national and international reporting (Kurz et al. 2002). Factors affecting forest C dynamics include natural (e.g. fire, insect outbreaks) and anthropogenic disturbances from land use management and change (e.g. harvest, reforestation, deforestation) (Kurz et al. 2002) as well as weather (Morgenstern et al. 2004), fertilization (Jassal et al. 2010), and stand age and species composition, which can be related to disturbance history, site edaphic characteristics, or management practices (e.g. silvicuture) (Humphreys et al. 2006; Krishnan et al. 2009).

To measure the net exchange of C between land ecosystems and the atmosphere (net ecosystem exchange, NEE, with -NEE referred to as net ecosystem productivity, NEP), a global network of over 400 eddy-covariance (EC) flux stations has been established across a range of ecosystems, building an extensive data record of NEP, in some cases spanning up to two decades. The majority of these towers are located on sites not undergoing major disturbances so the fluxes measured reflect the interaction of weather, vegetation composition, stand age, and seasonal phenology. EC flux-towers use micrometeorological equipment to take measurements at the canopy scale of the exchanges of gasses including CO₂, water vapor, sensible heat, and some are equipped to measure other trace gases (Baldocchi 2008). The source area contributing to measurements made by the instruments mounted on the EC flux-tower, the flux-tower footprint, is variable in size and shape depending on height of measurement, surface roughness length, wind speed and direction, and atmospheric stability (Leclerc and Thurtell 1990; Schmid 2002), and proper interpretation of EC flux-tower based measurements depends on the flux-footprint over which the fluxes are sampled (Chen et al. 2008). Early estimates modelled flux-footprints as simple ovals. Recently, estimates of flux-footprint climatology may demonstrate more complex geometries and continuous probability density surfaces to quantify the upwind distribution of weighting factors over long time periods (Chen et al. 2009).

Forest measurements can be taken to quantify forest C stocks (Dixon et al. 1994; Clark and Brown 2001; NFI 2008), and by measuring forest C stocks at two times with a sufficient interval between measurements, the C stock change (ΔC) can be determined and used to estimate the cumulative net uptake of C from the atmosphere to the forest for the period of measurement (ΣNEP) (Clark and Brown 2001). Measurements are typically made of individual plants, for example, all trees within a plot are measured for diameter, species, and height, then allometric relationships are used to determine tree mass (Ter-Mikaelian 1997). Similarly, woody debris pieces are measured along transects using line intersect sampling and allometric relationships are applied to find biomass (Brown 1971). Other components are made by direct measurements, for example, understory vegetation can be sampled, dried, and weighed (Bailey et al. 1998). To measure soil C, field samples are typically dried and weighed for calculation of bulk density, and sent for laboratory analysis of C using dry combustion (Janzen 2005). Soil C measurements are important because the underground processes driving soil respiration form a large component of ecosystem gas exchanges, but these are less well-understood and measured compared to above ground stocks and processes (Ehman et al. 2002). Changes in soil C storage over time have implications as a source or sink for the global C cycle, and as an indicator of environmental function and health (Janzen 2005).

One application of forest measurement data to provide accounting of forest management actions and subsequently inform forest management decisions is as an input to C budget models such as CBM-CFS3. CBM-CFS3 is a forest C budget model that utilizes growth and yield information for biomass and generates explicit simulation of dead organic matter dynamics (Kurz et al. 2009). In contrast, forest process models use approaches based on the understanding of processes such as photosynthesis; these process-based models have more potential to model changing conditions because they do not rely on historic growth and yield and therefore are better suited to simulate forest conditions under global change situations (Landsberg and Waring 1997). CBM-CFS3 was designed to meet reporting requirements for forest management, provide policy support, and a function as a C budgeting tool for operating foresters, thus functioning at a range of spatial scales, and estimating records of C stocks, transfers between pools, and emissions (Kurz et al. 2002; Kurz et al. 2009).

Comparing EC flux-tower cumulative NEP (ΣNEP) and inventory ground plot measurements of C stock changes (ΔC) over the same period can help inform forest C processes for several reasons. First, ground plot measurements can serve as an independent validation of ΣNEP measurements made at EC flux-towers. Second, inventory plot data can provide more detailed information about the stand structural changes that may be driving variation stand-level EC C exchange. Third, broad stand inventory measurement datasets are available for a greater spatial extent than the global EC flux-tower network. For example, the Canadian National Forest Inventory (NFI) samples approximately 22,000 photo plots locations across Canada with detailed ground plots measured at over 1000 locations (Gillis et al. 2005). However, comparing inventory plot changes and EC flux-tower measurements requires comparisons to be made across different spatial and temporal scales. Therefore, such comparisons pose methodological challenges beyond individual inventory plot and EC flux-tower measurements.

Comparisons between EC tower measurements of ΣNEP and inventory plot measurements of ΔC stocks have been completed at a number of sites globally (Schulze et al. 2000; Granier et al. 2000; Law et al. 2001; Barford et al. 2001; Ehman et al. 2002; Curtis et al. 2002; Miller et al. 2004; Black et al. 2005; Gough et al. 2008; Kominami et al. 2008; Yashiro et al. 2010; Gielen et al. 2013; Babst et al. 2014), and at two sites of the Canadian Carbon Program (CCP), Boreal Ecosystems Research and Monitoring Sites (BERMS) station, located in the boreal forest of northern Saskatchewan (Theede 2007). Over a ten-year interval (1994 – 2004), Theede (2007) found a convergence of 15.6 ± 4.0 MgC ha⁻¹ 10⁻¹ years (inventory plots) and 18.2 ± 8.09 MgC ha⁻¹ 10⁻¹ years (EC tower) at the Old Aspen site and 5.8 ± 2.0 MgC ha⁻¹ 10⁻¹ years (inventory plots) and 6.9 ± 1.6 MgC ha⁻¹ 10⁻¹ years (EC tower) at the Old Jack Pine site. However, the homogeneous stand structure contributing to canopy-level EC flux measurements and large fetch in the boreal forest reduced the need for consideration of spatial vegetation structure and footprint distributions.

These previous studies have contributed considerably to the understanding of forest C dynamics; however, at sites with more complex vegetation structure and topography, such as the forests found in coastal British Columbia, the spatial distribution of vegetation structure and the footprint distribution are important considerations for accurate comparisons of EC flux-tower and inventory plot measurements (Schmid and Lloyd 1999). More complex ecosystems, therefore, require more complex models and methods developed for accurate comparisons of EC flux-tower measurements and inventory-based measurements. For example, a limited number of studies, including Ehman et al. (2002), Chen et al. (2009), Ferster et al. (2011), and Gielen et al. (2013) applied weighting by EC flux-tower footprint climatology to inventory-based measurements, where areas of the footprint that contribute more to tower-measured flux are weighted more heavily than other areas, resulting in an overall improvement when compared to simple flux footprint geometries.

In this paper, the difference between 2002 and 2006 inventory plot measurements of ecosystem C stocks was determined and the ΔC stocks compared to CBM-CFS3 modelled and EC-flux-tower estimates of ΣNEP at the CCP coastal British Columbia Flux Station. Through a combination of variable topography and a complex disturbance history (Trofymow et al. 2008), the EC flux-tower sites possess fine-scale spatial variation in forest structure, thus presenting a challenge to the interpretation of EC flux data (Schmid and Lloyd 1999; Göckede et al. 2004). To allow the comparison between measurements from the inventory plots, C budget models, and EC flux-towers, our approach consisted of four steps: first, we found the change in C based on forest inventory ground plot measurements made in 2002 and 2006 (four year period); second, we developed an additional approach that utilized litterfall data and published decay equations, for a second inventory-based estimate of C stocks change; third, we defined stand attributes at a fine spatial resolution, stratified the site based on stand attributes, and weighted the two inventory plot estimates; and fourth, calculated ΔC. These were compared with CBM-CFS3 model-based estimates of ΣNEP and EC flux-tower-measured ΣNEP for four years (2003, 2004, 2005 and 2006) in the same period. Finally, we evaluated how the forest age and conditions at the different sites accounted for convergence or lack of convergence in estimates for the different methods, and discussed the relative advantages and constraints of each measurement method.

Methods

Study area

The coastal British Columbia Flux Station sites are located on the leeward central east coast of Vancouver Island (Figure 1), in the Coastal Western Hemlock biogeoclimatic zone (Green and Klinka 1994) with a mean annual air temperature of 10° Celsius. Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco) is the dominant species, with lesser amounts of western hemlock (Tsuga heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata Donn ex D. Donn), and red alder (Alnus rubra Bong.). Three eddy-covariance and meteorological towers were placed in stands of different ages (Table 1).

Figure 1

The Canadian Carbon Program Coastal British Columbia Flux Station consists of three sites with flux towers (DF49, HDF00, HDF88) as well as ancillary inventory ground plots (HDF90) on the central east coast of Vancouver Island. A 5 x 5 km area (red square) spanning the Oyster River area has been used for spatial modelling studies.

Table 1 Study site locations and characteristics

EC flux-tower footprints, which delineate the spatial distribution of upwind source-weighting factors, were calculated by Chen et al. (2009) to represent the size, shape, direction, and magnitude of the flux as a function of wind speed and direction, measurement height, and surface roughness. Chen et al. (2009) calculated hourly footprints for 10 m by 10 m cells, weighted by NEP, and averaged over the 2002 – 2006 measurement period to determine the footprint climatology for the sites. In this study, the flux-footprint climatology 85% cumulative flux probability boundary was selected as the maximum extent of analysis, since the majority of the probability density is concentrated within this area and the remainder extends across a large area (Figure 2).

Figure 2

Flux footprints and inventory plot locations for the A) DF49, B) HDF88, C) HDF00 sites. Background imagery is pan-sharpened 2004 Quickbird for A and C, and true-colour orthophoto for B.

Inventory ground plots

All ground plots were established and measured by the Canadian Forest Service at the four sites following Canadian National Forest Inventory (NFI) guidelines and protocols (NFI 2008) though they are not part of the primary NFI ground plot network. Plot locations were chosen to represent each site series class (based on air photo interpretation and field transects) within the preliminary flux footprint boundary estimated in 2002 (Figure 2). Establishment and initial measurements were made in September 2002 and re-measurements were made in September 2006, constituting a four-year measurement interval. NFI data-compilation software used to determine plot-level values through application of allometric equations to overstory trees and density values for woody debris NFI (2008). Since live overstory vegetation roots were not measured in 2002 or 2006, coarse and fine root C mass was estimated from overstory biomass using equations by Li et al. (2003).

Live biomass

Live trees greater than 1.3 m tall were measured in an 11.28 m radius circular plot. Where trees were very numerous (for example, at HDF88), a half plot or 5 m radius circular plot was measured following NFI guidelines (NFI 2008). Allometric equations are used in the NFI compiler to estimate stem, bark, branch, and foliar biomass based on work by Lambert et al. (2005). These biomass values were converted to C content by assuming 52, 56, 52, and 52% dry C concentration, respectively (Matthews 1993).

Understory vegetation from four 1 m × 1 m microplots was destructively sampled and sent to the laboratory for drying and weighing for determination of mass which was converted to C content assuming 50% dry C concentration.

Detritus

Woody debris, fallen woody material > 1 cm diameter, was measured for diameter, species, and decay class along four 15 m transects at each plot. Decay class was determined by field crews in the five-class ordinal rating system utilized by the NFI, from intact (Class 1) to highly decomposed (Class 5) (NFI 2008). NFI compilation routines utilize algorithms for line intercept sampling to determine plot level volume and a lookup table of woody debris density values by species and decay class to determine mass. C content was determined assuming 50% dry C concentration of woody debris mass. Fine woody debris, fallen woody material 2 mm to <1 cm diameter, was sampled from within the four 1 × 1 m understory microplots in each sample plot and sent for laboratory drying, weighing, and determination of mass, and C content assuming dry mass is composed of 50% C. Dead standing trees were measured for diameter, height, and species and decay class was observed and recorded by field crews. Stumps were measured for diameter and decay class in 2002 but were not remeasured in 2006.

Surface substrates were measured along four 15 m long transects to determine the average depth and percent area surface coverage in the entire plot. The organic litter, fibric, and humus layer (forest floor) was sampled and average depths measured from within four 20 × 20 cm templates (one at each microplot). These destructive samples were collected, sieved, dried, and weighed to determine bulk density and subsamples sent for laboratory analysis of C concentration.

Mineral soil

In 2002, <2 mm mineral soil was sampled for bulk density and C concentrations from 10 – 12 cm diameter excavated holes at three depth intervals 0–15 cm (4 samples), 15–35 cm (2 samples), and 35–55 cm (one sample) and % volumetric coarse fragments determined from one 55-cm-deep soil pit per ground plot. Data were scaled to the entire groundplot discounting for the area without mineral soil (i.e. exposed bedrock). In 2006, the 0–15 cm layer was re-measured for C content at four locations in each plot using a 2 cm diameter soil corer. Samples were stored in plastic bags at 2°C prior to laboratory processing. For a detailed description of soil sampling methods see NFI (2008). For 2002–2006 soil C stock changes, only the 0–15 cm layer was considered, soil C at deeper depths was assumed unchanged.

Changes in inventory ground plot C stocks

Changes in C stocks over the four-year interval were calculated from the difference between C stocks in 2002 and 2006. The total change in C stocks was calculated as:

$$ \varDelta \mathrm{C} = \varDelta {\mathrm{C}}_{\mathrm{L}} + \varDelta {\mathrm{C}}_{\mathrm{D}} + \varDelta {\mathrm{C}}_{\mathrm{S}} $$

where ∆C_L is the change in total live biomass C stocks including live trees, and shrubs, herbs, and bryophytes; ∆C_D is the total change in detrital C stocks including dead standing trees, woody debris, and the forest floor; and ∆C_S is the change in mineral soil C stocks (0–15 cm).

∆C_D2 Detritus annual accounting method

Preliminary examination of the ∆C_D values showed larger-than-expected changes in the measured C stocks over the four-year interval, especially for the forest floor component. To evaluate these measured changes, ∆C_D was estimated using a second method (referred to as ∆C_D2) that accounted for the annual inputs into the detrital pools from litterfall and mortality (i.e. the transfers from live biomass to forest floor and soil , and live trees to dead trees respectively) and annual losses and transfers from decomposition estimated using published equations and parameters (Smyth et al. 2010; Kurz et al. 2009).

Inputs to detrital C pools from 2002 – 2005 included overstory fine litterfall (needles, leaves, cones, twigs) which was collected quarterly in 3 0.189 m² conical-mesh litterfall traps located in each ground plot at DF49, HDF88 and HDF90. Annual litter fall masses were calculated and converted to C assuming 50% dry C content. Herbaceous material, which was not tall enough to be sampled using litterfall traps especially at the HDF00 site, was collected and measured at the sample plots in 2002 and 2006, linearly interpolated between measurement dates, and assumed to turn over at a rate of 80% annually. Shrubby litterfall material was divided into stem and branch and foliage components, and assumed to turn over at a rate of 3% and 95%, respectively similar to the fine branch wood and foliage for deciduous trees (Kurz et al. 2009). When trees died within the measurement interval, the biomass was transferred to the detrital pools half way through the measurement interval. Finally, the biomass of any trees (live or dead) that were measured in 2002 and not recorded by the field crew in 2006 was transferred to the detrital pools half-way through the measurement interval under the assumption that the tree fell during the measurement interval.

Losses through decomposition and transfer to the soil organic C of each detrital pool was calculated at an annual time step with an annual average temperature of 10°C using coefficients that were developed and validated nationally (Smyth et al. 2010) for the CBM-CFS3 model (Kurz et al. 2009). The transfers and coefficients used are presented in Figure 3. Soil C stock changes were assumed negligible, and ∆C₂ was calculated as: ∆C_{2 =} ∆C_L + ∆C_D2.

Figure 3

∆C₂ stocks inputs (+) and outputs (−). Detrital pools decomposition and transfer coefficients from Kurz et al. (2009) and Smyth et al. (2010).

Spatial distribution of C stocks and ∆C at the stand scale

To account for spatial heterogeneity due to stand structure and topography, the methodology developed by Ferster et al. (2011) was applied to the ground plot data at the flux tower sites. Following this approach, predictor variables from GIS, topography, and remote sensing data were used to impute measurements from the inventory ground plots across the forest site based on the three-most-similar-neighbours (K-MSN with K = 3) (Crookston and Finley 2008) using the Mahalanobis distance as a measure of similarity (Mahalanobis 1936). For each footprint cell, the detailed inventory target measurements were estimated as the weighted mean of the three nearest neighbours (inversely weighted by the Mahalanobis distance) and used for calculation of site-level means of ∆C. To assess the error of the model, the root mean squared difference (RMSD) was calculated based on leave-one-out cross validation (Stage and Crookston 2007). The procedure was completed using R version 2.11 (R Development Core Team 2011) and yaImpute version 1.0-10 (Crookston and Finley 2008). Predictor variables were selected if correlations with inventory-based C stocks in 2002 and 2006 were significant at the 95% confidence level and there was no significant co-linearity with other predictor variables. Following the calculation, the 2002 – 2006 cumulative NEP footprint probability density surface for each site (calculated by Chen et al. 2008 and Chen et al. 2009) was used to weight each 10 m by 10 m footprint cell for the calculation of the site level mean value.

Site-level estimates of ΔC and ΔC₂ stocks and were calculated three ways. First, inventory- based site estimates ΔC and ΔC₂ stocks were calculated as an arithmetic mean of ground plots at each site. Second, ΔC₂ stocks were calculated with the K-MSN prediction. Third, and finally, ΔC₂ stocks were calculated with the K-MSN prediction weighted by flux footprint probability density.

Flux tower ∑NEP

Measurement and calculations of flux tower annual NEP published by Black et al. (2008) and Krishnan et al. (2009) were summed to find the cumulative ∑NEP (2004–2006) used for this study. These values were gap filled for conditions at night when friction velocity was low (i.e. inadequate turbulent mixing) and corrected for energy-balance closure. Annual NEP values for 2003, 2004, 2005, and 2006 were summed for each site to estimate the ΣNEP total for the four-year period. Morgenstern et al. (2004) reported that uncertainty in the annual NEP measured using EC at the near-rotation stand may be as much as 0.9 MgC ha⁻¹ year⁻¹ (due to systematic error in the EC measurements). This suggests that the uncertainty in the estimate at the near-rotation stand may be ±3.6 MgC ha⁻¹ over the full four-year measurement interval. Detailed estimates of the uncertainty that propagates through the calculation based on the friction velocity thresholds (e.g. following Richardson and Hollinger 2007) was beyond the scope of this paper.

CBM-CFS3 ∑NEP

CBM-CFS3 is a forest C accounting model used for national reporting of annual C inventories for Canada’s managed forests (Kurz et al. 2009). This model uses growth and yield equations and allometric equations to estimate tree growth and stand net primary production (NPP). The C mass of overstory trees and overstory tree roots is use to estimate litterfall and root mortality transfers to the aboveground and belowground detritus C pools, respectively. A soil submodel estimates decomposition and heterotrophic respiration (Rh) based on the size of various detrital pools and mean annual temperature. The model tracks all major C pools to ensure closure and estimates annual NEP from NPP - Rh. Wang et al. (2011) modeled forest processes using the CBM-CFS3 model for the 5 × 5 km area spanning the Oyster River and encompassing the DF49, HDF90, and HDF00 sites (Figure 1). Model runs for the area were performed on 1 hectare grid cells of forest disturbance history data, forest cover data, growth and yield equations, and disturbance transition matrices from Trofymow et al. (2008). The modeled ∑NEP values for the DF49, HDF00 sites (from Wang et al. 2011) were calculated as site-level means unweighted and weighted by the flux probability distribution for the model cells in the footprint area; and at HDF90 as an unweighted mean of model cells over the spatial extent of the ground plots. Annual model values of NEP for 2003, 2004, 2005 and 2006 were summed to calculate the ∑NEP and mean annual NEP for the four year period (Figure 4).

Figure 4

Flux tower footprint isolines and CBM-CFS3 model grids for a) the regenerating clearcut (HDF00) b) the near-rotation stand (DF49). For main species, Fd = Douglas-fir (Pseudotsuga menziesii), Hw = western hemlock (Tsuga heterophylla), Cw = western redcedar (Thuja plicata), Dr = red alder (Alnus rubra), Ba = Amabilis Fir (Abies amabilis). Site index is an estimation of height of typical dominant and co-dominant trees in even-aged and undisturbed sites at 50 years age to indicate productivity.

Comparing convergence of ΔC, EC tower ΣNEP, and CBM-CFS3 ΣNEP

Comparisons were made among the estimates by first evaluating how each method ranked the stand ages. Second, estimates of EC tower ΣNEP and CBM-CFS3 ΣNEP were compared. Third, at each stand age from youngest to oldest, the estimates of ΔC and tower and CBM-CFS3 ΣNEP were compared for convergence by comparing the means, variance around the means indicated by the standard deviation, and spatial variance indicated by the RMSD. Comparisons of convergence among methods were made for the following estimates: ΔC₂ as an arithmetic mean of plots (ΔC₂ unweighted), ΔC₂ with K-MSN classification (ΔC₂ K-MSN), ΔC₂ with K-MSN classification and footprint weighting (ΔC₂ K-MSN FPW), CBM-CFS3 ΣNEP as an arithmetic mean of cells within the tower footprint, CBM-CFS3 ΣNEP as a footprint weighted mean (ΣNEP CBM-CFS3 FPW), and EC tower ΣNEP.

Results

Live biomass ΔC_L

For all sites, there was a positive change in the live C mass of overstory trees (Table 2). For the near-rotation stand and juvenile stands, the increase was larger than at the regenerating site. The increase in live biomass C was nearly equal at the near-rotation and juvenile stands. An increase in shrub, herb, and bryophyte understory ΔC was observed at the near rotation and regenerating stands, while a small decrease in understory ΔC stocks occured at the juvenile stands, possibly related to an expected increase in canopy closure as the overstory trees matured. Comparing the two juvenile stands (HDF88 and HDF90), HDF88 had higher live biomass C than HDF90 (for example, live stem C and total live C were more than 1 standard deviation larger).

Table 2 Inventory-based site C stocks and ∆C estimates using arithmetic plot means (standard deviation)

Detritus ΔC_D1

Dead standing tree C decreased in the near-rotation stand due to dead standing trees falling during the measurement interval, and there was a small increase at the other sites due to tree mortality (Table 2). Large woody debris increased in the near-rotation stand, and decreased at the other sites. There was considerable variability among plots, indicated by the high standard deviations compared to the means. Fine woody debris decreased at all sites, except one of the juvenile sites (where there was a small increase). The large decrease in fine woody debris at the regenerating clearcut was likely due to decomposition of debris from the recent harvest. Finally, forest floor material increased at all sites, with the smallest increase at the regenerating clearcut, and the largest increase at the juvenile stands.

Mineral soil ΔC_s

Soil C (0–15 cm) stocks was highest at the near-rotation site, followed by the two juvenile stands, and lowest at the regenerating clearcut (Table 2). Over the measurement interval, large increases in mineral soil (0–15 cm) ΔC_s were measured at all sites, with the highest at one of the juvenile stands (HDF88).

Detritus annual accounting method ∆C_D2

Measured litterfall was highest at the near-rotation stand, followed by the juvenile stands, and the regenerating stand had the lowest amount (Figure 5). The majority of litterfall at the near rotation stand and juvenile stands was needles followed by twigs. At the near-rotation stand, litterfall was highest in 2001–2002, decreased through 2004–2005, and slightly increased in 2005–2006. The other sites were relatively constant through time.

Figure 5

Measured annual litterfall 2001–2006. Needles include green and senescent foliage, cones include male and female cones, twigs include woody material such as twigs or small branches, and other includes miscellaneous pieces such as leaves, mosses, and lichens.

Results from the annual accounting method for detritus, showed a net gain for the sum of components at the near-rotation stand, and a decrease at the other sites (Table 3). The increase in C at the near-rotation stand was slightly larger using ∆C_D2 than the direct measurement method; however, the biggest difference was at the juvenile stand, which showed a net loss using ∆C_D2, and a large positive balance using the direct measurement method. For stumps, this method indicated that decomposition may have been greater than was initially expected (stumps were not re-measured due to an expectation of minimal decomposition), especially at the more recently disturbed juvenile and regenerating clearcut sites. Large woody debris increased at the near-rotation stand, and decreased at the other sites, due to less dead standing trees to fall and less branchfall. The increase in fine woody debris was estimated to be much lower using ∆C_D2, with the near-rotation stand having a net-accumulation of fine woody debris, the juvenile sites had a small loss, and the regeneration site experienced the largest decrease. This indicates that, given the available inputs to litterfall and expected rates of decomposition for fine woody debris, the direct measured stock changes in fine woody debris was very likely unrealistically high. Finally, for the forest floor layer, the near-rotation stand was estimated to have a large increase, while the other sites decreased, with the greatest decrease at the regenerating clearcut.

Table 3 Inventory-based site detrital C in 2002 and ∆C _D2 (2002–2006) using plot mean (standard deviations) litterfall, mortality and decomposition values

Since the plot values for ∆C_D2 were judged more reliable, all subsequent analyses to spatially scale the plot data to the site were made using ∆C₂.

Spatial distribution of C stocks and ∆C₂ stocks

Predictor variables (Table 4) were selected and used for the K-MSN procedure at the sites with flux towers (Figure 6). Slope was an important predictor for all sites as it was strongly correlated with soil C stocks (with less sloping plots having higher soil (0 – 15 cm) C stocks), and at the regenerating clearcut, slope was also correlated with living tree C stocks. At the near rotation site, forest-inventory mapping was a significant predictor variable for overstory C stocks. Due to the smaller footprints at the other sites, forest-inventory mapping varied little within the site; however, at all sites the other spatial predictor variables showed finer scale variation. At the regenerating clearcut, NDVI was a significant predictor for woody debris, forest floor, and mineral soil C stocks.

Table 4 Environmental predictor variables used to determine stratification units at DF49, HDF88, and HDF00

Figure 6

First most similar neighbour (MSN) sample plots demonstrate the patterns of spatial variability for A) the near rotation stand, B) juvenile stand, C) recent clearcut. Background imagery is pan-sharpened 2004 Quickbird for A and C, and true-colour orthophoto for B.

The Mahalanobis distance (Figure 7) demonstrates how comprehensively the inventory ground plots represent stand conditions across the footprint. Lower percentiles represent small distances that were relatively well represented by the ground plots. Areas with the smallest Mahalanobis distances and best ground plot representation were close to the towers in the areas contributing most to NEP measurements (within the 50% cumulative flux probability density isoline). Notable areas with large Mahalanobis distances included areas with different forest cover types (e.g., patches of hardwoods at the near-rotation stand), very few overstory trees (at the juvenile stand), and short steep slopes (between bench terraces at the regenerating clearcut).

Figure 7

Total Mahalanobis distance for K-MSN stratification units for A) the near rotation stand, B) juvenile stand, C) recent clearcut. Mahalanobis distance percentiles are for each site. Footprint cells with green tones are below the median Mahalanobis distance at site and are relatively well represented by inventory-based sample plots, while footprint cells with red tones are above the median Mahalanobis distance for each site and are less well represented. Background imagery is pan-sharpened true-colour 2004 Quickbird for A and C, and true-colour orthophoto for B.

For site level estimations of ΔC_L at the regenerating clearcut, K-MSN was more than 1 standard deviation lower than the unweighted mean, and K-MSN-FPW was 1 standard deviation lower than K-MSN. For ΔC_D2, K-MSN and K-MSN FPW were less than 1 standard deviation higher than the unweighted mean (Tables 3 and 5). The totals for ΔC₂ K-MSN and K-MSN FPW were less than 1 standard deviation higher than the arithmetic mean (Table 6).

Table 5 Site mean (± RMSD) changes in live biomass (∆C _L) and detrital (∆C _D2 ) C stocks calculated using three most similar neighbours classification (K-MSN), and K-MSN with footprint weighting means (K-MSN FPW) of plots for sites with flux-towers

Table 6 Comparisons of inventory-based site C stock changes (∆C ₂ ) determined from plot means (live biomass and detrital fluxes), K-MSN scaling mean (standard deviation), and K-MSN with footprint weighting (K-MSN FPW) mean (standard deviation) to estimates of ∑NEP from flux towers and ∑NEP from CBM-CFS3 for the (a) 4-year period 2003–2006 Mg C ha ⁻¹ and the (b) mean annual value for each method in Mg C ha ⁻¹ yr ⁻¹

For the juvenile stands, K-MSN and K-MSN FPW ΔC_L was less than 1 standard deviation higher than the unweighted mean. The K-MSN estimate was higher than the K-MSN FPW estimate. For ΔC_D2, K-MSN and K-MSN FPW estimates were slightly higher than the mean (less than 1 standard deviation). For the total ΔC₂, the estimates were less than 1 standard deviation different, and the standard deviations were very large compared to the means, indicating a large amount of variability at the site.

At the near-rotation stand, ΔC_L K-MSN and K-MSN FPW were slightly higher than the arithmetic mean (less than 1 standard deviation). The estimates using KMSN and KMSN-FPW were identical. For ΔC_D2, the unweighted mean was higher than K-MSN and K-MSN FPW, but less than 1 standard deviation different. For ΔC₂ total, all estimates were within 1 standard deviation.

Comparing the ΔC₂ and ΣNEP methods at the juvenile stands, HDF88 was compared with measurements from the flux tower installed at the site, and HDF90 was compared with CBM-CFS3 model output. While the two stands were similar in seral stage, species composition, and stand history, there were differences in stand composition. For example, HDF88 had higher live C stocks, in particular, live stem C, understory stem C, and total stem C. In addition, HDF88 had larger detritus C stocks including stumps C, large woody debris C, and total detritus C. The stand variability at HDF88 was also greater than at HDF90 indicated by the large standard deviations compared to the means. Given the large amount of spatial variability at HDF88, spatial scaling, and footprint weighting had a notable effect on the amount and sign at this site.

Comparing inventory ∆C₂, EC tower ΣNEP and CBM-CFS3 ΣNEP

All methods ranked the near-rotation stand as a sink for C (net C uptake), the juvenile stands as a weak sink or weak source (small C uptake or release), and the regenerating clearcut as a source (net C release). While the ranking was consistent for all methods, there were differences in the magnitude of values for each stand (Table 6).

At the regenerating clearcut, the CBM-CFS3 estimate of ΣNEP with footprint weighting converged more closely with the EC tower estimate of ΣNEP than the unweighted estimate, and both were within 1 standard deviation. Evaluating the yearly means (Figure 8a), demonstrated that footprint weighting improved convergence between the CBM-CFS3 estimate of ΣNEP and the EC flux ΣNEP. This was primarily due to less weighting for patches of uncut forest in the periphery of the flux footprint (Figure 4).

Figure 8

Yearly NEP estimates for CBM-CFS3 and EC flux tower (Black et al. 2008) for (A) the regenerating clearcut and (B) the near-rotation stand. CBM-CFS3 is shown as yearly mean of footprint cells with equal weighting (Fp equal) and weighting by NEP flux probability density distribution (Fp weight).

At the juvenile sites, CBM-CFS3 ΣNEP indicated that HDF90 was a weak sink (net C uptake) and EC tower ΣNEP indicated that HDF88 was a weak source (net C release) (Tables 3 and 6). At the near rotation stand, the estimates of ΣNEP from CBM-CFS3 were higher than the EC flux tower ΣNEP, with the footprint weighted estimate being the highest. Considering the means of ΣNEP from CBM-CFS3 on a yearly basis, footprint weighting decreased convergence with the EC flux tower measurement of ΣNEP (Figure 8b).

At the regenerating clearcut, unweighted ΔC₂ was within 1 standard deviation of the CBM-CFS3 estimates of ΣNEP, and closest to unweighted CBM-CFS3 ΣNEP. The unweighted mean of ΔC₂, ΔC₂ K-MSN, and ΔC₂ K-MSN FPW were greater than 1 standard deviation higher than the EC flux tower estimate of ΣNEP. ΔC₂ K-MSN was within 1 standard deviation of the CBM-CFS3 unweighted ΣNEP. Both ΔC₂ K-MSN and ΔC₂ K-MSN FPW were greater than 1 standard deviation higher than ΣNEP CBM-CFS 3 FPW (Table 6).

At the juvenile stands, all estimates ∆C₂ were within 1 standard deviation of the estimates of ΣNEP. At HDF90, the unweighted estimate of ∆C₂ matched ΣNEP CBM-CFS3. At HDF88, the unweighted mean of ground plots was the closest to the EC tower measurements of ΣNEP, while the ΔC₂ K-MSN had a different sign (indicating a small source), and the K-MSN ∆C₂ FPW had the same sign as ∆C₂ unweighted (Table 6).

At the near-rotation stand, all estimates of ∆C₂ were within 1 standard deviation of all estimates of ΣNEP. The estimates of ∆C₂ (and CBM-CFS3 ΣNEP) higher than the EC flux tower measurements of ΣNEP (Table 6).

Discussion

Several approaches can be used to compare forest inventory data with EC ΣNEP. For example, several authors have used measurements or estimates of soil and detritus respiration in combination with inventory measurements of overstory productivity (Law et al. 2001; Ehman et al. 2002; Kolari et al. 2004; Black et al. 2005). In the present study, comparisons were based on measurements of ∆C which were higher than EC ΣNEP. Black et al. (2005) also found a divergence between inventory measurements of ∆C (including ∆C_S measurements) and EC ΣNEP, with ∆C higher than EC ΣNEP; however, they also calculated ΣNEP using measurements of heterotrophic respiration combined with biometric estimates of NPP, and found that the bias was not observed. They concluded that ∆C systematically overestimated NEP due to unaccounted decomposition processes and uncertainties in ∆C_S. At DF49, Jassal et al. (2010) measured heterotrophic soil respiration in 2007, which was a major portion of total soil respiration. Including measurements of soil heterotrophic respiration with net primary production estimates from forest inventory may improve convergence with EC-flux tower measurements. However, a limitation of taking measurements of soil respiration is that the measurements require long-term collection using special equipment and are not typically collected in traditional forest inventory.

Kolari et al. (2004) found that management practices, such as clearcut harvesting had a large impact on the soil C balance, with clearcut sites soil C sources, while non-disturbed sites had more stable soil C balances. In this study, a trend in soil C stocks, with the highest soil C stocks at the oldest site, and the lowest at the recently disturbed site, was consistent with the trend by seral stage reported by Kolari et al. (2004). For the two similar-aged juvenile stands, since there was no flux tower installed at HDF90 and CBM-CFS3 model runs results were not available at HDF88, a direct comparison was not possible. In addition to the differences in measurement methods, the differences in ΣNEP between the two juvenile stands could also be partially attributed to site conditions. Notably, HDF88 had more detrital C stocks, likely leading to larger releases of C due to decomposition. This was reflected in the more negative ΔC_D2 at HDF88 compared to HDF90.

For the forest inventory measurement of ∆C and its components ∆C_L, ∆C_D, and ∆C_S, differences in procedure, measurement error by field and lab crews, and sample design may introduce error into inventory measurements that may be larger than the magnitude of change in C stocks over short measurement intervals, such as the 4 year period used in this study. Overall, ΔC_D1 and ΔC_S were much larger than reported in previous studies (e.g. Law et al. 2001), and was much larger than expected given the measured litterfall inputs and estimated decomposition rates. The NFI (2008) calls for 10-year remeasurement intervals, which may be appropriate for trees and detritus, but may be too short for ΔC_S. In addition, direct use of forest inventory data relies on subjective classification of decay class that can differ from one-field-crew to the next, so utilizing decay and transfer algorithms may reduce the potential effect of these classifications. Developing site-specific allometric equations is labour intensive, destructive, and costly, and therefore uncommon in most studies; however, use of existing relationships can lead to errors due to differences in tree architecture, and wood density. Further, the use of inappropriate allometric equations can be a significant source of error in forest productivity studies (Clark and Brown 2001). The NFI data compilation procedure provided regionally applicable allometric equations; however, the errors in this large pool of equations have not yet been estimated and therefore is a limitation.

The measured changes in ΔC_S were much larger than was expected for a four year period. For example, Law et al. (2001) estimated soil sequestration of 0.7 – 0.8 MgC ha⁻¹ year⁻¹ in a Douglas-fir stand and Gielen et al. (2013) found no significant change over an eight-year interval in a Scots pine forest. In other studies by Ehman et al. (2002), Kolari et al. (2004), Miller et al. (2004), Ohtsuka et al. (2007), and Granier et al. (2008) ΔC_S was assumed to be zero over the measurement interval. The large value for ΔC_S in this study may be explained by several factors related to measurement methodology. First, mineral soil bulk densities in the 0–15 cm layer were assumed not to change and thus the BD values measured in 2002 were also used in 2006. If the bulk density decreased then the C_S change would be overestimated. Second, at the regenerating clearcut and juvenile stands, an increase in fine root mass from 2002 to 2006 included in the <2 mm soil fraction could also have accounted for some of the increase in soil C. Third, samples in 2002 were taken by 10–12 cm diameter hole excavations, while samples in 2006 were collected using a 2 cm diameter soil corer, which, due to compression at the opening of the sample corer, may be biased to soil from shallower depths where C content is likely higher. Therefore, in this study, measurements of ∆C_S were deemed unreliable, due to the large positive changes in ∆C_S, and assumed to be negligible over the four year period and thus not used for subsequent calculations of site level ∆C and further comparisons. In future work, efforts to reduce any chances of variation in the way the samples are collected and analysed may result in more reliable measurements of C stock changes. For example, collecting and processing the samples using identical methods at nearly the same locations, or collecting a much larger number of samples to capture a wider range of spatial variation could reduce the effect of spatial variability on measurements of ∆C_S.

Several authors have applied footprint analysis in an effort to correct for non-homogenous environments, and horizontal advection, for example due to daily upslope or downslope winds, which may introduce bias into the eddy-covariance measurements (Baldocchi 2008). For example, Ehman et al. (2002) used footprints to calculate an estimate of sensor bias given using a vegetation index, and this gave confidence that the sensor measurements were representative of the inventory measurements. Gielen et al. (2013) used a more complex footprint model applied to the processing of NEP data, were fluxes originating outside of the area where inventory measurements were sampled were removed from the NEP estimates. In this study, footprint weighting by ecological attributes and footprint probability density was evaluated where it facilitated comparison between inventory-based measurements of ecosystem ΔC stocks and tower ΣNEP by accounting for fine scale spatial variability in forest structure. The effect of K-MSN classification and footprint weighting had an effect on the estimated site-level values; however, this effect was smaller than differences due to measurement methodology. Notably, applying footprint weighting to the estimation of CBM-CFS3 ΣNEP at the regenerating clearcut increased convergence with EC tower ΣNEP; however, footprint weighting did not have the same effect for ΔC₂. The largest differences observed among methods were seen at the regenerating clearcut, where the inventory approach indicated the site was a much lower source than the ΣNEP methods. A possible cause is that the inventory methods did not properly account for C losses from decomposition of coarse roots from stumps and large woody debris at recently harvested sites, which can be a substantial source of respiration (Janisch et al. 2005). Therefore, if an estimate of stump coarse root decomposition were included in ΔCD₂ it may have increased convergence with the ΣNEP estimates. Additionally, the original sample locations were designed to capture the range of conditions in near proximity of the tower (based on interpretation of aerial photographs), and the relationship was similar when scaled up to the site level. Future studies may seek to derive accurate footprint estimates prior to establishing sample plots to ensure they are representative of the broader site conditions. Another more costly alternative is to establish a much greater number of plots in a systematic basis around the tower, to ensure the spatial variation within the tower footprint is adequately captured.

This research highlights how different measurements and approaches for C accounting include different constraints, and as a result, computations using different data sources may not converge. Each type of measurement has “different strengths and weaknesses, but the combination of multiple measurements and modeling has the potential for refining estimates of C stocks and fluxes” (Law et al. 2004). For example, eddy-covariance data are temporally detailed, and provide information about daily and seasonal processes. However, the data are less spatially extensive and detailed than the forest inventory data, for example, only three sites were available in this study each covering an area with diverse vegetation within the tower footprint. Measurement error was estimated to be relatively low, on the order of 0.9 Mg C ha⁻¹ year ⁻¹ (Morgenstern et al. 2004). In contrast, forest inventory data are more spatially extensive and include detailed measurements of forest attributes. The main constraint for forest inventory data is that they are less temporally detailed (for example, measurements were available at only two points in time), and the error in these measurements may be greater than the expected rate of change over the measurement interval. The inventory plots in this study sampled a broad range of vegetation conditions and provided insight into the representativeness of flux-tower footprints. Estimates of error of 1.2 to 1.3 Mg C ha⁻¹ were identified relating to spatial representativeness of the inventory plots within the tower footprints. Finally, C budget models are important for understanding landscape scale C processes. Since these models depend on numerous coefficients, algorithms and assumptions, it is valuable to compare model outputs with EC flux and inventory measurements to assess their performance and evaluate their behaviour. Considering all of these approaches together, the relative strengths and limitation of each approach can be taken into account in evaluating the suitability of each type of measurement at varying spatial and temporal scales.

Conclusions

The comparison of ΔC from inventory measurements with ΣNEP from CBM-CFS3 and EC flux-tower measurements demonstrated an agreement among the methods in trends across stand seral stage; however, due to differences in the measurement approaches, there was some divergence in the results. Convergence amongst methods was closest for the juvenile sites (HDF90 and HDF88), which were in transition from C source- to sink. At the most recently harvested site (HDF00), the EC flux ΣNEP and CBM CFS3 ΣNEP indicated that the site was a greater source of C over the time period than ∆C from inventory methods. At the near-rotation site the inventory ∆C₂ method, unweighted CBM CFS3 ΣNEP and EC flux ΣNEP also converged.

Each of the measurement methods had advantages and limitations. For example, while EC flux towers provide temporally rich data, they had limited spatial coverage. Obtaining ΔC using forest inventory measurements included several challenges such as the consistency of data collected over long time intervals, small changes in ΔC compared to uncertainty in the measurements, and insufficient measurement of processes such as decomposition of all detrital pools. Forest inventory data are available over broad areas; therefore methods that combine measurements with estimates of decomposition are needed to provide information over a regional scale. Footprint analysis using the spatial predictors from remote sensing and topography can give confidence that sample locations represent broader area site conditions. Considering the advantages and constraints of each type of measurement is helpful in choosing appropriate measurement methods at varying spatial and temporal scales.