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

Colin J Ferster; JA (Tony) Trofymow; Nicholas C Coops; Baozhang Chen; Thomas Andrew Black

doi:10.1186/s40663-015-0038-3

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

Ferster, Colin J;Trofymow, JA (Tony);Coops, Nicholas C;Chen, Baozhang;Black, Thomas Andrew 2015-12-01 00:00:00 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. Keywords: Forest carbon; Micrometeorology; Biometry; Remote sensing; Geographic information systems * Correspondence: colin.ferster@ubc.ca Department of Forest Resources Management, University of British Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada Full list of author information is available at the end of the article © 2015 Ferster et al.; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. Ferster et al. Forest Ecosystems (2015) 2:13 Page 2 of 19 Background 2008), and by measuring forest C stocks at two times Forests are a large component of the global carbon (C) with a sufficient interval between measurements, the C stocks, containing an estimated 1146 Pg C (Dixon et al. stock change (ΔC) can be determined and used to esti- 1994). Forest processes, which may be influenced by forest mate the cumulative net uptake of C from the atmos- management, can therefore have a large impact on the phere to the forest for the period of measurement global C budget, either by storing or releasing C. It is (ΣNEP) (Clark and Brown 2001). Measurements are typ- therefore critical to understand forest C dynamics, includ- ically made of individual plants, for example, all trees ing forest C stock components and transfer mechanisms within a plot are measured for diameter, species, and in order to develop accurate forest C models such as the height, then allometric relationships are used to deter- Carbon Budget Model of the Canadian Forest Sector mine tree mass (Ter-Mikaelian 1997). Similarly, woody (CBM-CFS3) (Kurz et al. 2009), 3PG (Landsberg and debris pieces are measured along transects using line Waring 1997), and Ecosys (Grant et al. 2007), amongst intersect sampling and allometric relationships are ap- others, as well as to inform forest management policy and plied to find biomass (Brown 1971). Other components for national and international reporting (Kurz et al. 2002). are made by direct measurements, for example, under- Factors affecting forest C dynamics include natural (e.g. story vegetation can be sampled, dried, and weighed fire, insect outbreaks) and anthropogenic disturbances (Bailey et al. 1998). To measure soil C, field samples are from land use management and change (e.g. harvest, re- typically dried and weighed for calculation of bulk density, forestation, deforestation) (Kurz et al. 2002) as well as and sent for laboratory analysis of C using dry combustion weather (Morgenstern et al. 2004), fertilization (Jassal (Janzen 2005). Soil C measurements are important be- et al. 2010), and stand age and species composition, which cause the underground processes driving soil respiration can be related to disturbance history, site edaphic char- form a large component of ecosystem gas exchanges, but acteristics, or management practices (e.g. silvicuture) these are less well-understood and measured compared to (Humphreys et al. 2006; Krishnan et al. 2009). above ground stocks and processes (Ehman et al. 2002). To measure the net exchange of C between land eco- Changes in soil C storage over time have implications as a systems and the atmosphere (net ecosystem exchange, source or sink for the global C cycle, and as an indicator NEE, with -NEE referred to as net ecosystem productiv- of environmental function and health (Janzen 2005). ity, NEP), a global network of over 400 eddy-covariance One application of forest measurement data to provide (EC) flux stations has been established across a range of accounting of forest management actions and subse- ecosystems, building an extensive data record of NEP, in quently inform forest management decisions is as an in- some cases spanning up to two decades. The majority of put to C budget models such as CBM-CFS3. CBM-CFS3 these towers are located on sites not undergoing major is a forest C budget model that utilizes growth and yield disturbances so the fluxes measured reflect the inter- information for biomass and generates explicit simula- action of weather, vegetation composition, stand age, tion of dead organic matter dynamics (Kurz et al. 2009). and seasonal phenology. EC flux-towers use micro- In contrast, forest process models use approaches based meteorological equipment to take measurements at the on the understanding of processes such as photosyn- canopy scale of the exchanges of gasses including CO , thesis; these process-based models have more potential water vapor, sensible heat, and some are equipped to to model changing conditions because they do not rely measure other trace gases (Baldocchi 2008). The source on historic growth and yield and therefore are better area contributing to measurements made by the instru- suited to simulate forest conditions under global change ments mounted on the EC flux-tower, the flux-tower foot- situations (Landsberg and Waring 1997). CBM-CFS3 print, is variable in size and shape depending on height of was designed to meet reporting requirements for forest measurement, surface roughness length, wind speed and management, provide policy support, and a function as a C budgeting tool for operating foresters, thus function- direction, and atmospheric stability (Leclerc and Thurtell 1990; Schmid 2002), and proper interpretation of EC flux- ing at a range of spatial scales, and estimating records of tower based measurements depends on the flux-footprint C stocks, transfers between pools, and emissions (Kurz et al. 2002; Kurz et al. 2009). over which the fluxes are sampled (Chen et al. 2008). Early estimates modelled flux-footprints as simple ovals. Re- Comparing EC flux-tower cumulative NEP (ΣNEP) cently, estimates of flux-footprint climatology may dem- and inventory ground plot measurements of C stock changes (ΔC) over the same period can help inform for- onstrate more complex geometries and continuous probability density surfaces to quantify the upwind distri- est C processes for several reasons. First, ground plot bution of weighting factors over long time periods (Chen measurements can serve as an independent validation of ΣNEP measurements made at EC flux-towers. Second, et al. 2009). Forest measurements can be taken to quantify forest C inventory plot data can provide more detailed informa- stocks (Dixon et al. 1994; Clark and Brown 2001; NFI tion about the stand structural changes that may be Ferster et al. Forest Ecosystems (2015) 2:13 Page 3 of 19 driving variation stand-level EC C exchange. Third, broad determined and the ΔC stocks compared to CBM-CFS3 stand inventory measurement datasets are available for a modelled and EC-flux-tower estimates of ΣNEP at the greater spatial extent than the global EC flux-tower net- CCP coastal British Columbia Flux Station. Through a work. For example, the Canadian National Forest Inven- combination of variable topography and a complex dis- tory (NFI) samples approximately 22,000 photo plots turbance history (Trofymow et al. 2008), the EC flux- locations across Canada with detailed ground plots mea- tower sites possess fine-scale spatial variation in forest sured at over 1000 locations (Gillis et al. 2005). However, structure, thus presenting a challenge to the interpret- comparing inventory plot changes and EC flux-tower ation of EC flux data (Schmid and Lloyd 1999; Göckede measurements requires comparisons to be made across et al. 2004). To allow the comparison between measure- different spatial and temporal scales. Therefore, such com- ments from the inventory plots, C budget models, and parisons pose methodological challenges beyond individ- EC flux-towers, our approach consisted of four steps: ual inventory plot and EC flux-tower measurements. first, we found the change in C based on forest inventory Comparisons between EC tower measurements of ground plot measurements made in 2002 and 2006 (four ΣNEP and inventory plot measurements of ΔC stocks year period); second, we developed an additional ap- have been completed at a number of sites globally proach that utilized litterfall data and published decay (Schulze et al. 2000; Granier et al. 2000; Law et al. 2001; equations, for a second inventory-based estimate of C Barford et al. 2001; Ehman et al. 2002; Curtis et al. 2002; stocks change; third, we defined stand attributes at a fine Miller et al. 2004; Black et al. 2005; Gough et al. 2008; spatial resolution, stratified the site based on stand attri- Kominami et al. 2008; Yashiro et al. 2010; Gielen et al. butes, and weighted the two inventory plot estimates; 2013; Babst et al. 2014), and at two sites of the Canadian and fourth, calculated ΔC. These were compared Carbon Program (CCP), Boreal Ecosystems Research and with CBM-CFS3 model-based estimates of ΣNEP and Monitoring Sites (BERMS) station, located in the boreal EC flux-tower-measured ΣNEP for four years (2003, forest of northern Saskatchewan (Theede 2007). Over a 2004, 2005 and 2006) in the same period. Finally, we ten-year interval (1994 – 2004), Theede (2007) found a evaluated how the forest age and conditions at the dif- −1 −1 convergence of 15.6 ± 4.0 MgC ha 10 years (inventory ferent sites accounted for convergence or lack of conver- −1 −1 plots) and 18.2 ± 8.09 MgC ha 10 years (EC tower) at gence in estimates for the different methods, and −1 −1 the Old Aspen site and 5.8 ± 2.0 MgC ha 10 years (in- discussed the relative advantages and constraints of each −1 −1 ventory plots) and 6.9 ± 1.6 MgC ha 10 years (EC measurement method. tower) at the Old Jack Pine site. However, the homoge- neous stand structure contributing to canopy-level EC Methods flux measurements and large fetch in the boreal forest re- Study area duced the need for consideration of spatial vegetation The coastal British Columbia Flux Station sites are located structure and footprint distributions. on the leeward central east coast of Vancouver Island These previous studies have contributed considerably to (Figure 1), in the Coastal Western Hemlock biogeocli- the understanding of forest C dynamics; however, at sites matic zone (Green and Klinka 1994) with a mean annual with more complex vegetation structure and topography, air temperature of 10° Celsius. Douglas-fir (Pseudotsuga such as the forests found in coastal British Columbia, the menziesii var menziesii (Mirb.) Franco) is the dominant spatial distribution of vegetation structure and the foot- species, with lesser amounts of western hemlock (Tsuga print distribution are important considerations for accur- heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata ate comparisons of EC flux-tower and inventory plot Donn ex D. Donn), and red alder (Alnus rubra Bong.). measurements (Schmid and Lloyd 1999). More complex Three eddy-covariance and meteorological towers were ecosystems, therefore, require more complex models and placed in stands of different ages (Table 1). methods developed for accurate comparisons of EC flux- EC flux-tower footprints, which delineate the spatial tower measurements and inventory-based measurements. distribution of upwind source-weighting factors, were For example, a limited number of studies, including calculated by Chen et al. (2009) to represent the size, Ehman et al. (2002), Chen et al. (2009), Ferster et al. shape, direction, and magnitude of the flux as a function (2011), and Gielen et al. (2013) applied weighting by EC of wind speed and direction, measurement height, and flux-tower footprint climatology to inventory-based mea- surface roughness. Chen et al. (2009) calculated hourly surements, where areas of the footprint that contribute footprints for 10 m by 10 m cells, weighted by NEP, and more to tower-measured flux are weighted more heavily averaged over the 2002 – 2006 measurement period to than other areas, resulting in an overall improvement determine the footprint climatology for the sites. In this when compared to simple flux footprint geometries. study, the flux-footprint climatology 85% cumulative flux In this paper, the difference between 2002 and 2006 probability boundary was selected as the maximum ex- inventory plot measurements of ecosystem C stocks was tent of analysis, since the majority of the probability Ferster et al. Forest Ecosystems (2015) 2:13 Page 4 of 19 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. density is concentrated within this area and the remain- Live biomass der extends across a large area (Figure 2). Live trees greater than 1.3 m tall were measured in an 11.28 m radius circular plot. Where trees were very nu- merous (for example, at HDF88), a half plot or 5 m ra- Inventory ground plots dius circular plot was measured following NFI guidelines All ground plots were established and measured by the (NFI 2008). Allometric equations are used in the NFI Canadian Forest Service at the four sites following compiler to estimate stem, bark, branch, and foliar bio- Canadian National Forest Inventory (NFI) guidelines mass based on work by Lambert et al. (2005). These bio- and protocols (NFI 2008) though they are not part of mass values were converted to C content by assuming the primary NFI ground plot network. Plot locations 52, 56, 52, and 52% dry C concentration, respectively were chosen to represent each site series class (based (Matthews 1993). on air photo interpretation and field transects) within Understory vegetation from four 1 m × 1 m micro- the preliminary flux footprint boundary estimated in plots was destructively sampled and sent to the labora- 2002 (Figure 2). Establishment and initial measurements tory for drying and weighing for determination of mass were made in September 2002 and re-measurements were which was converted to C content assuming 50% dry C made in September 2006, constituting a four-year meas- concentration. urement interval. NFI data-compilation software used to determine plot-level values through application of allo- metric equations to overstory trees and density values for Detritus woody debris NFI (2008). Since live overstory vegetation Woody debris, fallen woody material > 1 cm diameter, roots were not measured in 2002 or 2006, coarse and fine was measured for diameter, species, and decay class along root C mass was estimated from overstory biomass using four 15 m transects at each plot. Decay class was deter- equations by Li et al. (2003). mined by field crews in the five-class ordinal rating system Ferster et al. Forest Ecosystems (2015) 2:13 Page 5 of 19 Table 1 Study site locations and characteristics mass. C content was determined assuming 50% dry C Site Location (latitude Site description concentration of woody debris mass. Fine woody debris, and longitude) fallen woody material 2 mm to <1 cm diameter, was sam- DF49 49.868797°, � Near-rotation stand established pled from within the four 1 × 1 m understory microplots in 1949. in each sample plot and sent for laboratory drying, weigh- (Figure 2a) - 125.33515° � 162 ha flux tower footprint. ing, and determination of mass, and C content assuming � 12 ground plots. dry mass is composed of 50% C. Dead standing trees were measured for diameter, height, and species and decay class � 260 to 470 m elevation. was observed and recorded by field crews. Stumps were � Harvested of old growth timber measured for diameter and decay class in 2002 but were in 1937, 1938, and 1943. not remeasured in 2006. � Broadcast burned in 1938 and 1943 following harvest. Surface substrates were measured along four 15 m long transects to determine the average depth and per- � Planted in 1949. cent area surface coverage in the entire plot. The organic HDF88 49.536655°, � Juvenile pole-sapling stand litter, fibric, and humus layer (forest floor) was sampled established in 1988. and average depths measured from within four 20 × 20 (Figure 2b) - 124.90146° � 35 ha flux tower footprint. cm templates (one at each microplot). These destructive � 6 ground plots. samples were collected, sieved, dried, and weighed to de- � 150 to 220 m elevation. termine bulk density and subsamples sent for laboratory � Harvested of second-growth analysis of C concentration. timber in 1987. � Broadcast burned following harvest. Mineral soil � Planted in 1988. In 2002, <2 mm mineral soil was sampled for bulk dens- HDF90 49.893666°, � Juvenile pole-sapling stand ity and C concentrations from 10 – 12 cm diameter ex- established in 1990. cavated holes at three depth intervals 0–15 cm (4 - 125.304415° � No flux tower. samples), 15–35 cm (2 samples), and 35–55 cm (one � 6 ground plots. sample) and % volumetric coarse fragments determined from one 55-cm-deep soil pit per ground plot. Data were � 175 m elevation. scaled to the entire groundplot discounting for the area � Harvested of second-growth timber in 1990. without mineral soil (i.e. exposed bedrock). In 2006, the 0–15 cm layer was re-measured for C content at four lo- � Broadcast burned following harvest. cations in each plot using a 2 cm diameter soil corer. Samples were stored in plastic bags at 2°C prior to la- � Planted in 1990. boratory processing. For a detailed description of soil � Similar in age and composition to HDF88. sampling methods see NFI (2008). For 2002–2006 soil C stock changes, only the 0–15 cm layer was considered, HDF00 49.872177°, � Regenerating clearcut stand established in 2000. soil C at deeper depths was assumed unchanged. (Figure 2c) - 125.29235° � 14 ha flux tower footprint. � 9 ground plots. Changes in inventory ground plot C stocks � 160 to 190 m elevation. Changes in C stocks over the four-year interval were cal- � Harvested from a second-growth culated from the difference between C stocks in 2002 stand in 2000. and 2006. The total change in C stocks was calculated � Logging debris were piled and as: burned. � Planted in 2000. ΔC ¼ ΔC þ ΔC þ ΔC L D S where ΔC is the change in total live biomass C stocks in- utilized by the NFI, from intact (Class 1) to highly decom- cluding live trees, and shrubs, herbs, and bryophytes; ΔC posed (Class 5) (NFI 2008). NFI compilation routines is the total change in detrital C stocks including dead utilize algorithms for line intercept sampling to determine standing trees, woody debris, and the forest floor; and plot level volume and a lookup table of woody debris ΔC is the change in mineral soil C stocks (0–15 cm). density values by species and decay class to determine S Ferster et al. Forest Ecosystems (2015) 2:13 Page 6 of 19 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. ΔC Detritus annual accounting method between measurement dates, and assumed to turn over D2 Preliminary examination of the ΔC values showed at a rate of 80% annually. Shrubby litterfall material larger-than-expected changes in the measured C stocks was divided into stem and branch and foliage compo- over the four-year interval, especially for the forest nents, and assumed to turn over at a rate of 3% and floor component. To evaluate these measured changes, 95%, respectively similar to the fine branch wood and ΔC was estimated using a second method (referred to foliage for deciduous trees (Kurz et al. 2009). When as ΔC ) that accounted for the annual inputs into the trees died within the measurement interval, the bio- D2 detrital pools from litterfall and mortality (i.e. the trans- mass was transferred to the detrital pools half way fers from live biomass to forest floor and soil , and live through the measurement interval. Finally, the biomass trees to dead trees respectively) and annual losses and of any trees (live or dead) that were measured in 2002 transfers from decomposition estimated using pub- and not recorded by the field crew in 2006 was trans- lished equations and parameters (Smyth et al. 2010; ferred to the detrital pools half-way through the meas- Kurz et al. 2009). urement interval under the assumption that the tree Inputs to detrital C pools from 2002 – 2005 included fell during the measurement interval. overstory fine litterfall (needles, leaves, cones, twigs) Losses through decomposition and transfer to the soil which was collected quarterly in 3 0.189 m conical- organic C of each detrital pool was calculated at an an- mesh litterfall traps located in each ground plot at nual time step with an annual average temperature of DF49, HDF88 and HDF90. Annual litter fall masses 10°C using coefficients that were developed and vali- were calculated and converted to C assuming 50% dry dated nationally (Smyth et al. 2010) for the CBM-CFS3 C content. Herbaceous material, which was not tall model (Kurz et al. 2009). The transfers and coefficients enough to be sampled using litterfall traps especially at used are presented in Figure 3. Soil C stock changes the HDF00 site, was collected and measured at the were assumed negligible, and ΔC was calculated as: sample plots in 2002 and 2006, linearly interpolated ΔC ΔC + ΔC 2= L D2. Ferster et al. Forest Ecosystems (2015) 2:13 Page 7 of 19 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 using R version 2.11 (R Development Core Team 2011) To account for spatial heterogeneity due to stand struc- and yaImpute version 1.0-10 (Crookston and Finley ture and topography, the methodology developed by 2008). Predictor variables were selected if correlations Ferster et al. (2011) was applied to the ground plot data at with inventory-based C stocks in 2002 and 2006 were the flux tower sites. Following this approach, predictor significant at the 95% confidence level and there was variables from GIS, topography, and remote sensing data no significant co-linearity with other predictor vari- were used to impute measurements from the inventory ables. Following the calculation, the 2002 – 2006 cumu- ground plots across the forest site based on the three- lative NEP footprint probability density surface for each most-similar-neighbours (K-MSN with K = 3) (Crookston site (calculated by Chen et al. 2008 and Chen et al. and Finley 2008) using the Mahalanobis distance as a 2009) was used to weight each 10 m by 10 m footprint measure of similarity (Mahalanobis 1936). For each cell for the calculation of the site level mean value. footprint cell, the detailed inventory target measure- Site-level estimates of ΔCand ΔC stocks and were ments were estimated as the weighted mean of the calculated three ways. First, inventory- based site esti- three nearest neighbours (inversely weighted by the mates ΔC and ΔC stocks were calculated as an arith- Mahalanobis distance) and used for calculation of site- metic mean of ground plots at each site. Second, ΔC level means of ΔC. To assess the error of the model, stocks were calculated with the K-MSN prediction. the root mean squared difference (RMSD) was calcu- Third, and finally, ΔC stocks were calculated with the lated based on leave-one-out cross validation (Stage K-MSN prediction weighted by flux footprint probability and Crookston 2007). The procedure was completed density. Ferster et al. Forest Ecosystems (2015) 2:13 Page 8 of 19 Flux tower ∑NEP Comparing convergence of ΔC, EC tower ΣNEP, and CBM- Measurement and calculations of flux tower annual CFS3 ΣNEP NEP published by Black et al. (2008) and Krishnan Comparisons were made among the estimates by first et al. (2009) were summed to find the cumulative evaluating how each method ranked the stand ages. Sec- ∑NEP (2004–2006) used for this study. These values ond, estimates of EC tower ΣNEP and CBM-CFS3 ΣNEP were gap filled for conditions at night when friction vel- were compared. Third, at each stand age from youngest ocity was low (i.e. inadequate turbulent mixing) and to oldest, the estimates of ΔC and tower and CBM-CFS3 corrected for energy-balance closure. Annual NEP ΣNEP were compared for convergence by comparing values for 2003, 2004, 2005, and 2006 were summed for the means, variance around the means indicated by the each site to estimate the ΣNEP total for the four-year standard deviation, and spatial variance indicated by the period. Morgenstern et al. (2004) reported that uncer- RMSD. Comparisons of convergence among methods tainty in the annual NEP measured using EC at the were made for the following estimates: ΔC as an arith- −1 near-rotation stand may be as much as 0.9 MgC ha metic mean of plots (ΔC unweighted), ΔC with K- 2 2 −1 year (due to systematic error in the EC measure- MSN classification (ΔC K-MSN), ΔC with K-MSN 2 2 ments). This suggests that the uncertainty in the esti- classification and footprint weighting (ΔC K-MSN −1 mate at the near-rotation stand may be ±3.6 MgC ha FPW), CBM-CFS3 ΣNEP as an arithmetic mean of cells over the full four-year measurement interval. Detailed within the tower footprint, CBM-CFS3 ΣNEP as a foot- estimates of the uncertainty that propagates through print weighted mean (ΣNEP CBM-CFS3 FPW), and EC the calculation based on the friction velocity thresholds tower ΣNEP. (e.g. following Richardson and Hollinger 2007) was be- yond the scope of this paper. Results Live biomass ΔC For all sites, there was a positive change in the live C CBM-CFS3 ∑NEP mass of overstory trees (Table 2). For the near-rotation CBM-CFS3 is a forest C accounting model used for na- stand and juvenile stands, the increase was larger than at tional reporting of annual C inventories for Canada’s the regenerating site. The increase in live biomass C was managed forests (Kurz et al. 2009). This model uses nearly equal at the near-rotation and juvenile stands. An growth and yield equations and allometric equations to increase in shrub, herb, and bryophyte understory ΔC estimate tree growth and stand net primary production was observed at the near rotation and regenerating (NPP). The C mass of overstory trees and overstory tree stands, while a small decrease in understory ΔC stocks roots is use to estimate litterfall and root mortality occured at the juvenile stands, possibly related to an ex- transfers to the aboveground and belowground detritus pected increase in canopy closure as the overstory trees C pools, respectively. A soil submodel estimates de- matured. Comparing the two juvenile stands (HDF88 composition and heterotrophic respiration (Rh) based and HDF90), HDF88 had higher live biomass C than on the size of various detrital pools and mean annual HDF90 (for example, live stem C and total live C were temperature. The model tracks all major C pools to en- more than 1 standard deviation larger). sure 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 Detritus ΔC D1 Oyster River and encompassing the DF49, HDF90, and Dead standing tree C decreased in the near-rotation HDF00sites (Figure1). Modelrunsfor theareawere stand due to dead standing trees falling during the meas- performed on 1 hectare grid cells of forest disturbance urement interval, and there was a small increase at the history data, forest cover data, growth and yield equations, other sites due to tree mortality (Table 2). Large woody and disturbance transition matrices from Trofymow debris increased in the near-rotation stand, and de- et al. (2008). The modeled ∑NEP values for the DF49, creased at the other sites. There was considerable vari- HDF00 sites (from Wang et al. 2011) were calculated ability among plots, indicated by the high standard as site-level means unweighted and weighted by the deviations compared to the means. Fine woody debris flux probability distribution for the model cells in the decreased at all sites, except one of the juvenile sites footprint area; and at HDF90 as an unweighted mean (where there was a small increase). The large decrease in of model cells over the spatial extent of the ground fine woody debris at the regenerating clearcut was likely plots. Annual model values of NEP for 2003, 2004, due to decomposition of debris from the recent harvest. 2005 and 2006 were summed to calculate the ∑NEP Finally, forest floor material increased at all sites, with and mean annual NEP for the four year period the smallest increase at the regenerating clearcut, and (Figure 4). the largest increase at the juvenile stands. Ferster et al. Forest Ecosystems (2015) 2:13 Page 9 of 19 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. Mineral soil ΔC Results from the annual accounting method for de- Soil C (0–15 cm) stocks was highest at the near-rotation tritus, showed a net gain for the sum of components at site, followed by the two juvenile stands, and lowest at the near-rotation stand, and a decrease at the other sites the regenerating clearcut (Table 2). Over the measure- (Table 3). The increase in C at the near-rotation stand ment interval, large increases in mineral soil (0–15 cm) was slightly larger using ΔC than the direct measure- D2 ΔC were measured at all sites, with the highest at one ment method; however, the biggest difference was at the of the juvenile stands (HDF88). juvenile stand, which showed a net loss using ΔC , and D2 a large positive balance using the direct measurement Detritus annual accounting method ΔC method. For stumps, this method indicated that decom- D2 Measured litterfall was highest at the near-rotation position may have been greater than was initially expected stand, followed by the juvenile stands, and the regener- (stumps were not re-measured due to an expectation of ating stand had the lowest amount (Figure 5). The ma- minimal decomposition), especially at the more recently jority of litterfall at the near rotation stand and juvenile disturbed juvenile and regenerating clearcut sites. Large stands was needles followed by twigs. At the near- woody debris increased at the near-rotation stand, and de- rotation stand, litterfall was highest in 2001–2002, de- creased at the other sites, due to less dead standing trees creased through 2004–2005, and slightly increased in to fall and less branchfall. The increase in fine woody deb- 2005–2006. The other sites were relatively constant ris was estimated to be much lower using ΔC ,with the D2 through time. near-rotation stand having a net-accumulation of fine Ferster et al. Forest Ecosystems (2015) 2:13 Page 10 of 19 Table 2 Inventory-based site C stocks and ΔC estimates using arithmetic plot means (standard deviation) DF49 HDF90 HDF88 HDF00 Live Foliage C 2002 12.7 (1.8) 2.7 (0.9) 2.6 (0.8) 0 (0.1) Foliage ΔC −0.1 (0.7) 1.2 (0.7) 1.5 (0.4) 0.4 (0.5) Branch C 2002 22.1 (4) 2.5 (1) 1.5 (0.4) 0 (0.1) Branch ΔC 0.7 (1.2) 1.3 (0.9) 1.8 (0.5) 0.1 (0.2) Bark C 2002 17.6 (3.8) 0.7 (0.1) 1.1 (0.4) 0 (0) Bark ΔC 0.9 (0.9) 0.7 (0.3) 1 (0.3) 0.1 (0.2) Stem C 2002 90.6 (23.4) 2.4 (0.4) 4.7 (2) 0.1 (0.3) Stem ΔC 5.7 (4.7) 3.3 (1.1) 4.4 (1.3) 0.6 (1.1) Roots C 2002 32.6 (8.3) 1.8 (0.3) 2.8 (1.3) 0.1 (0.1) Roots ΔC 1.6 (1.6) 1.8 (0.6) 3.9 (0.9) 0.7 (1.5) Understory C 2002 1.3 (0.6) 8.3 (1.9) 11.2 (4.7) 4.9 (1.8) Understory ΔC 1.3 (1.8) −0.2 (2.8) −2.6 (4.5) 1.7 (1.9) Total Live C 2002 177 (39.9) 18.5 (3.7) 24 (4.7) 5.2 (2) Total Live ΔC 10.2 (8.7) 8.1 (6) 10 (5.9) 3.7 (4.4) Detritus Standing Dead C 2002 21.8 (13.9) 0.2 (0.2) 0 (0.1) 0 (0.2) Standing Dead ΔC −3.4 (11.1) 0.5 (0.6) 0.1 (0.2) 0 (0.2) Stumps C 2002 6.4 (3.2) 5.7 (2.8) 11.7 (6.1) 11.9 (10.6) Stumps ΔC 0 (0) 0 (0) 0 (0) 0 (0) Large Woody Debris C 2002 34.9 (20.5) 12.3 (8.7) 43.3 (25.3) 14.4 (11.1) Large Woody Debris ΔC 1.6 (12.4) −3.5 (9.7) −1.5 (18.8) −0.9 (6.3) Fine Woody Debris C 2002 3.1 (1) 1.3 (0.7) 0.5 (0.6) 6.1 (2.9) Fine Woody Debris ΔC −0.9 (1) −0.7 (0.8) 0.8 (0.8) −4.5 (2.6) Forest Floor C 2002 13.8 (7.3) 12.1 (7.2) 18.8 (13.1) 18.2 (11.6) Forest Floor ΔC 11 (6.5) 1.9 (11.2) 24.8 (19.7) 0.1 (12.3) Total Detritus C 2002 79.9 (32.7) 31.6 (8.4) 74.4 (38.7) 50.6 (16.8) Total Detritus ΔC 8.2 (21) −1.9 (8.1) 24.2 (30.3) −5.4 (12.6) Soil Mineral Soil 0–15 cm C 2002 38.6 (21.9) 29 (8.9) 34.6 (16.1) 23.6 (11.3) Mineral Soil 0–15 cm ΔC 18.8 (24.1) 27.8 (27.7) 10.7 (32.4) 24.2 (15.2) Totals Total Live ΔC 18.3 (23.4) 6.2 (10.5) 34.2 (30.6) −1.7 (14.9) + Total Detritus ΔC Total Live ΔC 37.2 (34.7) 34 (29.8) 44.9 (41.3) 22.5 (14.2) + Total Detritus ΔC + Mineral Soil 0 - 15 cm ΔC −1 All units are Mg C ha . woody debris, the juvenile sites had a small loss, and the layer, the near-rotation stand was estimated to have a large regeneration site experienced the largest decrease. This in- increase, while the other sites decreased, with the greatest dicates that, given the available inputs to litterfall and ex- decrease at the regenerating clearcut. pected rates of decomposition for fine woody debris, the Since the plot values for ΔC were judged more reli- D2 direct measured stock changes in fine woody debris was able, all subsequent analyses to spatially scale the plot very likely unrealistically high. Finally, for the forest floor data to the site were made using ΔC . 2 Ferster et al. Forest Ecosystems (2015) 2:13 Page 11 of 19 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. Spatial distribution of C stocks and ΔC stocks sites the other spatial predictor variables showed finer Predictor variables (Table 4) were selected and used for scale variation. At the regenerating clearcut, NDVI was theK-MSN procedureatthe siteswithfluxtowers a significant predictor for woody debris, forest floor, (Figure6). Slopewas an importantpredictor forall sites and mineral soil C stocks. as it was strongly correlated with soil C stocks (with The Mahalanobis distance (Figure 7) demonstrates less sloping plots having higher soil (0 – 15 cm) C how comprehensively the inventory ground plots repre- stocks), and at the regenerating clearcut, slope was also sent stand conditions across the footprint. Lower per- correlated with living tree C stocks. At the near rota- centiles represent small distances that were relatively tion site, forest-inventory mapping was a significant well represented by the ground plots. Areas with the predictor variable for overstory C stocks. Due to the smallest Mahalanobis distances and best ground plot smaller footprints at the other sites, forest-inventory representation were close to the towers in the areas con- mapping varied little within the site; however, at all tributing most to NEP measurements (within the 50% Table 3 Inventory-based site detrital C in 2002 and ΔC (2002–2006) using plot mean (standard deviations) litterfall, D2 mortality and decomposition values DF49 HDF90 HDF88 HDF00 Detritus C Standing Dead 21.8 (13.9) 0.2 (0.2) 0 (0.1) 0 (0.2) ΔC Standing Dead −3.4 (11.1) 0.5 (0.6) 0.1 (0.2) 0 (0.2) D2 C Stumps 6.4 (3.2) 5.7 (2.8) 11.7 (6.1) 11.9 (10.6) ΔC Stumps ΔC −0.5 (0.2) −0.4 (0.2) −0.9 (0.4) −0.9 (0.8) D2 C Large Woody Debris C 34.9 (20.5) 12.3 (8.7) 43.3 (25.3) 14.4 (11.1) ΔC Large Woody Debris 3.9 (9.6) −1 (1) −4.3 (2.6) −1.4 (1) D2 C Fine Woody Debris 3.1 (1) 1.3 (0.7) 0.5 (0.6) 6.1 (2.9) ΔC Fine Woody Debris 1 (0.8) −0.4 (0.3) −0.1 (0.2) −2.3 (1.1) D2 C Forest Floor 13.8 (7.3) 12.1 (7.2) 18.8 (13.1) 18.2 (11.6) ΔC Forest Floor 8.2 (4.1) −4.8 (4.2) −7.4 (7.2) −9.3 (6.7) D2 C Total Detritus 79.9 (32.7) 31.6 (8.4) 74.4 (38.7) 50.6 (16.8) ΔC Total Detritus 9.2 (8.7) −6.1 (4.2) −12.5 (9.4) −13.8 (6.5) D2 −1 All units are Mg C ha . Ferster et al. Forest Ecosystems (2015) 2:13 Page 12 of 19 Table 4 Environmental predictor variables used to determine stratification units at DF49, HDF88, and HDF00 Selection Coverage Variable Units DF49 HDF88 HDF00 Forest-Inventory Site Index m selected Top Height m Disturbance History (Trofymow et al. 2008) 1st Harvest year 2nd Harvest year 1st Fire year 2nd Fire year Date Est. 1953 year Date Est. 2003 year 1st Fertilization year 2nd Fertilization year Fire Cause 1 nominal Fire Cause 2 nominal Topography aspect azimuth selected selected selected Elevation m asl selected Slope degrees selected selected selected SCOSA selected selected SSINA selected selected TSRAI selected selected selected −1 1999 Orthophoto Dominant Canopy Tree Density stems ha selected selected (Gougeon 1995) 2004 Multispectral NDVI NDVI selected selected Forest-Inventory Cover Species (Trofymow et al. 2008) Site Species Soil Survey of Canada (Jungen 1985) Most Common Soil Association CFS and Forest Companies Site Series selected selected (Trofymow et al. 2008) Site Index: Tree height at 50 years age at breast height (1.3 m). Topographic Variables: SCOSA = Slope × cos (Aspect), SSINA = Slope × sin (Aspect); TSRAI (Topographic Solar Radiation Aspect Index) = (1-cos((π/180°)(Aspect-30°)))/2 (Roberts and Cooper 1989). Topography from 2004 LiDAR survey at DF49 and HDF00 (Coops et al. 2007). Topography from 1:50,000 National Topographic Series map at HDF88. Multispectral data from 2004 Quickbird survey at DF49 and HDF00. Multispectral data from 2004 Landsat scene at HDF88. cumulative flux probability density isoline). Notable For the juvenile stands, K-MSN and K-MSN FPW ΔC areas with large Mahalanobis distances included areas was less than 1 standard deviation higher than the un- with different forest cover types (e.g., patches of hard- weighted mean. The K-MSN estimate was higher than woods at the near-rotation stand), very few overstory the K-MSN FPW estimate. For ΔC , K-MSN and K- D2 trees (at the juvenile stand), and short steep slopes (be- MSN FPW estimates were slightly higher than the mean tween bench terraces at the regenerating clearcut). (less than 1 standard deviation). For the total ΔC , the For site level estimations of ΔC at the regenerating estimates were less than 1 standard deviation different, clearcut, K-MSN was more than 1 standard deviation and the standard deviations were very large compared to lower than the unweighted mean, and K-MSN-FPW was the means, indicating a large amount of variability at the 1 standard deviation lower than K-MSN. For ΔC K- site. D2, MSN and K-MSN FPW were less than 1 standard devi- At the near-rotation stand, ΔC K-MSN and K-MSN ation higher than the unweighted mean (Tables 3 and 5). FPW were slightly higher than the arithmetic mean (less The totals for ΔC K-MSN and K-MSN FPW were less than 1 standard deviation). The estimates using KMSN than 1 standard deviation higher than the arithmetic and KMSN-FPW were identical. For ΔC , the un- D2 mean (Table 6). weighted mean was higher than K-MSN and K-MSN Ferster et al. Forest Ecosystems (2015) 2:13 Page 13 of 19 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. FPW, but less than 1 standard deviation different. For regenerating clearcut as a source (net C release). While ΔC total, all estimates were within 1 standard deviation. the ranking was consistent for all methods, there were Comparing the ΔC and ΣNEP methods at the juvenile differences in the magnitude of values for each stand stands, HDF88 was compared with measurements from (Table 6). the flux tower installed at the site, and HDF90 was com- At the regenerating clearcut, the CBM-CFS3 estimate pared with CBM-CFS3 model output. While the two of ΣNEP with footprint weighting converged more stands were similar in seral stage, species composition, closely with the EC tower estimate of ΣNEP than the and stand history, there were differences in stand com- unweighted estimate, and both were within 1 standard position. For example, HDF88 had higher live C stocks, deviation. Evaluating the yearly means (Figure 8a), dem- in particular, live stem C, understory stem C, and total onstrated that footprint weighting improved conver- stem C. In addition, HDF88 had larger detritus C stocks gence between the CBM-CFS3 estimate of ΣNEP and including stumps C, large woody debris C, and total de- the EC flux ΣNEP. This was primarily due to less tritus C. The stand variability at HDF88 was also greater weighting for patches of uncut forest in the periphery of than at HDF90 indicated by the large standard devia- the flux footprint (Figure 4). tions compared to the means. Given the large amount of At the juvenile sites, CBM-CFS3 ΣNEP indicated that spatial variability at HDF88, spatial scaling, and footprint HDF90 was a weak sink (net C uptake) and EC tower weighting had a notable effect on the amount and sign ΣNEP indicated that HDF88 was a weak source (net C at this site. release) (Tables 3 and 6). At the near rotation stand, the estimates of ΣNEP from CBM-CFS3 were higher than Comparing inventory ΔC , EC tower ΣNEP and CBM-CFS3 the EC flux tower ΣNEP, with the footprint weighted es- ΣNEP timate being the highest. Considering the means of All methods ranked the near-rotation stand as a sink for ΣNEP from CBM-CFS3 on a yearly basis, footprint C (net C uptake), the juvenile stands as a weak sink or weighting decreased convergence with the EC flux tower weak source (small C uptake or release), and the measurement of ΣNEP (Figure 8b). Ferster et al. Forest Ecosystems (2015) 2:13 Page 14 of 19 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. Table 5 Site mean (± RMSD) changes in live biomass (ΔC and detrital (ΔC ) C stocks calculated using three most L) D2 similar neighbours classification (K-MSN), and K-MSN with footprint weighting means (K-MSN FPW) of plots for sites with flux-towers DF49 HDF88 HDF00 Pool Component K-MSN K-MSN FPW K-MSN K-MSN FPW K-MSN K-MSN FPW Live ΔC Foliage 0.1 ± 1.7 0.2 (0.4) ± 1.7 1.3 ± 2 1.3 (0.3) ± 2 0.2 ± 0.9 0.1 ± 0.9 ΔC Branches 0.9 ± 1.6 1 (0.6) ± 1.6 1.9 ± 1.3 1.8 (0.6) ± 1.3 0.1 ± 1.1 0 ± 1.1 ΔC Bark 1.1 ± 1.5 1.1 (0.4) ± 1.5 1 ± 1.7 1 (0.3) ± 1.7 0.1 ± 0.8 0 ± 0.8 ΔC Stem 6.7 ± 1.4 6.1 (2.2) ± 1.4 4.4 ± 1.7 4.2 (1.2) ± 1.7 0.4 ± 0.9 0.1 ± 0.9 ΔC Understory 1.8 ± 1 2.5 (2.3) ± 1 −0.9 ± 0.8 −1.7 (5.2) ± 0.8 0.8 ± 1.3 0.6 ± 1.3 ΔC Roots 2 ± 1.5 1.9 (0.8) ± 1.5 4.3 ± 2 4.4 (0.6) ± 2 0.3 ± 1.1 0.1 ± 1.1 ΔC 12.8 ± 1.6 12.8 (3) ± 1.6 11.9 ± 0.8 10.9 (7.2) ± 0.8 1.8 ± 0.8 1 ± 0.8 LTotal Detritus ΔC Standing Dead −5.7 ± 1.2 −8.1 (16.1) ± 1.2 0.1 ± 1.3 −0.7 (0.4) ± 1.2 0 ± 1.4 0 ± 1.4 D2 ΔC Stumps −0.5 ± 1.1 −0.5 (0.1) ± 1.1 −0.6 ± 1.2 −3.5 (1.1) ± 1.9 −0.7 ± 1.3 −0.8 ± 1.3 D2 ΔC Large Woody Debris 3.5 ± 0.9 6.3 (13.7) ± 0.9 −3.5 ± 1.9 −0.1 (0.2) ± 1.1 −1.5 ± 0.3 −1.5 ± 0.3 D2 ΔC Find Woody Debris 0.6 ± 1.4 0.2 (0.8) ± 1.4 −0.1 ± 1.1 −6.9 (4.3) ± 1.5 −2 ± 1.3 −1.6 ± 1.3 D2 ΔC Litter 7.1 ± 0.9 8.7 (2.5) ± 0.9 −6.2 ± 1.5 −0.2 (11.2) ± 1.6 −6.1 ± 1.3 −4.8 ± 1.3 D2 ΔC 5 ± 1.5 −10.3 ± 1.6 1.3 (0.3) ± 2 −10.4 ± 1.3 −8.7 ± 1.3 D2Total Total ΔC = ΔC + ΔC 17.8 ± 1.2 19.5 (5.9) ± 1.5 1.6 ± 1.2 −0.2 ± 1.2 −8.6 ± 1.3 −7.7 ± 1.3 2 LTotal D2Total The Root Mean Square Difference (RMSD) is presented as an estimate of spatial estimation error. Ferster et al. Forest Ecosystems (2015) 2:13 Page 15 of 19 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 −1 −1 −1 2003–2006 Mg C ha and the (b) mean annual value for each method in Mg C ha yr a) Site Method DF49 HDF90 HDF88 HDF00 ΔC Mean 19.4 (7.2) 2.0 (8.5) −2.5 (13.3) −10.1 (8.2) ΔC K-MSN 17.8 (5.8) ± 1.2 na/ 1.6 (11.7) ± 1.2 −8.6 (5.0) ± 1.3 ΔC K-MSN FPW 19.5 (5.9) ± 1.2 na/ −0.2 (11.2) ± 1.2 −7.7 (3.5) ± 1.3 ΣNEP Flux Tower 13.6 ± 3.6 na/ −1.9 −20.1 ΣNEP CBM-CFS3* 18.9 2.0* na/ −11.5 ΣNEP CBM-CFS3 FPW 24.1 na/ na/ −17.9 b) Site Method DF49 HDF90 HDF88 HDF00 ΔC Mean 4.9 (1.8) 0.5 (2.1) −0.63 (3.3) −2.5 (2.0) ΔC K-MSN 4.5 (1.5) ± 1.2 na/ 0.4 (2.9) ± 1.2 −2.2 (1.2) ± 1.3 ΔC K-MSN FPW 4.9 (1.5) ± 1.2 na/ −0.1 (2.8) ± 1.2 −1.9 (1.9) ± 1.3 ΣNEP Flux Tower 3.4 ± 0.9 na/ −0.5 −5.0 ΣNEP CBM-CFS3* 4.7 0.5* na/ −2.9 ΣNEP CBM-CFS3 FPW 6.0 na/ na/ −4.47 *average of all grid cells encompassing the six plots 2003 – 2006. na – not available or not applicable. At the regenerating clearcut, unweighted ΔC was have used measurements or estimates of soil and de- within 1 standard deviation of the CBM-CFS3 estimates tritus respiration in combination with inventory mea- of ΣNEP, and closest to unweighted CBM-CFS3 ΣNEP. surements of overstory productivity (Law et al. 2001; The unweighted mean of ΔC , ΔC K-MSN, and ΔC Ehman et al. 2002; Kolari et al. 2004; Black et al. 2005). 2 2 2 K-MSN FPW were greater than 1 standard deviation In the present study, comparisons were based on mea- higher than the EC flux tower estimate of ΣNEP. surements of ΔC which were higher than EC ΣNEP. ΔC K-MSN was within 1 standard deviation of the Black et al. (2005) also found a divergence between in- CBM-CFS3 unweighted ΣNEP. Both ΔC K-MSN and ventory measurements of ΔC (including ΔC measure- 2 S ΔC K-MSN FPW were greater than 1 standard ments) and EC ΣNEP, with ΔC higher than EC ΣNEP; deviation higher than ΣNEP CBM-CFS 3 FPW however, they also calculated ΣNEP using measure- (Table 6). ments of heterotrophic respiration combined with bio- At the juvenile stands, all estimates ΔC were within 1 metric estimates of NPP, and found that the bias was standard deviation of the estimates of ΣNEP. At HDF90, not observed. They concluded that ΔC systematically the unweighted estimate of ΔC matched ΣNEP CBM- overestimated NEP due to unaccounted decomposition CFS3. At HDF88, the unweighted mean of ground plots processes and uncertainties in ΔC . At DF49, Jassal was the closest to the EC tower measurements of ΣNEP, et al. (2010) measured heterotrophic soil respiration in while the ΔC K-MSN had a different sign (indicating a 2007, which was a major portion of total soil respir- small source), and the K-MSN ΔC FPW had the same ation. Including measurements of soil heterotrophic sign as ΔC unweighted (Table 6). respiration with net primary production estimates from At the near-rotation stand, all estimates of ΔC were forest inventory may improve convergence with EC- within 1 standard deviation of all estimates of ΣNEP. flux tower measurements. However, a limitation of The estimates of ΔC (and CBM-CFS3 ΣNEP) higher taking measurements of soil respiration is that the mea- than the EC flux tower measurements of ΣNEP surements require long-term collection using special (Table 6). equipment and are not typically collected in traditional forest inventory. Discussion Kolari et al. (2004) found that management practices, Several approaches can be used to compare forest inven- such as clearcut harvesting had a large impact on the tory data with EC ΣNEP. For example, several authors soil C balance, with clearcut sites soil C sources, while Ferster et al. Forest Ecosystems (2015) 2:13 Page 16 of 19 than expected given the measured litterfall inputs and estimated decomposition rates. The NFI (2008) calls for 10-year remeasurement intervals, which may be appro- priate for trees and detritus, but may be too short for ΔC . In addition, direct use of forest inventory data re- lies 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 allomet- ric equations is labour intensive, destructive, and costly, and therefore uncommon in most studies; however, use of existing relationships can lead to errors due to differ- ences 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 pro- cedure provided regionally applicable allometric equa- tions; however, the errors in this large pool of equations have not yet been estimated and therefore is a limitation. The measured changes in ΔC 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 −1 −1 MgC ha year in a Douglas-fir stand and Gielen et al. (2013) found no significant change over an eight-year Figure 8 Yearly NEP estimates for CBM-CFS3 and EC flux tower (Black interval in a Scots pine forest. In other studies by Ehman et al. 2008) for (A) the regenerating clearcut and (B) the near-rotation et al. (2002), Kolari et al. (2004), Miller et al. (2004), stand. CBM-CFS3 is shown as yearly mean of footprint cells with equal Ohtsuka et al. (2007), and Granier et al. (2008) ΔC was weighting (Fp equal) and weighting by NEP flux probability density distribution (Fp weight). assumed to be zero over the measurement interval. The large value for ΔC in this study may be explained by several factors related to measurement methodology. non-disturbed sites had more stable soil C balances. In First, mineral soil bulk densities in the 0–15 cm layer this study, a trend in soil C stocks, with the highest soil were assumed not to change and thus the BD values C stocks at the oldest site, and the lowest at the recently measured in 2002 were also used in 2006. If the bulk disturbed site, was consistent with the trend by seral density decreased then the C change would be overesti- stage reported by Kolari et al. (2004). For the two mated. Second, at the regenerating clearcut and juvenile similar-aged juvenile stands, since there was no flux stands, an increase in fine root mass from 2002 to 2006 tower installed at HDF90 and CBM-CFS3 model runs included in the <2 mm soil fraction could also have results were not available at HDF88, a direct comparison accounted for some of the increase in soil C. Third, sam- was not possible. In addition to the differences in meas- ples in 2002 were taken by 10–12 cm diameter hole ex- urement methods, the differences in ΣNEP between the cavations, while samples in 2006 were collected using a two juvenile stands could also be partially attributed to 2 cm diameter soil corer, which, due to compression at site conditions. Notably, HDF88 had more detrital C the opening of the sample corer, may be biased to soil stocks, likely leading to larger releases of C due to de- from shallower depths where C content is likely higher. composition. This was reflected in the more negative Therefore, in this study, measurements of ΔC were ΔC at HDF88 compared to HDF90. deemed unreliable, due to the large positive changes in D2 For the forest inventory measurement of ΔC and its ΔC , and assumed to be negligible over the four year components ΔC , ΔC , and ΔC , differences in proced- period and thus not used for subsequent calculations of L D S ure, measurement error by field and lab crews, and site level ΔC and further comparisons. In future work, sample design may introduce error into inventory mea- efforts to reduce any chances of variation in the way the surements that may be larger than the magnitude of samples are collected and analysed may result in more change in C stocks over short measurement intervals, reliable measurements of C stock changes. For example, such as the 4 year period used in this study. Overall, collecting and processing the samples using identical ΔC and ΔC were much larger than reported in previ- methods at nearly the same locations, or collecting a D1 S ous studies (e.g. Law et al. 2001), and was much larger much larger number of samples to capture a wider range Ferster et al. Forest Ecosystems (2015) 2:13 Page 17 of 19 of spatial variation could reduce the effect of spatial vari- “different strengths and weaknesses, but the combination ability on measurements of ΔC . of multiple measurements and modeling has the potential Several authors have applied footprint analysis in an for refining estimates of C stocks and fluxes” (Law et al. effort to correct for non-homogenous environments, 2004). For example, eddy-covariance data are temporally and horizontal advection, for example due to daily up- detailed, and provide information about daily and seasonal slope or downslope winds, which may introduce bias processes. However, the data are less spatially extensive into the eddy-covariance measurements (Baldocchi and detailed than the forest inventory data, for example, 2008). For example, Ehman et al. (2002) used footprints only three sites were available in this study each covering to calculate an estimate of sensor bias given using a an area with diverse vegetation within the tower footprint. vegetation index, and this gave confidence that the sen- Measurement error was estimated to be relatively low, on −1 −1 sor measurements were representative of the inventory the order of 0.9 Mg C ha year (Morgenstern et al. measurements. Gielen et al. (2013) used a more complex 2004). In contrast, forest inventory data are more spatially footprint model applied to the processing of NEP data, extensive and include detailed measurements of forest at- were fluxes originating outside of the area where inven- tributes. The main constraint for forest inventory data is tory measurements were sampled were removed from that they are less temporally detailed (for example, mea- the NEP estimates. In this study, footprint weighting by surements were available at only two points in time), and ecological attributes and footprint probability density the error in these measurements may be greater than the was evaluated where it facilitated comparison between expected rate of change over the measurement interval. inventory-based measurements of ecosystem ΔC stocks The inventory plots in this study sampled a broad range and tower ΣNEP by accounting for fine scale spatial of vegetation conditions and provided insight into the rep- variability in forest structure. The effect of K-MSN clas- resentativeness of flux-tower footprints. Estimates of error −1 sification and footprint weighting had an effect on the of 1.2 to 1.3 Mg C ha were identified relating to spatial estimated site-level values; however, this effect was representativeness of the inventory plots within the tower smaller than differences due to measurement method- footprints. Finally, C budget models are important for un- ology. Notably, applying footprint weighting to the derstanding landscape scale C processes. Since these estimation of CBM-CFS3 ΣNEP at the regenerating models depend on numerous coefficients, algorithms and clearcut increased convergence with EC tower ΣNEP; assumptions, it is valuable to compare model outputs with however, footprint weighting did not have the same ef- EC flux and inventory measurements to assess their per- fect for ΔC . The largest differences observed among formance and evaluate their behaviour. Considering all of methods were seen at the regenerating clearcut, where these approaches together, the relative strengths and limi- the inventory approach indicated the site was a much tation of each approach can be taken into account in lower source than the ΣNEP methods. A possible cause evaluating the suitability of each type of measurement at is that the inventory methods did not properly account varying spatial and temporal scales. for C losses from decomposition of coarse roots from stumps and large woody debris at recently harvested Conclusions sites, which can be a substantial source of respiration The comparison of ΔC from inventory measurements (Janisch et al. 2005). Therefore, if an estimate of stump with ΣNEP from CBM-CFS3 and EC flux-tower measure- coarse root decomposition were included in ΔCD it ments demonstrated an agreement among the methods in may have increased convergence with the ΣNEP estimates. trends across stand seral stage; however, due to differences Additionally, the original sample locations were designed in the measurement approaches, there was some diver- to capture the range of conditions in near proximity of the gence in the results. Convergence amongst methods was tower (based on interpretation of aerial photographs), and closest for the juvenile sites (HDF90 and HDF88), which the relationship was similar when scaled up to the site were in transition from C source- to sink. At the most re- level. Future studies may seek to derive accurate footprint cently harvested site (HDF00), the EC flux ΣNEP and estimates prior to establishing sample plots to ensure they CBM CFS3 ΣNEP indicated that the site was a greater are representative of the broader site conditions. Another source of C over the time period than ΔC from inventory more costly alternative is to establish a much greater methods. At the near-rotation site the inventory ΔC number of plots in a systematic basis around the tower, to method, unweighted CBM CFS3 ΣNEP and EC flux ΣNEP ensure the spatial variation within the tower footprint is also converged. adequately captured. Each of the measurement methods had advantages This research highlights how different measurements and limitations. For example, while EC flux towers pro- and approaches for C accounting include different con- vide temporally rich data, they had limited spatial cover- straints, and as a result, computations using different data age. Obtaining ΔC using forest inventory measurements sources may not converge. Each type of measurement has included several challenges such as the consistency of Ferster et al. Forest Ecosystems (2015) 2:13 Page 18 of 19 data collected over long time intervals, small changes in Bailey JD, Mayrsohn C, Doescher PS, St Pierre E, Tappeiner JC (1998) Understory vegetation in old and young Douglas-fir forests of western Oregon. For Ecol ΔC compared to uncertainty in the measurements, and Manage 112:289–302 insufficient measurement of processes such as decom- Baldocchi DD (2008) “Breathing” of the terrestrial biosphere: lessons learned from position of all detrital pools. Forest inventory data are a global network of carbon dioxide flux measurement systems. Aust J Bot 56:1, doi:10.1071/BT07151 available over broad areas; therefore methods that com- Barford CC, Wofsy SC, Goulden ML, Munger JW, Pyle EH, Urbanski SP, Hutyra L, bine measurements with estimates of decomposition are Saleska SR, Fitzjarrald D, Moore K (2001) Factors controlling long- and needed to provide information over a regional scale. short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science 294:1688–1691, doi:10.1126/science.1062962 Footprint analysis using the spatial predictors from re- Black K, Bolger T, Davis P, Nieuwenhuis M, Reidy B, Saiz G, Tobin B, Osborne B mote sensing and topography can give confidence that (2005) Inventory and eddy covariance-based estimates of annual carbon sample locations represent broader area site conditions. sequestration in a Sitka spruce (Picea sitchensis (Bong.) Carr.) forest ecosystem. Eur J For Res 126:167–178, doi:10.1007/s10342-005-0092-4 Considering the advantages and constraints of each type Black TA, Fredeen A, Jassal R (2008) Carbon Sequestration in British Columbia’s of measurement is helpful in choosing appropriate meas- Forests and Management Options. University of Victoria, Victoria, BC, White urement methods at varying spatial and temporal scales. paper. Pacific Institute for Climate Solutions Brown JK (1971) A planar intersect method for sampling fuel volume and surface area. For Sci 17:96–102 Competing interests Chen B, Chen JM,Mo G,Black TA,WorthyDEJ (2008) Comparison of regional carbon The authors declare that they have no competing interests. flux estimates from CO 2 concentration measurements and remote sensing based footprint integration. Global Biogeochem Cycles. doi:10.1029/2007GB003024 Chen BB, Black TA, Coops NC, Hilker T, Trofymow JA, Morgenstern K, Black AT Authors’ contributions (2009) Assessing tower flux footprint climatology and scaling between CJF set research objectives, devised and undertook analysis, and composed the remotely sensed and eddy covariance measurements. Boundary-Layer manuscript. JAT set research objectives, directed the collection and synthesis of Meteorol 130:137–167, doi:10.1007/s10546-008-9339-1 field data, summarized output from CBM-CFS3 runs, advised analysis, assisted Clark D, Brown S (2001) Measuring net primary production in forests: concepts with manuscript composition, and provided extensive editorial work on the and field methods. Ecol Appl 11:356–370 manuscript. NCC advised on research objectives, advised analysis, and assisted Coops NC, Hilker T, Wulder MA, St-Onge B, Newnham GJ, Siggins A, Trofymow JA with manuscript composition. BC advised the footprint analysis, and provided (2007) Estimating canopy structure of Douglas-fir forest stands from discrete- commentary on manuscript. TAB advised on research objectives, advised the return LiDAR. Trees-Structure Funct 21:295–310. doi:10.1007/s00468-006-0119-6 analysis, and provide extensive review and commentary throughout the work. Crookston NL, Finley AO (2008) yaImpute: an R package for kNN imputation. All authors read and approved the final manuscript. J Stat Softw 23(10):1–16 Curtis P, Hanson P, Bolstad P, Barford C, Randolph J, Schmid H, Wilson K (2002) Acknowledgements Biometric and eddy-covariance based estimates of annual carbon storage in We thank Bob Ferris, Frank Eichel, and Glenda Russo of the Canadian Forest five eastern North American deciduous forests. Agric For Meteorol 113:3–19, Service and staff of B.A. Blackwell and Associates for their help processing and doi:10.1016/S0168-1923(02)00099-0 collecting National Forest-inventory-style ground plot data. We also thank Dixon RK, Solomon AM, Brown S, Houghton RA, Trexier MC, Wisniewski J (1994) François Gougeon, CFS, for determining canopy tree density values for the Carbon pools and flux of global forest ecosystems. Science 263:185–190, sites and to Graham Stinson, CFS, for the CBM-CFS3 output files used in doi:10.1126/science.263.5144.185 Wang et al. (2011). We thank the UBC Land and Food Systems Biometeorology Ehman JL, Schmid HP, Grimmond CSB, Randolph JC, Hanson PJ, Wayson CA, and Soil Physics Group, in particular Paul Jassal, Praveena Krishnan, Kai Cropley FD (2002) An initial intercomparison of micrometeorological and Morgenstern, and Elyn Humphreys for their work processing EC flux data, ecological inventory estimates of carbon exchange in a mid-latitude deciduous and Zoran Nesic, Dominic Lessard, Andrew Sauter, and Andrew Hum for forest. Glob Chang Biol 8:575–589, doi:10.1046/j.1365-2486.2002.00492.x their work running and maintaining the EC flux tower sites. We also thank Ferster CJ, Trofymow JA, Coops NC, Chen B, Black TA, Gougeon FA (2011) Mark Johnson for his comments and suggestions. Funding for this study Determination of ecosystem carbonstock distributions in the flux footprint of an was provided by the Canadian Forest Service Pacific Forestry Centre Graduate eddy-covariance tower in a coastal forest in British Columbia. Can J For Res Student Award, a CFCAS grant to the Canadian Carbon Program (CCP), and 41:1380–1393, doi:10.1139/X11-055 Natural Sciences and Engineering Research Council of Canada (NSERC) Gielen B, De Vos B, Campioli M, Neirynck J, Papale D, Verstraeten A, Ceulemans R, Discovery Grant to NCC. The LiDAR data were acquired by Benoit St-Onge Janssens IA (2013) Biometric and eddy covariance-based assessment of of the University of Quebec at Montreal as part of an ongoing collaborative decadal carbon sequestration of a temperate Scots pine forest. Agric For project with funds provided by NSERC and BIOCAP. Meteorol 174:135–143, doi:10.1016/j.agrformet.2013.02.008 Gillis M, Omule A, Brierley T (2005) Monitoring Canada’s forests: the national Author details forest inventory. For Chron 81:214–221 Department of Forest Resources Management, University of British Göckede M, Rebmann C, Foken T (2004) A combination of quality assessment Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada. Canadian Forest tools for eddy covariance measurements with footprint modelling for the Service (Pacific Forestry Centre), Natural Resources Canada, 506 Burnside characterisation of complex sites. Agric For Meteorol 127:175–188, Road West, Victoria, BC V8Z 1M5, Canada. Faculty of Land and Food doi:10.1016/j.agrformet.2004.07.012 Systems, University of British Columbia, 2357 Main Mall, Vancouver, BC V6T 4 Gougeon F (1995) A crown-following approach to the automatic delineation of 1Z4, Canada. Department of Biology, University of Victoria, PO Box 1700, individual tree crowns in high spatial resolution aerial images. Can J Remote Station CSC, Victoria, BC V8W 2Y2, Canada. Sens 21:274–284 Gough C, Vogel C, Schmid H, Su H, Curtis P (2008) Multi-year convergence of Received: 18 December 2014 Accepted: 21 April 2015 biometric and meteorological estimates of forest carbon storage. Agric For Meteorol 148:158–170, doi:10.1016/j.agrformet.2007.08.004 Granier A, Ceschia E, Damesin C, Dufrene E, Epron D, Gross P, Lebaube S, Le Dantec V, Le Goff N, Lemoine D, Lucot E, Ottorini JM, Pontailler JY, Saugier B References (2000) The carbon balance of a young Beech forest. Funct Ecol 14:312–325, Babst F, Bouriaud O, Papale D, Gielen B, Janssens I, Nikinmaa E, Ibrom A, Wu J, doi:10.1046/j.1365-2435.2000.00434.x Bernhofer C, Köstner B, Grünwald T, Seufert G, Ciais P, Frank D (2014) Above-ground woody carbon sequestration measured from tree rings is Granier A, Bréda N, Longdoz B, Gross P, Ngao J (2008) Ten years of fluxes and coherent with net ecosystem productivity at five eddy-covariance sites. New stand growth in a young beech forest at Hesse, North-eastern France. Ann Phytol 201:1289–1303, doi: 10.1111/nph.12589 For Sci 65:704–704, doi:10.1051/forest Ferster et al. Forest Ecosystems (2015) 2:13 Page 19 of 19 Grant RF, Black TA, Humphreys ER, Morgenstern K (2007) Changes in net Ohtsuka T, Mo W, Satomura T, Inatomi M, Koizumi H (2007) Biometric based ecosystem productivity with forest age following clearcutting of a coastal carbon flux measurements and net ecosystem production (NEP) in a Douglas-fir forest: testing a mathematical model with eddy covariance temperate deciduous broad-leaved forest beneath a flux tower. Ecosystems measurements along a forest chronosequence. Tree Physiol 27:115–131 10:324–334, doi:10.1007/s10021-007-9017-z Green RN, Klinka K (1994) A Field Guide to Site Identification and Interpretation R Development Core Team (2011) R: a language and environment for statistical for the Vancouver Forest Region. Land Management Handbook Number 28. computing. R Foundation for Statistical Computing, Vienna, Austria, URL British Columbia Ministry of Forests, Research Program, Victoria, BC http://www.R-project.org/ Humphreys ER, Black TA, Morgenstern K, Cai T, Drewitt GB, Nesic Z, Trofymow JA Richardson AD, Hollinger DY (2007) A method to estimate the additional (2006) Carbon dioxide fluxes in coastal Douglas-fir stands at different stages uncertainty in gap-filled NEE resulting from long gaps in the CO2 flux record. of development after clearcut harvesting. Agric For Meteorol 140:6–22, Agric For Meteorol 147:199–208, doi:10.1016/j.agrformet.2007.06.004 doi:10.1016/j.agrformet.2006.03.018 Roberts D, Cooper S (1989) Concepts and techniques of vegetation mapping. In: Janisch JE, Harmon ME, Chen H, Fasth B, Sexton J (2005) Decomposition of Ferguson D, Morgan P, Johnson FD (eds) Land classifications based on coarse woody debris originating by clearcutting of an old-growth conifer vegetation: applications for resource management. USDA Forest Service forest. Ecoscience 12:151–160, doi:10.2980/i1195-6860-12-2-151.1 General Technical Report: INT-257; USDA Forest Service Intermountain Forest and Range Experiment Station, Ogden, UT, USA, pp 90–96 Janzen HH (2005) Soil carbon: a measure of ecosystem response in a changing Schmid H (2002) Footprint modeling for vegetation atmosphere exchange world? Can J Soil Sci 85:467–480 studies: a review and perspective. Agric For Meteorol 113:159–183, Jassal RS, Black TA, Trofymow JA, Roy R, Nesic Z (2010) Soil CO and N Oflux 2 2 doi:10.1016/S0168-1923(02)00107-7 dynamics in a nitrogen-fertilized Pacific Northwest Douglas-fir stand. Geoderma Schmid HP, Lloyd CR (1999) Spatial representativeness and the location bias of 157:118–125, doi:10.1016/j.geoderma.2010.04.002 flux footprints over inhomogeneous areas. Agric For Meteorol 93:195–209, Jungen J (1985) Soils of Southern Vancouver Island. Victoria, BC doi:10.1016/S0168-1923(98)00119-1 Kolari P, Jukka P, Rannik U, Hannu I, Pertti H, Berninger F (2004) Carbon balance Schulze ED, Hogberg P, Van Oene H, Persson T, Harrison A, Read D, Kjoller A, of different aged Scots pine forests in Southern Finland. Glob Chang Biol Matteucci G (2000) Interactions Between the Carbon and Nitrogen Cycles 10:1106–1119, doi:10.1111/j.1365-2486.2004.00797.x and the Role of Biodiversity: A Synopsis of a Study Along a North–south Kominami Y, Jomura M, Dannoura M, Goto Y, Tamai K, Miyama T, Kanazawa Y, Transect Through Europe. In: Schulze ED (ed) Carbon and Nitrogen Cycling Kaneko S, Okumura M, Misawa N, Hamada S, Sasaki T, Kimura H, Ohtani Y in European Forest Ecosystems. Springer, Heidelberg, pp 468–491 (2008) Biometric and eddy-covariance-based estimates of carbon balance for Smyth C, Trofymow JA, Kurz W (2010) Decreasing Uncertainty in CBM-CFS3 a warm-temperate mixed forest in Japan. Agric For Meteorol 148:723–737, Estimates of Forest Soil Carbon Sources and Sinks Through use of Long-Term doi:10.1016/j.agrformet.2008.01.017 Data from the Canadian Intersite Decomposition Experiment. Information Krishnan P, Black TA, Jassal RS, Chen B, Nesic Z (2009) Interannual variability of Report BC-X-422. Natural Resources Canada, Canadian Forest Service, Pacific the carbon balance of three different-aged Douglas-fir stands in the Pacific Forestry Centre, Victoria, BC Northwest. J Geophys Res 114:G04011, doi:10.1029/2008JG000912 Stage AR, Crookston NL (2007) Partitioning error components for accuracy- Kurz W, Apps M, Banfield E, Stinson G (2002) Forest carbon accounting at the assessment of near-neighbor methods of imputation. For Sci 53:62–72 operational scale. For Chron 78:672–679 Ter-Mikaelian M (1997) Biomass equations for sixty-five North American tree Kurz WA, Dymond CC, White TM, Stinson G, Shaw CH, Rampley GJ, Smyth C, species. For Ecol Manage 97:1–24, doi:10.1016/S0378-1127(97)00019-4 Simpson BN, Neilson ET, Trofymow JA, Metsaranta J, Apps MJ (2009) CBM- Theede AD (2007) Biometric and Eddy-Covariance Estimates of Ecosystem Carbon CFS3: a model of carbon-dynamics in forestry and land-use change Storage at Two Boreal Forest Stands in Saskatchewan: 1994–2004. University implementing IPCC standards. Ecol Modell 220:480–504, doi:10.1016/ of Saskatchewan, MSc thesis j.ecolmodel.2008.10.018 Trofymow JA, Stinson G, Kurz WA (2008) Derivation of a spatially explicit 86-year Lambert M-C, Ung C-H, Raulier F (2005) Canadian national tree aboveground retrospective carbon budget for a landscape undergoing conversion from biomass equations. Can J For Res 35:1996–2018, doi:10.1139/x05-112 old-growth to managed forests on Vancouver Island, BC. For Ecol Manage Landsberg JJ, Waring RH (1997) A generalised model of forest productivity using 256:1677–1691, doi:10.1016/j.foreco.2008.02.056 simplified concepts of radiation-use efficiency, carbon balance and partitioning. Wang Z, Grant RF, Arain MA, Chen BN, Coops N, Hember R, Kurz WA, Price DT, For Ecol Manage 95:209–228, doi:10.1016/S0378-1127(97)00026-1 Stinson G, Trofymow JA, Yeluripati J, Chen Z (2011) Evaluating weather effects Law B, Thornton P, Irvine J (2001) Carbon storage and fluxes in ponderosa pine on interannual variation in net ecosystem productivity of a coastal temperate forests at different developmental stages. Glob Chang Biol 7:755–777 forest landscape: a model intercomparison. Ecol Modell 222:3236–3249, Law BE, Turner D, Campbell J, Sun OJ, Van Tuyl S, Ritts WD, Cohen WB (2004) doi:10.1016/j.ecolmodel.2011.06.005 Disturbance and climate effects on carbon stocks and fluxes across Western Yashiro Y, Lee N-YM, Ohtsuka T, Shizu Y, Saitoh TM, Koizumi H (2010) Biometric- Oregon USA. Glob Chang Biol 10:1429–1444, doi:10.1111/j.1365- based estimation of net ecosystem production in a mature Japanese cedar 2486.2004.00822.x (Cryptomeria japonica) plantation beneath a flux tower. J Plant Res 123:463–472, Leclerc M, Thurtell G (1990) Footprint prediction of scalar fluxes using a doi:10.1007/s10265-010-0323-8 Markovian analysis. Boundary-Layer Meteorol 52(3):247–258 Li Z, Kurz WA, Apps MJ, Beukema SJ (2003) Belowground biomass dynamics in the Carbon Budget Model of the Canadian Forest Sector: recent improvements and implications for the estimation of NPP and NEP. Can J For Res 33:126–136, doi:10.1139/x02-165 Mahalanobis P (1936) On the generalized distance in statistics. Proc Natl Inst Sci 12:49–55 Matthews G (1993) The carbon content of trees. Forestry Commission Technical Submit your manuscript to a Paper, vol. 4. Forestry Commission, London journal and beneﬁ t from: Miller SD, Goulden ML, Menton MC, da Rocha HR, de Freitas HC, Figueira AMES, de Sousa CA D (2004) Biometric and micrometeorological measurements of 7 Convenient online submission tropical forest carbon balance. Ecol Appl 14:114–126, doi:10.1890/02-6005 7 Rigorous peer review Morgenstern K, Andrew Black T, Humphreys ER, Griffis TJ, Drewitt GB, Cai T, Nesic Z, Spittlehouse DL, Livingston NJ (2004) Sensitivity and uncertainty of the 7 Immediate publication on acceptance carbon balance of a Pacific Northwest Douglas-fir forest during an El Niño/La 7 Open access: articles freely available online Niña cycle. Agric For Meteorol 123:201–219, doi:10.1016/ 7 High visibility within the ﬁ eld j.agrformet.2003.12.003 7 Retaining the copyright to your article NFI (2008) Canada’s National Forest Inventory - Ground Sampling Guidelines. Natural Resources Canada, Canadian Forest Service, Victoria, BC, [Electronic Resource]. Available: https://nfi.nfis.org/documentation/ground_plot/Gp_ Submit your next manuscript at 7 springeropen.com guidelines_v5.0.pdf. Accessed 13 April 2015 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png "Forest Ecosystems" Springer Journals http://www.deepdyve.com/lp/springer-journals/comparison-of-carbon-stock-changes-eddy-covariance-carbon-fluxes-and-RGDe4JLzcy

