Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Using inventory based field attributes to characterize carbon stocks and carbon stock changes within… Ferster, Colin Jay 2009

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2009_fall_ferster_colin.pdf [ 3.39MB ]
Metadata
JSON: 24-1.0067765.json
JSON-LD: 24-1.0067765-ld.json
RDF/XML (Pretty): 24-1.0067765-rdf.xml
RDF/JSON: 24-1.0067765-rdf.json
Turtle: 24-1.0067765-turtle.txt
N-Triples: 24-1.0067765-rdf-ntriples.txt
Original Record: 24-1.0067765-source.json
Full Text
24-1.0067765-fulltext.txt
Citation
24-1.0067765.ris

Full Text

Using inventory based field attributes to characterize carbon stocks and carbon stock changes within eddy-flux covariance tower footprints  by  Colin Jay Ferster  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR A DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Forestry)  The University of British Columbia (Vancouver) October 2009 © Colin Jay Ferster, 2009  Abstract Forests are an important part of the global carbon (C) cycle, and understanding and quantifying forest C dynamics is necessary for informed forest management decisions and accurate C budget accounting. This thesis contributes to the understanding of forest C dynamics by comparing biometric measurements of C stock changes (∆C) and eddy covariance (EC) flux-tower measurements of cumulative net ecosystem productivity (ΣNEP) at three sites at the Fluxnet Canada Research Network British Columbia Flux Station; a young (near-end-of-rotation) forest established in 1949 (DF1949), a pole sapling stand established in 1988 (HDF1988), and a recent clearcut established in 2000 (HDF2000). To address spatial variability in stand and tower footprint conditions, first, light detection and ranging (lidar) remote sensing data were used to quantify large tree and snag aboveground mass (TSAM) at DF1949 (r2=0.75, SEE=29.68 Mg/ha). Next, for all three sites, remote sensing estimates were combined with advanced GIS data to estimate the spatial distribution of C stocks that are more difficult to measure directly using only remote sensing data. The resulting spatial representation of forest structure was combined with published EC flux-tower footprint probability distributions to enable a comparison between biometric ∆C stocks and tower ΣNEP. Where biometric ∆C stocks were difficult to resolve, changes were modeled using parameters from the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3). The best agreement between biometric ∆C stocks and tower ΣNEP was at DF1949 (mean = 15.18 ± 7.94 MgC/ha/4 years ∆C; 13.63 MgC/ha/4 years ΣNEP). The other two sites followed the same rank but had larger  ii  divergences and differences of sign (HDF1988: mean = 11.52 ± 1.17 MgC/ha/4 years ∆C and -1.93 MgC/ha/4 years ΣNEP; HDF2000: 3.59 ± 1.59 MgC/ha/4 years ∆C and -20.08 MgC/ha/4 years ΣNEP). An unaccounted source of respiration at the more recently disturbed sites may be related to stump, coarse root, and loggingslash decomposition. This study develops an approach and methodology that may be applied at other EC flux-tower sites where researchers wish to compare biometric-based measurements with micrometeorological-based measurements.  iii  Table of Contents Abstract.................................................................................................................... ii Table of Contents ................................................................................................... iv List of Tables ......................................................................................................... vii List of Figures....................................................................................................... viii Acknowledgements ................................................................................................ ix Co-authorship Statement ....................................................................................... x 1. Introduction......................................................................................................... 1 1.1 The Carbon Cycle ................................................................................... 1 1.2 EC Flux-Towers...................................................................................... 2 1.3 Biometric Data ........................................................................................ 3 1.4 Spatial Scaling of Measurements............................................................ 4 1.5 Comparing Biometric and EC Flux-Tower Measurements of Fluxes .... 6 1.6 References............................................................................................. 10 2. Aboveground Large Tree Mass Estimation in a Coastal Forest in British Columbia using Plot Level Metrics and Individual Tree Detection from Lidar....................................................................................................................... 14 2.1 Introduction........................................................................................... 14 2.2 Methodology ......................................................................................... 16 2.2.1 Study Site and Groundplots ....................................................... 16 2.2.1 Individual Tree Detection from Lidar Data ............................... 18 2.2.3 Lidar Aboveground Mass Estimation ........................................ 19 2.3 Results................................................................................................... 20 2.4 Discussion ............................................................................................. 24 2.5 References............................................................................................. 27 3. Determination of Carbon Stock Distributions in the Flux Footprint of an Eddy-Covariance Tower in a Coastal Forest in British Columbia .................. 29 3.1 Introduction........................................................................................... 29 3.2 Materials and Methods.......................................................................... 34 3.2.1 Study Site ................................................................................... 34 3.2.2 GIS Data..................................................................................... 35 3.2.3 Ground plot Inventory Data ....................................................... 39 3.2.4 GIS Stratification Units Derivation ........................................... 39 3.2.5 Lidar Data .................................................................................. 41 3.2.6 Lidar Estimates of TSAM Spatial Distribution ......................... 41 3.2.7 EC Flux-Tower Footprint .......................................................... 42 3.2.8 Comparing Lidar and GIS Stratification Estimates of TSAM... 44 3.2.9 GIS Stratification Estimates of Other Component C Stocks ..... 45 3.3 Results................................................................................................... 45 3.3.1 GIS Stratification Units.............................................................. 45 3.3.2 GIS Stratification and Lidar Estimates of TSAM...................... 48 3.3.3 Stratification Estimates of C Stock Components....................... 51  iv  3.4 Discussion ............................................................................................. 53 3.4.1 GIS Stratification Units.............................................................. 54 3.4.2 GIS Stratification Estimate of TSAM........................................ 55 3.4.3 Lidar estimate of TSAM ............................................................ 55 3.4.4 Comparing Estimates of TSAM................................................. 56 3.4.5 Stratification Estimates of C Stock Components....................... 57 3.5 References............................................................................................. 60 4. Comparing Biometric Carbon Stock Changes with Eddy-Covariance FluxTower Based Cumulative NEP Measurements at Three Sites in a Coastal Forest in British Columbia................................................................................... 65 4.1 Introduction........................................................................................... 65 4.1.1 EC Flux-Towers..................................................................... 65 4.1.2 Biometric Measurements of Forest C Stocks......................... 66 4.1.3 Comparing EC Flux-Tower and Biometric Measurements ... 67 4.2 Methods................................................................................................. 70 4.2.1 Study Area ............................................................................. 70 4.2.2 DF1949 .................................................................................. 71 4.2.3 HDF1988................................................................................ 72 4.2.4 HDF2000................................................................................ 72 4.2.5 Site Stratification ................................................................... 73 4.2.6 Biometric C Stocks and Changes........................................... 74 4.2.7 Overstory Trees...................................................................... 75 4.2.8 Understory Vegetation ........................................................... 76 4.2.9 Woody Debris ........................................................................ 76 4.2.10 Forest Floor.......................................................................... 76 4.2.11 Mineral Soil ......................................................................... 77 4.2.12 Litterfall ............................................................................... 77 4.2.13 Modeled Changes in Fine Woody Debris, Forest Floor, and Mineral Soil C................................................................................. 77 4.2.14 Comparing Biometric Estimates of ∆C Stocks and EC Tower ΣNEP............................................................................................... 78 4.3 Results................................................................................................... 79 4.3.1 Site Stratification ................................................................... 79 4.3.2 Biometric C Stock Changes ................................................... 85 4.3.3 Overstory................................................................................ 87 4.3.4 Understory.............................................................................. 88 4.3.5 Coarse Woody debris............................................................. 88 4.3.6 Forest Floor............................................................................ 88 4.3.7 Mineral Soil ........................................................................... 89 4.3.8 Modeled Changes in Fine Woody Debris, Forest Floor, and Mineral Soil C................................................................................. 89 4.3.9 Comparing EC Tower ΣNEP and Biometric ∆C Stocks ....... 90 4.4 Discussion ............................................................................................. 91 4.5 Conclusion ............................................................................................ 94 4.6 References............................................................................................. 97  v  5. Conclusion ....................................................................................................... 101 5.1 Lidar Estimation of Biomass .............................................................. 101 5.2 Spatial Extrapolation of Ground Plot Data ......................................... 102 5.3 Comparing Biometric-Based and EC Flux-Tower-Based Measurements ........................................................................................... 103 5.4 References........................................................................................... 106  vi  List of Tables Table 2.1 Table 2.2 Table 3.1 Table 3.2 Table 3.3 Table 3.4 Table 4.1 Table 4.2 Table 4.3 Table 4.4  Matches of TreeVaW and ground-plot measured trees ...................... 21 Models for TSAM............................................................................... 23 DF1949 Environmental predictor variable suitability ........................ 37 DF1949 Environmental predictor variable selection .......................... 47 DF1949 TSAM estimates from lidar and GIS .................................... 51 DF1949 Unweighted and weighted mean C stocks ............................ 52 Three site environmental perdictor variable selection ........................ 81 Measured and modeled ∆C stocks ...................................................... 86 Modeled LFH, FWD, and mineral soil ∆C stocks.............................. 87 Comparing biometric ∆C stocks and tower ΣNEP ............................. 91  vii  List of Figures Figure 2.1 Figure 2.2 Figure 2.3 Figure 3.1 Figure 3.2 Figure 3.3 Figure 3.4 Figure 3.5 Figure 4.1 Figure 4.2 Figure 4.3 Figure 4.4  Lidar predictions of individual tree heights....................................... 21 Lidar predictor variables for above ground mass .............................. 22 Correlations for lidar above ground mass ......................................... 24 Coastal BC Station ............................................................................ 35 DF1949 EC flux-tower footprint estimates ....................................... 43 DF1949 Most similar neighbour classification ................................. 48 DF1949 lidar and GIS estimates of aboveground mass .................... 50 DF1949 site C stock extrapolation .................................................... 53 Coastal BC Station ............................................................................ 71 Three site EC flux footprint estimations ........................................... 82 Three site most similar neighbour classification ............................... 83 Three site Mahalanobis distance ....................................................... 84  viii  Acknowledgements Thanks to Bob Ferris, Eric Bol, Gurp Thandi, Frank Eichel, and Glenda Russo of the Canadian Forest Service (CFS) and staff of B.A. Blackwell and Associates for their help collecting and processing National Forest Inventory-style ground plot data. Sample chemical analyses were conducted by the Chemical Services laboratories at the CFS Pacific Forestry Centre and the BC Ministry of Forests, North Roads Lab. Thanks to the Biometeorology and Soil Physics Group in the Faculty of Land and Food Systems, UBC, who are responsible for measurement, analysis, and gap-filling the EC flux NEP data. Funding for this study was provided by the Canadian Forest Service Pacific Forestry Centre graduate student award, CFCAS grant to the Canadian Carbon Program (CCP), a National Sciences and Engineering Research Council of Canada (NSERC) discovery grant to Coops, University of British Columbia (UBC) and by the Natural Resource Canada PERD program. I would like to thank Dr. Nicholas Coops, Dr. Tony Trofymow, and Dr. Andy Black for their time, support, and inspiration. I would also like to thank the members of the IRSS lab and staff at the PFC for their company and support.  ix  Co-authorship Statement This thesis is a combination of three scientific papers of which I am lead author. For these scientific journal submissions, I performed primary research on data which was previously collected, data analysis, and prepared the final manuscripts. Co-authors provided advice and guidance on methodology and made valuable editorial comments as required. This research is an objective of the original Fluxnet Canada Research Network British Columbia British Columbia Flux Station project proposal. Dr Tony Trofymow first conceived the approach and offered methodological guidance throughout the process. Dr. Nicholas Coops, provided considerable direction toward the remote sensing and modelling aspects of this study. Dr. Andy Black provided guidance for the EC flux and other C cycle science aspects of the study.  x  1. Introduction 1.1 The Carbon Cycle The global carbon (C) biogeochemical cycle has important roles for ecosystem form and function. Carbon dioxide (CO2), is a greenhouse gas, and as a result of activities including fossil fuel consumption and land use change, atmospheric concentrations of CO2 have increased rapidly, raising concern about climate change (IPCC, 2007). In addition, plants depend on CO2 for photosynthesis and C composes about half of a plant's living mass. As a result, forests contain large C stocks, estimated to be approximately 1146 Pg C globally (Dixon et al., 1994), and 85.9 Pg C in Canada (Kurz et al. 1999). Forests release C into the atmosphere through decomposition (decay of forest materials) and respiration (use of carbohydrate to provide energy). Magnitudes of forest C fluxes are determined by a range of factors such as harvest (Humphreys et al., 2006), natural disturbance (Kurz et al., 2008), fertilization (Jassal et al., 2009) and weather (Krishnan et al., 2009; Morgenstern et al., 2006).  Informed forest management decisions may mitigate the release of CO2 to the atmosphere (Winjum et al., 1993), which requires a thorough understanding of forest C dynamics (Kurz et al., 2009). Results from in depth scientific studies of the C cycle are commonly integrated into advanced forest process models, which are in turn used to inform forest policy and for national and international reporting (Kurz et al., 2002). In Canada, the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS) integrates knowledge gained from many detailed studies of  1  forest C processes and has been used for reporting for international policy agreements, including the United Nations Framework Convention on Climate Change (UNFCC) and the Montreal Process (Kurz et al., 2009). Exercises, such as historical C budget modelling are used to test and fine tune the model in a variety of specific conditions (Trofymow et al., 2008), leading to incremental improvements in each version of the model.  1.2 EC Flux-Towers The global network of eddy-covariance (EC) flux-towers has done much to increase the understanding of forest C dynamics (Baldocchi, 2008). EC flux-towers use micrometeorological equipment to make stand scale continuous half-hourly measurements of CO2, sensible heat, and water vapour fluxes, and at some sites to measure other trace gasses. More than 400 sites have been established globally, covering a range of ecosystems, many operating continuously over time periods from years to decades. In Canada, the Canadian Carbon Program (CCP), has established 12 research stations and affiliated research stations across the country in a broad range of ecosystem types. The data record provided by these towers covers a range of management and weather conditions, thus increasing the understanding of the effects of climate variability and disturbance on Canada's forest C dynamics and reducing uncertainty in C budget model predictions (http://fluxnet-canada.ca).  2  However, several challenges exist for the measurement and analysis of CO2 fluxes using EC flux-tower methodology. First, at sites with complex topography and highly variable spatial distributions of vegetation, special considerations must be made to properly interpret EC flux data (Göckede, et al., 2004). Recent advances in flux footprint modelling, which evaluate the source area contributing to the measured fluxes, provide a framework to accurately evaluate EC flux-tower measurements over complex surfaces (Chen et al., 2009). Second, each site requires a large investment in equipment and expertise to operate, and as a result, the global network of EC flux sites still has a relatively limited spatial extent compared to other datasets provided by remote sensing, or on-the-ground biometric measurements.  1.3 Biometric Data In contrast to the network of EC sites, there exists a large network of biometric sample plots as part of national and regional forest inventories. For example, the Canadian National Forest Inventory (NFI) aims to establish sample plots at approximately 22,000 locations across Canada (Gillis et al., 2005). Forest C fluxes can be inferred from the difference between biometric measurements of forest C stocks at two points in time with a long interval (e.g. 4 – 10 years) in between (∆C) (Clark et al., 2001). An advantage of biometric methods of determining C fluxes is that information about the individual forest ecological elements contributing to fluxes is available since measurements are made of individual components at the plot scale. For example, trees and woody debris pieces are measured within plot  3  boundaries or along transects and allometric relationships are commonly applied to estimate C content. In addition, understory vegetation, forest floor material, and soil are destructively sampled and sent for laboratory analysis. Measurements of individual components are then summed to find plot-level totals.  Challenges to proper interpretation of biometric measurements also exist. Certain elements are difficult to measure, for example coarse root, and harvest slash which may decompose and contribute strongly to site-level respiration (Jansich et al. 2005). Since long intervals pass between measurements, different field crews may have different measurement biases, or slightly different methods may be used. Further, measurement inaccuracy may be high relative to small magnitude changes for some components. For example, soil C exhibits a large amount of spatial variability and has been a source of considerable uncertainty in previous studies (Curtis et al., 2002).  1.4 Spatial Scaling of Measurements Measurements from biometric sample plots must be spatially scaled to the extent of the landscape they are sampling. A common approach in forestry applications is stratified random sampling using expert interpretation of aerial photographs to define strata. In this approach, forest attributes are determined through the expert interpretation of aerial photographs, areas with common attributes are delimited as strata, and sample plots are measured within the strata (Haris and Dawson, 1979). This sampling scheme aims to increase sampling efficiency compared to simple  4  random sampling (Schreuder et al., 2004). A limitation of this approach is that it depends on the skills and judgement of expert airphoto interpreters and is highly labour intensive.  Modern remote sensing and geographic information systems (GIS) datasets provide quantitative spatial data sources more advanced than traditional aerial photographs. For example, light detection and ranging (lidar) has been used to reliably measure a number of forest attributes and has been used as part of Scandinavian forest inventories (Naesset, 2004), and a variety of lidar sensors have been used to make robust, regionally applicable, estimates of forest biomass (Lefsky et al., 2005). Multispectral imagery has been linked to vegetation composition and structure (Carlson and Ripley, 1997), and high spatial resolution optical data has been used to extract information about individual tree crown dimensions (Gougeon et al., 2005). As well, innovative GIS datasets have been developed to integrate historic inventory and disturbance records from a wide variety of sources (Trofymow et al., 2008).  Further, advanced classification algorithms have been developed to quantitatively classify landscapes using modern remote sensing and GIS datasets. One approach, the most-similar-neighbour classification, can use statistical distances to measure the similarity between landscape units (for example, a footprint probability density function cell) and measured inventory plots using environmental target variables (Crookston and Finley, 2004). Then the target biometric variables from the most similar biometric sample plot (determined by the shortest statistical distance for  5  predictor variables) are assigned to each landscape unit. The statistical distance measure, such as the Mahalanobis distance (Mahalanobis, 1936), is also useful to evaluate the representativeness of sample plots for a landscape. This approach is non-parametric and non-distributional, which is beneficial in landscapes which are not normally distributed.  1.5 Comparing Biometric and EC Flux-Tower Measurements of Fluxes Biometric and EC flux-tower measurements of fluxes are complementary, and comparisons are therefore valuable to our understanding of forest C dynamics. Biometric measurements can provide estimates of the component level changes driving stand level fluxes. In turn, EC flux-tower instrumentation can measure complete stand level fluxes and provide information about important fluxes that may be missed using biometric methods. Further, measurements made using EC flux-tower based instrumentation can provide finer temporal resolution data, as measurements made using biometric methods are limited to a long interval. Overall, estimating stand level fluxes using both independent methods may add confidence to C budget estimates.  Globally, comparative studies between biometric ∆C stocks and EC flux-tower ΣNEP have been completed at a limited number of sites (Gough et al., 2008; Kominami et al., 2007; Miller et al., 2004; Curtis et al., 2002; Ehman et al., 2002; Barford et al., 2001; Grainier et al., 2000; Schulze et al., 2000). As well, despite the many recent advances in flux-footprint models, footprint weighting has been  6  applied to compare EC flux-tower measurements at very few sites (Ehman 2002). At present, comparing biometric measurements with EC flux-tower measurements is a priority for the CCP. So far, this comparison has only been carried out at two sites at one Canadian flux station, the Boreal Ecosystems Monitoring Sites (BERMS) station in Northern Saskatchewan over 10 years with promising results (15.6 ± 4.0 (biometric) and 18.2 ± 8.09 MgC/ha (tower ∑NEP) at the Old Aspen site and 5.8 ± 2.0 MgC/ha (biometric) and 6.9 ± 1.6 MgC/ha (Tower NEP) at the Old Jack Pine site) (Theede, 2008). However, given the relatively ideal site conditions, spatial footprint probability density and stand structure distributions were not included in this study.  In contrast, the British Columbia Flux Station possesses complex stand structure and topography, requiring high resolution mapping of forest structure and analysis of the footprint probability density distribution (Chen et al., 2009). Despite complex site conditions, understanding C fluxes in this stand is critical, since it is representative of Canada's most productive forests. This thesis will compare biometric and EC flux-tower measurements of C flux at the CCP British Columbia Flux Station, answering two important research questions: 1.  How can ground plots be spatially extrapolated within the tower  footprints so that they are representative of the range of forest land characteristics that the tower is measuring? 2.  What is the change in biometric ground plot carbon stocks (∆C)  between 2002 and 2006 and how does the estimate of ∆C compare with  7  tower-based measurements of cumulative net ecosystem productivity (ΣNEP) for the same period?  These research questions were answered in three chapters:  In the first chapter, lidar remote sensing was utilized to estimate large tree and snag above ground mass (TSAM) distributions using plot level metrics and detection of individual trees. The TreeVaW algorithm (Popescu, 2007) was used to detect individual trees in lidar data, and for each plot, individual tree metrics were calculated (e.g. number of stems, mean height, and maximum height). In addition, for each plot, traditional lidar metrics of return height percentiles and return height densities were calculated. Two multiple linear regression models were developed to estimate plot level TSAM: one model using lidar metrics alone, and one model using plot level lidar metrics in combination with individual tree metrics. After this analysis, the two models were compared.  In the second chapter, a suite of GIS and non-lidar remote sensing coverages were used together with advanced classification algorithms to establish site footprint level estimates of C stock distributions. A correlation matrix was used to select significant environmental predictor variables for all measured 2002 C stock components. The selected variables were used in a nearest neighbour classification using the Mahalanobis distance, which is based on a transform of the covariance matrix of predictor variables to normalize the values. Spatial distributions of TSAM from the nearest neighbour site classification were compared with spatial  8  distributions of TSAM from the lidar measurement. Once confidence was gained from this independent comparison, the nearest neighbour classification was used to estimate spatial distributions of components that can not be measured directly using remote sensing, such as understory, woody debris, and soil C. The Mahalanobis distance, from the nearest neighbour classification, was used to assess the spatial representativeness of the measured ground sample plots.  In the third chapter, the difference between 2002 C stocks and 2006 C stocks was determined (∆C). For some difficult to measure components, ∆C was determined using a modeling approach. The site classification previously developed was used to determine the spatial distribution of ∆C, and footprint probability density distributions (from Chen et al., 2009) were used to calculate weighted means, allowing comparison with EC flux-tower ΣNEP, given the complex topography and vegetation distributions at the three sites.  9  1.6 References Baldocchi, D. 2008. Breathing of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems. Australian Journal of Botany. 56:1-26. Barford, C. C., Wofsy, S.C., Goulden, M.L., Munger, J.W., Pyle, E.H., Urbanski, S.P., Hutyra, L., Saleska, S.R., Fitzjarrald, D., and Moore, K. 2001. Factors controlling long- and short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science. 294:1688-1691. Carlson, T., and D. Ripley. 1997. On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sensing of Environment. 62:241-252. Chen, B., Black, T.A., Coops, N.C., Hilker, T., Trofymow, J.A. (Tony), and Morgenstern, K. 2009. Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Boundary-Layer Meteorology. 130:137-167. Clark, D., Brown, S., Kicklighter, D., Chambers, J., Thomlinson, J., and Ni, J. 2001. Measuring net primary production in forests: Concepts and field methods. Ecological Applications. 11:356-370. Crookston, N. L., and Finley, A.O. 2008. yaImpute: an R package for k-NN imputation. Journal of Statistical Software. 23. Curtis, P., Hanson, P., Bolstad, P., Barford, C., Randolph, J., Schmid, H., and Wilson, K. 2002. Biometric and eddy-covariance based estimates of annual carbon storage in five eastern North American deciduous forests. Agricultural and Forest Meteorology. 113:3-19. Dixon, R. K., Solomon, A.M., Brown, S., Houghton, R.A., Trexier, M.C., and Wisniewski, J. 1994. Carbon pools and flux of global forest ecosystems. Science. 263:185-190. Ehman, J. L., Schmid, H.P., Grimmond, C.S.B., Randolph, J.C., Hanson. P.J., Wayson, C., and Cropley, F.D. 2002. An initial intercomparison of micrometeorological and ecological inventory estimates of carbon exchange in a mid-latitude deciduous forest. Global Change Biology. 8:575-589. Gillis, M., Omule, A., and Brierley, T. 2005. Monitoring Canada's forests: The National Forest Inventory. Forestry Chronicle. 81:214-221.  10  Gougeon, F.A. 1995. A crown-following approach to the automatic delineation of individual tree crowns in high spatial resolution aerial images. Canadian Journal of Remote Sensing. 21:274-284. Gough, C., Vogel, C., Schmid, H., Su, H., and Curtis, P. 2008. Multi-year convergence of biometric and meteorological estimates of forest carbon storage. Agricultural and Forest Meteorology. 148:158-170. Göckede, M., Rebmann, C., and Foken, T. 2004. A combination of quality assessment tools for eddy covariance measurements with footprint modelling for the characterisation of complex sites. Agricultural and Forest Meteorology. 127:175-188. Granier, A., Ceschia, E., Damesin, C., Dufrene, E., Epron, D., Gross, P., Lebaube, S., Dantec, V.L., Goff, N.L., Lemoine, D., Lucot, E., Ottorini, J.M., Pontailler, J.Y., and Saugier, B. 2000. The carbon balance of a young beech forest. Functional Ecology. 14:312-325. Grant, R. F., Black, T. A., Humphreys, E. R., and Morgenstern, K. 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 Physiology 27:115-131. Harris, J.W.E., and Dawson, A.F. 1979. Evaluation of aerial forest pest damage survey techniques in British Columbia. Canandian Forest Service Information Report BC-X-198. Humphreys, E. R., Black, T.A., Morgenstern, K., Cai, T., Drewitt, G.B., Nesic, Z., and Trofymow, J.A. (Tony). 2006. Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of development after clearcut harvesting. Agricultural and Forest Meteorology. 140:6-22. Intergovernmental Panel on Climate Change (IPCC). 2007. Fourth assessment report. Cambridge University Press. New York, NY. Jassal, R. S., Black, T.A., Trofymow, J.A. (Tony), Roy, R., Nesic, Z. 2009. Forestfloor CO2 and N2O flux dynamics in a nitrogen-fertilized Pacific Northwest Douglas-fir stand. Soil Biology and Biochemistry, submitted. Janisch, J. E., M. E. Harmon, H. Chen, B. Fasth, and J. Sexton. 2005. Decomposition of coarse woody debris originating by clearcutting of an old-growth conifer forest. Ecoscience. 12:151–160. Kominami, Y., Jomura, M., Dannoura, M., Goto, Y., Tamai, K., Miyama, T., Kanazawa, Y., Kaneko, S., Okumura, M., Misawa, N., Hamada, S., Sasaki, T., Kimura, H., and Ohtani, Y. 2008. Biometric and eddy-  11  covariance-based estimates of carbon balance for a warm-temperate mixed forest in Japan. Agricultural and Forest Meteorology. 148:723-737. Kurz, W. A., and Apps, M.J. 1999. A 70-year retrospective analysis of carbon fluxes in the Canadian forest sector. Ecological Applications 9:526-547. Kurz, W., Apps, M., Banfield, E., and Stinson, G. 2002. Forest carbon accounting at the operational scale. Forestry Chronicle. 78:672-679. Kurz, W., and Apps, M. 2006. Developing Canada's national forest carbon monitoring, accounting and reporting system to meet the reporting requirements of the Kyoto Protocol. Mitigation and Adaptation Strategies for Global Change. 11:33-43. Kurz, W. A., Dymond, C.C., Stinson, G., Rampley, G.J., Neilson, E.T., Carroll, A.L., Ebata, T., and Safranyik, L. 2008. Mountain pine beetle and forest carbon feedback to climate change. Nature 452:987-990. Kurz, W., Dymond, C., White, T., Stinson, G., Shaw, C., Rampley, G., Smyth, C. Simpson, B., Neilson, E., Trofymow, J.A. (Tony), Metsaranta, J., and Apps, M. 2009. CBM-CFS3: A model of carbon-dynamics in forestry and land-use change implementing IPCC standards. Ecological Modelling. 220:480-504. Lefsky, M. A., Hudak, A.T., Cohen, W.B., and Acker, S. 2005. Geographic variability in lidar predictions of forest stand structure in the Pacific Northwest. Remote Sensing of Environment 95:532-548. Mahalanobis, P.C. 1936. On the generalized distance in statistics. Proceedings of the National Institute of Science of India. 12. pp. 49-55. Miller, S., Goulden, M., Menton, M., da Rocha, H., de Freitas, H., Figueira, A., and de Sousa, C. 2004. Biometric and micrometeorological measurements of tropical forest carbon balance. Ecological Applications. 14:114-126. Morgenstern, K., T. Black, T.A., Humphreys, E.R., Griffis, T.J., Drewitt, G.B., Cai, T., Nesic, Z., Spittlehouse, D.L., and Livingston, N.J. 2004. Sensitivity and uncertainty of the carbon balance of a Pacific Northwest Douglas-fir forest during an El Nino/La Nina cycle. Agricultural and Forest Meteorology. 123:201–219. Næsset, E. 2004. Practical large-scale forest stand inventory using a small-footprint airborne scanning laser. Scandinavian Journal of Forest Research 19:164. Popescu, S. 2007. Estimating biomass of individual pine trees using airborne lidar. Biomass and Bioenergy. Vol. 31. pp. 646-655.  12  Schreuder, H. T., Ernst, R., and Ramirez-Maldonado, H. 2004. Statistical techniques for sampling and monitoring natural resources. U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. RMRS-GTR-126:111. Schulze, E. D., Hogberg, P., Van Oene, H., Persson, T., Harrison, A.F., Read, D., Kjoller, A., and Matteucci, G. 2000. Interactions between the carbon and nitrogen cycle and the role of biodiversity: A synopsis of a study along a north-south transect through Europe. Ecological Studies:468–492. Trofymow, J. A., Stinson, G., and Kurz, W.A. 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. Forest Ecology and Management. 256, 1677-1691. Winjum, J., Dixon, R., and Schroeder, P. 1993. Forest management and carbon storage - an analysis of 12 key forest nations. Water Air and Soil Pollution 70:239-257.  13  2. Aboveground Large Tree Mass Estimation in a Coastal Forest in British Columbia using Plot Level Metrics and Individual Tree Detection from Lidar1 2.1 Introduction Aboveground forest carbon (C) stocks are closely linked to canopy height and stem size. Carbon composes approximately half of the live vegetation on a dry mass basis in forests (biomass), including trees (trunks, branches, foliage, and roots), shrubs, herbs, mosses, as well as decaying forest detritus (snags (ie. standing dead trees), woody debris, forest floor litter, and humus) (Chapin et al., 2006; Landsberg and Waring, 1997). National reporting requirements for the United Nations Framework Convention on Climate (UNFCC), as well as national and regional level policy decisions (Gillis et al., 2005), make the timely, accurate, and cost effective measurement of forest biomass a common goal. Light detection and ranging (lidar) is a relatively new remote sensing technology that is well suited to describing the vertical structure of forest canopies (Coops et al., 2007; Wulder et al., 2008) and has significant potential to measure aboveground forest biomass.  Lidar data is commonly acquired from an airborne platform equipped with a lidar instrument which emits pulses of infrared light in a scanning pattern, measuring the time for an emitted pulse to reflect from ground surfaces and return to the sensor. Using accurate positional information obtained from a GPS, and detailed data on aircraft orientation from an inertial measurement unit (IMU), the height and 1  A version of this chapter has been published. Ferster, C. J., Coops, N.C., and Trofymow, J.A. (Tony). 2009. Aboveground large tree mass estimation in a coastal forest in British Columbia using plot-level metrics and individual tree detection from lidar. Canadian Journal of Remote Sensing. 35:270 - 275.  14  position of features on the ground can be determined very accurately. Small footprint lidar devices often use beam footprints approximately 20 cm in diameter with very close spacing (~1.5 m between returns). Discrete return lidar devices may record up to five returns from an emitted pulse, with first returns typically incident from branches, foliage, and stems at the top of the canopy, subsequent returns originating from objects within the canopy, and the final return from the underlying terrain (Wulder et al., 2008). As part of national operational forest inventories, lidar data has been applied to measure plot-level, aboveground mass of trees including trunks, branches, and foliage of living and dead standing trees (Naesset, 2004). To predict plot-level aboveground biomass, inferential relationships are often applied using descriptive statistics from the point cloud of returns as predictor variables, including percentiles of return heights, and density of returns at height intervals within the canopy (Næsset, 2004; Lefsky et al., 2005; Næsset and Gobakken, 2008).  In contrast to approaches which use lidar return height and density metrics to estimate plot-level biomass, high density small footprint lidar data has also been used to identify individual tree crowns, measure individual tree height, and subsequently estimate individual tree biomass (Hyyppa et al., 2001; Popescu and Wynne, 2004; Popescu, 2007; Persson et al., 2002). For trees that are individually identified, height measurements from lidar data have accuracy approaching that of ground based measurements commonly made using vertex hypsometers, but covering much broader spatial extents (Andersen et al., 2006; Persson et al., 2002). As a result, biomass and volume have been accurately estimated using allometric  15  relationships using height and crown size derived from lidar data (Hyyppa et al., 2001; Persson et al., 2002; Popescu, 2007).  Thus, lidar data has successfully been used to measure forest biomass using two broad methods, first the estimation of plot-level biomass using return heights and densities summarised within a forest stand or plot, and secondly at the individual tree level where trees are identified and lidar height and crown dimensions used to estimate individual tree biomass. Despite these successes, however, few studies have combined these approaches, by first identifying individual trees and combining these tree level lidar results with plot-level lidar metrics of height and density to calculate plot-level aboveground biomass. In this short communication, we develop and apply this approach to estimate large tree and snag aboveground mass (total live and dead standing) (TSAM) in a forest in coastal British Columbia, and compare the predictions with measurements derived from Canadian National Forest Inventory (NFI)-style ground plots (NFI, 2004).  2.2 Methodology 2.2.1 Study Site and Groundplots The study site is located in the Oyster River watershed on the east coast of Vancouver Island, near Campbell River, British Columbia, Canada within the Coastal Western Hemlock very dry maritime biogeoclimatic subzone (CWHxm) (Pojar et al., 1991). The site is one of three flux-tower sites being studied as part of the Canadian Carbon Program (CCP) British Columbia Station (CCP, 2008). This second-growth stand was established by planting in 1949 following clearcut  16  harvesting and broadcast burning of the original old-growth stand. The overstory is dominated by Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco) with smaller amounts of western hemlock (Tsuga heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata Donn ex D. Donn), and red alder (Alnus rubra Bong.). Though the site is mapped as a single forest cover type and site index by forest companies operating in the area, the sloping topography (260 m to 470 m elevation) and disturbance history have resulted in finer scale variation in forest cover across the 33 ha study site.  In autumn of 2002, field measurements were made at twelve ground plots following NFI guidelines (NFI, 2004) to determine the range of stand characteristics in the study site. Plots were clustered in groups of three, systematically located at representative sites according to the range of mapped forest cover and ecosite series within the forest area estimated in 2002 to be contributing to flux-tower measurements (Humphreys et al., 2006). Live and dead trees (snags) greater than 9.0 cm diameter at breast height (DBH – measured 130 cm above ground) were measured for species, height, diameter, height to live crown, canopy class (dominant, codominant, intermediate and suppressed), and stem location in 11.28 m radius circular large tree plots (0.04 ha). Tree height measurements were made using a vertex hypsometer and stem mapping was completed using a compass and vertex hypsometer. In some cases, in dense homogeneous parts of the stand, half plots (0.02 ha) were measured when representative of stand conditions. All measured trees within the full extent of each plot were used for assessment of TreeVaW measurements of stem location and  17  height. For consistency in this study, all full plots were then divided into half plots and only half of the plot, selected at random, was used to estimate plot-level tree and snag mass. The same plot boundaries that were used for calculation of plotlevel tree and snag mass from biometric data were also used for calculation of plotlevel TreeVaW metrics and lidar return metrics. Small trees (0.5 -9.0 cm DBH) were measured on the large plot or in a nested 3.99 m radius small tree plot, but were not used in this analysis which focused on large trees only. Field data was submitted to the NFI database compilation program for calculation of large tree and snag aboveground mass (stem, bark, branches, and foliage of live trees; stem, bark, and branches of standing dead trees) using published national species-specific allometric equations (NFI, 2004). Estimates for each tree were then summed within each plot to determine plot-level large tree and snag aboveground mass (TSAM).  2.2.1 Individual Tree Detection from Lidar Data Lidar data was collected on June 8, 2004 with a beam footprint of 0.19m and ground point sampling density vendor-reported scene-wide average of 0.7 returns m −2 (Coops et al., 2007). Within the sample plots, actual return density ranged  between 1.8 and 7.3 returns m −2 (mean 5.3 returns m −2 , standard deviation 1.9 returns m −2 ). A 50cm digital elevation model (DEM) and canopy height model (CHM) were built using Fusion v2.6 software (USDA Forest Service Remote Sensing Applications Centre, Salt Lake City, Utah). Individual trees were detected from the canopy height model using TreeVaW v1.2 software (Popescu et al., 2003) which identifies high points in the CHM that are likely to be tree tops using local  18  maximum filtering, with a circular window sized proportionally with the relationship between tree height and crown diameter (Popescu and Wynne, 2004; Popescu et al., 2002). To parameterize TreeVaW for Douglas-fir, a regression model was developed from live tree height and crown diameter data available in Coops et al. (2007).  2.2.3 Lidar Aboveground Mass Estimation To estimate TSAM, two sets of predictor variables were used: 1) metrics from lidar returns within the plot, and metrics from measurements of individually identified trees. Lidar metrics calculated for each plot included height of percentiles of nonground returns (10, 20, 30, 40, 50, 60, 70, 80, 90, and 95 percentiles, P10 to P95 respectively); and 2) the canopy return density over a range of relative heights ie. proportion of laser returns in ten even height intervals (D0, D1, D2, D3, D4, D5, D6, D7, D8, D9) between 0.5 m (D0) and the 95th (D9) percentile of maximum height (Næsset, 2004; Næsset and Gobakken, 2008). Summary descriptive statistics for the TreeVaW identified trees (TVT) included maximum height (MaxHt), minimum height (MinHt), and variance in heights (HtCV), and stem density (Den) for trees ≥ 28 m height in each ground plot (the mean height of trees with codominant canopy status). All predictor and target variables were transformed using the natural logarithm to improve linearity. A linear multiple regression model was developed based on maximum r 2 selection. Variables were not transformed any further (eg. Principal Component Analysis was not used), making explanation of lidar measurements in terms of stand structure more straightforward. A  19  correlation matrix was used to test for models with combinations of variables with high collinearity. To reduce the bias often present in log-transformed allometric equations, a correction factor was applied to all variables when back-transforming to original units (following Sprugel, 1983). All statistical analysis was performed using the Statistical Analysis System (SAS) v9 (SAS Institute, Richmond Va.).  2.3 Results The relationship between stem locations identified by TreeVaW and live stem locations measured in the field is shown in Table 2.1. Initially a comparison of all trees, regardless of height, was undertaken. Secondly, a minimum height cutoff of 28 m was applied to both TreeVaW identified trees and live inventory measured trees and the relationship re-assessed. As expected, the results indicate that many of the smaller trees, lower in the canopy were not consistently detected using TreeVaW, due to the dense overstory canopy structure. Also, TreeVaW parameterization uses only live trees and only three of 117 dead standing trees were correctly identified. Restricting our analysis to live dominant and co-dominant trees (greater than 28 m) resulted in a greater proportion of correctly identified trees.  A strong linear relationship ( r 2 = 0.92 and SE E = 0.69m ) was observed for correctly identified trees greater than 28 m high between tree height identified by TreeVaW and height measured in ground plots (Figure 2.1), though the TreeVaW predicted heights deviated from a 1:1 line with a constant bias such that TVT heights were about 1 m less than measured GP tree heights.  20  Table 2.1. Matches between TreeVaW identified trees (TVT) and ground plot measured trees (GP). Matches occur when the crown radius of a tree located by TreeVaW encompasses a stem mapped tree in a ground plot. GP  1GP:1TVT  Many GP:1 TVT  Unmatched GP  0GP:1TVT  (# of stems)  (% of GP stems)  (% of GP stems)  (% of GP stems)  (# of unmatched TVT stems)  All trees  302  12  38  50  17  Trees ≥ 28m  105  41  27  32  7  Figure 2.1. For trees greater than 28 m high, tree heights measured using TreeVaW matched tree heights measured in ground plots with an r 2 = 0.92 and SE E = 0.69m , prediction 95% confidence interval in grey, 1:1 line in grey dots.  Correlations between TSAM and predictor variables from lidar (Figure 2.2A) and TreeVaW (Figure 2.2B) showed that the lidar return density at 60% maximum height ( D6 ) had the highest correlation with TSAM and that correlations with all TreeVaW variables were low.  21  a)  b)  Figure 2.2. a) Pearson correlation coefficients between: 1) TSAM and lidar plot metrics (height percentiles P10-P95 as black line with circles; height-densities D0D9 as black line with squares), and 2) height percentiles (P0-P10) and D6 (the most highly correlated plot-level metric) as grey line with diamonds. b) Pearson's correlation coefficient between: 1) TSAM and TreeVaW metrics as grey bars, and 2) between TreeVaW metrics and D6 (most highly correlated plot-level metric) as grey lines.  The two most applicable multiple linear regression models (highest r 2 and lowest SEE ) used either the best plot-level lidar variables alone (Model 1) or with the best TreeVaW plot summary individual tree variables in combination with plotlevel lidar variables (Model 2) (Table 2.2). For both models, a variable representing lidar canopy return density at 60% maximum plot height ( D6 ) was chosen. For the lidar plot-level model, a variable representing crown height was also chosen (P80, the 80th percentile of the height of non-ground returns). The Pearson's correlation coefficient between the two variables in the model was 0.49, indicating that collinearity was not a problem for this model. A plot of the Model 1 TSAM estimates demonstrated it deviated from the 1:1 line with variable bias such that it overestimated TSAM at low levels (Figure 2.3A). For model 2, in addition to D6, the second variable selected was TVTMaxHt (maximum height of TreeVaW  22  identified trees in each ground plot). The Pearson's correlation coefficient between the two variables in the model was 0.51, indicating that collinearity was not a problem for this model. A plot of the Model 2 predicted TSAM demonstrated it deviated from the 1:1 line with variable bias such that it overestimated TSAM at low levels (Figure 2.3A). The variable bias in both models is likely the result of the use of log transformed data for the regression.  Table 2.2. Multiple linear regression models for large tree and snag aboveground mass (TSAM). Models take the form ln Y = β 0 + β 1 * ln x1 + β 2 * ln x 2 + e . Where Y = AboveGroundMass . β0 β1 x1 β2 x2 SE E′ r2 Model 1 8.91**** 1.14***  D6  -0.86**  P80  0.75 29.68  Model 2 8.3***  D6  -0.64*  TVTMaxHt  0.72 31.95  1.09***  r 2 is the coefficient of determination; SE E′ is the standard error of the estimate in Mg/ha TSAM; Significance Level: * (<0.2); **(<.1); ***(<.05); ****(<0.001).  23  a)  b)  Figure 2.3. lidar multiple regression model estimates of large tree and snag aboveground mass (TSAM) versus ground plot measurements using (a) Model 1 or (b) Model 2, see Table 2.2 for model variables.  2.4 Discussion TreeVaW software was used to identify the locations and heights of live trees from lidar data collected over a Douglas-fir dominated stand on eastern Vancouver Island, BC Canada. The detection algorithm was initially developed for tree identification in plantation-like stands of loblolly pine (Pinus taeda) in the southeastern United States, often using high return-density lidar data. In these complex, multilayered stands of the Pacific Northwest, multiple trees often compose canopy units which are challenging to separate, resulting in multiple stems falling within a single lidar delimited crown. As well, in this stand, standing dead trees contribute more than 10% of the overall large tree mass; however, TreeVaW is not parameterized to measure these and few were accurately detected. When our analysis was restricted to taller live trees in the canopy, the TreeVaW algorithm was much more successful in identifying individual trees, and height was determined with considerable accuracy.  24  Measures of canopy return density were more successful than measures of height for predicting TSAM. On its own, the proportion of lidar returns at 60% of maximum tree height and higher ( D6 ) explained 65% of the variation in TSAM. Height percentile measures were less successful individually, with lower percentiles explaining more of the variation in TSAM, and were highly correlated with measures of canopy density. As a result, multiple linear regression models with two variables were examined in an attempt to explain more of the variation in TSAM.  Two multiple linear regression models for plot-level TSAM were developed. The first model selected predictor variables for lidar return height ( P80 ) and return density ( D6 ). The second model selected predictor variables for TreeVaW maximum height ( TVTMaxHt ) and lidar return density ( D6 ) and explained 72% of the variance. The 80th percentile return height performed slightly better than TVTMaxHt , explaining 75% of the variance. For both models canopy return density  was the most significant variable for determining TSAM with a less significant relationship between TVTMaxHt in plot and TSAM. For this study, individual tree measurements did not improve TSAM estimates. This is most likely due to a number of the TSAM components not being well captured using individual treebased approaches. A study by Falkowski et al. (2008) shows that in conifer forests with high canopy cover, individual tree detection algorithms commonly make errors of omission, such as missing subdominant stems or clumping multiple stems  25  together into single canopy objects. As a result of these trees not being measured properly, the accuracy of TSAM estimations in conifer forest with high canopy cover may be limited. For example, in this study, TSAM is distributed throughout the canopy in dead standing trees (not measured by TreeVaW), closely spaced shoulder height trees (as indicated by the success of the P80 lidar metric), and very tall live standing dominant and codominant trees. Individual tree-based approaches may be more successful for plot-level estimates of TSAM in stands where canopy units are more distinct, such as in conifer forests with less dense canopy cover. Given our success in estimating TSAM despite difficulties detecting individual suppressed or dead standing trees, for other studies seeking to estimate TSAM in stands with complex canopies, it may be more economical to collect lower-density lidar data and utilize plot-level lidar return metrics.  This study demonstrates that it is practical to apply published techniques for individual tree detection and canopy return density using small footprint discrete return lidar data and use the results to estimate large tree and snag aboveground mass.  26  2.5 References Andersen, H., Reutebuch, S., and McGaughey, R. 2006. A rigorous assessment of tree height measurements obtained using airborne lidar and conventional field methods. Canadian Journal of Remote Sensing. Vol 32: pp. 355-366. CCP. 2008. Canadian Carbon Program/Programme Canadien du Carbone. Retrieved October 3, 2008, from http://www.fluxnetcanada.ca/home.php?page=home&setLang=en. Chapin, F., Woodwell, G., Randerson, J., Rastetter, E., Lovett, G., Baldocchi, D., Clark, D., Harmon, M., Schimel, D., Valentini, R., Wirth, C., Aber, J., Cole, J., Goulden, M., Harden, J., Heimann, M., Howarth, R., Matson, P., McGuire, A., Melillo, J., Mooney, H., Neff, J., Houghton R., Pace, M., Ryan, M., Running, S., Sala, O., Schlesinger, W., and Schulze, E.. 2006. Reconciling carbon-cycle concepts, terminology, and methods. Ecosystems. Vol. 9. pp. 1041-1050. Coops, N., Hilker, T., Wulder, M., St-Onge, B., Newnham, G., Siggins, A., and Trofymow, JAT. 2007. Estimating canopy structure of Douglas-fir forest stands from discrete-return lidar. Trees - Structure and Function. Vol. 21. pp. 295-310. Falkowski, M. J., Alistair, M.S., Gessler, P.E., Hudak, A.T., Vierling, L.A., and Jeffrey, S.E. 2008. The influence of conifer forest canopy cover on the accuracy of two individual tree measurement algorithms using lidar data. Canadian Journal of Remote Sensing Vol. 34. pp 338-351. Gillis, M., Omule, A., and Brierley, T. 2005. Monitoring Canada's forests: The National Forest Inventory. Forestry Chronicle. Vol. 81. pp. 214-221. Humphreys, E. R., Black, T. A., Morgenstern, K., Cai, T., Drewitt, G. B., Nesic, Z., and Trofymow, JAT. 2006. Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of development after clearcut harvesting. Agricultural and Forest Meteorology. Vol. 140 pp. 6-22. Landsberg, J., and Waring, R. 1997. A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. Forest Ecology and Management. Vol. 95. pp. 209-228. Lefsky, M. A., Cohen, W. B., Acker, S.A., Parker, G.G., Spies, T.A., and Harding, D.. 1999. Lidar remote sensing of the canopy structure and biophysical properties of Douglas-fir western hemlock forests. Remote Sensing of Environment. Vol. 70 pp. 339-361.  27  Næsset, E. 2004. Practical large-scale forest stand inventory using a small-footprint airborne scanning laser. Scandinavian Journal of Forest Research. Vol 19. pp. 164-179. Næsset, E., and Gobakken, T. 2008. Estimation of above- and below-ground biomass across regions of the boreal forest zone using airborne laser. Remote Sensing of Environment. Vol. 112. pp. 3079-3090. NFI. 2004. Canada's National Forest Inventory - Ground Sampling Guidelines v4.1. Retrieved April 10, 2008, from https://nfi.nfis.org/documentation/ground_plot/Gp_guidelines_v4.1.pdf Persson, A., Holmgren, J., and Soderman, U. 2002. Detecting and measuring individual trees using an airborne laser scanner. Photogrammetric Engineering and Remote Sensing. Vol. 68. pp. 925-932. Pojar, J., Klinka, K., and Demarchi, D. 1991. Chapter 6: Coastal Western Hemlock Zone. In: Meidinger, D., Pojar, J. (Eds.). Ecosystems of British Columbia. BC Special Report Series. pp. 95-111 . Research Branch, Ministry of Forests. Popescu, S. 2007. Estimating biomass of individual pine trees using airborne lidar. Biomass and Bioenergy. Vol. 31. pp. 646-655. Popescu, S., and Wynne, R. 2004. Seeing the trees in the forest: Using lidar and multispectral data fusion with local filtering and variable window size for estimating tree height. Photogrammetric Engineering and Remote Sensing. Vol. 70. pp. 589-604. Popescu, S., Wynne, R., and Nelson, R. 2003. Measuring individual tree crown diameter with lidar and assessing its influence on estimating forest volume and biomass. Canadian Journal of Remote Sensing. Vol. 29. pp. 564-577. Popescu, S., Wynne, R., and Nelson, R. 2002. Estimating plot-level tree heights with lidar: local filtering with a canopy-height based variable window size. Computers and Electronics in Agriculture. Vol. 37. pp. 71-95. Sprugel, D. G. 1983. Correcting for bias in log-transformed allometric equations. Ecology Vol. 64. pp209-210. Wulder, M. A., Bater, C. W., Coops, N. C., Hilker, T., and White, J. C. 2008. The role of lidar in sustainable forest management. The Forestry Chronicle. Vol. 84. pp. 807-826.  28  3. Determination of Carbon Stock Distributions in the Flux Footprint of an Eddy-Covariance Tower in a Coastal Forest in British Columbia2 3.1 Introduction Global carbon (C) cycling is strongly influenced by forest C dynamics as forests cover 4.1 billion hectares of the Earth's surface and are estimated to contain 1146 Pg C (Dixon et al. 1994). As a result, management and opperations decisions, such as harvest and reforestation (Kurz et al. 2002), permanent deforestation (Houghton, 1999), afforestation (Nilsson and Schopfhauser, 1995), pest management (Kurz et al., 2008), and fire management (Amiro et al., 2001) can have potentially global  climate impacts (Houghton, 2003; IPCC, 2007). A combination of field-based observations and C budget modelling (Houghton et al., 1999; Kurz and Apps, 1999; Landsberg and Waring, 1997) is critical to meet national and international reporting requirements, inform policy decisions, and potentially mitigate anthropogenic impacts on climate (Kurz and Apps, 2006).  In recent years, in an effort to understand C fluxes between forest ecosystems and the atmosphere, a global network of eddy-covariance (EC) flux-towers has been established at approximately 400 sites across a broad range of ecosystems (Baldocchi, 2008). The EC flux methodology uses micrometeorological equipment to provide continuous ecosystem-scale measurements of CO2 exchange directly at the canopy-atmosphere interface over time scales of hours, days, and years. The  2  A version of this chapter has been submitted for publication. Ferster, C. J., Trofymow, J.A.(Tony), Coops, N.C., Chen, B., Black, T.A., and Gougeon, F.A. Determination of carbon stock distributions in the flux footprint of an eddy-covariance tower in a coastal forest in British Columbia.  29  tower footprint, the source area where fluxes measured by the tower originate, depends on wind speed and direction, surface roughness, and instrument height. The flux-tower footprint can be modelled using simple geometry; however, advances in footprint analysis have resulted in more representative footprint shapes and surfaces (Finnigan, 2004; Schmid, 2002). Footprints may be modelled on a half-hourly basis corresponding with EC flux measurements and accumulated over long time periods to represent the long-term average EC flux probability density (Chen et al., 2009).  At an ideal EC tower site, site conditions meet the assumptions of horizontal homogeneity, steady state, and non-advective flows; however, many forests have complex structures and patterns of vegetation, which violate these assumptions (Göckede et al., 2004). For example, in coastal British Columbia, forests are often in areas of complex topography and affected by small-area disturbances, such as windthrow, mass movement (such as landslides), and insect attack, resulting in a spatially heterogeneous stands with uneven distributions of species, age, diameter, canopy layering, canopy gaps, and understory vegetation (MacKinnon, 2003; Panel, 1995). Upslope (anabatic) and downslope (katabatic) winds often occur in sloping terrain. Finally, landscape position near large water bodies may affect daily weather patterns related to temperature differences between land surfaces and water bodies, causing sea breeze (daily onshore) and land breeze (nightly offshore) flows. In recent years, this global network of EC flux-towers is being used to answer important research questions; therefore, proper interpretation and understanding of EC flux-tower data depends on careful analysis of flux-footprints  30  and application of the footprint analysis to determine attributes of the forest within the flux-footprint.  Traditional biometric measurement of landscape C stocks allows us to directly account for changes in landscape C stocks over time and indirectly estimate fluxes to the atmosphere (Clark et al., 2001; Houghton, 2005). Numerous techniques, such as forest inventory (eg. Gillis et al., 2005), are available to measure landscape C stocks. Stands are sampled with one or more measured ground plots to determine various biometric attributes which can then be used to calculate live biomass and detritus mass components. Tree biomass is estimated by measuring attributes of individual trees and applying allometric equations, derived from destructive samples of representative trees (Brown, 2002). The mass of smaller understory shrubs, herbs, and bryophytes, as well as fine woody debris and forest floor material are typically determined through a destructive sample of one or more microplots (within the main plot) which is then dried, and weighed (NFI 2004A) and converted to a C basis either through laboratory determination of sample %C, or more typically through the assumption that biomass is composed of 50% C (Litton et al. 2007; Birdsey, 1992). Soil C stocks are typically determined through sampling one or more soil pits or soil cores to a fixed depth to determine soil type and coarse fragments and obtain samples of assessment of bulk density, texture, and the %C of fine soil and detritus fractions (NFI 2004A; Palmer et al. 2002). All of the individual estimates are combined to provide the overall value of ecosystem C stocks for a ground plot and then extrapolated in some manner to provide a stand level value.  31  Traditional forestry applications often employ a stratified sampling approach where environmental predictor attributes, potentially from expert interpretation of aerial photography, are used to define homogenous polygonal strata boundaries. Representative biometric measurements are taken from within the boundaries of each strata and the results are extrapolated to the full spatial extent of the area under study. This approach is advantageous because it reduces the redundancy in simple random sampling by strategically locating sample points (Schreuder et al., 2004). Modern computer technology has improved on this approach, now allowing multivariate extrapolation and quantification of the homogeneity within strata using spatial data environmental indicator coverages available from remote sensing and GIS forest land management data (Hudak et al. 2008; Crookston and Finley, 2008; Moeur and Stage, 1995).  A recent remote sensing technology, discrete return light detection and ranging (lidar), provides another effective way to measure tree biomass (Næsset and Gobakken, 2008; Wulder et al., 2009). Lidar instruments are mounted on an airborne platform with a very accurate global positioning system and inertial measurement system. Lidar sensors emit near infrared light at as many as 160,000 pulses per second in a scanning pattern, with footprints 20 - 30 cm in diameter spaced 1 m apart. The time taken for up to five discrete returns to reflect from the forest canopy, stems, branches, trunks, and the ground is measured to determine height. Initial returns originate from high in the canopy and late returns originate from the ground below the canopy (Wulder et al., 2008). Many studies have shown  32  lidar to be a robust tool for measuring tree biomass in forest landscapes (Lefsky et al., 2005; Næsset, 2004; Næsset and Gobakken, 2008). Furthermore, a very  accurate digital terrain model (DTM) is developed in the process, providing information about underlying landforms.  In this paper, we first compare two approaches to estimate the spatial distribution of aboveground mass in large trees and snags (greater than 9.0 cm diameter at breast height) (TSAM) for a EC flux-tower site: first, using a GIS based stratification method and second, using lidar remote sensing data based estimate. We then calculate and compare the site-level average estimate of TSAM for the two methods based on absolute spatial variability and flux-footprint-weighted variability. The lidar remote sensing estimate thus serves as an independent method to compare with the GIS based stratification estimate and establish how well we can scale from ground plots using the GIS based stratification technique. Once we are confident with the results for TSAM, we then use the GIS based stratification method to estimate spatial distribution of other C stock components at the site, including understory, detritus and soil - which are difficult to measure using remote sensing techniques, and derive a site-level average estimate of each measured component and total aboveground C stocks to the full extent of the footprint at the EC flux-tower site.  33  3.2 Materials and Methods 3.2.1 Study Site The Canadian Carbon Program (CCP)- Fluxnet Canada Research Network (FCRN), British Columbia Flux Station on east central Vancouver Island includes three EC flux-tower sites and other monitoring sites at which short-term and longterm measurements for ecosystem C balance are made, supporting several C budget modelling projects (Humphreys et al., 2006; Trofymow et al., 2009). This study examines one of these sites, DF1949, at which an EC flux-tower (49º52´7.8"N, 125º20´6.3"W) (Figure 3.1) was established in 1997 and has been making continuous measurements since establishment (Humphreys et al., 2006). The main stand at the DF1949 site is dominated by Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco) with lesser amounts of western hemlock (Tsuga heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata Donn ex D. Donn), and  red alder (Alnus rubra Bong.). This second-growth stand was established following harvest (1937, 1938, and 1943) of the original old-growth stand and subsequent broadcast burning (1939 and 1943). Due to limited regeneration, most of the site was planted with Douglas-fir in 1949, however stands in productive areas on the periphery of the site established naturally or were planted between 1942 and 1980. Forest inventory mapping by forest companies operating in the area show the site contains a few stands of different forest cover type and/or site index classes. However, differences in ecosite series due to topography (260 to 470 m elevation, slope 0 to 27˚) and disturbance history have resulted in finer scale variation in forest cover than that in the inventory mapping classes (Trofymow et al., 2008). Soils on the site are predominantly duric humo-ferric of morainal origin of over 1  34  m in depth and of a gravelly sandy loam texture (Jungen 1985). However a change in slope at the site means soils on the low slope eastern half are generally wetter than those on the steeper slopes of the western half.  Figure 3.1. General location of CCP/FCRN BC Flux Station study area .  3.2.2 GIS Data A library of GIS stratifying variables is available for the BC Flux Station Oyster River area which includes the DF1949 site (Trofymow et al., 2008) (Table 3.1). Operational forest inventory (forest cover species, site index, and date of establishment) and terrestrial ecosystem mapping of site series, the near seral climax vegetation type (Banner et al., 1996), was provided by Timberwest and Weyerhaeuser (now Island Timberlands), the two forest companies operating in the area. More finely detailed ecosystem site series mapping was collected by the Canadian Forest Service (CFS) near the tower, within the EC tower footprint boundaries as originally estimated in 2002 (Figure 3.2), using field transects and interpretation of aerial imagery acquired in 1999. The two site series coverages  35  were combined so that fine resolution site series mapping is used where it is available, and coarse mapping of site series is used elsewhere. As part of ongoing research at the study area, a comprehensive disturbance database was previously developed using historic maps and historic aerial imagery (Trofymow et al., 2008). Elevation was obtained from a digital terrain model (DTM) developed by Ferster et al. (2009) from lidar imagery acquired in 2004 (fully described in the next section)  and slope (in degrees), aspect, and topographic variables (following Hudak et al. 2008; Roberts and Cooper 1989) were calculated from the DTM. A 1999 orthophoto with 1 m x 1 m pixels was used to calculate individual tree crown (ITC) locations and subsequently calculate dominant canopy tree density. The ITC algorithm uses a valley-following approach to delimit tree crowns and identify bright pixels that are peaks in the canopy, and it is well suited to the conifer forests of Vancouver Island (Gougeon, 1995). Finally, a Quickbird high spatial resolution image (60 cm panchromatic, 2.4 m multispectral Blue = 450-520 nm, Green = 520600 nm, Red = 630-690 nm, and Near Infrared (NIR) = 760-900 nm) was acquired in 2004, and used to calculate the normalized difference vegetation index (NDVI) ( ( NIR − RED ) /( NIR + RED ) ).  36  Table 3.1. GIS and Remote Sensing data, variables, data sources, and suitability for use in most similar neighbour classification. Coverage  Variable  Data Level  Units  Forest Inventory  Site Index1  ratio  m  Disturbance History (Trofymow et al. 2008)  2004 Lidar Survey (Coops et al. 2007)  CFS and Forest Companies (Trofymow et al.2008)  37  29  36  Suitable for model yes  Top Height  ratio  m  0  31  24.7  28.7  yes  2nd Harvest  interval interval  year year  1937 0  1943 2008  1938 0  1943 0  yes no  1st Fire  interval  year  1939  1943  1939  1943  yes  2nd Fire Date Est. 1953  interval interval  year year  1943 1940  1943 1980  0 1945  0 1949  no yes  Date Est. 2003  interval  year  1940  1980  1945  1949  yes  1st Fertilization  interval  year  0  1994  1994  1994  no  2nd Fertilization  interval  year  0  0  0  0  no  359. 96 469. 75 27.1 6 27.0 5 25.9 5  13.7 8  314.25  yes  0.00  387.07  yes  1.48  10.48  yes  -7.73  7.55  yes  aspect  ratio  azimuth  Elevation  ratio  m asl  Slope  ratio  degrees  2  flat 262.44 0.00  ratio  -26.59  SSINA2  ratio  -23.23  -8.77  9.10  yes  TSRAI2  ratio ratio  stems/ha  0.50 57.96  0.51 254. 15  1.49 381.67  yes yes  Dominant Canopy Tree Density NDVI  1.50 504. 82  ratio  NDVI  0.5  0.78  0.64  0.7  yes  Cover Species3  nominal  Site Species3  nominal  Dr, Fd, NP  Fd  no  Fire Cause 14 Fire Cause 24  nominal nominal  M-PB, SB SB  SB null  no no  Most Common Soil Association  nominal  Quimper (morainal), Piggot (colluvial)  Quimper (morainal)  no  Site Series5  nominal  1,3,4,5,7,12  1, 5, 12  yes, as wet/dry type  Forest Inventory (Trofymow et al.2008)  Disturbance History (Trofymow et al. 2008) Soil Survey of Canada (Jungen, 1985)  15  Ground Plots MIN MAX  1st Harvest  SCOSA  1999 Orthophoto (Gougeon, 1995) 2004 Multispectral (Quickbird) Nominal Variables  Footprint MIN MAX  Unique Values CwFdHw, Dr, DrCwFd, Fd, FdCw, FdCwHw, FdDr, FdHwBa, FdHwCw, NP, SW  Unique Values Fd, FdCw  no  37  Footnotes for Table 3.1: 1  Site Index: Tree height (m) at 50 years age. Topographic Variables: SCOSA = Slope * cos (Aspect), SSINA = Slope * sin (Aspect); TSRAI (Topographic Solar Radiation Aspect Index) = (1cos((pi/180)(Aspect-30)))/2 (Roberts and Cooper, 1989) 3 Species codes for Cover Species and Site Species: Fd (Pseudotsuga menziesii var menziesii), Dr (Alnus Rubra), Cw (Thuja plicata), Hw (Tsuga heterophylla), Sw (Swamp), and NP (Not Productive) 4 Fire cause codes for Fire Cause 1 and Fire Cause 2: M-PB (escaped broadcast burn into standing timber, partial disturbance) and SB (broadcast slash burn). 5 Site Series codes: 1 (Western Hemlock/Douglas-fir – Kindebergia), 3 (Douglasfir/Western Hemlock – Salal), 4 (Douglas-fir – Swordfern), 5 (Western Red Cedar – Sword Fern), 7 (Western Red Cedar – Foamflower), 12 (Western Red Cedar – Skunk Cabbage). Dry site series (1, 3, 4, and 5). Wet site series (7 and 12). 2  38  3.2.3 Ground plot Inventory Data In 2002 and 2006, ecological ground plots were established and measured by the CFS following National Forest Inventory guidelines (NFI 2004a). Ground plots were established in four clusters of three individual ground plots, in a triangular arrangement 30 m appart, resulting in 12 individual inventory plots systematically located to sample the range of mapped forest cover and ecosite series variation within the forest area contained by the estimated flux-tower footprint in 2002 (Humphreys et al. 2006) (Figure 3.2). Ground plots were measured for large trees (>9.0 cm diameter at breast height, 130 cm above ground; measurements of height, DBH, species, status, stem location), small trees (<9.0 cm DBH, >1.3 m height; measurements of height, DBH, species, status), understory vegetation and fine woody debris (in four 1 m2 microplots), woody debris (four 15 m transects), forest floor (four 20 cm by 20 cm samples and depth measurements along the four 15 m transects), and soil (four 10 cm diameter excavation cores 0 – 15 cm, three 15 – 35 cm, one 35 – 55 cm used for bulk density and soil chemistry, and one 0 – 55 cm depth 50 x 50 cm soil profile pit used for pedogenic descriptions and coarse fragment determination). Field and laboratory measurement data were formatted to the National Forest Inventory standards and submitted to the NFI database for data compilation (e.g. calculation of biomass from NFI allometric equations for tree species, DBH and height) (NFI, 2004b) and total and component ecosystem C stocks summarized per plot.  3.2.4 GIS Stratification Units Derivation The approach of this study to delineate GIS stratification units and estimate the spatial distribution of target biometric variables is explained as follows:  39  In the first step, GIS and remote sensing predictor variables were selected from the full suite of available data based on their predictive capability for the target biometric variables using the following criteria: 1.  The range of values sampled by the ground plots representative of the range  and spatial scale of values observed across the footprint extent. 2.  A statistically significant relationship with the target variables as assessed  using a correlation matrix with Pearson’s product-moment coefficient (R) (for continuous variables) and the point biserial coefficient (for nominal data). 3.  Unique information (no collinearity with other variables).  Groups of variables with statistically significant R or point biserial coefficients (at  α = 0.1) were ranked according to the strength of the relationship and then examined for collinearity with other predictor variables, indicated by an R value 0.9 or higher (following Hudak et al., 2008). Where collinearity existed within groups of predictor variables, the variable with the highest R value for the target variables was selected for the model and the other variables removed to minimize redundancy.  In the next step, the selected predictor variables were then used in a nearest neighbour approach to assign target biometric variables from representative ground plots to each footprint cell. Values for each cell were calculated using Most Similar Neighbour (v2) software (USDA, Moscow Id.) measured with the Mahalanobis distance (Mahalanobis, 1936), a statistical distance based on a transformation of  40  the covariance matrix of the predictor variables. The measured biometric attributes from the representative ground plot clusters are then applied to each footprint cell.  Finally to evaluate the representativeness of the ground plots within the footprint extent, the Mahalanobis distance was mapped to illustrate the spatial representativeness of the ground plots for the predictor variables.  3.2.5 Lidar Data Discrete return light detection and ranging (lidar) data with mean beam footprint size of 19 cm and average 5.3 non-ground returns per square meter return density was collected in 2004 by Terra Systems (Sydney, British Columbia) (Coops et al. 2007). Ground returns were separated from non-ground returns by the lidar vendor, a digital terrain model was built using Fusion v2.0 software (USDA Forest Service), and height above ground of non-ground returns was determined by subtracting non-ground return height from the corresponding ground height at the same location on the DEM. Further processing of lidar data was done using SAS software (Richmond, Virginia).  3.2.6 Lidar Estimates of TSAM Spatial Distribution Lidar data provides an alternative approach to measure above ground tree mass from living and dead trees across the landscape at the study site. Ferster et. al (2009) calculated TSAM for the 10 m footprint cells at the site using metrics for return heights, and density of return heights at regular intervals through the canopy,  41  following Næsset and Gobakken (2008). These lidar metrics include height of percentiles of non-ground returns (10, 20, 30, 40, 50, 60, 70, 80, 90, 95, P10 – P95, respectively) and canopy return density ie. proportion of returns at 10 regular intervals (D0, D1, D2, D3, D4, D5, D6, D7, D8, D9) between 0.5 m (minimum nonground height, D0) and the height of the 95 percentile (maximum height, D9). All variables were log transformed, to improve linearity, a maximum r2 selection was used to find combinations of variables that were highly correlated with TSAM, and a correlation matrix was used to reject models with a high amount of collinearity. A multiple linear regression model was built to estimate biomass using one lidar metric of percentiles of return height, and one lidar metric of return density (ln TSAM = 8.91 + 1.14 * ln D6 + -0.86 * P80). Finally, to reduce the bias often present in log-transformed biometric allometric equations, a correction factor was applied when back transforming to original units (following Sprugel, 1983).  3.2.7 EC Flux-Tower Footprint A detailed footprint analysis for the DF1949 site was completed by Chen et al., (2009) using 30-minute EC tower-measured net ecosystem productivity (NEP) flux-footprints averaged over a five-year period from 2002 - 2006 (Figure 3.2). Footprint data for 10 m by 10 m cells for probability density and cumulative probability density were spatially referenced for this analysis. The flux density from each cell rapidly diminishes with distance from the tower and the 99% fluxfootprint is spread over a very large area (816 ha) extending beyond the bounds of the GIS and lidar data, therefore, it was necessary to choose a cut-off value that  42  encompassed between 80 and 90% of the cumulative flux. For this study, the 85% cumulative flux isoline was chosen as the footprint boundary for analysis and the individual cell probability densities were recalculated by renormalizing to the 85% footprint area (162 ha).  Figure 3.2. Mean 5-year NEP flux-footprint cumulative isolines and 85% flux boundary at the DF1949 site (from Chen et. al 2009) and original (2002) estimated footprint boundary . Background image 2004 Quickbird panchromatic.  43  3.2.8 Comparing Lidar and GIS Stratification Estimates of TSAM To assess the difference between the GIS stratification and lidar large tree above ground mass estimates of TSAM, comparisons were made of the spatial distribution of differences within the footprint, the stratified and unweighted means, and the stratified and unweighted means. To compare the spatial distribution of high and low values, all footprint cells for both methods were normalized to the mean of each respective surface by dividing the cell value by the mean surface value for that method. Subsequently, corresponding pixels for each method were subtracted from each other. To determine the stratified and unweighted site-level mean, all footprint cells for both estimates were averaged, and the standard deviations calculated. To determine the stratified and footprint weighted mean, all footprint cells for both estimates were averaged with weighting by the footprint NEP flux probability density, and the unbiased weighted sample standard deviation was calculated. Root mean square difference (RMSD) (Stage and Crookston, 2007) was used to calculate error for each variable associated with the model. For the GIS stratification based estimate, a cross validation procedure was utilized, due to the nearest neighbour classification resulting in no residuals (each plot is its own nearest neighbour). Following this cross validation approach, one ground plot was removed from the model at a time and the biometric values imputed from the pool of other ground plots, and these differences were used to calculate the RMSD. For the lidar estimation, the standard error of the estimate (SEE) was calculated for the residuals of the ground plots used to build the model.  44  3.2.9 GIS Stratification Estimates of Other Component C Stocks Using the nearest neighbour classification described above for the GIS based stratification units, C stocks for components of live biomass carbon, detritus mass carbon, and mineral soil carbon were extrapolated for each footprint cell for the full extent of the footprint. The stratified and unweighted, and stratified and footprint weighted mean and standard deviation, as well as model RMSD, was calculated for each component.  3.3 Results 3.3.1 GIS Stratification Units Seven GIS and remote sensing environmental predictor variables were selected to estimate the spatial distribution of TSAM and other C stock biometric variables within the flux-tower footprint (Table 3.2). Using the selected predictor variables, the nearest neighbour classifications for each footprint cell and representative ground plots are mapped in Figure 3.3a. The nearest neighbour Mahalanobis distance is mapped in Figure 3.3b, to demonstrate how well the footprint is represented by the ground plots spatially in this classification. Footprint cells with small Mahalanobis distances are relatively well represented by the ground plots, while cells with large Mahalanobis distances indicate that the sample plots have relatively dissimilar predictor attributes. The figure shows that near the tower, Mahalanobis distances were relatively small, with small patches with higher values, while the largest Mahalanobis distances were near the southeast and southwest corners of the footprint, areas where ground plots had not been installed, as these were beyond the area considered to be the tower footprint in 2002. The south east sector of the footprint is a shallow sloping bench upland of the tower,  45  incised by a stream channel with steep banks sloping southwest away from the tower that are classified as non-productive in the forest inventory GIS data. Areas in the southwest corner with large Mahalanobis distances are in a gently sloping area where water collects in wetlands and seasonal water channels. In addition, a band of higher Mahanobis distances is also present in the along the western edge of the footprint, an upland area far from the tower where regeneration occurred much later (1980) than the rest of the stand.  46  Table 3.2. Significant predictor variables for ground plots at DF1949 (rows) and Pearson product-moment correlation coefficients (R) for biometric target variables (columns). 1. Tree and snag above ground mass (TSAM) (Mg/ha), 2. Large and small live trees (Tree) (MgC/ha) 3. Understory shrubs, herbs, and bryophytes (U. St.) (MgC/ha); 4. Standing dead (Snag) (MgC/ha) 5. Woody debris (WD) (MgC/ha) 6. Stumps (Stmp) (MgC/ha) 7. Litter, fibric, and humus (LFH), (MgC/ha) 8. Fine woody debris (FWD) (MgC/ha) 9. Mineral soil 0-55 cm (Soil). (MgC/ha) 10. Maximum Pearson product-moment correlation coefficient (Max R). 11. Selection decision (Selection), selected variables are in bold, or explanation given for noninclusion. TSAM  Tree  U. St.  Snag WD Stmp LFH FWD Soil Max R  Site index  -0.15  -0.39  0.85  -0.2  Top height  -0.15  -0.39  0.85  -0.2  First Harvest  0.15  0.39  -0.85  0.2  Fire year 2  0.15  0.39  -0.85  0.2  Establishment 1953  -0.15  -0.39  0.85  -0.2  Establishment 2003  -0.15  -0.39  0.85  -0.2  Site Series (wet/dry) 1  -0.85  -0.72  0.09  -0.43  Elevation  0.15  0.35  -0.73  0.27  NDVI ITCDen Aspect Slope TSRAI  0.24 0.07 -0.33 0.15 -0.43  0.41 -0.20 -0.45 0.25 -0.53  -0.68 0.70 0.31 -0.38 0.16  0.2 0.15 -0.01 0.29 -0.05  Selection  -0.16 0.19 -0.64 0.29 0.44 0.85 Selected Collinear -0.16 0.19 -0.64 0.29 0.44 0.85 with Site Index Collinear 0.16 -0.19 0.64 -0.29 -0.44 0.85 with Site Index Collinear 0.16 -0.19 0.64 -0.29 -0.44 0.85 with Site Index Collinear -0.16 0.19 -0.64 0.29 0.44 0.85 with Site Index Collinear -0.16 0.19 -0.64 0.29 0.44 0.85 with Site Index -0.40 -0.58 0.19 -0.26 0.19 0.85 Selected Collinear 0.29 -0.19 0.65 -0.20 -0.44 0.73 with Site Index 0.11 -0.29 0.71 -0.21 -0.06 0.71 Selected 0.43 -0.03 -0.25 0.47 0.06 0.70 Selected -0.62 -0.16 -0.26 0.15 0.26 0.62 Selected 0.30 -0.31 0.40 0.13 -0.62 0.62 Selected -0.56 0.05 -0.01-0.02 0.58 0.58 Selected  Highlighted values are significant at p<0.1. Note that no predictors are significant for snags, but when combined with large live trees for TSAM, there are two significant predictors. 1 Point biserial coefficient used for nominal variable.  47  Figure 3.3. A) Nearest neighbour (NN) classification from selected environmental predictor variables, and spatial distribution for TSAM within footprint pixels. B) Mahalanobis distances (Dist.) (unit-less) for nearest neighbour classification of footprint pixels.  3.3.2 GIS Stratification and Lidar Estimates of TSAM The GIS stratification method estimate of TSAM from each ground plot extrapolated to the footprint extent is shown in Figure 3.4a, the lidar method estimate of TSAM is shown in Figure 3.4b, and the difference between the two estimates is shown in Figure 3.4c. Agreement is good near the tower, with scattered pixels both higher and lower. Lidar based estimates were higher than GIS stratification based estimates in wetland patches north of the tower. Lidar-based estimates were lower than GIS stratification based estimates in areas with canopy gaps and roadways. Large contiguous areas of disagreement are present at the far southwest and southeast parts of the footprint. In the southwest corner, the lidarbased estimate was consistently higher than the GIS stratification based estimate, except for an exposed stream bank, where lidar estimates were lower. In the southeast corner, the lidar estimates were greater than GIS stratification estimates in a large wetland area.  48  For TSAM, the unweighted average of 12 ground plots and unweighted and footprint density surface weighted means for both methods are listed in Table 3.3. The lowest estimate was the unweighted GIS stratification method estimate (282.93 Mg/ha, standard deviation = 35.21 Mg/ha, RMSD = 71.80), the unweighted mean of 12 ground plots was about 15 Mg/ha higher (298.42 Mg/ha, standard deviation = 64.75 Mg/ha), and the lidar method estimate was 37 Mg/ha higher (320.32 Mg/ha, standard deviation = 62.38 Mg/ha, SEE = 29.68). Using footprint probability density surface weighting increased the GIS stratification estimate of TSAM to 290.08 Mg/ha (standard deviation = 34.83 Mg/ha), 7.13 Mg/ha or 2.52% higher than the unweighted stratified mean. Using footprint density surface weighting decreased the lidar method estimate of TSAM to 318.80 Mg/ha (standard deviation = 45.21 Mg/ha), 1.52 Mg/ha or 0.47% lower than the unweighted stratified mean. The difference between standardized estimates of TSAM from GIS stratification and lidar methods is mapped in Figure 3.4. The footprint area contributing to 85% of the flux probability density is within +/- 25% relative agreement between the two methods, while 75% of the total footprint area is within +/- 25% agreement between the two methods of TSAM estimation. This indicates that the areas of the flux probability density surface contributing the most to long-term tower measurements of NEP are in better agreement than areas that are far from the tower and contribute less.  49  Figure 3.4. A) Large tree and snag above ground mass (TSAM) extrapolated from ground plots to footprint extent using the GIS stratification method (Strat.) (Mg/ha). B) TSAM estimated using lidar (Mg/ha). C) Spatial distribution of the difference between methods in estimates of TSAM (normalized to mean values).  50  Table 3.3. Tower site range of values, mean, and standard deviation (SD) for large tree and snag aboveground mass (TSAM) Mg/ha from 12 ground plots, Stratification (all footprint cells), and lidar (all footprint cells) methods, unweighted or weighted (W.) by the footprint flux probability density, and the percentage difference relative to unweighted mean. Standard error of the mean (SEM) is provided for the lidar model, and the root mean squared difference (RMSD) is provided for the Stratification model. Source  12 Ground Plots Lidar  Min  Max  Unweighted Mean (SD)  Weighted by Flux Probability Density Mean (SD)  Difference  124.45 8.83  371.48 691.69  298.42 (64.75) 320.32 (62.38)  N/A 318.80 (45.21)  -0.47%  124.45  371.48  282.93 (35.21)  290.08 (34.75)  2.52%  SEM = 29.68 Stratification RMSD = 71.80  3.3.3 Stratification Estimates of C Stock Components The components of live biomass, detritus mass, and mineral soil (0-55 cm) ecosystem C stocks, some of which can not be measured directly with remote sensing, were extrapolated to the footprint extent using the GIS stratification methodology (Figure 3.5 and Table 3.4). For total live biomass C stocks, the highest values are found in the south central part of the footprint where stem density is very high and in patches in the eastern part of the footprint. For total detritus mass C stocks, values are high in the south central part of the footprint, as well as in patches west of the tower. Total aboveground (live and detritus) C stock distributions follow the sum of the components, with consistently high values uphill of the tower to the west, and scattered very high values throughout the stand. Mineral soil (0-55 cm) C stocks are highest in flat areas east of the tower and also uphill of the tower to the west. For these C stock components, the greatest magnitude difference due to footprint weighting was 28.97 % higher than the stratified unweighted mean estimate for fine woody debris. While other C stock  51  components also changed due to footprint weighting, the total aboveground C stocks were only -0.12% lower than the unweighted stratified mean.  Table 3.4. Tower site range, means, standard deviations (SD), and model root mean squared difference (RMSD) for carbon stocks components in MgC/ha from stratification, unweighted or weighted (W.) by the footprint flux probability density and difference in weighted relative to unweighted mean. Min  Max  Unweighted Mean (SD)  Live Tree RMSD = 40.86 Understory shrub, herb, and briophytes RMSD = 0.66  69.81  188.97  0.29  2.08  132.42 (17.05) 1.69 (0.44)  Total Live RMSD = 40.53 Coarse Woody Debris RMSD = 29.86 Snags RMSD = 14.56 Stumps RMSD = 3.42 Litter, fibric, and humus RMSD = 17.49  71.23  190.35  8.91  Weighted by Flux Probability Density Mean (SD) 132.88 (14.90) 1.85 (0.33)  Difference  0.35% 10.00%  134.74 (14.74)  0.47%  67.08  134.11 (16.80) 44.09 (13.32)  41.70 (15.43)  -5.44%  2.66  55.75  18.53 (13.10)  20.47 (16.78)  10.51%  0.49  14.37  6.55 (1.60)  6.84 (1.39)  3.80  37.80  15.39 (7.98)  13.74 (8.59)  -10.73%  Fine Woody Debris RMSD = 1.31 Total Detritus RMSD = 46.08 Total Above Ground (Live + Detritus) RMSD = 70.79  1.96  4.81  3.20 (1.31)  4.13 (1.03)  28.97%  39.57  164.88  87.76 (29.05)  86.88 (37.49)  -1.01%  110.79  309.44  221.87 (34.28)  221.62 (37.77)  -0.12%  Mineral Soil (0 – 55 cm) RMSD = 44.42 Total Above Ground + Mineral Soil RMSD = 83.85  39.20  157.50  80.61 (29.92)  68.38 (31.72)  -15.17%  209.09  414.34  302.48 (49.75)  289.99 (63.06)  -4.13%  4.46%  52  Figure 3.5. Carbon (C) stock components (MgC/ha) extrapolated to the full footprint extent. A) Total Live (Live) B) Total Detritus (Detritus) C) Total Aboveground (Total Live + Total Detritus) (ABVG.) D) Mineral Soil (0-55 cm) (Soil)  3.4 Discussion Detailed analysis of EC flux-tower data will inform resource management decisions that are intended to limit the impact of anthropogenic activities on the global carbon cycle. The global EC flux databank holds records of C exchange over a broad range of regions, ecosystems, and time scales. Proper understanding of EC flux data records, however, is dependant on spatial scale, as EC flux-towers measure at the canopy scale, while biometric measurements are often made at the finer detailed plot scale. In this paper, we presented a methodology to scale plot-  53  level biometric measurements to the extent of an EC tower site footprint and then weight the spatial distribution based on mean five-year EC flux-tower NEP footprint probability density distribution, allowing comparisons to be made across scales over long time periods.  3.4.1 GIS Stratification Units Our results indicate that a range of environmental predictor variables can be used to stratify an EC flux-tower footprint and produce homogeneous areas represented by plots located strategically within the footprint. High spatial resolution Quickbird multispectral data was used to calculate NDVI, which was an important environmental predictor variable for this classification, and the spatial resolution (2.4 m by 2.4 m) was appropriate for comparison with lidar data. At sites where Quickbird data is not available, or is prohibitively expensive, multispectral Landsat data is freely available and is appropriate for calculating NDVI, though the spatial resolution is much coarser (30 m by 30 m).  Following nearest neighbour classification, the Mahalanobis distance was a useful tool to quantify and evaluate how well footprint cells are represented by the ground plots: small Mahalanobis distances represent good representation and large Mahalanobis distances indicate poor representation. This is demonstrated within the DF1949 site footprint, where outside of the area included in the original footprint estimation, there are some areas with large Mahalanobis distances, since few ground plots were established here. Therefore, this study is an advancement  54  over other studies, which seldom thoroughly assess how representative sample plots are of the study landscape.  3.4.2 GIS Stratification Estimate of TSAM The GIS based stratification was used to extrapolate above ground C stocks to the full extent of the EC flux-tower footprint. As expected, the GIS stratification estimate extrapolation technique maintained the range of values measured in the biometric ground plots, while more knowledge of the spatial distribution of values improved our estimate of the mean, especially where the ground plots are well representative of the footprint. In areas where ground plots are not representative of the ground cover, for example on roads or eroded stream banks, neither this technique nor the unweighted unstratified estimate yields accurate estimates of ground cover. As expected, the RMSD was high compared to the mean because it includes sources of error beyond the pure error of the model. For example, by applying a cross validation methodology, the difference between predictor observations (which were selected to be unique) is larger than the difference between target observations and representative predictor observations (Stage and Crookston 2007). Increasing the number of sample plots and creating sample plot redundancy in each ecosystem type would lower the RMSD for each estimate.  3.4.3 Lidar estimate of TSAM Lidar is a new tool available to make highly detailed measurements of stand structure (Wulder et al. 2008). Measurements from small footprint discrete return  55  lidar have been used to characterize vertical canopy distribution (Coops et al. 2007) and robust measurements of aboveground biomass from individual trees (Popescue, 2007) and at the plot-level (Naesset and Gobakken, 2008). In this study we use lidar to make a highly detailed and continuous estimate of the distribution of TSAM.  As expected, the lidar estimate captures more fine scaled details than the stratification-based estimate. As a result, fine positive and negative variations, related to canopy gaps and clusters of trees are scattered throughout the footprint in the TSAM estimation. In certain parts of the footprint, there are systematic trends due to measurement methods. For example, lidar estimates are low where there are canopy gaps on roads or exposed areas, and lidar estimates are high where there is a very closed canopy, though due to the intercept value and the nature of linear regression models, these estimate may be over or under- estimations of TSAM. As well, the lidar estimate includes a wider range of values than measured in the ground plots and we are less certain of values outside of this range. Finally, the Mahalanobis distance from the above section is also a useful way to assess the representativeness of the ground plots used to develop this lidar-based model.  3.4.4 Comparing Estimates of TSAM Agreement between the lidar and stratification based estimates was highest near the tower, where ground plots are most representative and EC flux-footprint probability density is highest. The areas which are in less agreement are near the  56  edges of the footprint, where EC flux-footprint probability density is lower. The agreement between the stratification and lidar estimates is higher with footprint weighting as compared to the unweighted estimate. This gives confidence to our use of GIS stratification to estimate the spatial distribution of plot-level biometric measurements to the extent of the EC flux-tower footprint.  3.4.5 Stratification Estimates of C Stock Components In order to describe the spatial distribution of C stock components other than TSAM, which can be measured directly using lidar, we utilized the GIS stratification methodology. The GIS stratification predictor variables that were chosen are correlated with multiple components of stand structure, and these correlations have a biologic basis. For example, site index is related to soil composition and nutrient availability (Carmean, 1975). Topographic variables, such as slope and aspect, are also linked to vegetation structure and distribution (McGarigal et al. 2009; Franklin, 1995; Granger and Schulze, 1977). Therefore, given these relationships, and our confirmation using lidar measurements of TSAM, these GIS variables are suitable to estimate the spatial distribution of C stock components that are difficult to measure using remote sensing.  Flux-tower weighting made a difference to the total C stocks, and changed the relative contribution of each component to the total. For example, footprint weighting increased above ground tree mass, while it decreased mineral soil C  57  stocks. Mineral C stocks were correlated with slope, where flat areas have accumulated high mineral soil C stocks.  Spatial heterogeneity within a flux-tower footprint was examined by Schmidt and Lloyd (1999) in a "tiger bush" landscape, composed of patches of vegetation (trees and shrubs) and bare sand, finding that changes in tower position can cause bias in terms of which parts of a landscape are measured by EC flux instrumentation. This finding underlines the importance of considering flux-footprints, and the research in the present study applies the flux-footprint concept to examine spatial heterogeneity in a much more complex landscape. As well, Chen et al. (2009) used footprint weighting to compare EC tower measurements with remote sensing based measurements of GPP at the same study site, and the present research builds on this approach, utilizing the flux-tower footprint to analyze finely detailed biometric measurements and characterize forest heterogeneity within an EC flux-tower footprint.  Previous work by Humphreys et al. (2006) described the stands at the coastal BC flux station by the range of biometric values observed in ground plots. With the methodology in the present study, we are able to estimate site-level spatially distributed and footprint flux probability density weighted means and standard deviations of biometric values observed in ground plots. In future research, we will directly compare changes in biometric measurements over a time interval with long-term tower NEP measurements over the same time interval using this methodology, providing an independent validation of EC tower based  58  measurements, and insight into the detailed ecological processes driving canopy scale C fluxes.  This type of spatial scaling is applicable globally to flux-towers at sites with nonuniform distributions of vegetation. By validating EC flux-tower measurements of C fluxes, and by quantifying the C stock change in the biologic components contributing to stand level fluxes, a greater understanding of forest C dynamics can be achieved. This understanding will improve process models, which are used to inform management and policy decisions.  59  3.5 References Amiro, B., Todd, J., Wotton, B., Logan, K., Flannigan, M., Stocks, B., Mason, J., Martell, D., and Hirsch, K. 2001. Direct carbon emissions from Canadian forest fires, 1959-1999. Canadian Journal of Forest Research. 31, 512-525. Baldocchi, D. 2008. Breathing of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems. Australian Journal of Botany. 56, 1-26. Banner, A., Meidinger, D.V., Lea, E.C., Maxwell, R.E., and Sacken, B.C. 1996. Ecosystem mapping methods for British Columbia. Environmental Monitoring and Assessment. 39, 97-117. Birdsey, R. A. 1992. Carbon storage and accumulation in United States forest ecosystems. U.S. Dept. of Agriculture, Forest Service, Washington, D.C. Brown, S. 2002. Measuring carbon in forests: current status and future challenges. Environmental Pollution116:363-372. Carmean, W. H. 1975. Forest site quality evaluation in the United States. Advances in Agronomy. 27, 209-269. Clark, D. A., Brown, S., Kicklighter, D.W., Chambers, J.Q., Thomlinson, J.R., and Ni, J. 2001. Measuring net primary production in forests: Concepts and field methods. Ecological Applications. 11, 356-370. Coops, N., Hilker, T., Wulder, M., St-Onge, B., Newnham, G., Siggins, A., and Trofymow, J.A. 2007. Estimating canopy structure of Douglas-fir forest stands from discrete-return lidar. Trees - Structure and Function. 21:295310. Chen, B., Black, T.A., Coops, N.C., Hilker, T., Trofymow, J.A. (Tony), and Morgenstern, K. 2009. Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Boundary-Layer Meteorology. 130:137-167. Crookston, N., and Finley, A. 2008. yaImpute: An R package for kNN imputation. Journal of Statistical Software. 23. Dixon, R. K., Solomon, A.M., Brown, S., Houghton, R.A., Trexier, M.C., and Wisniewski, J. 1994. Carbon pools and flux of global forest ecosystems. Science. 263, 185-190.  60  Ferster, C.J., Coops, N.C., Trofymow, J.A. (Tony). 2009. Above ground large tree mass estimation in a coastal forest in British Columbia using plot level metrics and individual tree detection. Canadian Journal of Remote Sensing. 35:270 – 275. Finnigan, J. 2004. The footprint concept in complex terrain. Agricultural and Forest Meteorology. 127, 117-129. Franklin, J. 1995. Predictive vegetation mapping: geographic modelling of biospatial patterns in relation to environmental gradients. Progress in Physical Geography. 19, 474-499. Gillis, M., Omule, A., and Brierley, T. 2005. Monitoring Canada's forests: The National Forest Inventory. Forestry Chronicle. 81, 214-221. . Göckede, M., Rebmann, C., and Foken, T. 2004. A combination of quality assessment tools for eddy covariance measurements with footprint modelling for the characterisation of complex sites. Agricultural and Forest Meteorology. 127, 175-188. Gougeon, F. 1995. A crown-following approach to the automatic delineation of individual tree crowns in high spatial resolution aerial images. Canadian Journal of Remote Sensing 21:274-284. Granger, J. E., and Schulze, R.E. 1977. Incoming solar radiation patterns and vegetation response: examples from the Natal Drakensberg. Plant Ecology. 35, 47-54. Houghton, R. A. 1999. The annual net flux of carbon to the atmosphere from changes in land use. Tellus B. 51, 298-313. Houghton, R. A. 2003. Revised estimates of the annual net flux of carbon to the atmosphere from changes in land use and land management. Tellus B. 55, 378-390. Houghton, R. A. 2005. Aboveground forest biomass and the global carbon balance. Global Change Biology. 11, 945-958. Houghton, R. A., Hackler, J.L., and Lawrence, K.T. 1999. The US carbon budget contributions from land-use change. Science. 285, 574. Hudak, A. T., Crookston, N.L., Evans, J.S., Hall, D.E., and Falkowski, M.J. 2008. “Nearest neighbor imputation of species-level, plot-scale forest structure attributes from lidar data.” Remote Sensing of Environment 112, 22322245.  61  Humphreys, E. R., Black, T.A., Morgenstern, K., Cai, T., Drewitt, G.B., Nesic, Z., and Trofymow, J.A. 2006. Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of development after clearcut harvesting. Agricultural and Forest Meteorology. 140, 6-22. IPCC. 2007. Fourth Assessment Report: Working Group I "The Physical Science Basis". Retrieved April 21, 2008, from http://www.ipcc.ch/ipccreports/ar4wg1.htm. Jungen, J. R. 1985. Soils of southern Vancouver Island. British Columbia Soil Survey #44, MOE Technical Report #17, Ministry of Environment, Victoria, BC. Kurz, W., Apps, M., Banfield, E., and Stinson, G. 2002. Forest carbon accounting at the operational scale. Forestry Chronicle. 78, 672-679. Kurz, W., and Apps, M. 2006. Developing Canada's national forest carbon monitoring, accounting and reporting system to meet the reporting requirements of the Kyoto Protocol. Mitigation and Adaptation Strategies for Global Change. 11, 33-43. Kurz, W. A., and Apps, M.J. 1999. A 70-year retrospective analysis of carbon fluxes in the Canadian forest sector. Ecological Applications. 9, 526-547. Kurz, W. A., Dymond, C.C., Stinson, G., Rampley, G.J., Neilson, E.T., Carroll, A.L., Ebata, T., and Safranyik, L. 2008. Mountain pine beetle and forest carbon feedback to climate change. Nature. 452, 987-990. Landsberg, J., and Waring, R.. 1997. A generalized model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. Forest Ecology and Management. 95, 209-228. Lefsky, M. A., Hudak, A.T., Cohen, W.B., and Acker, S. 2005. Geographic variability in lidar predictions of forest stand structure in the Pacific Northwest. Remote Sensing of Environment. 95, 532-548. Litton, C. M., Raich, J.W., and Ryan, M.G. 2007. Carbon allocation in forest ecosystems. Global Change Biology 13:2089-2109. MacKinnon, A. 2003. West coast, temperate, old-growth forests. Forestry Chronicle. 79, 475-484. Mahalanobis, P.C. 1936. On the generalized distance in statistics. Proceedings of the National Institute of Science of India. 12. pp. 49-55.  62  McGarigal, K., Tagil, S., and Cushman, S.A. 2009. Surface metrics: An alternative to patch metrics for the quantification of landscape structure. Landscape Ecology. 24:433-450. Moeur, M., and Stage, A.R. 1995. Most similar neighbor: An improved sampling inference procedure for natural resource planning. Forest Science. 41, 337359. Næsset, E. 2004. Practical large-scale forest stand inventory using a small-footprint airborne scanning laser. Scandinavian Journal of Forest Research. 19, 164179. Næsset, E., and Gobakken, T. 2008. Estimation of above- and below-ground biomass across regions of the boreal forest zone using airborne laser. Remote Sensing of Environment. 112, 3079-3090. NFI. 2004a. Canada's National Forest Inventory - Ground Sampling Guidlines. Retrieved April 10, 2008, from https://nfi.nfis.org/documentation/ground_plot/Gp_guidelines_v4.1.pdf. NFI. 2004b. Canada's National Forest Inventory: National Standards for Ground Plots Compilation Procedures. Retrieved October 3, 2008, from https://nfi.nfis.org/documentation/ground_plot/Gp_compilation_v1.5.pdf. Nilsson, S., and Schopfhauser, W. 1995. The carbon-sequestration potential of a global afforestation program. Climatic Change. 30, 267-293. Palmer, C. J., Smith, W.D., and Conkling, B.L. 2002. Development of a protocol for monitoring status and trends in forest soil carbon at a national level. Environmental Pollution. 116, 209-219. Panel, (Carmanah Sound Science Panel) 1995. Sustainable Ecosystem Management in Clayoquot Sound: Planning and Practices. Gov. of British Columbia, Victoria, BC. Pojar, J., Klinka, K., and Demarchi, D. 1991. Chapter 6: Coastal Western Hemlock Zone. In: Meidinger, D., Pojar, J. (Eds.). Ecosystems of British Columbia. BC Special Report Series. Pages 95-111 . Research Branch, Ministry of Forests. Roberts, D. W., and Cooper, S.V. 1989. Concepts and techniques of vegetation mapping, Land classifications based on vegetation: Applications for resource management. GTR-INT-257. USDA Forest Service Intermountain Research Station, Ogden, UT. pp. 90-96.  63  Schmid, H. P. 2002. Footprint modeling for vegetation atmosphere exchange studies: a review and perspective. Agricultural and Forest Meteorology. 113, 159-183. Schmid, H. P., and Lloyd, C.R. 1999. Spatial representativeness and the location bias of flux footprints over inhomogeneous areas. Agricultural and Forest Meteorology. 93, 195-209. Schreuder, H. T., Ernst, R., and Ramirez-Maldonado, H. 2004. Statistical Techniques for Sampling and Monitoring Natural Resources. U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. RMRS-GTR-126:111. Stage, A. R., and Crookston, N.L. 2007. Partitioning error components for accuracy-assessment of near-neighbor methods of imputation. Forest Science. 53, 62-72. Trofymow, J. A., Stinson, G., and Kurz, W.A. 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. Forest Ecology and Management. 256, 1677-1691. Wulder, M. A., Bater, C.W., Coops, N.C., Hilker, T., and White, J.C. 2008. The role of lidar in sustainable forest management. The Forestry Chronicle 84, 803-826.  64  4. Comparing Biometric Carbon Stock Changes with Eddy-Covariance Flux-Tower Based Cumulative NEP Measurements at Three Sites in a Coastal Forest in British Columbia3 4.1 Introduction Forests are a large component of the global carbon (C), containing an estimated 1146 PgC (Dixon et al., 1994). It is therefore critical that we understand forest C dynamics, including forest carbon stock components and transfer mechanisms in order to develop accurate forest C models such as the Carbon Budget Model of the Canadian Forest Sector (CBM-CFS3) (Kurz et al., 2009), 3PG (Landsberg and Waring, 1997), and ecosys (Grant et al., 2007), amongst others, to inform forest management policy and for national and international reporting (Kurz and Apps 2006). Factors affecting forest C dynamics include natural (e.g. fire, insect outbreaks) and anthropogenic disturbances such as land use management and change (eg. harvest, reforestation, deforestation) (Kurz et al. 2002) as well as weather (Morgenstern et al. 2006), fertilization (Jassal et al. 2009) and species composition, which is related to site stand age, disturbance history and site edaphic characteristics (Humphreys et al., 2006, Krishnan et al., 2009).  4.1.1 EC Flux-Towers To measure the net exchange of C between land ecosystems and the atmosphere (NEE with –NEE referred to as net ecosystem production, NEP), a global network of over 400 eddy-covariance (EC) flux stations has been established across a range 3  A version of this chapter will be submitted for publication. Ferster, C.J., Trofymow, J.A.(Tony), Coops, N.C., Chen, B., Black, T.A. Comparing biometric carbon stock changes with eddycovariance flux-tower based cumulative NEP measurements at three sites in a coastal forest in British Columbia.  65  of ecosystems, building an extensive data record of NEP. The majority of these towers are located on sites not undergoing major disturbances and so the fluxes measured reflect the interaction of weather, vegetation composition, and seasonal phenology. EC flux-towers use micrometeorological equipment to make continuous half-hourly measurements at the canopy scale of CO2, water vapour, sensible heat, and some are equipped to measure the exchange of other trace gases (Baldocchi 2008). Some EC flux-tower sites have records spanning up to almost two decades. The source area contributing to measurements made by the instruments mounted on the EC flux-tower, the flux-tower footprint, is variable in size and shape depending on wind speed and direction, surface roughness, and instrument height (Finnigan, 2004; Schmid, 2002). Early estimates modelled fluxfootprints as simple ovals. More recent estimates of flux-footprints may demonstrate more complex geometries and continuous probability density surfaces (Chen et al. 2009).  4.1.2 Biometric Measurements of Forest C Stocks Biometric measurements can be used to quantify forest C stocks (Dixon et al., 1994; Clark et al., 2004; NFI, 2009), and by measuring forest carbon stocks at two times with a sufficient interval between measurements, the C stock change can be determined and used to estimate the net exchange between the atmosphere and the forest (Clark et al. 2001; Smith, 2004). Measurements are typically made of individual plants, for example, all trees within a plot are measured for diameter, species, and height, then allometric relationships are used to determine tree mass.  66  Similarly, woody debris pieces are measured along transects using line intersect sampling and allometric relationships are applied to find biomass. Other components are made by direct measurements. For example, understory vegetation can be sampled, dried, and weighed, and measurements of individual components summed for plot-level totals. By placing sample plots across a landscape using an established sampling scheme, such as stratified random, landscape level inferences of ecosystem C stocks can be made (Schreuder et al., 2004).  4.1.3 Comparing EC Flux-Tower and Biometric Measurements Comparing EC flux-tower cumulative NEP (ΣNEP) and biometric measurements of C stock changes (∆C stocks) over the same period is useful for a number of reasons. First, as an independent validation, biometric measurements should confirm NEP measurements made at EC flux-towers. Second, comparisons with biometric data provide more specific information about the ecological changes that are driving stand canopy-level carbon exchange. Third, broad biometric measurement datasets are available with a greater spatial extent than the global EC flux-tower network. For example, the Canadian National Forest Inventory (NFI) samples approximately 22,000 locations across Canada (Gillis et al, 2005). Data from EC flux-towers used in combination with biometric data will aid the proper interpretation of biometric datasets, in particular for components that are difficult to measure, such as woody debris, snag, stump, and stump root decomposition following harvest (Janisch et al., 2005). However, comparing biometric and EC flux-tower measurements requires comparisons to be made across spatial and  67  temporal scales. Therefore, these comparisons pose methodological challenges beyond individual biometric and EC flux-tower measurements.  Comparisons between EC tower measurements of cumulative NEP and biometric measurements of C stock changes have been completed at a limited number of sites globally (Gough et al., 2008; Kominami et al., 2007; Miller et al., 2004; Curtis et al., 2002; Ehman et al., 2002; Barford et al., 2001; Grainier et al., 2000; Schulze et al., 2000), and at two sites of the Canadian Carbon Program (CCP), Boreal  Ecosystems Research and Monitoring Sites (BERMS) station, located in the boreal forest of northern Saskatchewan (Theede, 2007). Over a ten-year interval (1994 – 2004) Theede (2007), found reasonable agreement of 15.6 ± 4.0 (biometric) and 18.2 ± 8.09 MgC/ha (tower NEP) at the Old Aspen site and 5.8 ± 2.0 MgC/ha (biometric) and 6.9 ± 1.6 MgC/ha (tower NEP) at the Old Jack Pine site. In addition, this study provided detailed information about the biological components contributing to canopy-level measurements of flux. The homogeneous stand structure and large fetch in the boreal forest reduced the need for consideration of spatial vegetation structure and footprint distributions.  These previous studies have provided valuable information about forest C dynamics, however, at sites with more complex vegetation structure, such as coastal BC forests, the spatial distribution of vegetation structure and the footprint distribution are important considerations for accurate comparisons of EC fluxtower and biometric measurements (Schmidt and Lloyd, 1999). More complex ecosystems, therefore, require more complex models and methods developed for  68  accurate comparisons of EC flux-tower measurements and biometric measurements. A limited number of studies, including Ehman et al. (2002) and Ferster et al. (2009) applied EC flux-tower footprint weighting to biometric measurements, where areas of the footprint area that contribute more to towermeasured flux are weighted more heavily than other areas.  In this paper, we will compute the difference between 2002 and 2006 biometric measurements of ecosystem C stocks in ground plots at the CCP British Columbia Flux Station. This flux station measures C flux in Canada’s most productive forests, and therefore, understanding of this ecosystem's role in the C cycle depends on proper interpretation of the measurements at this station. Proper interpretation of fluxes, however, presents unique challenges since coastal forests are typically affected by small area disturbance such as windthrow and mass movement (landslides) (MacKinnon, 2003). Through a combination of variable topography and a complex disturbance history (Trofymow et al. 2008), the EC flux-tower sites possess fine-scale spatial variation in forest structure, thus presenting a challenge to the interpretation of EC flux data (Göckede et al., 2004; Schmid and Lloyd, 1999). To allow the comparison between biometric and EC flux-tower measurements, our approach consisted of four steps: first, we found the difference in biometric C stocks from ground plot measurements made in 2002 and 2006 (four year period); second, we modelled C stocks changes for components that are difficult to measure directly; third, we described stand structure at a fine spatial resolution, stratified the site based on forest structure, and applied footprint probability density weighting; and fourth, biometric ∆C stocks as unweighted and weighted means (with  69  measured and modelled values) were calculated and then were compared with EC flux-tower-measured ΣNEP over the same period.  4.2 Methods 4.2.1 Study Area British Columbia Flux Station sites are located on the leeward central east coast of Vancouver Island (Figure 4.1), in the Coastal Western Hemlock biogeoglimatic zone (Green and Klinka, 1994). The mean annual air temperature is 10˚ Celsius. Douglas-fir (Pseudotsuga menziesii var menziesii (Mirb.) Franco) is the dominant species, with lesser amounts of western hemlock (Tsuga heterophylla (Raf.) Sarg.), western redcedar (Thuja plicata Donn ex D. Donn), and red alder (Alnus rubra Bong.). Three towers were placed in a chronosequence of stands: a young forest (near end of rotation) established in 1949 (DF1949), a pole sapling stand established in 1988 (HDF1988), and a recent clearcut established in 2000 (HDF2000).  70  Figure 4.1. The Canadian Carbon Program’s British Columbia Flux Station consists of three sites on the central east coast of Vancouver Island.  4.2.2 DF1949 The Oyster River young forest (near-end-of-rotation) stand at DF1949 was harvested of its original old growth timber in 1937, 1938, and 1943. Following harvest, broadcast slash burning was applied in 1938 and 1939. In some places slash burning escaped into standing timber, which was subsequently harvested and broadcast burned again. Due to limited regeneration, the site was planted in 1949, leading to the establishment of the current stand (Trofymow et al., 2008). Elevation ranges from 260 m to 470 m above sea level, with gentler slopes and generally wetter soils in the eastern half of the footprint. The soil is Quimper gravelly sandy loam overlying a dense compact till at the 1 m depth (Jungen, 1985). As a result of the relatively complex disturbance history and topography, the stand is diverse in terms of stand structure with a much finer scale variation than is denoted in operational forest inventory mapping. For a detailed description of site  71  characteristics, see Trofymow et al. (2008). To cover the variation in stand conditions within the tower footprint area at this site (162 ha), a total of 12 NFIstyle groundplots were established and measured in 2002 and 2006 (NFI plot numbers 2000111, 2000112, 2000113, 2000121, 2000122, 2000123, 2000141, 2000142, 2000144, 2000181, and 2000182).  4.2.3 HDF1988 The Buckley Bay pole-sapling stand at HDF1988 was harvested from secondgrowth timber in 1987. The second-growth stand was composed of Douglas-fir, western redcedar, and western hemlock in the upland parts of the site, and Douglasfir, western redcedar and red alder in the downslope areas. Following harvest, logging slash was broadcast burned, and the site was replanted. Elevations at the site range between 150 m and 220 m, gently sloping down toward the east. The soils are deep (> 1 m depth) Reegan gravelly loam (Jungen 1985). Since the fluxtower was shorter than at DF1949, the footprint area at this site was smaller (35 ha) and a total of 6 NFI style groundplots were established and measured in 2002 and 2006 (NFI plot numbers 2000611, 2000613, 2000614, 2000621, 2000622, and 2000624).  4.2.4 HDF2000 The Oyster River clearcut stand at HDF2000 was harvested from a second-growth stand in 2000. Logging slash was piled and burned, and the site was replanted in 2000. Elevations at the site range from 160 m to 190 m, and the site is on generally flat upper and lower terraced benches with the east end of the site wetter than the west end. The soils are Hawarth very-gravelly loamy sand. (Jungen, 1985). The  72  short flux-tower meant this site had the smallest fetch (14.32 ha) and a total of 9 NFI style groundplots were established and measured in 2002 and 2006 (NFI plot numbers 2000311, 2000313, 2000314, 2000321, 2000323, 2000324, 2000331, 2000333, and 2000334).  4.2.5 Site Stratification To account for spatial heterogeneity in stand structure, the three sites were stratified following Ferster et al. (2009). This methodology uses environmental predictor variables from GIS and remote sensing data (Table 4.1) to extrapolate measurements from sample plots across the forest site. A nearest-neighbour classification was used (Crookston and Finlay, 2009), utilizing the Mahalanobis distance (Mahalanobis, 1937), a dimensionless statistical transform of the covariance matrix to identify which ground plot was most similar to each footprint cell, based on environmental predictor variables. Then, for each footprint cell, the detailed biometric target measurements from the most representative ground plot are assumed to be representative, and these assigned target variables at each footprint cell are used for calculation of site-level means. As an estimate of model error for the target biometric variables, the root mean square difference (RMSD) between a plot and its nearest neighbour is commonly used. However, this is not analogous to the standard error of a parametric estimate since it contains other sources of error beyond the error of the mean (Stage and Crookston, 2007). Regression estimators are designed to explicitly minimize this error, while nearest neighbour classification schemes are designed for spatial representativeness and  73  maintaining distribution shape. Furthermore, stratified sampling strategies strive to reduce sampling redundancy among sample plots (Schreuder et al., 2004) and thus sites with well designed and unique sample plots will have a very large RMSD, despite the reliable estimates.  In this study, correlations between all measured C stocks components (overstory, understory, detritus, and soil) both in 2002 and 2006 and the environmental predictor variables are used to select meaningful environmental predictors. Environmental predictor variables with correlations with biometric C stocks in either year that were significant at the 95 % confidence level were selected for use in the nearest neighbour classification and stratification of the site. Following stratification, the 2002 – 2006 cumulative NEP footprint probability density surface for each site (Chen et al. 2008) was used to weight each 10 m by 10 m footprint cell.  4.2.6 Biometric C Stocks and Changes Ground sample plots were established and measured by the Canadian Forest Service at the three sites following Canadian National Forest Inventory (NFI) guidelines (NFI 2009). Initial measurements were made in September 2002 and remeasurements were made in September 2006. Measurements were submitted to the NFI database for data compilation, which includes the calculation of plot-level values from individual measurements and application of allometric equations to overstory tree and woody debris measurements (NFI 2009). Estimates of biomass were converted to C composition by applying an approximate conversion  74  coefficient of 0.5 (Mathews, 1993). Estimates of aboveground tree C changes were made using national allometric equations from the NFI. Biometric measurements of other ecosystem C stock components, such as woody debris, forest floor organic matter and mineral soil were also made on the two dates, however, given the spatial variability in these components and small relative changes in mass over 4 years due to decomposition, estimates of changes in these components will be more uncertain (Curtis et al. 2002). Furthermore, live overstory coarse roots were not measured, so their C mass and changes due to growth and decomposition were estimated by modelling (see below).  4.2.7 Overstory Trees Overstory trees were measured in an 11.28 m radius circular plot. Where trees were very numerous (ie. HDF1988), a half plot or 5 m radius circular plot was measured due to logistical constraints and to reduce redundancy. Allometric equations from the NFI were used to estimate stem, bark, branch, and foliar biomass. These biomass values were converted to C values by applying the coefficients of 0.52, 0.56, 0.52, and 0.52, respectively (Mathews, 1993). When finding overstory C stock differences, corrections were not made for ingrowth or mortality (e.g., Clark et al. 2001), since these components were measured elsewhere in the plot (e.g.,  seedlings in understory, fallen trees in coarse woody debris). Coarse and fine root biomass C for live overstory trees was estimated using equations by Li et al. (2003).  75  4.2.8 Understory Vegetation Understory vegetation from four 1 m x 1 m microplots was destructively sampled and sent to the laboratory for drying and weighing for determination of mass which was converted to C content assuming 50% C.  4.2.9 Woody Debris Woody debris was measured for diameter, species, and decay class along four 15 m transects at each plot. NFI compilation routines utilize algorithms for line intercept sampling to determine plot level volume and a lookup table of woody debris density values by species and decay class to determine mass. C content was determined assuming 50% C content of woody debris mass. Fine woody debris (2 mm to <1 cm diameter) was sampled from within the four 1 × 1 m understory microplots in each sample plot and sent for laboratory drying, weighing, and determination of mass, and C content assuming dry mass is composed of 50% C.  4.2.10 Forest Floor Surface substrates were measured along four 15 m long transects to determine the average depth and percent area surface coverage in the entire plot. The forest floor organic litter, fibric, and humus (LFH) layer was sampled and average depths measured from within four 20 × 20 cm templates (one at each transect). These destructive samples were collected, sieved, dried, and weighed to determine bulk density and subsamples sent for laboratory analysis of %C. There was a difference in laboratory processing, where >8mm woody pieces that were excluded in the bulk density calculation and %C determination in 2002 were included in 2006.  76  4.2.11 Mineral Soil In 2002, mineral soil was sampled for bulk density and C concentrations in three depth intervals (0-15 cm, 15-35 cm, and 35-55 cm) from one 55 cm deep soil pit per ground plot and scaled to the entire groundplot discounting for the area without mineral soil (i.e. exposed bedrock). In 2002, samples were taken from a complete excavation. In 2006, only the 0-15 cm layer was re-measured for C content, but at four locations in each plot. In 2006, a 1 cm diameter soil corer was used to collect samples. Samples were stored in plastic bags at 2˚ Celsius prior to laboratory processing. For a detailed description of soil sampling methods see NFI (2009). For total C stock changes, only the 0-15 cm layer was considered.  4.2.12 Litterfall From 2002 – 2006, overstory fine litterfall (needles, leaves, cones, twigs) were collected quarterly in 3 0.189 m2 conical-mesh litterfall traps located in each ground plot at DF1949 and HDF1988. Annual litter fall masses were calculated and converted to C assuming 50% C.  4.2.13 Modeled Changes in Fine Woody Debris, Forest Floor, and Mineral Soil C In addition to the measured changes in fine woody debris, forest floor, and mineral soil (0 – 15 cm), the changes in these components was modeled using algorithm components and parameters from CBM-CFS3, a forest C accounting model used for national reporting of annual C inventories for Canada's managed forests (Kurz et al. 2009). This model uses the C mass of overstory trees and overstory tree roots  to estimate litterfall and root mortality transfers to the aboveground and  77  belowground detritus C pools respectively. For this study, the full model was not applied, but specific published algorithms and parameters were used to estimate LFH, fine woody debris, and mineral soil (0 – 15 cm) C stock changes. Where available, annual aboveground litterfall data collected by the Canadian Forest Service over 2002 - 2006 at the DF1949 and HDF1988 sites was used directly as an input into the detritus pools. Where litterfall data was not available (at HDF2000 due to the low height of trees), the default CBM-CFS3 relationships with tree biomass was used for litterfall inputs to the detritus pools. CBM-CFS3 parameters were used to model mortality and turnover in fine and coarse roots to estimate transfers to the detritus pools. CBM-CFS3 calculates decomposition rates of the different detritus pools as a function of mean annual temperature which determines the amount of C released to the atmosphere or transferred to a very slowly decaying soil C pool. The sub-model used in this study was run on an annual time step to calculate decomposition and transfers between the detrital and soil C pools.  4.2.14 Comparing Biometric Estimates of ∆C Stocks and EC Tower ΣNEP Calculations of tower NEP available in Krishnan et al. (2009) and Black et al. (2008) were used for this study. These values are gap filled for high stability conditions and corrected for lack of energy-balance closure. Annual NEP values for 2003, 2004, 2005, and 2006 were summed for each site to estimate the ΣNEP total for the four-year period.  Biometric-based site estimates of mean ∆C stocks were calculated as unweighted ground-plots and weighted by site classification and flux footprint using measured  78  and modeled component values. For unweighted means, the 95% confidence interval was used as an estimate of measurement error. For weighted by site classification and footprint means, the RMSD was used as an estimate of error.  4.3 Results 4.3.1 Site Stratification Selected environmental predictor variables are presented in Table 4.1, and maps of the sites with footprints and stratification units from most similar neighbour classification are presented in Figures 4.2 and 4.3 (a-c), respectively. The individual correlations used to select environmental predictor variables for the site classification are described as follows: Slope was an important predictor for all sites as it was strongly correlated with soil C stocks (flatter plots having higher soil C stocks), and at HDF2000, slope was also correlated with living tree C stocks. Due to the broad spatial extent at DF1949, forest inventory mapping captured some variation between plots, and at this site it was a significant predictor for overstory C stocks. At the other sites, forest inventory mapping did not indicate any variation within the stand, and at all sites the other predictor variables showed finer scale variation than shown in the forest inventory mapping. At HDF2000, NDVI was an important predictor for woody debris, forest floor, and mineral soil C stocks. Another remote sensing variable, individual tree crown density (Gougeon, 1995), was an important predictor of overstory C stocks at HDF1988, and understory C stocks at DF1949. Finally, at DF1949 and HDF1988, site series mapping was an important predictor for overstory and live C stock components.  79  Mahalanobis distance is mapped in Figure 4.2 (a-c). Lower percentiles represent relatively small distances which are well sampled and represented by the biometric ground plots. Areas with Mahalanobis distances less than the median value at each site are shown in green. Areas with high percentiles had relatively large distances, and should be considered less well represented by the biometric ground-plots. Areas with Mahalanobis distances greater than the median at each site are shown in red. Areas with the smallest Mahalanobis distances and best ground-sample-plot representation are close to the towers in the areas contributing most to NEP measurements (within the 50% isoline). Notable areas with large Mahalanobis distances include areas with different forest cover types (e.g., patches of hardwoods at DF1949), very few overstory trees (at HDF1988), and short steep slopes (between bench terraces at HDF2000).  80  Table 4.1. Environmental predictor variables used to determine stratification units at DF1949, HDF1988, and HDF2000. Selection Coverage Forest Inventory Disturbance History (Trofymow et al. 2008)  Topography3  Variable  Units  DF1949  Site Index  m  yes  Top Height 1st Harvest  m year  1  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  aspect  azimuth  Elevation  m asl  Slope  degrees  SCOSA SSINA  Forest Inventory (Trofymow et al.2008) Soil Survey of Canada (Jungen, 1985) CFS and Forest Companies (Trofymow et al.2008)  yes  yes  yes  yes  yes  yes  yes  yes  yes  yes  yes  stems/ha  yes  yes  -  NDVI  yes  2  Dominant Canopy Tree Density NDVI  yes  HDF2000  yes  2  TSRAI2 1999 Orthophoto (Gougeon, 1995) 2004 Multispectral4  yes  HDF1988  yes  Cover Species Site Species Most Common Soil Association Site Series  yes  yes  1  Site Index: Tree height at 50 years age. Topographic Variables: SCOSA = Slope * cos (Aspect), SSINA = Slope * sin (Aspect); TSRAI (Topographic Solar Radiation Aspect Index) = (1cos((pi/180)(Aspect-30)))/2 (Roberts and Cooper, 1989) 3 Topography from 2004 lidar survey at DF1949 and HDF2000 (Coops et al. 2008). Topography from 1:50,000 National Topographic Series map at HDF1988. 4 Multispectral data from 2004 Quickbird survey at DF1949 and HDF2000. Multispectral data from 2004 Landsat scene at HDF1988. 2  81  1 2 3  Figure 4.2 Site flux footprints for A) DF1949, B) HDF1988, C) HDF2000. Background imagery is pan-sharpened 2004 Quickbird for A and C, and true color orthophoto for B. 82  1 2  Figure 4.3 Most similar neighbour (MSN) classification of stratification units for A) DF1949, B) HDF1988, C) HDF2000.  83  1 2 3 4 5  Figure 4.4 Total Mahalanobis distance for MSN stratification units for A) DF1949, B) HDF1988, C) HDF2000. 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 biometric sample plots, while footprint cells with red tones are above the median Mahalanobis distance for each site and are less well represented. 84  Spatial weighting resulted in between 650% and 1% difference for individual plot weighting towards the site classification and footprint weighted site mean compared to even weighting for the unweighted site mean. As a result, for C stock mean differences, weighted standard deviations were generally lower than unweighted standard deviations.  4.3.2 Biometric C Stock Changes Biometric C stock changes (∆C) from are presented in Table 4.2 as the absolute difference between ground plot measurements. Modeled estimates of ∆C in fine woody debris (FWD), forest floor (LFH), and mineral soil are presented in Table 4.3. Individual components are described in detail below.  85  Table 4.2 Site mean and sample standard deviations (SD) for measured or modelled (denoted by *) biometric ∆C stocks in groundplots from 2002 to 2006. Error is represented by the 95% confidence interval (CI) for unweighted means, and root mean square difference (RMSD) for weighted means for each respective site. All units are in MgC/ha/4 years. Unweighted Component Stem Bark Branch Foliage Live Tree Understory Aboveground Live Coarse Woody Debris Standing Dead Stumps LFH Fine Woody Debris Total Detritus Total Aboveground Coarse Roots* Fine Roots* Total Roots Mineral Soil (0-15 cm) Total Belowground Total Difference  DF1949 Mean ± CI (SD) 5.71 ± 2.99 (4.70) 0.93 ± 0.58 (0.91) 0.67 ± 0.76 (1.19) -0.11 ± 0.44 (0.69) 7.2 ± 4.64 (7.31) 1.33 ± 1.13 (1.78) 8.53 ± 4.54 (7.14) 1.56 ± 7.85 (12.35) -3.42 ± 7.08 (11.15) 0.00 ± 0.00 (0.00) 18.13 ± 12.25 (19.27) -0.92 ± 0.65 (1.02) 15.35 ± 18.07 (28.43) 23.88 ± 18.97 (29.86) 1.56 ± 1.00 (1.58) 0.04 ± 0.03 (0.05) 1.61 ± 1.03 (1.62) 19.13 ± 15.69 (24.7) 20.73 ± 15.89 (25.01) 44.61 ± 23.63 (37.19)  HDF1988 Mean ± CI (SD) 4.41 ± 1.31 (1.25) 1.01 ± 0.33 (0.31) 1.81 ± 0.54 (0.52) 1.49 ± 0.46 (0.44) 8.72 ± 2.47 (2.36) -2.58 ± 4.73 (4.51) 6.14 ± 6.40 (6.09) -1.49 ± 19.76 (18.83) 0.12 ± 0.23 (0.22) 0 ± 0.00 (0.00) 88.77 ± 72.23 (68.83) 0.8 ± 0.86 (0.82) 88.19 ± 72.36 (68.95) 94.33 ± 67.88 (64.69) 1.57 ± 0.41 (0.39) 0.69 ± 0.11 (0.11) 2.25 ± 0.51 (0.49) 10.07 ± 32.72 (31.18) 12.32 ± 32.7 (31.16) 106.65 ± 87.9 (83.75)  HDF2000 Mean ± CI (SD) 0.71 ± 0.71 (1.12) 0.15 ± 0.15 (0.23) 0.1 ± 0.12 (0.19) 0.38 ± 0.34 (0.54) 1.33 ± 1.27 (2.00) 1.74 ± 1.21 (1.91) 3.08 ± 2.02 (3.18) -1.42 ± 3.99 (6.29) -0.06 ± 0.1 (0.15) 0 ± 0.00 (0.00) 9.59 ± 17.04 (26.81) -4.62 ± 1.67 (2.63) 3.49 ± 17.84 (28.09) 6.57 ± 18.11 (28.5) 0.33 ± 0.31 (0.49) 0.22 ± 0.21 (0.33) 0.55 ± 0.52 (0.82) 21.84 ± 8.55 (13.46) 22.38 ± 8.49 (13.37) 28.95 ± 19.04 (29.97)  Weighted Component (RMSD for DF1949, HDF1988, and HDF2000, respectively) Stem (4.31, 1.63, 1.39) Bark (0.84, 0.41, 0.28) Branch (1.18, 0.84, 0.16) Foliage (0.78, 0.55, 0.65) Live Tree (6.64, 3.21, 2.47) Understory (1.9, 6.3, 2.56) Aboveground Live (6.38, 9.08, 4.27) Coarse Woody Debris (15.41, 16.98, 6.9) Standing Dead (10.37, 0.34, 0.18) Stumps (0, 0, 0) LFH (41.97, 105.09, 13.12) Fine Woody Debris (1.1, 1.05, 3.4) Total Detritus (45.09, 103.29, 15.36) Total Aboveground (47.11, 95.77, 16.02) Coarse Roots (1.43, 0.57, 0.57) Fine Roots (0.04, 0.13, 0.38) Total Roots (1.47, 0.68, 0.96) Mineral Soil (0-15 cm) (29.08, 53.74, 14.99) Total Belowground (1.66, 0.93, 1.12) Total Difference (41.42, 135.45, 13.06)  DF1949 Mean (SD) 4.81 (4.48) 0.86 (0.96) 0.71 (1.31) 0.03 (0.79) 6.41 (7.54) 2.19 (0.66) 8.60 (6.67) -3.54 (3.73) -4.54 (5.55) 0.00 (0.00) 20.67 (6.42) -2.37 (0.88) 10.22 (9.30) 18.82 (14.46) 1.39 (1.63) 0.03 (0.05) 1.42 (1.67) 8.65 (7.79) 10.07 (9.36) 28.89 (21.92)  HDF1988 Mean (SD) 4.62 (0.24) 1.06 (0.07) 1.84 (0.03) 1.65 (0.05) 9.17 (0.34) -2.95 (0.95) 6.22 (4.47) -8.59 (4.09) 0.06 (0.07) 0.00 (0.00) 75.07 (9.75) 0.92 (0.22) 67.46 (11.27) 73.69 (12.12) 1.59 (0.06) 0.73 (0.02) 2.32 (0.08) 18.12 (21.31) 20.44 (21.38) 94.13 (32.26)  HDF2000 Mean (SD) 0.37 (0.31) 0.08 (0.07) 0.06 (0.04) 0.22 (0.15) 0.73 (0.56) 2.14 (1.6) 2.87 (3.13) 0.33 (4.25) 0.00 (0.00) 0.00 (0.00) 14.09 (9.11) -6.37 (2.11) 8.05 (13.11) 10.91 (11.92) 0.25 (0.22) 0.18 (0.16) 0.44 (0.37) 18.18 (3.47) 18.62 (3.73) 29.53 (12.65)  86  Table 4.3 Means and sample standard deviations for modeled ∆C stocks over 2002 – 2006. Modeled components (denoted by *) are combined with NFI measured differences for total differences and summed totals (in bold). Model error is represented by the 95% confidence interval (CI) for unweighted means, and root mean square difference (RMSD) for the weighted means presented for each respective site. All units are in MgC/ha/4 years. Unweighted DF1949 Mean ± CI (SD)  HDF1988 Mean ± CI (SD)  HDF2000 Mean ± CI (SD)  7.64 ± 0.53 (0.84)  2.23 ± 0.3 (0.29)  0.29 ± 0.25 (0.4)  Total Aboveground Mineral Soil (0-15 cm) *  12.75 ± 7.57 (11.91) 3.11 ± 0.3 (0.47)  8.48 ± 6.47 (6.16) 0.83 ± 0.23 (0.22)  3.31 ± 2.13 (3.36) 0.1 ± 0.08 (0.13)  Total Belowground  4.72 ± 1.18 (1.86)  3.09 ± 0.67 (0.64)  0.64 ± 0.6 (0.95)  Total Difference  17.47 ± 8.31 (13.07)  11.57 ± 7.08 (6.75)  3.95 ± 2.6 (4.09)  Component (SEM) LFH and FWD *  Weighted Component (for DF1949, HDF1988, and HDF2000, repectively) LFH and FWD (0.96, 0.4, 0.49) *  DF1949 Mean (SD) 6.77 (0.27)  HDF1988 Mean (SD) 2.14 (0.06)  HDF2000 Mean (SD) 0.21 (0.18)  Total Aboveground (10.96, 9.22, 4.55) Mineral Soil (0-15 cm) (0.58, 0.31, 0.17) *  10.82 (6.45) 2.93 (0.14)  8.43 (1.11) 0.77 (0.03)  3.08 (1.46) 0.07 (0.06)  Total Belowground (1.66, 0.93, 1.12)  4.35 (1.80)  3.09 (0.08)  0.51 (0.43)  Total Difference (12.05, 10.08, 5.53)  15.18 (7.94)  11.52 (1.17)  3.59 (1.59)  4.3.3 Overstory Overstory trees were an important driver of site biometric ∆C stock changes. Livetree increment was greatest at HDF1988, slightly less at DF1949, and relatively small at HDF2000. Site classification and footprint weighting made a small difference in means at each site. Standing dead trees showed a small decrement in biometric ∆C stocks at DF1949. This was due to lower snag heights and stem numbers, from stem breakage during wind events exceeding the mortality rate, and resulting in transfers to the coarse and fine woody debris C pool. At HDF1988, standing dead C stocks increased slightly due to tree mortality for the pole sapling stand. Stump C changes were not measured, since it would be difficult to detect any change since decomposition would be minimal over four years (Table 4.2).  87  4.3.4 Understory Understory vegetation is a small component of the total ecosystem C stocks. A small increment was observed at DF1949 and HDF2000, while a small decrement was observed at HDF1988. Site classification and footprint weighting resulted in larger increments for DF1949 and HDF2000, while making little difference at HDF1988 (Table 4.2).  4.3.5 Coarse Woody debris Coarse woody debris C stocks decreased slightly over the time interval at HDF1988 and HDF2000, indicating that decomposition or transfers to other C pools took place over the measurement interval. This was a result of field crews classifying woody debris pieces into higher decay classes (more decayed) than in 2002. At DF1949, there was a small increase in coarse woody debris C stocks (Table 4.2).  4.3.6 Forest Floor Forest floor litter, fibric, and humus (LFH) increments measured in the ground plots were exceptionally large over the time interval. A large increase in was noticed in measured values of bulk density (BD) and attributed to mineral soil being included in the samples due to subjective differences in interpretation of the mineral/organic soil interface. Further, the inclusion of >8 mm woody material in 2006 BD calculation but not in 2002 BD inflated the forest floor C increment. Recalculation of the 2002 forest floor BD and %C with the >8mm woody material included made a small reduction in changes for forest floor values, however, C increments remained very high. Forest floor ∆C stocks were very spatially  88  variable, and therefore, footprint weighting had a large effect on the calculation of site-level means (Table 4.2).  4.3.7 Mineral Soil Mineral soil C changes in the 0-15cm layer were much higher than expected, in particular, at HDF1988. Mineral soil bulk densities in the 0-15cm layer were assumed not to change and thus the BD values measured in 2002 were also used in 2006. If BD decreased then the C change would be overestimated. At HDF00 and HDF1988, an increase in fine root mass from 2002 to 2006 included in the <2 mm soil fraction could also have accounted for some of the increase in soil C. Further, samples in 2002 were complete excavations, while samples in 2006 were collected using a soil corer, which may be biased to shallower depths in rocky soil. Mineral soil C stock changes had a high amount of spatial variability, and therefore, footprint weighting had a large effect on the calculation of site-level means (Table 4.2).  4.3.8 Modeled Changes in Fine Woody Debris, Forest Floor, and Mineral Soil C Modeled values of CBM-CFS3 using measured aboveground litterfall or estimated root mortality inputs were much lower than the measured values for changes in forest floor and mineral soil (0 – 15 cm). Modeled changes for LFH and fine woody debris were between 2% and 44% of the measured change. Modelled mineral soil changes were between 0.45% and 36% of measured changes. These large differences indicate that the measured values were likely overestimates, given the available inputs for transfers to this C stock pool. Measured fine woody debris changes were comparable with modeled values. In addition, modeled estimates had  89  less spatial variability, decreasing the effect of footprint classification and weighting for C stock changes between 2002 and 2006 (Table 4.2).  4.3.9 Comparing EC Tower ΣNEP and Biometric ∆C Stocks Between 2002 and 2006, EC flux-tower measurements of ΣNEP indicate that DF1949 was a sink (13.63 MgC/ha), HDF1988 was a small source (-1.93 MgC/ha), and HDF2000 was a large source (-20.08 MgC/ha) (Krishnan et al. 2009; Black et al. 2008; Humphreys et al. 2006). Morgenstern et al. (2004) reported that  uncertainty in the annual NEP measured using EC at DF1949 may be as much as 0.9 MgC/ha/year (due to systematic error in the EC measurement). This suggests that the uncertainty in the estimate at DF1949 may be ±3.6 MgC/ha over the fouryear measurement interval. The best agreement between biometric measurements of ecosystem ∆C stocks and tower ΣNEP between 2002 and 2006 was found at the young forest (near end-of-rotation) site (DF1949), the most mature stand at the station, using a site classification and footprint weighted mean of ground plots (Table 4.4), and the tower ΣNEP value also falls within the 95% confidence interval for the unweighted mean ∆C stocks for that site. For HDF1988 and HDF2000, the biometric ∆C stocks followed the same trend in magnitude, though negative biometric ∆C stock values were not found at these sites.  90  Table 4.4 Comparison of biometric ∆C stocks and tower ΣNEP 2002 - 2006. For unweighted biometric changes, the mean, 95% confidence interval, and standard error of the mean (SEM) are presented. For weighted biometric changes, the weighted mean and root mean squared difference (RMSD) are presented. All units are in MgC/ha/4 years. Estimate  DF1949  HDF1988  HDF2000  Tower Cumulative NEP  13.63 ± 3.5  -1.93  -20.08  44.61 ± 23.63 (10.74)  106.65 ± 83.75 (34.19)  28.95 ± 29.97 (8.65)  28.89 (41.42) 17.47 ± 8.31 (3.77)  94.13 (135.45) 11.57 ± 7.08 (2.75)  29.53 (13.06) 3.95 ± 2.60 (4.09)  15.18 (7.94)  11.52 (1.17)  3.59 (1.59)  Biometric Unweighted Change (SEM) Weighted Change (RMSD) Unweighted Change with Modeled Values (SEM) Weighted Change with Modeled Values (RMSD)  4.4 Discussion Footprint weighting by ecological attributes and footprint probability density significantly affected values of the site-level means of measured changes in forest floor and mineral soil C stocks likely due to high spatial variability. Modeled ∆C stock values for forest floor and mineral soil were more realistic, and these estimates had less spatial variability, thus reducing the effect of footprint weighting on site-level means. As well as for other ecosystem components, while the distribution of C stocks across the EC flux-tower footprint was highly variable (Ferster et al., 2009), the spatial distribution of biometric ∆C was less variable. Site classification and footprint weighted means and unweighted means of ground plots were similar and within the estimates of error for each method.  An important aspect of this study is the use of Mahalanobis distance to assess how representative the ground plots are of the site. Site sampling could be improved had ground sample plots been established where Mahalanobis distance are relatively high. Future studies wishing to compare biometric measurements with EC flux-  91  tower measurements may seek to derive accurate footprint estimates and use similar remote sensing and GIS environmental predictor variables prior to establishing sample plots and so ensure they are representative of site conditions. Another more costly alternative is to establish a much greater number of plots in a systematic basis around the tower, to ensure the spatial variation within the tower footprint is adequately captured.  Developing site-specific allometric is labour intensive, destructive, and costly, and therefore uncommon in most studies; however, use of existing relationships will lead to errors due to differences in tree architecture, and wood density. Further, the use of inappropriate allometric equations can be a significant source of error in forest productivity studies (Clark et al., 2002). The NFI data compilation procedure was a valuable source for regionally applicable allometric equations, however the error of this large pool of equations is not available may warrant further research. Losses due to decomposition of coarse woody debris, snags, stump and snag roots, are difficult to measure or model, but may represent a major source of respiration at the recently disturbed sites HDF1988 and HDF2000, since a relatively short period of time has elapsed since harvest. The estimate of understory ∆C stocks does not include losses due to herbivory, which are likely minimal. However, the observed decrease in aboveground understory C at HDF1988 may be a result of increasing canopy closure.  Slight decreases in woody debris C stocks, due to pieces of woody debris being classified into higher decay classes (more decomposed), are somewhat subjective  92  given that they depend on field crew judgement of decay class. However they may be a result of decomposition occurring since the last measurement. Fine woody debris showed a strong decrement at HDF2000, likely as leftover fine logging slash continued to decompose. A small increment was observed at HDF1988 likely as litterfall inputs from the shrubs and growing trees.  Measured forest floor ∆C stocks were much larger than expected. Further investigation is warranted to explain these differences in terms of sample collection, laboratory analysis, and data compilation. Therefore, we feel that at present the modeled estimate is more realistic. A potential limitation of the model, however, is that the model does not include contributions from understory vegetation. At HDF2000, understory vegetation forms a greater proportion of total C stocks, therefore, we expect that the true LFH increment is slightly larger than predicted by the model at HDF2000, though still much less than the measured difference. Another potential limitation of this approach is that litterfall samples were measured after a period of decomposition inside the collection traps (up to three months between measurements), and therefore litterfall inputs into the detritus pools may have been slightly larger, resulting in a slight underestimation of the LFH increment. Modeled changes in fine woody debris and forest floor C stocks were near to, or greater than, the overstory increment at all sites. Modeled mineral soil C stock changes were much lower than measured differences in the groundplots, and therefore systematic differences in sampling, processing, or data compilation may have resulted in higher C content in 2006, and therefore, further investigation is warranted. Another consideration is that spatial variability in soil C  93  stocks within the plot may be greater than the change over the time period. Modeled changes show a much smaller potential for increases than was measured in the samples.  Agreement between EC flux-tower-measured ΣNEP and biometric ∆C stocks was within the bounds of error estimates for the most mature stand, DF1949. At the other sites, EC flux-tower-measured NEP and biometric C stocks followed the same trend as tower measurements of NEP, but were of different magnitude and sign - at all sites, the change in biometric C stocks was positive, though HDF1988 and HDF2000 had negative EC flux-tower-measured NEP values. However, a major source of respiration was not included in the biometric measurements. Further investigation is needed to explain what this source of respiration might be. One possible explanation is that the biometric C stock inventories did not account for coarse root or woody debris decomposition, and following harvest, large C stocks are left to decompose as stumps, coarse roots, and harvest slash (Janisch et al., 2005). We hypothesize that the discrepancy between the biometric ∆C stocks  and tower ΣNEP at the recently harvested sites may be due to this continued decomposition.  4.5 Conclusion The global network of EC flux-towers is an important source of information, aiding our understanding of the global C cycle, and improving management decisions that mitigate the release of C to the atmosphere, and maximize C sequestration in the  94  forest. Comparing biometric ∆C stocks and tower ΣNEP at the same sites and time intervals are complementary sources of information that can increase our understanding and interpretation of the global EC flux and biometric data records. For example, measurement of biometric ∆C stocks helps to inform researchers of the individual component ∆C stock changes that determine canopy-level fluxes. As well, tower ΣNEP measurements may indicate that stump, coarse root, and harvest slash decomposition following harvest may be an important source of respiration that is not accounted for using biometric methods.  Biometric ∆C stocks in certain ecosystem components, such as forest floor (LFH) and mineral soil are challenging to measure over short intervals, since the spatial variation, potential for subjective differences in measurement, or measurement error may be larger than the change. In these cases, applying C budget models may provide a realistic estimate of the change. As well, decomposition from stumps, coarse roots, and harvest slash is likely an important driver of stand level C fluxes in stands following harvest. This is difficult to measure and difficult to model. Therefore, better estimates of decomposition and respiratory losses from these components are topics which warrant further research.  For sites with complex vegetation distributions, detailed footprint analysis is important. However, very few studies comparing EC flux measurements with biometric measurements consider the spatial distribution of vegetation within the footprint and the footprint probability density distribution over the stand being measured. In this study, classification of vegetation structure within the footprint  95  and footprint probability-density weighting was used to test the effects of spatial variability in stand structure and footprint distribution. This consideration allowed comparison of biometric ∆C stocks and tower ΣNEP at sites with complex topography and stand composition with greater confidence. Other studies comparing biometric ∆C stocks and tower ΣNEP may find remote sensing and GIS datasets useful to spatially characterize stand conditions, establish efficient and representative sampling designs for biometric plots, and used in combination with footprint probability distributions, increase our understanding of the effect of spatial variability on EC flux measurements.  96  4.6 References Baldocchi, D. 2008. Breathing of the terrestrial biosphere: lessons learned from a global network of carbon dioxide flux measurement systems. Australian Journal of Botany. 56:1-26. Barford, C. C., Wofsy, S.C., Goulden, M.L., Munger, J.W., Pyle, E.H., Urbanski, S.P., Hutyra, L., Saleska, S.R., Fitzjarrald, D., and Moore, K. 2001. Factors controlling long- and short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science. 294:1688-1691. Black, T.A., Jassal, R., and A. Fredeen, A. 2008. Carbon sequestration in british columbia's forests and management options. Report for the: Pacific Institute for Climate Solutions. Chen, B., Black, T.A., Coops, N.C., Hilker, T., Trofymow, J.A. (Tony), and Morgenstern, K. 2009. Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Boundary-Layer Meteorology. 130:137-167. Clark, D., Brown, S., Kicklighter, D., Chambers, J., Thomlinson, J., and Ni, J. 2001. Measuring net primary production in forests: Concepts and field methods. Ecological Applications. 11:356-370. Crookston, N. L., and Finley, A.O. 2008. yaImpute: an R package for k-NN imputation. Journal of Statistical Software. 23. Curtis, P., Hanson, P., Bolstad, P., Barford, C., Randolph, J., Schmid, H., and Wilson, K. 2002. Biometric and eddy-covariance based estimates of annual carbon storage in five eastern North American deciduous forests. Agricultural and Forest Meteorology. 113:3-19. Dixon, R. K., Solomon, A.M., Brown, S., Houghton, R.A., Trexier, M.C., and Wisniewski, J. 1994. Carbon pools and flux of global forest ecosystems. Science. 263:185-190. Ehman, J. L., Schmid, H.P., Grimmond, C.S.B., Randolph, J.C., Hanson. P.J., Wayson, C., and Cropley, F.D. 2002. An initial intercomparison of micrometeorological and ecological inventory estimates of carbon exchange in a mid-latitude deciduous forest. Global Change Biology. 8:575-589. Ferster, C., J., Trofymow, J.A. (Tony), Coops, N.C., Chen, B., Black, T.A., and Gougeon, F.A. 2009. Determination of carbon stock distributions in the flux footprint of an eddy-covariance tower in a coastal forest in British Columbia. Ecological Applications. In Review.  97  Finnigan, J. 2004. The footprint concept in complex terrain. Agricultural and Forest Meteorology. 127:117-129. Gillis, M., Omule, A., and Brierley, T. 2005. Monitoring Canada's forests: The National Forest Inventory. Forestry Chronicle. 81:214-221. Gougeon, F.A. 1995. A crown-following approach to the automatic delineation of individual tree crowns in high spatial resolution aerial images. Canadian Journal of Remote Sensing. 21:274-284. Gough, C., Vogel, C., Schmid, H., Su, H., and Curtis, P. 2008. Multi-year convergence of biometric and meteorological estimates of forest carbon storage. Agricultural and Forest Meteorology. 148:158-170. Göckede, M., Rebmann, C., and Foken, T. 2004. A combination of quality assessment tools for eddy covariance measurements with footprint modelling for the characterisation of complex sites. Agricultural and Forest Meteorology. 127:175-188. Granier, A., Ceschia, E., Damesin, C., Dufrene, E., Epron, D., Gross, P., Lebaube, S., Dantec, V.L., Goff, N.L., Lemoine, D., Lucot, E., Ottorini, J.M., Pontailler, J.Y., and Saugier, B. 2000. The carbon balance of a young beech forest. Functional Ecology. 14:312-325. Grant, R. F., Black, T. A., Humphreys, E. R., and Morgenstern, K. 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 Physiology 27:115-131. Green and Klinka. 1994. A field guide to site identification and interpretation for the vancouver forest region/by RN Green and K. Klinka. Ministry of Forests, Research Program. Humphreys, E. R., Black, T.A., Morgenstern, K., Cai, T., Drewitt, G.B., Nesic, Z., and Trofymow, J.A. (Tony). 2006. Carbon dioxide fluxes in coastal Douglas-fir stands at different stages of development after clearcut harvesting. Agricultural and Forest Meteorology. 140:6-22. Janisch, J. E., M. E. Harmon, H. Chen, B. Fasth, and J. Sexton. 2005. Decomposition of coarse woody debris originating by clearcutting of an old-growth conifer forest. Ecoscience. 12:151–160. Jassal, R. S., Black, T.A., Trofymow, J.A. (Tony), Roy, R., Nesic, Z. 2009. Forestfloor CO2 and N2O flux dynamics in a nitrogen-fertilized Pacific Northwest Douglas-fir stand. Soil Biology and Biochemistry, submitted.  98  Jungen, J. R. 1985. Soils of southern Vancouver Island. Page 209. MOE technical report 17, Ministry of Environment, Victoria, BC. Kominami, Y., Jomura, M., Dannoura, M., Goto, Y., Tamai, K., Miyama, T., Kanazawa, Y., Kaneko, S., Okumura, M., Misawa, N., Hamada, S., Sasaki, T., Kimura, H., and Ohtani, Y. 2008. Biometric and eddy-covariance-based estimates of carbon balance for a warm-temperate mixed forest in Japan. Agricultural and Forest Meteorology. 148:723-737. Krishnan, P., Black, T. A., Jassal, R. S., Chen, B., and Nesic, Z. 2009. Interannual variability in the carbon balance of three different-aged Douglas-fir stands in the Pacific Northwest. Journal of Geophysical Research (In press). Kurz, W., and Apps, M. 2006. Developing Canada's national forest carbon monitoring, accounting and reporting system to meet the reporting requirements of the Kyoto Protocol. Mitigation and Adaptation Strategies for Global Change. 11:33-43. Kurz, W., Dymond, C., White, T., Stinson, G., Shaw, C., Rampley, G., Smyth, C. Simpson, B., Neilson, E., Trofymow, J.A. (Tony), Metsaranta, J., and Apps, M. 2009. CBM-CFS3: A model of carbon-dynamics in forestry and landuse change implementing IPCC standards. Ecological Modelling. 220:480504. Landsberg, J., and Waring, R.. 1997. A generalised model of forest productivity using simplified concepts of radiation-use efficiency, carbon balance and partitioning. Forest Ecology and Management. 95:209-228. Li, Z., Kurz, W.A., Apps, M.J., and Beukema, S.J. 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. Canadian Journal of Forest Research. 33:126–136. MacKinnon, A. 2003. West coast, temperate, old-growth forests. Forestry Chronicle. 79:475-484. Mahalanobis, P.C. 1936. On the generalized distance in statistics. Proceedings of the National Institute of Science of India. 12. pp. 49-55. Matthews, G., and Great Britain. Forestry Commission. 1993. The carbon content of trees / George Matthews. Pages iv, 21 p.. Forestry Commission, Edinburgh. Miller, S., Goulden, M., Menton, M., da Rocha, H., de Freitas, H., Figueira, A., and de Sousa, C. 2004. Biometric and micrometeorological measurements of tropical forest carbon balance. Ecological Applications. 14:114-126.  99  Morgenstern, K., T. Black, T.A., Humphreys, E.R., Griffis, T.J., Drewitt, G.B., Cai, T., Nesic, Z., Spittlehouse, D.L., and Livingston, N.J. 2004. Sensitivity and uncertainty of the carbon balance of a Pacific Northwest Douglas-fir forest during an El Nino/La Nina cycle. Agricultural and Forest Meteorology. 123:201–219. National Forest Inventory Taskforce 2009. Canada's National Forest Inventory Ground Sampling Guidelines. Retrieved from: https://nfi.nfis.org/documentation/ground_plot/Gp_guidelines_v4.1.pdf. Schmid, H. P. 2002. Footprint modeling for vegetation atmosphere exchange studies: a review and perspective. Agricultural and Forest Meteorology. 113:159-183. Schmid, H. P., and Lloyd, C.R. 1999. Spatial representativeness and the location bias of flux footprints over inhomogeneous areas. Agricultural and Forest Meteorology. 93:195-209. Schreuder, H. T., Ernst, R., and Ramirez-Maldonado, H. 2004. Statistical techniques for sampling and monitoring natural resources. U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station. RMRSGTR-126:111. Schulze, E. D., Hogberg, P., Van Oene, H., Persson, T., Harrison, A.F., Read, D., Kjoller, A., and Matteucci, G. 2000. Interactions between the carbon and nitrogen cycle and the role of biodiversity: A synopsis of a study along a North-south transect through Europe. Pages 468–492 in Shulz, E. D. (ed.) Carbon and nitrogen cycling in European forest ecosystems. Springer Verlag, Heidelberg, Heidelberg. Stage, A. R., and Crookston, N.L. 2007. Partitioning error components for accuracy-assessment of near-neighbor methods of imputation. Forest Science. 53:62-72. Theede, A. D. 2007. Biometric and eddy-covariance estimates of ecosystem carbon storage at two boreal forest stands in Saskatchewan: 1994-2004. MSc Thesis: University of Saskatchewan. Trofymow, J. A., Stinson, G., and Kurz, W.A. 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. Forest Ecology and Management. 256, 1677-1691.  100  5. Conclusion A thorough understanding of the temporal and spatial forest C dynamics depends on compiling complementary data from a variety of sources, since each approach has its own advantages and disadvantages. Lidar remote sensing data offers the ability to accurately estimate the spatial distribution of aboveground live tree and snag mass, and this thesis provides some practical considerations for this task. To estimate the spatial distribution of other components of forest C stocks, remote sensing data may be used together with GIS data. Spatial scaling exercises such as the most-similar-neighbour classification of the EC flux-tower footprint are necessary in order to be able to compare biometric and EC flux-tower data sources over the same spatial scale and extent. An advantage of analysing biometric data sources is the ability to provide estimates of changes by component, while an advantage of analysing EC flux-tower-based data is the ability to measure the complete stand level flux. Where individual components of forest C stocks are difficult to measure accurately with biometric methods, modelling approaches may be useful to make informed estimates.  5.1 Lidar Estimation of Biomass Tree and snag above ground mass (TSAM) was estimated at the DF1949 study area using a combination of individual tree and plot-level lidar metrics. Most other biomass studies use plot-level lidar metrics for stand level estimates (e.g. Naesset, 2004; Lefsky et al., 2005), or individual tree metrics alone for measurement of individual tree biomass (eg. Popescue, 2007), but do not combine the two. In this  101  thesis, using individual tree metrics did not improve plot-level estimates compared to plot-level lidar metrics alone.  A study by Falkowski et al., (2008) shows that for individual tree detection studies in stands with dense canopies, errors of omission are common, including nonidentification of subdominant stems, and misidentification of multiple stems as single units. Further, the individual tree detection algorithm used was not calibrated for the detection of snags, and therefore did not improve the estimate. Accurate detection of individual trees, requires high return density lidar data (Anderson et al., 2006), which is more expensive to collect. Therefore, it may be more  economical for studies seeking stand level estimates to collect lower return density lidar data and focus on obtaining adequate return density for accurate estimates of target parameters, such as TSAM, rather than detection of individual trees.  5.2 Spatial Extrapolation of Ground Plot Data Advanced remote sensing and GIS predictor variables were used to estimate the spatial distribution of a variety of biometric variables across the extent of the footprint estimate using a most similar neighbour classification at DF1949. Since lidar provides a reliable estimate of TSAM, it was compared to the GIS estimate as an independent validation. An advantage of this approach was the Mahalanobis distance associated with the most-similar-neighbour classification was a useful indicator of how representative ground plots were of site conditions. Where Mahalanobis distances were low, ground plots were most representative of site  102  conditions. Agreement between the lidar and GIS models for TSAM was highest near the EC flux-tower, where plots were most representative of site conditions, and areas with less agreement were close to the edges of the footprint, where the footprint flux probability is also lower. This gives us confidence using the GIS stratification-based estimate of the spatial distribution of C stocks. The use of footprint weighting made a large difference in the calculation of site-level C stocks.  5.3 Comparing Biometric-Based and EC Flux-Tower-Based Measurements In contrast to the C stock extrapolations, footprint weighting made less of a difference for estimated changes in C over the measurement interval (∆C stocks), since ∆C stocks was lower and less spatially variable than the original C stock distribution. Certain elements of forest C stocks were difficult to measure consistently (for example, forest floor, mineral soil, and roots), and in these cases, modelling provided an important tool to achieve a reliable estimate of the change in C stocks over the time period. We hypothesise that an important source of respiration is responsible for the difference between biometric ∆C stocks and cumulative EC flux-tower-based measurements of NEP (ΣNEP) at the more recently disturbed sites which was not measured in this study, and an important future research goal is to identify the reason for this difference.  This thesis shows that a thorough analysis of forest C stocks and C stock changes requires integrating multiple approaches, including remote sensing, biometric inventory, C budget modelling, and EC flux-tower-based observations. This is also  103  illustrated in other C accounting efforts, such as Canada's National Forest Carbon Monitoring, Accounting, and Reporting System (NFCMARS), where one national forest inventory measurement is used with CBM-CFS modeling of change (Kurz and Apps, 2006). When future re-measurements of NFI sample plots are available, Canadian National forest C accounting may rely more heavily on biometric remeasurement data, and where difficulties are encountered, due to measurement inconsistencies between different field crews or spatial variability within plots, a combined measurement and modelling approach may be utilized, similar to the methods in this thesis.  When combining measurements from multiple sources, such as EC flux-towers, biometric inventory plots, and regional modelling efforts, spatial and temporal scaling must be considered so that comparison can be made across the same temporal and spatial extents. In the NFCMARS example, spatial scaling is accomplished in the analysis of National Forest Inventory photo plots, where land and forest cover is sampled and extrapolated to the landscape scale (Gillis et al. 2006). In the present study, spatial scaling was accomplished at a very high level of spatial detail using advanced remote sensing and GIS data sources and classification algorithms. At the British Columbia Flux Station, which is a highly intensive study area, the availability of a combination of data from remote sensing, GIS, EC flux-tower instrumentation, and biometric inventory plots makes it possible to try and test innovative approaches from multiple data sources. Information from these intensive studies may potentially be considered for broader regional projects.  104  Comparing biometric ∆C stocks and tower ΣNEP is currently a priority for the Canadian Carbon Program (CCP) (Trofymow et al., 2008), with significant attention to quantifying variability of stand structure at EC flux-tower sites. This concern is addressed at the British Columbia Flux Station in this thesis by identifying environmental predictor variables from GIS and remote sensing data sources that are correlated with biometric target variables, and using statistical distances, such as the Mahalanobis distance, to quantify how representative sample plots are of stand structure. At other sites, researchers may also be successful with similar data sources and methodologies.  Estimates of C fluxes from EC flux-tower-based data and biometric-based data provided complementary insights into forest C processes. For example, data from biometric sample plots provided data about changes by component, while EC fluxtower-based data suggest that an important component of respiration is not being measured at the more recently disturbed sites. Direct measurement of woody debris, stump, snag, and root decomposition are highly research intensive (e.g. Jansich et al. 2005), and thus very costly. Given the currently available range of allometric equations, it is difficult to estimate belowground coarse root mass for stumps. If these values were available, decomposition estimates and transfers from CBM-CFS3 could be used to model this source of respiration, which is a missing component in this study. Future research reducing uncertainty for biometric estimation of dead coarse root C fluxes, may improve this estimate, and furthermore, may have direct forest harvest practice policy implications.  105  5.4 References Andersen, H., Reutebuch, S., and McGaughey, R. 2006. A rigorous assessment of tree height measurements obtained using airborne lidar and conventional field methods. Canadian Journal of Remote Sensing. 32:355-366. Falkowski, M. J., Smith, A.M.S., Gessler, P.E., Hudak, A.T., Vierling, L.A., and Evans, J.S. 2008. The influence of conifer forest canopy cover on the accuracy of two individual tree measurement algorithms using lidar data. Canadian Journal of Remote Sensing. 34:1-13. Gillis, M., Omule, A., and Brierley, T. 2005. Monitoring Canada's forests: The National Forest Inventory. Forestry Chronicle. 81:214-221. Janisch, J. E., Harmon, M.E., Chen, H., Fasth, B., and Sexton, J. 2005. Decomposition of coarse woody debris originating by clearcutting of an old-growth conifer forest. Ecoscience. 12:151–160. Kurz, W., and Apps, M. 2006. Developing Canada's National Forest Carbon Monitoring, Accounting and Reporting System to Meet the Reporting Requirements of the Kyoto Protocol. Mitigation and Adaptation Strategies for Global Change. 11:33-43. Lefsky, M. A., Hudak, A.T., Cohen, W.B., and Acker, S. 2005. Geographic variability in lidar predictions of forest stand structure in the Pacific Northwest. Remote Sensing of Environment. 95:532-548. Næsset, E. 2004. Practical large-scale forest stand inventory using a small-footprint airborne scanning laser. Scandinavian Journal of Forest Research 19:164. Popescu, S. C. 2007. Estimating biomass of individual pine trees using airborne lidar. Biomass and Bioenergy. 31:646-655. Trofymow, J. A., Stinson, G., and Kurz, W.A. 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. Forest Ecology and Management. 256, 1677-1691.  106  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0067765/manifest

Comment

Related Items