Automated approaches for snow and ice cover monitoring using optical remote sensing

Snow and ice cover exhibits a high degree of spatial and temporal variability. Data from multispectral optical remote sensing instruments such as Landsat are an underutilized resource that can extend our ability for mapping these phenomena. High resolution imagery is used to demonstrate that even at finer spatial resolutions (below 100 m), pixels with partial snow cover are common throughout the year and nearly ubiquitous during the meltout period. This underscores the importance of higher spatial resolution datasets for snow cover monitoring as well as the utility of fractional snow covered area (fSCA) monitoring approaches. Landsat data are used to develop a fully automated approach for mapping persistent ice and snow cover (PISC). This approach relies on the availability of numerous Landsat scenes, an improved technique for automated cloud cover mapping, and a series of automated postprocessing routines. Validation at 12 test sites suggest that the automated PISC mapping approach provides a good approximation of debris-free glacier extent across the Arctic. The PISC mapping approach is then used to produce the first single-source, temporally well-constrained (2010-2014) map of PISC across the conterminous western U.S. The Landsat-derived PISC map is more accurate than both a previously published dataset based on aerial photography acquired during the 1960s, 1970s and 1980s and the National Land Cover Database (NLCD) 2011 extent of perennial snow and ice cover. Further analysis indicates differences between the newly developed Landsat-derived PISC dataset and the previously published glacier dataset can likely be attributed to changes in the extent of PISC over time. Finally, in order to map mean annual snow cover persistence across the entire landscape, we implement a novel canopy adjustment approach designed to improve the accuracy of Landsat-derived fSCA in forested areas. In situ observations indicate canopy-adjusted snow covered area calculated from all available Landsat scenes can provide an accurate estimate of mean annual snow cover duration. The work presented here lays the groundwork for addressing scientific questions regarding the spatial and temporal variability of snow cover, snow accumulation and ablation processes, and the impact of changes in snow cover on physical and ecological systems.