Loading next page...

References (60)

M Ter-Mikaelian (1997)
Biomass equations for sixty-five North American tree species
For Ecol Manage, 97
C Rebmann M Göckede (2004)
A combination of quality assessment tools for eddy covariance measurements with footprint modelling for the characterisation of complex sites
Agric For Meteorol, 127
D Turner BE Law (2004)
Disturbance and climate effects on carbon stocks and fluxes across Western Oregon USA
Glob Chang Biol, 10
M Apps W Kurz (2002)
Forest carbon accounting at the operational scale
For Chron, 78
JK Brown (1971)
A planar intersect method for sampling fuel volume and surface area
For Sci, 17
Black TA Humphreys ER (2006)
Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of development after clearcut harvesting
CC Dymond WA Kurz (2009)
CBM-CFS3: a model of carbon-dynamics in forestry and land-use change implementing IPCC standards
Ecol Modell, 220
T Andrew Black K Morgenstern (2004)
Sensitivity and uncertainty of the carbon balance of a Pacific Northwest Douglas-fir forest during an El Niño/La Niña cycle
Agric For Meteorol, 123
citation_title=R (2011)
R: a language and environment for statistical computing
S Brown D Clark (2001)
Measuring net primary production in forests: concepts and field methods
Ecol Appl, 11
HH Janzen (2005)
Soil carbon: a measure of ecosystem response in a changing world?
Can J Soil Sci, 85
Jungen J (1985)
Soils of Southern Vancouver Island
G Stinson JA Trofymow (2008)
Derivation of a spatially explicit 86-year retrospective carbon budget for a landscape undergoing conversion from old-growth to managed forests on Vancouver Island, BC
For Ecol Manage, 256
M Jomura Y Kominami (2008)
Biometric and eddy-covariance-based estimates of carbon balance for a warm-temperate mixed forest in Japan
Agric For Meteorol, 148
P Mahalanobis (1936)
On the generalized distance in statistics
Proc Natl Inst Sci, 12
N Bréda A Granier (2008)
Ten years of fluxes and stand growth in a young beech forest at Hesse, North-eastern France
Ann For Sci, 65
G Matthews (1993)
The carbon content of trees. Forestry Commission Technical Paper, vol. 4
, 4
O Bouriaud F Babst (2014)
Above-ground woody carbon sequestration measured from tree rings is coherent with net ecosystem productivity at five eddy-covariance sites
New Phytol, 201
Chen JM Chen B (2008)
Comparison of regional carbon flux estimates from CO 2 concentration measurements and remote sensing based footprint integration
ME Harmon JE Janisch (2005)
Decomposition of coarse woody debris originating by clearcutting of an old-growth conifer forest
Ecoscience, 12
DY Hollinger AD Richardson (2007)
A method to estimate the additional uncertainty in gap-filled NEE resulting from long gaps in the CO2 flux record
Agric For Meteorol, 147
Trofymow JA Ferster CJ (2011)
Determination of ecosystem carbonstock distributions in the flux footprint of an eddy-covariance tower in a coastal forest in British Columbia
G Thurtell M Leclerc (1990)
Footprint prediction of scalar fluxes using a Markovian analysis
Boundary-Layer Meteorol, 52
H Schmid (2002)
Footprint modeling for vegetation atmosphere exchange studies: a review and perspective
Agric For Meteorol, 113
P Thornton B Law (2001)
Carbon storage and fluxes in ponderosa pine forests at different developmental stages
Glob Chang Biol, 7
NL Crookston AR Stage (2007)
Partitioning error components for accuracy-assessment of near-neighbor methods of imputation
For Sci, 53
W Mo T Ohtsuka (2007)
Biometric based carbon flux measurements and net ecosystem production (NEP) in a temperate deciduous broad-leaved forest beneath a flux tower
Ecosystems, 10
T Bolger K Black (2005)
Inventory and eddy covariance-based estimates of annual carbon sequestration in a Sitka spruce (Picea sitchensis (Bong.) Carr.) forest ecosystem
Eur J For Res, 126
RF Grant Z Wang (2011)
Evaluating weather effects on interannual variation in net ecosystem productivity of a coastal temperate forest landscape: a model intercomparison
Ecol Modell, 222
AD Theede (2007)
Biometric and Eddy-Covariance Estimates of Ecosystem Carbon Storage at Two Boreal Forest Stands in Saskatchewan: 1994–2004
E Ceschia A Granier (2000)
The carbon balance of a young Beech forest
Funct Ecol, 14
ML Goulden SD Miller (2004)
Biometric and micrometeorological measurements of tropical forest carbon balance
Ecol Appl, 14
F Gougeon (1995)
A crown-following approach to the automatic delineation of individual tree crowns in high spatial resolution aerial images
Can J Remote Sens, 21
WA Kurz Z Li (2003)
Belowground biomass dynamics in the Carbon Budget Model of the Canadian Forest Sector: recent improvements and implications for the estimation of NPP and NEP
Can J For Res, 33
AM Solomon RK Dixon (1994)
Carbon pools and flux of global forest ecosystems
Science, 263
G Matteucci (2000)
Interactions Between the Carbon and Nitrogen Cycles and the Role of Biodiversity: A Synopsis of a Study Along a North–south Transect Through Europe
TA Black RS Jassal (2010)
Soil CO2 and N2O flux dynamics in a nitrogen-fertilized Pacific Northwest Douglas-fir stand
Geoderma, 157
AO Finley NL Crookston (2008)
yaImpute: an R package for kNN imputation
J Stat Softw, 23
A Omule M Gillis (2005)
Monitoring Canada’s forests: the national forest inventory
For Chron, 81
TA Black P Krishnan (2009)
Interannual variability of the carbon balance of three different-aged Douglas-fir stands in the Pacific Northwest
J Geophys Res, 114
N-YM Lee Y Yashiro (2010)
Biometric-based estimation of net ecosystem production in a mature Japanese cedar (Cryptomeria japonica) plantation beneath a flux tower
J Plant Res, 123
P Jukka P Kolari (2004)
Carbon balance of different aged Scots pine forests in Southern Finland
Glob Chang Biol, 10
S Cooper (1989)
Concepts and techniques of vegetation mapping
RH Waring JJ Landsberg (1997)
A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning
For Ecol Manage, 95
citation_title=Canada’s National Forest Inventory - Ground Sampling Guidelines, citation_publication_date= (2008)
Canada’s National Forest Inventory - Ground Sampling Guidelines
DD Baldocchi (2008)
“Breathing” of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems
Aust J Bot, 56
CR Lloyd HP Schmid (1999)
Spatial representativeness and the location bias of flux footprints over inhomogeneous areas
Agric For Meteorol, 93
R Jassal (2008)
Carbon Sequestration in British Columbia’s Forests and Management Options
C Mayrsohn JD Bailey (1998)
Understory vegetation in old and young Douglas-fir forests of western Oregon
For Ecol Manage, 112
K Klinka (1994)
A Field Guide to Site Identification and Interpretation for the Vancouver Forest Region. Land Management Handbook Number 28
B Vos B Gielen (2013)
Biometric and eddy covariance-based assessment of decadal carbon sequestration of a temperate Scots pine forest
Agric For Meteorol, 174
C-H Ung M-C Lambert (2005)
Canadian national tree aboveground biomass equations
Can J For Res, 35
C Vogel C Gough (2008)
Multi-year convergence of biometric and meteorological estimates of forest carbon storage
Agric For Meteorol, 148
TA Black BB Chen (2009)
Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements
Boundary-Layer Meteorol, 130
TA Black RF Grant (2007)
Changes in net ecosystem productivity with forest age following clearcutting of a coastal Douglas-fir forest: testing a mathematical model with eddy covariance measurements along a forest chronosequence
Tree Physiol, 27
SC Wofsy CC Barford (2001)
Factors controlling long- and short-term sequestration of atmospheric CO2 in a mid-latitude forest
Science, 294
T Hilker NC Coops (2007)
Estimating canopy structure of Douglas-fir forest stands from discrete-return LiDAR
Trees-Structure Funct, 21
P Hanson P Curtis (2002)
Biometric and eddy-covariance based estimates of annual carbon storage in five eastern North American deciduous forests
Agric For Meteorol, 113
HP Schmid JL Ehman (2002)
An initial intercomparison of micrometeorological and ecological inventory estimates of carbon exchange in a mid-latitude deciduous forest
Glob Chang Biol, 8
W Kurz (2010)
Decreasing Uncertainty in CBM-CFS3 Estimates of Forest Soil Carbon Sources and Sinks Through use of Long-Term Data from the Canadian Intersite Decomposition Experiment. Information Report BC-X-422

