Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Remotely sensed image time series to inform on forest structure and post-disturbance recovery over Canadian… Frazier, Ryan James 2016

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

Item Metadata

Download

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

Full Text

REMOTELY SENSED IMAGE TIME SERIES TO INFORM ON FOREST STRUCTURE AND POST-DISTURBANCE RECOVERY OVER CANADIAN BOREAL FORESTS by  Ryan James Frazier  B.A. Clark University 2006 M.A. Clark University 2007  A THESIS SUBMITTED IN PARTIAL FULFILLMENT  OF THE REQUIREMENTS FOR THE DEGREE OF   Doctor of Philosophy in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (FORESTRY)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)   August 2016  ©  Ryan James Frazier, 2016 ii  Abstract Over the last century Canadian boreal forests have warmed by 2-3° C, causing growing seasons to lengthen and alterations to annual productivity, which result in numerous responses from boreal tree species. Both disturbance and recovery cycles are affected, although change in northern Canadian boreal forests is difficult to detect, since they remain non-inventoried and lack detailed spatially explicit descriptive data. Through the use of Landsat time series, remote sensing offers the ability to map and monitor large forested areas over time to provide valuable information about boreal forests. The overall objective of this dissertation is to assess the capacity of remotely-sensed spectral time series to characterize forest recovery following disturbance in Canadian boreal forests Major findings produced from the research presented in this dissertation show: • Boreal forest attributes are better estimated with Landsat time series metrics than single date information, and the inclusion of recovery metrics substantially improves accuracy • Choice of spectral index to monitor recovery is important, and the use of multiple spectral indexes can provide better and meaningful insights into forest recovery • The East/West division of the Boreal Shield ecozone is reinforced due differing spectral forest recovery trajectories that are suggestive of distinct recovery processes in each region.  • Forest recovery rates are not fixed across the Boreal and Taiga Shield ecozones, with Taiga Shield spectral forest recovery rates showing a consistent positive trend, possibly indicating forests are recently recovering at an accelerated rate.  iii  The research presented in this dissertation advances the use of remote sensing to detect post-disturbance recovery in boreal forest ecosystems. Monitoring large forested areas such as the boreal for change is increasingly important as climatic conditions alter, and the spectral time series methods shown herein provide new tools to observe change in boreal forests. Future research directions are identified around first lengthening time series across longer periods of time, then extending these spectral time series approaches across jurisdictional lines in the pan-boreal region, and finally incorporating the data generated from these methods to be incorporated into carbon accounting frameworks. iv  Preface The research questions and objectives of this dissertation were conceived by me and are the direct result of discussions between myself and my supervisory committee. Portions of the research presented in this dissertation appear as co-authored publications in peer reviewed journals. For those portions that are published, I conducted the research, data analysis and subsequent interpretation of results, and prepared the final submitted manuscript:  • Chapter 4: Frazier, R. J., Coops, N. C., Wulder, M. A., & Kennedy, R. (2014). Characterization of aboveground biomass in an unmanaged boreal forest using Landsat temporal segmentation metrics. ISPRS Journal of Photogrammetry and Remote Sensing, 92, 137-146. • Chapter 5: Frazier, R. J., Coops, N. C., & Wulder, M. A. (2015). Boreal Shield forest disturbance and recovery trends using Landsat time series. Remote Sensing of Environment, 170, 317-327. • Chapter 6: Frazier, R. J., Coops, N. C., & Wulder, M. A., Hermosilla, T., White, J.C. (2016). Spatial and Temporal Analysis of Boreal Forest Spectral Recovery Trends from Dense Landsat Time Series.  R. J. F. Faculty of Forestry The University of British Columbia v  Table of Contents Abstract .......................................................................................................................................... ii  Preface ........................................................................................................................................... iv  Table of Contents ...........................................................................................................................v List of Tables ................................................................................................................................ ix List of Figures ............................................................................................................................... xi List of Abbreviations ................................................................................................................. xiv Acknowledgements .................................................................................................................... xvi Dedication ................................................................................................................................. xviii  Introduction, Research Needs, Questions, and Dissertation Overview...................1 Chapter 1:1.1 Introduction ..................................................................................................................... 1 1.2 Research Needs ............................................................................................................... 7 1.3 Research Questions ......................................................................................................... 8 1.4 Dissertation Overview .................................................................................................... 9  Study Area & Data .....................................................................................................13 Chapter 2:2.1 The Canadian Boreal Zone ........................................................................................... 13 2.2 Ecozone Study Areas .................................................................................................... 15 2.3 Landsat Data ................................................................................................................. 19  Canadian Boreal Change as Assessed with Remotely Sensed Datasets ................26 Chapter 3:3.1 Coarse Spatial Resolution Monitoring .......................................................................... 26 3.2 Landsat Based Medium Spatial Resolution Monitoring ............................................... 30  Spectral Time Series Metrics for Forest Structure Prediction ..............................39 Chapter 4:4.1 Introduction ................................................................................................................... 39 vi  4.2 Methods......................................................................................................................... 43 4.2.1 Study Area ................................................................................................................ 43 4.2.2 Datasets ..................................................................................................................... 44 4.2.2.1 LiDAR Data ...................................................................................................... 44 4.2.2.2 Ancillary Environmental Data .......................................................................... 46 4.2.2.3 Landsat Data ..................................................................................................... 47 4.2.3 Image Processing ...................................................................................................... 49 4.2.3.1 Temporal Segmentation .................................................................................... 49 4.2.3.1.1 Cloud Masking ............................................................................................ 49 4.2.3.1.2 Radiometric Correction and Normalization ................................................ 50 4.2.3.1.3 Yearly Cloud Free Composites ................................................................... 50 4.2.3.1.4 Tasseled Cap Components .......................................................................... 51 4.2.3.1.5 LandTrendr .................................................................................................. 53 4.2.3.1.6 LandTrendr Spectral Time Series Metrics .................................................. 55 4.2.3.1.7 Metric Grouping .......................................................................................... 56 4.2.3.2 Spatial Segmentation ........................................................................................ 58 4.2.3.3 Random Forests ................................................................................................ 59 4.3 Results ........................................................................................................................... 60 4.4 Discussion ..................................................................................................................... 64 4.4.1 Understanding Variance and RMSE in Context with Previous Studies ................... 64 4.4.2 Metric Complexity .................................................................................................... 66 4.4.3 Variable Importance.................................................................................................. 67  Spectral Recovery Trajectories Illuminate Differences in Forest Recovery ........69 Chapter 5:vii  5.1 Introduction ................................................................................................................... 69 5.2 Study Area .................................................................................................................... 71 5.3 Methods......................................................................................................................... 73 5.3.1 Data ........................................................................................................................... 73 5.3.1.1 Landsat .............................................................................................................. 73 5.3.2 Image Processing ...................................................................................................... 75 5.3.2.1 Image Selection, Preprocessing, and Compositing ........................................... 76 5.3.2.2 Spectral Indices and Recovery .......................................................................... 78 5.3.2.3 Temporal Segmentation .................................................................................... 80 5.3.2.4 Cohen’s d statistic ............................................................................................. 81 5.4 Results ........................................................................................................................... 83 5.4.1 Wetness Recovery Trajectories and Cohen’s d Values ............................................ 83 5.4.2 Greenness Recovery Trajectories and Cohen’s d Values ......................................... 86 5.4.3 Brightness Recovery Trajectories and Cohen’s d Values ......................................... 88 5.4.4 Areal Disturbance and Recovery Trends .................................................................. 90 5.5 Discussion ..................................................................................................................... 93 5.5.1 Temporal Trajectory Considerations & Recovery Characterization ........................ 94 5.5.2 Case Study: Partitioning of the East and West Boreal Shield .................................. 97  Spectral Recovery Rates are Not Fixed over Time and Space .............................101 Chapter 6:6.1 Introduction ................................................................................................................. 101 6.2 Study Area .................................................................................................................. 105 6.3 Data & Methods .......................................................................................................... 108 viii  6.3.1 Spectral Data ........................................................................................................... 108 6.3.2 Methods................................................................................................................... 109 6.3.2.1 Data Stratification and Filtering...................................................................... 109 6.3.2.2 Spectral Forest Recovery Metrics ................................................................... 110 6.3.2.3 Trend Detection .............................................................................................. 118 6.4 Results ......................................................................................................................... 119 6.4.1 Non-Recovery Trends, Fire Frequency and Extent ................................................ 119 6.4.2 100 km Grid Spectral Forest Recovery Rate Trend ................................................ 122 6.4.3 Ecozone-Level Analysis ......................................................................................... 129 6.5 Discussion ................................................................................................................... 131 6.5.1 Spectral Forest Recovery Rate Trends .................................................................... 131 6.5.2 Spectral Forest Recovery Metrics ........................................................................... 131 6.5.3 The Need for Continued Boreal Monitoring ........................................................... 134  Conclusions ...............................................................................................................136 Chapter 7:7.1 Research Innovations .................................................................................................. 136 7.2 Methods for Charactering Boreal Forest Recovery with Spectral Time Series .......... 136 7.3 Application of Spectral Forest Recovery for Data Generation in Boreal Forests....... 139 7.4 Insights into Canadian Boreal Forest Dynamics ......................................................... 141 7.5 Limitations .................................................................................................................. 143 7.6 Future Research & New Questions ............................................................................. 146 References ...................................................................................................................................148 ix  List of Tables Table 2.1: Summary of biotic and abiotic factors of Boreal Cordillera, Boral Shield, and Taiga Shield ecozones. ............................................................................................................................ 17 Table 2.2: Spatial, spectral  and temporal characteristics of Landsat sensors. ............................. 23 Table 3.1: Coarse spatial resolution remote sensing based vegetation studies summarized and pertinent Canadian forest results highlighted ............................................................................... 28 Table 3.2: Landsat based time series analysis of forests summarized with recovery aspects highlighted .................................................................................................................................... 35 Table 3.3: Frequently used spectral vegetation index calculation equations. ............................... 37 Table 3.3: Frequently used spectral vegetation index calculation equations. ............................... 38 Table 4.1: Explained variance (R2) and Root Mean Square Error (RMSE) values of LiDAR transect predicted forest attributes from the 25x25 meter plots database. .................................... 46 Table 4.2: Landsat image dates and years used to construct spectral time series. ........................ 48 Table 4.3: Variable groups used in Random Forest analysis to predict forest attributes .............. 57 Table 4.4: General characteristics of forest units produced from multilevel object oriented segmentation ................................................................................................................................. 59 Table 4.5: Explained variance (R2) and Root Mean Square Error % (RMSE%) values of random forests models used to predict aboveground ................................................................................. 61 Table 4.6: Mean variable importance values for random forests models ..................................... 62 Table 4.7: Explained variance (R2) and Root Mean Square Error % (RMSE%) output values of random forests models. ................................................................................................................. 63 Table 5.1: Summary statistics about Landsat scenes used to create spectral time series. ............ 75 x  Table 5.2: Summarization of Cohen’s d recovery trajectories based on Brightnes, Greennes, Wetness and comparison trajectories. ........................................................................................... 85 Table 6.1: Boreal and Taiga Shield East and West study area in respect to the Canadian boreal zone as defined by Brandt et al., (2009), .................................................................................... 106 Table 6.2: Global Moran’s I spatial autocorrelation results indicating that dispersed, random, or clustered patterns exist ................................................................................................................ 129  xi  List of Figures Figure 1.1: Conceptual diagram of boreal forest dynamics and spectral time series. .................... 9 Figure 2.1 The boreal zone of Canada (Brandt 2009) highlighted against a false color imagery composite ...................................................................................................................................... 14  Figure 2.2 Map of Canadian ecozones with management areas displayed; highlighted zones indicate where research was performed ........................................................................................ 16 Figure 2.3: Timeline of Landsat satellites and sensors illustrating the temporal depth of the archive. .......................................................................................................................................... 21  Figure 4.1: Map showing study area location and other characteristics near Watson Lake, YK. 44 Figure 4.2: Spectral time series, temporal segmentation, and derived metrics for a hypothetical forest unit. ..................................................................................................................................... 54  Figure 4.3: Aboveground Biomass explained variance (R2) and Root Mean Square Error % (RMSE%) for random forests models using scenarios 1-10. ........................................................ 61 Figure 4.4: Explained variance (R2) and Root Mean Square Error % (RMSE%) for random forests models estimating each forest attribute using scenario 11. ............................................... 64 Figure 5.1: Boreal Shield East and West ecozones and Canadian forest management zones within the eight study scene portions. ...................................................................................................... 73 Figure 5.2: Distribution of dates and years of imagery for each Landsat footprint portion used for study in Chapter 5. ........................................................................................................................ 77 Figure 5.3: Boreal Shield East and West Tasseled Cap Wetness recovery and Cohen’s d trajectories. .................................................................................................................................... 84 Figure 5.4: Boreal Shield East and West Tasseled Cap Greenness recovery and Cohen’s d trajectories. .................................................................................................................................... 87 xii  Figure 5.5: Boreal Shield East and West Tasseled Cap Brightness recovery and Cohen’s d trajectories. .................................................................................................................................... 89 Figure 5.6: Annual amounts of disturbed, and continuously unabated recovering pixels for a) eastern study areas and b) western study areas. ............................................................................ 91 Figure 5.7: Boreal Shield East study areas spatial disturbance and recovery trends shown as a percent of the total forest area in the study areas. ......................................................................... 92 Figure 5.8: Boreal Shield West study areas spatial disturbance and recovery trends shown as a percent of the total forest area in the study areas. ......................................................................... 93 Figure 6.1: The Boreal and Taiga Shield East and West study area is displayed with a false color Landsat composite (R; Band 5, G: Band 4, B: Band 3). ............................................................. 106 Figure 6.2: Example of spectral data showing pre-disturbance NBR values, Fire Year, and dNBR values for three sites in the Taiga Shield West. .......................................................................... 110 Figure 6.3: NBR time series from the year of disturbance through five years post disturbance. 112 Figure 6.4: The Recovery Indicator (RI), adapted from Kennedy et al., (2012). ....................... 114 Figure 6.5: The Ratio 80% of Pre-disturbance (R80P) metric, adapted from Pickell et al., (2016)...................................................................................................................................................... 116  Figure 6.6: The Year on Year Average (YrYr) metric. .............................................................. 117 Figure 6.7: Sample of spectral forest recovery rate metrics calculated for the same areas shown in Figure 6.2 and 6.3. ...................................................................................................................... 118 Figure 6.8: Presence of significant trends over time in either pre-disturbance values or disturbance severity in the 100 km grid cells.............................................................................. 120 Figure 6.9: The number of years a cell experienced a moderate or high severity fire disturbance within the time series. ................................................................................................................. 121 xiii  Figure 6.10: The percentage of each grid cell disturbed by a moderate or high severity fire within the time series. ............................................................................................................................ 122 Figure 6.11: Significantly trending grid cells and Theil-Sen slopes for the Recovery Indicator (RI) metric time series. ............................................................................................................... 123 Figure 6.12: Significantly trending grid cells and Theil-Sen slopes for the Ratio 80% of Pre-disturbance (R80P) metric time series. ....................................................................................... 125 Figure 6.13: Significantly trending grid cells and Theil-Sen slopes for the Year on Year Average (YrYr) metric time series. ........................................................................................................... 127 Figure 6.14: Number of significant trends in a grid cell, across all spectral forest recovery rate metric time series Mann-Kendall results. ................................................................................... 128 Figure 6.15: Ecozone analysis level spectral forest recovery rate trend detection results. ......... 130    xiv  List of Abbreviations AERO: Aerosol Band AK: Alaska AVHRR: Advanced Very High Resolution Radiometer BAP: Best Available Pixel CMFDA: Continuous Monitoring of Forest Disturbance Algorithm DEM: Digital Elevation Model dNBR: Difference Normalized Burn Ratio DOY: Day of Year EOSD: Earth Observation for Sustainable Development of Forests ETM+: Enhanced Thematic Mapper Plus Fmask: Function of Mask LEDAPS: Landsat Ecosystem Disturbance Adaptive Processing System LiDAR: Light Detecting and Ranging MADCAL: Multivariate Alteration Detection Calibration Algorithm  ME: Maine MSS: Multi Spectral Scanner NBR: Normalized Burn Ratio NBR y-1: Normalized Burn Ratio values one year before disturbance NBRpre: Normalized Burn Ratio values before disturbance NBRy+1: Normalized Burn Ratio values on year after disturbance NBRy+2: Normalized Burn Ratio values two years after disturbance NBRy+3: Normalized Burn Ratio values three years after disturbance NBRy+4: Normalized Burn Ratio values four years after disturbance NBRy+5: : Normalized Burn Ratio values five year after disturbance NBRy0: Normalized Burn Ratio values in the year of disturbance NBRy-2: Normalized Burn Ratio values two years before disturbance ∆NBRdisturbance: The magnitude change in NBR value due to disturbance ∆NBRregrowth: The magnitude change in NBR values after disturbance NDMI: Normalized Difference Moisture Index NDVI: Normalized Difference Vegetation Index NFI: National Forest Inventory NIR: Near Infrared NWT: Northwest Territories OLI: Operational Land Imager PAN: Panchromatic QC: Quebec R80P: Ratio of Eighty Percent of Pre-Disturbed Values RI: Recovery Indicator RMSE%: Root Mean Square Error Percentage RMSE: Root Mean Square Error SE: Southeast SPOT: System Pour le Observation de la Terre  SVI(s): Spectral Vegetation Index(es) xv  SWIR: Shortwave Infrared TCB: Tasseled Cap Brightness TCG: Tasseled Cap Greenness TCW: Tasseled Cap Wetness TIR: Thermal Infrared TM: Thematic Mapper USA: United States of America USGS: United States Geological Survey UTM: Universal Transverse Mercator VCT: Vegetation Change Tracker VGT: Vegetation sensor onboard System Pour le Observation de la Terre  WRS-2: Worldwide Reference System -2 YK: Yukon Territory YrYr: Year on Year Average   xvi  Acknowledgements  Many thanks go to my committee - I’d be much poorer without their input. Nicholas Coops is singled out for thanks due to his exceptional guidance throughout my PhD career; at every phase he knew what should be done, and always offered a good bit of advice. Special thanks go to Mike Wulder for many useful conversations about Landsat and boreal forests; his expert opinions on scientific writing were invaluable. None of this research would have been possible without the University of British Columbia benefiting me a Four Year Fellowship. Additionally, the United States Geological Survey is thanks for maintaining the Landsat program and the data created by Landsat sensors since 1972; without that data this research would have been impossible.  All past and present IRSS lab members deserve to receive thanks as well, as they either lent a hand when asked for help, or an ear when asked for support! Special thanks go to Txomin Hermosilla for his mastery of and insights into Landsat data processing, but for also being a reliable friend. Doug Bolton is also worthy of thanks for being a positive force in the lab. A dependable friend, Wiebe Nijland is thanked for starting and then entrusting me to carry the lunch walk tradition on in his absence. I will always appreciate the good camaraderie-fostering conversation that we all had as lab mates on lunch walk. Many thanks go to Trevor Jones- from Worcester to Vancouver, HERO to the IRSS, you always have been a great source of advice! xvii  Many thanks go to my previous teachers, coworkers, instructors, and professors; they provided me a solid base of knowledge from which I could delve deep into research. Enduring gratitude goes to John Rogan for introducing me to Landsat 5 and inspiring with its capabilities! Most importantly, thanks must go to my family and friends who have supported me throughout this endeavor - without their help and encouragement, none of this would have been possible. Thanks to Beauregard for his calming nature and never ending positive demeanor. My mothers, father, and brother are all thanked for their endless positive encouragement, which I have relied upon for so long. Deserving of never ending thanks, my wife made the long hard road of completing this dissertation much easier. Without her continual support, encouragement, and loving care, I would have never finished this work, and likely never applied myself fully. Simply put, I could not have asked for a better partner to journey through life with! Finally, I thank my son, Henry, for always inspiring me to ask more questions.   xviii  Dedication        to my two loves, Henry and Eynat1   Introduction, Research Needs, Questions, and Dissertation Chapter 1:Overview 1.1 Introduction Boreal forest comprise 27% of the world’s forest cover (FAO, 2001), with approximately 88% of the North American boreal located in Canada (Brandt, 2009). Latitudinal location and climatic conditions limit annual growth of the slow-growing coniferous tree species that are present and interspersed with varieties of cold tolerant deciduous tree species (Hytteborn et al., 2006; Kimball et al., 2004; Jarvis et al., 2000). The importance of Canadian boreal forests lies not only in their extent, but also in services and goods they provide. The managed boreal forests of Canada store an estimated 28 Pg of carbon (Kurz et al., 2013), and the pan-boreal forests hold an estimated 32% of world carbon stocks (Pan et al., 2011). The timber in Canadian boreal forests is a highly valued resource and as a result, Canada is an important producer of roundwood and sawnwood (FAO 2011). Additionally, products derived from the boreal forests encompass 13% of gross export earnings in Canada (Nilsson 2000). Additionally, Canadian boreal forests provide essential habitats for a diverse range of animal and plant species (Kerr et al., 2003). Furthermore, boreal forests also play an important role: 1) affecting volume,  timing and purification of the freshwater that discharges into the Arctic Ocean (Powers et al., 2013; Ruckstuhl et al., 2008; Woo et al., 2008); 2) global aboveground and below ground carbon storage (Kurz et al., 2013; Shvidenko et al., 2006); 3) global and regional climate regulation through the absorption of incoming radiation and warming of local temperatures (Brandt et al., 2013; Brandt 2009; Bonan 2008). 2  Boreal forests will be altered extensively by a changing climate (Gauthier et al., 2014; Price et al., 2013; Brandt et al., 2013), with average annual temperatures already rising 2-3° C over the last century throughout the region and contributing towards broad environmental changes (Stocks et al., 1998). These new climatic conditions will modify disturbance and growth dynamics. Changes in growing season length (Luckman 2004, Xu et al., 2013), seasonal phenology (Lemprière et al., 2008;  Colombo 1998), annual productivity (Nemani et al., 2003;  Jarvis et al., 2000), and disturbance characteristics (Flannigan et al., 2005; Flannigan et al., 2000; Stocks et al., 1998) are all underway throughout Canadian boreal forests. Fire-dominated disturbance cycles that are typical of boreal forest ecosystems are also affected by the increase in temperature. Climatic changes can alter the frequency, severity, and extent of boreal fires (Gillet et al., 2004; Amiro et al., 2001), and extend fire seasons beyond their usual length (Stocks et al., 2002). This changing climate has already produced mixed effects on annual forest growth and post-disturbance recovery. Research has shown increases in productivity and annual growth in boreal forests are occurring unevenly in both magnitude and spatial distribution (Boisvenue et al., 2006). A number of widespread boreal tree species have shown a variety of responses to a changing climate. Aspen forests in Saskatchewan have experienced earlier leaf emergence due to warmer mean annual temperatures (Chen et al., 1999). Black spruce (Picea mariana) dominated systems showed an increase in annual growth under cooler and wetter conditions, while jack pine (Pinus banksiana) ecosystems displayed increased growth with warmer temperatures and greater spring precipitation (Brooks et al., 1998). Additionally, white spruce (Picea glauca) trees have reacted positively to warmer spring temperatures with increased annual growth, but also 3  negatively to hotter summers with decreased annual growth (Wilmking et al., 2004) likely caused by a lack of available moisture (Barber et al., 2000). Nevertheless, it is difficult to more comprehensively assess the widespread effects of a changing climate on Canadian boreal forests due to a lack of spatial information describing northern forests (Price et al., 2013). Canadian boreal forests are subject to varying levels of management and inventory practices (Wulder et al., 2007), with southern boreal forests typically subject to inventory on a periodic basis in support of sustainable forest management and planning objectives. These southern managed forests are generally where most of the human population resides, road access is greater, and vegetative productivity is higher. In contrast, the less productive northern Canadian boreal forests are more remote and not typically utilized for harvest, and as a result less emphasis is placed on inventory and monitoring of those regions. Those unmanaged northern boreal forests are consequently not as well characterized as their southern counterparts (Andrew et al., 2012). Despite the paucity of data describing the condition of up to half of the Canadian boreal forests, international and internal reporting requirements must still be met. Forest inventory methods driven by in situ sample collection are expensive due to the high costs associated with establishing plots in remote areas that lack accessibility (Wulder et al., 2004). In support of establishing some data for those unmanaged areas, the National Forest Inventory (NFI) conducts limited inventory at a number of locations, and often using remotely-sensed imagery to develop some forest inventory attributes. For those northern boreal forests that lack spatially explicit data, a more efficient and cost effective method of acquiring comprehensive data must be employed. The need for disturbance mapping and revegetation monitoring is accentuated by increased 4  quantity and severity of wildfire disturbances in drier western boreal areas (Flannigan et al., 2005) and increased productivity in temperature limited predominantly eastern boreal areas (Boisvenue et al.,  2006), both due to a changing climate.  An ideal data source to generate new and useful information about large forested areas over long time periods is remotely-sensed imagery (Kennedy et al., 2014; Hansen et al., 2013; Zhu et al., 2012a; Kennedy et al., 2010; Huang et al., 2010; Kennedy et al., 2007;   ). The value of remotely-sensed data lies with its systematic observations that have captured changes over large areas, which are comparable through time. Remotely-sensed data can offer information over spatial scales that may not be otherwise achievable with in situ methods or other data sources. Remotely sensed Landsat imagery has proven effective for a wide variety of land surface applications, especially forest monitoring (Wulder et al., 2008a; Woodcock et al., 2008; Hall et al., 2006; Wulder et al., 2004b; Skakun et al., 2003). A beneficial balance between spatial (30 m), spectral (visible, near and shortwave infrared spectra), radiometric (cross comparable data), and temporal (approximately 16 day revisit period) resolutions struck by the Landsat sensor series creates an ideal data source to observe change in boreal forests (Kennedy et al., 2014), especially in the absence of any other extant descriptive forest data.  As the need to understand landscape dynamics has increased, the analysis of Landsat data has also shifted towards multi-temporal approaches to capture change over time (Hansen et al.,  2012; Wulder et al., 2011; Kennedy et al., 2010). To date there has been a marked emphasis on characterizing forest disturbance events using multi-temporal Landsat imagery (such as fire, harvest and other non-stand replacing disturbances) (Zhu et al., 2012a; Thomas et al., 2011;  Frolking et al., 2009; Huang et al., 2009;  Wulder et al., 2009;  Olsson 2009; French et al., 2008;  5  Masek et al., 2008;  Cohen et al., 2002 ). For example, The Canadian National Fire Database tracks fire across Canada and in part characterizes fire events using remotely-sensed datasets (Stocks et al., 2002). Insect infestations (Goodwin et al., 2008) and harvest (Masek et al., 2008; Schroeder et al., 2011) have also been mapped and characterized with remotely-sensed datasets. Thus, the main types of boreal forest disturbances have been well characterized with remotely-sensed datasets.  In direct contrast to the focus on boreal forest disturbance, relatively less attention has been placed on observing post-disturbance forest recovery spectrally (Griffiths et al., 2014; Kennedy et al., 2012; Gitas et al., 2012; Schroeder et al., 2011;  Chen et al., 2011). Thus there is a knowledge gap in observing and understanding how the forest recovers spectrally. Forest recovery is typically defined as the reestablishment of forest biomass or canopy structure following disturbance (Frolking et al., 2009; Oliver et al., 1996), making it a process rather than a state and therefore observable over time (Bartels et al., 2016; Hermosilla et al., 2015b; Kennedy et al., 2012;  Zhu et al., 2012a; Huang et al., 2010). The current knowledge of boreal forest dynamics describes forests capable of withstanding and requiring periodic disturbances, and largely recovering through multiple means of self-replacement. However observable forest recovery may be with remotely-sensed datasets, data describing forest recovery is virtually non-existent in the non-inventoried northern boreal forests and only available on a periodic basis in the southern managed forests; the data divide between northern and southern forests represents a gap in information and knowledge about Canadian boreal forests. Current theory suggests that the regrowth of forests after disturbance is a key process that must occur if the previous function and structure of forests are to be restored, and also provide 6  services critical to the ecosystem in which they inhabit. Information about forest recovery is important to collect, as it can be utilized to: 1) assess ecosystem health (Gauthier et al., 2015); 2) be used to create more accurate carbon accounting models (Kurz et al., 2013); and, 3) make better forest management decisions (Shatford et al., 2007). Moreover, the presence of forest recovery data for all Canadian boreal forests would facilitate greater understanding of the effects caused by a changing climate. For example , Kasischke et al., (2006) used fire disturbance data and showed that across the North American boreal forest between 1959 and 1999 annual burnt area doubled, large fires became more frequent, and that fire occurrence shifted more towards the end the growing season all thought to be driven by a changing climate. Although this information is valuable, of equal importance and much less explored is how those fire-affected forests have reacted to those changing aspects of the disturbance cycle. Information about forest recovery, available across the Canadian boreal, would enable a better understanding of how forests are adapting to the new set of factors presented by a changing climate.  The recovery of post-disturbed boreal forests could be monitored effectively with Landsat-based spectral time series. Spectral signatures are theorized to show that vegetation has a unique signal across the collected spectral bands of Landsat, and that over time these signatures can change with the vegetation type and age. For example, the spectral characterization of forest recovery through time can be used to monitor the progression of a stand towards successful reestablishment and maturation (Schroeder et al., 2011). This valuable information is important for forest managers and decision makers alike, (Frokling et al., 2009; Shatford et al., 2007) offering insights into how different regions are potentially changing, both in terms of pressures on disturbance regimes as well as the subsequent recovery (Powers et al., 2015; Pickell et al., 7  2014; Xu et al., 2013; Beck et al., 2011; Berner et al., 2011). Monitoring Canadian boreal forests with remotely-sensed datasets: 1) offers the capacity to examine large areas over time; 2) can reveal insights into forest recovery not evident with other spatially disparate datasets; and, 3) production of information relevant for researchers and decision makers alike. Lastly, forest recovery information is vital, because it is not enough to understand the effects of disturbance, but full comprehension of the ensuing recovery of Canadian boreal forests is important in order to unravel the full effects of a changing climate.   1.2 Research Needs Unfettered access to analysis ready remotely sensed imagery, when combined with new time series analysis techniques, offers the ability to improve upon existing forest data and current techniques and provide new insights into disturbance and recovery of boreal forests consistently over large areas through time. Recent research has established the ability of remotely sensed time series to characterize disturbance events and their impact on forests. However the nature of how spectral recovery derived from satellite data, such as from Landsat, relates to forest recovery is still a new field of active research, especially in boreal environments. Since remote sensing datasets are limited in their time span, space for time substitution is often used in a similar ecological manner as established in situ plots, but to what degree spectral recovery can report on forest recovery is unclear. Spectral signatures theory suggests that as a forest recovers, it should have a changing spectral signature. As a result, five research needs will be addressed in this dissertation: (1) improved understanding of the relationship between forest recovery and spectral recovery; (2) new methods for characterizing and comparing spectral recovery, (3) application of newly developed spectral recovery methods across larges areas and over time; and (4) generation 8  of spatially explicit spectral forest recovery data; and finally, (5)  derivation of  insights from the comparison of forest recovery over space and time to better understand forest conditions in Canadian boreal forests. 1.3 Research Questions The overall objective of this dissertation is to assess the capacity of remotely-sensed spectral time series to characterize forest recovery following disturbance in Canadian boreal forests. Three questions are posed to meet this objective: 1. Can remotely-sensed spectral time series be used to inform on boreal forest attributes? 2. Do spectral differences in recovery trajectories reflect distinct post-disturbed forest recovery patterns? 3. Are spectral forest recovery rates fixed across space and time?   These questions give a unique perspective on the use of remotely sensed data to monitor boreal forests, and how to characterize post-disturbance forest recovery with spectral time series as shown in Figure 1.1. The first question investigates how spectral time series metrics can be used to estimate forest attributes, especially utilizing the entire spectral time series. Only post-disturbance recovery is examined and characterized by Question 2 using spectral time series. Finally, Question 3 focuses on the initial period of stand reestablishment to inform on the rate at which that forest recovery occurs throughout a time series.    9   Figure 1.1: Conceptual diagram of boreal forest dynamics and spectral time series. Research Question 1 addresses aspects of the entire spectral trajectory and all forest conditions. Research Question 2 examines the post disturbance recovery spectral trajectory and the forested states that occur after disturbance. Research Question 3 then examines a small portion of the spectral recovery trajectory and likely the initial forest recovery in the stand initiation period 1.4 Dissertation Overview The overarching theme of this dissertation is to provide improved information about boreal forest condition and recovery using remotely sensed imagery, innovative methods, and offering unique analysis. The 2008 opening of the Landsat archive has provided a revolution in terms of data availability and quality (Wulder et al., 2012c) for examining past and present boreal forest conditions. Additionally, the use of the remotely-sensed Landsat imagery archive offers the ability to systematically observe and report on very large and remote areas with less cost and logistic complexity than in situ forest mensuration methods. The structure of this dissertation is as follows: 10   Chapter 1 provides the background information, research needs, research questions and overview. In this chapter the motivating factors around the research are described and then an outline is provided.  Chapter 2 contextualizes the study area and details the common dataset used throughout this dissertation. Canadian boreal forests are discussed in general, and then the three ecozones which form the basis of the research are detailed. The Landsat satellite series, its sensors and legacy are also described in this chapter.  Chapter 3 provides a review of coarse and moderate scale remote sensing studies focused on vegetation trends in the boreal, writ large. Exceptions to the focus on boreal studies occur where the application of remote sensing to monitor forest recovery warrants review.   Chapter 4 addresses the first research question and examines how a time series of remotely-sensed imagery can be used to  estimate a range of boreal forest attributes, in an unmanaged northern Canadian boreal forest. To do so, a Landsat data time series over a single study site in a northern Canadian boreal forest was temporally segmented. A number of spectral time series metrics were calculated to predict a range of forest attributes derived from spatially coincident Light Detection and Ranging (LiDAR) data. Outcomes from this research demonstrate that spectral time series metrics can be used to estimate forest attributes, with their predictive accuracy varying markedly depending on the forest attribute that is estimated. The research also highlights that complex spectral time series metrics are not inherently better than simpler metrics for estimating forest attributes.  11  Chapter 5 addresses the second research question, by: first extending the analysis in Chapter 4 from a single Landsat scene, to a much more comprehensive sample of sites across the Canadian Boreal Shield ecozone, and; second focusing on the spectral recovery to examine differences in forest recovery. First, East and West Boreal Shield post-disturbed recovery trajectories extracted from Landsat time series are compared. Those results show that differences in trajectories between the East and West mirror differences in forest recovery. Specifically, the spectral recovery trajectories detected in the Western Boreal Shield forest are driven by initial conifer establishment, whereas the Eastern forests experience a broadleaf regeneration first which is then succeeded by conifers.  Findings also indicate that selection of the spectral index from which to derive the time series is critical, as it greatly affects which recovery properties are detected and thus the rate of spectral recovery from the time of disturbance to the return of pre-disturbed spectral conditions. A discussion of how spectral and forest recovery is also reviewed in this chapter.  The final research Chapter 6 takes the improved understanding of forest attribute estimation from spectral trajectories from chapter 4 and analysis of the spectral indices behavior from Chapter 5 and examines the rates of spectral recovery over all stand replacing (fire) disturbances within the Taiga and Boreal Shield ecozones. This chapter investigates if spectral forest recovery rates have varied over the past 30 years. To do so, five year post-disturbance recovery trajectories are calculated and used to assess spectral recovery rates at regional scales. Critical results from this final research chapter indicate that spectral recovery rates are not static over time or space across the Canadian Boreal. Furthermore, recovery rates in both the East and West Taiga Shield ecozone divisions show an increasing trend over time.  12   Chapter 7 concludes with a review of the major findings of the research, research innovations, limitations, and possible directions for future research.    13   Study Area & Data Chapter 2:2.1 The Canadian Boreal Zone The extent of Canadian boreal forests (Figure 2.1) are often defined by their high latitudinal location and constraints placed on annual growth by cold temperatures, which result in short growing seasons and frequent long lasting snow cover (Gauthier et al., 2015). Brandt (2009) further defined the boreal zone as dominated by forests, and stretching from Newfoundland across northern Canada into Alaska; northern boundaries are determined by where trees will grow, while southern boundaries are variable depending on local climates, topography and soil conditons (Brandt 2009). Boreal forest stands are typically dominated by coniferous tree species such as white spruce (Picea glauca), black Spruce (Picea mariana), larch (Larix laricina), balsam fir (Abies balsamea), and jack pine (Pinus banksiana), but broadleaf species trembling aspen (Populus tremuloides), balsam poplar (Populus balsamifera), and white birch (Betula papyrifera) can also occur either in pure stands covering large areas or more frequently in mixed conifer stands (Brandt 2009). Those tree species cover a wide range of climatic and topographic conditions with a number of understory species and configurations, and the current distribution is still being relatively novel, as almost all of Canada was covered by ice sheets until 16,000-15,000 years before present (Brandt et al., 2009; Dyke et al., 2003). Each area`s post glacial history and current tree species tolerance for very cold climates combine to explain their present distribution (Sakai et al., 1987, Brandt 2009). Thus, boreal tree species are thought to be wide ranging generalists with wide environmental tolerances.  14    Figure 2.1 The boreal zone of Canada (Brandt 2009) highlighted against a false color imagery composite   Boreal forest changes are largely driven by natural disturbances which vary in spatial extent, intensity, and through time (Volney et al., 2005, Brassard et al., 2008). Fire, insect outbreaks, harvests, and disease all cause disturbances in the boreal forest (Volney et al., 2005); however, fire is the most significant disturbance with respect to severity and extent (Bond-Lamberty et al., 2007; Weber et al., 1997). Fire intensity and behavior largely determine the subsequent recovery and regrowth of vegetation (Wofsy et al., 2002; Harmon et al., 1990). Insect infestations are also 70°W80°W90°W100°W110°W120°W60°N50°N40°NCanadianBorealZoneUnitedCanadaStatesCanadian Boreal ZoneProvinces15  a key disturbance agent in boreal forests that can result in widespread impacts to annual forest stand growth, and in some cases may result in mortality (Flemming et al., 2000).  Harvest disturbances occur most often in the form of clear cuts, and more often in the southern managed boreal forest, with replanting occurring in some cases. However widespread or severe the impact, disturbances are essential for boreal forest ecosystems as they are critical for maintaining the patchworks of species composition, ages, and stand structures that characterize this forest. Thus, Canadian boreal forests involve complex interaction between disturbance types, tree species, climatic conditions, and topography combining to form the present non-static state of the boreal. 2.2 Ecozone Study Areas The boreal forest mosaic can be sub-divided along similar biotic and abiotic conditions, such that homogenous areas of climate, topography, species, and ecosystem processes are grouped together at multiple scales. Ecozones reside at a high level within the spatial framework, and are “area[s] where organisms and physical environment endure as a system” (Wiken 1986). This spatial framework is useful when examining large boreal forest areas, as they delineate similar areas where forest processes such as disturbance and recovery cycles are driven by similar agents and time scales across the landscape.    16    Figure 2.2 Map of Canadian ecozones with management areas displayed; highlighted zones indicate where research was performed Three ecozones (Figure 2.2) comprise the focus of this dissertation; the Boreal Cordillera, Boreal Shield, and the Taiga Shield. In Chapter 4, a subset of the Boreal Cordillera is used to examine the applicability of temporal trajectory metrics for estimating forest attributes. The research in Chapter 5 occurs across the Boreal Shield ecozone, where spectral forest recovery is compared across an East/West division. In Chapter 6, both the Taiga Shield and Boreal Shield ecozones are 20°W40°W60°W80°W120°W140°W160°W 100°W60°N50°N40°NEcozonesUnmanaged Boreal ForestManaged Boreal ForestProvincesBorealCordilleraTaiga  Shield    WestBoreal Shield   WestBoreal ShieldEastTaiga ShieldEast17  used to investigate spectral forests rates over time. A summary of the biotic and abiotic factors in the three highlighted zones is provided in Table 2.1. Table 2.1: Summary of biotic and abiotic factors of Boreal Cordillera, Boral Shield, and Taiga Shield ecozones. Boreal Cordillera The Boreal Cordillera extends from northern British Columbia to southern Yukon and is considered the midsection Canada’s Western Cordillera. Mountain ranges with high peaks and extensive plateaus permeate the landscape, separated by wide valley lowlands. Climatically the area experiences cold winters with short warm summers. Large rivers wind through the landscape, enabling hydroelectricity development and forestry in additional to the extraction of mineral resources.  Open and closed canopy forests white and black spruce, alpine fir (Abies lasiocarpa), and lodgepole pine (Pinus contorta) are common, though terrain and elevation moderate species and density with valleys, steep slopes and plateaus all containing differing assemblages (Wiken 1986; Yukon Ecoregions Working Group 2004).    Mean Temperature Summer 10° C Winter -18° C Mean Annual Precipitation 300 – 1500 mm, elevation dependent    140°0'0"W 40°0'0"W60°0'0"W80°0'0"W100°0'0"W100°0'0"W120°0'0"W160°0'0"W60°0'0"N60°0'0"N50°0'0"N50°0'0"N40°0'0"N40°0'0"NEcozonesBoreal ForestProvinces18  Boreal Shield The largest of all ecozones where much of the land remains inaccessible and in a wilderness condition. This ecozone forms a horseshoe-shaped area between northern Saskatchewan and Newfoundland. Vast expanses of trees, water, and bedrock characterize the Boreal Shield landscape with over 80% of the characteristic broad rolling uplands forested and interspersed with wetlands. Its continental climate comprises long, cold winters and short, warm summers. Eleven and half percent of Canada’s population occupies the Boreal Shield, with human settlement developing around natural resource bases. Forests are dominated by closed stands of white and black spruce that pervade the landscape with more broadleaf species such as white birch, trembling aspen, and balsam poplar mixing in the south (Wiken 1986; Ecological Stratification Working Group 1996).    Mean Temperature Summer 13 C Winter -10 C Mean Annual Precipitation 400 - 1000 mm, higher on coastal areas    140°0'0"W 40°0'0"W60°0'0"W80°0'0"W100°0'0"W100°0'0"W120°0'0"W160°0'0"W60°0'0"N60°0'0"N50°0'0"N50°0'0"N40°0'0"N40°0'0"NEcozonesBoreal ForestProvinces19  Taiga Shield Divided in half by Hudson Bay, the Taiga Shield is a patchwork of numerous lakes, wetlands and open forests. Bounded in the east by the Atlantic coast of Labrador through Quebec to Hudson Bay, it continues on the opposite shore of Hudson Bay through northern Manitoba, Saskatchewan, an eastern portion of the Northwest Territories, and the southernmost reaches of Nunavut. While the latitudinal center contains stunted trees, the northern edge is the latitudinal limit of tree growth. Winters are long and cold with short, cool summers. The zone’s rolling topography is a mosaic of treed uplands and wetlands. Forest are dominated by stunted black spruce and jack pine stand, with mixed wood stands of balsam fir, trembling aspen, and white birch commonly occurring (Wiken 1986, Ecological Stratification Working Group 1996).     Mean Temperature Summer 8.5° C Winter -18° C Annual Mean Precipitation 200 – 800 mm, higher in coastal areas  2.3 Landsat Data Remotely-sensed datasets can provide detailed spatially explicit data on Canadian boreal forests (Wulder et al., 2004). The unique characteristics of Landsat data result in a dataset that is well suited for many applications, and particularly forest monitoring  and reporting, which have become a central of focus (Wulder et al., 2012a; Wulder et al., 2008a; Hall et al., 2006; Wulder et al., 2004b; Skakun et al., 2003).  140°0'0"W 40°0'0"W60°0'0"W80°0'0"W100°0'0"W100°0'0"W120°0'0"W160°0'0"W60°0'0"N60°0'0"N50°0'0"N50°0'0"N40°0'0"N40°0'0"NEcozonesBoreal ForestProvinces20  Landsat data is frequently used for monitoring applications due to its cost free availability, synoptic capture of landscapes, and frequent sixteen-day revisit period . While these characteristics are important, the consistent radiometry, distribution of analysis ready data, and the temporally deep imagery archive add to the value of Landsat data (Wulder et al., 2012c). The historic imagery archive is an important feature of the data, allowing detection of changes in reflected radiation corresponding to changes in ground conditions; i.e. multiple snapshots of landscapes over a forty year time period are captured in the imagery archive (White et al., 2013;).   21  The history of the Landsat series of sensors and missions are summarized in Figure 2.3 and Table 2.2, and displays the length of the over 40 year archive of imagery available to users. The Multi-Spectral Scanner (MSS) era started with Landsat 1 in 1972. The MSS sensor collected data across four bands of the electromagnetic spectrum, (green 0.5-0.6 µm, red 0.6-0.7µm, and near infrared 0.7-0.8 and 0.8-1.1µm) at 59 by 79m spatial resolution (distributed now as 60 by 60m), and was also onboard Landsat 2 and Landsat 3. Although the imagery captured by the MSS sensor represents some of the earliest snapshots of an area, these data are often not included in studies due to their coarser spatial resolution, and lack of shortwave infrared bands, compared to later Landsat sensors.   Figure 2.3: Timeline of Landsat satellites and sensors illustrating the temporal depth of the archive. The launch of Landsat 4 in 1982 ushered in a new era of land remote sensing with the Thematic Mapper (TM). The TM sensor offered a number of improvements over the MSS sensor, which 22  was also onboard Landsat 4, and included increased spatial resolution (30m) and data collection from an additional portion of the electromagnetic spectrum: the shortwave infrared. Thermal infrared (longwave infrared) data was also collected by the TM sensor, which was also collected by Landsat 3 for a very short period of time before that channel failed. As a follow on mission, Landsat 5 launched in 1984 carried both MSS and TM sensors, with TM data collected operationally until 2011. Landsat 5’s MSS sensor collected data from 1984 until 1999, and then again in limited bursts between 2012 and 2013. The legacy of Landsat 5 is that it ensured imagery capture was continuous despite Landsat 6’s failure to reach orbit in 1993. Landsat 7 was launched in 1999 carrying the Enhanced Thematic Mapper Plus (ETM+) instrument, continuing data collection in the same bands as TM sensor but dividing the thermal spectrum into two separate bands acquired at an increased spatial resolution of 60m. The Enhanced Thematic Mapper Plus instrument also added panchromatic information, which had the highest resolution (15m) of any Landsat data (Loveland & Dwyer 2012). However, in 2003 the sensor suffered an issue where 22% of the data within an image is not collected, with areas near the edges most affected.  Continuing in the same tradition as Landsat 7, Landsat 8 is the latest addition to the program, carrying a reduced-spatial resolution (100m) Thermal Infrared Sensor and the Operational Land Imager (OLI). The OLI sensor continues to collect information in the same bands as Landsats 5 and 7 and at the same spatial resolution, which importantly maintains consistency with previous Landsat sensors. In addition, two new spectral bands were added to the sensors capabilities: an ultra-blue band for coastal applications and an infrared band for cirrus cloud detection. Using much of the same technology as Landsat 8, Landsat 9’s launch is slated for 2023 and will 23  continue imaging in the same historic spectral bands and spatial resolutions as previous sensors, possibly pushing the archive past the 50 year mark.  The success story of the Landsat program has also been beset with troubles and failures. After Landsat 5 was launched, the program underwent commercialization, which had disastrous effects for its users. As the program was transferred to a commercial entity, data costs skyrocketed and data redistribution was curtailed; many previous users were locked out of using Landsat data for their applications. The Landsat program was retaken by the US federal government in 1992, but the still commercially owned Landsat 6 failed to reach orbit in 1993. The lack of a new satellite in orbit created the fear of a data-gap, however Landsat 5 stayed operational until 2012, even through the launch of Landsat 7 in 1999. Landsat 7 suffered its own imaging error 2003, which underscored the need of a new satellite in orbit, which required fourteen years to fund, build, launch and then make operational. Thus, although the Landsat program is largely a success story, it is not without its problems throughout its timeline.  Table 2.2: Spatial, spectral  and temporal characteristics of Landsat sensors. Aero: Aerosol/Coastal; NIR: Near Infrared; PAN: Panchromatic; SWIR: Shortwave Infrared; TIR: Thermal Infrared. Sensor Satellites Years Active Spatial Resolution (m) Band Spectral Resolution (µm) Temporal Resolution Multi-Spectral Scanner Landsat 1 – Landsat 5 1972 – 1999, 2012-2013 60 Green 0.5-0.6 18  days Red 0.6-0.7 NIR1 0.7-0.8 NIR2 0.8-1.1 Thematic Mapper Landsat 4, Landsat 5 1984 -2011 30 Blue 0.45–0.52 16 days Green 0.52–0.60 Red 0.63–0.69 NIR 0.76–0.90 SWIR1 1.55–1.75 SWIR2 2.08–2.35 120 TIR2 10.40–12.50 24  Sensor Satellites Years Active Spatial Resolution (m) Band Spectral Resolution (µm) Temporal Resolution Enhanced Thematic Mapper Landsat 6 1993,  Did not achieve orbit 30 Blue 0.45–0.52 16 days Green 0.52–0.60 Red 0.63–0.69 NIR 0.76–0.90 SWIR1 1.55–1.75 SWIR2 2.08–2.35 120 TIR2 10.40–12.50 15 Pan 0.52–0.90 Enhanced Thematic Mapper + Landsat 7 1999 - present 30 Blue 0.45–0.52 16 days Green 0.52–0.60 Red 0.63–0.69 NIR 0.76–0.90 SWIR1 1.55–1.75 SWIR2 2.08–2.35 120 TIR2 10.40–12.50 15 Pan 0.52–0.90 Operational Land Imager & Thermal Infrared Sensor Landsat 8 2013 - present 30 Aero 0.43–0.45 16 days Blue 0.45–0.51 Green 0.53–0.59 Red 0.64–0.67 NIR 0.85–0.88 Cirrus 1.36–1.38 SWIR1 1.57–1.65 SWIR2 2.11–2.29 100 TIR1 10.60–11.19 120 TIR2 11.50–12.51 15 PAN 0.50–0.68  Landsat sensors have been collecting imagery over the globe since 1972, continually adding to the archive of imagery. After the imagery in the archive was made freely available, the use of Landsat data for scientific inquiry became more widespread, as the number of images distributed annually increased dramatically after the archive was opened. In their research on the effect of opening archive, Wulder et al., (2012c) stated: 25  “With nearly 40 years of continuous observation… the Landsat program has benefited from insightful technical specification, robust engineering, and the necessary infrastructure for data archive and dissemination. … [The spectral and spatial] resolutions have proven of broad utility and have remained largely stable over the life of the program. ... The new [open access] data policy is revolutionizing the use of Landsat data, spurring the creation of robust standard products and new science and applications approaches, …[and] also increased international collaboration to meet the Earth observing needs of the 21st century.”   26   Canadian Boreal Change as Assessed with Remotely Sensed Chapter 3:Datasets  3.1 Coarse Spatial Resolution Monitoring The large expanse of the boreal has been observed with a variety of remote sensors and analysis techniques, however these observations and models have often undertaken with coarse resolution data (≥ 1km2) due to their long time series and analysis ready data. Research is often performed at the continental, pan-boreal, or global scales, to detect Normalized Difference Vegetation Index (NDVI) trends over time that are representative of changes in photosynthetic capacity, or vegetation vigor and health. Many researchers have found negative NDVI trends, suggestive of less photosynthetically active vegetation, or decreasing quantity, or vigor (Beck et al., 2011; Goetz et al., 2005; Bunn et al., 2009; Zhou et al., 2001) as shown in Table 3.1. However, others using similar datasets have found positive trends, suggesting increased vegetation photosynthetic activity and vigor (Myneni et al., 1997; Nemani et al., 2003; Slayback et al., 2003; Olthof et al., 2007).   Research has acknowledged there is more complexity to these negative and positive NDVI trends demonstrating there can be multiple and divergent trends simultaneously occurring throughout the Canadian boreal, both spatially and temporally. For example, both increasing and decreasing trends were found by Pouliot et al., (2009) using NDVI time series. In that research they noted that southeastern Canadian boreal forests showed a decreasing trend, while the remainder showed an increasing forest trend. Furthermore, Angert et al., (2007) found that the Canadian boreal forests first increased in NDVI during one time period (1985-1991) and then decreased throughout a second period (1994-2002). These conflicting trends suggest that boreal 27  forests undergo complex changes over time and they may not be best observed with such coarse spatial scale data. Recent research has now shown that the preprocessing steps and inter-sensor calibration of course grained data can bias the resultant time series towards one direction or another (Alcaraz-Segura et al., 2010). Moreover, most of those coarse scale studies improperly handle the inherent spectral trends associated with boreal disturbance and recovery. By not analyzing them separately, could lead to erroneous trends to be detected. Lastly, coarse datasets likely cannot capture the fine scale mosaicked nature that defines boreal landscapes.28  Table 3.1: Coarse spatial resolution remote sensing based vegetation studies summarized and pertinent Canadian forest results highlighted   Extent Time Period Time Step/ Resolution/ Sensor Analysis Results Forest Trend Source Pan Boreal 1981 – 1991 Daily, 8 km, AVHRR Regression of NDVI Time Series Photosynthetic activity increased  Increasing Myneni et al., 1997 Northern Hemisphere 1981 -1999 Annual, 8km, AVHRR Linear regression of NDVI time series Canadian boreal forests generally decrease in NDVI  Decreasing Zhou et al., 2001 Worldwide 1982 – 1999 Annual, 8km, AVHRR Linear regression of NPP time series derived from NDVI North American Boreal forests NPP increased Increasing Nemani et al., 2003 Pan-Boreal 1982-1999 Annual, 8km, AVHRR Non-linear regression of NDVI time series Forest areas increase in photosynthetic activity Increasing Slayback et al., 2003 North American Boreal 1981-2003 15-Day, 8 km, AVHRR Auto-Regression of NDVI Time Series Forest areas decrease in photosynthetic activity Decreasing Goetz et al., 2005 Northern Hemisphere 1982 – 2002 Monthly, 8km, AVHRR Linear regression of spring NPP time series derived from NDVI North American Boreal forests increase then decrease in NDVI over 2 periods Both increasing and decreasing Angert et al., 2007 Canadian Boreal 1998 – 2004 Daily, 1km,  SPOT-VGT Linear Regression of NDVI Time Series Forests areas increase in NDVI Increasing Olthof et al., 2007 Pan-Boreal 1982 – 2003 Annual,    64 km,  AVHRR NDVI derived photosynthetic capacity analyzed with trend detection tests North American Boreal forests decrease in photosynthetic capacity Decreasing Bunn et al., 2009 29    Extent Time Period Time Step/ Resolution/ Sensor Analysis Results Forest Trend Source Canadian Boreal 1985 – 2006 10-day, 1km, AVHRR Thiel-Sen Regression of NDVI Time Series Most Forest areas increase in NDVI, SE decreasing Both increasing and decreasing Pouliot et al., 2009 North American Boreal 1982 – 2008 15-day, 8km , AVHRR Linear Regression of NDVI Time Series Forest areas decrease in photosynthetic activity Decreasing Beck et al., 2011 30  3.2 Landsat Based Medium Spatial Resolution Monitoring Despite issues surrounding the use of coarse spatial resolution datasets to monitor boreal forests, they were still representative of the best available data at the time of analysis. In contrast to the studies performed at coarser spatial resolution in boreal forests, medium spatial resolution (10 m2 - 1km2) based research has historically focused on using fewer observations throughout their time series, but with inherently greater spatial detail (Table 3.2). For example, Jin  et al., (2005) used multitemporal Landsat data and image differencing techniques to determine that coniferous forest clear-cuts  are still detectable five years after they occurred. Using an image pair with 28 years lapsing between them, a study in northern Quebec and Northwest Territories, Canada found that boreal forests did not expand beyond their original extent observed in the first image (Masek et al., 2001). Forest disturbance rates across North America were determined by Masek et al., (2012) using decadal imagery pairs across a large sample of sites distributed across North America. In that study they used an image differencing technique and found that areal disturbance rates of 2-3% were common for many forests across the United States and Canada. The information about forests that these studies generated is crucial and represents a step forward from using single date of imagery analysis techniques towards observing the boreal over longer time periods. However research using imagery pairs are prone to bias introduced from the selection of anniversary dates and the length of time between images, both of which can critically affect the results (Gillanders et al., 2008; Coppin et al., 2008).   A new era in multitemporal Landsat research began when the Landsat imagery archive was opened in 2008, making imagery freely available, easy to obtain, and distributed in an analysis ready state (Wulder et al., 2012c; Woodcock et al., 2008). As a result newer applications are not 31  limited by imagery acquisition costs and preprocessing burdens, thus many novel methods are available to exploit the temporally deep and analysis ready Landsat archive. The imagery present in the archive represents the condition of landscapes over time, providing multiple views into the condition of landscapes at multiple points in time. By utilizing the archive, past and current conditions of landscapes can be analyzed using new time series methods. The resultant time series data can then be used to inform on trends at local, regional and national scales.  New Landsat-designed, automated time series analysis algorithms have allowed large spatial areas to be analyzed, importantly dividing a spectral time series into periods of no change, disturbance, and recovery. These novel temporal segmentation methods have allowed the quantification of disturbance as well as recovery, and typically at annual time steps with a number of spectral vegetation indices (Table 3.3). For example, some algorithms composite multiple images within a single growing season into a single annual observation at the pixel level, and then detect change as departures from observations taken in previous years. One temporal segmentation method that detects change as departures from stable observations is the Vegetation Change Tracker (VCT) algorithm, which examines annual spectral values in reference to statistical properties of nearby forested pixels to divide a spectral time series (Huang et al., 2010). Similarly, Kennedy et al., (2010) developed the LandTrendr algorithm, which emphasized the ability to extract information about disturbance and recovery. Lastly and a departure from annual observations, the Continuous Monitoring of Forest Disturbance Algorithm (CMFDA) can ingest and segment twice monthly observations to determine change between and within years, but is only viable where many observations occur regularly throughout a year (i.e. the continental U.S.) (Zhu et al., 2012a).  32  These three algorithms (VCT, CMFDA, and LandTrendr) have been applied to a number of forested ecosystems to detect disturbance and recovery over time (Table 3.2). Examining disturbance and recovery trends in the Pacific Northwest United States after a forest harvest policy change, Kennedy et al., (2012) showed that restrictions on timber harvest did reduce the number of disturbances on public land, but private land saw a corresponding increase in clear cutting. Equally important are the recovery trends found by that study, which indicated that moist regions produced quicker forest recovery than drier areas (Kennedy et al., 2012). Forest recovery as observed by a Landsat time series proved an important variable when estimating stand biomass in eastern Oregon (Pflugmacher et al., 2012). Similarly Powell et al., (2010) found at sites Minnesota and New Mexico, USA that biomass estimates after a disturbance event were more accurate when spectral forest recovery information was used in the prediction . Thus, temporally-segmented Landsat time series can produce data that informs on vegetation succession and structural development.  In addition to examining change over time, open access to data makes cross-jurisdiction large area examination possible, which can reveal broad forest trends that are not evident at the spatial scale of single jurisdictions. Kennedy et al., (2012) examined Landsat imagery over three US states and was able to cross multiple state and ownership boundaries to arrive at their disturbance and recovery findings, which showed different harvesting patterns resulting from policy decisions.  Over a larger area and encompassing 20 separate jurisdictions, the entirety of European Russia was used as a study area to detect increasing forest cover and recovery in abandoned agricultural fields (Potapov et al., 2015). Griffiths et al., (2014) examined Landsat time series across Eastern Europe and noted that forest recovery differed by national boundaries 33  and over time. Using a similar study area, Turubanova et al., (2015) found that 85% of disturbed forests showed signs of recovery, which may not have been evident if national harvest and disturbance summary statistics were examined in isolation.  Due to the frequent disturbances that spur the subsequent regrowth, Canadian boreal forests are an ideal place to examine forest recovery using Landsat time series. The recovery of a forest disturbed by fire in the Yukon was examined by Epting et al., (2005) who showed that spectral indices plateau after 10-15 years of recovery. Interestingly, spectral indices based on the shortwave infrared bands were more important for tracking post disturbance recovery since they report more on forest structure than canopy cover (Schroder et al., 2011). Similarly, Fraser et al., (2011) related trends in the Tasseled Cap Greenness and Wetness spectral indices to post disturbance forest recovery in northern Manitoba. Extending that concept over a larger forested area in the Northwest Territories, the same authors were able to show that the spectral indices continued to increase until a change in forest successive state occurred (Fraser et al., 2014). A more detailed examination of forest recovery conducted over multiple time periods revealed that spectral recovery could vary over time, with some periods showing faster or slower regrowth rates (Pouliot et al., 2009). Monitoring forest recovery through spectral time series can also generate valuable information that can be used for land cover change classification. Data on spectral forest recovery was used to correctly classify harvest and fire disturbances in central Saskatchewan (Schroeder et al., 2011). Also working in the forested portion of Saskatchewan, Hermosilla et al., (2015a, 2015b) demonstrated that multiple disturbance agents, with varying severities, could be accurately identified using Landsat time series. In that work, they also showed the utility of quantifying 34  spectral forest recovery by indicating that post disturbance trends are important to understand current forest conditions. For example, Landsat time series analysis demonstrated forest recovery became more spatially dispersed over time in Albertan forests (Gómez et al., 2015).   Spectral time series and temporal segmentation methods benefit from having more data to draw upon to form gap free annual composites, obtain multiple observations within a growing season, or to extend a spectral time series beyond the limitations of one sensor. Although Landsat data has the long term historic and open access imagery archive, other sensors can complement spectral time series. For example, some SPOT sensors have been observing the Earth since 1986 with some compatibility between spectral bands and spatial resolution. Most promising is the recent launch of the Sentinel 2 satellites, as their spatial and spectral resolutions offer a level of data interchangeability previously unavailable. With the current Landsat and Sentinel satellites in orbit, it is possible to observe a region multiple times every seven days.  In addition, the imagery collected by Sentinel 2 satellites will be distributed freely, in a similar manner and policy as that of Landsat data.  In conclusion, coarse spatial resolution datasets have conventionally been used to monitor boreal forests; however, there remains disagreement on the vegetation trends observed by those studies. Most studies have confirmed that coarse spatial resolution datasets are poorly designed to examine the fine scaled change that is typical of the boreal forest ecosystem, suggesting that a finer spatial scale dataset would be more useful. Landsat data are ideal for this monitoring task and have a proven record of providing useful forest information, which is often generated from spectral time series.  35  Table 3.2: Landsat based time series analysis of forests summarized with recovery aspects highlighted Extent Time Period & Step Analysis Results Forest Recovery Aspect Source 2-37,500  km2 areas in Northern NWT and QC 1972 – 2000, 28 year interval NDVI change detection of tree lines Forest expansion did not occur Boreal forests exhibit a stable distribution through the Landsat imagery record Masek 2001 170 km2 boreal forest fire on AK/YK border 1986 – 2002, 3 year interval Post fire recovery displayed by NBR and NDVI for burned sites Post fire recovery dependent on fire severity Spectral indices peak 10-15 years post fire; trees dominantly self-replace Epting et al., 2005 5000  km2 in ME 1988-1993 & 2000-2002,  Multiple intervals Image differencing between two dates of TCW and NDMI Clear cuts could be delineated with 1,2 and 5  intervals between images None Jin et al., 2005 23-35,000 km2  regions in US 1990 – 2000, 10 year interval SVI based differencing with thresholds Disturbance rates of 2-3% are common across the US and CA Forest recovery hard to detect without an annual time step Masek et al., 2008 2 – 34,225 km2 areas in northern QC and NWT/YK border 1986 - 2006, 6 year interval Linear regression of NDVI time series Positive NDVI trends occur Post disturbance NDVI recovery rates can change over time Pouliot et al., 2009 Multiple U.S. locations 1984-2012, 2 year interval Detection of disturbance and linear segmentation of spectral time series VCT segmentation algorithm displayed and accurate over a number of forest environments Forest recovery tied to similarity of neighboring undisturbed forests Huang et al., 2010 4-35,000 km2 in WA and OR 1984 – 2009, 1 year interval Linear segmentation of time series for disturbance and recovery detection LandTrendR segmentation algorithm displayed  Spectral recovery tied to canopy cover thresholds Kennedy et al., 2010 36  Extent Time Period & Step Analysis Results Forest Recovery Aspect Source 2-35,000 km2 regions in MN & NM 1985 – 2006,2 year interval   LandTrendR based temporal segmentation of a time series then related to biomass change over time Biomass trends were fit to spectral trends  Better tracking of forest regrowth leads to better post disturbance biomass estimates Powell et al., 2010 34,225 km2  total area in central SK 1986 – 2008, 1 year interval Temporal segmentation of spectral time series, post disturbance NDVI and NBR tracked Harvest and Fire disturbances types can be determined from recovery trajectories SWIR based indexes outperform NIR based indexes for recovery trajectory tracking Schroeder et al., 2011 11,476 km2 in northern MT 194-2009,  1 year interval Non-parametric trend detection of Tasseled Cap Indices and NDVI time series Increasing Greenness, Wetness and NDVI trends, decreasing Brightness trend Increasing Greenness and Wetness trends related to post disturbance regeneration Fraser et al., 2011 13,818 km2 in AB 1973 – 2008,  3 year interval Lagrange interpolation of Tasseled Cap Angle time series Areal size of regions undergoing change increased near the end of the time series Forest recovery was more dispersed at the end of the time series  Gómez et al., 2011 229900  km2 area in WA, OR and CA 1985 – 2008, 1 year interval LandTrendR algorithm applied to a large area Disturbance and recovery enumerated for federal, state, and private landholders Recovery Indicator introduced, used to show more moist regions recover more quickly Kennedy et al., 2012 830 km2 in OR 1972 – 2010, 1 year interval LandTrendR based temporal segmentation of a time series, metrics used to predict forest structure Time series metrics better predict forest structure than single date Integration of spectral recovery information better helped predict biomass than other metrics Pflugmacher et al., 2012 3,600 km2 in GA 2001-2004, 16 day Sine wave fitting of bi-weekly reflectance for disturbance detection CMFDA segmentation algorithm displayed, showing accurate mapping of disturbances None Zhu et al., 2012a 37  Extent Time Period & Step Analysis Results Forest Recovery Aspect Source 300,000 km2 area in NWT 1985-2011,  1 year interval Non-parametric trend detection of Brightness, Greenness, and Wetness time series Trend composites are useful to identify change agents Over time post-fire Brightness decreases, Greenness and Wetness increase, until a next successive state is reached Fraser et al., 2014 390,000  km2 in Eastern Europe 1985 -2010, 5 year interval Disturbance detection based on post classification comparison, spectral recovery tracked  20% of forests experienced a stand replacing disturbance, a small net increase in forest area occurred Forest recovery differed by national boundaries Griffiths et al., 2014 31 781 km2  in Mongolia 1995 – 2014, 1 year interval Post fire recovery in Larch dominated forests monitored with two SVIS Full recovery of canopy characteristics in 30-45 years, correlated with burn severity A first index used to monitor canopy cover, a second used to monitor structure  Chu et al., 2015 378,818 km2 forested area of SK 1984- 2012,  1 year interval Temporal segmentation of NBR time series  Annul change detected, agents reliably attributed using trajectory metrics Post disturbance spectral recovery rates quantized  for a large area Hermosilla et al., 2015a, 2015b 1,138,825  km2 , in Eastern Europe 1985 – 2012, 1 year interval Multi-temporal metrics used to map forest cover change over time Forest cover increased by 4.7% over the region, though average annual loss was .41% 85% of disturbed areas show recovery Turubanova et al., 2015 Canada and Alaska 1984 – 2012, 1 year interval Linear regression of peak NDVI time series Greening trends for northeastern Canada and  browning trends for northern AB and SK NDVI typically rises post fire, not evident in coarser data Ju et al., (2016) Table 3.3: Frequently used spectral vegetation index calculation equations. Band references are made in terms of the portion of the electromagnetic spectrum that is covered, and in reference to Table 2.2   38  Table 3.4: Frequently used spectral vegetation index calculation equations. Band references are made in terms of the portion of the electromagnetic spectrum that is covered, and in reference to Table 2.2    Spectral Index Acronymn Equation Normalized Difference Vegetation Index NDVI  =(NIR − RED)(NIR + RED) Normalized Burn Ratio NBR  =(NIR − SWIR2)(NIR + SWIR2) Tasseled Cap Brightness TCB  = (	 ∗ 	0.3037) + ( ∗ 	0.2793) + ( ∗ 	0.4743) + ( ∗ 	0.5585)+ (1 ∗ 0.5082	) + (2 ∗ 	0.1863) Tasseled Cap Greenness TCG  = (	 ∗ 	−0.2484) + ( ∗ −0.2435	) + ( ∗ 	−0.5436) + (∗ 	0.7243) + (1 ∗ 	0.0840) + (2 ∗ 	−0.1800) Tasseled Cap Wetness TCW  = (	 ∗ 	0.1509) + ( ∗ 	0.3279) + ( ∗ 0.3279	) + (∗ 	0.3406) + (1 ∗ 	−0.7112) + (2 ∗ 	−0.4572) 39   Spectral Time Series Metrics for Forest Structure Prediction Chapter 4:4.1 Introduction  Boreal forests comprise 27% of the world’s forest cover (FAO 2001), and in northern latitudes represent areas where human intervention has been minimal and ecosystem processes generally remain undisturbed (FAO 2006). With the exception of Scandinavia, much of the global boreal forest does not have forest inventory level data or is not regularly monitored (Wulder et al., 2007). While most of the global boreal forests are located in Russia, Canada contains a substantial portion of the North American boreal (Brandt 2009). In Canada, provincial and territorial resource management agencies hold forest stewardship and management responsibilities, which include implementation of inventory and monitoring programs. Boreal stands that are located in southern Canadian boreal forest are typically subject to inventory on a periodic basis (often on a 10 year cycle) to support forest management and planning goals to meet a range of sustainable forest management objectives. These managed forests are generally where most of the human population resides, road access is often greater, and vegetative productivity is higher. In contrast, Canada’s northern boreal forests are subject to more limited and spatially focused management activities and as a consequence, inventory and monitoring programs receive less emphasis. As a result, the condition and state of the unmanaged northern forests are not as well characterized and understood when compared to the southern managed forests (Andrew et al., 2012).  Managed and unmanaged areas of Canadian boreal forests are subject to a range of disturbances. Fire, harvests, insect outbreaks, and disease are all agents of boreal forest disturbance (Volney & 40  Hirsch 2005), each with variability in intensity and scale through time and by region (Brassard & Chen 2008). Of these, fire is the most notable in severity and extent over boreal forests (Bond-Lamberty et al., 2007). Therefore, in boreal ecosystems fire intensity is largely responsible for determining the subsequent recovery and regrowth of vegetation (Johnstone et al., 2006a, Johnstone et al., 2006b, Johnstone et al., 2005, Payette et al., 2002). Although the location and extent of fires is tracked by a number of agencies (de Groot et al., 2007; Stocks et al., 2002), the subsequent recovery from these disturbances is only monitored through the re-occurring inventory practices in southern Canadian boreal forests, but not in northern Canadian boreal forests.  Despite their remoteness, information on northern Canadian boreal forests is still necessary, both to meet regional information needs and to support national and international reporting requirements (Wulder et al., 2004). In response to this need and in partnership with the provinces and territories, Canada has implemented a National Forest Inventory (NFI). The NFI is a sample based system based upon a national network of 2 x 2 km aerial photo plots on a 20 km grid. The photo plots are calibrated against, and supported by, an integrated network of ground plots (Gillis et al., 2005). For some larger Canadian jurisdictions, the NFI can provide a sufficient sample to provide reliable statistical estimation of a variety of important forest attributes (e.g. species, age, crown closure, disturbance agents and year, and other detailed information) used for monitoring and reporting purposes across the managed boreal in ten year intervals.  To augment this paucity of spatial information in the northern boreal, remote sensing techniques offer potential solutions to the remoteness, difficulty of access, and high costs of establishing and then maintaining ground plots (Wulder et al., 2007). This type of remote sensing based 41  characterization includes: (1) the capture of changes in forest stand condition over time and; (2) the prediction of aboveground biomass and its components for a range of applications, including the replication of some forest inventory attributes to improve understanding of and quantification of national carbon dynamics.  There is a long history of estimating aboveground biomass in forests using single date Landsat imagery. Fazakas et al., (1999) utilized Landsat Thematic Mapper bands 1 -5, and 7 to estimate aboveground biomass in Sweden. Hall et al., (2006) used Landsat Thematic Mapper bands 3,4,5, and 7 to derive estimates of canopy cover and stand height, which were subsequently used to predict biomass in Alberta, Canada. Labrecque et al., (2006) utilized Landsat bands, numerous vegetation indices, band ratios, and Tasseled Cap components to estimate aboveground biomass in Newfoundland, Canada, demonstrating that Tasseled Cap components were the most correlated all of the variables tested. Due to the use of many different spectral variables, there is no broad consensus on which bands or indices should be used to estimate aboveground biomass and other boreal forest attributes from Landsat data.  Previous studies have shown that temporal trajectory methods aid in forest attribute estimation over conventional methods which rely on a single date of optical imagery (e.g.,  Pflugmacher et al., 2012; Powell et al., 2010). Powell et al., (2010) used a large number of spectral indices and environmental variables to estimate aboveground forest biomass in Arizona and Minnesota over a 20+ year time series. Results showed that aboveground biomass change estimation was improved when using temporal trajectories as opposed to using only single date imagery. Pflugmacher et al., (2012) found that the use of temporal trajectory metrics produced higher 42  correlations than variables created from single date imagery. They found that temporal trajectory metrics could reliably estimate relative aboveground biomass in Douglas-fir forests in Oregon.  The primary objective of this chapter is to measure the capacity of temporal trajectory metrics for estimating selected forest attributes in a northern Canadian boreal forest. A secondary goal was to investigate the importance of different types of spectral time series metrics when estimating these forest attributes. To do so we: (1) acquired, normalized, and composited a dense stack of Landsat images of a northern boreal study area; (2) subjected the composited Landsat stack to temporal segmentation using the Landtrendr algorithm; (3) calculated a number of simple and complex temporal trajectory metrics and constructed variable sets based on each Tasseled Cap components; and 4) created random forests models of LiDAR derived above ground biomass with the variable sets. To develop an improved understanding of which temporal trajectory metrics had the most predictive power we then compared R2, RMSE% and importance values from the developed models for aboveground biomass. By selecting and then combining the most important variables from previous models, a new and final variable group was created which was then used to estimate aboveground biomass, stem biomass, stem volume, basal area, and Lorey’s mean height. We review the R2, RMSE% and mean importance values from these results and lastly we conclude with a discussion on the ability of temporal trajectory metrics and the importance of types of spectral time series metrics for forest attribute estimation.  43  4.2 Methods 4.2.1 Study Area The 29,067 km2 rectangular study area is defined by a Landsat image footprint (path 54, row 18 of Worldwide Reference System II) in northeast British Columbia and southeast Yukon Territory which is entirely within the Boreal Cordillera ecozone and substantially divided by portions of four ecoregions (Figure 4.1). The Liard Basin ecoregion covers the central 77% of the study area, and is characterized by low rolling hills (elevation 600-1200m) separated by wide flat areas. Rainfall is evenly distributed throughout the year totaling approximately 500mm, with a drier period between February and May that partially coincides with the cold season. Winter and summer temperatures average -25 °C and 12 °C, respectively. Surrounding the Liard Basin are three ecoregions with mountainous terrain, high plains, and plateaus that can reach elevations up to 2000 meters (Ecological Stratification Working Group 1996, Yukon Ecoregions Working Group, 2004). Dominant vegetation types in these three ecoregions are relatively even aged boreal Black Spruce (Picea mariana) and White Spruce (Picea glauca) forests, with Paper Birch (Betula papyrifera) and Lodgepole Pine (Pinus contorta) occurring where fire has recently affected the land. In higher elevations, White Spruce and Subalpine fir (Abies lasiocarpa) occur more frequently than Black Spruce, which occurs in more moist conditions (Meidinger et al., 1991). The main driver of change across all ecoregions in the study area is fire, often ignited by summer lightning storms (Yukon Ecoregions Working Group, 2004). 44   Figure 4.1: Map showing study area location and other characteristics near Watson Lake, YK. Elevation, ecoregions, provincial boundaries, LiDAR flight tracks, and the town of Watson Lake are displayed. 4.2.2 Datasets 4.2.2.1 LiDAR Data  Light Detection and Ranging (LiDAR) data are increasingly being used to predict and report a range of forest attributes (Wulder et al., 2012a; Næsset et al., 2009), however acquiring wall-to-wall LiDAR data remains costly, with logistical difficulties in remote environments (Wulder et al., 2012b). The Canadian Consortium for LiDAR Environmental Applications Research in coordination with the Canadian Forest Service acquired multiple flight paths of LiDAR data to 45  sample boreal forests on a continental scale transect acquisition plan in the summer of 2010. Thirty-four separate flight lines totaling over 25,000 km in length were flown spanning 13 UTM zones. All LiDAR data were collected using a discrete return LiDAR system (Optech ALTM 3100). LiDAR data were specified to be collected at 1,200m above ground level, with scan angles between +/- 15°, pulse repetition frequency of 70 kilohertz, and a pulse density of approximately 2.8 returns per a square meter with a swath width of 800 m. LiDAR metrics were calculated from the returns that fell within 25x25 meter (625m2) plots. Bater et al., (2011) calculated forest attributes for these LiDAR sampling flights by creating multiple linear regression models relating the height, volume, and density of LiDAR returns to coincident field and aerial photo based samples of forest inventory attributes estimated from published allometrics developed by Ung et al., (2008). These attributes include total aboveground biomass and its components (wood, bark, and branch biomass), tree height, dominant tree height, Lorey’s mean tree height, basal area, and gross stem volume. These attributes were predicted and stored in a database attached to plot centroids. Table 4.1 shows the R2 and RMSE values of each LiDAR predicted forest attribute. Generally, LiDAR derived models explained between 65% and 83% of the variance. Root mean square error values were variable, but all were under 25% of the mean. In total five LiDAR transects covered the study area representing 5% of the total distance of the entire continental transect, and approximately 3% of the land within the study area.   46   Table 4.1: Explained variance (R2) and Root Mean Square Error (RMSE) values of LiDAR transect predicted forest attributes from the 25x25 meter plots database. Metric Name Amount of Variance Explained R2  RMSE (Root Mean Square Error ) Unit of RMSE Total Aboveground Biomass 76% 1353.21 kg Stem Biomass 78% 1125.35 kg Stem Volume 80% 2.74 m3 Basal Area 65% 0.30 m2 Lorey`s Mean Height 83% 1.34 m   4.2.2.2 Ancillary Environmental Data Environmental data that describes elevation and land cover have proven useful in the classification and stratification of diverse landscapes (Rogan et al., 2003; Richetti 2000; Fazakas et al., 1999; Franklin et al., 1995). Two digital elevation models were combined from GeoBC (2012) and Geomatics Yukon (2012) to create a seamless 30 m DEM across the provincial boundaries for the study area. Land cover was derived from the Earth Observation for Sustainable Development of Forests (EOSD) dataset which was compiled from Landsat imagery representative circa year 2000 (Wulder et al., 2008a). The land cover classification of the study area was dominated by four broad categories: water, deciduous forest, coniferous forest, and mixed forest. The forest categories are further divided into three density categories (dense, sparse, and open), but were collapsed into single density categories for each forest class for this study. These land cover data represent the most detailed information available describing the unmanaged northern forests in our study area, and often throughout Canada, (Wulder et al., 2008a).  47  4.2.2.3 Landsat Data For the years 1984 to 2012, multiple Landsat images per year were selected and downloaded from the USGS based on three criteria. First, images closest to an anniversary date of July 19th (Julian day 200) were selected as they represented peak phenological conditions. Images were then selected based on estimated cloud cover amount (< 75 %). Finally cloud location within an image was considered as an additional criterion. Up to four images per year were selected, with the goal of creating annual cloud free composites as described below in the compositing section. Images used in the study are shown in Table 4.2.    48  Table 4.2: Landsat image dates and years used to construct spectral time series. Mean day of year shows composite’s average day of the year based on the images used to create that annual composite. Image dates that are italicized and bolded were dropped from further use before compositing due to poor radiometric calibration during image processing.YEAR 1st Image Day of Year 2cd Image Day of Year 3rd Image Day of Year 4th Image Day of Year Annual Mean Day 1984 161 177 193 257 193 1985 211    211 1986 166 198   182 1987 201 217 249  222 1988 188 220 236  215 1989 190 238 254  227 1990 161 193 209  188 1991 212 228   220 1992 167 183 231 247 215 1993 185 201 217 233 209 1994 188 204 220  196 1995 159 223 239  207 1996 194 210 223  209 1997 212 228 260  233 1998 183 199 215 263 215 1999 194 218   224 2000 261 277   235 2001 167 183 215 231 199 2002 218 250   250 2003 269 277   273 2004 176 200    2005 162 178 186  175 2006 165 173 197  178 2007 216 240 256  237 2008 179 219   199 2009 165 213   165 2010 168 200 208 264 210 2011 187 219 227 243 219 2012 198 230     214  Mean Annual Day of Every Year 211 49  4.2.3 Image Processing The Landsat imagery time series was subject to cloud masking, relative radiometric normalization, yearly compositing, and temporal segmentation following methods described in detail in Kennedy et al., (2010) and detailed below. Temporal trajectory metrics were then extracted and grouped by Tasseled Cap component and complexity, ranging from simple (such as individual spectral index values) to more complex (such as the amount of change in a spectral index relative to pre-change values). A cloud free image from 2010 matching the time of the LiDAR overpass was spatially segmented using homogeneity criteria to produce a set of forest units that resemble southern boreal forest inventory polygons. LiDAR derived forest attributes, ancillary environmental variables, and temporal trajectory metrics were then averaged within each of the forest units. Variable sets were constructed based on the Tasseled Cap components and complexity of the derived temporal trajectories and then used in random forests models to estimate aboveground biomass, stem biomass, stem volume, basal area, and Lorey’s mean height. 4.2.3.1 Temporal Segmentation 4.2.3.1.1 Cloud Masking For each Landsat image, the Function of mask (Fmask) cloud masking algorithm was used to create masks defining clear sky, water, snow, ice, cloud, and cloud shadow pixels (Zhu & Woodcock 2012b). The algorithm uses rules based on cloud physics and Landsat’s optical bands to first delineate clouds. Cloud shadows are located by comparing the cloud shapes to potential shadow locations on the ground and then the similarity between them to determine if that area is a shadow.  50  4.2.3.1.2 Radiometric Correction and Normalization Normalizing each image in the Landsat time series to a reference image was implemented to minimize variations in atmospheric conditions (Kennedy et al., 2010). A high quality cloud free image (3rd August 1998) close to midsummer and the center of the time series was used as the radiometric reference point. The radiometric reference image was calibrated using Landsat Ecosystem Disturbance Adaptive Process (LEDAPS). LEDAPS corrects the image to surface reflectance using known atmospheric conditions for that date and given image characteristics (Masek et al., 2006). All Landsat images in the time series were radiometrically normalized to the LEDAPS corrected image using the Multivariate Alteration Detection Calibration Algorithm (MADCAL) (Canty et al., 2004). MADCAL identifies invariant pixels between pairs of images and then applies orthogonal regression, thus performing a relative normalization. Of the available 84 images selected from the archive for this study, nine were excluded from further study due to poor normalization caused by unmasked haze (see images noted in bold and italicized text in Table 4.2).  4.2.3.1.3 Yearly Cloud Free Composites Image compositing followed methods described by Kennedy et al., (2010). From the four images within each year that were selected from the archive, image composites were created by selecting an observation which was closest to the anniversary date, and cloud free. Table 4.2 shows the mean Julian day for each year`s image using this compositing approach and relates the variation of the composites around the Julian day 200 anniversary date  51  4.2.3.1.4 Tasseled Cap Components The yearly composited and normalized image stack was then transformed to Brightness, Greenness, and Wetness using the Tasseled Cap transformation (Crist et al., 1984a). The Brightness component is mostly related to variability in soil reflectivity weighted particularly to the blue band. The Greenness component is defined as the perpendicular axis to Brightness, and contrasts the visible bands with the infrared bands. A third perpendicular axis is defined as Wetness and contrasts the visible and near-infrared bands with the short-wave infrared bands (Crist et al., 1984b). Of the three components, Tasseled Cap Wetness has most often been used to assess changes in vegetation and soil, as well as to forest structure (Skakun et al., 2003; Franklin et al., 2002; Franklin et al., 2000; Collins et al., 1996). Franklin et al., (2002) found that Wetness could detect large, stand-replacing disturbance and minor, non-stand replacing forest disturbance. Healey et al., (2006) also found that Wetness was one of the most useful spectral indices for forest structure characterization and structure change across a range of Douglas-fir (Pseudotsuga menziesii) and Ponderosa pine (Pinus ponderosa) forests in Washington. In the same study, Healey et al., (2006) was successful both measuring forest change and at classifying mortality and harvest related disturbances into intensity categories. Likewise, in lodgepole pine and white spruce forests located in British Columbia, Wulder et al., (2004b) found that Wetness showed a positive relationship to stand age, and was a strong indicator of forest growth. Additionally, Tasseled Cap Wetness in particular has been shown to be sensitive to structural changes in forest vegetation (Healey at el., 2006). For this study, we utilized Tasseled Cap Wetness as our principal vector of change due to its success in previous studies; however we also assess the 52  predictive power of the Brightness and Greenness Tasseled Cap components to provide a contrast and assess their potential use in a boreal forest environment.    53  4.2.3.1.5 LandTrendr This study made use of the LandTrendr algorithm (Landsat based detection of trends in disturbance and recovery) as described by Kennedy et al., (2010). LandTrendr examines the history of each pixel, referred to as a temporal trajectory (e.g. spectral time series) and consisting of yearly spectral index values. Up to six straight-line segments that statistically capture the overall shape of the trajectory are fit to spectral index values. The possible seven anchor index values are referred to as vertices in the fitted temporal trajectory and describe the start and end of each linear segment. By fitting these linear segments, changes or trends in a single spectral index can be better described by reducing noise that is inherent in time series datasets (Figure 4.2). Once vertices are identified and linear segments fit in the target spectral index (here Tasseled Cap Wetness), the fitted solution can then be applied to other spectral indices and the minimization of year-to-year noise across any other spectral index. In this study, we use the terms temporal segmentation and temporal segments to refer to the fitting process and the resultant segments, respectively.    54   Figure 4.2: Spectral time series, temporal segmentation, and derived metrics for a hypothetical forest unit. A number of metrics can be derived and calculated from the data describing a unit’s temporal trajectory, as illustrated. A number of metrics were calculated describing each pixel’s temporally segmented trajectory (Figure 4.2) including Tasseled Cap Brightness, Greenness and Wetness values at each vertex, as well as the change between those vertices. Vertex years were also identified, and the time between vertices defining the duration. Additionally, year to year change was calculated from the yearly fitted Brightness, Greenness, and Wetness values.    55  4.2.3.1.6 LandTrendr Spectral Time Series Metrics Individual temporal segments can be classified into broad categories of recovery, disturbance and no change. We defined recovery segments as those with an uninterrupted increase in a Tasseled Cap component, whether or not it was preceded by an initial decrease in the same index. A decrease in a Tasseled Cap component over any length of time was defined as disturbance. Thus variables describing the greatest numeric increase or decrease in Wetness were defined as the “Greatest Recovery” or “Greatest Disturbance”, respectively. Variables that describe the second-greatest numeric increase or decrease were defined as the “Second Greatest Recovery” or “Second Greatest Disturbance respectively.” The “Longest Recovery” or “Longest Disturbance” was defined by determining which segments of recovery or disturbance had the largest number of years the segment occupies in the temporal trajectory (Kennedy et al., 2010). Consecutive segments of like trends were combined to describe long term trends of different magnitudes. The “Last Monotonic Trend” was the union of the last temporal segment and any preceding segments that have the same directional trend, allowing grouping of temporally adjacent and consistently recovering or disturbed segments (Pflugmacher et al., 2012).   Once temporal segments were labeled as recovery, disturbance, or no change, the timing and change in those segments were used to calculate detailed statistics. The previously calculated segment duration and magnitude were used to calculate a simple rate of change, and from those rates year to year recovery and disturbance amounts were derived. To normalize this change in context to the initial values, “Relative Magnitude” was calculated as the amount of change divided by the initial vertex value of the temporal segment, and “Weighted Magnitude” was calculated as the magnitude multiplied by the duration of the temporal segment.  56  Many metrics describe individual segments within a temporal trajectory, but metrics that examine and quantify the entire trajectory can also be useful. Totals of recovery and disturbance can be quantified by summing the year to year disturbance or recovery change values from each temporal segment. Those totals can then be used to calculate the weighted magnitude of year to year disturbance or recovery and a simple ratio of recovery to disturbance. A more complex and previously unused metric named Normalized Disturbance and Recovery Ratio was calculated by subtracting the total recovery from the total disturbance and then dividing that by the sum of recovery and disturbance.  4.2.3.1.7 Metric Grouping We grouped temporal segmentation metrics based on Tasseled Cap component and complexity (Table 4.3). The first group represented Tasseled Cap component values for 2012. The next three groups (2 – 4) consisted of simple segment descriptor statistics for each Tasseled Cap component; the next three groups (5-7) utilized only complex temporal trajectory metrics organized around each Tasseled Cap component. Groups (8-10) utilized the combination of the simple segment descriptor statistics and the more complex temporal trajectory metrics, again separated by Tasseled Cap component and included land cover and elevation data.   57  Table 4.3: Variable groups used in Random Forest analysis to predict forest attributes, listed by scenario with temporal trajectory metrics and group name. Scenario Variable Group Name Temporal Trajectory Metrics Used 1 Tasscap Brightness, Greenness, and Wetness from a 2012 image 2 Simple Brightness LandTrendr Brightness Vertices, Temporal Brightness Segment Magnitudes & Duration, & Year Vertices 3 Simple Greenness LandTrendr Greenness Vertices, Temporal Greenness Segment Magnitudes & Duration, & Year Vertices 4 Simple Wetness LandTrendr Wetness Vertices, Temporal Wetness Segment Magnitudes & Duration, & Year Vertices 5 Complex Brightness Metrics Greatest Recovery, Second Greatest Recovery, Longest Recovery, Year to Year Recovery & Calculated Recovery Metrics, Greatest Disturbance, Second Greatest Disturbance, Year to Year Disturbance & Calculated Disturbance Metrics, Last Monotonic Trend & Calculated Last Monotonic Trend Metrics of Brightness 6 Complex Greenness Metrics Greatest Recovery, Second Greatest Recovery, Longest Recovery, Year to Year Recovery & Calculated Recovery Metrics, Greatest Disturbance, Second Greatest Disturbance, Year to Year Disturbance & Calculated Disturbance Metrics, Last Monotonic Trend & Calculated Last Monotonic Trend Metrics of Greenness 7 Complex Wetness Metrics Greatest Recovery, Second Greatest Recovery, Longest Recovery, Year to Year Recovery & Calculated Recovery Metrics, Greatest Disturbance, Second Greatest Disturbance, Year to Year Disturbance & Calculated Disturbance Metrics, Last Monotonic Trend & Calculated Last Monotonic Trend Metrics of Wetness 8 All Brightness Metrics Metrics from scenarios 2 and 5, Landcover and Elevation data 9 All Greenness Metrics Metrics from scenario 3 and 6, Landcover and Elevation data 10 All Wetness Metrics Metrics from Scenario 4 and 7, Landcover and Elevation data 11 All Important Brightness, Greenness, Wetness Metrics and Ancillary Variables Brightness, Greenness, and Wetness Vertex values 1 and 2, Brightness, Greenness and Wetness Last Monotonic Trend pre event value, Wetness Last monotonic trend post event value, Brightness Longest Recovery post event value, Wetness Longest Recovery pre event value, Brightness Greatest Recovery post event value, Greenness and Wetness Greatest Recovery pre event value, Landcover, and Elevation.   58  Once the initial models (scenarios 1-10) were developed, we then selected metrics with importance values (see Random Forests section) greater than 20 from the model results. These high importance metrics were placed into a new group (scenario 11) named Important All Landtrendr Brightness, Greenness, Wetness Metrics, and Derived Metrics, and Ancillary Variables. This group represents the important variables from each Tasseled Cap component’s model utilizing simple and complex temporal trajectory metrics and ancillary elevation and land cover data and was used to then estimate all forest attributes.  4.2.3.2 Spatial Segmentation Spatial segmentation was used to partition the imagery into discrete image objects (Powers et al., 2012). The resultant image objects allow for generalization of features and based upon the segmentation settings used are similar to manually delineated stands of a forest inventory (Hay et al., 2005). Segmentation was based on a cloud free image acquired on September 21st 2010 (Julian day 264) using the following parameters: scale 10, shape = 0.3, and compactness of 0.5 (settings after Wulder et al., (2003)). Comparable to forest inventory polygons, forest units were created from the segmentation results. Forest units were on average 6750 m2 and in size ranging between 900 m2 and 77,400m2 as described in Table 4.4. Although not a focus of the study, we acknowledge different spatial segmentation parameters will produce a range of forest units which is discussed in detail in Wulder et al., (2003). All LiDAR derived forest attributes, temporal trajectory metrics, and ancillary environmental variables were averaged into these forest units and used in the remaining analysis.   59  Table 4.4: General characteristics of forest units produced from multilevel object oriented segmentation of a Landsat image.       4.2.3.3 Random Forests The machine-learning algorithm “random forests” (Brieman 2001) was used to predict forest attributes as calculated from the LiDAR data using predictor variables extracted from the spectral time series metrics. Random forests have been used successfully in other forest attribute modeling studies (Falkowski et al., 2009; Hudak et al., 2008; Brieman 2001) aided by the inherent favorable characteristics of the algorithm including: 1) higher accuracies than single classification trees; 2) ability to utilize categorical as well as continuous data; 3) robust predictions are created by relying on multiple trees; 4) reliable internal error estimates; 5) the nonparametric approach that makes no assumption about the data distribution; 6) over fitting avoidance through the drawing of randomized  predictor variable subsets for each tree; and 7) the ease with which intercorrelated variables can be incorporated (Cutler et al., 2007; Prasad et al., 2006; Brieman 2001). Random forests has also been used to predict current and future tree distributions under climate change scenarios (Prasad et al., 2006) and forest attributes (Falkowski et al., 2009). For our implementation, we set the training/testing split at 70/30%, the maximum number of trees at 100, and the maximum levels per tree at 100. Importance, as calculated by Number of Forest Units: 56540     Characteristic Average Amount Unit Area 6750 m2 Total Aboveground Biomass 8870 Kg Stem Biomass 3805 Kg Stem Volume 4.65 m3 Basal Area 1.15 m2 Lorey’s Mean Height 13.17 m 60  summing the added improvement in accuracy gained at each node for all input variables was used to gauge the influence of variables in random forest models, and was used to establish which temporal trajectory metrics had the most predictive power.  We approached the analysis in two phases. We first focused on testing the predictive capacity and importance of the spectral time series metrics to estimate aboveground biomass, as compared to single date Tasseled Cap components. Using the results from the aboveground biomass estimates, high importance scoring metrics were selected from scenarios 8-10, creating a new group that was used to estimate all selected forest attributes. 4.3 Results Models of total aboveground biomass were developed using the listed sets of metrics (scenarios 1-10) in Table 4.3. A baseline model (scenario 1) using only Tasseled Cap component values from 2012 explained 13% of the variance with a high RMSE% value of 74% (Table 4.5, Figure 4.3). The next three groups (scenarios 2-4) utilized simple temporal trajectory metrics  generated with each Tasseled Cap component. Of those models, generally metrics based on Wetness led to higher R2 (38%) and lower RMSE% (62%) values. The use of more complex metrics (scenarios 5-7) provided modest additional explanatory power in estimations of total aboveground biomass (Table 4.5 & Figure 4.3), with Wetness based variables reporting the highest R2 (42%) and lowest RMSE% (61%) values. Combining simple and complex metrics by Tasseled Cap component with ancillary environmental variables (scenarios 8-10) again increased the models’ variance on average by an additional 12% of variance and reducing RMSE% by 6.3%. The models that used simple, complex, and environmental variables explained the most variance and 61  had the lowest RMSE% values for each Tasseled Cap component’s models. The aboveground biomass All Wetness Metrics model was the best of all models that utilized metrics based only on one Tasseled Cap component. That model (scenario 10) explained the most variance of scenarios 1-10 with an R2 of 50% and an RMSE% of 56% (Table 4.5, Figure 4.3).  Figure 4.3: Aboveground Biomass explained variance (R2) and Root Mean Square Error % (RMSE%) for random forests models using scenarios 1-10. Table 4.5: Explained variance (R2) and Root Mean Square Error % (RMSE%) values of random forests models used to predict aboveground biomass listed by scenarios 1-10 and group name. Forest Attribute Scenario Variable Group R2% RMSE% Aboveground Biomass 1 Tasscap 13% 74% 2 Simple Brightness Metrics 20% 71% 3 Simple Greenness Metrics 24% 69% 4 Simple Wetness Metrics 38% 62% 5 Complex Brightness Metrics 18% 72% 6 Complex Greenness Metrics 34% 65% 7 Complex Wetness Metrics 42% 61% 8 All Brightness Metrics & Ancillary Variables 32% 66% 9 All Greenness Metrics & Ancillary Variables 48% 57% 10 All Wetness Metrics & Ancillary Variables 50% 56%  62  Mean importance values for scenarios 8-10 are shown by variable type in Table 4.6. The values represent the relative importance of similar metrics, such as all Greatest Disturbance metrics, or all segment magnitudes. The mean importance values show the dominance of environmental variables in the random forests models, with a mean importance value (68) three times higher than the next ranked variable type (Table4.6). Simple vertex values for each Tasseled Cap component were ranked second most important variable type, though its value (17) was much lower than the first ranked variable type. The third highest ranked variable type was the Last Monotonic trend metrics, which reported a mean importance value (9) that was almost half of the second ranked variable type.  Table 4.6: Mean variable importance values for random forests models using scenarios 8-10, shown by variable type and metric. Scenarios 8-10 Model's Mean Importance Rank for Variable Type Mean Importance Value Ancillary Variables 68 Vertex Values 17 Last Monotonic Trend Metrics 9 Year to Year Recovery Metrics 7 Greatest Recovery Metrics 7 Year to Year Disturbance Metrics 6 Longest Recovery Metrics 6 Segment Durations 6 Disturbance ,Recovery, and No Change Years 5 Year to Year Recovery 4 Year to Year Disturbance 4 Greatest Disturbance Metrics 4 Disturbance and Recovery Ratios 4 Vertex Years 3 Segment Magnitudes 3 Second Greatest Recovery Metrics 1 Second Greatest Disturbance Metrics 1  63  Aboveground biomass, stem biomass, stem volume, basal area, and Lorey’s mean height were modeled using scenario 11’s variable set with R2 and RMSE% values reported in Table 4.7 and depicted in Figure 4.4. The amount of explained variance ranged from 62-65% for all forest attributes, however the RMSE% values differed by as much as 31%. The aboveground biomass model was the least accurate (RMSE% = 49%), with stem biomass and stem volume, reporting modestly lower RMSE% values of 44% and 42%. Basal Area had the second lowest RMSE% value of 34% and Lorey’s mean height reported the lowest value of 18%.  Table 4.7: Explained variance (R2) and Root Mean Square Error % (RMSE%) output values of random forests models. Models were used to predict aboveground biomass, stem biomass, stem volume, basal area, and Lorey’s mean height listed by scenario 11 and group name.    Forest Attribute Scenario Variable Group R2% RMSE% Aboveground Biomass 11 All Important  Brightness, Greenness, Wetness  Metrics and Ancillary Variables 62% 49% Stem Biomass 11 All Important  Brightness, Greenness, Wetness  Metrics and Ancillary Variables 63% 44% Stem Volume 11 All Important  Brightness, Greenness, Wetness  Metrics and Ancillary Variables 64% 42% Basal Area 11 All Important  Brightness, Greenness, Wetness  Metrics and Ancillary Variables 65% 34% Lorey`s Mean Height 11 All Important  Brightness, Greenness, Wetness  Metrics and Ancillary Variables 63% 18% 64   Figure 4.4: Explained variance (R2) and Root Mean Square Error % (RMSE%) for random forests models estimating each forest attribute using scenario 11. 4.4 Discussion 4.4.1 Understanding Variance and RMSE in Context with Previous Studies Reporting root mean square error as a percent of the mean (RMSE%) is the most often reported statistic for aboveground biomass estimation using optical remotely-sensed data . By normalizing with the mean, it is easier to compare studies in different areas, using different methods, with differing biomass amounts. Labrecque et al., (2006) reported RMSE values between 44% and 56.1%. Powell et al., (2010) reported RMSE%s of 68-87% for an Arizona study area and 61-69% for a Minnesota study area. Fazakas et al., (1999) reported RMSE% between 66.5% and 78.5%.  Our best aboveground biomass estimate reported a RMSE% of 49%, placing it within the bounds of previous work of aboveground biomass estimates using single date and time series of Landsat imagery.  65  The amount of variance explained in percentage (R2) is also often reported. The models of Hall et al., (2006) explained 65% of the variance in their biomass estimates using multivariate regression. Labrecque et al., (2006) explained 8.4% of the variance using direct radiometric relationships. Relying on land cover lookup values, Wulder et al., (2008b) reported a baseline value of 54% of the variance when compared to field measurements. Pflugmacher et al., (2010) found that some temporal trajectory variables explained as much as 69% of the aboveground biomass variance. Our best model explained 62% of the LiDAR estimated aboveground biomass which is in the range of reported values for aboveground biomass estimation using optical data.  Our model results using Tasseled Cap component values from 2012 (scenario 1) confirm the known difficulties in estimating aboveground biomass from single date imagery (Figure 4.3). However, the results also show that the inclusion of variables derived from temporal trajectories can greatly improve aboveground biomass estimates as indicated by higher R2 values and lower RMSE% values for scenarios 2-10. Unexpectedly, the use of only more complex variables in models (scenarios 5-7) did not substantially improve upon R2 values and RMSE% values when compared to models that used simpler metrics (scenarios 2-4). Models of stem biomass, stem volume, basal area, and Lorey’s mean height were created using the variable group (scenario 11) that utilized only high importance scoring metrics from scenarios 8-10.  Interestingly, best models for stem biomass, stem volume, basal area, and Lorey’s mean height reported similar R2 values (63%, 64%, 65%, and 63%). However, despite the R2 values staying relatively consistent between forest attribute models, the RMSE% values differed by as much as 31%. Lorey’s mean height reported the lowest RMSE% values (18%), considerably lower than next closest forest attribute’s RMSE%, basal area. Interestingly, basal 66  area is used in the calculation of Lorey’s mean height and so was used by Bater et al., (2011)in thier calculation of Lorey’s mean height from the field data, which was then used to estimate forest attributes from the LiDAR data in multiple linear regression models. It is worth noting that the Bater et al., (2011) regression’s R2 values were higher for Lorey’s mean height than basal area. Our estimates of Lorey’s mean height and basal area may reflect this difference in the LiDAR data, which may explain why our estimates of Loreys mean height are more precise than our estimates of basal area.  4.4.2 Metric Complexity  Spectral time series (e.g. temporal trajectory) metrics can vary in their complexity. Some metrics are intuitive and are the result of the temporal trajectory fitting process. Other metrics isolate key points in a spectral time series i.e. the Greatest Disturbance or Greatest Recovery, the Last Monotonic Trend. Still other metrics aim to sum the whole time series i.e. ratios of total disturbance to total recovery. As indicated by the R2 and RMSE% values in Table 4.5, the use of more complex temporal trajectory metrics provide modest, but not substantial increases in explained variance and error estimates when compared to simpler metrics. When simple and complex metrics from all Tasseled Cap components were combined, such as in scenario 11 (Table 4.7), the highest R2 and lowest RMSE% values were reported for aboveground biomass estimates, thus suggesting that the combination of simple and complex temporal trajectory metrics from each Tasseled Cap component are the best estimators of aboveground biomass. However, our study cannot confirm this, since those models using scenario 11’s variable set also made use of ancillary environmental variables.  67  4.4.3 Variable Importance Variable importance can be a useful metric for understanding the ability and power of individual or groups of metrics for forest attribute estimation. Firstly, environmental data describing the type of land cover and elevation are of paramount importance. This may have occurred because the mean importance value of many spectral time series metric groups in Table 4.6 simply does not capture the variability of importance values for individual metrics within that group. For example, the mean importance value of Wetness vertices 1-7 for scenario 10 was 27. However the previously unreported Wetness vertex importance values for each individual vertex 1-6 were 100, 64, 10, 7, 5, 3 and 1 in order. The first two vertices have much higher importance values than the later vertices, but when all the vertices were averaged together as a whole set in their own group, the mean importance value is much lower than the first two vertices. Thus the low mean importance value for the whole group obscures the high importance values of the first two vertices. Secondly, simpler metrics such as the vertex values may be better for estimating forest attributes than the more complex metrics. As a single set of variables the individual vertex values (Brightness, Greenness, and Wetness values at vertices 1-7) do not have temporal information. The vertex values describe the key points in time for each forest unit, highlighting periods of stability or change, regardless if that change is a disturbance or recovery. However those metrics describe the entire spectral time series of a forest unit, which is potentially useful since these metrics are much easier to calculate and account for disturbance and recovery. The relatively high rank of Year to Year Recovery Metrics and Year to Year Disturbance Metrics (Table 4.6) lends support to the idea that metrics describing the whole trajectory may be better for estimating forest attributes than more complex metrics that single out a period in time. 68  The Year to Year Disturbance and Recovery variables when combined essentially describe year to year change, whether that be positive or negative. When the two variable types are used together, they combination describes the same trends as the vertex values. However the key difference is that the vertex values describe key points in the history of forest unit, whereas the Year to Year Disturbance and Year to Year Recovery variables describe annual change.    69   Spectral Recovery Trajectories Illuminate Differences in Forest Chapter 5:Recovery  The previous chapter investigated how time series of remotely-sensed imagery can be used to estimate a range of boreal forest attributes. One aspect from the results of that work showed that metrics describing the recovery portion of a temporal trajectory were useful to estimate forest attributes. As a result, this prompted further interest and research into the post disturbance spectral recovery trajectory, which then lead to examining different methods to characterize spectral forest recovery trajectories. 5.1 Introduction The Canadian boreal can be stratified into ecozones which provide a spatial framework of ecologically homogenous units delineating distinct areas based on biophysical factors (Ecological Stratification Working Group 1996; Wiken 1986). The Boreal Shield ecozone is the largest unit stretching from Newfoundland into northern Saskatchewan with much of the area inaccessible and in a wilderness state. Despite the reported homogeneity of this ecozone, studies have suggested that the ecozone could be split along an east / west divide based on a number of distinct factors. Kurz et al., (1992) suggested the ecozone be split due to differing ecoclimatic conditions. Kull et al., (2006) split the ecozone citing colder and drier climes in the western section than the east. To characterize fire in the Canadian boreal forest, the ecozone was divided into eastern and western sections by Stocks et al., (2002). Likewise, citing differing climatic patterns and forest processes, Andrew et al., (2012) divided the Boreal Shield into two separate sections to aid their analysis of protected areas in Canada’s boreal zone. 70  As a result of differing levels of accessibility (Andrew et al., 2012), large areas of the Canadian Boreal Shield are not subject to direct anthropogenic influences with up to 35% and 40% of the forested East and West Boreal Shield ecozone sections not subject to any forest management practices (Wulder et al., 2007; Shvidenko et al., 2006). Thus the Boreal Shield can be further divided into two management zones: a generally southern managed portion that enjoys a relative wealth of descriptive data, and a generally northern unmanaged section that is less well characterized. The different management zones create a divide in forest information, which leaves a large part of the ecozone uncharacterized. The lack of spatially explicit forest information for more than a third of the ecozone is a problem when assessing both current and future climate change (Price et al., 2013). Increased productivity in temperature limited predominantly eastern boreal areas (Boisvenue et al., 2006), and increased quantity and severity of wildfire disturbances in drier western boreal areas (Flannigan et al., 2005) is expected, highlighting the need for disturbance mapping and then monitoring of the subsequent return of vegetation. Forest recovery has been monitored using image based time series using a variety of approaches. The most common, used with varying success, is to develop a time series of a given spectral index which is then related to biophysical parameters at a range of spatial scales (Chu et al., 2013). Griffiths et al., (2014) assessed variability in forest recovery both spatially and temporally across political jurisdictions in the Carpathian by defining spectral recovery as a combination of pre-disturbance spectral index values and the spectral magnitude of the disturbance. Kennedy et al., (2012) defined spectral recovery as a ratio of the magnitude of the disturbance event to spectral conditions five years post disturbance and illustrated recovery differences between 71  public and privately owned lands in the Pacific Northwest of the United States. In the same study Kennedy et al., (2012) showed that drier areas experienced slower recovery as opposed to more moist areas across management regimes, ownership, and political lines. Thus meaningful information about forest recovery across large areas can be derived from remotely sensed image based time series, providing useful data to managers and researchers alike and facilitating decisions about forest recovery.  Landsat time series analysis offers an effective approach to monitor large forested areas (Hermosilla et al., 2015b; Zhu et al., 2012a; Huang et al., 2010) like the Boreal Shield ecozone for disturbances and can enable characterization of the subsequent vegetative recovery (Kennedy et al., 2010). In this chapter we compare and contrast spectral forest recovery trajectories following forest disturbance across the east and west sections of the Boreal Shield ecozone by developing and then applying a methodology using Landsat time series imagery. To meet this objective we examine stand replacing disturbances and the subsequent spectral signals of recovery by constructing spectral trajectories of recovery from the post-disturbed spectral time series. These trajectories are then compared to nominally undisturbed forested signals, and the extent and recovery rate trends examined and compared across the ecozone.   5.2 Study Area The Boreal Shield ecozone (Figure 5.1) is generally characterized by rolling and hilly topography with many small lakes, streams and rocky outcrops. A precipitation gradient exists with higher amounts (1000mm) in the coastal east and less (400mm) in the more continental west (Urquizo et al., 2000; Ecological Stratification Working Group 1996). July average 72  maximum temperatures are 13°C for both sections, however the east typically has less severe winters with January average minimum temperatures of -1°C compared to the west with -20°C (Urquizo et al., 2000; Ecological Stratification Working Group 1996). The Boreal Shield ecozone is dominated by forests of black (Picea mariana) and white spruce (Picea glauca) with the more southerly portions having a wider mix of broadleaf treed vegetation such as white birch (Betula papyrifera), trembling aspen (Populus tremuloides) as well as an array of needle leaf species such as white (Pinus strobus), and red (Pinus resinosa) and jack pine (Pinus banksiana) (Ecological Stratification Working Group 1996).  As discussed, fire and harvest are the primary agents of stand replacing forest disturbances (Bergeron et al., 2004; Bergeron et al., 2000), with fires occurring more often over larger areas in the west than the east (Stocks et al., 2002). Less common are disturbances caused by insect infestations, wind and storm related damage, and disease (Urquizo et al., 2000).For the purpose of this study we follow the division of the Boreal Shield ecozone established by Kull et al., (2006), who cited the colder and drier climate in the west resulting in differing forest processes. As depicted in Figure 5.1 the ecozone can be divided by southern managed zone, where commercial forestry activities are present, and a northern largely unmanaged (non-commercial) zone.   73   Figure 5.1: Boreal Shield East and West ecozones and Canadian forest management zones within the eight study scene portions. 5.3 Methods 5.3.1 Data 5.3.1.1 Landsat The Landsat World Referencing System II (WRS-2) acquisition grid (185 x 185km) was used to create a stratified random non overlapping sample of Landsat scene footprints within the Boreal Shield ecozone (Masek et al., 2013; Kennedy et al., 2010; Healey et al., 2006; Wulder et al., 2001; Tucker et al., 2000). Images were required to be at least 50% within the Boreal Shield ecozone, with portions outside of the ecozone clipped from further analysis. The Earth 74  Observation for Sustainable Development of Forests (EOSD) land cover dataset was used to mask water bodies and to define areas of forest cover within the selected Landsat scenes (Wulder et al., 2008a). Deciduous, coniferous, and mixed forested type classes from the EOSD were collapsed into one category to create a binary forest layer, which was then used as a second filter to ensure that scenes were dominated by forest cover (> 80%). The scenes were further filtered by forest management zone, removing any scene that was not at least 50% within one management zone or another; if a scene covered both management zones, it was clipped to the portion within the largest management zone. Scenes were then randomly selected from the final candidate pool, while minimizing differences in total area where possible. Out of a possible 185 scenes covering the Boreal Shield ecozone, a total of 8 Landsat scenes were selected covering over 200,000 km2 or 25 % of the ecozone’s forested area (Table 5.1).  75  Table 5.1: Summary statistics about Landsat scenes used to create spectral time series. Management and inventory area, and boreal forested area are shown by individual path/row from which scene portions were used to construct the Landsat time series. The amount of area experiencing each type of management paradigm in each Boreal Shield section and as a percent of the total forested area in that Boreal Shield section and management zone is also displayed in Table 5.1. Table 5.1 additionally enumerates the total area under study per a Boreal Shield section and as a percent of the total forested area within that Boreal Shield section. The values reflected in Table 5.1 show that we have divided our study areas as evenly as possible between zones of management and Boreal Shield East and West sections.  5.3.2 Image Processing In summary, sets of Landsat imagery from each  scene was subject to radiometric correction, cloud masking, yearly compositing, and then temporal segmentation using the LandTrendr algorithm (Kennedy et al., 2012). Temporal trajectory metrics were then extracted from the temporal segmentation results and used to characterize recovery following stand replacing disturbances and compared between the east and west sections. Ecozone Name Boreal Forest Managed/ Unmanaged WRS2 Path and Row Portion Boreal Forest Area in Scene (km2) Study Zone Area per Management (km2)  Study Zone Area as a % of Ecozone per Manage-ment Total Study Zone Area in Ecozone (km2) Total Study Zone Area as a % of Ecozone Boreal Shield East Managed 17/26 31649 62592 7.60% 97865 11.89% 12/25 30943 Unmanaged 13/23 11571 35272 4.29% 16/24 23701 Boreal Shield West Managed 31/21 15742 46078 5.88% 105485 13.47% 27/24 30336 Unmanaged 26/23 29305 59406 7.58% 38/20 30101 76  5.3.2.1 Image Selection, Preprocessing, and Compositing Once selected, the imagery was extracted from the Landsat archive based on two criteria: 1) an acquisition date within the typical boreal zone growing season (which corresponds to Julian days 152-243 and the calendar months of June, July, and August); and, 2) cloud cover < 75%. A total of 1220 images were selected from the archive. The distribution of imagery over the years for each scene is displayed in Figure 5.2.  All imagery was ortho-rectified and converted to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Process (LEDAPS) (Masek et al., 2006). A cloud mask was produced for each image using the Function of mask (Fmask) algorithm (Zhu et al., 2012b). Pixels delineated as cloud, cloud shadow, ice and snow were masked and not used in the remainder of the analysis.    77   Figure 5.2: Distribution of dates and years of imagery for each Landsat footprint portion used for study in Chapter 5.Eastern Boreal Shield path row scene footprints are shown in blue while western Boreal Shield path row scene footprints are shown in orange. Diamond shaped symbols indicate images located in the management zone, while circle shaped symbols indicate images within the non-management zone. The straight black line indicates DOY 200, or mid-growing season. A total of 1220 images were used divided between the 8 scene footprints.   78  An annual time series for each of the eight Landsat scenes was constructed from the multiple images per year using rules developed by Kennedy et al., (2010): 1) for each year a base image was selected based on its proximity to July 19th (Julian day 200), coinciding with the mid-growing season; and, 2) all other images in order of their proximity to July 19th within the same year were then used to fill in any areas that had been masked in the base image. This compositing process was repeated for each of the 29 years in the analysis period.  5.3.2.2 Spectral Indices and Recovery Recent studies have highlighted differences in spectral recovery depending on the spectral index examined (Pickell et al., 2016; Banskota et al., 2014; Cuevas-Gonzalez et al., 2009; Epting et al., 2005). Frequently used in Landsat time series studies, the Tasseled Cap components are the result of a linear transformation of Landsat’s spectral bands into three orthogonal axes of Brightness, Greenness, and Wetness using published coefficients (Crist et al., 1984b). Previous studies have found that Tasseled Cap Brightness provides an indication of overall pixel albedo; Greenness provides an indication of vegetation photosynthetic condition; and Wetness is sensitive to changes in forest structure (Healey et al., 2006; Skakun et al., 2003; Franklin et al., 2002; Franklin et al., 2000; Collins et al., 1996; Cohen et al., 1995).  After a stand replacing disturbance, the Tasseled Cap components will have differing recovery trajectories responding to different properties of forest stands as the recovery process progresses: Brightness relates to the general albedo (Crist et al., 1984b) of the land surface and values are typically high after a stand replacing disturbance and remain high due to the lack of absorption from vegetation cover and 79  increased exposure of soil and non photsynthetic vegetation (Banskota et al., 2014). Throughout forest stand recovery Brightness values will reduce with time due to increasing vegetation cover and subsequent increased absorption of visible and NIR wavelengths. This reduction in Brightness is typically observed as the stand reinitiates and growing space becomes increasingly occupied. As a result Brightness is uniquely suited to detect if a disturbed forest stand is recovering and if vegetation is re-establishing post disturbance. Greenness contrasts the visible with the near and shortwave infrared wavelengths. The absorption of visible radiation by vegetation, and scattering driven by leaf by structure in the near infrared combine to produce high Greenness values as dense leafy forest canopies develop (Crist et al., 1984b). Low Greenness values are typically due to lack of photosynthetic vegetation, or vegetation that is under considerable stress. After a stand replacing disturbance, Greenness would be anticipated to increase quickly due to its sensitivity to green vegetation (Pickell et al., 2015).  Wetness was named for its response to moisture content in vegetation and soils, and relies heavily upon the shortwave infrared bands. High values of Wetness indicate an abundance of moist plant matter, while lower values indicate reduced vegetation (Crist et al., 1984b). Changes in forest structure have been related to changes in Wetness (Hansen et al., 2001; Cohen et al., 80  1995; Cohen et al., 1992), similar to the changes a forest stand undergoes during the recovery process. As a forest stand recovers, structure becomes more complex, creating a denser, moister, and multifaceted canopy with crown shadows and branching, all of which are reflected by increased Wetness values. Thus Wetness is well suited to observe forest recovery due to its ability to detect structural changes over time.  Throughout the recovery process forest stands undergo many physical changes, which can be observed using Landsat based Tasseled Cap time series. However, the capability to detect change is limited over time by the asymptotical nature of reflected radiation from a closed and dense forest canopy, as passive optical remote sensors can mostly detect changes on the topmost reflective surface, i.e. the canopy.  5.3.2.3 Temporal Segmentation The LandTrendr temporal segmentation algorithm was used to reduce year to year noise over the time series and segment the Wetness spectral values into distinct phases of disturbance, recovery, and no change (Frazier et al., 2014; Kennedy et al., 2010; Powell et al., 2010). LandTrendr divides each pixel’s time series using an iterative linear fitting model, establishing breakpoints that define periods of no change, losses, and gains within a temporal trajectory. Temporal trajectory metrics were then calculated from the resultant temporal segments (Frazier et al., 2014; Pflugmacher et al., 2012; Kennedy et al., 2010). To detect and isolate stand replacing disturbance events within the time series, we followed Kennedy et al., (2012) and applied thresholds to the largest (>700 change magnitude) and most rapid (<3 years) changes in Wetness to identify severe and short duration fire and harvest stand replacing disturbances. 81  Once disturbance events were identified using the Wetness trajectory, recovery trajectories that follow a disturbance event were isolated and only selected for further study if they exhibited a positive slope from the disturbance event to the end of the time period, thus ensuring a second disturbance did not occur. Using the year of disturbance, all recovery trajectories were then normalized by time since disturbance (see Schroeder et al., 2011). The mean and standard deviation of the three Tasseled Cap components was computed annually for all recovery trajectories in the east and west sections. For comparison, a stable, no change, forest mean threshold and standard deviation were also calculated from 18,652,249 no change pixels in the eastern section and 19,443,793 pixels in the western section. Finally, every recovering pixel in each year was classified using the annual rate of recovery into one of five Wetness recovery classes: far above average recovery, above average recovery, average recovery, below average recovery, and far below average recovery. The far above average recovery represented values 2+ standard deviations above the ecozone section mean while above average values were between 1 - 2 standard deviations above the average. The average recovery category was bounded by the ± 1 standard deviation from the recovery mean. Below average recovery was bounded by 1-2 standard deviations below the average. Pixels experiencing recovery that was more than 2 standard deviations below the mean were classified as far below average recovery. The amount of area for each disturbance and recovery category was also calculated as a percent of the total forest area under observation. 5.3.2.4 Cohen’s d statistic Traditional statistical significance tests using very large sample sizes can be misleading, as any small difference would be significant at some level and often unimportant in terms of causation 82  (Johnson 1999). We required a statistic that did not take into account our sample size and could determine meaningful levels of difference between our stable forest signals and recovery trajectory, and then between both East and West Boreal Shield recovery mean sets. The Cohen’s d statistic (Cohen 1988) can indicate similarity or difference in datasets regardless of sample sizes. Cohen’s d is a standardized statistic similar to z-scores, allowing the similarity or difference between two datasets to be compared. Although not often used in remote sensing applications, Cohen’s d has been applied in the biological and ecological fields to determine likeness or difference between samples. Studies have shown Cohen’s d ability to determine differences between datasets, when traditional significance testing may have shown less than optimal results (Cavin et al., 2013; Wesner et al., 2012; Nakagawa et al., 2007).  To calculate a d value, the means of two separate groups are subtracted and then divided by the pooled standard deviation as shown in equation 5.1:  Equation 5.1:  =̅ − ̅ where ̅ and ̅ are the means of two groups, and  is the pooled standard deviation and d is a Cohen’s d value. As described by Cohen (1988) and demonstrated in further studies (Cavin et al., 2013; Nakagawa et al., 2007), predefined difference thresholds set at d values of 0.2, 0.5, and 0.8  indicate small, moderate, or large differences between two datasets respectively. We utilized Cohen’s d to compare recovery trends across the Boreal Shield ecozone and to establish the degree of difference between stable and recovering vegetation, and difference between recovering vegetation in each Boreal Shield section. 83  5.4 Results 5.4.1 Wetness Recovery Trajectories and Cohen’s d Values After a stand replacing disturbance, Wetness has been shown to increase with time, and often corresponds to increases in forest structure (Schroeder et al., 2010). A positive Wetness recovery trajectory occurs in both the eastern and western Boreal Shield sections (Figure 5.3a). The eastern Boreal Shield initial Wetness value post disturbance is far below the stable range of no change vegetation, but after 15 years of recovery Wetness values have increased to be within the stable range.  Wetness values then continue to increase until the end of the Landsat record when they resemble the stable, no change vegetation mean. The western trajectory shows a similar trend, reaching the stable range by year 16, however not exceeding the stable mean during the analysis period. 84   Figure 5.3: Boreal Shield East and West Tasseled Cap Wetness recovery and Cohen’s d trajectories. a) Observed Wetness recovery trajectory of each Boreal Shield section. Shading for each section shows +/-1 standard deviation of the recovery means per a year. Straight horizontal lines indicate the stable means and +/-1 standard deviation range for each Boreal Shield section. b) Cohen’s d Wetness recovery trajectories for each Boreal Shield are shown with the predetermined levels of difference. Also shown is the Cohen’s d trajectory comparing the eastern Boreal Shield trajectory means to the western Boreal Shield trajectory means  The eastern and western Cohen’s d Wetness recovery trajectories (Figure 5.3b) indicate that both sections are initially very different, but end the analysis period with relatively little or no difference from their stable means. The d statistic for the eastern trajectory (2.09) indicates a large difference from its stable mean. After 25 years however the d value recovers to indicate no difference. The western trajectory’s initial d value of 1.3 also indicates a large difference from its 85  stable mean, but again after 22 years the trajectory crosses into the zone of no change forest. The trajectory ends at a d value of 0.37 indicating that the differences in Wetness of recovering vegetation and stable vegetation is small. When comparing the eastern and western trajectories to each other (Figure 5.3b) the Cohen’s d trajectory has an initial value of 0.01 indicating no initial difference between the two trajectories immediately following disturbance. After 20 years the recovery trajectories are markedly different from each other, and the comparison d trajectory ends at a value of 0.90 indicative of a large difference between the two Boreal Shield sections. A summary of the Wetness, Greenness, Brightness, and associated Cohen’s d trajectories are shown in Figures 5.3b-5.6b and Table 5.2. Each trajectory’s start, midpoint and end difference levels and d values are presented for the eastern and western recovery trajectories, as well as the comparison d trajectory. Table 5.2: Summarization of Cohen’s d recovery trajectories based on Brightnes, Greennes, Wetness and comparison trajectories. Values are shown at the start, midpoint, and end for each d trajectory comparing east and west Boreal Shield recovery means to their stable means, and comparing eastern Boreal Shield recovery means to western Boreal Shield recovery means.  East Recovery Means vs East No Change Stable Forest Mean West Recovery Means vs West No Change Stable Forest Mean East to West Comparison of Recovery Means Recovery Year 1 14 28 1 14 28 1 14 28 Brightness -0.84 -0.74 -0.49 -.036 -0.45 -0.48 0.51 0.62 0.21 Greenness 0.49 -0.60 -104 0.90 0.39 0.03 0.50 1.51 1.31 Wetness 2.09 0.98 0.02 1.31 0.97 0.37 0.01 0.61 0.90           Levels of Difference: Large Moderate Little None 86   5.4.2 Greenness Recovery Trajectories and Cohen’s d Values Like the wetness trajectories Greenness recovery is expected to increase as vegetation returns to the stand, and photosynthetic capability of the stand increases. Both East and West Boreal Shield Greenness trajectories conform to the expected positive trend following disturbance (Figure 5.4a).  The eastern trajectory begins below  the stable mean and within the stable range, but then quickly increases above the stable mean in the first 5 years. By year 19 recovery Greenness means are above the stable range, with the trajectory stable at that level until the end of the time series. The western Greenness trajectory has a very similar trend also starting below its stable range and increasing to higher than stable mean values after 20 years of recovery, but not exceeding the stable range by the end of the trajectory 87   Figure 5.4: Boreal Shield East and West Tasseled Cap Greenness recovery and Cohen’s d trajectories.a) Observed Greenness recovery trajectory of each Boreal Shield section. Shading for each section shows +/-1 standard deviation of the recovery means per a year. Straight horizontal lines indicate the stable means and +/-1 standard deviation range for each Boreal Shield section. b) Cohen’s d Greenness recovery trajectories for each Boreal Shield are shown with the predetermined levels of difference. Also shown is the Cohen’s d trajectory comparing the eastern Boreal Shield trajectory means to the western Boreal Shield trajectory means Cohen’s d Greenness recovery trajectories (Figure 5.4b) show similar shapes for the eastern and western Boreal Shield, but have different d ranges indicating similar recovery rates but different relationships to their no change stable means. The eastern d trajectory’s initial value of 0.49 indicates little difference from its stable mean, but after 5 years is not different and after 18 years is within the negative large difference range, where the trajectory ends with a d value of -1.04 88  indicating the recovery mean has exceeded the stable mean. The western d trajectory’s initial d value of 0.90 indicates a large difference from its stable mean, decreasing by year 17 to be within the no difference zone, remaining there until the end of the time series with a d value of 0.03. This suggests that the Greenness signal from recovering vegetation is indistinguishable from stable mature forest vegetation after 18 years in the western Boreal Shield, however the eastern recovery trajectory shows more complicated trends.  The eastern and western Greenness Cohen’s d comparison trajectory initial value of 0.50 indicates moderate difference (Figure 5.4b) between the two Greenness recovery trajectories. The d values quickly increase and by year 5 are within the large difference zone, remaining there until the end at a d value of 1.31. The high comparison d values (>0.8) from 4 years  post disturbance until the end of the comparison trajectory indicate that Greenness recovery is largely different between both Boreal Shield sections.  5.4.3 Brightness Recovery Trajectories and Cohen’s d Values As a forest recovers from a stand replacing disturbance Brightness will likely decrease representing an increase in canopy cover and a reduction in overall albedo. The eastern Brightness trajectory conforms to this expectation, however the western trajectory shows no trend (Figure 5.5a). The eastern trajectory starts above its stable range and ends above the stable mean. The western trajectory starts within the stable range and remains there throughout its entire trajectory finishing at similar initial Brightness values. Differences between the eastern and western trajectories are mirrored in their Cohen’s d trajectories (Figure5.5b). The eastern d trajectory’s initial value of -0.84 indicates a large difference from the stable mean. The trajectory then crosses into the moderate difference zone after 13 years and ends within the little difference 89  zone with a d value of -0.49. The western d trajectory’s initial value of -0.36 is within the little difference zone and stays within that zone for the remainder of the analysis period with a final d value of -0.48.   Figure 5.5: Boreal Shield East and West Tasseled Cap Brightness recovery and Cohen’s d trajectories.  a) Observed Brightness recovery trajectory of each Boreal Shield section. Shading for each section shows +/-1 standard deviation of the recovery means per a year. Straight horizontal lines indicate the stable means and +/-1 standard deviation range for each Boreal Shield section. b) Cohen’s d Brightness recovery trajectories for each Boreal Shield are shown with the predetermined levels of difference. Also shown is the Cohen’s d trajectory comparing the eastern Boreal Shield trajectory means to the western Boreal Shield trajectory means. 90  The Brightness comparison Cohen’s d trajectory decreases over time with an initial value of 0.51 within the moderate difference zone (Figure 5.5b). After 20 years the trajectory crosses into the little difference zone and ends with a d value of 0.21 near the lower bounds of the same zone signalling that recovering vegetation is always somewhat different between the two sections, however becoming more similar as  time since disturbance time increases.  5.4.4 Areal Disturbance and Recovery Trends  The total yearly area disturbed in the eastern Boreal Shield sample shows a constant pattern of disturbances between 1985-2000 (Figure 5.6). Between 2001-2012 there is a marked increase in disturbed forest area. The western Boreal Shield shows a variable pattern of disturbance from 1985-2000 when compared to the eastern section, and again shows an increase in annual disturbance from 2001-2012.   91   Figure 5.6: Annual amounts of disturbed, and continuously unabated recovering pixels for a) eastern study areas and b) western study areas.  The breakdown by trend rate for eastern and western Boreal Shield sections are shown in Figures 5.7 and 5.8 respectively as a percent of total forest area. A relatively stable disturbance trend in the eastern Boreal Shield from 1985 to 2000 is apparent from the yearly totals (Figure 5.7). The increasing annual percentage of disturbance is also evident from 2001-2012. The cumulative and classified recovery categories shown in Figure 5.7 demonstrates the variability in recovery rate after stand replacing disturbances. Most apparent is the marked relative increases in the areas 92  that are classified as far above average recovery rate, defined as being 2 or more standard deviations above mean Wetness recovery.   Figure 5.7: Boreal Shield East study areas spatial disturbance and recovery trends shown as a percent of the total forest area in the study areas. Disturbance is calculated annually. Recovery is cumulative and the categories of recovery are classified into five categories based on the recovery trajectory means and standard deviations from the Wetness trajectories in Figure 5.3a.  Variable disturbance trends characterize the western Boreal Shield between 1985-2000 as is evident from the yearly totals (Figure 5.8). An increasing amount of disturbance from 2001-2012 is evident as well as an increase in the area undergoing far above average recovery, similar to the eastern section. Throughout the same period, the amount of area classified as having an average 93  rate of recovery remains stable from 2001-2012,  at a level of about 0.5% of total forested area of the Boreal Shield.   Figure 5.8: Boreal Shield West study areas spatial disturbance and recovery trends shown as a percent of the total forest area in the study areas. Disturbance is calculated annually.Recovery is cumulative and the categories of recovery are classified into five categories based on the recovery trajectory means and standard deviations from the Wetness trajectories in Figure 5.3a. 5.5 Discussion In this research we examined the utility of Landsat time series analysis to characterize disturbance and recovery trends across the boreal forests of the Canadian Boreal Shield ecozone. We discuss the results under two broad headings: temporal trajectory considerations and 94  recovery characterization, and then our case study partitioning of the East and West Boreal Shield. 5.5.1 Temporal Trajectory Considerations & Recovery Characterization Forest recovery refers to the reestablishment of key forest biophysical variables following a disturbance event, and is a process not a state (Frolking et al., 2009). Forest recovery rates can differ due to severity, frequency and type of disturbance, site characteristics, and climate (Bolton et al., 2015; Frolking et al., 2009, Oliver et al., 1996), Post disturbance forest recovery can take multiple pathways including: 1) no recovery (e.g. land use conversion), 2) recovery towards a predefined stocking or structure (e.g. plantation forestry), or 3) natural recovery. Critical to any study of forest disturbance and recovery, regardless of location, is the rate that forest vegetation is returning from the disturbance event. However, the multiple ways in which forest recovery can be defined can lead to confusion when linking recovery to spectral time series. Chief among the concerns is that ecological and silvicultural definitions of forest recovery are not dependent on the actual presence of a mature forest, but more so predicated on provision of ecosystem services or emergent characteristics regarding treed vegetation height and canopy cover requirements.  Landsat time series can be used to characterize forest disturbance and forest recovery as shown in this work. The information generated from the Landsat time series aids in better understanding forest processes. Recovery extent, timing, and magnitude can now be reported and compared across large areas to inform on regional or temporally varying trends.  The terminology used by Oliver et al., (1996) when describing forest stand dynamics is particularly useful when describing spectral forest recovery trajectories, as it grounds spectral 95  forest recovery in easier to understand terms that are more generalizable and easier understood by a wider audience than spectral indices.  Oliver et al., (1996) detail a generally applicable model of forest stand dynamics briefly summarized here with additional insight on how each Tasseled Cap component is affected. A stand replacing disturbance opens growing space for new individuals (Oliver et al., 1996) causing an increase in albedo that in turn is reflected by increased in Brightness values. Individuals compete for resources until the canopy is closed, which causes large increases in Greenness and moderate increases in Wetness values, and a sharp decrease in Brightness. The stand reduces in density over time and the canopy increases in height (Oliver et al., 1996) causing further increases in Wetness over time. An understory redevelops as growing space and resources are released by the canopy (Oliver and Larson 1996), which is mostly undetectable to each Tasseled Cap Index since it is occurring below the canopy. Shade tolerant individuals grow from the understory into the canopy supplanting dominant canopy species (Oliver and Larson 1996) potentially causing further change all three Tasseled Cap indices. The stand is then governed by gap dynamics caused by small disturbances (Oliver and Larson 1996 until a stand replacing disturbance occurs, which will restart the recovery process (Oliver and Larson 1996).  While the spectral recovery trajectories described above are simplified archetypes, the reality will be different as the combined elements that drive recovery (including initial disturbance severity, proximity of regenerative source seed bank or undisturbed stand, and competition) can vary through a stands recovery, and thus effect the signals of recovery observed using Landsat time series. For example, the apparent linearity of our recovery trajectories is due to the selection for recovering stands that experience no further disturbances within our time series. In reality the 96  recovery process will vary based upon the nature and amount of vegetation present over the course of forest recovery over time and could result in more complex and subtlety varied recovery trajectories. The within pixel averaging of spectral vegetation conditions will also serve to mute the year-on-year variability in spectral vegetation recovery. Likewise, forest recovery monitoring using passive optical remote sensing are, by definition, limited to observing changes in reflected light from forest canopies. Previous efforts have demonstrated the difficulty inherent in assessing successional state from remotely-sensed datasets (Song et al., 2007; Schroeder et al., 2006; Song et al., 2003). As a result, the spectral forest recovery trajectories do not necessarily directly correspond to changes in forest successional stages or may only be able to report on a small amount of variation in biophysical variables. Rather spectral trajectories essentially track the reestablishment of vegetation towards canopy closure or when the spectral signals of a recovering forest stand resemble that of nearby undisturbed vegetation or its previous state before the disturbance occurred (Kennedy et al., 2012; Buma 2012; Schroeder et al., 2007;).  In this research we found that less than 1% of disturbed forests in both Boreal Shield sections showed no signs of recovery, with nearly all disturbed forests showing signs of recovery within 5 years. This implies that Boreal Shield forests overwhelmingly initiate recovery processes after stand replacing disturbances. This is in agreement with previous descriptions of disturbance and subsequent recovery as key processes that shape the patchwork of forest structures, types, and ages critical to maintaining the diversity of the Canadian boreal zone (Bonan et al., 1989, Heinselman 1981).  The spectral recovery trajectories constructed in this chapter depend upon changes in Wetness through time. The shortwave infrared bands are weighted heavily in the calculation of Wetness when compared to the visible and near infrared bands (Crist et al., 97  1984b). As a result it is well suited to detecting changes from disturbance events such as fire. If a disturbance or recovery event produces an observable change in other portions of the spectrum other than the shortwave infrared, it is possible that this change may not be detected by the temporal segmentation algorithm. Thus care should be taken when selecting a spectral index and knowledge of the disturbance types of interest and their characteristics is useful.  It is important to recall this link between the spectral channels used in a given index (or weight given in a combined index) with what is portrayed as recovery. An index will likely cross a non-disturbed forest threshold sooner with a visible and near-infrared driven index than an index driven by longer wavelengths. In this research, we have been conservative by focusing on the SWIR-driven Wetness index which has been shown to reflect changes in structural complexity and shadowing and is therefore more sensitive to the presence of trees, than early successional herbs and shrubs (Franklin et al., 2000; Collins et al., 1996; Cohen et al., 1995). Alternatively, initial increases in vegetated cover are well characterized with shorter visible and near-infrared wavelengths (as used in Greenness). This more rapid saturation of the visible and near-infrared driven indices (i.e. Greenness) offers initial insight that vegetation is present but has less capacity for providing information on increasing structural complexity (Healey et al., 2006; Skakun et al., 2003; Franklin et al., 2002; Franklin et al., 2000).   5.5.2 Case Study: Partitioning of the East and West Boreal Shield  We were able to characterize forest recovery using Landsat time series, focusing our work on the known east/west division of the Boreal Shield ecozone. One of the rationales for splitting the Boreal Shield ecozone into two distinct features was that the regions experienced different ecoclimatic conditions, disturbance regimes, and forest processes. Using our temporally 98  segmented datasets describing the disturbance and recovery trends in the both Boreal Shield sections, we suggest that characteristics of disturbance and recovery for the east and west Boreal Shield study areas are also different, further reinforcing the division of this ecozone.  The annual area disturbed shows different patterns between the east and the west sections. From 1985-2000 there is a relatively constant area disturbed in the east. In contrast the west has a much more variable pattern of disturbance over the same time period (Figure 5.6). Interestingly, both the east and west show the same increasing disturbance trends in the 2001-2012 periods, but at differing amounts with more disturbances in the east than west.  The recovery trends of most eastern and western spectral indices follow the expected trend of an initial low value rising over time. However, Brightness trajectories are anticipated to begin with an initial high value and decrease over time. When the east and west Brightness section’s trajectories are compared, large differences become apparent; the eastern trajectory follows a decreasing trend, while the western trajectory show no clear trend. We interpret this difference as related to the type and characteristics of disturbance that is most likely occurring. Bergeron et al., (2004) demonstrated that fires are less frequent in the eastern than the western Boreal Shield and fire size is on average much larger in the west (Stocks et al., 2002), which can lead to more patchiness and mixed severity effects. The clear decreasing pattern of Brightness in the East implies disturbances that are homogeneous entirely clearing a forest stand and leaving few if any remnant patches resulting in a uniform recovery signal. The lack of clear Brightness recovery trajectory in the West implies a more heterogeneous situation with mixed fire severities leaving more remnant patches that diminish the Brightness signal through the presence of more vegetation and decreased reflectivity in the visible wavelengths.  99  Spectral recovery trajectory patterns relative to their stable mean also highlight differences in the East and West ecozone sections. For example most eastern trajectories increase to their stable mean values (e.g. Wetness) and in some cases exceed their stable mean (e.g. Greenness). In contrast the western trajectories rarely exceed their stable mean, with most trajectories remaining below the mean throughout the analysis period (e.g. Wetness). We suggest these spectral recovery differences also relate known patterns of forest stand recovery between the two sections. After a stand replacing disturbance, eastern Boreal Shield forest stands have been shown to initiate with regenerating broad leaf shade intolerant treed vegetation, which have higher Wetness, and Greenness levels than the shade tolerant conifers that often replace as the stand matures (Bergeron et al., 2007; Chen et al., 2002; Bergeron et al., 2001; Bergeron et al., 2000). After a stand replacing disturbance, western Boreal Shield stands typically initiate from the surviving seedbank and mature towards even aged shade tolerant conifers or more complex and layered conifer stands if another stand replacing disturbance does not occur (Brassard et al., 2006; Bergeron et al., 2000). This pattern is shown in the western Boreal Shield in the spectral recovery trajectory with lower Greenness, and Wetness values corresponding to the coniferous dominated re-initiating stand.  The classifications of recovery rates show differences between the east and west (Figures 5.7 and 5.8). While both areas show an increase in the above average recovery rate towards the end of their time series, this is more evident in the east than the west. Also the amount of average recovery in the western study areas remains relatively stable from year 14, the midpoint of our analysis period, until the end of the time series, while the eastern Boreal Shield study areas average recovery increases over the same period of time.  100  The Cohen’s d statistic enabled the direct comparison of recovery patterns between the eastern and western Boreal Shield (Figures 5.3b-5.6b, Table 5.2). Three of four Cohen’s d east to west comparison trajectories commence within zones of small or moderate difference and end with large difference between the two sections, indicating that their recovery is observably different. The difference between eastern and western Wetness recovery is especially important because it is sensitive to changes in forest structure. Noticeably, the Wetness comparison d trajectory shows the most change from its initial point of no difference (d value 0.01) ending with large differences (d value 0.90). This essentially means that although stand clearing disturbances may reduce the signal in disturbed areas in the east and west to similar Wetness levels, they recover from that level in different manners, as indicated by their Wetness recovery trajectories and their Cohen’s d trajectories. Finally the comparison d trajectory shows that by the end of the time series, that east and west Boreal Shield recovery is largely (d>0.8) different.  It is important to recognize the factors that drive the variability in recovery between these two sections of the Boreal Shield ecozone. A known precipitation gradient exists, along with differences in winter average minimum temperatures, contributing to differing climates (Urquizo et al., 2000; Ecological Stratification Working Group 1996) and as a result observed differences in vegetation recovery rates. Differences in disturbance frequency and size are already well known (Stocks et al., 2002) for the two Boreal Shield sections which will also impact the observed spectral forest recovery. Additionally, soil conditions could also affect forest recovery rates (Lecomte et al., 2005). The status and variability of all these factors combine to drive the observed differences in the east and west spectral forest recovery trajectories.  101   Spectral Recovery Rates are Not Fixed over Time and Space Chapter 6:The previous chapter investigated how to characterize post disturbance recovery trajectories, which then provided insights into how forest recovery processes differ between and East/West division of the Boreal Shield ecozone. Based on the resultant differing spectral forest recovery trends, this prompted further research into the stability or variable nature of spectral forest recovery rates.  6.1 Introduction Climate and disturbance are the two of the most important factors that combine to shape the Canadian boreal landscape over time (Brandt et al., 2013). First, climate primarily controls where and which tree species grow and adapt in forested areas, with annual tree growth limited by short growing seasons and severe winters (Gauthier et al., 2015). Second, disturbances drive change in the boreal, particularly fires, which occur frequently, over large areas, and are critical for ecosystem functions (Stocks et al., 2002). As a result of the interplay between climate and disturbance, Canadian boreal forest ecosystems are a mosaic of patches with varying age, structure, biodiversity, productivity and species composition (Weber et al., 1997). Boreal tree species are wide ranging generalists and have adapted to frequent disturbance. These same species utilize multiple post-disturbance recovery methods that occur over long periods of time. This post-disturbance recovery period is a critical establishment period to monitor the health of boreal forests (Gauthier et al., 2015) and is highly susceptible to changes in climate, which in turn can impact the essential ecosystem goods and services provided by boreal forests (Gauthier et al., 2015; Gauthier et al., 2014; Brandt et al., 2013).  102  Boreal forests will be altered extensively by a changing climate (Gauthier et al., 2014; Brandt et al., 2013; Price et al., 2013). Currently, boreal ecosystems are undergoing alterations in phenology (Lemprière et al., 2008; Colombo 1998), productivity (Nemani et al., 2003; Jarvis et al., 2000), and disturbance (Flannigan et al., 2005; Flannigan et al., 2000; Stocks et al., 1998), all attributed to changes in climate. With respect to disturbance, a changing climate has had mixed effects on forest recovery. For example, three key tree species within the Canadian boreal forest have shown divergent responses to a changing climate. Black spruce (Picea mariana) dominated systems showed an increase in growth under cooler and wetter conditions, while jack pine (Pinus banksiana) ecosystems displayed increased growth with warmer temperatures and increased spring precipitation (Brooks et al., 1998). Additionally, white spruce (Picea glauca) trees have reacted positively to warmer spring temperatures with increased annual growth, but also negatively to warmer summers with decreased annual growth (Wilmking et al., 2004) likely caused by a lack of available moisture (Barber et al., 2000).  Spatially extensive temporal trends describing increasing or conversely decreasing vegetation quantities and vigor have been found occurring across the Canadian boreal using coarse and medium (100-30 ) grained remotely-sensed data. Importantly, those remote sensing-based trends derived from spectral indices in those studies have previously been tied to physical properties of forests (Fraser et al., 2014; McManus et al., 2012; Pouliot et al., 2009) with those same indices able to represent changes to forest canopy cover or structure (Schroeder et al., 2011). While some have found a widespread trend of increasing vegetative over time that is associated with a thickening canopy (Olthof et al., 2007; Slayback et al., 2003; Myneni et al., 1997), Goetz et al., (2005) found that Canadian boreal forests experienced a negative trend between 1981-2003 103  related to reductions in vegetation abundance (Bi et al., 2013; Beck et al., 2011). For example changes over time in Normalized Difference Vegetation Index (NDVI), the spectral index often used in coarse grained studies, have been linked to active photosynthetic vegetation quantity (Turubanova et al., 2015; McManus et al., 2012; Cuevas-Gonzalez et al., 2009). Used more often in fire severity studies, the Normalized Burn Ratio (NBR) spectral index and its change over time has been linked to forest structure properties (Epting et al., 2005; Wulder 1998), and may be more suitable than NDVI for forest recovery tracking (Pickell et al., 2016). Spectral trends found at coarse scale have been shown to be influenced by the inter-instrument calibration, sensor drift, and preprocessing steps (Alcaraz-Segura et al., 2010). Moreover, lower spatial resolution studies can fail to capture the fine scale mosaiced nature of the boreal forest landscape and are often in disagreement with research done at finer spatial scales (Fraser et al., 2011; Alcaraz-Segura et al., 2010; Wulder et al., 2010). In contrast to known issues with inter-instrument calibration of broad-scale sensors (Alcaraz-Segura et al., 2010), finer spatial resolution Landsat data (30m) has fewer inter-sensor calibration issues (Vogelmann et al., 2016; Markham et al., 2012), although issues regarding inter-sensor calibration of certain spectrum wavelenthgs are noted to persist (Sulla-Menashe et al., 2016; Ju and Masek 2016). Detailed and annual forest recovery information is now extractable from fine spatial scale wall-to-wall remotely-sensed datasets, principally due to free and open data access and data storage and processing capacity (Kennedy et al., 2014). The United States Geological Survey Landsat data archive provides an ideal data set for boreal-wide research that offers a well calibrated data record initiated in 1972 with complete spatial coverage in fine detail (Wulder et al., 2012c). Larger areas and longer time periods of spectral data can now be used to inform on the finely 104  detailed process of forest recovery, especially after wildfire-caused stand replacing disturbances that are typical of the Canadian boreal. For example, Pickell et al., (2016) used Landsat time series analysis across a range of forested boreal bioclimatic zones and observed differing spectral forest recovery rates for each zone. Likewise, Frazier et al., (2015) examined spectral forest recovery using Landsat time series data and found spectral recovery differences between two ecozones were related to distinct forest recovery processes. Research has indicated that boreal forest recovery rates vary spatially (Bartels et al., 2015) and have been suggested to vary temporally (Chu et al., 2016), but post-disturbance forest recovery rates and how they may vary over time have been little explored (Ju and Masek 2016, Goetz et al., 2006). In this chapter we examine boreal forest recovery following fire under a period of changing climate to determine how recovery rates have changed over the past 26 years, both temporally and spatially. We focus on the Boreal and Taiga Shield ecozones as representing nearly half of Canadian boreal zone, ensuring a wide array of boreal forest conditions and fire severities are considered. The 2,976,480 km2 study area and wall-to-wall mapping affords the ability to compare all observed forest fire disturbances and their subsequent recovery. To do so we: 1) use multiple spectral forest recovery metrics derived from Landsat time series data that inform on different aspects of post-fire spectral recovery; 2) detect trends in post-fire spectral forest recovery rates first at a 100 km spatial scale to visualize detailed rate trend patterns and then second at the ecozone analysis unit scale to understand broad environmental change; and, (3) determine if statistically significant spectral recovery temporal trends exist within the time series. Significant results are then discussed in relation to the disturbance and recovery regimes 105  of Canadian boreal forests, as well as the benefits and limitations of spectral forest recovery metric approaches to assess forest recovery used in this study.  6.2 Study Area Our study area covers the forested boreal areas with the Canadian Boreal Shield and Taiga Shield ecozones (Figure 6.1; Ecological Stratification Working Group 1996; Brandt 2009). For scientific analyses, the Boreal Shield and Taiga Shield ecozones are often divided into eastern and western sections as a result of differences in climate and disturbance regimes (Bolton et al., 2015; Andrew et al., 2012). Eastern sections generally experience a less harsh winter and more annual precipitation than their western counterparts (Frazier et al., 2015; Kull et al., 2006; Kurz et al., 1992). Both western ecozone sections experience more forest fires (Stocks et al., 2002); forest harvesting occurs in all ecozones, and insect infestations play a larger role in eastern than western disturbance regimes (Brandt et al., 2013). When combined together, these two ecozones account for 49% of the Canadian Boreal zone (Brandt 2009) (Table 6.1).    106  Table 6.1: Boreal and Taiga Shield East and West study area in respect to the Canadian boreal zone as defined by Brandt et al., (2009),  and ecozones as defined by Ecological Stratification Working Group (1996). Boreal Area in Ecozone Total Area of Ecozone % Boreal of Total Ecozone % Boreal Area of CDN Boreal Area Boreal Shield East 952673 km² 1236647 km² 77.04% 15.51% Boreal Shield West 818726 km² 818980 km² 99.97% 13.33% Taiga Shield East 692545 km² 789149 km² 87.76% 11.28% Taiga Shield West 550051 km² 598086 km² 91.97% 8.96% Total 3013995 km² 3442862 km² 87.54% 49.07%   Figure 6.1: The Boreal and Taiga Shield East and West study area is displayed with a false color Landsat composite (R; Band 5, G: Band 4, B: Band 3). The four ecozones under study are displayed: the East and West sections of the Boreal and Taiga Shield. The featured ecozone areas are also restricted to the Boreal ecosystem as defined by Brandt (2009).   The West and East Boreal Shield ecozones are dominated by forests and characterized by many small lakes and streams interspersed between rocky outcrops with rolling and hilly topography. 70°W80°W90°W100°W60°N50°NTaiga   Shield       WestBoreal Shield WestTaiga ShieldEastBoreal Shield EastUnitedCanadaStates107  A precipitation gradient exists varying from higher (1000mm per year) amounts in the coastal east, and lesser (400mm per year) in the continental west. January mean temperatures are -20°C and -1°C in the west and east showing a pronounced colder west to warmer east temperature gradient (Urquizo 2000), and the mean annual temperature is approximately 3°C. Boreal Shield forests are dominated by black  and white spruce  stands, with southerly portions having a wider mix of broadleaf and coniferous species, i.e. white birch (Betula papyrifera), trembling aspen (Populus tremuloides), white (Pinus strobus), red (Pinus resinosa) and jack pine (Ecological Stratification Working Group 1996).  Located north of the Boreal Shield Ecozones, the Taiga Shield Ecozones are physically divided by Hudson Bay. The topography is marked by rolling uplands punctuated by rocky outcrops, and glacial moraines and eskers. Differences in temperature and precipitation reflect the coastal eastern and continental western sections; the eastern portion receives 500-800mm per year of precipitation and experiences a milder winter with mean January temperatures near -11°C. Annual precipitation in the Taiga Shield West ranges 200mm lower than that of the eastern section, and average January temperatures are -25°C. The northern edge of the Taiga Shield delineates the edge of where tree growth is possible; however, common throughout the ecozone are forest stands dominated by short black spruce and jack pine, with white spruce dominated mixed wood stands including balsam fir (Abies balsamea), trembling aspen, balsam poplar (Populus balsamifera), and white birch (Ecological Stratification Working Group 1996).  108  6.3 Data & Methods 6.3.1 Spectral Data Annual, cloud-free, and seamless surface reflectance image composites and change products derived from Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) data were utilized for this analysis (Hermosilla et al., 2015a, Hermosilla et al., 2015b). These data were generated using the Composites-to-Change (C2C) algorithm, which is briefly described here. First, annual best-available-pixel (BAP) spectral composites were produced using all peak growing season imagery (White et al., 2014) and then the Normalized Burn Ratio spectral index (NBR) was calculated (Key et al., 1999). A temporal segmentation algorithm was then applied to the NBR time series to determine and quantify change over time and to infill data gaps with proxy spectral values (Hermosilla et al., 2015a). Importantly, the NBR spectral index ranges between 1 and -1, with low values indicating low or no vegetation present, while high values are known to represent dense healthy vegetation. NBR has been successfully used to: 1) assess burn severity of fires (French et al., 2008); 2) detect and classify forest disturbances (Wulder et al., 2009), and; 3) inform on forest attributes (Frazier et al., 2014; Pflugmacher et al., 2012). Lastly, change objects were geometrically, temporally and spectrally characterized, and then classified by change type. Hermosilla et al., (2015b) found these C2C methods to be accurate over the forested area of Saskatchewan and reported user’s and producer’s accuracy for wildfire of 97.8% and 96.7%, respectively.  109  6.3.2 Methods 6.3.2.1 Data Stratification and Filtering  Analysis of trends in spectral forest recovery rates was first undertaken at a 100 km grid cell scale to ascertain if there were any spatial patterns in forest recovery trends within the ecozones themselves. Then each ecozone analysis unit was analyzed as a whole to determine if spectral forest recovery rates showed a consistent temporal trend. In order to better capture stand initiating forest recovery characteristics, the Landsat time series data were filtered such that only spectral recovery trajectories from moderate and high severity fires (typically stand replacing fires) were considered (Figure 6.2). Fires were selected using a threshold value of NBR change (dNBR) that represented a level of severity high enough to initiate the growth of a new stand. Hall et al., (2008) examined dNBR values for fire disturbances in the Boreal Shield West determining that moderate and severe fires corresponded to dNBR values greater than 0.284.  110   Figure 6.2: Example of spectral data showing pre-disturbance NBR values, Fire Year, and dNBR values for three sites in the Taiga Shield West. The first column is associated with site 1, the second column site 2 and the third column site 3; this pattern is repeated in Figures 6.3 & 6.7, and each site represents an area affected by a fire that occurred in the early, middle or later periods of the time series respectively. The top row shows the pre-disturbance NBR values, while second row shows which year fire occurred. The bottom row shows the dNBR values for each fire year.  6.3.2.2 Spectral Forest Recovery Metrics  Three spectral forest recovery rate metrics were calculated on a yearly per pixel basis using a five year recovery window with the knowledge that the initial first few years of tree species reestablishment after disturbance is critical for forest recovery monitoring (Figure 6.3) (Kennedy 100°W110°W120°W55°N90°WNBR-1-0.250.010.20.350.451.0Fire YeardNBR0.00.160.320.50.660.841.0SITE 1 SITE 2 SITE 3200619861996SITE 1SITE 2SITE 3111  et al., 2012; Goetz et al., 2006). Three metrics were used due to their ability to gauge different aspects of spectral recovery as described below. 112   Figure 6.3: NBR time series from the year of disturbance through five years post disturbance. Time since disturbance increases from the top row and downward, with each site corresponding to the sites shown in Figure 6.2. Generally, NBR increases with time since disturbance.  SITE 1 SITE 2 SITE 3NBR Y0-1-0.250.010.20.350.451.0NBR Y+1-1-0.250.010.20.350.451.0NBR Y+2-1-0.250.010.20.350.451.0NBR Y+3-1-0.250.010.20.350.451.0NBR Y+4-1-0.250.010.20.350.451.0NBR Y+5-1-0.250.010.20.350.451.0113  The first metric calculated was the Recovery Indicator (RI) as proposed by Kennedy et al., (2012), which compares the spectral disturbance magnitude to the recovery magnitude five years post- disturbance (Figure 6.4)  shown in Equation 6.1:.  Equation 6.1:   =∆	∆ where ∆NBRregrowth is five year total change in NBR after disturbance and ∆NBRdisturbance is the change in NBR due to disturbance. We modified the RI to work with our unfitted post-disturbance spectral data by allowing ∆NBRregrowth to use the maximum NBR value from four or five year post-disturbance shown in Equation 6.2:  Equation 6.2: ∆	 = 		(, 	) − 	 where Max(NBRy+5, NBRy+4) is the maximum of either NBR value four or five years after disturbance, and NBRy0 is the NBR value the year disturbance occurred. We modified ∆NBRdisturbance to use the average of NBR values from two years before disturbance to allow for our unfitted spectral data shown in Equation 6.3:  Equation 6.3:  ∆ = 		 − 	 where NBRpre is the average NBR value from two year pre-disturbance and shown in Equation 6.4: 114  Equation 6.4:    	2 where NBRy-2 and NBRy-1 are the NBR values from two years each before disturbance. Unique to the RI metric is the scaling that occurs as spectral recovery is based on the amount of disturbance. This scaling can be helpful in areas such as the boreal that have diverse and heterogeneous signals before and after disturbance. This metric has been used to detect differences in spectral forest recovery across forest ownerships, climate, and substrates (Kennedy et al., 2012). A value of zero indicates that no spectral forest recovery has taken place, while one indicates that an equal amount of spectral forest recovery and disturbance has occurred. Values greater than one indicate that more spectral forest recovery has occurred than spectral forest disturbance, and RI values practically range between zero and two.   Figure 6.4: The Recovery Indicator (RI), adapted from Kennedy et al., (2012). We adapted this metric to, utilize the average of two years pre-disturbance NBR values to determine a pre-disturbance condition, which was then used to calculate ∆NBRdisturbance. The maximum NBR value from four or five years post-disturbance and the year of disturbance were used to calculate ∆NBRregrowth.  Ratio of Eighty Percent (R80P), adapted from Pickell et al., (2016) was the second metric calculated, and represented the number of years required for a disturbed pixel to return to 80% of 115  the pre-disturbance NBR value (Figure 6.5). For this study, the metric was modified to be a ratio of the spectral value five years post-disturbance compared to 80% of the pre-disturbance average of the NBR values from the two preceding years, and represents the amount of spectral recovery relative to the pre-disturbance condition shown in Equation 6.5: Equation 6.5:  80 = 		(, 	) ∗ 	 .8 where the maximum NBR value from four or five year post-disturbance is used. This recovery rate metric is different from the RI in that it gauges the amount of spectral recovery by the spectral signal that existed before disturbance. This scaling can contextualize the spectral recovery by relating the post-disturbance conditions to pre- disturbance conditions. The R80P metric as used by Pickell et al., (2016) was able to distinguish spectral recovery differences between biogeoclimatic zones, suggesting either differing pre-disturbance conditions or forest recovery processes occurring during their time series. R80P values are similar to the RI, with a value of zero indicating that spectral forest recovery did not occur, while a value of one indicates that NBR values have recovered to be equivalent to 80% of the pre- disturbance values. Values greater than one indicate that spectral forest recovery has exceeded that of the 80% of pre- disturbance threshold.  116   Figure 6.5: The Ratio 80% of Pre-disturbance (R80P) metric, adapted from Pickell et al., (2016). Additionally, the NBRpre denominator was calculated as the NBR average of two year prior to disturbance and the maximum NBR value from four or five years after disturbance is used in the numerator. The final metric was Year on Year Average (YrYr) which is the average of spectral values from each year in the reference period to the prior year for five post-disturbance years (Figure 6.6) shown in Eqution 6.6:  Equation 6.6: Year on Year Average spectral recovery rate metric.          5 and yearly differences are calculated for a,b,c,d, and e shown in Equations 6.7-6.11: Equation 6.7:   	 			 Equation 6.8:  	  	 	 Equation 6.9:    	 	 Equation 6.10:    	 	 Equation 6.11:    	 	 117   where NBR0, NBR+1, NBR+2, NBR+3, NBR+4,and NBR+5 are NBR values from the year disturbance occurred through five years post-disturbance as denoted by their subscript. The YrYr metric only considers spectral recovery and is neither referenced to the disturbance magnitude, nor the pre-disturbance values. This metric is important because it offers a view into annual spectral recovery, instead of a using a single point to mark recovery, unlike the RI and R80P metrics. This metric provides an insight on the average annual change after disturbance and can help capture annual post-disturbance growth than a single measurement. A YrYr value of zero indicates that for the five year average, zero spectral forest recovery had occurred, while positive values indicate the average gain in NBR in five years that occurred.   Figure 6.6: The Year on Year Average (YrYr) metric. This metric was created for this study, and is the average of yearly change in a spectral index for five years after disturbance. Year to year differences a, b, c, d, and e are calculated by subtracting spectral index values of one year from the previous year. Areas that experienced a disturbance before 1986 or after 2006 were removed; either an insufficient amount of spectral data existed before disturbance to establish a spectral trend or not enough spectral data existed after disturbance to assess spectral recovery in our five year window (Figure 6.7). The annual pre-disturbance values, dNBR values, and the spectral forest recovery 118  rates were then combined to create five time series, which were then each subjected to statistical testing for trend detection.  Figure 6.7: Sample of spectral forest recovery rate metrics calculated for the same areas shown in Figure 6.2 and 6.3. Note that the YrYr values are lower than the other metrics, and class divisions are adjusted accordingly. 6.3.2.3 Trend Detection A Thiel-Sen slope estimator was used to determine the average yearly change of the pre-disturbance values, dNBR fire severity, and spectral forest recovery rates over time. A Mann-Kendall test was then used to determine whether a time series exhibited a significant monotonic increasing or decreasing trend (p < 0.05). Non-significant trends were not considered important and were discarded from further analysis, although their Theil-Sen slopes were still displayed in the results for completeness. Only Mann-Kendall significant Thiel-Sen slope results are reported YrYr00.050.10.150.20.25-0.05RI00.10.20.30.40.5-0.05R80P00.10.20.30.40.5-0.05SITE 1 SITE 2 SITE 3119  on in the results. Additionally, time series with less than six fires were determined to have an insufficient quantity of data to test for trends. Non-significant results for any of the spectral forest recovery metric time series only indicate that the rate at which forests are spectrally recovering in that cell has not changed consistently over time. A non-significant result cannot indicate that forests are not recovering after disturbance nor does it imply that the spectral forest recovery rate is static. However, a non-significant result does indicate that either the rates of spectral forest recovery do not change greatly over time or that those rates have changed erratically and do not show a consistent trend.  Annual average pre-disturbance average NBR value and NBR disturbance magnitude time series were tested to determine if spectral recovery rates might be increasing as a result of an overall NBR increase over time or increased disturbance magnitudes, which could affect the spectral recovery rates. Spectral forest recovery rate time series were then tested for significant trends. A significant trend in spectral forest recovery rates indicates that the rate at which vegetation is recovering within the cell is consistently changing over time, with direct impacts for forest productivity and sustained health. Results from both statistical tests were first assessed at the 100 km grid scale and then by each ecozone analysis unit. 6.4 Results 6.4.1 Non-Recovery Trends, Fire Frequency and Extent  Average pre-disturbance NBR values and NBR disturbance severity time series trend detection results are shown in Figure 6.8. Approximately 85% of cells had no significant trend in NBR pre-disturbance values or disturbance severity over time. The remaining cells either exhibited a 120  trend in pre-disturbance values (8.1%) or a trend in disturbance severity (5.7%) or very rarely both (.1%). Interestingly, the cells that reported a disturbance severity trend average a very small Thiel-Sen slope of -0.000399 per a year, indicating that NBR assessed severity has declined over time. The cells that produced a trend in pre-disturbance average values conversely contain a very small positive Thiel-Sen slope of 0.000191, suggesting that a very small increase in NBR values over time has occurred in those cells. These results should not alter the recovery rate metrics time series as they both represent very small change over time. When compared to the possible range of NBR, the dNBR fire severity and pre-disturbance value Theil-Sen slopes represent 0.05% and less than 0.01% average yearly change, respectively.   Figure 6.8: Presence of significant trends over time in either pre-disturbance values or disturbance severity in the 100 km grid cells. A trend in pre-disturbance values indicates that NBR values before disturbance either increase or decrease consistently, while a trend in disturbance severity indicates that dNBR values have either increased or decreased over time. 60°W70°W80°W90°W100°W110°W60°N50°NNone Pre-DisturbanceValuesBothEcozone BoundsInsufficient DataDisturbanceSeverity121  The number of fire-affected years and percent area disturbed by fire throughout the time series for each 100 km grid cell is presented in Figure 6.9 and Figure 6.10. As expected there are more moderate and severe fires and more area burnt in western than eastern ecozone analysis units. Eighty percent of grid cells exhibited 10 or more years of fire, and overall 35% of all cells experienced between 10 and 15 years of fire in their time series. Although 97% of all grid cells experienced fire, 78% of the grid cells had less than 10% of their area burnt, and only 1% of all cells recorded more than 30% burnt area.   Figure 6.9: The number of years a cell experienced a moderate or high severity fire disturbance within the time series. 60°W70°W80°W90°W100°W110°W60°N50°NFire Disturbance Years01515 1021Ecozone Bounds122   Figure 6.10: The percentage of each grid cell disturbed by a moderate or high severity fire within the time series. 6.4.2 100 km Grid Spectral Forest Recovery Rate Trend  The cell-based Thiel-Sen regression slopes and significant Mann-Kendall test results for the RI metric time series are shown in Figure 6.11. Overall 13% of cells with sufficient data (more than five years of fire disturbances occurring between 1986 and 2006) report a significant trend, and 69% of all significant cells contain a positive Thiel-Sen slope. The largest proportion of significant trending cells (60%) contains a slope between 0.0-0.025, with an average slope among them of 0.0178. This indicates that RI values, which measure the rate of spectral recovery, were larger at the end of the time series than the beginning. Since RI can only practically range between 0 and 2, the average RI rate is increasing by 0.89% yearly for those positive and significant trending cells. Interestingly, 20% of all the cells in the Taiga Shield East are significant and positively trending, which is the highest proportion of any ecozone.  60°W70°W80°W90°W100°W110°W60°N50°N% Cell Distrubed by Fire20 0515 103040Ecozone Bounds123   Figure 6.11: Significantly trending grid cells and Theil-Sen slopes for the Recovery Indicator (RI) metric time series. Cells with significant Mann-Kendall results are outlined in thick gold borders, indicating a trend over time is present for that metric. Grey areas indicate that an insufficient amount of data was available to calculate a trend. A histogram of significant and non-significant slopes is also displayed. 124  For the R80P time series (Figure 6.12), only nine percent of cells with sufficient data reported a significant recovery rate trend and 70% of those significant cells contain positive Theil-Sen slopes. The majority (55%) of significant slopes range between 0.0 and 0.025, with an average yearly increase of 0.01565. This shows that R80P values among those cells increase by 0.01565 annually throughout the time series, indicating that the rate of spectral forest recovery is increasing with time for those positive and significant cells. Similar to the spatial pattern shown for the RI results, the Taiga Shield East contains the highest proportion of significant trending cells.  125   Figure 6.12: Significantly trending grid cells and Theil-Sen slopes for the Ratio 80% of Pre-disturbance (R80P) metric time series. Cells with significant Mann-Kendall results are outlined in thick gold border, indicating a trend over time is present for that metric. Grey areas indicate that an insufficient amount of data was available to calculate a trend. A histogram of significant and non-significant slopes is also displayed. 126  YrYr grid cell time series results show significant trends (Figure 6.13) in 13% of all cells, with 64% of significant grid cells having a positive slope. All positive significant slopes average 0.0021 NBR units annually, indicating the difference between NBR recovery rates are on average 0.0021 higher from year to year for those cells. Again, a contiguous pattern of significant grid cells is located in the Taiga Shield East, with that ecozone containing the majority (32%) significant cells, and significant cells occupying 21% of all cells in the ecozone.  127   Figure 6.13: Significantly trending grid cells and Theil-Sen slopes for the Year on Year Average (YrYr) metric time series. Mann-Kendall significant trending grid cells are outlined in a thick gold border, indicating a trend over time is present for that metric. Grey areas indicate that an insufficient amount of data was available to calculate a trend. A histogram of significant and non-significant slopes is also displayed. 128  The combination of all significant trend results across the three metrics shows two general areas where significant increasing trends tend to spatially cluster (Figure 6.14). First, as previously shown in Figures 6.11, 6.12 and 6.13, a contiguous block of significant cells is located in the Taiga Shield East near 70°W and 55°N. Grid cells reporting at least one significant trend occupy 26% of all cells in that ecozone. Interestingly, 25% of grid cells in the Taiga Shield West ecozone report at least one significant trend; however, no cells report significant trends for all three time series of spectral forest recovery rates. Moran’s I spatial autocorrelation results show that the significant metrics are clustered over the study area, but when each ecozone is solely considered a clear pattern may not be evident (Table 6.2). Both Eastern ecozone analysis units had a significant clustered pattern, but the Western ecozones had non-significant and random patterns of significant recovery rate trends.  Figure 6.14: Number of significant trends in a grid cell, across all spectral forest recovery rate metric time series Mann-Kendall results. 60°W70°W80°W90°W100°W110°W60°N50°NEcozone BoundsInsufficient DataTotal SignifcantMetric Trends13 2 0129  Table 6.2: Global Moran’s I spatial autocorrelation results indicating that dispersed, random, or clustered patterns exist in the total number of significantly trending metrics displayed in Figure 6.14. Global Moran's I Z Score P Value Result Boreal Shield East 0.16 2.85 0.00 Clustered Boreal Shield West -0.05 -0.64 0.52 Random Taiga Shield East 0.29 4.00 0.00 Clustered Taiga Shield West 0.09 1.23 0.22 Random All Cells 0.24 6.74 0.00 Clustered  6.4.3 Ecozone-Level Analysis  The annual NBR pre-disturbance average and disturbance magnitude time series for each the ecozone analysis units showed no significant trend over time, except for the Boreal Shield West. The Boreal Shield West pre-disturbance value time series had a significant trend with a very small negative Thiel-Sen slope (-0.00063 NBR units). That small negative slope value is 0.03% of the whole range of NBR and could affect R80P or RI time series, but not the YrYr recovery metric time series. Neither the Boreal Shield East nor West ecozone analysis units had a significant trend for any of the three spectral recovery rate metrics (Figure 6.15). In contrast, the Taiga Shield East had a significant and positive trend for the R80P rates time series (Figure 6.15). Theil-Sen analysis confirmed that the R80P recovery rate values increased by a total of 18% over the 26 year period. Likewise, the Taiga Shield West ecozone showed significant positive trends for the RI and YrYr recovery rate time series. Theil-Sen analysis of the RI time series showed a 9% increase in spectral forest recovery rates when start and end values were compared. The YrYr Theil-Sen slope value indicated that the spectral recovery rate was 12% greater at the end of the time series than the beginning. Most importantly, these significant results for the Taiga Shield 130  West and East indicated that the rate of spectral forest recovery consistently increased over time for the ecozone.  Figure 6.15: Ecozone analysis level spectral forest recovery rate trend detection results. Recovery Indicator (Panel A), Ratio 80 Percent (Panel B), and Year on Year Average (Panel C) is displayed by ecozone, as well the Thiel-Sen regression displaying a significant Mann-Kendall trend denoted by a thicker line. Panel D shows the number of pixels per year for each ecozone that were used in the subsequent five spectral recovery period 0.40.50.60.70.80.040.060.080.100.12Recovery Indicator  (RI)Year on Year Average (YrYr)Boreal Shield East Recovery Rate Metric RatesBoreal Shield West Recovery Rate Metric RatesTaiga Shield East Recovery Rate Metric RatesTaiga Shield West Recovery Rate Metric RatesNon-Significant-Mann-KendallTheil-Sen TrendSignificant-Mann-KendallTheil-Sen TrendNumber of ObservationsBoreal Shield EastBoreal Shield WestTaiga Shield WestTaiga Shield East0.40.50.60.70.80.91.0Ratio of Eighty Percent  (R80P)1985 1990 1995 2000 2005Years1985 1990 1995 2000 2005Years0200000400000600000800000100000012000001400000Number of Observations16000001985 1990 1995 2000 2005Years1985 1990 1995 2000 2005Yearsa)b)c)d)131  6.5 Discussion 6.5.1 Spectral Forest Recovery Rate Trends We examined three different time series of spectral forest recovery rates to determine the existence of trends, that is, if spectral recovery rates showed a consistent year to year change. The results showed that annual spectral forest recovery rates consistently increased for a portion of the 100 km cell locations and also for two of four ecozone analysis units over time. A change in spectral forest recovery rates can be crucial to future sustained boreal forest health because changes in forest recovery rates can have important impacts on annual growth, productivity, essential ecosystem services provided by Canadian boreal forests (Gauthier et al., 2015).  The 100 km grid cell results display an array of significant and non-significant spectral forest recovery rate trends, confirming that spectral forest recovery rates are not fixed over time or across space. Those significant trends suggest an acceleration of spectral recovery over time for some areas of the boreal forest. In contrast to 100 km cell results, the Boreal Shield ecozones did not report any trend in spectral forest recovery rates. However, the Taiga Shield ecozone analysis units collectively had significant and positive trends for all three spectral forest recovery rate time series. The positive trend detected in the Taiga Shield East R80P time series and the positive trends detected in the RI and YrYr time series indicate that spectral forest recovery rates have increased over time.  6.5.2 Spectral Forest Recovery Metrics Spectral forest recovery metrics aim to quantify the spectral change after disturbance by using a spectral vegetation index that indirectly assesses the quality and quantity of vegetation (Chu et 132  al., 2015). The spectral index used in this study has been used successfully to detect and classify disturbances (Hermosilla et al., 2015b; Schroeder et al., 2011), as well as inform on the age and structure of recovering forests (Pflugmacher et al., 2012). However caution should be used when interpreting spectral forest recovery metrics, as they may not be directly related to some physical forest recovery indicators (Bartels et al., 2015; Kennedy et al., 2012; Frolking et al., 2009).  The three spectral forest recovery metrics used in this chapter address different possible methods to spectrally measure forest recovery. For instance, the RI compares the spectral recovery magnitude to the spectral disturbance magnitude, e.g. the metric is anchored by disturbance amounts. That disturbance magnitude is quantified in terms of dNBR, which can be related to fire severity (Hall et al., 2006). Likewise, recovery magnitude has been shown to be related to forest structure throughout the recovery period (Pflugmacher et al., 2012). In contrast to the RI, recovery using the R80P metric is related to pre-disturbance values. Because NBR values have been used to predict canopy cover (Kennedy et al., 2010) and forest structure (Frazier et al., 2014), those pre-disturbance NBR values are then a critical measure of forest recovery.  As opposed to evaluating recovery based on pre-disturbance conditions, the YrYr metric is calculated as the average year-to-year NBR change for five years after disturbance and importantly gave each post-disturbance year separate and equal consideration. Critical to this metric, greater NBR values are generally related to increasing forest structure and canopy cover (Wulder et al., 2009). These three metrics contextualize the measured spectral forest recovery differently, providing separate insights into post-disturbance boreal forest recovery.  133  A five-year post-disturbance period was used to examine spectral forest recovery in this study. Although the temporal window used to observe recovery in our research was relatively short, it allowed more data to be examined and a longer time span to be analyzed; a ten-year post-disturbance window would limit the last recovery rate to be collected in 2002.  Kennedy et al., (2012) proposed a five-year window stating that it was critical to monitor the initiation of recovery, while Pickell et al., (2016) noted that the majority of boreal forests recover spectral values quickly. Moreover, a meta-analysis by Bartels et al., (2106) noted that forest height and canopy cover meet forest recovery definitional criteria between 5 and 10 years after fire or harvest disturbances in most areas of the Canadian boreal forest. Those findings indicate that observational windows longer than about five years are not necessarily required for observing boreal forest recovery. Furthermore, Gauthier et al., (2015) identified the initial phase of forest recovery as an important indicator of future vitality and this would be captured better with a relatively short observational window. Our spectral recovery results build upon previous research efforts that have aimed to assess the recovery of boreal forests using Landsat time series. Fraser et al., (2014) examined recovery resulting from fire disturbances in Canadian boreal forests and found that multiple spectral vegetation indices increase with time since disturbance. Similarly, Ju et al., (2016) noted in their North American study of Landsat-based NDVI that areas recovering from fires primarily showed an increase in NDVI related to accumulating recovering vegetation. Furthermore, Pouliot et al., (2009) found that spectral recovery rates varied over time in Canadian boreal forests, although their study only applied to a limited area. These results compare favorably to our results that indicated spectral forest recovery rates are not fixed over time or across space.  134  6.5.3 The Need for Continued Boreal Monitoring  The boreal forests of the Taiga Shield and Boreal Shield ecozones are both dominated by stands of black and white spruce (Ecological Stratification Working Group 1996). Fire is an integral part of these forests, and is often severe enough to be considered stand replacing (Weber et al., 1997), which is followed by regeneration in most cases (Pickell et al., 2016; Frazier et al., 2015). This primarily fire-based disturbance cycle has created the resultant patchwork mosaic of ages, compositions, structures, biodiversity, and productivity that comprise the current boreal conditions (Brandt 2009). Adding to this complexity is that post-disturbance recovery can follow multiple paths that are affected by many localized factors, such as fire severity, pre-disturbance species composition and proximity to undisturbed patches (Chen et al., 2002). Thus, the disturbance and recovery cycles of the boreal forest ensure that any landscape is not static. Continued monitoring is required to enable a comprehensive understanding of ecosystem change over time.  The initial forest reestablishment and regrowth period is critical to monitor to assess the continued health of the boreal forest in the face of a changing climate (Gauthier et al., 2015). Our research examined that initial five year period of recovery over a 26 year time series and found that spectral forest recovery rates increased over time for the Taiga Shield East and West ecozones. Recovery was assessed by using the NBR spectral index, which primarily responds to chlorophyll content and moisture in leaves (Key 2006). Thus, higher NBR levels indicate greater quantities of vegetation (Wulder et al., 2009). The increasing NBR-based recovery rates found in this study could be a result of the already noted productivity changes occurring throughout the boreal, which would correspond to increased amounts of photosynthetic vegetation (Nemani et 135  al., 2003; Jarvis et al., 2000). In addition, growing season length alterations directly affect seedling reestablishment and annual sapling growth (Lemprière et al., 2006; Colombo 1998). The capability of species to re-establish in a disturbed area is directly impacted by any change in disturbance characteristics (Johnstone et al., 2010; Flannigan et al., 2005; Flannigan et al., 2000; Stocks et al., 1998). Most importantly, tree species capable of exploiting the new environmental conditions will thrive, while others will decline (Gauthier et al., 2015); however, gains in annual growth due to the surmised favorable conditions can be offset by declines caused by another factor (Brandt et al., 2013; Brooks et al., 1998).  136   Conclusions Chapter 7:7.1 Research Innovations The research undertaken in this dissertation provides advances in the spectral characterization of forest recovery in the Canadian boreal. New methods were developed and then applied to landscapes and regions that result in novel insights about boreal forests. These advances and insights include:  • New methods for characterizing and comparing spectral recovery and forest recovery; • A time series approach to estimate forest attributes for a non-inventoried boreal forest; • Tools to assess forest spectral Forest recovery over large areas and through time; • Examination of forest recovery patterns across the division of the Boreal Shield ecozone between East and West sections; and  • A better understanding of broad environmental change occurring in the Boreal and Taiga Shield ecozones 7.2 Methods for Charactering Boreal Forest Recovery with Spectral Time Series The application of Landsat spectral time series to monitor forest recovery over large areas across the entire archive is still a relatively new endeavor. To date most coarse scale and AVHRR-based efforts did not emphasize forest recovery, instead focusing on undisturbed vegetation trends or mapping the spatial extent of forest disturbance (fires). The detection and characterization of forest recovery after a disturbance has more recently been a focus of study using Landsat data, which importantly recognizes that forest recovery and spectral recovery are not always synonymous (Kennedy et al., 2012). Thus the link between spectral recovery and forest recovery can be dependent upon the user, with the definition of forest recovery variable among forest 137  information consumers (Bartels et al., 2016). The research presented in chapters 4, 5, and 6 acknowledged this difficulty and encouraged grounding spectral change in terms defined by what a sensor could reliably detect, i.e. spectral signals and their change over time. By focusing on these signals, I have added to a better understanding of the spectral and forest recovery.  Spectral recovery and forest recovery are not mutually exclusive, and care must be taken to define what spectral recovery may mean in terms of forest recovery. The generalized model of stand dynamics outlined by Oliver and Larson (1996) becomes especially useful when discussing spectral forest dynamics. For example, forest disturbance as described by Oliver and Larson (1996) is an event that causes a disruption of forest function or structure in time and space. They also go on to describe recovery using structure and canopy characteristics as important indicators of forest recovery. This disturbance and recovery language has been widely adopted by the remote sensing community to describe spectral time series that show change. Rooting those spectral changes in a physical manner provides a reliable method for interpreting spectral changes in forests as shown in Chapters 4, 5, and 6. The research in Chapter 5 specifically details how the model presented in Oliver and Larson (1996) can be used to more comprehensively interpret spectral change by describing physical and spectral forest recovery at the same time aiding our understanding of the dynamic nature of boreal forests. This adds to the current understanding of spectral recovery by adding a lens to view what spectral recovery may mean for forest recovery.  Many previous efforts to study forest recovery using remotely sensed image time series confined their research to a single spectral vegetation index. For example, many of the coarse resolution studies highlighted in Chapter 3 make exclusive use of NDVI. Although, the use of NDVI in 138  those studies may have been imposed by technological constraints because, the AVHRR sensor lacked bands in the longer spectral wavelengths. However, the longer wavelengths represented in Landsat TM and ETM+ data allow for more varied use of spectral vegetation indexes, as shown in Chapter3. The Normalized Difference Vegetation Index, Normalized Burn Ratio, Brightness, Greenness, Wetness, Tasseled Cap Angle, Disturbance Index, Integrated Forestry Z-scores (exclusively used by the VCT algorithm), and Forest Recovery Index (a reciprocal of the Integrated Forestry Z-scores index) have all been used singly (Table 3.2) , but rarely more than one are used to compare forest recovery.  In Chapter 5, I compared Tasseled Cap Brightness, Greenness and Wetness  spectral recovery trajectories. Results showed that the use of multiple spectral indices can provide unique insights into forest recovery by informing on distinct recovery processes occurring in different sections of the Boreal Shield ecozone. This adds new knowledge to how spectral recovery and forest recovery can relate, specifically how the recovery of multiple spectral indices may indicate different complementary things about a forest recovery. In addition, Chapter 4 showed that use of multiple spectral vegetation indices added more explanatory power to forest attribute prediction models, which is important in the unmanaged boreal forest context, where little information exists describing the forests.  New methods were applied to observe spectral forest recovery through the course of this disertation. Firstly, building on work by Kennedy et al., (2010) and Pflugmacher et al., (2012), new metrics describing spectral recovery and the whole spectral trajectory were extracted from temporally segmented time series in Chapter 4. The new metrics I developed for that research quantified both annual change in spectral indices, and the generalized shape of a given spectral 139  time series. Both of these new metrics were shown to have a moderate amount of importance in forest attribute estimation models.    When traditional statistical techniques lacked the ability to determine meaningful difference between two datasets, I utilized a method that had previously been used in biological research, but is novel to remote sensing research: the Cohen’s d effect size statistic. This effect size metric allowed my research to make statements about similarity or difference between two datasets, regardless of sample size. Without using this type of metric, many other traditional statistical techniques would have declared even small differences between datasets as significant due to the very large sample sizes used in Chapter 5. Thus, the innovative use of Cohen’s d allowed me to compare spectral recovery to gain meaningful information about forest recovery.  Lastly, a new spectral time series metric was created to describe spectral recovery rates: the Year to Year Average metric. I used this metric in Chapter 6 of this dissertation to quantify annual spectral recovery over a five year moving window throughout a time series. I created this metric to capture the pure spectral recovery, without reference to the disturbance magnitudes, like that of the Recovery Indicator, or pre disturbance conditions, like that of the Ratio of Pre-Disturbance 80% metric used in Chapter 6. I also modified the Ratio of Pre-Disturbance 80% metric for the recovery rate calculation in Chapter 6, importantly setting the five year recovery window as the critical time period to monitor boreal forest reestablishment after disturbance.  7.3 Application of Spectral Forest Recovery for Data Generation in Boreal Forests Non-inventoried boreal forests lack detailed spatial data describing their condition. Most often these forests remain non-inventoried due to their remoteness, lack of merchantable timber, and 140  cost of establishing and carrying out in situ forest attribute mensuration. The work presented in Chapters 4, 5, and 6 of this dissertation provide new data about non-inventoried and inventoried Canadian boreal forests. Chapter 4 demonstrated the application of spectral time series can reliably estimate forest attributes. Chapter 5 displayed the utility of monitoring forest recovery with spectral time series, while Chapter 6 applied spectral time series to determine forest recovery rates over 20 years.   Spectral time series metrics quantify periods of no change, disturbance, and recovery using concepts such as magnitude and duration, which is critical to understanding the change that had occurred in each particular segment. Chapter 4 explored the use of temporal trajectory metrics in an unmanaged boreal forest context. A highlight from that research showed that temporal trajectory metrics are much better for estimating forest attributes than single image date techniques showing a nearly 30% increase in r2 values. Previous to the effort in Chapter 4, the only spatially comprehensive information that existed describing the unmanaged forests in that location had been a land cover dataset, which was also created using Landsat data and was nearly 15 years old at the time of this study.  Forest recovery has typically been tracked using techniques which compare classifications of forest recovery, however the research presented in Chapter 5 utilized spectral recovery trajectories to show forest recovery over time. This application of spectral recovery trajectories is unique in that it collapsed and averaged all recovery trajectories through the time series into one annual average recovery trajectory, providing the ability to quantify and track forest recovery from stand initiation through later successive stages of growth using spectral datasets. I demonstrated in Chapter 5 that spectral recovery can be calculated for a large area, regardless of 141  inventory status, and was able to quantify the average annual post disturbed spectral recovery for East and West Boreal Shield. This allowed me to compare nearly 30 years of expected spectral recovery across two different regions while utilizing multiple spectral vegetation indices to gain new insights into forest recovery across the regions.  Lastly, I highlighted the application of spectral forest recovery rates across a large area. The research performed in Chapter 6 is wholly unique, as the application of spectral time series has allowed me to take a census of the Boreal and Taiga Shield ecozones using repeatable and consistent methods. This is especially important where inventory data might not exist and across provincial and management jurisdictions. By applying the knowledge about spectral forest recovery gained in Chapter 5, I was able to calculate five year recovery rates based on consecutive non-overlapping epochs starting from the year of disturbance, which were then used to compare forest recovery rates over time and between ecozones. This application of spectral time series aided in determining how forests are recovering from disturbance through time, and also enabled the detection of trends across time indicating broad environmental change had been occurring.  7.4 Insights into Canadian Boreal Forest Dynamics  Up to half of the Canadian boreal lacks detailed spatial data describing the condition of forests. The research described in Chapters 4, 5, and 6 generated valuable information about boreal forests. Forest attribute estimation using time series generated metrics results demonstrated that simpler metrics summarizing the entire time series were more important than complex metrics that described one temporal segment in detail. The importance of those simple metrics lends 142  support to the idea that spectral time series metrics do provide information that would be unavailable with other datasets. In addition, metrics describing the recovery of disturbed forest stands were important, which formed a key insight: knowledge of past conditions is as important as current conditions for estimating forest attributes.   The research presented in Chapter 5 indicated that spectral recovery can provide key insights into forest recovery across large areas. The results in that chapter added further justification to the division of the Boreal Shield ecozone by displaying differences in forest recovery between different regions. Moreover, the differing patterns between the East/West division observed in the spectral datasets mirrored physical forest recovery differences from smaller areas within each ecozone. The interpretations of the spectral time series indicated that the West section experienced a coniferous dominated recovery, while the East section underwent an initial state dominated by broadleaf vegetation then succeeded by a coniferous recovery. Furthermore and previously unreported, the area undergoing faster than average recovery increased disproportionally to other recovery rate categories. This is a valuable insight into forest recovery, as it may not have been observable across such large areas without the use of remotely sensed data sets. Lastly, annual areal disturbance calculations highlighted differences between both sections reflected the already known spatial disparity in disturbance regimes within both East and West ecozone sections.  The disturbance and recovery cycles of Canadian boreal forests are generally well understood, with fire as the dominant change agent across the landscape. Equally important as mapping fire occurrence is the subsequent recovery of forest vegetation. As these cycles have been working synchronously for millennia, a changing climate exerts new pressures on both disturbance and 143  recovery processes. Forest recovery rates were examined in Chapter 6 using spectral time series. This examination of spectral recovery rates over time provided new insights into how forest recovery is being affected by a changing climate. I was able to show that the Taiga Shield ecozone has experienced an increase in spectral forest recovery rates over a 26 year period. This result indicates that boreal forests are recovering more rapidly now than they were 26 years ago, and is also suggestive of a change in annual productivity. In addition, spectral forest recovery rates were found to be fixed across neither space nor time, with many regions showing highly variable rates, but with no readily discernable trend. Interestingly, separate East and West sections for the Taiga and Boreal Shield ecozone showed different spatial patterns of spectral recovery rate trends. Western sections displayed a dispersed pattern, while eastern sections had a clustered pattern, indicating that concentrated hot-spots of change exist in the Boreal and Taiga Shield regions.   7.5 Limitations The limitations of the research presented in this dissertation mostly revolve around temporal segmentation, and the use of remotely sensed imagery time series to detect forest recovery. Although forest disturbances are well characterized with remotely sensed time series, research into using the same datasets to characterize subsequent forest recovery is still ongoing. As such, uncertainties around 1) the use of one spectral index over another, 2) the time step to observe change, or  3) the length of the time series still remain.  The use of remotely sensed time series to characterize forest disturbance and recovery is sometimes undertaken due to the prohibitive costs of other more direct methods. The low cost, 144  seamless coverage, and frequent revisit cycle are all unique strengths of the Landsat program and the temporally deep remotely sensed imagery archive that was utilized in this dissertation. As easily as the imagery archive lends itself to research inaccessible areas, those same areas typically suffer from a lack of detailed data describing current conditions. In fact, if any other data exists describing those remote places it was likely generated using Landsat data or another remotely-sensed dataset. Data on the cover and structure of vegetation would prove very useful for validating remote sensing based boreal studies where no in situ data exists. The lack of in situ measured forest attribute data can therefore hinder research relating spectral conditions to physical recovery conditions over time and space. Moreover, establishing, maintaining, and collecting in situ data from plots in remote areas is critical for forest monitoring, remote sensing based or not. The missed opportunity to collect that type of data now becomes a lack of information that can be perpetuated far into the future, and hamper future monitoring activities, be they in situ or remote sensing based.  Caution is urged when interpreting spectral recovery trajectories, as they may not always represent a definition of forest recovery acceptable to all users of forest information (Kennedy et al., 2010). Additionally, the definition of forest recovery, and a “recovering” or “recovered” forest can change dependent upon who is observing the forest and for what purpose (Bartels et al., 2016). The research in this dissertation did not create a spectral threshold to define a “recovered forest” state, though that threshold has been created and used in other studies by correlating spectral values to a canopy cover thresholds (Kennedy et al., 2012).  A number of factors need to be carefully considered when temporally segmenting a spectral time series. Chapter 5 showed that the selection of an appropriate spectral vegetation index is one 145  factor that merits careful consideration. Spectral vegetation indices that rely on shorter wavelengths like that of NDVI, or Tasseled Cap Greenness may saturate very quickly after disturbance. In contrast, indices that use longer wavelengths like Tasseled Cap Wetness and NBR show less sensitivity to initial post disturbance conditions. Thus the selection of which spectral vegetation index to use should be guided by what phenomenon is most important to observe.  The length of time examined or number and spacing of observations within a time span is a final critical factor that affects spectral time series. The longest possible time series of remotely sensed observations is currently nearly 45 years long, as represented by the Landsat imagery archive. However, some observations over this time span may be limited to a single annual or biennial occurrence, which may not represent ideal change detection conditions to produce usable data. Additionally, the Landsat imagery archive may include over forty years of data, but most Landsat time series studies do not include the previous MSS era data or the current OLI data (Vogelmann et al., 2016). Data from the MSS sensors are typically unused due to the low degree of imagery preprocessing (Pflugmacher et al., 2012), and poor geolocation (Gillanders et al., 2008) and limited spectral bands available; MSS data does not contain the longer wavelengths that are better indicators of forest recovery. Thus Landsat time series are typically limited to the TM and ETM+ eras, with data first available in 1984 and regularly available through 2016. Thus the lack of viable data before 1984 hampers efforts to report on change throughout the entire 40 plus year Landsat record.  146  7.6 Future Research & New Questions Landsat time series continually prove their ability to provide new insight into landscape dynamics. The research presented in this dissertation demonstrates how valuable the Landsat imagery archive can be when used for forest change detection. The use of Landsat time series as presented in this dissertation becomes especially important when little or no information exists for remote forests, and especially for unmanaged and non-inventoried Canadian boreal forests.  The methods used to perform the research are easily extendable and offer new ways of characterizing forest recovery spectrally. As explored in Chapter 6, spectral forest recovery rates offer a view into how forest recovery can change over time. This type of research could be used to verify reforestation claims made under international reporting requirements or further help inform carbon budgeting models.  I leave the reader with three directions to further the use of Landsat time series in analyzing forested landscapes change over time. These new directions incorporate Landsat time series to generate information about boreal forest conditions through time and can provide novel insight into how the pan-boreal region may be changing.  1. Backwards extension of Landsat time series into the MSS era and inclusion of OLI data through the use of inter-sensor calibration and a common spectral vegetation index.  2. Extension of these approaches across the entire pan-boreal to monitor worldwide boreal forest health.  3. Incorporation of spectrally assessed forest recovery into carbon accounting models.  147  The first new direction aims to integrate the older MSS era and newer OLI data into the current TM and ETM+ time series analysis paradigm. Methods for incorporating those data must be found, as a longer time period will reveal information that was previously unknown; these data can report on an earlier time period that remains unexamined. Forest recovery processes occur over long time periods, and a longer spectral time series may reveal new information beyond the nearly 30 year time period reported on in Chapters 5 and 6.  The second new direction seeks to apply the use of spectral time series across boreal forests worldwide, in service of creating a better understanding of the regional response to a changing climate. As previously noted, the initial post-disturbance forest regrowth period is a critical time to monitor boreal forests for signs of degradation. Monitoring pan-boreal conditions through the Landsat record could provide many agencies and users with information that was not previously accessible, and is easily comparable, as opposed the locationally piecemeal fashion in which many studies are presently conducted. For example, the response of Canadian and Russian boreal forests to a changing climate could be better assessed, with specific forest recovery information providing a improved insight into forest extent.  The third new direction would demonstrate that the incorporation of spectral forest recovery information into carbon accounting models can lead to more accurate information. In the absence of detailed spatial information, averages are often drawn from small samples over very large areas. The use of spectral forest recovery data at a fine spatial scale in carbon accounting models would provide a more precise information about recently disturbed areas and allow for a more accurate enumeration of carbon flux in boreal forests over time.  148  References  Adams, P., Babooram, A., Dewar, H., St. Lawrence, J., Wang, J. (2008) Climate change in Canada. Statistics Canada. Human activity and the environment: annual statistics 2007 and 2008. Ottawa. 16-201-X, 1-159.  Alcaraz-Segura, D., Chuvieco, E., Epstein, H. E., Kasischke, E. S., & Trishchenko, A. (2010). Debating the greening vs. browning of the North American boreal forest: differences between satellite datasets. Global Change Biology, 16(2), 760-770.  Amiro, B. D., Barr, A. G., Barr, J. G., Black, T. A., Bracho, R., Brown, M., Chen, J. Clark, K. L. , Davis, K. J.,  Desai, A. R. , Dore, S., Engel, V., Fuentes, J. D. , Goldstein, A. H. , Goulden, M. L. , Kolb, T. E.,  Lavigne, M. B.,  Law, B. E.,  Margolis, H. A.,  Martin, T., McCaughey, J. H.,  Mission, L., Montes-Helu, M., Noorments, A., Randerson, J. T.,  Starr, G., & Xiao, J. (2010). Ecosystem carbon dioxide fluxes after disturbance in forests of North America. Journal of Geophysical Research, 115(G4), 1-13.  Andrew, M. E., Wulder, M. A., & Coops, N. C. (2012). Identification of de facto protected areas in boreal Canada. Biological Conservation, 146(1), 97-107.  Angert, A., Biraud, S., Bonfils, C., Henning, C. C., Buermann, W., Pinzon, J., Tucker, C. J., & Fung, I. (2005). Drier summers cancel out the CO2 uptake enhancement induced by warmer springs. Proceedings of the National Academy of Sciences of the United States of America, 102(31), 10823-10827.  Banskota, A., Kayastha, N., Falkowski, M. J., Wulder, M. A., Froese, R. E., & White, J. C. (2014). Forest monitoring using landsat time series data: a review. Canadian Journal of Remote Sensing, 40(5), 362-384.  Barber, V. A., Juday, G. P., & Finney, B. P. (2000). Reduced growth of Alaskan white spruce in the twentieth century from temperature-induced drought stress. Nature, 405(6787), 668-673.  Bartels, S. F., Chen, H. Y., Wulder, M. A., & White, J. C. (2016). Trends in post-disturbance recovery rates of Canada’s forests following wildfire and harvest. Forest Ecology and Management, 361, 194-207.  Bater, C. W., Wulder, M. A., Coops, N. C., Hopkinson, C., Coggins, S. B.,  Arsenault, E., Beaudoin A., Guindon, L., Hall, R. J. , Villemaire, P., & Woods, M. (2011). Model development for the estimation of aboveground biomass using a lidar based sample of Canada’s boreal fire. Silvilaser 2011. Oct 16-19, 2011. Hobart, Tasmania, Australia.   Baumann, M., Ozdogan, M., Kuemmerle, T., Wendland, K. J., Esipova, E., & Radeloff, V. C. (2012). Using the Landsat record to detect forest-cover changes during and after the 149  collapse of the Soviet Union in the temperate zone of European Russia. Remote Sensing of Environment, 124, 174-184.  Beck, P. S., & Goetz, S. J. (2011). Satellite observations of high northern latitude vegetation productivity changes between 1982 and 2008: ecological variability and regional differences. Environmental Research Letters, 6(4), 045501.  Bergeron, Y., Drapeau, P., Gauthier, S., & Lecomte, N. (2007). Using knowledge of natural disturbances to support sustainable forest management in the northern Clay Belt. The Forestry Chronicle, 83(3), 326-337.  Bergeron, Y., Flannigan, M., Gauthier, S., Leduc, A., & Lefort, P. (2004). Past, current and future fire frequency in the Canadian boreal forest: implications for sustainable forest management. AMBIO: A Journal of the Human Environment,33(6), 356-360.  Bergeron, Y., Gauthier, S., Kafka, V., Lefort, P., & Lesieur, D. (2001). Natural fire frequency for the eastern Canadian boreal forest: consequences for sustainable forestry. Canadian Journal of Forest Research, 31(3), 384-391.  Bergeron, Y. (2000). Species and stand dynamics in the mixed woods of Quebec's southern boreal forest. Ecology, 81(6), 1500-1516.  Berner, L. T., Beck, P. S., Bunn, A. G., Lloyd, A. H., & Goetz, S. J. (2011). Highlatitude tree growth and satellite vegetation indices: Correlations and trends in Russia and Canada (1982–2008). Journal of Geophysical Research: Biogeosciences (2005–2012), 116(G1).  Bi, J., Xu, L., Samanta, A., Zhu, Z., & Myneni, R. (2013). Divergent arctic-boreal vegetation changes between North America and Eurasia over the past 30 years. Remote Sensing, 5(5), 2093-2112.  Blaishke, T. (2010). Object based image analysis for remote sensing. Photogrammetry and Remote Sensing. 65(1), 2-16.  Boisvenue, C., & Running, S. W. (2006). Impacts of climate change on natural forest productivity–evidence since the middle of the 20th century. Global Change Biology, 12(5), 862-882.  Bolton, D. K., Coops, N. C., & Wulder, M. A. (2015). Characterizing residual structure and forest recovery following high-severity fire in the western boreal of Canada using Landsat time-series and airborne lidar data. Remote Sensing of Environment, 163, 48-60.  Bolton, D. K., Coops, N. C., & Wulder, M. A. (2013). Measuring forest structure along productivity gradients in the Canadian boreal with small-footprint Lidar. Environmental Monitoring and Assessment, 185(8), 6617-6634. 150    Bonan, G. B. (2008). Forests and climate change: forcings, feedbacks, and the climate benefits of forests. science, 320(5882), 1444-1449.  Bonan, G. B., & Shugart, H. H. (1989). Environmental factors and ecological processes in boreal forests. Annual Review of Ecology and Systematics, 1-28.  Bond-Lamberty, B., Peckham, S. D.,  Ahl, D. E.,  & Gower, S.T. (2007). Fire as the dominant driver of central Canadian boreal forest carbon balance. Nature. 450, 89-92.   Brandt, J. P., Flannigan, M. D., Maynard, D. G., Thompson, I. D., & Volney, W. J. A. (2013). An introduction to Canada’s boreal zone: ecosystem processes, health, sustainability, and environmental issues 1. Environmental Reviews, 21(4), 207-226.  Brandt, J. P. (2009). The extent of the North American boreal zone. Environmental Reviews, 17(NA), 101-161.  Brassard, B. W.,  & Chen, Y.H. (2008). Effects of Forest Type and Disturbance on Diversity of Coarse Woody Debris in Boreal Forest. Ecosystems. 11(7), 1078-1090.  Brassard, B. W., & Chen, H. Y. (2006). Stand structural dynamics of North American boreal forests. Critical Reviews in Plant Sciences, 25(2), 115-137.  Breiman, L. (2001). Random forests. Machine learning, 45(1), 5-32.  Breiman, L., Friedman, J., Stone, C. J., & Olshen, R. A. (1984). Classification and regression trees. Chapman & Hall/CRC.  Brooks, J. R., Flanagan, L. B., & Ehleringer, J. R. (1998). Responses of boreal conifers to climate fluctuations: indications from tree-ring widths and carbon isotope analyses. Canadian Journal of Forest Research, 28(4), 524-533.  Buma, B. (2012). Evaluating the utility and seasonality of NDVI values for assessing post-disturbance recovery in a subalpine forest. Environmental Monitoring and Assessment, 184(6), 3849-3860.  Bunn, A. G., & Goetz, S. J. (2006). Trends in satellite-observed circumpolar photosynthetic activity from 1982 to 2003: the influence of seasonality, cover type, and vegetation density. Earth Interactions, 10(12), 1-19.  Bunn A. G., Goetz, S. J., Fiske, G. J. (2005) Observed and predicted responses of plant growth to climate across Canada. Geophysical Research Letters. 32:L16710  151  Canty, M. J., Nielsen, A. A., & Schmidt, M. (2004). Automatic radiometric normalization of multitemporal satellite imagery. Remote Sensing of Environment. 91(3-4), 441−451.  Cavin, L., Mountford, E. P., Peterken, G. F., & Jump, A. S. (2013). Extreme drought alters competitive dominance within and between tree species in a mixed forest stand. Functional Ecology, 27(6), 1424-1435.  Chen, H. Y., & Popadiouk, R. V. (2002). Dynamics of North American boreal mixedwoods. Environmental Reviews, 10(3), 137-166.  Chen, X., Vogelmann, J. E., Rollins, M., Ohlen, D., Key, C. H., Yang, L., Huang, C.,  & Shi, H. (2011). Detecting post-fire burn severity and vegetation recovery using multitemporal remote sensing spectral indices and field-collected composite burn index data in a ponderosa pine forest. International Journal of Remote Sensing, 32(23), 7905-7927.  Chen, W. J., Black, T. A., Yang, P. C., Barr, A. G., Neumann, H. H.,  Nesic, Z., Blanken, P. D., Novak, M. D.,  Eley, J. Ketler, R. J.,  & Cuenca, R. (1999) Effect of climatic variability on the annual carbon sequestration by a Boreal Aspen forest. Global Change Biology, 5(1), 41–53.  Chu, T., & Guo, X. (2013). Remote sensing techniques in monitoring post-fire effects and patterns of forest recovery in boreal forest regions: a review. Remote Sensing, 6(1), 470-520.  Cohen, W. B., & Goward, S. N. (2004). Landsat's role in ecological applications of remote sensing. Bioscience, 54(6), 535–545.  Cohen, W. B., Fiorella, M., Gray, J., Helmer, E., & Anderson, K. (1998). An efficient and accurate method for mapping forest clearcuts in the Pacific Northwest using Landsat imagery. Photogrammetric Engineering and Remote Sensing, 64(4), 293-299.  Cohen, W. B., Spies, T. A., & Fiorella, M. (1995). Estimating the age and structure of forests in a multi-ownership landscape of western Oregon, USA. International Journal of Remote Sensing, 16(4), 721-746.   Cohen, W. B., & Spies, T. A. (1992). Estimating structural attributes of Douglas-fir/western hemlock forest stands from Landsat and SPOT imagery. Remote Sensing of Environment, 41(1), 1-17.  Cohen, J. (1988) Statistical power analysis for the behavioral sciences. 2nd Ed. Hove: Lawrence Erlbaum Associates.  Collins, J. B., & Woodcock, C. E. (1996). An assessment of several linear change detection techniques for mapping forest mortality using multitemporal Landsat TM data. Remote Sensing of Environment, 56(1), 66−77. 152   Colombo, S. J. (1998). Climatic warming and its effect on bud burst and risk of frost damage to white spruce in Canada. The Forestry Chronicle, 74(4), 567-577.  Crist, E. P., & Cicone, R. C. (1984a). A physically-based transformation of Thematic Mapper data---The TM Tasseled Cap. IEEE Transactions on Geoscience and Remote Sensing, IEEE Transactions On, (3), 256-263.  Crist, E. P., & Cicone, R. C. (1984b). Application of the tasseled cap concept to simulated thematic mapper data. Photogrammetric Engineering and Remote Sensing, 50(3), 327 – 331.  Cuevas-González, M., Gerard, F., Balzter, H., & Riano, D. (2009). Analysing forest recovery after wildfire disturbance in boreal Siberia using remotely sensed vegetation indices. Global Change Biology, 15(3), 561-577.  Cutler, R. D., Edwards, T. C., Jr., Beard, K. H., Cutler, A., Hess, K. T., Gibson, J., & Lawler, J. J. (2007). Random Forests for classification in ecology. Ecology, 88(11), 2783−2792.  De Groot, W. J., Landry, R., Kurz, W. A., Anderson, K. R., Englefield, P., Fraser, R. H., Hall, R. J.,  Banfield, E., Raymond, D. A.,  Decker, V., Lynham, T. J.,  Pritchard, J. M. (2007). Estimating direct carbon emissions from Canadian wildland fires. International Journal of Wildland Fire, 16(5), 593-606.  Duncanson, L. I., Niemann, K. O.,  & Wulder, M.A.. (2010). Integration of GLAS and Landsat TM data for aboveground biomass estimation. Canadian Journal of Remote Sensing, 36(2), 129-141.  Ecological Stratification Working Group. (1996). A national ecological framework for Canada. Agriculture and Agri-Food Canada and Environment Canada: Ottawa, ON. Available from: http://sis.agr.gc.ca/cansis/publications/ecostrat/cad_report.pdf [cited on May 12th 2015]   Eidenshink, J., Schwind, B., Brewer, K., Zhu, Z., Quayle, B., & Howard, S. (2007). A project for monitoring trends in burn severity. Fire Ecology 3 (1): 3-21.Fire Ecology Special Issue Vol, 3, 4.   Epting, J., Verbyla, D., & Sorbel, B. (2005). Evaluation of remotely sensed indices for assessing burn severity in interior Alaska using Landsat TM and ETM+. Remote Sensing of Environment, 96(3), 328-339.  Falkowski, M. J., Evans, J. S., Martinuzzi, S., Gessler, P. E., & Hudak, A. T. (2009). Characterizing forest succession with Lidar data: An evaluation for the inland Northwest USA. Remote Sensing of Environment, 113(1), 946−956  153  Flannigan, M. D., Logan, K. A., Amiro, B. D., Skinner, W. R., & Stocks, B. J. (2005). Future area burned in Canada. Climatic Change, 72(1-2), 1-16  Flannigan, M. D., Stocks, B. J., & Wotton, B. M. (2000). Climate change and forest fires. Science of the total environment, 262(3), 221-229.  Food and Agriculture Organization of the United Nations [FAO] (2001). Global Forest Resources Assessment 2000 Main Report. FAO Forestry Paper 140. Rome. 479 p. fao.org/forestry/site/24691/en/  Food and Agriculture Organization of the United Nations [FAO](2006). Global Forest Resources Assessment 2005 Progress towards sustainable forest management. FAO Forestry Paper 147. Rome.  Franklin, S. E.,  Lavigne, M. B.,  Wulder, M. A.,  & Stenhouse, G.B. (2002). Change Detection and Landscape structure mapping using remote sensing. The Forestry Chronicle. 78(5), 618-625  Franklin, S. E., Moskal, L. M., Lavigne, M. B., & Pugh, K. (2000). Interpretation and classification of partially harvested forest stands in the Fundy Model Forest using multitemporal Landsat TM. Canadian Journal of Remote Sensing, 26, 318−333.  Franklin, J. (1995). Predictive vegetation mapping: geographic modeling of biospatial patterns in relation to environmental gradients. Progress in Physical Geography. 19 (4), 474-499.  Fraser, R. H., Olthof, I., Kokelj, S. V., Lantz, T. C., Lacelle, D., Brooker, A., Wolfe, S.,  & Schwarz, S. (2014). Detecting Landscape Changes in High Latitude Environments Using Landsat Trend Analysis: 1. Visualization. Remote Sensing, 6(11), 11533-11557.  Fraser, R. H., Olthof, I., Carrière, M., Deschamps, A., & Pouliot, D. (2011). Detecting long-term changes to vegetation in northern Canada using the Landsat satellite image archive. Environmental Research Letters, 6(4), 045502.  Fazakas, Z., Nilsson, M., & Olsson, H. (1999). Regional forest biomass and wood volume estimation using satellite data and ancillary data. Agricultural and Forest Meteorology. 98-99, 417-425.  Frazier, R. J., Coops, N. C., & Wulder, M. A. (2015). Boreal Shield forest disturbance and recovery trends using Landsat time series. Remote Sensing of Environment, 170, 317-327.  Frazier, R. J., Coops, N. C., Wulder, M. A., & Kennedy, R. (2014). Characterization of aboveground biomass in an unmanaged boreal forest using Landsat temporal segmentation metrics. ISPRS Journal of Photogrammetry and Remote Sensing, 92, 137-146. 154   French, N. H., Kasischke, E. S., Hall, R. J., Murphy, K. A., Verbyla, D. L., Hoy, E. E., & Allen, J. L. (2008). Using Landsat data to assess fire and burn severity in the North American boreal forest region: an overview and summary of results. International Journal of Wildland Fire, 17(4), 443-462.  Frolking, S., Palace, M. W., Clark, D. B., Chambers, J. Q., Shugart, H. H., & Hurtt, G. C. (2009). Forest disturbance and recovery: A general review in the context of spaceborne remote sensing of impacts on aboveground biomass and canopy structure. Journal of Geophysical Research: Biogeosciences (2005–2012), 114(G2).  Gamache, I., Payette, S..(2004). Height growth response of tree line black spruce to recent warming across the forest-tundra of eastern Canada. Journal of Ecology. 92, 835-845.   Gauthier, S., Bernier, P., Kuuluvainen, T., Shvidenko, A. Z., & Schepaschenko, D. G. (2015). Boreal forest health and global change. Science, 349(6250), 819-822.  Gauthier, S., Bernier, P., Burton, P. J., Edwards, J., Isaac, K., Isabel, N., Jayen, K., Le Goff, H., & Nelson, E. A. (2014). Climate change vulnerability and adaptation in the managed Canadian boreal forest. Environmental Reviews, 22(3), 256-285.  GeoBC. (2012) Gridded DEM. http://archive.ilmb.gov.bc.ca/crgb/products/imagery/gridded.htm. Accessed 7/31/2012.  Geomatics Yukon. (2012) 30 Meter DEM. ftp://ftp.geomaticsyukon.ca/DEMs/30m. Accessed 7/31/2012  Gillanders, S. N., Coops, N. C., Wulder, M. A., Gergel, S. E., & Nelson, T. (2008). Multitemporal remote sensing of landscape dynamics and pattern change: describing natural and anthropogenic trends. Progress in Physical Geography, 32(5), 503-528.  Gillis, M. D.,  Omule, A. Y.,  & Brierley, T. (2005). Monitoring Canada’s forests: The National Forest Inventory. The Forestry Chronicle 81(2): 214-221.  Gitas, I., Polychronaki, A., Mitri, G., & Veraverbeke, S. (2012). Advances in remote sensing of post-fire vegetation recovery monitoring-a review. INTECH Open Access Publisher.  Goetz, S. J., Fiske, G. J., & Bunn, A. G. (2006). Using satellite time-series data sets to analyze fire disturbance and forest recovery across Canada. Remote Sensing of Environment, 101(3), 352-365.  Goetz, S. J., Bunn, A. G., Fiske, G. J., & Houghton, R. A. (2005). Satellite-observed photosynthetic trends across boreal North America associated with climate and fire disturbance. Proceedings of the National Academy of Sciences of the United States of America, 102(38), 13521-13525. 155   Gómez, C., White, J. C., & Wulder, M. A. (2011). Characterizing the state and processes of change in a dynamic forest environment using hierarchical spatio-temporal segmentation. Remote Sensing of Environment, 115(7), 1665-1679.  Goodwin, N. R., Coops, N. C., Wulder, M. A., Gillanders, S., Schroeder, T. A., & Nelson, T. (2008). Estimation of insect infestation dynamics using a temporal sequence of Landsat data. Remote sensing of environment, 112(9), 3680-3689.  Griffiths, P., Kuemmerle, T., Baumann, M., Radeloff, V. C., Abrudan, I. V., Lieskovsky, J., Munteanu, C., Ostapowicz, K., & Hostert, P. (2014). Forest disturbances, forest recovery, and changes in forest types across the Carpathian ecoregion from 1985 to 2010 based on Landsat image composites. Remote Sensing of Environment, 151, 72-88.  Guindon, L., Bernier, P. Y., Beaudoin, A., Pouliot, D., Villemaire, P., Hall, R. J., Latifovic, R.,  & St-Amant, R. (2014). Annual mapping of large forest disturbances across Canada’s forests using 250 m MODIS imagery from 2000 to 2011.Canadian Journal of Forest Research, 44(12), 1545-1554.  Hall, R. J., Freeburn, J. T., De Groot, W. J., Pritchard, J. M., Lynham, T. J., & Landry, R. (2008). Remote sensing of burn severity: experience from western Canada boreal fires. International Journal of Wildland Fire, 17(4), 476-489.  Hall, R. J.,  Skakun, R. S.,  Arsenault, E. J.,  Case, B.S. (2006). Modeling forest stand structure attributes using Landsat ETM+ data: Application to mapping aboveground biomass and stand volume. Forest Ecology and Management. 225(1) 378-390.  Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova, S. A., Tyukavina, A., Thau,  D., Stehman, S. V., Goetz, S. J., Loveland, T. R., Kommareddy, A., Egorov, A., Chini, L., Justice, C. O., & Townshend, J. R. G. (2013). High-resolution global maps of 21st-century forest cover change. Science, 342(6160), 850-853.  Hansen, M. C., & Loveland, T. R. (2012). A review of large area monitoring of land cover change using Landsat data. Remote sensing of Environment, 122, 66-74.  Hansen, M. J., Franklin, S. E., Woudsma, C., & Peterson, M. (2001). Forest structure classification in the North Columbia mountains using the Landsat TM Tasseled Cap Wetness component. Canadian Journal of Remote Sensing, 27(1), 20-32.  Harmon, M. E., Ferrell, W. K., & Franklin, J. F. (1990). Effects on carbon storage of conversion of old- growth forests to young forests. Science, 247(4943), 699−702.  Hay, G. J., Castilla, G., Wulder, M. A., & Ruiz, J. R. (2005). An automated object-based approach for the multiscale image segmentation of forest scenes. International Journal of Applied Earth Observation and Geoinformation, 7(4), 339-359. 156   Healey, S. P., Yang, Z., Cohen, W. B., & Pierce, D. J. (2006). Application of two regression-based methods to estimate the effects of partial harvest on forest structure using Landsat data. Remote Sensing of Environment, 101(1), 115-126.  Healey, S. P., Cohen, W. B., Zhiqiang, Y., & Krankina, O. N. (2005). Comparison of Tasseled Cap-based Landsat data structures for use in forest disturbance detection. Remote Sensing of Environment, 97(3), 301-310.  Heinselman, M. L. (1981). Fire and succession in the conifer forests of northern North America. In Forest succession (pp. 374-405). Springer New York.  Hermosilla, T., Wulder, M. A., White, J. C., Coops, N. C., & Hobart, G. W. (2015a). Regional detection, characterization, and attribution of annual forest change from 1984 to 2012 using Landsat-derived time-series metrics. Remote Sensing of Environment, 170, 121-132.  Hermosilla, T., Wulder, M. A., White, J. C., Coops, N. C., & Hobart, G. W. (2015b). An integrated Landsat time series protocol for change detection and generation of annual gap-free surface reflectance composites. Remote Sensing of Environment, 158, 220-234.  Hogg, E. H., Brandt, J. P., Kochtubajda, B. (2005) Factors affecting interannual variation in growth of western Canadian aspen forests during 1951–2000. Canadian Journal of Forest Research, 35(3), 610–622.  Horler, D.N. H.,  Ahern, F.J. (1986). Forestry information content of Thematic Mapper data. International Journal of Remote Sensing. (7) 5449-5464.  Houghton, R. A. (1999). The annual net flux of carbon to the atmosphere from changes in land use 1850-1990. Tellus, 51B, 298−313.  Huang, C., Goward, S. N., Masek, J. G., Thomas, N., Zhu, Z., & Vogelmann, J. E. (2010). An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks. Remote Sensing of Environment, 114(1), 183-198.  Huang, C., Wylie, B., Yang, L., Homer, C., & Zylstra, G. (2002). Derivation of a Tasselled Cap transformation based on Landsat 7 at-satellite reflectance. International Journal of Remote Sensing, 23(8), 1741-1748.  Hudak, A. T., Crookston, N. L., Evans, J. S., Hall, D. E., & Falkowski, M. J. (2008). Nearest neighbor imputation modeling of species-level, plot-scale structural attributes from lidar data. Remote Sensing of Environment, 112(5), 2232−2245.  157  Huete, A. R., Liu, H. Q., Batchily, K., & Van Leeuwen, W. J. D. A. (1997). A comparison of vegetation indices over a global set of TM images for EOS-MODIS. Remote Sensing of Environment, 59(3), 440-451.  Iverson, L., & Prasad, A. (1998). Predicting abundance of 80 tree species following climate change in the eastern United States. Ecological Monographs. 68(4), 465-485.  Jarvis, P., & Linder, S. (2000). Botany: constraints to growth of boreal forests. Nature, 405(6789), 904-905.  Jin, S., & Sader, S. A. (2005). Comparison of time series tasseled cap wetness and the normalized difference moisture index in detecting forest disturbances. Remote Sensing of Environment. 94(3), 364–372  Johnson, D. H. (1999). The insignificance of statistical significance testing. The Journal of Wildlife Management, 763-772.  Johnson, E. A. (1996). Fire and vegetation dynamics: studies from the North American boreal forest. Cambridge University Press.  Johnstone, J. F., Hollingsworth, T. N., Chapin, F. S., & Mack, M. C. (2010). Changes in fire regime break the legacy lock on successional trajectories in Alaskan boreal forest. Global Change Biology, 16(4), 1281-1295.  Johnstone, J. F., & Chapin III, F. S. (2006a). Effects of soil burn severity on post-fire tree recruitment in boreal forest. Ecosystems, 9(1), 14-31.  Johnstone, J. F., & Chapin III, F. S. (2006b). Fire interval effects on successional trajectory in boreal forests of northwest Canada. Ecosystems, 9(2), 268-277.  Johnstone, J.F.; Kasischke, E.S. (2005). Stand-level effects of soil burn severity on postfire regeneration in a recently-burned black spruce forest. Canadian Journal of Forest Research. 35: 2151–2163  Ju, J., & Masek, J. G. (2016). The vegetation greenness trend in Canada and US Alaska from 1984–2012 Landsat data. Remote Sensing of Environment, 176, 1-16.  Kasischke, E. S., Loboda, T., Giglio, L., French, N. H., Hoy, E. E., de Jong, B., & Riano, D. (2011). Quantifying burned area for North American forests: Implications for direct reduction of carbon stocks. Journal of Geophysical Research: Biogeosciences (2005–2012), 116(G4).  Kasischke, E. S., & Turetsky, M. R. (2006). Recent changes in the fire regime across the North American boreal region—spatial and temporal patterns of burning across Canada and Alaska. Geophysical Research Letters, 33(9). 158   Kennedy, R. E., Andréfouët, S, Cohen, W. B., Gomez, C., Griffiths, P., Hais, M., Healey, S. P., Helmer, E. H., Hostert, P., Lyons, M. B., Meigs, G. W., Pflugmacher, D., Phinn, S. R., Powell, S. L., Scarth, P., Sen, S., Schroeder, T. A., Schneider, A., Sonnenschein, R., Vogelmann, J. E., Wulder, M. A.; Zhu, Z. (2014). Bringing an ecological view of change to Landsat based remote sensing. Frontiers in Ecology and the Environment, 12(6), 339-346.  Kennedy, R. E., Yang, Z., Cohen, W. B., Pfaff, E., Braaten, J., & Nelson, P. (2012). Spatial and temporal patterns of forest disturbance and regrowth within the area of the Northwest Forest Plan. Remote Sensing of Environment, 122, 117-133.  Kennedy, R. E., Yang, Z., & Cohen, W. B. (2010). Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sensing of Environment, 114(12), 2897-2910.   Kennedy, R. E., Townsend, P. A., Gross, J. E., Cohen, W. B., Bolstad, P., Wang, Y. Q., & Adams, P. (2009). Remote sensing change detection tools for natural resource managers: Understanding concepts and tradeoffs in the design of landscape monitoring projects. Remote Sensing of Environment. 113(7), 1382–1396  Kennedy, R.E, Cohen, W. B.,  Schroeder, T.A. (2007). Trajectory based change detection for automated characterization of forest disturbance dynamics. Remote Sensing of Environment. 110(3), 370-386.  Key, C. H. (2006). Ecological and sampling constraints on defining landscape fire severity. Fire Ecology, 2(2), 34-59.  Key, C. H., & Benson, N. C. (1999). The Normalized Burn Ratio (NBR): A Landsat TM radiometric measure of burn severity. United States Geological Survey.  Kull, S. J., Kurz, W. A., Rampley, G. J., Banfield, G. E., Schivatcheva, R. K., & Apps, M. J. (2006). Operational-scale carbon budget model of the Canadian forest sector (CBM-CFS3) version 1.0: user's guide. 2006.  Kurz, W. A., Shaw, C. H., Boisvenue, C., Stinson, G., Metsaranta, J., Leckie, D., ... & Neilson, E. T. (2013). Carbon in Canada’s boreal forest—A synthesis 1. Environmental Reviews, 21(4), 260-292.  Kurz, W. A. (2010). An ecosystem context for global gross forest cover loss estimates. Proceedings of the National Academy of Sciences, 107(20), 9025-9026.  Kurz, W.A.; Apps, M.J.; Webb, T.M.; McNamee, P.J. (1992).   The carbon budget of the Canadian forest sector: Phase   I. For. Can., Northwest Reg., Edmonton, AB. Inf. Rep.   NOR-X-326 159   Labrecque, S., Fournier, R. A.,  Luther, J. E.,  & Piercey, D.(2006) A comparison of four methods to map biomass from Landsat-TM and inventory data in western Newfoundland. Forest Ecology and Management, 226(1-3), 129-144.  LeComte, N.,  Bergeron, Y. (2005). Successional pathways on different surficial deposits in the coniferous boreal forest of the Quebec Clay Belt. Canadian Journal of Forest Research, 2005, vol. 35, no 8, p. 1984-1995.  Lemprière, T. C.,  Bernier, P. Y., Carroll, A. L., Flannigan, M. D., Gilsenan, R. P., McKenney, D. W., Hogg, E. H., Pedlar,  J. H.,  & Blain, D. (2008). The importance of forest sector adaptation to climate change (Vol. 416). Northern Forestry Centre.  Lentile, L. B., Holden, Z. A., Smith, A. M., Falkowski, M. J., Hudak, A. T., Morgan, P., Lewis, S. A., Gessler,  P. E.,  & Benson, N. C. (2006). Remote sensing techniques to assess active fire characteristics and post-fire effects. International Journal of Wildland Fire, 15(3), 319-345.  Lloyd, A. H., Bunn, A. G.,  (2007) Responses of the circumpolar boreal forest to the 20th century climate variability. Environmental Research Letters. (2), 1-13  Loboda, T. V., Zhang, Z., O'Neal, K. J., Sun, G., Csiszar, I. A., Shugart, H. H., & Sherman, N. J. (2012). Reconstructing disturbance history using satellite-based assessment of the distribution of land cover in the Russian Far East. Remote Sensing of Environment, 118, 241-248.  Luckman, B. H., Youngblut, D. K., Wilson R. J. S., Watson, E., Colenutt, M.E.,  & St. George, R. S. (2004). Dendroclimatology in the Canadian Cordillera. In Meeting on Tree Rings and Climate: Sharpening the Focus. Marriott University Park, Tucson, AZ.  Markham, B. L. and Helder, D. L. (2012). Forty-year calibrated record of earth-reflected radiance from Landsat: A review. Remote Sensing of Environment. 122, 30-40.  Masek, J. G., Goward, S. N., Kennedy, R. E., Cohen, W. B., Moisen, G. G., Schleeweis, K., & Huang, C. (2013). United States forest disturbance trends observed using Landsat time series. Ecosystems, 16(6), 1087-1104.  Masek, J. G., Huang, C., Wolfe, R., Cohen, W., Hall, F., Kutler, J., & Nelson, P. (2008). North American forest disturbance mapped from a decadal Landsat record. Remote Sensing of Environment, 112(6), 2914-2926.  Masek, J. G.,  E.F. Vermote, N. Saleous, R. Wolfe, F.G. Hall, Huemmrich, K. F.,   Gao, J. Kutler, and T.K. Lim, (2006), A Landsat surface reflectance data set for North America, 1990-2000, Geoscience and Remote Sensing Letters, 3, 68-72.  160  Masek, J. G. (2001). Stability of boreal forest stands during recent climate change: evidence from Landsat satellite imagery. Journal of Biogeography, 28(8), 967-976.  McManus, K. M., Morton, D. C., Masek, J. G., Wang, D., Sexton, J. O., Nagol, J. R., Ropars, P., & Boudreau, S. (2012). Satellitebased evidence for shrub and graminoid tundra expansion in northern Quebec from 1986 to 2010. Global Change Biology, 18(7), 2313-2323.  Mehrabi, Z., Slade, E. M., Solis, A., & Mann, D. J. (2014). The importance of microhabitat for biodiversity sampling. PloS one, 9(12), e114015.  Meidinger, D., & Pojar, J. (1991). Ecosystems of British Columbia. Victoria, B.C.: B.C. Ministry of Forests.  Myneni, R. B.,  Keeling, C. D.,  Tucker, C. J.,  Asrar, G., Memani, R.R. (1997) Increased Plant growth in the northern high latitudes from 1981-1991. Nature. 386, 698-701  Næsset, E., Gobakken T., & Nelson, R. (2009). Sampling and mapping forest volume and biomass using airborne LIDARs. In Proceedings of the eight annual forest inventory and analysis symposium, October 16-19, 2006, Monterey, CA, USA. General Technical Report WO-79, United States Department of Agriculture, Forest Service, Washington D.C. 297-301.  Nakagawa, S., & Cuthill, I. C. (2007). Effect size, confidence interval and statistical significance: a practical guide for biologists. Biological Reviews,82(4), 591-605.  Nemani, R. R., Keeling, C. D., Hashimoto, H., Jolly, W. M., Piper, S. C., Tucker, C. J., Myneni, R.,  & Running, S. W. (2003). Climate-driven increases in global terrestrial net primary production from 1982 to 1999. Science, 300(5625), 1560-1563  Oliver, C. D., & Larson, B. C. (1990). Forest stand dynamics. McGraw-Hill, Inc..   Olsson, H. (2009). A method for using Landsat time series for monitoring young plantations in boreal forests. International Journal of Remote Sensing, 30(19), 5117-5131.  Olthof, I., & Latifovic, R. (2007). Short term response of arctic vegetation NDVI to temperature anomalies. International Journal of Remote Sensing, 28(21), 4823-4840.  Pacala, S., Birdsey, R. A., Bridgham, S. D., Conant, R. T., Davis, K., & Hales, B. (2007). The North American carbon budget past and present. In A. W. King (Ed.), The first State of the Carbon Cycle Report (SOCCR), Asheville, N.C.: NOAA National Climate Data Center.  Pan, Y., Birdsey, R. A., Fang, J., Houghton, R., Kauppi, P. E., Kurz, W. A., Phillips, O. L., Shvidenko, A., Lewis, S. L., Canadell, J. G., Ciais, P., Jackson, R. B., Pacala, S. W., 161  McGuire, A. D., Piao, S., Rautianen, A., Sitch, S., & Hayes, D. (2011). A large and persistent carbon sink in the world’s forests. Science, 333(6045), 988-993.  Paradella, W. R.,  Da Silva, M.F. F.,  De N., Rosa, A. & Kushigbor, C.A. (1994). A geobotanical approach to the tropical rain forest environment of the Carajás mineral province Amazon region, Brazil, based on digital TM-Landsat and DEM data. International Journal of Remote Sensing, 15(8), 1633-1648.  Payette, S. (1992). Fire as a controlling process in the North American Boreal Forest. In Shugart, H. H.,  Leemans, R., Bonan, G. B.,  Eds. A systems analysis of the global boreal forests. Cambridge: Cambridge University press, 144-169.  Pflugmacher, D., Cohen, W.B.,  & Kennedy, R.E. (2012). Using Landsat derived disturbance history (1972-2012) to predict current forest structure. Remote Sensing of Environment, 122, 146-165.  Pickell, P. D., Hermosilla, T., Frazier, R. J.,  Coops, N. C., & Wulder, M. A. (2016). Forest recovery trends derived from Landsat time series for North American boreal forests. International Journal of Remote Sensing, 37(1), 138-149.  Potapov, P. V., Turubanova, S. A., Tyukavina, A., Krylov, A. M., McCarty, J. L., Radeloff, V. C., & Hansen, M. C. (2014). Eastern Europe's forest cover dynamics from 1985 to 2012 quantified from the full Landsat archive. Remote Sensing of Environment. 159 (2015): 28-43.  Potapov, P., Hansen, M. C., Stehman, S. V., Loveland, T. R., & Pittman, K. (2008). Combining MODIS and Landsat imagery to estimate and map boreal forest cover loss. Remote Sensing of Environment, 112(9), 3708-3719.  Pouliot, D., Latifovic, R., & Olthof, I. (2009). Trends in vegetation NDVI from 1 km AVHRR data over Canada for the period 1985–2006. International Journal of Remote Sensing, 30(1), 149-168.  Powell, S. L., Cohen, W. B., Healey, S. P., Kennedy, R. E., Moisen, G. G., Pierce, K. B., & Ohmann, J. L. (2010). Quantification of live aboveground forest biomass dynamics with Landsat time series and field inventory data: A comparison of empirical modeling approaches. Remote Sensing of Environment, 114, 1053−1068.  Powers, R. P., Hay, G. J., & Chen, G. (2012). How wetland type and area differ through scale: A GEOBIA case study in Alberta's Boreal Plains. Remote Sensing of Environment, 117, 135-145.  Powers, R. P., Hermosilla, T., Coops, N. C., & Chen, G. (2015). Remote sensing and object-based techniques for mapping fine-scale industrial disturbances. International Journal of Applied Earth Observation and Geoinformation, 34, 51-57. 162   Prasad, A. M., Iverson, L. R., & Liaw, A. (2006). Newer classification and regression tree techniques: Bagging and random forests for ecological prediction. Ecosystems, 9, 181−199.  Price, D. T., Alfaro, R. I., Brown, K. J., Flannigan, M. D., Fleming, R. A., Hogg, E. H., Girardin, M. P.,  Lakusta, T., Johnston, M., McKenney, D. W.,  Pedlar, J., Stratton, T., Sturrock, R., Thompson, I., Trofymow, JA., Venier, L. A., . (2013). Anticipating the consequences of climate change for Canada’s boreal forest ecosystems 1. Environmental Reviews, 21(4), 322-365.  Richetti, E. (2000). Multispectral satellite image and ancillary data integration for geological classification. Photogrammetric Engineering and Remote Sensing, 66(4), 429-435.  Rogan, J., Miller, J., Stow, D., Franklin, J., Levien, L., & Fisher, C. (2003). Land-Cover Change Monitoring with Classification Trees Using Landsat TM and Ancillary Data. Photogrammetric Engineering & Remote Sensing, 69(7), 793-804.  Ruckstuhl, K. E., Johnson, E. A., & Miyanishi, K. (2008). Introduction. The boreal forest and global change. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 363(1501), 2243-2247.  Roy, D. P., Ju, J., Kline, K., Scaramuzza, P. L., Kovalskyy, V., Hansen, M. C., Loveland, T. R., Vermote, E. & Zhang, C. (2010). Web-enabled Landsat Data (WELD): Landsat ETM+ composited mosaics of the conterminous United States. Remote Sensing of Environment, 114(1), 35–49.  Sader, S. A., Waide, R. B., Lawrence, W. T., & Joyce, A. T. (1989). Tropical forest biomass and successional age class relationships to a vegetation index derived from Landsat TM data. Remote Sensing of Environment, 28, 143-198.   Schroeder, T. A., Wulder, M. A., Healey, S. P., & Moisen, G. G. (2011). Mapping wildfire and clearcut harvest disturbances in boreal forests with Landsat time series data. Remote Sensing of Environment, 115(6), 1421-1433.  Schroeder, T. A., Cohen, W. B., & Yang, Z. (2007). Patterns of forest regrowth following clearcutting in western Oregon as determined from a Landsat time-series. Forest Ecology and Management, 243(2), 259-273.  Schroeder, T. A., Cohen, W. B., Song, C., Canty, M. J., & Yang, Z. (2006). Radiometric correction of multi-temporal Landsat data for characterization of early successional forest patterns in western Oregon. Remote Sensing of Environment, 103(1), 16-26.  Shatford, J. P. A., Hibbs, D. E., & Puettmann, K. J. (2007). Conifer regeneration after forest fire in the Klamath-Siskiyous: How much, how soon?. Journal of Forestry, 105(3), 139-146. 163   Shvidenko, A., & Apps, M. (2006). The International Boreal Forest Research Association: understanding boreal forests and forestry in a changing world. Mitigation and Adaptation Strategies for Global Change, 11(1), 5-32.  Skakun, R. S., Wulder, M. A., & Franklin, S. E. (2003). Sensitivity of the Thematic Mapper enhanced wetness difference index to detect mountain pine beetle red-attack damage. Remote Sensing of Environment, 86(4), 433−443.  Slayback, D. A., Pinzon, J. E., Los, S. O., & Tucker, C. J. (2003). Northern hemisphere photosynthetic trends 1982–99. Global Change Biology, 9(1), 1-15.  Song, C., Schroeder, T. A., & Cohen, W. B. (2007). Predicting temperate conifer forest successional stage distributions with multitemporal Landsat Thematic Mapper imagery. Remote Sensing of Environment, 106(2), 228-237.  Song, C., & Woodcock, C. E. (2003). Monitoring forest succession with multitemporal Landsat images: factors of uncertainty. IEEE Transactions on Geoscience and Remote Sensing, 41(11), 2557-2567.  Steyaert, L. T., Hall, F. G., & Loveland, T. R. (1997). Land cover mapping, fire regeneration, and scaling studies in the Canadian boreal forest with 1 km AVHRR and Landsat TM data. Journal of Geophysical Research: Atmospheres (1984–2012), 102(D24), 29581-29598.  Stocks, B. J.,  Mason, J. A.,  Todd, J. B.,  Bosch, E. M.,  Wotton, B. M.,  Amiro, B. D.,  Flannigan, M. D.,  Hirsch, K. G.,  Logan, K. A.,  Martell, D. L.,  Skinner, W.R. (2002). Large forest fires in Canada, 1959–1997. Journal of Geophysical Research: Atmospheres (1984–2012), 107(D1), FFR-5.  Stocks, B. J., Fosberg, M. A., Lynham, T. J., Mearns, L., Wotton, B. M., Yang, Q., Jin, J. Z.,  Lawrence, K., Hartley, G. R.,  Mason, J. A.,   & McKenney, D. W. (1998). Climate change and forest fire potential in Russian and Canadian boreal forests. Climatic change, 38(1), 1-13.  Sulla-Menashe, D., Friedl, M. A., & Woodcock, C. E. (2016). Sources of bias and variability in long-term Landsat time series over Canadian boreal forests. Remote Sensing of Environment, 177, 206-219.  Taylor, A. R.,  & Chen, H.Y.H. (2011). Multiple successional pathways of boreal forest stands in central Canada. Ecography, 34(2), 208-219.  Turubanova, S., Potapov, P., Krylov, A., Tyukavina, A., McCarty, J. L., Radeloff, V. C., & Hansen, M. C. (2015). Using the Landsat data archive to assess long-term regional 164  forest dynamics assessment in Eastern Europe, 1985-2012. The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 40(7), 531.  Ung, C., Benrier, P., Guo X. (2008). Canadian national biomass equations: new parameter estimates that include British Columbia data. Canadian Journal of Forest Research. 38, 1123–1132.  Urquizo, N., Bastedo, J., Brydges, T., & Shear, H. (2000). Ecological assessment of the boreal shield ecozone. Environment Canada, Ottawa, Ont. ISBN-662-28679-0.  Viglas, J. N.,  Brown, C.D, Johnstone, J.F. (2013). Age and size effects on seed productivity of northern black spruce. Canadian Journal of Forest Research. 43(6), 534-543  Vogelmann, J. E., Gallant, A. L., Shi, H., & Zhu, Z. (2016). Perspectives on monitoring gradual change across the continuity of Landsat sensors using time-series data. Remote Sensing of Environment.  Volney, W. J.,  Hirsh, K.G. (2005). Disturbing forest disturbances. The Forestry Chronicle. 81(5), 662-668.  Weber, M. G., & Flannigan, M. D. (1997). Canadian boreal forest ecosystem structure and function in a changing climate: impact on fire regimes. Environmental Reviews, 5(3-4), 145-166.  Wesner, J. S., Billman, E. J., & Belk, M. C. (2012). Multiple predators indirectly alter community assembly across ecological boundaries. Ecology, 93(7), 1674-1682.  White, J. C., Wulder, M. A., Hobart, G. W., Luther, J. E., Hermosilla, T., Griffiths, P., Coops, N. C.,  Hall, R. J.,  Hostert, P., Dyk, A., & Guindon, L. (2014). Pixel-based image compositing for large-area dense time series applications and science. Canadian Journal of Remote Sensing, 40(3), 192-212.  White, J. C., & Wulder, M. A. (2014). The Landsat observation record of Canada: 1972–2012. Canadian Journal of Remote Sensing, 39(6), 455-467. Wiken, E.B. (1986). Terrestrial ecozones of Canada.  Ecological Land Classification Series No. 19. Environment Canada, Hull, Quebec.  Wilmking, M., Juday, G. P., Barber, V. A., & Zald, H. S. (2004). Recent climate warming forces contrasting growth responses of white spruce at treeline in Alaska through temperature thresholds. Global Change Biology, 10(10), 1724-1736.  Wofsy, S. C., & Harris, R. C. (2002). The North American Carbon Program (NACP): Report of the NACP Committee of the U.S. Interagency Carbon Cycle Science Program. Washington, D.C.: U.S. Global Change Research Program.  165  Woo, M. K., Thorne, R., Szeto, K., & Yang, D. (2008). Streamflow hydrology in the boreal region under the influences of climate and human interference. Philosophical Transactions of the Royal Society of London B: Biological Sciences, 363(1501), 2249-2258.  Woodcock, C. E.,  R. Allen, M. Anderson, A. Belward, R. Bindschadler, W. Cohen, F. Gao, S.N. Goward, D. Helder, E. Helmer, R. Nemani, L. Oreopoulos, J. Schott, P.S. Thenkabail, E.F. Vermote, J. Vogelmann, M.A. Wulder, and R. Wynne. (2008). Free access to Landsat imagery. Science. 302(5879), 1011.   Wulder, M. A.,  White, J. C., Bater, C. W., Coops, N. C., Hopkinson, C., & Chen, G. (2012a). LiDAR plots - A new large-area data collection option: context, concepts, and case study. Canadian Remote Sensing Journal, 38(5):600–618.  Wulder, M. A.,  White, J. C., Nelson R. F.,  Næsset E., Orka, H. O.,  Coops, N. C.,  Hilker T., Bater, C.W. & Gobakken, T. (2012b). Lidar sampling for large area forest characterization: A Review. Remote Sensing of Environment, 121, 196-209.  Wulder, M. A., Masek, J. G., Cohen, W. B., Loveland, T. R., & Woodcock, C. E. (2012c). Opening the archive: How free data has enabled the science and monitoring promise of Landsat. Remote Sensing of Environment, 122, 2-10.  Wulder, M. A., White, J. C., Masek, J. G., Dwyer, J., & Roy, D. P. (2011). Continuity of Landsat observations: Short term considerations. Remote Sensing of Environment, 115(2), 747–751 Wulder, M. A.; White, J. C.; Gillis, M. D.; Walsworth, N.; Hansen, M. C., and Potapov, P. (2010). Multiscale satellite and spatial information and analysis framework in support of a large-area forest monitoring and inventory update. Environmental Monitoring and Assessment, 170:417-433.  Wulder, M. A., White, J. C., Alvarez, F., Han, T., Rogan, J., & Hawkes, B. (2009). Characterizing boreal forest wildfire with multi-temporal Landsat and LIDAR data. Remote Sensing of Environment, 113(7), 1540-1555.  Wulder, M. A.,  White, J. C.,  Cranny, M. M.,  Hall, R. J.,  Luther, J. E.,  Beaudoin, A. R.,  Goodenough, D. G.,  & Dechka, J.A. (2008a). Monitoring Canada's forests. Part 1: Completion of the EOSD Land Cover Project. Canadian Journal of Remote Sensing, 34(6), 549−562.  Wulder, M. A., Ortlepp, S. M., White, J. C., & Maxwell, S. (2008b). Evaluation of Landsat-7 SLC-off image products for forest change detection. Canadian Journal of Remote Sensing, 34(2), 93-99.  166  Wulder, M. A.,  Campbell, C., White, J. C.,  Flannigan, M., & Campbell, I.D. (2007). National circumstances in the international circumboreal community. The Forestry Chronicle. 83(4), 539-556.  Wulder, M. A., Skakun, R. S., Kurz, W. A., & White, J. C. (2004). Estimating time since forest harvest using segmented Landsat ETM+ imagery. Remote Sensing of Environment, 93(1-2), 179−187.  Wulder, M. A.,  & Seemann, D. (2003). Forest inventory height update through the integration of LiDAR data with segmented Landsat imagery. Canadian Journal of Remote Sensing, 29(5), 536–543.  Wulder, M. A., & Seemann, D. (2001). Spatially partitioning Canada with the Landsat worldwide referencing system. Canadian Journal of Remote Sensing, 27(3), 225-231.  Wulder, M. (1998). Optical remote-sensing techniques for the assessment of forest inventory and biophysical parameters. Progress in physical Geography, 22(4), 449-476.  Xu, L., Myneni, R. B.,  Chapin III, F. S.,  Callaghan, T. V.,  Pinzon, J. E.,  Tucker, C. J.,  Zhu, Z., Bi, J., Ciais, P., Tommervik, H., Euskirchen, E. S.,  Forbes, B. C.,  Piao, S. L.,  Anderson, B. T., Ganguly, S., Nemani, R. R.,  Goet, S. J.,  Beck, P.S. A.,  Bunn, A. G.,  Cao, C., Stroeve, J.C. (2013). Temperature and vegetation seasonality diminishment over northern lands. Nature Climate Change Letters. 1836, 1-6   Yukon Ecoregions Working Group. (2004). Ecoregions of the Yukon Territory: Biophysical Properties of Yukon Landscapes. Agriculture and Agri-Food Canada, Summerland, British Columbia.  Zhang, Q., Pavlic, G., Chen, W., Latifovic, R., Fraser, R., & Cihlar, J. (2004). Deriving stand age distribution in boreal forests using SPOT Vegetation and NOAA AVHRR imagery. Remote Sensing of Environment, 91(3), 405-418.  Zhang, X., Brown, R., Vincent, L., Skinner, W., Feng, Y. and Mekis, E. (2011). Canadian climate trends, 1950-2007. Canadian Biodiversity: Ecosystem Status and Trends 2010, Technical Thematic Report No. 5. Canadian Councils of Resource Ministers. Ottawa, ON.  Zhou, L., Tucker, C. J., Kaufmann, R. K., Slayback, D., Shabanov, N. V., & Myneni, R. B. (2001). Variations in northern vegetation activity inferred from satellite data of vegetation index during 1981 to 1999. Journal of Geophysical Research: Atmospheres, 106(D17), 20069-20083.  Zhu, Z., Woodcock, C. E., & Olofsson, P. (2012a). Continuous monitoring of forest disturbance using all available Landsat imagery. Remote Sensing of Environment, 122, 75-91.  167  Zhu, Z., & Woodcock, C. E. (2012b). Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sensing of Environment, 118, 83-94.    

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-0307371/manifest

Comment

Related Items