AUTOMATED APPROACHES FOR SNOW AND ICE COVER MONITORING USING OPTICAL REMOTE SENSING by David James Selkowitz A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Geography College of Social and Behavioral Science The University of Utah August 2017 Copyright © David James Selkowitz 2017 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of David James Selkowitz has been approved by the following supervisory committee members: Richard Forster , Chair 2/3/2017 Phillip Dennison , Member Date Approved 2/3/2017 Date Approved Simon Brewer , Member 2/3/2017 Date Approved Mitchell Power , Member 2/3/2017 Date Approved Dennis Dye , Member 2/3/2017 Date Approved and by Andrea Brunelle the Department of Geography and by David B. Kieda, Dean of The Graduate School. , Chair of ABSTRACT Snow and ice cover exhibits a high degree of spatial and temporal variability. Data from multispectral optical remote sensing instruments such as Landsat are an underutilized resource that can extend our ability for mapping these phenomena. High resolution imagery is used to demonstrate that even at finer spatial resolutions (below 100 m), pixels with partial snow cover are common throughout the year and nearly ubiquitous during the meltout period. This underscores the importance of higher spatial resolution datasets for snow cover monitoring as well as the utility of fractional snow covered area (fSCA) monitoring approaches. Landsat data are used to develop a fully automated approach for mapping persistent ice and snow cover (PISC). This approach relies on the availability of numerous Landsat scenes, an improved technique for automated cloud cover mapping, and a series of automated postprocessing routines. Validation at 12 test sites suggest that the automated PISC mapping approach provides a good approximation of debris-free glacier extent across the Arctic. The PISC mapping approach is then used to produce the first single-source, temporally well-constrained (2010-2014) map of PISC across the conterminous western U.S. The Landsat-derived PISC map is more accurate than both a previously published dataset based on aerial photography acquired during the 1960s, 1970s and 1980s and the National Land Cover Database (NLCD) 2011 extent of perennial snow and ice cover. Further analysis indicates differences between the newly developed Landsat-derived PISC dataset and the previously published glacier dataset can likely be attributed to changes in the extent of PISC over time. Finally, in order to map mean annual snow cover persistence across the entire landscape, we implement a novel canopy adjustment approach designed to improve the accuracy of Landsat-derived fSCA in forested areas. In situ observations indicate canopy-adjusted snow covered area calculated from all available Landsat scenes can provide an accurate estimate of mean annual snow cover duration. The work presented here lays the groundwork for addressing scientific questions regarding the spatial and temporal variability of snow cover, snow accumulation and ablation processes, and the impact of changes in snow cover on physical and ecological systems. iv This dissertation is dedicated to Suka, superstar Siberian husky and connoisseur of fine snow conditions from Alaska to Utah and everywhere in between. You will be missed. TABLE OF CONTENTS ABSTRACT ................................................................................................................................................... iii LIST OF TABLES ...................................................................................................................................... viii LIST OF FIGURES ........................................................................................................................................ x ACKNOWLEDGEMENTS ...................................................................................................................... xiv Chapters 1. INTRODUCTION ................................................................................................................................... 1 2. PREVALENCE OF PURE VERSUS MIXED SNOW COVER PIXELS ACROSS SPATIAL RESOLUTIONS IN ALPINE ENVIRONMENTS .......................................................................... 6 2.1 Introduction ............................................................................................................................ 8 2.2 Study Area and Methods ................................................................................................... 9 2.3 Results .................................................................................................................................... 21 2.4 Discussion ............................................................................................................................. 30 2.5 Conclusions .......................................................................................................................... 33 3. AN AUTOMATED APPROACH FOR MAPPING PERSISTENT ICE AND SNOW COVER OVER HIGH LATITUDE REGIONS .............................................................................. 38 3.1 Introduction ......................................................................................................................... 39 3.2 Study Regions ..................................................................................................................... 41 3.3 Data ......................................................................................................................................... 42 3.4 Methods ................................................................................................................................. 44 3.5 Results .................................................................................................................................... 48 3.6 Discussion ............................................................................................................................. 54 3.7 Conclusions .......................................................................................................................... 56 4. AUTOMATED MAPPING OF PERSISTENT ICE AND SNOW COVER ACROSS THE WESTERN US WITH LANDSAT .......................................................................................................... 60 4.1 Introduction ......................................................................................................................... 61 4.2 Study Area and Methods ................................................................................................ 62 4.3 Results .................................................................................................................................... 68 4.4 Discussion ............................................................................................................................. 71 4.5 Conclusions .......................................................................................................................... 74 5. THE USGS LANDSAT SNOW COVERED AREA PRODUCTS: METHODS AND PRELIMINARY VALIDATION ....................................................................................................... 76 5.1 Introduction ......................................................................................................................... 77 5.2 Study Areas and Methods .............................................................................................. 81 5.3 Results .................................................................................................................................... 97 5.4 Discussion .......................................................................................................................... 107 5.5 Conclusions ....................................................................................................................... 112 5.6 Acknowledgements………………………………………………………………………… 114 5.7 References……………………………………………………………………………………… 114 6. CONCLUSIONS ................................................................................................................................. 120 6.1 References………………………………………………………………………………………..125 vii LIST OF TABLES 2.1. Geographic and climatic characteristics for each study area. Climatic characteristics are derived from the PRISM Climate Dataset………………………...... 10 2.2 Study area, image type, date, snow cover fraction, 0.5 m binary SCA classification accuracy, and recent snowfall history (in terms of snow water equivalent) for each WorldView (panchromatic) or WorldView-2 (3 band) image strip included in the analysis……………………………………………………………………….... 13 2.3 Pixel sizes used to compute frequency of mixed pixels and corresponding pixel sample sizes for each ISA………………………………………………….………………………...... 15 2.4 Study area location and site characteristics for 60 x 60 m grid cell footprints where in situ monitoring was conducted..……..…………………….………………………........... 19 2.5 Full, partial, and total snow cover days, as well as snow-covered to snow-free transition period metrics for 60 x 60 m footprints.………………….………………………...... 27 2.6 Snow cover fraction uncertainty metrics for each 60 m grid cell..……..……...... 29 3.1 Study area locations and characteristics.……..………………………………………....…....42 3.2 Comparison of cloud-free and shadow-free surface views (CFSFSV) for each study area using the original CFmask and the revised cloud masking approach…………………………..………………………………………………….………………………...... 51 4.1 F metric for USGS Topographic Maps glacier layer, NLCD 2011 perennial snow and ice cover, initial Landsat-derived PISC and revised Landsat-derived PISC........... 68 4.2 Accuracy, precision, recall and F metrics for 1000 randomly selected validation points……………………………..………………………………………………….………………………........... 72 4.3 Area of PISC mapped by the USGS DRG dataset, the automated Landsat dataset, and the NLCD 2011 dataset by region..……..………………………….………………………...... 73 4.4 PISC area mapped by semi-automated Landsat approach for 1987-1988 and 2008-2010 for three test regions..……..………………………………….………………………...... 73 5.1 Landsat scenes used by path row.……..…………………………………………………...... 84 5.2 30 x 30 km subsets used for detailed analysis of the Landsat mean annual snow cover duration product….……..………………………………………….………………………............... 84 5.3 Criteria for identification of surrogate pixels in the surrounding 11 x 11 or 31 x 31 pixel neighborhoods……..………………………………………………….………………………......... 90 5.4 Accuracy metrics for TMSCAG, canopy adjusted TMSCAG, and neighborhood canopy adjusted TMSCAG….………………………………………………….………………………...... 98 5.5 Accuracy metrics for mean annual snow cover duration (days) calculated using unadjusted TMSCAG and canopy adjusted TMSCAG relative to mean annual snow cover duration calculated from SNOTEL sites……………………………………………...... 106 5.6 Accuracy metrics for mean annual snow cover days calculated using adjusted TMSCAG for periods 1991-1995, 1996-2000, 2001-2005, 2006-2011, and 20112015……………………..……..………………………………………………….………………………...... 107 ix LIST OF FIGURES 2.1. Study area locations in the Western United States………………………………............ 11 2.2 (a) Oregon Cascades ISA and (b) Rocky Mountain NP ISA. …….…….……....… 12 2.3 Semivariograms for elevation derived from 10 m DEM for (a) the Oregon Cascades and Rocky Mountain NP ISAs and, (b) high- and low-relief subsets from the Oregon Cascades and Rocky Mountain NP ISAs…………………………………….……....….. 16 2.4 Locations for 60 x 60 m grid cells instrumented with arrays of temperature data loggers at (a) the Cinnamon Pass FSA and (b) the Niwot Ridge FSA……..….......... 18 2.5 Schematic diagram of the arrangement of temperature data loggers at each 60 x 60 m grid cell footprint. ……………………..…………………………….………………………........... 19 2.6 WorldView imagery (top row) and high resolution binary snow-covered area maps (bottom row) for the Oregon Cascades ISA…………………….……..………………...... 22 2.7 WorldView imagery (top row) and high-resolution binary snow-covered area maps (bottom row) for the Rocky Mountain NP ISA………………………………………........... 23 2.8 Examples of 0.5 m binary snow-covered area classifications (right side) based on WorldView 2 imagery (left side)………………………………………….…………………...... 24 2.9 Fraction of mixed (partially snow-covered) pixels for pixel resolutions between 1 m and 500 m for (a) the Oregon Cascades imagery study area, and (b) the Rocky Mountain NP imagery study area. …………………………….………………………............ 24 2.10 (a) Fraction of mixed pixels for Oregon Cascades ISA high-relief subset, (b) fraction of mixed pixels for the Rocky Mountain NP ISA high-relief subset, (c) fraction of mixed pixels for the Oregon Cascades ISA low-relief subset, and (d) fraction of mixed pixels for the Rocky Mountain NP ISA low-relief subset………...... 25 2.11 Relationship between pixel resolution and total study area snow cover fraction derived from binary SCA for (a) Oregon Cascades ISA and (b) Rocky Mountain NP imagery study area.……………………………………………………….….…………………………...... 26 2.12 Spatial distribution of 400 x 400 m blocks from the Rocky Mountain NP imagery study area where absolute differences between binary and fractional SCA were low (< 0.05), medium (0.05-0.10), high (0.10-0.15), and highest (>0.15)…...... 26 2.13 Daily 60 m grid cell snow cover fraction time series from (a) grid cells 1, 5 and 6 at Cinnamon Pass FSA (b) grid cells 3 and 4 at Cinnamon Pass FSA, and (c) grid cells 8, 9, and 10 at Niwot Ridge FSA..……..……………………………………………….….....….. 28 2.14 Results from sensitivity analysis indicating 60 m grid cell snow cover fraction using 0.5 °C, 1 °C, 2 °C, and 3 °C temperature thresholds for (a) an example site with partially snow-covered conditions throughout the winter and spring (site #2) and (b) for an example site with consistent fully snow-covered conditions throughout the winter and spring (site #9)..……..………………………………………………….……......…….. 29 3.1 Study area locations across the circumpolar Arctic..…….…….…....….. 42 3.2 Diagram of processing flow for classification of Persistent Ice and Snow Cover (PISC) for a single pixel..……..…………………………………………………………..…………....….. 44 3.3 Comparison of original CFmask and revised CFmask for Bylot Island….........….. 45 3.4 Agreement with the Randolph Glacier Inventory (RGI) glacier extent for each full study area and validation study area. ..……..…………………………………………....…….. 48 3.5 Agreement with VHRI-derived glacier extent for each validation study area…………………………………………………...……..……………………………………………….…........ 49 3.6 Agreement between PISC and glaciers mapped with RGI for the Trollaskagi Peninsula, Iceland, with area of detail showing Landsat imagery and areas of false positives for PISC outlined in red...……..……………………..…………………………………...... 50 3.7 Cloud-free and shadow-free views for a subset of the Brooks Range study area using the original CFmask and the revised cloud masking approach..……………........... . 51 3.8 Accuracy (agreement with the RGI) for each study area as a function of late summer snow cover days threshold using the original CFmask and the revised cloud masking approach……………………………...……..…………………………………………………...... 51 3.9 Effect of fraction of Days With Ice and Snow Cover (fDISC) threshold and Normalized Difference Snow Index (NDSI) threshold on PISC map accuracy (defined as agreement with RGI glacier area)..……..……………………………………………………...... 52 3.10 Effect of number of cloud and shadow free views on precision (user's accuracy for the ice covered class) and recall (producer's accuracy for the ice-covered class) for RGI validation study areas and VHRI validation study areas………………………...... 53 xi 4.1 Overview map of the western United States, including glacier regions defined by the USGS DRG dataset, very high resolution imagery (VHRI) validation subsets, change analysis regions, and locations for other figures…………………………………...... 63 4.2 Diagram of processing flow for classification of Persistent Ice and Snow Cover (PISC)…………………………………………………………………………..…………………………….......….. 65 4.3 An example of differences between cloud cover mapped by the CFmask algorithm and the revised cloud masking approach for the Wind River Range of Wyoming…………………………………………………………………….………………………………......... 66 4.4 An example of differences between the number of available cloud-free shadowfree views resulting from use of the original CFmask algorithm and the revised cloud masking approach for the North Cascades of Washington………………………………...... 66 4.5 Agreement between PISC mapped using late summer WorldView 2 imagery and PISC datasets for the Mammoth Glacier (Wyoming) VHRI subset.…….……….......... 69 4.6 Agreement between PISC mapped using late summer WorldView 2 imagery and PISC mapped by the automated Landsat approach plus ancillary data for all 13 VHRI subsets……………………………………………………………….………………………………...... 70 4.7 Relationship between fraction of validation subset with PISC and Accuracy and F metrics……………………………..……………………………………….………………………………...... 71 4.8 Accuracy, precision, recall and F metrics for 1000 points randomly selected from areas on or near (within 500 m) of previously mapped glaciers and perennial snow cover patches.…………….……………………………………….………………………………...... 71 4.9 Relationship between the number of cloud-free, shadow-free views and F metric for all pixels from VHRI validation subsets….…….………………………………...... 72 4.10 Area of PISC mapped by the USGS DRG dataset and by the automated Landsatderived dataset for 8 regions of the western U.S…………...………………………………...... 72 5.1 Study area locations in the western U.S., including Landsat path/rows used, 30 x 30 km analysis subsets, locations of SNOTEL sites used, and locations of in situ fSCA sites………….……………………………………………………………………………….………...... 83 5.2 Flow chart for canopy adjustment processes used for adjustment of scenebased Landsat fSCA. ………………………………………………………………..……….…………...... 85 5.3 Identification of surrogate pixels meeting criteria listed in Table 5.2 in the 11 x 11 pixel neighborhood surrounding target pixel.………….…………………………….……...... 91 xii 5.4 Calculation of the fraction of days in June with snow cover for the period 19862015 for a single 30 m grid cell. ………….……………………………………………….………...... 93 5.5 TMSCAG canopy-adjusted fSCA compared to in situ fSCA calculated from temperature data logger arrays. ………………………………………………………….………...... 99 5.6 Demonstration of canopy adjustment for an area in the northern Sierra Nevada on April 20, 2009.………….…………………………………………………………………….………...... 100 5.7 Mean annual snow cover duration for three 30 x 30 km subsets. ……………..... 101 5.8 Additional snow cover days added using the neighborhood canopy adjustment approach for three 30 x 30 km subsets……..………………………………………….………...... 102 5.9 Percent of all cloud-free pixels where snow cover was added via canopy adjustment, shown by month.………….……………………………………………………………...... 103 5.10 Canopy adjustment model failure frequency for three 30 x 30 km subsets………………………………………………………………………………………………………....... 104 5.11 Comparison between mean annual snow cover duration calculated from SNOTEL data and mean annual snow cover duration calculated from Landsat….... 106 xiii ACKNOWLEDGEMENTS I wish to acknowledge all members of my committee for all they have taught me. I also wish to acknowledge the USGS Land Remote Sensing Program for funding much of the research contained within this dissertation. Finally, I wish to thank Melanie Cota, who encouraged and supported me throughout the entire long process of earning a PhD. CHAPTER 1 INTRODUCTION Glaciers and seasonal snow cover serve as a crucial water resource across many regions of the world. Changes in glaciers and seasonal snow cover also serve as a key climate indicator across the globe and are particularly valuable where long term in situ measurements are unavailable. The quantity and timing of seasonal snow exerts a strong (and in some cases dominant) influence on a wide range of arctic, alpine, montane and boreal ecosystem processes, ranging from tree seedling establishment to ungulate travel routes and habitat selection. At the broadest spatial scales (i.e., regional to global), existing remotely sensed and modeled data products provide a clear picture of the status and variability in snow and ice cover and have been used extensively in assessments of global and regional climate. The spatial resolution of the existing datasets (such as those derived from MODIS), however, is often insufficient to resolve much of the finer scale variability in snow and ice cover that impacts regional water supplies and ecological processes. For instance, a collection of late lying snow patches approximately 1 ha in size that would not typically be resolved in most regional to global scale snow cover products can provide substantial summer runoff or serve as 2 crucial habitat for caribou seeking refuge from mosquitos. Consequently, the ability to map and monitor snow and ice cover at finer spatial resolutions has the potential to provide major benefits for science and society. The Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) have provided a nearly continuous time series of optical remote sensing data appropriate for snow cover mapping across many regions of the globe since the mid 1980s. Landsat data have often been overlooked as a tool for snow cover monitoring because the 16-day interval between scene acquisitions is insufficient for many snow cover monitoring applications requiring a higher temporal resolution. However, the high spatial resolution (nominally 30 m) and extensive archive of scenes stretching back to the mid 1980s allow for a unique contribution to snow cover monitoring and snow hydrology. For many regions, the Landsat archive can be used to compute snow cover duration metrics (e.g., monthly and for the entire year for periods of 5 years or more) as well as the extent of glaciers and perennial snow cover. In addition to this introduction and a concluding chapter, this dissertation is composed of four individual chapters which have been published or submitted as peer-reviewed journal manuscripts which address the potential for using data from Landsat or similar multispectral instruments for mapping and monitoring snow and ice cover at 30 m spatial resolution. Chapter 2, "Prevalence of pure versus mixed snow cover pixels across spatial resolutions in alpine environments," uses a combination of very high spatial resolution spaceborne imagery and in situ measurements to document the 3 frequency of partially snow-covered pixels in mountainous environments. The data presented in this chapter indicate that in mountain regions, pure snow cover pixels are rare at the scale of most regional to global datasets such as the MODIS snow cover products, and quite uncommon even at the Landsat spatial scale. The data show that even in locations where deep and spatially contiguous snowpacks accumulate during the winter, heterogeneity in accumulation and ablation processes results in extended periods during the spring and summer when Landsat scale pixels are partially snow covered. This underscores the utility of higher spatial resolution datasets, as well as the benefits of remote sensing approaches that provide estimates of fractional snow covered area for each pixel over more traditional approaches that only provide binary snow covered area. Chapter 3, "An automated approach for mapping persistent ice and snow cover over high latitude regions," documents an approach developed for automated classification of glaciers and perennial snow cover across high latitude regions. A key development presented in this manuscript is a revised version of the CFMask algorithm for cloud masking optimized for use in mountainous regions where snow, ice, and rock surfaces are commonly located in close proximity to one another. Mixed areas of snow, ice and rock are frequently misclassified by the original version of the CFMask algorithm included with Landsat surface reflectance products. The revised cloud masking algorithm incorporates the original CFMask cloud cover classification but substantially reduces errors of commission for cloud cover. 4 Chapter 4, "Automated mapping of persistent ice and snow cover across the western US with Landsat," expands upon the work presented in Chapter 2. In this manuscript, the approach developed in the previous chapter is adapted for use at lower latitudes and applied to the entire western conterminous United States. Validation of the persistent ice and snow cover (PISC) map using high spatial resolution imagery indicates the new dataset is more accurate than the 2011 National Land Cover Database snow/ice cover map and more accurate than a USGS atlas of glacier outlines compiled from topographic maps based on aerial photography. While the higher accuracy relative to the NLCD snow/ice cover class can be attributed to differences in mapping methods, the differences between the Landsat-derived PISC dataset and the dataset compiled from topographic maps appear to be due primarily to decreases in the extent of persistent ice and snow cover over time. Chapter 5, "The USGS Landsat snow covered area products: methods and preliminary validation," describes the development and validation of a new set of Landsat-derived snow covered area products. These products are now available for production on demand by user request and will eventually be included as standard products available alongside raw Landsat imagery and surface reflectance data. The Landsat snow covered area products include scene-based fractional snow covered area and canopy-adjusted fractional snow covered area as well as mean annual snow cover duration computed over the period 1986-2015. While viewable fractional snow covered area is validated in a separate publication, canopy adjusted fractional snow covered area for individual Landsat scenes is validated using a 5 network of in situ sensor arrays in the Sierra Nevada. Mean annual snow cover duration for the period 1986-2015, as well as for shorter 5-year periods, is validated using data from the SNOTEL network in California, Washington, and Wyoming. Results indicate that the RMSE for scene-based canopy-adjusted fSCA is 0.21, while the RMSE for 30-year mean annual snow cover duration is 14.7 days. The canopy adjustment approach introduced in this manuscript substantially improves accuracy and reduces bias for both scene-based fSCA and mean annual snow cover duration. While the research presented in Chapters 2-5 focuses on the remote sensing and image processing approaches necessary for development and production of Landsat-derived snow and ice cover datasets, Chapter 6 briefly delves into the broader scientific questions that can be addressed using these datasets. CHAPTER 2 PREVALENCE OF PURE VERSUS MIXED SNOW COVER PIXELS ACROSS SPATIAL RESOLUTIONS IN ALPINE ENVIRONMENTS Selkowitz, D. J., Forster, R. R., & Caldwell, M. K. (2014). Prevalence of pure versus mixed snow cover pixels across spatial resolutions in alpine environments. Remote Sensing, 6(12), 12478-12508. Published by MDPI 2014. This work was produced as part of the lead author's fulfillment of official government duties and is therefore considered to be in the public domain and not subject to copyright protection. 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 CHAPTER 3 AN AUTOMATED APPROACH FOR MAPPING PERSISTENT ICE AND SNOW COVER OVER HIGH LATITUDE REGIONS Selkowitz, D. J., & Forster, R. R. (2015). An automated approach for mapping persistent ice and snow cover over high latitude regions. Remote Sensing, 8(1), 16. Published by MDPI 2015. This work was produced as part of the lead author's fulfillment of official government duties and is therefore considered to be in the public domain and not subject to copyright protection. 39 Article An Automated Approach for Mapping Persistent Ice and Snow Cover over High Latitude Regions David J. Selkowitz 1,2, * and Richard R. Forster 2 Received: 30 June 2015; Accepted: 21 December 2015; Published: 25 December 2015 Academic Editors: Daniel J. Hayes, Santonu Goswami, Guido Grosse, Benjamin Jones, Richard Gloaguen and Prasad S. Thenkabail 1 2 * U.S. Geological Survey, Alaska Science Center, 4210 University Drive, Anchorage, AK 99508, USA Department of Geography, University of Utah, 260 S. Central Campus Dr., Room 270, Salt Lake City, UT 84112, USA; [e-mail] Correspondence: [e-mail]; Tel.: [phone]; Fax: [phone] Abstract: We developed an automated approach for mapping persistent ice and snow cover (glaciers and perennial snowfields) from Landsat TM and ETM+ data across a variety of topography, glacier types, and climatic conditions at high latitudes (above ~65˝ N). Our approach exploits all available Landsat scenes acquired during the late summer (1 August-15 September) over a multi-year period and employs an automated cloud masking algorithm optimized for snow and ice covered mountainous environments. Pixels from individual Landsat scenes were classified as snow/ice covered or snow/ice free based on the Normalized Difference Snow Index (NDSI), and pixels consistently identified as snow/ice covered over a five-year period were classified as persistent ice and snow cover. The same NDSI and ratio of snow/ice-covered days to total days thresholds applied consistently across eight study regions resulted in persistent ice and snow cover maps that agreed closely in most areas with glacier area mapped for the Randolph Glacier Inventory (RGI), with a mean accuracy (agreement with the RGI) of 0.96, a mean precision (user's accuracy of the snow/ice cover class) of 0.92, a mean recall (producer's accuracy of the snow/ice cover class) of 0.86, and a mean F-score (a measure that considers both precision and recall) of 0.88. We also compared results from our approach to glacier area mapped from high spatial resolution imagery at four study regions and found similar results. Accuracy was lowest in regions with substantial areas of debris-covered glacier ice, suggesting that manual editing would still be required in these regions to achieve reasonable results. The similarity of our results to those from the RGI as well as glacier area mapped from high spatial resolution imagery suggests it should be possible to apply this approach across large regions to produce updated 30-m resolution maps of persistent ice and snow cover. In the short term, automated PISC maps can be used to rapidly identify areas where substantial changes in glacier area have occurred since the most recent conventional glacier inventories, highlighting areas where updated inventories are most urgently needed. From a longer term perspective, the automated production of PISC maps represents an important step toward fully automated glacier extent monitoring using Landsat or similar sensors. Keywords: remote sensing of glaciers; snow and ice; Landsat; arctic 1. Introduction Glaciers have been identified as one of the most sensitive indicators of changes in climate [1,2] and have been identified as an essential climate variable that should be monitored globally [3]. Glaciers not only respond to changes in climate, but can also drive changes in the earth climate system through changes in albedo and contribution to sea level rise [4-7]. From a more local to regional Remote Sens. 2016, 8, 16; doi:10.3390/rs8010016 www.mdpi.com/journal/remotesensing 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 CHAPTER 4 AUTOMATED MAPPING OF PERSISTENT ICE AND SNOW COVER ACROSS THE WESTERN US WITH LANDSAT Selkowitz, D. J., & Forster, R. R. (2016). Automated mapping of persistent ice and snow cover across the western US with Landsat. ISPRS Journal of Photogrammetry and Remote Sensing, 117, 126-140. Published by Elsevier 2016. This work was produced as part of the lead author's fulfillment of official government duties and is therefore considered to be in the public domain and not subject to copyright protection. 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 CHAPTER 5 THE USGS LANDSAT SNOW COVERED AREA PRODUCTS: METHODS AND PRELIMINARY VALIDATION David J. Selkowitz 1,2 Thomas H. Painter3 Karl Rittger4 Gail Schmidt5 Richard Forster2 1 U.S. Geological Survey, Alaska Science Center 2 University of Utah, Department of Geography 3 NASA/Jet Propulsion Laboratory 4 National Snow and Ice Data Center, University of Colorado 5 Stinger Ghaffarian Technologies (SGT), Contractor to the U.S. Geological Survey (USGS) 77 5.1 Introduction Seasonal snow cover is vitally important to Earth's climate, ecology, and hydrology. Streamflow is generated primarily by snow cover runoff in many regions, and approximately one sixth of the world's population depends on snow cover for their water supply (Barnett et al., 2005). The timing and duration of seasonal snow cover is one of the key drivers for both short term fluctuations and long term changes in Earth's albedo, and therefore impacts climate dynamics at the global scale (Cohen & Entekhabi, 1999; Groisman et al., 1994). Snow cover insulates soil from cold winter temperatures (Groffman et al., 2001; Zhang, 2005) but can also inhibit thawing when air temperatures rise above the freezing point in the spring, thus altering drainage characteristics (Quinton et al., 2009). The duration of seasonal snow cover is often the dominant factor controlling the distribution of arctic and alpine plant species (Billings & Bliss, 1959; Walker et al., 1993) and can also impact the configuration of forests and meadows at the alpine treeline and below (Bekker, 2005; Hessl & Baker, 1997; Magee & Antos, 1992). Snow cover also impacts animal movement and habitat distribution (Aubry et al., 2007; Stenseth et al., 2004; Sweeney & Sweeney, 1984). The influence of snow cover on all of these crucial hydrological, climatological, and ecological processes underscores the importance of monitoring the spatial and temporal variability of snow cover across the Earth's surface at a variety of scales. Remote sensing is one of the most effective approaches for regular, spatially comprehensive snow cover monitoring. For many applications, the fine to moderate scale (10 m to 1 km) spatial distribution of snow cover is important for 78 understanding scales of controls and spatial variability (Deems et al., 2006; Tinkham et al., 2014). Several studies (Anderton et al., 2002; Luce & Tarboton, 1998) have demonstrated that in areas where heterogeneous seasonal snow covers develop, explicit representation of the spatial variability of snow cover is essential for accurate simulation of snowmelt runoff unless these parameters can be effectively accounted for in subgrid parameterization schemes. For snow simulation models, fine to moderate resolution snow cover patterns retrieved from remote sensing also provide an additional source of validation data besides runoff. Unlike runoff, however, remotely sensed snow covered area (SCA) at fine to moderate spatial resolutions can be used to assess the representation of individual processes in the model (e.g., wind redistribution of snow cover) (Bloschl et al., 1991). Lundquist and Dettinger (2005) demonstrate that the spatial heterogeneity of snow cover plays a key role in determining diurnal streamflow variations in larger basins. In colder climates, hillslope drainage is largely controlled by fine scale patterns of snow covered area because high latitude soils overlain by snow cover typically remain frozen and inhibit subsurface flow (Quinton et al., 2009). Finally, fine to moderate scale heterogeneity of snow cover controls the distribution and abundance of many plant species in arctic and alpine environments (Beck et al., 2005; Billings & Bliss, 1959; Walker et al., 1993), influences animal habitat selection (Eastland et al., 1989; LaPerriere & Lent, 1977), and impacts predator prey interactions (Huggard, 1993). The availability of standardized, freely distributed SCA products derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) has vastly 79 improved our ability to monitor and understand regional to global scale SCA patterns and variability. Remotely sensed SCA from the Landsat TM, ETM+, and OLI sensors, however, offers tremendous potential for extending this monitoring to a finer spatial scale. Landsat sensors are well suited for mapping SCA at 30 m spatial resolution. This has been demonstrated in numerous studies where Landsat data were used for snow cover mapping prior to the availability of MODIS data (Bronge† & Bronge†*, 1999; Dozier, 1984; Fily et al., 1997; Hall et al., 1989; Klein & Isacks, 1999; Rosenthal & Dozier, 1996; Winther & Hall, 1999), to assist with development and validation of MODIS snow algorithms (Hall et al., 1995; Klein et al., 1998; Painter et al., 2009; Rittger et al., 2013; Salomonson & Appel, 2004), for validation of spatially explicit snow cover models (Bernhardt et al., 2010; Bernhardt & Schulz, 2010; Fily et al., 1999; Letsinger & Olyphant, 2007), and for reconstruction of peak snow water equivalent (Cline et al., 1998; Durand et al., 2008; Margulis et al., 2016; Martinec & Rango, 1981; Molotch, 2009). Effective, snow cover remote sensing across most regions requires an approach that can reliably detect snow cover beneath forest canopies. Optical remote sensing approaches, however, typically map the viewable fraction of snow cover not obscured by forest canopy, rather than the true fraction of snow covering the ground. Although a number of algorithms have been developed with the aim of monitoring ground snow cover fraction via optical remote sensing (Klein, Hall & Riggs, 1998; Moloch & Margulis, 2008; Vikhamar & Solberg, 2003), this problem remains an active area of research. 80 The most commonly used approach for adjusting sensor-viewable fSCA to reflect the true in situ fSCA in areas with forest canopy is to assume that the viewable snow cover fraction for a given pixel will be identical (or at least similar) to the hidden (canopy-obscured) snow cover fraction for the same pixel (Coons et al., 2014; Durand & Molotch, 2008; Molotch & Margulis, 2008; Raleigh et al., 2013). Approaches assuming similar snow cover fractions for the canopy-obscured and canopy-free portions of a pixel are usually reasonably accurate as long as a snow cover fraction > 0 is retrieved at each snow covered pixel. They are, however, ineffective for identifying snow-covered pixels where the viewable snow cover fraction is less than the snow cover detection limit for the snow mapping algorithm.. There is therefore a need for a canopy adjustment approach that can identify snow cover that is frequently missed by optical remote sensing approaches. Snow covered area is the most basic measurement that can be made from optical remote sensing and serves as a key input for remote sensing or combined remote sensing/modeling approaches that endeavor to provide more complex snow metrics such as SWE. Fractional snow covered area (fSCA) provides more information per-pixel than binary SCA and is particularly useful in mountainous environments where 25-93% of all pixels at the Landsat spatial resolution are mixed pixels composed of two or more land surface types (Selkowitz et al., 2014). In order to meet the need for a standardized, analysis-ready Landsat snow cover dataset, the US Geological Survey is now producing a Landsat scene-based snow cover product based on Painter et al. (2009) that provides 30 m resolution fSCA, canopy-adjusted fSCA, and a cloud mask optimized for use in mountainous 81 environments. The first iteration of the scene-based Landsat snow cover product is available on demand for nearly any Landsat TM or ETM+ scene available in the archives stretching back to 1984. For users more interested in characterizing patterns of snow cover duration and potential changes in snow and ice cover over decades, the Landsat snow cover duration product will provide the mean annual snow cover duration days at 30 m resolution for periods as short as 5 years and as long as the full Landsat 5/Landsat 7 period of record (1984 to the present). The snow cover duration product will incorporate canopy-adjusted fSCA and cloud mask data from the Landsat scenebased snow products covering the period of interest. Initial production of the snow cover duration product will focus on 30-year mean annual snow cover duration across key mountain ranges in the western U.S., with areas outside of the region available upon request. The goals of this publication are: (1) to provide a detailed description of the methods used for production of both the scene-based fSCA and mean annual snow cover duration products, and (2) provide limited validation for each of the two products. 5.2. Study Area and Methods 5.2.1 Study Area Locations Validation of fSCA from individual Landsat scenes was conducted using in situ sensor arrays in the Sierra Nevada at sites covered by Landsat path rows 43/33, 43/34, 42/34, 42/35, and 41/35, while validation of mean annual snow cover 82 duration was conducted using SNOTEL sites in the Cascades of Washington and Oregon (path rows 45/27 and 45/28), the Sierra Nevada of California (path row 43/33), and the Rocky Mountains of Wyoming, Montana, and Idaho (path rows 38/29 and 38/30) (Figure 5.1, Table 5.1). Detailed analysis of the mean annual snow cover duration products was also conducted at one 30 x 30 km subset from each of the three regions (Figure 5.1, Table 5.2). The three regions were selected to represent the variation in snow climate regimes (Mock & Birkeland, 2000; Trujillo & Molotch, 2014) and vegetation types present across the western U.S. 5.2.2 Overview of Methods We retrieved the visible snow cover fraction for each 30 m Landsat pixel using the TMSCAG (Thematic Mapper Snow Covered Area and Grain Size) model (Painter et al., 2003; Painter et al., in review), a spectral unmixing approach. We then applied a series of adjustments (Figure 5.2) designed to produce fSCA values that more closely matched the fraction of snow covered ground (including rock, soil, low-growing vegetation, and woody debris) in areas with forest canopy. The first adjustment handled pixels where retrieved fSCA was > 0 but likely underrepresented the fraction of snow covered ground. The second adjustment handled pixels where retrieved fSCA was 0 but a combination of ancillary data and retrieved viewable fSCA from nearby pixels suggested snow cover was likely present. Finally, cloud cover was identified and masked using the revised CFMask approach described in Selkowitz and Forster (2015). These steps resulted in three layers: (1) viewable fSCA (computed directly from TMSCAG), (2) canopy-adjusted 83 Figure 5.1. Study area locations in the western U.S, including Landsat path/rows used, 30 x 30 km analysis subsets, locations of SNOTEL sites used, and locations of in situ fSCA sites. 84 Table 5.1. Landsat scenes used by path row. Landsat Type(s) of Analysis Path/Row 45/27 SNOTEL comparison, Cascades subset analysis 45/28 SNOTEL comparison 43/33 In situ sensor array comparison, SNOTEL comparison, Sierra subset analysis 43/34 In situ sensor array comparison 42/34 In situ sensor array comparison 42/35 In situ sensor array comparison 41/35 In situ sensor array 38/29 SNOTEL comparison 38/30 SNOTEL comparison, Gros Ventre subset analysis All Path Rows Scenes TM ETM+ Total 370 494 864 369 502 370 511 871 881 5 8 7 7 369 365 5 8 7 7 879 864 510 499 877 2516 4393 Table 5.2. 30 x 30 km subsets used for detailed analysis of the Landsat mean annual snow cover duration product. Forest cover indicates the percentage of pixels from the subset where the National Land Cover Database (NLCD) forest canopy layer indicates > 15% canopy cover. Mean canopy indicates the mean NLCD forest canopy for all pixels from the subset. Subset Landsat PR Elev. Range (m) Forest Cover Mean Canopy Cascades Path 45 Row 27 410-2459 84% 52.6% Sierra Nevada Path 43 Row 33 1332-3162 33% 34.6% Gros Ventre Path 38 Row 30 1751-3233 55% 26.7% 85 Figure 5.2. Flow chart for canopy adjustment processes used for adjustment of scene-based Landsat fSCA. fSCA, and (3) a cloud mask. Mean annual snow cover duration was calculated using all available Landsatderived fSCA layers available for the period of record. For each pixel, we determined the fraction of cloud-free, valid pixels with snow cover for each month (e.g., the fraction of cloud free, valid pixels imaged during the month of June over the period 1986-2015). The monthly fractions were then averaged and multiplied by 365 to calculate mean annual snow cover duration. 86 5.2.3 Production of Scene-based Landsat fSCA Product 5.2.3.1 Datasets We obtained 30 m resolution Landsat Climate Data Record (CDR) top-of- atmosphere (TOA) and surface reflectance (SR) products (Masek et al., 2006, available at http://earthexplorer.usgs.gov) for a total of 4383 Landsat scenes acquired between 1986 and 2016 in mountainous regions of the conterminous western United States (Table 5.2). We acquired the 30 m resolution National Land Cover Database (NLCD) 2011 land cover and percent forest canopy datasets for the conterminous U.S. (available from http://www.mrlc.gov) and then extracted subsets covering our study areas. Each subset was reprojected to the UTM projection associated with Landsat scenes within the study area. We obtained 1/3 arc second (approximately 10 m resolution) digital elevation models (DEMs) covering each of our study areas from the U.S. 3D Digital Elevation (3DEP) program (available from https://viewer.nationalmap.gov/basic/). Individual 1° tiles were mosaicked together, reprojected to the UTM projection associated with Landsat scenes within the study area, and then aggregated from 10 m to 30 m spatial resolution to correspond with the resolution of the Landsat scenes and ancillary data. 5.2.4.1 Image Processing We retrieved fSCA for each pixel from each Landsat scene using TMSCAG, a spectral mixture analysis model that evolved from the original MEMSCAG algorithm 87 that calculated fSCA and snow grain size from optical imaging spectrometer data (Painter et al., 1998; Painter et al., 2003). The TMSCAG model is similar to the MODSCAG model (Painter et al., 2009), which works with multispectral MODIS data and has been widely used for retrieval of fSCA; the key difference is that the TMSCAG model is configured to handle radiometric saturation in spectral bands 1-4. A more detailed description of the TMSCAG model as well as validation of model performance is provided in Painter et al. (in review). Potential solar radiation grids were calculated at 30 m spatial resolution using the r.sun algorithm available in GRASS GIS and the resampled 30 m DEM covering the study area. For each study area, we calculated potential solar radiation for every 10th day starting with day of year 274 (October 1 for regular years, September 30 for leap years). Potential solar radiation was interpolated for days in between, resulting in a daily time series of potential solar radiation. We then calculated cumulative potential solar radiation since October 1 for each day of the year. We used a modification of the CFmask cloud masking approach (Zhu & Woodcock, 2012) described in Selkowitz and Forster (2015) to identify cloudcovered pixels in each Landsat scene. While the original CFmask has been demonstrated to consistently classify certain landscape patches containing snow and ice or a mixture of snow/ice and rock as cloud cover, this problem is minimized in the revised version. In high mountain areas, accuracy for the original CFmask algorithm is 66%, while accuracy for the revised CFmask approach is 88% (Selkowitz & Forster, 2015). 88 Figure 5.2 provides a flow chart that describes the full canopy adjustment process. In cases where retrieved fSCA is > 0 and < 1, we incorporate the retrieved fSCA value as well as the forest canopy fraction value from the National Land Cover Database (NLCD) to calculate an adjusted fSCA value. In previous work, several authors (Coons et al., 2014; Durand & Molotch, 2008; Molotch & Margulis, 2008; Raleigh et al., 2013) have used a similar approach that normalizes retrieved fSCAv by the noncanopy fraction of the pixel (1-Fc) to calculate an adjusted fSCAadj value. This approach is defined in Equation 1: fSCA!"# = 𝑚𝑖𝑛 !"#$! !!!! , 1.0 (1) In our approach, we use the NLCD canopy percent value as Fc. In addition, we add 0.35 Fc to the result calculated using equation 1. We added this term to equation 1 because the NLCD canopy dataset tends to underestimate forest canopy by an average of 9.7% (and by as much as 23.4% in the Sierra Nevada) (Nowak & Greenfield, 2010), and because the viewable snow fraction in areas with forest cover tends to be poorly illuminated due to shading from the canopy, often leading to underestimation of viewable fSCA. Our canopy adjustment approach is defined in Equation 2: fSCA!"# = 𝑚𝑖𝑛 !"#$! !!!! + 0.35 𝐹𝑐, 1.0 (2) 89 In cases where retrieved fSCA is 0, we implement an additional approach designed to detect snow-covered pixels where forest canopy or a combination of forest canopy and shading would make snow cover detection otherwise impossible using optical remote sensing. We refer to this approach as the neighborhood canopy adjustment approach. In many cases where snow cover is present but not initially detected by TMSCAG, the viewable snow cover fraction is below the TMSCAG detection threshold of approximately 0.15. The neighborhood canopy adjustment approach (defined below) relies on the examination of at least 10 surrogate pixels located near the target pixel. The ratio of surrogate pixels with snow cover to total surrogate pixels determines whether snow is classified at the target pixel. This process is conducted in two separate phases, with the second phase only implemented if necessary. In the first phase, surrogate pixels are identified within a 11 x 11 pixel neighborhood centered on the target pixel. To qualify as a surrogate pixel, NLCD forest canopy percent must be lower than at the target pixel, potential solar radiation must be greater than at the target pixel, and elevation must be no more than 75 m greater than at the target pixel (Table 5.3, Figure 5.3). In addition, to qualify as a surrogate pixel a pixel must also contain valid data (i.e., not be within a Landsat 7 scan line gap) and be cloudfree, according to the revised cloud cover mask described above. If at least 10 surrogate pixels can be identified in the 11 x 11 local window, we compute the ratio of surrogate pixels with snow cover (fSCA > 0) to the total number of surrogate pixels. If this ratio exceeds 0.3, the target pixel is labeled as snow covered and given an fSCA value of 0.15. In some cases, however, less than 10 90 Table 5.3. Criteria for identification of surrogate pixels in the surrounding 11 x 11 or 31 x 31 pixel neighborhoods. Criteria Rule for surrogate eligibility NLCD canopy < target pixel AND < 60% Potential solar radiation >= target pixel Elevation < target pixel + 75 Cloud cover Must be cloud-free, and valid (not in SLC gap) surrogate pixels can be identified within the 11 x 11 pixel window. This can occur when the target pixel has unique canopy or topographic aspects when compared to nearby pixels, when all nearby pixels are covered by dense forest, or when the number of available surrogate pixels is reduced due to the presence of scan line corrector gaps or cloud cover. In cases where < 10 surrogate pixels can be identified within the 11 x 11 local window, we initiate a second phase that examines a larger 31 x 31 pixel window with a larger pool of potential surrogate pixels. The criteria for identification of surrogate pixels in the second phase are the same as the criteria used in the first phase (Table 5.3). If at least 15 surrogate pixels can be identified in the 31 x 31 pixel window, we compute the ratio of snow covered surrogate pixels to total surrogate pixels in the same manner as the first phase. For the larger 31 x 31 pixel window, if the ratio of snow covered pixels to total surrogate pixels exceeds 0.45, the target pixel is labeled as snow covered and given an fSCA value of 0.15. The size of neighborhoods used for identification of potential surrogate pixels was chosen to balance the need for a sufficient sample of surrogate pixels for decision making with the need to constrain potential surrogate pixels to those pixels 91 Figure 5.3. Identification of surrogate pixels meeting criteria listed in Table 5.2 in the 11 x 11 pixel neighborhood surrounding target pixel. The same process for identification of surrogate pixels is used for the larger 31 x 31 pixel window used in phase 2 (if necessary). with similar climatic characteristics to the target pixel. This excludes solar radiation, which can vary substantially within the 11 x 11 or 31 x 31 pixel neighborhoods, but is explicitly accounted for. Initially, we attempt to identify at least 10 surrogate pixels within 150 m of the target pixel because, given the solar radiation and elevation constraints, the closest pixels are likely to exhibit snow cover conditions most similar to the target pixel. When insufficient surrogate pixels are available within the smaller neighborhood, we examine the larger neighborhood for surrogate pixels. The algorithm is designed to be conservative when identifying missed snow cover, and to be especially conservative when missed snow cover is identified using the larger 31 x 31 pixel neighborhood. For this reason, we require a 92 higher ratio of snow-covered to total surrogate pixels to identify additional snow covered pixels in the larger 31 x 31 pixel neighborhood. Once the first, and, if necessary, second phases of neighborhood canopy adjustment have been implemented, standard canopy adjustment described in Equation 2 is applied to the resulting fSCA value of 0.15 (if snow cover is determined to be present) at the target pixel. 5.2.5 Calculation of Mean Annual Snow Cover Duration For the Landsat snow cover duration product, we exploit the historical Landsat archive by incorporating fSCA and cloud cover calculated for all scenes acquired during the period of interest. Using these data, we compute the ratio of snow covered days for all cloud-free surface views to the total number of cloud-free surface views for each calendar month (e.g., all cloud-free surface views acquired during the month of June over the period 1986-2015) for each 30 m pixel (Figure 5.4). For this calculation, all fSCA values > 0 are counted as snow cover. The monthly ratio of snow covered days to total cloud-free days is then weighted by number of days in the month to compute the fraction of days with snow cover for the entire year, which is then multiplied by 365 to convert to mean annual snow cover duration in units of days. Calculation of individual monthly ratios which are then converted to annual snow cover days is preferable to simply computing the ratio of snow-covered days to total cloud-free days for the entire period because at many locations, more cloud- 93 Figure 5.4. Calculation of the fraction of days in June with snow cover for the period 1986-2015 for a single 30 m grid cell. The matrix indicates 1 of 4 potential outcomes for each day in the period 1986-2015. Light grey squares indicate days when no Landsat data were acquired, dark grey squares indicate days when Landsat was acquired, but was not used for calculation of snow cover duration statistics due to cloud cover or missing data, red squares indicate snow free land, and light blue squares indicate snow-cover. The ratio of snow-covered days to snow-covered and snow-free days (total cloud-free days) is used to calculate the snow cover days fraction. 94 free views are available during the summer months than during the winter months. Consequently, computing the ratio of snow cover days to total cloud-free days for the entire period would result in underestimation of snow cover days because a disproportionate number of cloud-free views usually come from the summer months where snow cover is much less common. 5.2.6 Validation Approach We use separate datasets to provide an accuracy assessment for the scene- based canopy-adjusted fSCA products and the snow cover duration product. Validation data for canopy-adjusted fSCA from individual Landsat scenes consisted of in situ snow cover fraction data collected at sites across the Sierra Nevada of California, while validation data for the snow cover duration datasets consisted of data from SNOTEL sites maintained by the Natural Resources Conservation Service (NRCS) at locations in California, Oregon, Washington, and Wyoming (Figure 5.1). The two separate datasets were chosen for validation of the two products because the in situ snow cover fraction data allowed for accuracy assessment of fractional snow covered area, which was not possible with SNOTEL data. The in situ snow cover fraction dataset, however, covered only the period 2014-2016 and thus could not be used for validation of mean annual snow cover duration calculated over longer periods. The SNOTEL dataset was selected for validation of the snow cover duration product because this was the only dataset available covering a 30-year period of record with sites covering multiple mountain ranges in the western United States. 95 5.2.6.1 Validation of Individual Scenes Using In Situ Sensor Arrays We deployed 100 x 100 m temperature data arrays at 27 forest-covered sites across the Sierra Nevada of California. Sites ranged in elevation from 1860 m to 2930 m, included a wide range of forest canopy densities (4-83%), and included flat, gentle, and moderately steep (up to 23°) slopes on all aspects. Arrays consisted of 6 x 6 sensors spaced at 20 m intervals (2014-2015) or 5 x 5 sensors spaced at 25 m intervals (2015-2016). At each site, an array of temperature data loggers were buried 2-5 cm below the soil surface. Sensors were set to record temperature at 1.5 hour intervals. Temperature data loggers were deployed between August and November each year and retrieved between April and August of each year. We used the algorithm designed by Raleigh et al. (2013) and adapted by Selkowitz et al. (2014) to convert hourly or 1.5 hourly temperature time series from individual temperature data loggers to daily snow cover fraction for each 100 m grid cell footprint we monitored. This algorithm classifies snow cover if temperature varies by less than 1°C at an individual temperature data logger over two consecutive 24-hour periods. A more detailed discussion of this approach can be found in Raleigh et al. (2013), Selkowitz et al. (2014), and Lundquist and Lott (2008). The number of temperature data loggers used to compute daily snow cover fraction varied from 15-34. Although either 25 or 36 data loggers were installed at each site, some data loggers malfunctioned or stopped recording due to insufficient battery power and some data loggers could not be located (often as a result of ground disturbance by marmots or ground squirrels). 96 In order to verify that in situ temperature data loggers could accurately monitor the presence/absence of snow cover, we observed the presence/absence of snow cover at the location of the deployed temperature data loggers at four sites for a total of 5 days (one site was surveyed twice) in the spring and early summer of 2016. We used a GPS unit that provided real-time accuracy of +- 1.5 m or better to navigate to the location of temperature data loggers and record snow cover presence or absence above each temperature data logger. We collected 106 visual snow cover presence/absence observations for comparison with snow cover presence/absence classified using hourly temperature data from the data loggers using the algorithm described above. 5.2.6.2 Validation of Mean Annual Snow Cover Duration for 30-year and 5-year Periods Using SNOTEL Data In order to assess the accuracy of the snow cover duration product, we used data from 72 SNOTEL sites from California, Oregon, Washington, and Wyoming. While SNOTEL pillow measurements of snow water equivalent cover < 2 m2, a much smaller area than the 900 m2 covered by the nominal Landsat pixel size, they are one of the only long term measurements of snow cover in mountains regions of the western U.S. Therefore, despite the mismatch in area monitored by a SNOTEL pillow compared to a Landsat pixel, SNOTEL sites still represent the best source of data for validation of mean annual snow cover duration over decades. For our analysis, we calculated the mean annual snow cover duration (days with SWE > 0) observed at each SNOTEL station for the 30-year period 1986-2015 97 and compared this to the mean annual snow cover duration calculated from Landsat using the methods described above. We also conducted the same comparison between SNOTEL and Landsat-derived mean annual snow cover duration for the five year periods 1991-1995, 1996-2000, 2001-2005, 2006-2010, and 2011-2015. The period 1986-1990 was excluded because relatively few Landsat scenes were available during this period. 5.3 Results 5.3.1 Validation of In Situ Temperature Data Logger Snow Cover Monitoring Approach Assessment of snow cover classification using in situ temperature data loggers compared to visual observations of snow cover in May and June of 2016 indicated agreement in 102 out of 106 cases, with one false positive and three false negatives (96% accuracy). 5.3.2 Landsat-derived fSCA Compared to In Situ fSCA Comparison of TMSCAG and TMSCAG canopy adjusted fSCA to in situ fSCA from temperature data logger arrays indicated the canopy adjustment approaches used here substantially improved agreement between Landsat-derived and in situ measured fSCA (Table 5.4, Figure 5.5). While the standard canopy adjustment approach alone resulted in a substantial increase in accuracy, reducing RMSE from 0.49 to 0.25, the neighborhood adjustment approach resulted in further 98 Table 5.4. Accuracy metrics for TMSCAG, canopy adjusted TMSCAG, and neighborhood canopy adjusted TMSCAG. Metric TMSCAG Adjusted Adjusted + Neighborhood Canopy Adjustment RMSE 0.49 0.25 0.20 Mean Error (Bias) -0.38 -0.11 -0.07 Binary Accuracy 0.94 0.94 0.96 Binary Precision 0.99 0.99 0.99 Binary Recall 0.95 0.95 0.97 Binary F 0.97 0.97 0.98 improvement, reducing RMSE to 0.20. The mean error (bias) was initially -0.38 for the unadjusted TMSCAG data but was reduced to -0.07 when both adjustment approaches were applied. In the 139 instances where in situ fSCA was compared with TMSCAG and canopy-adjusted TMSCAG fSCA, there was only a single false positive snow cover result. There were no false positives resulting from addition of snow cover via the neighborhood canopy adjustment approach. It is important to note that the full impact of the neighborhood canopy adjustment approach is less evident in this comparison because fSCA from 9 Landsat pixels is compared to in situ fSCA across a 100 x 100 m grid, reducing the impact of individual 30 m Landsat pixels where snow cover is missed on the overall accuracy. An example of TMSCAG fSCA compared to fSCA adjusted using both the standard canopy adjustment approach and neighborhood canopy adjustment approach is shown in Figure 5.6, which also maps the spatial distribution of additional snow cover pixels added using the neighborhood canopy adjustment approach. 99 Figure 5.5. TMSCAG canopy-adjusted fSCA compared to in situ fSCA calculated from temperature data logger arrays. Both the standard canopy adjustment approach and the neighborhood canopy adjustment approach were applied to produce canopyadjusted fSCA. 5.3.3 Mean Annual Snow Cover Duration Mean annual snow cover duration computed from canopy-adjusted Landsat- derived fSCA for three 30 x 30 km subsets and three corresponding areas of detail are shown in Figure 5.7. The effect of elevation on snow cover duration is evident at the broader scale for the 30 x 30 km subsets, while the effects of both elevation and topographic position are evident at the finer scale shown for the 5 x 5 km areas of detail. The spatial distribution of additional snow cover days added via the canopy adjustment algorithm is shown for the three 30 x 30 km subsets in Figure 5.8. The three maps of additional snow cover days added via neighborhood canopy adjustment and the corresponding cumulative histograms demonstrate that impact of the canopy adjustment algorithm was largest in the Cascades subset, where 30 or more days were added to 27% of all pixels. The impact of the canopy adjustment 100 Figure 5.6. Demonstration of canopy adjustment for an area in the northern Sierra Nevada on April 20, 2009: (a) Landsat surface reflectance (bands 7-4-2), (b) TMSCAG fSCA, (c) canopy adjusted TMSCAG fSCA, and (d) areas of snow cover added using the neighborhood canopy adjustment approach. Figure 5.7. Mean annual snow cover duration for three 30 x 30 km subsets: (a) Cascades (Washington), (b) Sierra Nevada (California), and (c) Gros Ventre (Wyoming). 101 102 Figure 5.8. Additional snow cover days added using the neighborhood canopy adjustment approach for three 30 x 30 km subsets: (a) Cascades (Washington), (b) Sierra Nevada (California), and (c) Gros Ventre (Wyoming). The dashed line corresponds to 30 days added via the neighborhood canopy adjustment approach. algorithm was substantially less for the Sierra Nevada subset, where 30 or more days were added to just 9% of pixels, and for the Gros Ventre subset, where 30 more days were added to only 5% of all pixels. The percent of all cloud-free pixels where snow cover was added via neighborhood canopy adjustment, shown by month (Figure 5.9) tells a similar story, with snow cover added for substantially more instances in the Cascades subset than in the Sierra Nevada subset or Gros Ventre subset. Figure 5.9 also indicates how the impact on mapped SCA varies seasonally. For all three subsets, the percentage of 103 Figure 5.9. Percent of all cloud-free pixels where snow cover was added via canopy adjustment, shown by month. snow cover pixels added is high during December and January and low during the summer and early fall months. In the Cascades, and to a lesser extent in the Sierra Nevada, a secondary peak in the percentage of added snow cover pixels occurs in April and May. The neighborhood canopy adjustment approach occasionally fails due to insufficient surrogate pixels within both the 11 x 11 and 31 x 31 pixel neighborhoods. The incidence of failure is determined ultimately by not only the prevalence of surrogate pixels within the local neighborhood, but, on a scene-byscene basis, whether or not those potential surrogate pixels are obscured by clouds or scan-line correction failure gaps. The spatial distribution of failure frequency (as a fraction of total valid cloud-free days at each pixel) is shown in Figure 5.10, along with cumulative histograms that indicate the percent of pixels affected by various model failure rates. These data suggest that failure of the neighborhood canopy adjustment approach is more common in the Cascades subset than in the other two subsets, where 17% of pixels experienced failure in more than 10% of instances. By 104 Figure 5.10. Canopy adjustment model failure frequency for three 30 x 30 km subsets: (a) Cascades (Washington), (b) Sierra Nevada (California) and, (c) Gros Ventre (Wyoming). contrast, in the Sierra Nevada and Gros Ventre subsets, only about 2% of pixels experienced model failure in more than 10% of instances. Comparison between mean annual snow cover duration calculated from SNOTEL stations and from the original TMSCAG and canopy-adjusted TMSCAG image data allow for quantification of the improvement in accuracy achieved via the neighborhood canopy adjustment approach. This comparison is based on binary snow cover retrievals, where the total fraction of days with fSCA > 0 is compared to the number of days where the corresponding SNOTEL site recorded SWE > 0. It is important to note that the standard canopy adjustment approach described in Equation 2 therefore has no impact on accuracy in this comparison. This is because 105 the standard canopy adjustment approach only adjusts pixels where a snow cover fraction > 0 has been detected. Consequently, differences between the unadjusted and adjusted results for mean annual snow cover duration are due entirely to the implementation of the neighborhood canopy adjustment approach. While the unadjusted results are reasonably accurate for sites with canopy cover < 50%, mean annual snow cover duration is severely underestimated at many sites with canopy cover > 50% (Figure 5.11a). When fSCA from individual scenes is adjusted using the neighborhood canopy adjustment approach, however, agreement between mean annual snow cover duration from SNOTEL sites and from Landsat improves substantially at all but one site (Figure 5.11b). The local window canopy adjustment approach reduces RMSE from 22.6 days to 14.7 days and essentially eliminates the negative bias in mean snow cover duration (Table 5.5). Figures 5.11c and 5.11d also indicate that neighborhood canopy adjustment has a much larger impact at sites in the Cascades than at sites in the Sierra Nevada or Rocky Mountains. Agreement between TMSCAG canopy-adjusted snow cover duration and snow cover duration calculated from SNOTEL sites for 5-year periods is slightly lower than for the full 30 1986-2015 period, with RMSE ranging from 15.4 - 20.7 (Table 5.6). However, the improvement in accuracy resulting from the neighborhood canopy adjustment approach is still evident, and the negative bias for the canopy adjusted results is < 3 days for 4 of the 5 periods considered. 106 Figure 5.11. Comparison between mean annual snow cover duration calculated from SNOTEL data and mean annual snow cover duration calculated from Landsat. (a) unadjusted snow cover duration, with colors indicating forest canopy density (b) adjusted snow cover duration, with colors indicating forest canopy density, (c) unadjusted snow cover duration, with colors indicating region, and (d) adjusted snow cover duration, with colors indicating region. Table 5.5. Accuracy metrics for mean annual snow cover duration (days) calculated using unadjusted TMSCAG and canopy adjusted TMSCAG relative to mean annual snow cover duration calculated from SNOTEL sites. Metric TMSCAG Canopy Adjusted RMSE 22.6 14.7 Mean Error (Bias) -10.4 0.6 107 Table 5.6. Accuracy metrics for mean annual snow cover days calculated using adjusted TMSCAG for periods 1991-1995, 1996-2000, 2001-2005, 2006-2011, and 2011-2015. Metric 1986 - 1991 - 1996 20012006 - 20112015 1995 2000 2005 2010 2015 RMSE 15.1 20.0 15.4 18.1 19.3 20.7 Mean Error 0.5 -1.3 -5.0 -2.3 0.0 -0.4 (Bias) 5.4 Discussion TMSCAG has been demonstrated effective for retrieval of visible fSCA across a wide range of snow cover conditions, topography, vegetation types, and solar illumination conditions (Painter et al., in review). When significant forest canopy is present, however, the difference between the retrieved viewable fSCA and in situ fSCA beneath the canopy can be significant. In cases where snow cover is missed entirely, this can also impact mean annual snow cover duration calculated using all available scenes for a period of record. The results presented here indicate that TMSCAG fSCA mapping can be extended to allow for effective retrievals of in situ fSCA under forest canopies across the forests of the western conterminous U.S. mountains, and likely in other regions as well. The canopy adjustment approach presented here not only adjusts viewable fSCA values > 0, but also adds snow to forested pixels where snow cover was initially not retrieved if conditions at surrounding pixels with similar characteristics suggest that snow cover was likely missed. The adjustment of initial fSCA values > 0 108 has been applied previously in several studies (Coons et al., 2014; Durand & Molotch, 2008; Molotch & Margulis, 2008; Raleigh et al., 2013). Addition of snow cover to forested pixels initially identified as snow-free, however, is a novel approach critical for effectively monitoring snow cover conditions in forested regions at the Landsat spatial resolution. Comparison between TMSCAG fSCA and in situ fSCA in forested areas of the Sierra Nevada indicated that, as expected, the viewable snow cover fraction retrieved from TMSCAG is usually substantially lower than the in situ fSCA. While the standard canopy adjustment approach improved agreement between in situ fSCA and Landsat-derived fSCA in many instances, the neighborhood canopy adjustment approach resulted in further improvement in cases where snow cover was initially missed at some of the 9 Landsat pixels covering the 100 x 100 m in situ grid. The lack of false positives for snow cover resulting from application of the neighborhood canopy adjustment approach suggest that this is a relatively conservative approach to the problem of missed snow cover beneath forest canopy. In fact, the algorithm is structured so that snow cover can only be added if a substantial fraction of pixels within 450 m of the target pixel have been identified as snow covered by TMSCAG. Checks are also in place to prohibit the addition of snow cover when nearby snow cover is only identified at pixels with a lower cumulative solar radiation load or at substantially higher elevations. Despite these conditions and the lack of false positives for snow cover identified in the in situ fSCA 109 comparison dataset, we acknowledge that some false positives will inevitably result from this canopy adjustment approach. Our results suggest the importance of the neighborhood canopy adjustment varies both by region and over time, with frequent instances of adjustment to add snow cover initially missed concentrated in both clusters of pixels and specific times of year. The substantially higher amount of snow covered pixels added in the Cascades from neighborhood canopy adjustment relative to the Sierra Nevada and Gros Ventre subsets can be explained primarily by the higher forest density and greater prevalence of dense forest in the Cascades. The canopy adjustment approach may also be less necessary in mid-winter in the colder continental climate of the Gros Ventre subset, where snow cover is often retained for long periods in the canopy (Hedstrom & Pomeroy, 1998), resulting in a higher viewable fSCA and thus fewer missed snow cover pixels. The higher percentage of added snow cover pixels occurring in the months of November, December, and January for all three subsets can be explained by the relatively poor solar illumination conditions resulting from higher solar zenith angles during these months. Painter et al. (in review) found that TMSCAG is more likely to underestimate snow cover under poor illumination conditions, and this is likely exacerbated by forest canopy that leads to a further reduction in illumination at the snow surface. Under these conditions, which result in a higher frequency of pixels where snow cover is not initially identified by TMSCAG, the neighborhood canopy adjustment approach is employed more frequently. The secondary peak of added snow cover pixels during the months of April and May in the Cascades and 110 Sierra Nevada is likely due to the higher prevalence of partially snow covered pixels during this period, which typically corresponds with snowmelt across much of these two study areas. The density of forest canopy cover necessary to result in errors of omission declines substantially as the ground snow cover fraction declines. The comparison between mean annual snow cover duration calculated from SNOTEL data and from Landsat-derived fSCA confirms the effectiveness of the neighborhood canopy adjustment approach. Without an approach that considers nearby pixels or some sort of ancillary data, mean annual snow cover duration calculated from Landsat-derived fSCA is significantly underestimated at many pixels with moderate to dense forest canopy. While the canopy adjustment approach described here is not effective for eliminating all errors of omission at all Landsat pixels in all instances, it significantly reduces underestimation of mean snow cover duration at many forested pixels. Perhaps the largest limitation for the current approach to production of both scene-based and mean annual snow cover duration products is the inability of the neighborhood canopy adjustment approach to effectively correct all pixels. While instances of failure due to insufficient surrogate pixels are relatively rare, they are typically concentrated in both space and time and can therefore have a notable impact on results in certain areas and for certain periods. It may be possible to reduce the number of cases where the neighborhood canopy adjustment approach fails by extending the size of the neighborhood for identification of surrogate pixels. A completely effective solution to this limitation, however, will likely require the 111 inclusion of different remote sensing data such as lidar or a physically-based modeling approach. Another limitation of both the scene-based canopy-adjusted fSCA and the snow cover duration products is that the accuracy of canopy-adjusted snow covered area depends on accurate and consistent forest canopy information. Presently, this information is provided by the NLCD 2011 forest canopy layer, which has been shown to underestimate forest canopy by an average of 9.7% nationally and by 23.4% in the Sierra Nevada (Nowak & Greenfield, 2010). Perhaps more importantly, the NLCD canopy layer represents only a brief period of time and is likely to be incorrect in instances where land cover change has occurred either before or after the publication of the dataset. For example, a large fire might remove most or all of the canopy for a patch of pixels. If the reduction in canopy from the fire is not reflected in the canopy layer, canopy adjustment will be incorrectly applied to this patch of pixels, resulting in an overestimation of fSCA and possibly the generation of false positive snow cover pixels (although only if other nearby pixels are snow covered). A possible solution would be to generate new canopy layers annually or possibly even a new canopy layer for every Landsat scene processed in order to reduce the possibility of these types of errors. The relatively long 16-day interval between scene acquisitions from Landsat also limits our ability to produce a snow cover duration product for time periods shorter than about 5 years. The sporadic occurrence of cloud cover and the 16-day repeat interval for the Landsat 5 and 7 spacecraft have resulted in an irregular frequency of cloud-free surface views. Over periods longer than approximately 5 112 years, this irregular availability of cloud-free surface views amounts to a semirandom sample of cloud-free surface views from the early, middle and late portions of each month and from above-average, average, and below-average snow cover years. For shorter periods, however, the impact of the timing of each cloud-free surface view can skew the calculation of mean annual snow cover duration. Our results indicate that accuracy of the snow cover duration product is lower for 5-year periods than it is for the full 1986-2015 period. Future iterations of the Landsat snow cover products will incorporate data from Landsat 8 and possibly from the Sentinel 2A instrument, potentially allowing for generation of snow cover duration products based on as little as one year of data. 5.5 Conclusions Results presented here indicate that while the unadjusted Landsat snow cover products underestimate fSCA for individual scenes and underestimate mean annual snow cover duration calculated from all available scenes, this underestimation is substantially reduced for the canopy-adjusted versions of these products. The incorporation of a canopy adjustment approach that considers the snow cover status of nearby pixels allows for accurate estimation of scene-based fSCA in many cases even when the initially retrieved fSCA value is zero. When combined with a cloud mask optimized for use in mountainous environments, the resulting canopy-adjusted Landsat fSCA data can be used to provide an accurate estimate of mean annual snow cover duration at 30 m spatial resolution for the entire Landsat period of record as well as for temporal subsets as short as 5 years. 113 Data from the SNOTEL network indicate the mean annual snow cover duration for the period 1986-2015 calculated using this approach has an RMSE of 14.7 days and a bias of +0.6 days over the range of 90-270 days of annual snow cover. While the impact of the canopy adjustment approach varies by region due primarily to differences in forest cover, it can have a large impact on local snow cover duration estimates for pixels with forest canopy even within regions where forest cover is relatively sparse or mostly absent. The primary limitation of the canopy adjustment approach is that it is ineffective in cases where an insufficient number of surrogate pixels can be located within the local neighborhood, such as when a target pixel is surrounded by large tracts of contiguous, high density forest cover. Despite this limitation, the canopy adjustment approach substantially increases accuracy of fSCA maps in forested and partially forested regions. Together, the improved accuracy of the scene-based fSCA product and the approach developed to incorporate all Landsat data acquired during a period of record to calculate mean annual snow cover duration enables the production of an accurate snow cover duration product at a higher spatial resolution than has previously been available. The relatively high spatial resolution of both the Landsat scene-based snow cover product and the Landsat snow cover duration product can help illuminate relatively fine scale snow cover patterns common in rugged topography that would be obscured by the coarser spatial resolution of sensors like MODIS or VIIRS. This will result in enhanced understanding and possibly new insights into snow cover patterns and processes, particularly in regions with complex topography that consistently feature a high degree of fine scale snow cover variability. 114 5.6 Acknowledgements Funding for this research was provided by the Land Remote Sensing Program of the US Geological Survey. The authors gratefully acknowledge assistance in the field provided by Melanie Cota, Mike Stockwell and Erin Orozco. Part of this work was performed at the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA 5.7 References Anderton, S.P., White, S.M., & Alvera, B. (2002). Micro-scale spatial variability and the timing of snow melt runoff in a high mountain catchment. Journal of Hydrology, 268, 158-176. Aubry, K.B., Mckelvey, K. S., & Copeland, J. P. (2007). Distribution and broadscale habitat relations of the wolverine in the contiguous United States. The Journal of Wildlife Management, 71(7), 2147-2158. Barnett, T.P., Adam, J.C., & Lettenmaier, D.P. (2005). Potential impacts of a warming climate on water availability in snow-dominated regions. Nature, 438(17), 303-309. Beck, P. S., Kalmbach, E., Joly, D., Stien, A., & Nilsen, L. (2005). Modelling local distribution of an Arctic dwarf shrub indicates an important role for remote sensing of snow cover. Remote Sensing of Environment, 98(1), 110-121. Bekker, M. F. (2005). Positive feedback between tree establishment and patterns of subalpine forest advancement, Glacier National Park, Montana, USA. Arctic, Antarctic, and Alpine Research, 37(1), 97-107. Bernhardt, M., Liston, G.E., Strasser, U., Zangl, G., & Schulz, K. (2010). High resolution modeling of snow transport in complex terrain using downscaled MM5 wind fields. The Cryosphere, 4, 99-113. Bernhardt, M., & Schulz, K. (2010). SnowSlide: A simple routine for calculating gravitational snow transport. Geophysical Research Letters, 37(11). Bloschl, G., Kirnbauer, R., & Gutknecht, D. (1991). Distributed snowmelt simulations in an alpine catchment 1. Model evaluation on the basis of snow cover patterns. Water Resources Research, 27(12), 3171-3179. 115 Billings, W.D., & Bliss, L.C. (1959). An alpine snowbank environment and its effects on vegetation, plant development, and productivity. Ecology, 40(3), 388-397. Bronge†, L.B., & Bronge†*, C. (1999). Ice and snow-type classification in the Vestfold Hills, East Antarctica, using Landsat-TM data and ground radiometer measurements. International Journal of Remote Sensing, 20(2), 225-240. Cline, D.W., Bales, R.C., & Dozier, J. (1998). Estimating the spatial distribution of snow in mountain basins using remote sensing and energy balance modeling. Water Resources Research, 34(5), 1275-1285. Cohen, J., & Entekhabi, D. (1999). Eurasian snow cover variability and Northern Hemisphere climate predictability. Geophysical Research Letters, 26(3), 345-348. Coons, L. P., Nolin, A.W., Gleason, K.E., Mar, E.J., Rittger, K., Roth, T.R., & Painter, T.H. (2014). Seeing the snow through the trees: Toward a validated canopy adjustment for satellite snow-covered Area. In V. Lakshmi (Ed.) Remote Sensing of the Terrestrial Water Cycle. Hoboken, NJ: John Wiley and Sons. Deems, J.S., Fassnacht, S.R., & Elder, K. J. (2006). Fractal distribution of snow depth from LiDAR data. Journal of Hydrometeorology, 7(2), 285-297. Dozier, J. (1984). Snow reflectance from LANDSAT-4 Thematic Mapper. IEEE Transactions on Geoscience and Remote Sensing, GE-22(3), 323-328. Durand, M., Molotch, N.P., & Margulis, S.A. (2008). Merging complementary remote sensing datasets in the context of snow water equivalent reconstruction. Remote Sensing of Environment, 112: 1212-1225. Eastland, W.G., Bowyer, R.T., & Fancy, S.G. (1989). Effects of snow cover on selection of calving sites by caribou. Journal of Mammology, 70(4), 824-828. Fily, M., Bourdelles, B., Dedieu, J.P., & Sergent, C. (1997). Comparison of in situ and Landsat Thematic Mapper derived snow grain characteristics in the Alps. Remote Sensing of Environment, 59, 452-460. Fily, M., Dedieu, J.P., & Durand, Y. (1999). Comparison between the results of a snow metamorphism model and remote sensing derived snow parameters in the Alps. Remote Sensing of Environment, 68, 254-263. Groffman, P.M., Driscoll, C.T., Fahey, T.J., Hardy, J.P., Fitzhugh, R.D., & Tierney, G.L. (2001). Colder soils in a warmer world: A snow manipulation study in the northern hardwood forest ecosystem. Biogeochemistry, 56, 135-150. 116 Groisman, P.Y, Karl, T.R, & Knight, R.W. (1994). Observed impact of snow cover on the heat balance and the rise of continental spring temperatures. Science, 263, 198200. Hall, D.K., Chang, A.T.C., Foster, J.L., Benson, C.S., & Kovalick, W.M. (1989). Comparison of in situ and Landsat derived reflectance of Alaskan glaciers. Remote Sensing of Environment, 28, 23-31. Hall, D.K., Riggs, G.A., & Salomonson, V.V. (1995). Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer data. Remote Sensing of Environment, 54, 127-140. Hedstrom, N. R., & Pomeroy, J. W. (1998). Measurements and modelling of snow interception in the boreal forest. Hydrological Processes, 12(1011), 1611-1625. Hessl, A.E., & Baker, W.L. (1997). Spruce and fir regeneration and climate in the forest-tundra ecotone of Rocky Mountain National Park, Colorado, USA. Arctic and Alpine Research, 29(2), 173-183. Huggard, D.J. (1993). Effects of snow depth on predation and scavenging by gray wolves. The Journal of Wildlife Management, 57(2): 382-388. Klein, A.G., Hall, D.K., & Riggs, G.A. (1998). Improving snow cover mapping in forests through the use of a canopy reflectance model. Hydrological Processes, 12: 17231744. Klein, A.G., & Isacks, B.L. (1999). Spectral mixture analysis of Landsat thematic mapper images applied to the detection of the transient snowline on tropical Andean glaciers. Global and Planetary Change, 22, 139-154. Laperriere, A.J., & Lent, P.C. (1977). Caribou feeding sites in relation to snow characteristics in northeastern Alaska. Arctic, 10, 101-108. Letsinger, S.L., & Olyphant, G.A. (2007). Distributed energy-balance modeling of snow-cover evolution and melt in rugged terrain: Tobacco Root Mountains, Montana, USA. Journal of Hydrology, 336(1-2), 48-60. Luce, C.H., & Tarboton, D.G. (1998). The influence of the spatial distribution of snow on basin-averaged snowmelt. Hydrological Processes, 12, 1671-1683. Lundquist, J.D., & Dettinger, M.D., (2005). How snowpack heterogeneity affects diurnal streamflow timing. Water Resources Research, 41, W05007. Lundquist, J.D., & Lott, F. (2008). Using inexpensive temperature sensors to monitor the duration and heterogeneity of snow-covered areas. Water Resources Research, 44(4). 117 Magee, T.K., & Antos, J.A. (1992). Tree invasion into a mountain-top meadow in the Oregon Coast Range, USA. Journal of Vegetation Science, 3(4), 485-494. Margulis, S.A., Cortés, G., Girotto, M., & Durand, M. (2016). A Landsat-Era Sierra Nevada snow reanalysis (1985-2015). Journal of Hydrometeorology,17(4), 12031221. Martinec, J., & Rango, A. (1981). Areal distribution of snow water equivalent evaluated by snow cover monitoring. Water Resources Research, 5, 1480-1488. Masek, J.G., Vermote, E.F., Saleous, N.E., Wolfe, R., Hall, F.G., Huemmrich, K. F., & Lim, T. K. (2006). A Landsat surface reflectance dataset for North America, 19902000. IEEE Geoscience and Remote Sensing Letters, 3(1), 68-72. Mock, C. J., & Birkeland, K. W. (2000). Snow avalanche climatology of the western United States mountain ranges. Bulletin of the American Meteorological Society, 81(10), 2367. Molotch, N.P. (2009). Reconstructing snow water equivalent in the Rio Grande headwaters using remotely sensed snow cover data and a spatially distributed snowmelt model. Hydrological Processes, 23(7), 1076-1089. Molotch, N.P., & Margulis, S.A. (2008). Estimating the distribution of snow water equivalent using remotely sensed snow cover data and a spatially distributed snowmelt model: A multi-resolution, multi-sensor comparison. Advances in Water Resources, 31(11), 1503-1514. Nowak, D.J., & Greenfield, E.J. (2010). Evaluating the National Land Cover Database tree canopy and impervious cover estimates across the conterminous United States: A comparison with photo-interpreted estimates. Environmental Management, 46(3), 378-390. Painter, T.H., Rittger, K., & Selkowitz, D. (in review). Retrieval of fractional snowcovered area from Landsat Thematic Mapper data. Painter, T.H., Roberts, D.A., Green, R.O., & Dozier, J. (1998). The effect of grain size on spectral mixture analysis of snow-covered area from AVIRIS data. Remote Sensing of Environment, 65(3), 320-332. Painter, T.H., Dozier, J., Roberts, D.A., Davis, R.E., & Green, R.O. (2003). Retrieval of subpixel snow-covered area and grain size from imaging spectrometer data. Remote Sensing of Environment, 85(1), 64-77. Painter, T.H., Rittger, K., McKenzie, C., Slaughter, P., & Davis, R.E. (2009). Retrieval of subpixel snow covered area, grain size, and albedo from MODIS. Remote Sensing of Environment, 113(4), 868-879. 118 Quinton, W.L., Bemrose, R.K., Zhang, Y., & Carey, S.K., (2009). The influence of spatial variability in snowmelt and active layer thaw on hillslope drainage for an alpine tundra hillslope. Hydrological Processes, 23(18), 2628-2639. Rittger, K., Painter, T.H., & Dozier, J. (2013). Assessment of methods for mapping snow cover from MODIS. Advances in Water Resources, 51, 367-380. Rosenthal, W., & Dozier, J. (1996). Automated mapping of montane snow cover at subpixel resolution from the Landsat Thematic Mapper. Water Resources Research, 32(1), 115-130. Salomonson, V.V., & Appel, I., (2004). Estimating fractional snow cover from MODIS using the normalized difference snow index. Remote Sensing of Environment, 89(3), 351-360. Selkowitz, D. J., Forster, R. R., & Caldwell, M. K. (2014). Prevalence of pure versus mixed snow cover pixels across spatial resolutions in alpine environments. Remote Sensing, 6(12), 12478-12508. Selkowitz, D.J., & Forster, R.R. (2015). An automated approach for mapping persistent ice and snow cover over high latitude regions. Remote Sensing, 8(1), 16. Stenseth, N.C., Shabbar, A., Chan, K.S., Boutin, S., Rueness, E.K., Ehrich, D., & Jakobsen, K. S. (2004). Snow conditions may create an invisible barrier for lynx. Proceedings of the National Academy of Sciences of the United States of America, 101(29), 1063210634. Sweeney, J.M., & Sweeney, J.R. (1984). Snow depths influencing winter movements of elk. Journal of Mammalogy, 65(3), 524-526. Tinkham, W. T., Smith, A. M., Marshall, H. P., Link, T. E., Falkowski, M. J., & Winstral, A. H. (2014). Quantifying spatial distribution of snow depth errors from LiDAR using Random Forest. Remote Sensing of Environment, 141, 105-115. Trujillo, E., & Molotch, N.P. (2014). Snowpack regimes of the western United States. Water Resources Research, 50(7), 5611-5623. Vikhamar, D., & Solberg, R. Snow-cover mapping in forests by constrained linear spectral unmixing of MODIS data. (2003). Remote Sensing of Environment, 88, 309323. Walker, D.A., Halfpenny, J.C., Walker, M.D., & Wessman, C.A., (1993). Long-term studies of snow-vegetation interactions. Bioscience, 43(5), 287-301. Winther, J.G., & Hall, D.K. (1999). Satellite-derived snow coverage related to hydropower production in Norway: Present and future. International Journal of Remote Sensing, 20(15-16), 2991-3008. Zhang, T., (2005). Influence of the seasonal snow cover on the ground thermal regime: An overview. Reviews of Geophysics, 43: RG4002. 119 CHAPTER 6 CONCLUSIONS Work presented in Chapters 2-5 demonstrates the utility of Landsat and other moderate to high spatial resolution multispectral instruments for mapping both seasonal snow cover and persistent ice and snow cover such as glaciers and perennial snowfields. While the approaches described in Chapters 3-5 for mapping persistent ice and snow cover, snow covered area in forested areas, and mean snow cover duration have been validated, additional innovation and adjustments will likely lead to further improvements in product accuracy. There are two specific avenues of research with strong potential for improving Landsat-derived snow cover data products. First, snow covered area mapping algorithms should be extended to work with data from similar multispectral sensors such as the Operational Land Imager (onboard Landsat 8) and the European Space Agency's Sentinel-2 instrument. Extension of the basic algorithms, including the TMSCAG spectral unmixing algorithm, the forest canopy adjustment approach, and the approach for mapping persistent ice and snow cover will be relatively straightforward. However, comprehensive assessment of differences between products resulting from differences in the spectral and spatial resolutions of these sensors will be necessary. In particular, the reduced potential for radiometric 121 saturation in the visible bands provided by Landsat 8 and Sentinel-2 relative to Landsat 5/7 may have a substantial impact on the retrieval of snow covered area and could bias change analysis that incorporates data from both sensors if not explicitly addressed. Extension of the algorithms presented here to sensors with similar spatial resolutions and spectral bands will be essential for maximizing the available data for analysis, particularly in regions where cloud-free views of the earth surface are relatively scarce. Ideally, snow cover products will eventually incorporate data from multiple sensors with different capabilities. Combining Landsat-derived snow cover with data from optical remote sensing instruments with significantly different spatial resolutions (e.g., MODIS, VIIRS) offers the potential for providing snow covered area products with better temporal resolution. Finally, combination of data from Landsat, other types of instruments such as radar or lidar, and physicallybased snow cover modeling enables the estimation of snow water equivalent, perhaps the most sought-after snow metric. The science data products described in Chapters 3-5 have a wide array of potential uses across a variety of disciplines. These datasets can improve our understanding of basic snow processes as well as the variability of snow and ice cover in the recent past and into the future. While a comprehensive discussion of potential scientific questions these datasets could help answer is beyond the scope of this work, several of the most pressing questions and research applications that could benefit from Landsat-derived snow cover data are outlined below. 122 Landsat-derived snow cover datasets are particularly valuable for providing a comprehensive inventory of snow cover across the full range of elevations, slopeaspect combinations and vegetation types present throughout a region such as an individual mountain range. Remotely sensed snow cover data are also crucial for monitoring snow cover above the treeline, where in situ observations are typically sparse or nonexistent. While other remotely sensed snow cover products, such as the MODIS snow products, can also provide this type of comprehensive inventory, Landsat's higher spatial resolution provides a unique opportunity to assess the relationship between snow cover duration and landscape characteristics such as slope, aspect, topographic position, and vegetation type and density. The higher spatial resolution snow cover data can be used to assess changes in snow cover duration for specific landscape types or positions over time. For example, recent research suggests that under future warming scenarios, topographic effects will have a strong impact on the timing of snowmelt (and resulting streamflow), with areas subject to topographic shading potentially more resistant to earlier snowmelt brought on by warmer temperatures (Lundquist & Flint, 2006). Landsat-derived snow cover duration data have the potential to provide a detailed assessment of this hypothesis. In another example, assessment of the spatial patterns of snow cover duration can also provide insight into the physical processes impacting snow accumulation, redistribution, and accumulation. Arctic and alpine environments typically experience substantial redistribution of snow by wind transport, resulting in substantial snow cover heterogeneity at scales < 100 m (Liston, 1998; Pomeroy, 2004). In recent years, understanding and modeling physical processes like wind 123 redistribution of snow has been a top priority for snow researchers. Landsatderived patterns of snow cover duration can be used to validate and improve physically-based snow evolution models. These models will in turn be useful for forecasting changes in snow cover that will accompany forecasted changes in temperature and precipitation. The comprehensive, high spatial resolution snow cover duration datasets derived from Landsat can also be used to identify areas of persistent ice and snow cover, such as glaciers and perennial snow cover, as demonstrated in Chapters 3 and 4. While automated identification of areas of PISC for a single time period represents a significant step forward, the next step is to monitor changes in persistent ice and snow cover over time. While numerous efforts have already explored changes in glacier area for various regions and over various periods of time, the use of automated techniques that exploit the full Landsat data archive will allow this process to be standardized and extended to regions where previous analysis has not been conducted. The creation of regional 30 m resolution datasets at regular temporal intervals can provide further insight into the processes responsible for changes in glaciers by examining the distribution of changes in relation to topography. For example, glaciers with accumulation zones in protected cirque basins (ideal for both enhanced accumulation and reduced insolation) are often less responsive to regional climate signals (Hoffman et al., 2007). Landsat-derived PISC datasets can be used to test this and other hypotheses. Insights derived from this type of analysis can be applied to improve predictions of change for individual glaciers and snowfields. 124 Landsat-derived snow cover duration datasets can also be used to assess the impact of changes in snow cover duration on the distribution of plant communities. Research has demonstrated feedback loops between snow cover and vegetation often play a role in the establishment of shrubs in arctic tundra (Sturm et al., 2001; Sturm et al., 2005) and trees in alpine tundra (Bekker, 2005; Moir et al., 1999). Since these types of land cover conversion typically occur incrementally and begin as changes isolated to small patches < 100 x 100 m in size, changes occurring since the establishment of satellite remote sensing programs are best observed at finer spatial resolutions such as the 30 m resolution of Landsat. Finally, the detailed patterns of snow cover duration available from Landsat provide the opportunity to assess the impact of snow cover patterns on animal movement, habitat preferences, and reproductive success. For example, ungulates often search out landscape patches with the shallowest snow cover (Ball et al., 2001). Snow cover duration combined with physically-based snow modeling can be used to reconstruct depth and snow water equivalent (Molotch, 2009), over the course of a winter and identify these areas. In another example, caribou often seek out late lying snow patches (easily identified from a Landsat-derived snow cover duration product) for protection from mosquitos during the spring calving season. The individual questions and scientific applications addressed here represent only a limited subset of those that can be explored using Landsat snow cover data. In summary, the high spatial resolution, wall-to-wall coverage, and relatively long period of record of the Landsat sensors have the potential to provide insight into changing snow and ice cover conditions across arctic, alpine, and 125 montane systems that would not be possible using limited in situ observations or coarser resolution remote sensing data. 6.1 References Ball, J.P., Nordengren, C., & Wallin, K. (2001). Partial migration by large ungulates: Characteristics of seasonal moose Alces alces ranges in northern Sweden. Wildlife Biology, 7(1), 39-47. Bekker, M.F. (2005). Positive feedback between tree establishment and patterns of subalpine forest advancement, Glacier National Park, Montana, USA. Arctic, Antarctic, and Alpine Research, 37(1), 97-107. Essery, R., & Pomeroy, J. (2004). Vegetation and topographic control of wind-blown snow distributions in distributed and aggregated simulations for an Arctic tundra basin. Journal of Hydrometeorology, 5(5), 735-744. Hoffman, M.J., Fountain, A.G., & Achuff, J. M. (2007). 20th-century variations in area of cirque glaciers and glacierets, Rocky Mountain National Park, Rocky Mountains, Colorado, USA. Annals of Glaciology, 46(1), 349-354. Liston, G.E., & Sturm, M. (1998). A snow-transport model for complex terrain. Journal of Glaciology, 44(148), 498-516. Lundquist, J. D., & Flint, A.L. (2006). Onset of snowmelt and streamflow in 2004 in the western United States: How shading may affect spring streamflow timing in a warmer world. Journal of Hydrometeorology, 7(6), 1199-1217. Molotch, N.P. (2009). Reconstructing snow water equivalent in the Rio Grande headwaters using remotely sensed snow cover data and a spatially distributed snowmelt model. Hydrological Processes, 23(7), 1076-1089. Moir, W.H., Rochelle, S.G., & Schoettle, A.W. (1999). Microscale patterns of tree establishment near upper treeline, Snowy Range, Wyoming, USA. Arctic, Antarctic, and Alpine Research, 379-388. Sturm, M., Holmgren, J., McFadden, J.P., Liston, G.E., Chapin III, F.S., & Racine, C.H. (2001). Snow-shrub interactions in Arctic tundra: A hypothesis with climatic implications. Journal of Climate, 14(3), 336-344. 126 Sturm, M., Schimel, J., Michaelson, G., Welker, J.M., Oberbauer, S.F., Liston, G.E., ... & Romanovsky, V. E. (2005). Winter biological processes could help convert arctic tundra to shrubland. Bioscience, 55(1), 17-26.