Publisher: Springer Journals
Copyright: 2015 Ferster et al.; licensee Springer.
eISSN: 2197-5620
DOI: 10.1186/s40663-015-0038-3
Publisher site: See Article on Publisher Site

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. Keywords: Forest carbon; Micrometeorology; Biometry; Remote sensing; Geographic information systems * Correspondence: colin.ferster@ubc.ca Department of Forest Resources Management, University of British Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada Full list of author information is available at the end of the article © 2015 Ferster et al.; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. Ferster et al. Forest Ecosystems (2015) 2:13 Page 2 of 19 Background 2008), and by measuring forest C stocks at two times Forests are a large component of the global carbon (C) with a sufficient interval between measurements, the C stocks, containing an estimated 1146 Pg C (Dixon et al. stock change (ΔC) can be determined and used to esti- 1994). Forest processes, which may be influenced by forest mate the cumulative net uptake of C from the atmos- management, can therefore have a large impact on the phere to the forest for the period of measurement global C budget, either by storing or releasing C. It is (ΣNEP) (Clark and Brown 2001). Measurements are typ- therefore critical to understand forest C dynamics, includ- ically made of individual plants, for example, all trees ing forest C stock components and transfer mechanisms within a plot are measured for diameter, species, and in order to develop accurate forest C models such as the height, then allometric relationships are used to deter- Carbon Budget Model of the Canadian Forest Sector mine tree mass (Ter-Mikaelian 1997). Similarly, woody (CBM-CFS3) (Kurz et al. 2009), 3PG (Landsberg and debris pieces are measured along transects using line Waring 1997), and Ecosys (Grant et al. 2007), amongst intersect sampling and allometric relationships are ap- others, as well as to inform forest management policy and plied to find biomass (Brown 1971). Other components for national and international reporting (Kurz et al. 2002). are made by direct measurements, for example, under- Factors affecting forest C dynamics include natural (e.g. story vegetation can be sampled, dried, and weighed fire, insect outbreaks) and anthropogenic disturbances (Bailey et al. 1998). To measure soil C, field samples are from land use management and change (e.g. harvest, re- typically dried and weighed for calculation of bulk density, forestation, deforestation) (Kurz et al. 2002) as well as and sent for laboratory analysis of C using dry combustion weather (Morgenstern et al. 2004), fertilization (Jassal (Janzen 2005). Soil C measurements are important be- et al. 2010), and stand age and species composition, which cause the underground processes driving soil respiration can be related to disturbance history, site edaphic char- form a large component of ecosystem gas exchanges, but acteristics, or management practices (e.g. silvicuture) these are less well-understood and measured compared to (Humphreys et al. 2006; Krishnan et al. 2009). above ground stocks and processes (Ehman et al. 2002). To measure the net exchange of C between land eco- Changes in soil C storage over time have implications as a systems and the atmosphere (net ecosystem exchange, source or sink for the global C cycle, and as an indicator NEE, with -NEE referred to as net ecosystem productiv- of environmental function and health (Janzen 2005). ity, NEP), a global network of over 400 eddy-covariance One application of forest measurement data to provide (EC) flux stations has been established across a range of accounting of forest management actions and subse- ecosystems, building an extensive data record of NEP, in quently inform forest management decisions is as an in- some cases spanning up to two decades. The majority of put to C budget models such as CBM-CFS3. CBM-CFS3 these towers are located on sites not undergoing major is a forest C budget model that utilizes growth and yield disturbances so the fluxes measured reflect the inter- information for biomass and generates explicit simula- action of weather, vegetation composition, stand age, tion of dead organic matter dynamics (Kurz et al. 2009). and seasonal phenology. EC flux-towers use micro- In contrast, forest process models use approaches based meteorological equipment to take measurements at the on the understanding of processes such as photosyn- canopy scale of the exchanges of gasses including CO , thesis; these process-based models have more potential water vapor, sensible heat, and some are equipped to to model changing conditions because they do not rely measure other trace gases (Baldocchi 2008). The source on historic growth and yield and therefore are better area contributing to measurements made by the instru- suited to simulate forest conditions under global change ments mounted on the EC flux-tower, the flux-tower foot- situations (Landsberg and Waring 1997). CBM-CFS3 print, is variable in size and shape depending on height of was designed to meet reporting requirements for forest measurement, surface roughness length, wind speed and management, provide policy support, and a function as a C budgeting tool for operating foresters, thus function- direction, and atmospheric stability (Leclerc and Thurtell 1990; Schmid 2002), and proper interpretation of EC flux- ing at a range of spatial scales, and estimating records of tower based measurements depends on the flux-footprint C stocks, transfers between pools, and emissions (Kurz et al. 2002; Kurz et al. 2009). over which the fluxes are sampled (Chen et al. 2008). Early estimates modelled flux-footprints as simple ovals. Re- Comparing EC flux-tower cumulative NEP (ΣNEP) cently, estimates of flux-footprint climatology may dem- and inventory ground plot measurements of C stock changes (ΔC) over the same period can help inform for- onstrate more complex geometries and continuous probability density surfaces to quantify the upwind distri- est C processes for several reasons. First, ground plot bution of weighting factors over long time periods (Chen measurements can serve as an independent validation of ΣNEP measurements made at EC flux-towers. Second, et al. 2009). Forest measurements can be taken to quantify forest C inventory plot data can provide more detailed informa- stocks (Dixon et al. 1994; Clark and Brown 2001; NFI tion about the stand structural changes that may be Ferster et al. Forest Ecosystems (2015) 2:13 Page 3 of 19 driving variation stand-level EC C exchange. Third, broad determined and the ΔC stocks compared to CBM-CFS3 stand inventory measurement datasets are available for a modelled and EC-flux-tower estimates of ΣNEP at the greater spatial extent than the global EC flux-tower net- CCP coastal British Columbia Flux Station. Through a work. For example, the Canadian National Forest Inven- combination of variable topography and a complex dis- tory (NFI) samples approximately 22,000 photo plots turbance history (Trofymow et al. 2008), the EC flux- locations across Canada with detailed ground plots mea- tower sites possess fine-scale spatial variation in forest sured at over 1000 locations (Gillis et al. 2005). However, structure, thus presenting a challenge to the interpret- comparing inventory plot changes and EC flux-tower ation of EC flux data (Schmid and Lloyd 1999; Göckede measurements requires comparisons to be made across et al. 2004). To allow the comparison between measure- different spatial and temporal scales. Therefore, such com- ments from the inventory plots, C budget models, and parisons pose methodological challenges beyond individ- EC flux-towers, our approach consisted of four steps: ual inventory plot and EC flux-tower measurements. first, we found the change in C based on forest inventory Comparisons between EC tower measurements of ground plot measurements made in 2002 and 2006 (four ΣNEP and inventory plot measurements of ΔC stocks year period); second, we developed an additional ap- have been completed at a number of sites globally proach that utilized litterfall data and published decay (Schulze et al. 2000; Granier et al. 2000; Law et al. 2001; equations, for a second inventory-based estimate of C Barford et al. 2001; Ehman et al. 2002; Curtis et al. 2002; stocks change; third, we defined stand attributes at a fine Miller et al. 2004; Black et al. 2005; Gough et al. 2008; spatial resolution, stratified the site based on stand attri- Kominami et al. 2008; Yashiro et al. 2010; Gielen et al. butes, and weighted the two inventory plot estimates; 2013; Babst et al. 2014), and at two sites of the Canadian and fourth, calculated ΔC. These were compared Carbon Program (CCP), Boreal Ecosystems Research and with CBM-CFS3 model-based estimates of ΣNEP and Monitoring Sites (BERMS) station, located in the boreal EC flux-tower-measured ΣNEP for four years (2003, forest of northern Saskatchewan (Theede 2007). Over a 2004, 2005 and 2006) in the same period. Finally, we ten-year interval (1994 – 2004), Theede (2007) found a evaluated how the forest age and conditions at the dif- −1 −1 convergence of 15.6 ± 4.0 MgC ha 10 years (inventory ferent sites accounted for convergence or lack of conver- −1 −1 plots) and 18.2 ± 8.09 MgC ha 10 years (EC tower) at gence in estimates for the different methods, and −1 −1 the Old Aspen site and 5.8 ± 2.0 MgC ha 10 years (in- discussed the relative advantages and constraints of each −1 −1 ventory plots) and 6.9 ± 1.6 MgC ha 10 years (EC measurement method. tower) at the Old Jack Pine site. However, the homoge- neous stand structure contributing to canopy-level EC Methods flux measurements and large fetch in the boreal forest re- Study area duced the need for consideration of spatial vegetation The coastal British Columbia Flux Station sites are located structure and footprint distributions. on the leeward central east coast of Vancouver Island These previous studies have contributed considerably to (Figure 1), in the Coastal Western Hemlock biogeocli- the understanding of forest C dynamics; however, at sites matic zone (Green and Klinka 1994) with a mean annual with more complex vegetation structure and topography, air temperature of 10° Celsius. Douglas-fir (Pseudotsuga such as the forests found in coastal British Columbia, the menziesii var menziesii (Mirb.) Franco) is the dominant spatial distribution of vegetation structure and the foot- species, with lesser amounts of western hemlock (Tsuga print distribution are important considerations for accur- heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata ate comparisons of EC flux-tower and inventory plot Donn ex D. Donn), and red alder (Alnus rubra Bong.). measurements (Schmid and Lloyd 1999). More complex Three eddy-covariance and meteorological towers were ecosystems, therefore, require more complex models and placed in stands of different ages (Table 1). methods developed for accurate comparisons of EC flux- EC flux-tower footprints, which delineate the spatial tower measurements and inventory-based measurements. distribution of upwind source-weighting factors, were For example, a limited number of studies, including calculated by Chen et al. (2009) to represent the size, Ehman et al. (2002), Chen et al. (2009), Ferster et al. shape, direction, and magnitude of the flux as a function (2011), and Gielen et al. (2013) applied weighting by EC of wind speed and direction, measurement height, and flux-tower footprint climatology to inventory-based mea- surface roughness. Chen et al. (2009) calculated hourly surements, where areas of the footprint that contribute footprints for 10 m by 10 m cells, weighted by NEP, and more to tower-measured flux are weighted more heavily averaged over the 2002 – 2006 measurement period to than other areas, resulting in an overall improvement determine the footprint climatology for the sites. In this when compared to simple flux footprint geometries. study, the flux-footprint climatology 85% cumulative flux In this paper, the difference between 2002 and 2006 probability boundary was selected as the maximum ex- inventory plot measurements of ecosystem C stocks was tent of analysis, since the majority of the probability Ferster et al. Forest Ecosystems (2015) 2:13 Page 4 of 19 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. density is concentrated within this area and the remain- Live biomass der extends across a large area (Figure 2). Live trees greater than 1.3 m tall were measured in an 11.28 m radius circular plot. Where trees were very nu- merous (for example, at HDF88), a half plot or 5 m ra- Inventory ground plots dius circular plot was measured following NFI guidelines All ground plots were established and measured by the (NFI 2008). Allometric equations are used in the NFI Canadian Forest Service at the four sites following compiler to estimate stem, bark, branch, and foliar bio- Canadian National Forest Inventory (NFI) guidelines mass based on work by Lambert et al. (2005). These bio- and protocols (NFI 2008) though they are not part of mass values were converted to C content by assuming the primary NFI ground plot network. Plot locations 52, 56, 52, and 52% dry C concentration, respectively were chosen to represent each site series class (based (Matthews 1993). on air photo interpretation and field transects) within Understory vegetation from four 1 m × 1 m micro- the preliminary flux footprint boundary estimated in plots was destructively sampled and sent to the labora- 2002 (Figure 2). Establishment and initial measurements tory for drying and weighing for determination of mass were made in September 2002 and re-measurements were which was converted to C content assuming 50% dry C made in September 2006, constituting a four-year meas- concentration. urement interval. NFI data-compilation software used to determine plot-level values through application of allo- metric equations to overstory trees and density values for Detritus woody debris NFI (2008). Since live overstory vegetation Woody debris, fallen woody material > 1 cm diameter, roots were not measured in 2002 or 2006, coarse and fine was measured for diameter, species, and decay class along root C mass was estimated from overstory biomass using four 15 m transects at each plot. Decay class was deter- equations by Li et al. (2003). mined by field crews in the five-class ordinal rating system Ferster et al. Forest Ecosystems (2015) 2:13 Page 5 of 19 Table 1 Study site locations and characteristics mass. C content was determined assuming 50% dry C Site Location (latitude Site description concentration of woody debris mass. Fine woody debris, and longitude) fallen woody material 2 mm to <1 cm diameter, was sam- DF49 49.868797°, � Near-rotation stand established pled from within the four 1 × 1 m understory microplots in 1949. in each sample plot and sent for laboratory drying, weigh- (Figure 2a) - 125.33515° � 162 ha flux tower footprint. ing, and determination of mass, and C content assuming � 12 ground plots. dry mass is composed of 50% C. Dead standing trees were measured for diameter, height, and species and decay class � 260 to 470 m elevation. was observed and recorded by field crews. Stumps were � Harvested of old growth timber measured for diameter and decay class in 2002 but were in 1937, 1938, and 1943. not remeasured in 2006. � Broadcast burned in 1938 and 1943 following harvest. Surface substrates were measured along four 15 m long transects to determine the average depth and per- � Planted in 1949. cent area surface coverage in the entire plot. The organic HDF88 49.536655°, � Juvenile pole-sapling stand litter, fibric, and humus layer (forest floor) was sampled established in 1988. and average depths measured from within four 20 × 20 (Figure 2b) - 124.90146° � 35 ha flux tower footprint. cm templates (one at each microplot). These destructive � 6 ground plots. samples were collected, sieved, dried, and weighed to de- � 150 to 220 m elevation. termine bulk density and subsamples sent for laboratory � Harvested of second-growth analysis of C concentration. timber in 1987. � Broadcast burned following harvest. Mineral soil � Planted in 1988. In 2002, <2 mm mineral soil was sampled for bulk dens- HDF90 49.893666°, � Juvenile pole-sapling stand ity and C concentrations from 10 – 12 cm diameter ex- established in 1990. cavated holes at three depth intervals 0–15 cm (4 - 125.304415° � No flux tower. samples), 15–35 cm (2 samples), and 35–55 cm (one � 6 ground plots. sample) and % volumetric coarse fragments determined from one 55-cm-deep soil pit per ground plot. Data were � 175 m elevation. scaled to the entire groundplot discounting for the area � Harvested of second-growth timber in 1990. without mineral soil (i.e. exposed bedrock). In 2006, the 0–15 cm layer was re-measured for C content at four lo- � Broadcast burned following harvest. cations in each plot using a 2 cm diameter soil corer. Samples were stored in plastic bags at 2°C prior to la- � Planted in 1990. boratory processing. For a detailed description of soil � Similar in age and composition to HDF88. sampling methods see NFI (2008). For 2002–2006 soil C stock changes, only the 0–15 cm layer was considered, HDF00 49.872177°, � Regenerating clearcut stand established in 2000. soil C at deeper depths was assumed unchanged. (Figure 2c) - 125.29235° � 14 ha flux tower footprint. � 9 ground plots. Changes in inventory ground plot C stocks � 160 to 190 m elevation. Changes in C stocks over the four-year interval were cal- � Harvested from a second-growth culated from the difference between C stocks in 2002 stand in 2000. and 2006. The total change in C stocks was calculated � Logging debris were piled and as: burned. � Planted in 2000. ΔC ¼ ΔC þ ΔC þ ΔC L D S where ΔC is the change in total live biomass C stocks in- utilized by the NFI, from intact (Class 1) to highly decom- cluding live trees, and shrubs, herbs, and bryophytes; ΔC posed (Class 5) (NFI 2008). NFI compilation routines is the total change in detrital C stocks including dead utilize algorithms for line intercept sampling to determine standing trees, woody debris, and the forest floor; and plot level volume and a lookup table of woody debris ΔC is the change in mineral soil C stocks (0–15 cm). density values by species and decay class to determine S Ferster et al. Forest Ecosystems (2015) 2:13 Page 6 of 19 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. ΔC Detritus annual accounting method between measurement dates, and assumed to turn over D2 Preliminary examination of the ΔC values showed at a rate of 80% annually. Shrubby litterfall material larger-than-expected changes in the measured C stocks was divided into stem and branch and foliage compo- over the four-year interval, especially for the forest nents, and assumed to turn over at a rate of 3% and floor component. To evaluate these measured changes, 95%, respectively similar to the fine branch wood and ΔC was estimated using a second method (referred to foliage for deciduous trees (Kurz et al. 2009). When as ΔC ) that accounted for the annual inputs into the trees died within the measurement interval, the bio- D2 detrital pools from litterfall and mortality (i.e. the trans- mass was transferred to the detrital pools half way fers from live biomass to forest floor and soil , and live through the measurement interval. Finally, the biomass trees to dead trees respectively) and annual losses and of any trees (live or dead) that were measured in 2002 transfers from decomposition estimated using pub- and not recorded by the field crew in 2006 was trans- lished equations and parameters (Smyth et al. 2010; ferred to the detrital pools half-way through the meas- Kurz et al. 2009). urement interval under the assumption that the tree Inputs to detrital C pools from 2002 – 2005 included fell during the measurement interval. overstory fine litterfall (needles, leaves, cones, twigs) Losses through decomposition and transfer to the soil which was collected quarterly in 3 0.189 m conical- organic C of each detrital pool was calculated at an an- mesh litterfall traps located in each ground plot at nual time step with an annual average temperature of DF49, HDF88 and HDF90. Annual litter fall masses 10°C using coefficients that were developed and vali- were calculated and converted to C assuming 50% dry dated nationally (Smyth et al. 2010) for the CBM-CFS3 C content. Herbaceous material, which was not tall model (Kurz et al. 2009). The transfers and coefficients enough to be sampled using litterfall traps especially at used are presented in Figure 3. Soil C stock changes the HDF00 site, was collected and measured at the were assumed negligible, and ΔC was calculated as: sample plots in 2002 and 2006, linearly interpolated ΔC ΔC + ΔC 2= L D2. Ferster et al. Forest Ecosystems (2015) 2:13 Page 7 of 19 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 using R version 2.11 (R Development Core Team 2011) To account for spatial heterogeneity due to stand struc- and yaImpute version 1.0-10 (Crookston and Finley ture and topography, the methodology developed by 2008). Predictor variables were selected if correlations Ferster et al. (2011) was applied to the ground plot data at with inventory-based C stocks in 2002 and 2006 were the flux tower sites. Following this approach, predictor significant at the 95% confidence level and there was variables from GIS, topography, and remote sensing data no significant co-linearity with other predictor vari- were used to impute measurements from the inventory ables. Following the calculation, the 2002 – 2006 cumu- ground plots across the forest site based on the three- lative NEP footprint probability density surface for each most-similar-neighbours (K-MSN with K = 3) (Crookston site (calculated by Chen et al. 2008 and Chen et al. and Finley 2008) using the Mahalanobis distance as a 2009) was used to weight each 10 m by 10 m footprint measure of similarity (Mahalanobis 1936). For each cell for the calculation of the site level mean value. footprint cell, the detailed inventory target measure- Site-level estimates of ΔCand ΔC stocks and were ments were estimated as the weighted mean of the calculated three ways. First, inventory- based site esti- three nearest neighbours (inversely weighted by the mates ΔC and ΔC stocks were calculated as an arith- Mahalanobis distance) and used for calculation of site- metic mean of ground plots at each site. Second, ΔC level means of ΔC. To assess the error of the model, stocks were calculated with the K-MSN prediction. the root mean squared difference (RMSD) was calcu- Third, and finally, ΔC stocks were calculated with the lated based on leave-one-out cross validation (Stage K-MSN prediction weighted by flux footprint probability and Crookston 2007). The procedure was completed density. Ferster et al. Forest Ecosystems (2015) 2:13 Page 8 of 19 Flux tower ∑NEP Comparing convergence of ΔC, EC tower ΣNEP, and CBM- Measurement and calculations of flux tower annual CFS3 ΣNEP NEP published by Black et al. (2008) and Krishnan Comparisons were made among the estimates by first et al. (2009) were summed to find the cumulative evaluating how each method ranked the stand ages. Sec- ∑NEP (2004–2006) used for this study. These values ond, estimates of EC tower ΣNEP and CBM-CFS3 ΣNEP were gap filled for conditions at night when friction vel- were compared. Third, at each stand age from youngest ocity was low (i.e. inadequate turbulent mixing) and to oldest, the estimates of ΔC and tower and CBM-CFS3 corrected for energy-balance closure. Annual NEP ΣNEP were compared for convergence by comparing values for 2003, 2004, 2005, and 2006 were summed for the means, variance around the means indicated by the each site to estimate the ΣNEP total for the four-year standard deviation, and spatial variance indicated by the period. Morgenstern et al. (2004) reported that uncer- RMSD. Comparisons of convergence among methods tainty in the annual NEP measured using EC at the were made for the following estimates: ΔC as an arith- −1 near-rotation stand may be as much as 0.9 MgC ha metic mean of plots (ΔC unweighted), ΔC with K- 2 2 −1 year (due to systematic error in the EC measure- MSN classification (ΔC K-MSN), ΔC with K-MSN 2 2 ments). This suggests that the uncertainty in the esti- classification and footprint weighting (ΔC K-MSN −1 mate at the near-rotation stand may be ±3.6 MgC ha FPW), CBM-CFS3 ΣNEP as an arithmetic mean of cells over the full four-year measurement interval. Detailed within the tower footprint, CBM-CFS3 ΣNEP as a foot- estimates of the uncertainty that propagates through print weighted mean (ΣNEP CBM-CFS3 FPW), and EC the calculation based on the friction velocity thresholds tower ΣNEP. (e.g. following Richardson and Hollinger 2007) was be- yond the scope of this paper. Results Live biomass ΔC For all sites, there was a positive change in the live C CBM-CFS3 ∑NEP mass of overstory trees (Table 2). For the near-rotation CBM-CFS3 is a forest C accounting model used for na- stand and juvenile stands, the increase was larger than at tional reporting of annual C inventories for Canada’s the regenerating site. The increase in live biomass C was managed forests (Kurz et al. 2009). This model uses nearly equal at the near-rotation and juvenile stands. An growth and yield equations and allometric equations to increase in shrub, herb, and bryophyte understory ΔC estimate tree growth and stand net primary production was observed at the near rotation and regenerating (NPP). The C mass of overstory trees and overstory tree stands, while a small decrease in understory ΔC stocks roots is use to estimate litterfall and root mortality occured at the juvenile stands, possibly related to an ex- transfers to the aboveground and belowground detritus pected increase in canopy closure as the overstory trees C pools, respectively. A soil submodel estimates de- matured. Comparing the two juvenile stands (HDF88 composition and heterotrophic respiration (Rh) based and HDF90), HDF88 had higher live biomass C than on the size of various detrital pools and mean annual HDF90 (for example, live stem C and total live C were temperature. The model tracks all major C pools to en- more than 1 standard deviation larger). sure 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 Detritus ΔC D1 Oyster River and encompassing the DF49, HDF90, and Dead standing tree C decreased in the near-rotation HDF00sites (Figure1). Modelrunsfor theareawere stand due to dead standing trees falling during the meas- performed on 1 hectare grid cells of forest disturbance urement interval, and there was a small increase at the history data, forest cover data, growth and yield equations, other sites due to tree mortality (Table 2). Large woody and disturbance transition matrices from Trofymow debris increased in the near-rotation stand, and de- et al. (2008). The modeled ∑NEP values for the DF49, creased at the other sites. There was considerable vari- HDF00 sites (from Wang et al. 2011) were calculated ability among plots, indicated by the high standard as site-level means unweighted and weighted by the deviations compared to the means. Fine woody debris flux probability distribution for the model cells in the decreased at all sites, except one of the juvenile sites footprint area; and at HDF90 as an unweighted mean (where there was a small increase). The large decrease in of model cells over the spatial extent of the ground fine woody debris at the regenerating clearcut was likely plots. Annual model values of NEP for 2003, 2004, due to decomposition of debris from the recent harvest. 2005 and 2006 were summed to calculate the ∑NEP Finally, forest floor material increased at all sites, with and mean annual NEP for the four year period the smallest increase at the regenerating clearcut, and (Figure 4). the largest increase at the juvenile stands. Ferster et al. Forest Ecosystems (2015) 2:13 Page 9 of 19 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. Mineral soil ΔC Results from the annual accounting method for de- Soil C (0–15 cm) stocks was highest at the near-rotation tritus, showed a net gain for the sum of components at site, followed by the two juvenile stands, and lowest at the near-rotation stand, and a decrease at the other sites the regenerating clearcut (Table 2). Over the measure- (Table 3). The increase in C at the near-rotation stand ment interval, large increases in mineral soil (0–15 cm) was slightly larger using ΔC than the direct measure- D2 ΔC were measured at all sites, with the highest at one ment method; however, the biggest difference was at the of the juvenile stands (HDF88). juvenile stand, which showed a net loss using ΔC , and D2 a large positive balance using the direct measurement Detritus annual accounting method ΔC method. For stumps, this method indicated that decom- D2 Measured litterfall was highest at the near-rotation position may have been greater than was initially expected stand, followed by the juvenile stands, and the regener- (stumps were not re-measured due to an expectation of ating stand had the lowest amount (Figure 5). The ma- minimal decomposition), especially at the more recently jority of litterfall at the near rotation stand and juvenile disturbed juvenile and regenerating clearcut sites. Large stands was needles followed by twigs. At the near- woody debris increased at the near-rotation stand, and de- rotation stand, litterfall was highest in 2001–2002, de- creased at the other sites, due to less dead standing trees creased through 2004–2005, and slightly increased in to fall and less branchfall. The increase in fine woody deb- 2005–2006. The other sites were relatively constant ris was estimated to be much lower using ΔC ,with the D2 through time. near-rotation stand having a net-accumulation of fine Ferster et al. Forest Ecosystems (2015) 2:13 Page 10 of 19 Table 2 Inventory-based site C stocks and ΔC estimates using arithmetic plot means (standard deviation) DF49 HDF90 HDF88 HDF00 Live Foliage C 2002 12.7 (1.8) 2.7 (0.9) 2.6 (0.8) 0 (0.1) Foliage ΔC −0.1 (0.7) 1.2 (0.7) 1.5 (0.4) 0.4 (0.5) Branch C 2002 22.1 (4) 2.5 (1) 1.5 (0.4) 0 (0.1) Branch ΔC 0.7 (1.2) 1.3 (0.9) 1.8 (0.5) 0.1 (0.2) Bark C 2002 17.6 (3.8) 0.7 (0.1) 1.1 (0.4) 0 (0) Bark ΔC 0.9 (0.9) 0.7 (0.3) 1 (0.3) 0.1 (0.2) Stem C 2002 90.6 (23.4) 2.4 (0.4) 4.7 (2) 0.1 (0.3) Stem ΔC 5.7 (4.7) 3.3 (1.1) 4.4 (1.3) 0.6 (1.1) Roots C 2002 32.6 (8.3) 1.8 (0.3) 2.8 (1.3) 0.1 (0.1) Roots ΔC 1.6 (1.6) 1.8 (0.6) 3.9 (0.9) 0.7 (1.5) Understory C 2002 1.3 (0.6) 8.3 (1.9) 11.2 (4.7) 4.9 (1.8) Understory ΔC 1.3 (1.8) −0.2 (2.8) −2.6 (4.5) 1.7 (1.9) Total Live C 2002 177 (39.9) 18.5 (3.7) 24 (4.7) 5.2 (2) Total Live ΔC 10.2 (8.7) 8.1 (6) 10 (5.9) 3.7 (4.4) Detritus Standing Dead C 2002 21.8 (13.9) 0.2 (0.2) 0 (0.1) 0 (0.2) Standing Dead ΔC −3.4 (11.1) 0.5 (0.6) 0.1 (0.2) 0 (0.2) Stumps C 2002 6.4 (3.2) 5.7 (2.8) 11.7 (6.1) 11.9 (10.6) Stumps ΔC 0 (0) 0 (0) 0 (0) 0 (0) Large Woody Debris C 2002 34.9 (20.5) 12.3 (8.7) 43.3 (25.3) 14.4 (11.1) Large Woody Debris ΔC 1.6 (12.4) −3.5 (9.7) −1.5 (18.8) −0.9 (6.3) Fine Woody Debris C 2002 3.1 (1) 1.3 (0.7) 0.5 (0.6) 6.1 (2.9) Fine Woody Debris ΔC −0.9 (1) −0.7 (0.8) 0.8 (0.8) −4.5 (2.6) Forest Floor C 2002 13.8 (7.3) 12.1 (7.2) 18.8 (13.1) 18.2 (11.6) Forest Floor ΔC 11 (6.5) 1.9 (11.2) 24.8 (19.7) 0.1 (12.3) Total Detritus C 2002 79.9 (32.7) 31.6 (8.4) 74.4 (38.7) 50.6 (16.8) Total Detritus ΔC 8.2 (21) −1.9 (8.1) 24.2 (30.3) −5.4 (12.6) Soil Mineral Soil 0–15 cm C 2002 38.6 (21.9) 29 (8.9) 34.6 (16.1) 23.6 (11.3) Mineral Soil 0–15 cm ΔC 18.8 (24.1) 27.8 (27.7) 10.7 (32.4) 24.2 (15.2) Totals Total Live ΔC 18.3 (23.4) 6.2 (10.5) 34.2 (30.6) −1.7 (14.9) + Total Detritus ΔC Total Live ΔC 37.2 (34.7) 34 (29.8) 44.9 (41.3) 22.5 (14.2) + Total Detritus ΔC + Mineral Soil 0 - 15 cm ΔC −1 All units are Mg C ha . woody debris, the juvenile sites had a small loss, and the layer, the near-rotation stand was estimated to have a large regeneration site experienced the largest decrease. This in- increase, while the other sites decreased, with the greatest dicates that, given the available inputs to litterfall and ex- decrease at the regenerating clearcut. pected rates of decomposition for fine woody debris, the Since the plot values for ΔC were judged more reli- D2 direct measured stock changes in fine woody debris was able, all subsequent analyses to spatially scale the plot very likely unrealistically high. Finally, for the forest floor data to the site were made using ΔC . 2 Ferster et al. Forest Ecosystems (2015) 2:13 Page 11 of 19 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. Spatial distribution of C stocks and ΔC stocks sites the other spatial predictor variables showed finer Predictor variables (Table 4) were selected and used for scale variation. At the regenerating clearcut, NDVI was theK-MSN procedureatthe siteswithfluxtowers a significant predictor for woody debris, forest floor, (Figure6). Slopewas an importantpredictor forall sites and mineral soil C stocks. as it was strongly correlated with soil C stocks (with The Mahalanobis distance (Figure 7) demonstrates less sloping plots having higher soil (0 – 15 cm) C how comprehensively the inventory ground plots repre- stocks), and at the regenerating clearcut, slope was also sent stand conditions across the footprint. Lower per- correlated with living tree C stocks. At the near rota- centiles represent small distances that were relatively tion site, forest-inventory mapping was a significant well represented by the ground plots. Areas with the predictor variable for overstory C stocks. Due to the smallest Mahalanobis distances and best ground plot smaller footprints at the other sites, forest-inventory representation were close to the towers in the areas con- mapping varied little within the site; however, at all tributing most to NEP measurements (within the 50% Table 3 Inventory-based site detrital C in 2002 and ΔC (2002–2006) using plot mean (standard deviations) litterfall, D2 mortality and decomposition values DF49 HDF90 HDF88 HDF00 Detritus C Standing Dead 21.8 (13.9) 0.2 (0.2) 0 (0.1) 0 (0.2) ΔC Standing Dead −3.4 (11.1) 0.5 (0.6) 0.1 (0.2) 0 (0.2) D2 C Stumps 6.4 (3.2) 5.7 (2.8) 11.7 (6.1) 11.9 (10.6) ΔC Stumps ΔC −0.5 (0.2) −0.4 (0.2) −0.9 (0.4) −0.9 (0.8) D2 C Large Woody Debris C 34.9 (20.5) 12.3 (8.7) 43.3 (25.3) 14.4 (11.1) ΔC Large Woody Debris 3.9 (9.6) −1 (1) −4.3 (2.6) −1.4 (1) D2 C Fine Woody Debris 3.1 (1) 1.3 (0.7) 0.5 (0.6) 6.1 (2.9) ΔC Fine Woody Debris 1 (0.8) −0.4 (0.3) −0.1 (0.2) −2.3 (1.1) D2 C Forest Floor 13.8 (7.3) 12.1 (7.2) 18.8 (13.1) 18.2 (11.6) ΔC Forest Floor 8.2 (4.1) −4.8 (4.2) −7.4 (7.2) −9.3 (6.7) D2 C Total Detritus 79.9 (32.7) 31.6 (8.4) 74.4 (38.7) 50.6 (16.8) ΔC Total Detritus 9.2 (8.7) −6.1 (4.2) −12.5 (9.4) −13.8 (6.5) D2 −1 All units are Mg C ha . Ferster et al. Forest Ecosystems (2015) 2:13 Page 12 of 19 Table 4 Environmental predictor variables used to determine stratification units at DF49, HDF88, and HDF00 Selection Coverage Variable Units DF49 HDF88 HDF00 Forest-Inventory Site Index m selected Top Height m Disturbance History (Trofymow et al. 2008) 1st Harvest year 2nd Harvest year 1st Fire year 2nd Fire year Date Est. 1953 year Date Est. 2003 year 1st Fertilization year 2nd Fertilization year Fire Cause 1 nominal Fire Cause 2 nominal Topography aspect azimuth selected selected selected Elevation m asl selected Slope degrees selected selected selected SCOSA selected selected SSINA selected selected TSRAI selected selected selected −1 1999 Orthophoto Dominant Canopy Tree Density stems ha selected selected (Gougeon 1995) 2004 Multispectral NDVI NDVI selected selected Forest-Inventory Cover Species (Trofymow et al. 2008) Site Species Soil Survey of Canada (Jungen 1985) Most Common Soil Association CFS and Forest Companies Site Series selected selected (Trofymow et al. 2008) Site Index: Tree height at 50 years age at breast height (1.3 m). Topographic Variables: SCOSA = Slope × cos (Aspect), SSINA = Slope × sin (Aspect); TSRAI (Topographic Solar Radiation Aspect Index) = (1-cos((π/180°)(Aspect-30°)))/2 (Roberts and Cooper 1989). Topography from 2004 LiDAR survey at DF49 and HDF00 (Coops et al. 2007). Topography from 1:50,000 National Topographic Series map at HDF88. Multispectral data from 2004 Quickbird survey at DF49 and HDF00. Multispectral data from 2004 Landsat scene at HDF88. cumulative flux probability density isoline). Notable For the juvenile stands, K-MSN and K-MSN FPW ΔC areas with large Mahalanobis distances included areas was less than 1 standard deviation higher than the un- with different forest cover types (e.g., patches of hard- weighted mean. The K-MSN estimate was higher than woods at the near-rotation stand), very few overstory the K-MSN FPW estimate. For ΔC , K-MSN and K- D2 trees (at the juvenile stand), and short steep slopes (be- MSN FPW estimates were slightly higher than the mean tween bench terraces at the regenerating clearcut). (less than 1 standard deviation). For the total ΔC , the For site level estimations of ΔC at the regenerating estimates were less than 1 standard deviation different, clearcut, K-MSN was more than 1 standard deviation and the standard deviations were very large compared to lower than the unweighted mean, and K-MSN-FPW was the means, indicating a large amount of variability at the 1 standard deviation lower than K-MSN. For ΔC K- site. D2, MSN and K-MSN FPW were less than 1 standard devi- At the near-rotation stand, ΔC K-MSN and K-MSN ation higher than the unweighted mean (Tables 3 and 5). FPW were slightly higher than the arithmetic mean (less The totals for ΔC K-MSN and K-MSN FPW were less than 1 standard deviation). The estimates using KMSN than 1 standard deviation higher than the arithmetic and KMSN-FPW were identical. For ΔC , the un- D2 mean (Table 6). weighted mean was higher than K-MSN and K-MSN Ferster et al. Forest Ecosystems (2015) 2:13 Page 13 of 19 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. FPW, but less than 1 standard deviation different. For regenerating clearcut as a source (net C release). While ΔC total, all estimates were within 1 standard deviation. the ranking was consistent for all methods, there were Comparing the ΔC and ΣNEP methods at the juvenile differences in the magnitude of values for each stand stands, HDF88 was compared with measurements from (Table 6). the flux tower installed at the site, and HDF90 was com- At the regenerating clearcut, the CBM-CFS3 estimate pared with CBM-CFS3 model output. While the two of ΣNEP with footprint weighting converged more stands were similar in seral stage, species composition, closely with the EC tower estimate of ΣNEP than the and stand history, there were differences in stand com- unweighted estimate, and both were within 1 standard position. For example, HDF88 had higher live C stocks, deviation. Evaluating the yearly means (Figure 8a), dem- in particular, live stem C, understory stem C, and total onstrated that footprint weighting improved conver- stem C. In addition, HDF88 had larger detritus C stocks gence between the CBM-CFS3 estimate of ΣNEP and including stumps C, large woody debris C, and total de- the EC flux ΣNEP. This was primarily due to less tritus C. The stand variability at HDF88 was also greater weighting for patches of uncut forest in the periphery of than at HDF90 indicated by the large standard devia- the flux footprint (Figure 4). tions compared to the means. Given the large amount of At the juvenile sites, CBM-CFS3 ΣNEP indicated that spatial variability at HDF88, spatial scaling, and footprint HDF90 was a weak sink (net C uptake) and EC tower weighting had a notable effect on the amount and sign ΣNEP indicated that HDF88 was a weak source (net C at this site. release) (Tables 3 and 6). At the near rotation stand, the estimates of ΣNEP from CBM-CFS3 were higher than Comparing inventory ΔC , EC tower ΣNEP and CBM-CFS3 the EC flux tower ΣNEP, with the footprint weighted es- ΣNEP timate being the highest. Considering the means of All methods ranked the near-rotation stand as a sink for ΣNEP from CBM-CFS3 on a yearly basis, footprint C (net C uptake), the juvenile stands as a weak sink or weighting decreased convergence with the EC flux tower weak source (small C uptake or release), and the measurement of ΣNEP (Figure 8b). Ferster et al. Forest Ecosystems (2015) 2:13 Page 14 of 19 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. Table 5 Site mean (± RMSD) changes in live biomass (ΔC and detrital (ΔC ) C stocks calculated using three most L) D2 similar neighbours classification (K-MSN), and K-MSN with footprint weighting means (K-MSN FPW) of plots for sites with flux-towers DF49 HDF88 HDF00 Pool Component K-MSN K-MSN FPW K-MSN K-MSN FPW K-MSN K-MSN FPW Live ΔC Foliage 0.1 ± 1.7 0.2 (0.4) ± 1.7 1.3 ± 2 1.3 (0.3) ± 2 0.2 ± 0.9 0.1 ± 0.9 ΔC Branches 0.9 ± 1.6 1 (0.6) ± 1.6 1.9 ± 1.3 1.8 (0.6) ± 1.3 0.1 ± 1.1 0 ± 1.1 ΔC Bark 1.1 ± 1.5 1.1 (0.4) ± 1.5 1 ± 1.7 1 (0.3) ± 1.7 0.1 ± 0.8 0 ± 0.8 ΔC Stem 6.7 ± 1.4 6.1 (2.2) ± 1.4 4.4 ± 1.7 4.2 (1.2) ± 1.7 0.4 ± 0.9 0.1 ± 0.9 ΔC Understory 1.8 ± 1 2.5 (2.3) ± 1 −0.9 ± 0.8 −1.7 (5.2) ± 0.8 0.8 ± 1.3 0.6 ± 1.3 ΔC Roots 2 ± 1.5 1.9 (0.8) ± 1.5 4.3 ± 2 4.4 (0.6) ± 2 0.3 ± 1.1 0.1 ± 1.1 ΔC 12.8 ± 1.6 12.8 (3) ± 1.6 11.9 ± 0.8 10.9 (7.2) ± 0.8 1.8 ± 0.8 1 ± 0.8 LTotal Detritus ΔC Standing Dead −5.7 ± 1.2 −8.1 (16.1) ± 1.2 0.1 ± 1.3 −0.7 (0.4) ± 1.2 0 ± 1.4 0 ± 1.4 D2 ΔC Stumps −0.5 ± 1.1 −0.5 (0.1) ± 1.1 −0.6 ± 1.2 −3.5 (1.1) ± 1.9 −0.7 ± 1.3 −0.8 ± 1.3 D2 ΔC Large Woody Debris 3.5 ± 0.9 6.3 (13.7) ± 0.9 −3.5 ± 1.9 −0.1 (0.2) ± 1.1 −1.5 ± 0.3 −1.5 ± 0.3 D2 ΔC Find Woody Debris 0.6 ± 1.4 0.2 (0.8) ± 1.4 −0.1 ± 1.1 −6.9 (4.3) ± 1.5 −2 ± 1.3 −1.6 ± 1.3 D2 ΔC Litter 7.1 ± 0.9 8.7 (2.5) ± 0.9 −6.2 ± 1.5 −0.2 (11.2) ± 1.6 −6.1 ± 1.3 −4.8 ± 1.3 D2 ΔC 5 ± 1.5 −10.3 ± 1.6 1.3 (0.3) ± 2 −10.4 ± 1.3 −8.7 ± 1.3 D2Total Total ΔC = ΔC + ΔC 17.8 ± 1.2 19.5 (5.9) ± 1.5 1.6 ± 1.2 −0.2 ± 1.2 −8.6 ± 1.3 −7.7 ± 1.3 2 LTotal D2Total The Root Mean Square Difference (RMSD) is presented as an estimate of spatial estimation error. Ferster et al. Forest Ecosystems (2015) 2:13 Page 15 of 19 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 −1 −1 −1 2003–2006 Mg C ha and the (b) mean annual value for each method in Mg C ha yr a) Site Method DF49 HDF90 HDF88 HDF00 ΔC Mean 19.4 (7.2) 2.0 (8.5) −2.5 (13.3) −10.1 (8.2) ΔC K-MSN 17.8 (5.8) ± 1.2 na/ 1.6 (11.7) ± 1.2 −8.6 (5.0) ± 1.3 ΔC K-MSN FPW 19.5 (5.9) ± 1.2 na/ −0.2 (11.2) ± 1.2 −7.7 (3.5) ± 1.3 ΣNEP Flux Tower 13.6 ± 3.6 na/ −1.9 −20.1 ΣNEP CBM-CFS3* 18.9 2.0* na/ −11.5 ΣNEP CBM-CFS3 FPW 24.1 na/ na/ −17.9 b) Site Method DF49 HDF90 HDF88 HDF00 ΔC Mean 4.9 (1.8) 0.5 (2.1) −0.63 (3.3) −2.5 (2.0) ΔC K-MSN 4.5 (1.5) ± 1.2 na/ 0.4 (2.9) ± 1.2 −2.2 (1.2) ± 1.3 ΔC K-MSN FPW 4.9 (1.5) ± 1.2 na/ −0.1 (2.8) ± 1.2 −1.9 (1.9) ± 1.3 ΣNEP Flux Tower 3.4 ± 0.9 na/ −0.5 −5.0 ΣNEP CBM-CFS3* 4.7 0.5* na/ −2.9 ΣNEP CBM-CFS3 FPW 6.0 na/ na/ −4.47 *average of all grid cells encompassing the six plots 2003 – 2006. na – not available or not applicable. At the regenerating clearcut, unweighted ΔC was have used measurements or estimates of soil and de- within 1 standard deviation of the CBM-CFS3 estimates tritus respiration in combination with inventory mea- of ΣNEP, and closest to unweighted CBM-CFS3 ΣNEP. surements of overstory productivity (Law et al. 2001; The unweighted mean of ΔC , ΔC K-MSN, and ΔC Ehman et al. 2002; Kolari et al. 2004; Black et al. 2005). 2 2 2 K-MSN FPW were greater than 1 standard deviation In the present study, comparisons were based on mea- higher than the EC flux tower estimate of ΣNEP. surements of ΔC which were higher than EC ΣNEP. ΔC K-MSN was within 1 standard deviation of the Black et al. (2005) also found a divergence between in- CBM-CFS3 unweighted ΣNEP. Both ΔC K-MSN and ventory measurements of ΔC (including ΔC measure- 2 S ΔC K-MSN FPW were greater than 1 standard ments) and EC ΣNEP, with ΔC higher than EC ΣNEP; deviation higher than ΣNEP CBM-CFS 3 FPW however, they also calculated ΣNEP using measure- (Table 6). ments of heterotrophic respiration combined with bio- At the juvenile stands, all estimates ΔC were within 1 metric estimates of NPP, and found that the bias was standard deviation of the estimates of ΣNEP. At HDF90, not observed. They concluded that ΔC systematically the unweighted estimate of ΔC matched ΣNEP CBM- overestimated NEP due to unaccounted decomposition CFS3. At HDF88, the unweighted mean of ground plots processes and uncertainties in ΔC . At DF49, Jassal was the closest to the EC tower measurements of ΣNEP, et al. (2010) measured heterotrophic soil respiration in while the ΔC K-MSN had a different sign (indicating a 2007, which was a major portion of total soil respir- small source), and the K-MSN ΔC FPW had the same ation. Including measurements of soil heterotrophic sign as ΔC unweighted (Table 6). respiration with net primary production estimates from At the near-rotation stand, all estimates of ΔC were forest inventory may improve convergence with EC- within 1 standard deviation of all estimates of ΣNEP. flux tower measurements. However, a limitation of The estimates of ΔC (and CBM-CFS3 ΣNEP) higher taking measurements of soil respiration is that the mea- than the EC flux tower measurements of ΣNEP surements require long-term collection using special (Table 6). equipment and are not typically collected in traditional forest inventory. Discussion Kolari et al. (2004) found that management practices, Several approaches can be used to compare forest inven- such as clearcut harvesting had a large impact on the tory data with EC ΣNEP. For example, several authors soil C balance, with clearcut sites soil C sources, while Ferster et al. Forest Ecosystems (2015) 2:13 Page 16 of 19 than expected given the measured litterfall inputs and estimated decomposition rates. The NFI (2008) calls for 10-year remeasurement intervals, which may be appro- priate for trees and detritus, but may be too short for ΔC . In addition, direct use of forest inventory data re- lies 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 allomet- ric equations is labour intensive, destructive, and costly, and therefore uncommon in most studies; however, use of existing relationships can lead to errors due to differ- ences 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 pro- cedure provided regionally applicable allometric equa- tions; however, the errors in this large pool of equations have not yet been estimated and therefore is a limitation. The measured changes in ΔC 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 −1 −1 MgC ha year in a Douglas-fir stand and Gielen et al. (2013) found no significant change over an eight-year Figure 8 Yearly NEP estimates for CBM-CFS3 and EC flux tower (Black interval in a Scots pine forest. In other studies by Ehman et al. 2008) for (A) the regenerating clearcut and (B) the near-rotation et al. (2002), Kolari et al. (2004), Miller et al. (2004), stand. CBM-CFS3 is shown as yearly mean of footprint cells with equal Ohtsuka et al. (2007), and Granier et al. (2008) ΔC was weighting (Fp equal) and weighting by NEP flux probability density distribution (Fp weight). assumed to be zero over the measurement interval. The large value for ΔC in this study may be explained by several factors related to measurement methodology. non-disturbed sites had more stable soil C balances. In First, mineral soil bulk densities in the 0–15 cm layer this study, a trend in soil C stocks, with the highest soil were assumed not to change and thus the BD values C stocks at the oldest site, and the lowest at the recently measured in 2002 were also used in 2006. If the bulk disturbed site, was consistent with the trend by seral density decreased then the C change would be overesti- stage reported by Kolari et al. (2004). For the two mated. Second, at the regenerating clearcut and juvenile similar-aged juvenile stands, since there was no flux stands, an increase in fine root mass from 2002 to 2006 tower installed at HDF90 and CBM-CFS3 model runs included in the <2 mm soil fraction could also have results were not available at HDF88, a direct comparison accounted for some of the increase in soil C. Third, sam- was not possible. In addition to the differences in meas- ples in 2002 were taken by 10–12 cm diameter hole ex- urement methods, the differences in ΣNEP between the cavations, while samples in 2006 were collected using a two juvenile stands could also be partially attributed to 2 cm diameter soil corer, which, due to compression at site conditions. Notably, HDF88 had more detrital C the opening of the sample corer, may be biased to soil stocks, likely leading to larger releases of C due to de- from shallower depths where C content is likely higher. composition. This was reflected in the more negative Therefore, in this study, measurements of ΔC were ΔC at HDF88 compared to HDF90. deemed unreliable, due to the large positive changes in D2 For the forest inventory measurement of ΔC and its ΔC , and assumed to be negligible over the four year components ΔC , ΔC , and ΔC , differences in proced- period and thus not used for subsequent calculations of L D S ure, measurement error by field and lab crews, and site level ΔC and further comparisons. In future work, sample design may introduce error into inventory mea- efforts to reduce any chances of variation in the way the surements that may be larger than the magnitude of samples are collected and analysed may result in more change in C stocks over short measurement intervals, reliable measurements of C stock changes. For example, such as the 4 year period used in this study. Overall, collecting and processing the samples using identical ΔC and ΔC were much larger than reported in previ- methods at nearly the same locations, or collecting a D1 S ous studies (e.g. Law et al. 2001), and was much larger much larger number of samples to capture a wider range Ferster et al. Forest Ecosystems (2015) 2:13 Page 17 of 19 of spatial variation could reduce the effect of spatial vari- “different strengths and weaknesses, but the combination ability on measurements of ΔC . of multiple measurements and modeling has the potential Several authors have applied footprint analysis in an for refining estimates of C stocks and fluxes” (Law et al. effort to correct for non-homogenous environments, 2004). For example, eddy-covariance data are temporally and horizontal advection, for example due to daily up- detailed, and provide information about daily and seasonal slope or downslope winds, which may introduce bias processes. However, the data are less spatially extensive into the eddy-covariance measurements (Baldocchi and detailed than the forest inventory data, for example, 2008). For example, Ehman et al. (2002) used footprints only three sites were available in this study each covering to calculate an estimate of sensor bias given using a an area with diverse vegetation within the tower footprint. vegetation index, and this gave confidence that the sen- Measurement error was estimated to be relatively low, on −1 −1 sor measurements were representative of the inventory the order of 0.9 Mg C ha year (Morgenstern et al. measurements. Gielen et al. (2013) used a more complex 2004). In contrast, forest inventory data are more spatially footprint model applied to the processing of NEP data, extensive and include detailed measurements of forest at- were fluxes originating outside of the area where inven- tributes. The main constraint for forest inventory data is tory measurements were sampled were removed from that they are less temporally detailed (for example, mea- the NEP estimates. In this study, footprint weighting by surements were available at only two points in time), and ecological attributes and footprint probability density the error in these measurements may be greater than the was evaluated where it facilitated comparison between expected rate of change over the measurement interval. inventory-based measurements of ecosystem ΔC stocks The inventory plots in this study sampled a broad range and tower ΣNEP by accounting for fine scale spatial of vegetation conditions and provided insight into the rep- variability in forest structure. The effect of K-MSN clas- resentativeness of flux-tower footprints. Estimates of error −1 sification and footprint weighting had an effect on the of 1.2 to 1.3 Mg C ha were identified relating to spatial estimated site-level values; however, this effect was representativeness of the inventory plots within the tower smaller than differences due to measurement method- footprints. Finally, C budget models are important for un- ology. Notably, applying footprint weighting to the derstanding landscape scale C processes. Since these estimation of CBM-CFS3 ΣNEP at the regenerating models depend on numerous coefficients, algorithms and clearcut increased convergence with EC tower ΣNEP; assumptions, it is valuable to compare model outputs with however, footprint weighting did not have the same ef- EC flux and inventory measurements to assess their per- fect for ΔC . The largest differences observed among formance and evaluate their behaviour. Considering all of methods were seen at the regenerating clearcut, where these approaches together, the relative strengths and limi- the inventory approach indicated the site was a much tation of each approach can be taken into account in lower source than the ΣNEP methods. A possible cause evaluating the suitability of each type of measurement at is that the inventory methods did not properly account varying spatial and temporal scales. for C losses from decomposition of coarse roots from stumps and large woody debris at recently harvested Conclusions sites, which can be a substantial source of respiration The comparison of ΔC from inventory measurements (Janisch et al. 2005). Therefore, if an estimate of stump with ΣNEP from CBM-CFS3 and EC flux-tower measure- coarse root decomposition were included in ΔCD it ments demonstrated an agreement among the methods in may have increased convergence with the ΣNEP estimates. trends across stand seral stage; however, due to differences Additionally, the original sample locations were designed in the measurement approaches, there was some diver- to capture the range of conditions in near proximity of the gence in the results. Convergence amongst methods was tower (based on interpretation of aerial photographs), and closest for the juvenile sites (HDF90 and HDF88), which the relationship was similar when scaled up to the site were in transition from C source- to sink. At the most re- level. Future studies may seek to derive accurate footprint cently harvested site (HDF00), the EC flux ΣNEP and estimates prior to establishing sample plots to ensure they CBM CFS3 ΣNEP indicated that the site was a greater are representative of the broader site conditions. Another source of C over the time period than ΔC from inventory more costly alternative is to establish a much greater methods. At the near-rotation site the inventory ΔC number of plots in a systematic basis around the tower, to method, unweighted CBM CFS3 ΣNEP and EC flux ΣNEP ensure the spatial variation within the tower footprint is also converged. adequately captured. Each of the measurement methods had advantages This research highlights how different measurements and limitations. For example, while EC flux towers pro- and approaches for C accounting include different con- vide temporally rich data, they had limited spatial cover- straints, and as a result, computations using different data age. Obtaining ΔC using forest inventory measurements sources may not converge. Each type of measurement has included several challenges such as the consistency of Ferster et al. Forest Ecosystems (2015) 2:13 Page 18 of 19 data collected over long time intervals, small changes in Bailey JD, Mayrsohn C, Doescher PS, St Pierre E, Tappeiner JC (1998) Understory vegetation in old and young Douglas-fir forests of western Oregon. For Ecol ΔC compared to uncertainty in the measurements, and Manage 112:289–302 insufficient measurement of processes such as decom- Baldocchi DD (2008) “Breathing” of the terrestrial biosphere: lessons learned from position of all detrital pools. Forest inventory data are a global network of carbon dioxide flux measurement systems. Aust J Bot 56:1, doi:10.1071/BT07151 available over broad areas; therefore methods that com- Barford CC, Wofsy SC, Goulden ML, Munger JW, Pyle EH, Urbanski SP, Hutyra L, bine measurements with estimates of decomposition are Saleska SR, Fitzjarrald D, Moore K (2001) Factors controlling long- and needed to provide information over a regional scale. short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science 294:1688–1691, doi:10.1126/science.1062962 Footprint analysis using the spatial predictors from re- Black K, Bolger T, Davis P, Nieuwenhuis M, Reidy B, Saiz G, Tobin B, Osborne B mote sensing and topography can give confidence that (2005) Inventory and eddy covariance-based estimates of annual carbon sample locations represent broader area site conditions. sequestration in a Sitka spruce (Picea sitchensis (Bong.) Carr.) forest ecosystem. Eur J For Res 126:167–178, doi:10.1007/s10342-005-0092-4 Considering the advantages and constraints of each type Black TA, Fredeen A, Jassal R (2008) Carbon Sequestration in British Columbia’s of measurement is helpful in choosing appropriate meas- Forests and Management Options. University of Victoria, Victoria, BC, White urement methods at varying spatial and temporal scales. paper. Pacific Institute for Climate Solutions Brown JK (1971) A planar intersect method for sampling fuel volume and surface area. For Sci 17:96–102 Competing interests Chen B, Chen JM,Mo G,Black TA,WorthyDEJ (2008) Comparison of regional carbon The authors declare that they have no competing interests. flux estimates from CO 2 concentration measurements and remote sensing based footprint integration. Global Biogeochem Cycles. doi:10.1029/2007GB003024 Chen BB, Black TA, Coops NC, Hilker T, Trofymow JA, Morgenstern K, Black AT Authors’ contributions (2009) Assessing tower flux footprint climatology and scaling between CJF set research objectives, devised and undertook analysis, and composed the remotely sensed and eddy covariance measurements. Boundary-Layer manuscript. JAT set research objectives, directed the collection and synthesis of Meteorol 130:137–167, doi:10.1007/s10546-008-9339-1 field data, summarized output from CBM-CFS3 runs, advised analysis, assisted Clark D, Brown S (2001) Measuring net primary production in forests: concepts with manuscript composition, and provided extensive editorial work on the and field methods. Ecol Appl 11:356–370 manuscript. NCC advised on research objectives, advised analysis, and assisted Coops NC, Hilker T, Wulder MA, St-Onge B, Newnham GJ, Siggins A, Trofymow JA with manuscript composition. BC advised the footprint analysis, and provided (2007) Estimating canopy structure of Douglas-fir forest stands from discrete- commentary on manuscript. TAB advised on research objectives, advised the return LiDAR. Trees-Structure Funct 21:295–310. doi:10.1007/s00468-006-0119-6 analysis, and provide extensive review and commentary throughout the work. Crookston NL, Finley AO (2008) yaImpute: an R package for kNN imputation. All authors read and approved the final manuscript. J Stat Softw 23(10):1–16 Curtis P, Hanson P, Bolstad P, Barford C, Randolph J, Schmid H, Wilson K (2002) Acknowledgements Biometric and eddy-covariance based estimates of annual carbon storage in We thank Bob Ferris, Frank Eichel, and Glenda Russo of the Canadian Forest five eastern North American deciduous forests. Agric For Meteorol 113:3–19, Service and staff of B.A. Blackwell and Associates for their help processing and doi:10.1016/S0168-1923(02)00099-0 collecting National Forest-inventory-style ground plot data. We also thank Dixon RK, Solomon AM, Brown S, Houghton RA, Trexier MC, Wisniewski J (1994) François Gougeon, CFS, for determining canopy tree density values for the Carbon pools and flux of global forest ecosystems. Science 263:185–190, sites and to Graham Stinson, CFS, for the CBM-CFS3 output files used in doi:10.1126/science.263.5144.185 Wang et al. (2011). We thank the UBC Land and Food Systems Biometeorology Ehman JL, Schmid HP, Grimmond CSB, Randolph JC, Hanson PJ, Wayson CA, and Soil Physics Group, in particular Paul Jassal, Praveena Krishnan, Kai Cropley FD (2002) An initial intercomparison of micrometeorological and Morgenstern, and Elyn Humphreys for their work processing EC flux data, ecological inventory estimates of carbon exchange in a mid-latitude deciduous and Zoran Nesic, Dominic Lessard, Andrew Sauter, and Andrew Hum for forest. Glob Chang Biol 8:575–589, doi:10.1046/j.1365-2486.2002.00492.x their work running and maintaining the EC flux tower sites. We also thank Ferster CJ, Trofymow JA, Coops NC, Chen B, Black TA, Gougeon FA (2011) Mark Johnson for his comments and suggestions. Funding for this study Determination of ecosystem carbonstock distributions in the flux footprint of an was provided by the Canadian Forest Service Pacific Forestry Centre Graduate eddy-covariance tower in a coastal forest in British Columbia. Can J For Res Student Award, a CFCAS grant to the Canadian Carbon Program (CCP), and 41:1380–1393, doi:10.1139/X11-055 Natural Sciences and Engineering Research Council of Canada (NSERC) Gielen B, De Vos B, Campioli M, Neirynck J, Papale D, Verstraeten A, Ceulemans R, Discovery Grant to NCC. The LiDAR data were acquired by Benoit St-Onge Janssens IA (2013) Biometric and eddy covariance-based assessment of of the University of Quebec at Montreal as part of an ongoing collaborative decadal carbon sequestration of a temperate Scots pine forest. Agric For project with funds provided by NSERC and BIOCAP. Meteorol 174:135–143, doi:10.1016/j.agrformet.2013.02.008 Gillis M, Omule A, Brierley T (2005) Monitoring Canada’s forests: the national Author details forest inventory. For Chron 81:214–221 Department of Forest Resources Management, University of British Göckede M, Rebmann C, Foken T (2004) A combination of quality assessment Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada. Canadian Forest tools for eddy covariance measurements with footprint modelling for the Service (Pacific Forestry Centre), Natural Resources Canada, 506 Burnside characterisation of complex sites. Agric For Meteorol 127:175–188, Road West, Victoria, BC V8Z 1M5, Canada. Faculty of Land and Food doi:10.1016/j.agrformet.2004.07.012 Systems, University of British Columbia, 2357 Main Mall, Vancouver, BC V6T 4 Gougeon F (1995) A crown-following approach to the automatic delineation of 1Z4, Canada. Department of Biology, University of Victoria, PO Box 1700, individual tree crowns in high spatial resolution aerial images. Can J Remote Station CSC, Victoria, BC V8W 2Y2, Canada. Sens 21:274–284 Gough C, Vogel C, Schmid H, Su H, Curtis P (2008) Multi-year convergence of Received: 18 December 2014 Accepted: 21 April 2015 biometric and meteorological estimates of forest carbon storage. Agric For Meteorol 148:158–170, doi:10.1016/j.agrformet.2007.08.004 Granier A, Ceschia E, Damesin C, Dufrene E, Epron D, Gross P, Lebaube S, Le Dantec V, Le Goff N, Lemoine D, Lucot E, Ottorini JM, Pontailler JY, Saugier B References (2000) The carbon balance of a young Beech forest. Funct Ecol 14:312–325, Babst F, Bouriaud O, Papale D, Gielen B, Janssens I, Nikinmaa E, Ibrom A, Wu J, doi:10.1046/j.1365-2435.2000.00434.x Bernhofer C, Köstner B, Grünwald T, Seufert G, Ciais P, Frank D (2014) Above-ground woody carbon sequestration measured from tree rings is Granier A, Bréda N, Longdoz B, Gross P, Ngao J (2008) Ten years of fluxes and coherent with net ecosystem productivity at five eddy-covariance sites. New stand growth in a young beech forest at Hesse, North-eastern France. Ann Phytol 201:1289–1303, doi: 10.1111/nph.12589 For Sci 65:704–704, doi:10.1051/forest Ferster et al. Forest Ecosystems (2015) 2:13 Page 19 of 19 Grant RF, Black TA, Humphreys ER, Morgenstern K (2007) Changes in net Ohtsuka T, Mo W, Satomura T, Inatomi M, Koizumi H (2007) Biometric based ecosystem productivity with forest age following clearcutting of a coastal carbon flux measurements and net ecosystem production (NEP) in a Douglas-fir forest: testing a mathematical model with eddy covariance temperate deciduous broad-leaved forest beneath a flux tower. Ecosystems measurements along a forest chronosequence. Tree Physiol 27:115–131 10:324–334, doi:10.1007/s10021-007-9017-z Green RN, Klinka K (1994) A Field Guide to Site Identification and Interpretation R Development Core Team (2011) R: a language and environment for statistical for the Vancouver Forest Region. Land Management Handbook Number 28. computing. R Foundation for Statistical Computing, Vienna, Austria, URL British Columbia Ministry of Forests, Research Program, Victoria, BC http://www.R-project.org/ Humphreys ER, Black TA, Morgenstern K, Cai T, Drewitt GB, Nesic Z, Trofymow JA Richardson AD, Hollinger DY (2007) A method to estimate the additional (2006) Carbon dioxide fluxes in coastal Douglas-fir stands at different stages uncertainty in gap-filled NEE resulting from long gaps in the CO2 flux record. of development after clearcut harvesting. Agric For Meteorol 140:6–22, Agric For Meteorol 147:199–208, doi:10.1016/j.agrformet.2007.06.004 doi:10.1016/j.agrformet.2006.03.018 Roberts D, Cooper S (1989) Concepts and techniques of vegetation mapping. In: Janisch JE, Harmon ME, Chen H, Fasth B, Sexton J (2005) Decomposition of Ferguson D, Morgan P, Johnson FD (eds) Land classifications based on coarse woody debris originating by clearcutting of an old-growth conifer vegetation: applications for resource management. USDA Forest Service forest. Ecoscience 12:151–160, doi:10.2980/i1195-6860-12-2-151.1 General Technical Report: INT-257; USDA Forest Service Intermountain Forest and Range Experiment Station, Ogden, UT, USA, pp 90–96 Janzen HH (2005) Soil carbon: a measure of ecosystem response in a changing Schmid H (2002) Footprint modeling for vegetation atmosphere exchange world? Can J Soil Sci 85:467–480 studies: a review and perspective. Agric For Meteorol 113:159–183, Jassal RS, Black TA, Trofymow JA, Roy R, Nesic Z (2010) Soil CO and N Oflux 2 2 doi:10.1016/S0168-1923(02)00107-7 dynamics in a nitrogen-fertilized Pacific Northwest Douglas-fir stand. Geoderma Schmid HP, Lloyd CR (1999) Spatial representativeness and the location bias of 157:118–125, doi:10.1016/j.geoderma.2010.04.002 flux footprints over inhomogeneous areas. Agric For Meteorol 93:195–209, Jungen J (1985) Soils of Southern Vancouver Island. Victoria, BC doi:10.1016/S0168-1923(98)00119-1 Kolari P, Jukka P, Rannik U, Hannu I, Pertti H, Berninger F (2004) Carbon balance Schulze ED, Hogberg P, Van Oene H, Persson T, Harrison A, Read D, Kjoller A, of different aged Scots pine forests in Southern Finland. Glob Chang Biol Matteucci G (2000) Interactions Between the Carbon and Nitrogen Cycles 10:1106–1119, doi:10.1111/j.1365-2486.2004.00797.x and the Role of Biodiversity: A Synopsis of a Study Along a North–south Kominami Y, Jomura M, Dannoura M, Goto Y, Tamai K, Miyama T, Kanazawa Y, Transect Through Europe. In: Schulze ED (ed) Carbon and Nitrogen Cycling Kaneko S, Okumura M, Misawa N, Hamada S, Sasaki T, Kimura H, Ohtani Y in European Forest Ecosystems. Springer, Heidelberg, pp 468–491 (2008) Biometric and eddy-covariance-based estimates of carbon balance for Smyth C, Trofymow JA, Kurz W (2010) Decreasing Uncertainty in CBM-CFS3 a warm-temperate mixed forest in Japan. Agric For Meteorol 148:723–737, Estimates of Forest Soil Carbon Sources and Sinks Through use of Long-Term doi:10.1016/j.agrformet.2008.01.017 Data from the Canadian Intersite Decomposition Experiment. Information Krishnan P, Black TA, Jassal RS, Chen B, Nesic Z (2009) Interannual variability of Report BC-X-422. Natural Resources Canada, Canadian Forest Service, Pacific the carbon balance of three different-aged Douglas-fir stands in the Pacific Forestry Centre, Victoria, BC Northwest. J Geophys Res 114:G04011, doi:10.1029/2008JG000912 Stage AR, Crookston NL (2007) Partitioning error components for accuracy- Kurz W, Apps M, Banfield E, Stinson G (2002) Forest carbon accounting at the assessment of near-neighbor methods of imputation. For Sci 53:62–72 operational scale. For Chron 78:672–679 Ter-Mikaelian M (1997) Biomass equations for sixty-five North American tree Kurz WA, Dymond CC, White TM, Stinson G, Shaw CH, Rampley GJ, Smyth C, species. For Ecol Manage 97:1–24, doi:10.1016/S0378-1127(97)00019-4 Simpson BN, Neilson ET, Trofymow JA, Metsaranta J, Apps MJ (2009) CBM- Theede AD (2007) Biometric and Eddy-Covariance Estimates of Ecosystem Carbon CFS3: a model of carbon-dynamics in forestry and land-use change Storage at Two Boreal Forest Stands in Saskatchewan: 1994–2004. University implementing IPCC standards. Ecol Modell 220:480–504, doi:10.1016/ of Saskatchewan, MSc thesis j.ecolmodel.2008.10.018 Trofymow JA, Stinson G, Kurz WA (2008) Derivation of a spatially explicit 86-year Lambert M-C, Ung C-H, Raulier F (2005) Canadian national tree aboveground retrospective carbon budget for a landscape undergoing conversion from biomass equations. Can J For Res 35:1996–2018, doi:10.1139/x05-112 old-growth to managed forests on Vancouver Island, BC. For Ecol Manage Landsberg JJ, Waring RH (1997) A generalised model of forest productivity using 256:1677–1691, doi:10.1016/j.foreco.2008.02.056 simplified concepts of radiation-use efficiency, carbon balance and partitioning. Wang Z, Grant RF, Arain MA, Chen BN, Coops N, Hember R, Kurz WA, Price DT, For Ecol Manage 95:209–228, doi:10.1016/S0378-1127(97)00026-1 Stinson G, Trofymow JA, Yeluripati J, Chen Z (2011) Evaluating weather effects Law B, Thornton P, Irvine J (2001) Carbon storage and fluxes in ponderosa pine on interannual variation in net ecosystem productivity of a coastal temperate forests at different developmental stages. Glob Chang Biol 7:755–777 forest landscape: a model intercomparison. Ecol Modell 222:3236–3249, Law BE, Turner D, Campbell J, Sun OJ, Van Tuyl S, Ritts WD, Cohen WB (2004) doi:10.1016/j.ecolmodel.2011.06.005 Disturbance and climate effects on carbon stocks and fluxes across Western Yashiro Y, Lee N-YM, Ohtsuka T, Shizu Y, Saitoh TM, Koizumi H (2010) Biometric- Oregon USA. Glob Chang Biol 10:1429–1444, doi:10.1111/j.1365- based estimation of net ecosystem production in a mature Japanese cedar 2486.2004.00822.x (Cryptomeria japonica) plantation beneath a flux tower. J Plant Res 123:463–472, Leclerc M, Thurtell G (1990) Footprint prediction of scalar fluxes using a doi:10.1007/s10265-010-0323-8 Markovian analysis. Boundary-Layer Meteorol 52(3):247–258 Li Z, Kurz WA, Apps MJ, Beukema SJ (2003) Belowground biomass dynamics in the Carbon Budget Model of the Canadian Forest Sector: recent improvements and implications for the estimation of NPP and NEP. Can J For Res 33:126–136, doi:10.1139/x02-165 Mahalanobis P (1936) On the generalized distance in statistics. Proc Natl Inst Sci 12:49–55 Matthews G (1993) The carbon content of trees. Forestry Commission Technical Submit your manuscript to a Paper, vol. 4. Forestry Commission, London journal and beneﬁ t from: Miller SD, Goulden ML, Menton MC, da Rocha HR, de Freitas HC, Figueira AMES, de Sousa CA D (2004) Biometric and micrometeorological measurements of 7 Convenient online submission tropical forest carbon balance. Ecol Appl 14:114–126, doi:10.1890/02-6005 7 Rigorous peer review Morgenstern K, Andrew Black T, Humphreys ER, Griffis TJ, Drewitt GB, Cai T, Nesic Z, Spittlehouse DL, Livingston NJ (2004) Sensitivity and uncertainty of the 7 Immediate publication on acceptance carbon balance of a Pacific Northwest Douglas-fir forest during an El Niño/La 7 Open access: articles freely available online Niña cycle. Agric For Meteorol 123:201–219, doi:10.1016/ 7 High visibility within the ﬁ eld j.agrformet.2003.12.003 7 Retaining the copyright to your article NFI (2008) Canada’s National Forest Inventory - Ground Sampling Guidelines. Natural Resources Canada, Canadian Forest Service, Victoria, BC, [Electronic Resource]. Available: https://nfi.nfis.org/documentation/ground_plot/Gp_ Submit your next manuscript at 7 springeropen.com guidelines_v5.0.pdf. Accessed 13 April 2015

Journal

"Forest Ecosystems" – Springer Journals

Published: Dec 1, 2015

Keywords: Ecology; Ecosystems; Forestry

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

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

Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

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

References (60)

Abstract

Journal

Recommended Articles

There are no references for this article.

Our policy towards the use of cookies