Multispectral Classification and Reflectance of Glaciers: In situ data collection, satellite data algorithm development, and application in Iceland & Svalbard Allen J. Pope Trinity College & Scott Polar Research Institute University Of Cambridge This dissertation is submitted for the degree of Doctor of Philosophy 2013 ii iii Declaration This dissertation is the result of my own work and includes nothing which is the outcome of work done in collaboration except where specifically indicated in the text and Acknowledgements. It has not been submitted in whole or in part for a degree at any other university. This thesis is 146 pages long, including front matter, text, figure, tables and references. Allen Pope 28 July, 2013 iv vAcknowledgements My doctoral studies were supported by a graduate studentship from Trinity College, Cambridge as well as by the National Science Foundation Graduate Research Fellowship Programme under Grant No. DGE-1038596. Further research support came from UK Natural Environment Research Council’s Field Spectroscopy Facility, ARCFAC (the European Centre for Arctic Environmental Research), Trinity College Cambridge, Sigma Xi, the Norwegian Marshall Fund, the Explorers Club, the National Geographic Society Young Explorers Program, the Scott Polar Research Institute, the Cambridge University Geography Department, the Cambridge University Department of Anglo-Saxon, Norse, and Celtic Studies, and the Cambridge University Worts Fund. Thanks to Maria Pearman for her help sorting out all the finances. There are too many indviduals to thank for their help over the past four years, but I will try to do my best. Thank you to my primary advisor, Gareth Rees, for being a great guiding influence with a flexible but steady hand, and to my unofficial secondary supervisor, Ian Willis, for keeping me on the straight and narrow. For providing data, thank you to Finnur Pálsson, Neil Arnold, Florian Tolle, Alasdair Mac Arthur, and the Norwegian Meteorological Institute (via yr.no). For their thoughful comments as part of this thesis was submitted for publication, thank you to Marie Dumont, Jack Kohler, and three anonymous reviewers. For their help in talking through problems and finding solutions, thanks to my colleagues at SPRI and BAS, in particular: Toby, Steve, Liz, Kelly, Ruth, Sophie, Katya, Will, Georgina, Renuka, Willow, Alison, Nanna, Christine, Kat, Kate, Nick, Joe, Lauren, Jackie, Craig, Martin, Jon, Zuza, and Anna-Maria. Special thanks to Evan Miles for all his hard work on getting LandCor and 6S working. For feeding my fieldwork and outdoor learning addictions, thanks to the people of Tarfala Research Station, the Juneau Icefield Research Program, the University Centre in Svalbard, and the McCarthy International Summer School in Glaciology. For more fully developing me as a graduate student, thanks to Jenny Baeseman and all my colleagues at the Association of Polar Early Career Scientists and the UK Polar Network - you know who you are! Thanks to my Tweeps for discussion and distraction (#IcePhoto, #FieldPhotoFriday). For keeping me (mostly) sane, thanks to Alaine, Ali, Ali, Amy, Chris, Ellie, Enya, James, Jij, Katie, Meghan, Moritz, Nelly, Reed, Rob, Ryan, Sam, and Seb. Thanks to my brothers for keeping my studies in perspective, to my Dad for his formative guidance, and to my Mom for relinquishing her son to cold and snowy places but always being happy to talk about them. Last, but not least, I would like to acknowledge those that made the fieldwork possible for this thesis, including Fiona Danks and the staff of Norsk Polar’s Sverdrup Station in Ny-Ålesund, Addi Hermannsson and Ice, Ltd. for logistical support on Langjökull, and my field assistants Alexandra Giese, Richard Morris, Jon Peatman, and Fjord the Reindeer. vi vii Abstract Glaciers and ice caps (GIC) are central parts of the hydrological cycle, are key to understanding regional and global climate change, and are important contributors to global sea level rise, regional water resources and local biodiversity. Multispectral (visible and near-infrared) remote sensing has been used for studying GIC and their changing characteristics for several decades. Glacier surfaces can be classified into a range of facies, or zones, which can be used as proxies for annual mass balance and also play a significant role in understanding glacier energy balance. However, multispectral sensors were not designed explicitly for snow and ice observation, so it is not self-evident that they should be optimal for remote sensing of glaciers. There are no universal techniques for glacier surface classification which have been optimized with in situ reflectance spectra. Therefore, the roles that the various spectral, spatial, and radiometric properties of each sensor play in the success and output of resulting classifications remain largely unknown. Therefore, this study approaches the problem from an inverse perspective. Starting with in situ reflectance spectra from the full range of surfaces measured on two glaciers at the end of the melt season in order to capture the largest range of facies (Midtre Lovénbreen, Svalbard & Langjökull, Iceland), optimal wavelengths for glacier facies identification are investigated with principal component analysis. Two linear combinations are produced which capture the vast majority of variance in the data; the first highlights broadband albedo while the second emphasizes the difference in reflectance between blue and near-infrared wavelengths for glacier surface classification. The results confirm previous work which limited distinction to snow, slush, and ice facies. Based on these in situ data, a simple, and more importantly completely transferrable, classification scheme for glacier surfaces is presented for a range of satellite multispectral sensors. Again starting with in situ data, application of relative response functions, scaling factors, and calibration coefficients shows that almost all simulated multispectral sensors (at certain gain settings) are qualified to classify glacier accumulation and ablation areas but confuse classification of partly ash- covered glacier surfaces. In order to consider the spatial as well as the spectral properties of multispectral sensors, airborne data are spatially degraded to emulate satellite imagery; while medium-resolution sensors (~20-60 m) successfully reproduce high-resolution (2 m) observations, low-resolution sensors (i.e. 250 m+) are unable to do so. These results give confidence in results from current sensors such as ASTER and Landsat ETM+ as well as ESA’s upcoming Sentinel-2 and NASA’s recently launched LDCM. In addition, images from the Landsat data archive are used to classify glacier facies and calculate the albedo of glaciers on the Brøgger Peninsula, Svalbard. The time series is used to observe seasonal and interannual trends and investigate the role of melt-albedo feedback in thinning of Svalbard glaciers. The dissertation concludes with recommendations for glacier surface classification over a range of current and future multispectral sensors. Application of the classification schemes suggested should help to improve the understanding of recent and continuing change to GIC around the world. viii ix Contents Declaration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .ix List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii Acronyms and Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Rationale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Research Aims and Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Thesis Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2. Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review . . . . . . . . . . . . . . . . . . . . . . 5 2.1 Glacier Facies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 Linking Glacier Facies and Mass Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Glacier Facies and Energy Balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Snow, Ice, and Visible to Shortwave Infrared Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.1 Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.2 Modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4.3 Spectral Response and Facies Identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 Multispectral Remote Sensing of Glaciers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5.1 Choosing a Sensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.6 Classification and Glacier Facies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.6.1 Classification Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.6.2 Impact of Spatial Resolution on Multispectral Classification . . . . . . . . . . . . . . . . . . . 17 2.6.3 Multispectral Classification of Glacier Facies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6.4 Alternatives to Multispectral Imagery for Classification . . . . . . . . . . . . . . . . . . . . . . . 19 2.7 Albedo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.7.1 Estimation of Albedo from Landsat Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3. Field Sites and Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.1 Langökull, Iceland . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 Brøgger Peninsula, Svalbard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3 Landsat Imagery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.1 Obtaining Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3.2 Pre-Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 x3.4 Airborne Thematic Mapper Campaigns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.5 In Situ Spectral Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5.1 Midtre Lovénbreen data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.5.2 Langjökull Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.6 In Situ Facies Classification with Time-lapse Imagery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4. Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1 Spectral Response Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.2 Multispectral Sensor Radiometric Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.3 Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.4 Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.5 Statistical Classification Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.6 Software Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5. Using In Situ Spectra to Explore Landsat Classification of Glacier Facies . . . . . . . . . . . . . . . . . . . . . . 49 5.1 Landsat ETM+ Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2 Qualitative Spectral Analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.3 Full-Spectrum Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.4 Spectrum Clustering Using Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 5.5 Building Transferrable Linear Combinations for Classification . . . . . . . . . . . . . . . . . . . . . . . . 56 5.6 Section Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6. Using Field Spectroscopy and ATM Imagery to Assess the Role of the Spatial, Spectral, and Radio- metric Properties of Multispectral Imagers in Glacier Surface Classification . . . . . . . . . . . . . . . . . . . . . . 61 6.1 Impact of Spectral and Radiometric Properties on Glacier Surface Classification . . . . . . . . . 62 6.1.1 Linear Combinations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 6.1.2 Midtre Lovénbreen Clustering and Classification Analysis . . . . . . . . . . . . . . . . . . . . . 65 6.1.3 Langjökull Clustering and Classification Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.1.4 Concluding Considerations on Spectral and Radiometric Properties . . . . . . . . . . . . 73 6.2 Impact of Spatial Resolution on Glacier Surface Classification . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.2.1 Experimental Strategy and Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.2.2 Independent Resolution Assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.2.3 Down-Sampled Resolution Assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.2.4 Up-Sampled Resolution Assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.2.5 Concluding Considerations on Spatial Resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 6.3 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 7. Glacier Albedo on the Brøgger Peninsula, Svalbard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7.1 Data Preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.1.1 Atmospheric Correction of Landsat Imagery. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.1.2 Albedo Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.1.3 Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7.2 Seasonal Variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 7.3 Interannual Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 xi 7.5 Section Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 8. Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 8.1 Synthesis of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 8.1.1 Using In Situ Spectra to Explore Landsat Classification of Glacier Facies . . . . . . . . . 95 8.1.2 Impact of Spatial and Spectral Resolution on Glacier Surface Classification . . . . . . 96 8.1.3 Glacier Albedo on the Brøgger Peninsula, Svalbard . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 8.2 Achievement of Research Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 8.3 Transferability of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 8.4 Potential for Glaciology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 8.4.1 Mass Balance Proxies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 8.5 Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 8.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 9. References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 xii xiii List of Figures Figure 2.1: Cross section of a glacier showing idealised glacier facies. After Williams et al. (1991). . . . . 5 Figure 2.2: (a) Where mass balance is more negative than the previous year, the equilibrium line is above the firn line. (b) Where mass balance is less negative than the previous year, snow completely covers firn from previous years’ accumulation. After De Ruyter De Wildt et al. (2002). . . . . . . . . . . . . . . 6 Figure 2.3: Snow and ice spectral reflectance from the ASTER Spectral Library 2.0 (Baldridge et al. 2009). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Figure 2.4: Modelled reflectance of snow as a function of grain size (0 – 1 mm) and wavelength (0.4 – 2.5 μm). From Nolin & Dozier (1993). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Figure 2.5: Comparison of the spectral bands of major multispectral sensors. Background spectrum is the same as ‘Fine Snow’ in Figure 2.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Figure 3.1: Regional maps and false-colour Landsat ETM+ images of (a) the Brøgger Peninsula’s Austre Brøggerbreen (AB), Midtre Lovénbreen (ML), and Austre Lovénbreen (AL) and (b) Langjökull, Ice- land. Landsat bands 4 (NIR), 3 (red), and 2 (green) are used for red, green, and blue respectively in the false-colour composite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 3.2: False-colour ATM image of Midtre Lovénbreen collected on 9 August, 2003. ATM bands 4, 3, and 2 are used for red, green, and blue respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 3.3: Measurement setup for collecting surface reflectance spectra as seen on Langjökull. Note the FieldSpec in the backpack, controlling laptop on a belly board, sensor attached to the end of a pole, and Spectralon white-reference panel on the tripod in the foreground. Angles indicate the relative position of the operator for multiple measurements of the same sample. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 3.4: In situ surface reflectance spectra from (a) Midtre Lovénbreen and (b) Langjökull. The thickness of each line in (a) includes error bars based on the combination of many spectra in each class; error bars are spaced every 100 nm for (b). Data in the ranges 1350-1460 nm and 1790-1960 nm have been removed due to confounding water absorption in these wavelengths. Horizontal bars in (a) are an example of weighted averages of the continuous data to mimic bands of Landsat 7’s ETM+ (see Section 4.1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 xiv Figure 3.5: Locations of field spectra collected on (a) Midtre Lovénbreen and (b) Langjökull. Contours for (a) were derived from a 2005 LiDAR digital elevation model and for (b) from a 2004 SPOT-derived digital elevation model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Figure 4.1: The relative spectral response functions for Bands 1-10 of the Airborne Thematic Mapper. Note Band 1 (heavy line), which has an anomalous second peak; all other bands have one well-defined peak. The grey bars are the bands as defined by the manufacturer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Figure 5.1: Unsupervised classification of Landsat ETM+ imagery of glaciers on (a) Svalbard’s Brøgger Peninsula and (c) Langjökull Ice Cap, Iceland and the reflectance of each class (b, and d, respectively). Basic geographical and spectral information can be used to recognize some classes, but identification of intermediate facies remains unreliable. While SWIR wavelengths are useful for defining glaciated and non-glaciated areas, absolute brightness and the slope of spectra in the VNIR region are more useful for identifying glacier facies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 Figure 5.2: Principal component analysis of in situ spectra from Midtre Lovénbreen and Langjökull. Note the similarity of PCs between field sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Figure 5.3: Clustering biplots of in situ spectra from Midtre Lovénbreen calculated using full-spectrum information: (a) PC1 vs. PC2 and (b) PC1 vs. PC3 as shown in Figure 5.2 and calculated using weighted averages of restricted spectra which mimic the bands of the ETM+, (c) LC1 vs. LC2 and (d) LC1 vs. LC3 as shown in Equations 5.7-5.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Figure 5.4: Clustering biplots of in situ spectra from Langjökull calculated using full-spectrum informa- tion: (a) PC1 vs. PC2 and (b) PC1 vs. PC3 as shown in Figure 5.2 and calculated using weighted aver- ages of restricted spectra which mimic the bands of the ETM+, (c) LC1 vs. LC2 and (d) LC1 vs. LC3 as shown in Equations 5.7-5.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Figure 6.1: Biplots of the first and second linear combinations of in situ spectral data from Midtre Lovénbreen modified to mimic the spectral and radiometric properties of a range of multispectral im- agers. See Table 6.1 for the details of all individual plots, (a) through (p). The ellipses in Figure 6.1a also show the three groups into which the data were classified. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Figure 6.2: Biplots of the first and second linear combinations of in situ spectral data from Langjökull modified to mimic the spectral and radiometric properties of a range of multispectral imagers. See Table 6.2 for the details of all individual plots, (a) through (p). The ellipse in Figure 6.2a highlights the clean ice spectra which were classified as distinct from the rest. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Figure 6.3: Midtre Lovénbreen surface classification using ATM imagery at 2 m resolution (a) and the same image degraded to 20 m (b), 30 m (c), 60 m (d), and 250 m pixels (e). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 xv Figure 7.1: Surface classifications of Austre Lovénbreen using oblique photography (a-c) and Landsat imagery (d-f) on 20 August, 2010 (a, d), 23 June, 2011 (b, e), and 17 August, 2011 (c, f). Stripes in the Landsat images are the result of the Landsat 7 SLC error. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Figure 7.2: Estimated Transient Accumulation Area Ratio (TAAR) for Austre Lovénbreen, Midtre Lovénbreen, and Austre Brøggerbreen for years spanning 1976 to 2011 calculated using an unsuper- vised classification of non-shadowed areas of Landsat images. The strong seasonal transition from snow-covered to ice-covered glacier is visible across all years, as well as onset of autumn in 1988. Error bars on TAAR estimates are on the order of 5-10% (see Section 8.4.1). . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Figure 7.3: Albedo estimates for the study area displayed seasonally, broken up by glacier as well as fa- cies. Different years are indicated by colour while different glaciers are indicated by symbol. (a) and (b) show albedo for entire glaciers, (c) and (d) show albedo for accumulation facies, and (e) and (f) show albedo for ablation facies. (a), (c), and (e) present averages of the data in (b), (d), and (f). As discussed above, from earlier studies, all measurements should be considered to have error bars on the scale of ~0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Figure 7.4: Long-term albedo estimates from Landsat imagery for three study glaciers in Svalbard; dif- ferent glaciers are indicated by the symbols in the legend in (a). Measurements are averaged over the whole glacier (a), accumulation facies (b), or ablation facies (c). As discussed above, from earlier stud- ies, all measurements should be considered to have error bars on the scale of ~0.1. . . . . . . . . . . . . . . . . 90 Figure 7.5: Long-term albedo estimates from Landsat imagery for three study glaciers in Svalbard. The data from Figure 7.4a have been corrected for seasonal effects using the 3rd order polynomial fit in Figure 7.3a. All measurements should be considered to have error bars of at least 0.1. . . . . . . . . . . . . . 91 xvi xvii List of Tables Table 2.1: Properties of popular and prominent multispectral sensors and their VIR bands. . . . . . . . . 13 Table 3.1: Landsat images used for surface classification and albedo timeseries of glaciers on the Brøg- ger Peninsula, Svalbard (see Chapter 7 for use). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Table 3.2: Descriptions of the 7 surface classes measured on Midtre Lovénbreen and 16 surface classes measured on Langjökull. See Figure 3.4 for BHR (reflectance) spectra. The numbered classes on Langjökull are the result of classes that were observed to be similar in the field but had sufficiently dis- tinct spectra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Table 3.2 (cont.): Images of the surface classes defined in the earlier table, presented in the same order. Photos show people, tripods, shoes, and a fieldwork mascot to provide a sense of scale. . . . . . . . . . . . . 32 Table 4.1: Radiometric calibration values for a range of multispectral imagers. QMIN and QMAX are the minimum and maximum DNs which can be recorded by a sensor (i.e. radiometric resolution). LMIN and LMAX are the minimum and maximum radiance measureable by a sensor (i.e. radiometric range) and are reported in units of W m-2 sr-1 μm-1. ESUN is the mean solar exoatmospheric solar irradiance and is recorded in units of W m-2 μm-1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Table 4.2: Descriptive terms for classification based upon Cohen’s K values, following Monserud and Leemans (1992). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 Table 5.1: Contingency table of in situ spectra groupings as observed in the field (rows) from Midtre Lovénbreen as compared to (a) clustered classes of a merged k-means clustering of PC1 and PC2, full spectra and (b) clustered classes of a merged k-means clustering of LC1 and LC2, Landsat-band restricted spectra (columns). Values are percentage of total spectra measured during the Svalbard field campaign. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Table 6.1: For Midtre Lovénbreen spectra: sensors, radiometric resolution, gain settings, K (Cohen’s Kappa) of clustering success, description of classification success (Monserud & Leemans 1992), and corresponding panel in Figure 6.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 Table 6.2: For Langjökull spectra: sensors, radiometric resolution, gain settings, K (Cohen’s Kappa) of clustering success, description of classification success (Monserud & Leemans 1992), and correspond- ing panel in Figure 6.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 xviii Table 6.3: Comparison of the results of the classifications presented in Figure 6.3 by considering relative area of each class as well as their aggregation into glaciologically meaningful groups of accumulation and ablation area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Table 6.4: Comparison of the results of the classifications presented in Figure 6.3 by down-sampling results with majority filtering to lower resolutions. A is overall agreement and K is Cohen’s Kappa, both presented in Section 4.5; perfect agreement would lead to a value of 1 for both A and K. . . . . . . . . . . . 77 Table 6.5: Comparison of the results of the classifications presented in Figure 6.3 by downscaling results to higher resolutions. A is overall accuracy agreement and K is Cohen’s Kappa, both described in Sec- tion 4.5; perfect agreement would lead to a value of 1 for both A and K. . . . . . . . . . . . . . . . . . . . . . . . . . 78 Table 7.1: Mean and standard deviation of BRDF correction factors calculated using Equation 7.1 and parameters for the suite of Landsat images used in this study. The range of reflectances was determined by the range for which empirical constants were determined to be valid. . . . . . . . . . . . . . . . . . . . . . . . . . 84 Table 7.2: Contingency table comparing in situ photography-based and Landsat-based surface classifi- cations for Austre Lovénbreen on 20 August, 2010 (a), 23 June, 2011 (b), and 17 August, 2011 (c). A is the overall classification agreement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 xix Acronyms and Abbreviations AAR – Accumulation Area Ratio AB – Austre Brøggerbreen AL – Austre Lovénbreen ARSF – Airborne Research and Survey Facility a.s.l – above sea level ASTER – Advanced Spaceborne Thermal Emission and reflection Radiometer ATM – Airborne Thematic Mapper AVHRR – Advanced Very High Resolution Radiometer BHR – Bi-Hemispherical Reflectance BRDF – Bidirectional Reflectance Distribution Function DEM – Digital Elevation Model DHR – Directional-Hemispherical Reflectance DN – Digital Number DOY – Day Of Year ELA – Equilibrium Line Altitude ETM+ – Enhanced Thematic Mapper Plus FSF – Field Spectroscopy Facility FWHM – Full-Width Half-Maximum GIC – Glaciers and Ice Caps HDR – Hemispherical-Directional Reflectance InGaAs – Indium Gallium Arsinide ISODATA – Iterative Self-Organising DATa Analysis LC – Linear Combination LDCM – Landsat Data Continuity Mission (Landsat 8) ML – Midtre Lovénbreen MMU – Minimal Mapping Unit MODIS – MODerate resolution Imaging Spectroradiometer MSS – MultiSpectral Scanner NDSI – Normalised Difference Snow Index NERC – Natural Environment Research Council NIR – Near InfraRed (wavelengths) NTB – Narrow To Broadband OLI – Operational Land Imager PC – Principal Component PCA – Principal Component Analysis xx Variables A – overall classification agreement A* – random chance classification agreement AG – glacier surface area AAR – Accumulation Area Ratio AARL – Lower bound on AAR AARU – Upper bound on AAR α – albedo b – annual area-averaged mass balance c – surface area of accumulation facies crc – number of observations in a contingency matrix, row r and column c c4 & c5 – empirical constants d – Earth-Sun distance DNBAND – digital number of a field spectrum emulating a multispectral sensor’s band DOY – Day of Year EC – error of commission EO – error of omission ESUN – mean exoatmospheric solar irradiance f – proportionality factor between narrowband reflectance and narrowband albedo k – proportionality constant for mass balance proxies K – Cohen’s Kappa λ – wavelength L – radiance Lλ – spectral radiance LMIN – a sensor’s minimum measureable radiance LMAX – a sensor’s maximum measureable radiance PLDA – Probabilistic Linear Discriminant Analysis SAR – Synthetic Aperture Radar SLE – Sea Level Equivalent SWIR – ShortWave InfraRed (wavelengths) TAAR – Transient Accumulation Area Ratio TM – Thematic Mapper UK – United Kingdom VIR – Visible to InfraRed (wavelengths) VNIR – Visible and Near InfraRed (wavelengths) w.e. – Water Equivalent xxi m – empirical constant N – total number of observations (for classification agreement assessment) Q – digital number QMIN – minimum recordable digital number QMAX – maximum recordable digital number ρλ – at-sensor spectral reflectance ρMIN – minimum measurable at-sensor spectral reflectance ρMAX – maximum measurable at-sensor spectral reflectance r – reflectance r(λ) – spectral reflectance rSUB – reflectance for a particular SUBset of the electromagnetic spectrum (e.g. nb to indicate a narrow spectral range in general, or Landsat2 to indicate a Landsat TM or ETM+ Band 2) R(λ) – relative spectral response θz – solar zenith angle xxii 1Chapter 1: Introduction 1. Introduction 1.1 Rationale While significant emphasis has been placed on a comprehensive understanding of the continental ice sheets because of their potential to raise global sea level by over 70 m (Anderson 1999; IPCC 2007), in the present day 60% of sea level rise is traced to significantly smaller glaciers and ice caps (GIC) which respond much more quickly to shifting climate (Meier et al. 2007). Indeed, GIC around the world have an estimated volume ranging from 0.35 ± 0.07 m sea level equivalent (SLE; Grinsted 2013) to 0.43 ± 0.06 m SLE (Huss & Farinotti 2012) or 0.60 ± 0.07 m SLE (Radić & Hock 2010), and they are projected to contribute 0.1240 ± 0.037 m to global sea level by 2100 (Radić & Hock 2011). GIC are central parts of the world’s hydrological cycle and are key to understanding regional and global climate change (e.g. Oerlemans 1994). Glaciers contribute to local biodiversity (Jacobsen et al. 2012) and mediate the hydrology and flooding of some mountain systems (Dahlke et al. 2012). GIC provide crucial water resources for large populations of the world necessary for drinking, agriculture, industry, and power generation (Hopkinson & Demuth 2006; Björnsson & Pálsson 2008; Baraer et al. 2011; Barry 2011; Bolch et al. 2012). The changing state of GIC has contributed significant geomorphological hazards to mountain populations, particular in the Andes (e.g. Hubbard et al. 2005; Harrison et al. 2006) and Himalaya (e.g. Reynolds 1998; Bajracharya & Mool 2009; Benn et al. 2012). Through annual mass balance measurements, glacier volume change studies, and analysis of ice core records, GIC provide information about past and present local climate variability. The mass balance or mass budget is “the change in the mass of a glacier, or part of the glacier, over a stated span of time” (Cogley et al. 2011). It integrates information about a given year’s temperature, insolation, and precipitation into one metric which can then be tracked across many years. In the case of significant change, extrapolation of the mass balance to a common reference geometry can be made (Huss et al. 2012). While there is growing acknowledgement of incompletely understood drivers of glacier advance and retreat (Roe 2011), GIC continue to be crucial to understanding climate variability. Glacier surface properties are integral to the behaviour of GIC. The equilibrium line altitude (ELA) and accumulation area ratio (AAR, see Section 2.2; Cogley et al. 2011) can be used as proxies for glacier mass balance (Braithwaite 1984; Dyurgerov 1996). The division of GIC into accumulation and ablation areas is just the beginning of classification of glacier facies, or zones (see Chapter 2; Benson 1960; Williams et al. 1991). In addition, the glacier surface controls much of a glacier’s energy balance (Cuffey & Patterson 2010). Energy balance models both assimilate remotely sensed data about glacier surfaces to improve their results (Machguth et al. 2009; van Angelen et al. 2012), as well as validate their results (De Woul et al. 2006; Braun et al. 2007). Although fieldwork is crucial to glaciology, the mass balance of only a small and biased portion of GIC is monitored in situ (Zemp et al. 2009) due to the extensive logistical costs involved in maintaining a monitoring programme. Remote sensing tools provide ideal solutions, yielding data at various levels Chapter 1: Introduction 2 of spatial and spectral detail for otherwise inaccessible areas (Campbell 2006). Many different and complementary satellite and airborne techniques provide a breadth of information about glaciers, but multispectral imagery (see Section 2.5) is often the ideal tool for studying glacier surfaces (Pellikka & Rees 2010). Reflectance information over a range of wavelengths, good spatial resolution, frequent repeat imaging, an extensive image archive, and often cost-free data access are what make multispectral imagers ideal tools for glaciologists. However, multispectral sensors were not designed by glaciologists or with remote sensing of snow and ice in mind. Multispectral satellites like the original Landsat were (and continue to be) designed for a range of tasks including agricultural, oceanographic, and atmospheric monitoring (Markham & Helder 2012). Therefore, it is not self-evident that they should be optimal for imaging GIC. While many studies have measured and modelled the reflectance of snow and ice (see Section 2.4), and there are many good algorithms for (automatically) classifying different snow and ice facies (see Section 2.6), there are no universal techniques for glacier surface classification which have been optimized with in situ reflectance spectra. Thus, the roles that the various spectral, spatial, and radiometric properties of each sensor play in the success and output of resulting classifications remain unquantified. The efficacy of glacier surface classification has implications for many related questions within glaciology. Efficient surface classification can facilitate widespread mass balance measurements (e.g. Rabatel et al. 2005; Rabatel et al. 2008) and validation of glacier melt and hydrology models (De Woul et al. 2006; Braun et al. 2007). Glacier facies classification and albedo calculation of classified zones are also important in building improved surface energy balance models and therefore understanding the role of melt-albedo feedback on glacier surfaces. GIC are quickly changing, resulting in local and global impacts. Multispectral imagers provide a wealth of information about glacier surface behaviour, and there are many well-developed and established methods for classifying glacier surface types. However, these methods have not been built up from in situ data, and various multispectral sensors have different spatial, spectral, and radiometric properties. The combination of these factors leads to questions about the transferability of both methods and results between satellite imagers, with large implications for how best glaciologists should harness multispectral data. 1.2 Research Aims and Approach The introduction above frames the importance of and raises many interesting questions about multispectral classification of glacier surfaces. Specifically, it points towards a focused analysis of glacier surface classification using in situ spectral reflectance data. Multispectral imagers are powerful tools, and with the increasing availability of a range of high quality multispectral data including the recently launched Landsat 8 (known pre-operation as the Landsat Data Continuity Mission or LDCM, Irons et al. 2012) and upcoming Sentinel 2 (Drusch et al. 2012), it is increasingly crucial that they are fully understood. This thesis therefore takes an inverse perspective; it aims to start with in situ data to investigate the extent of information available from full-spectrum data and what that means for efficient and consistent application across multispectral sensors with different band capabilities and combinations. Specifically, this investigation is broken into three components: 3Chapter 1: Introduction • First, starting with in situ reflectance spectra from the full range of surfaces measured at the end of the melt season, in order to capture the largest range of facies, what are the optimal wavelengths for glacier facies identification and classification? Is it possible to distil this information simply and with repeatability? How do in situ measured classes compare with surface types measured in Landsat imagery (chosen for its ubiquity) and grouped with a proven classification scheme? Are there limitations to the new classification scheme depending on the glacier or ice cap? • Second, every step towards full automation and transferability is beneficial for many reasons, including the increasing amount of multispectral data as well as the importance of research repeatability. From the knowledge of what wavelengths and wavelength combinations are important for glacier classification, what are the best band combinations for a range of commonly used multispectral sensors? How do the spectral and radiometric properties of these sensors limit or enhance their performance in glacier classification? And because sensors are characterised by both spatial and spectral properties, how does the spatial resolution of these various sensors impact the resultant surface classification? What does this mean for glaciological applications? • Third, the classification scheme used above with in situ spectral data is used on satellite imagery to investigate a standing glaciological question. It has been suggested that enhanced thinning in upper regions of Svalbard glaciers may be related to melt-albedo feedback in the accumulation area (James et al. 2012), similar to what has been seen in Greenland (Box et al. 2012). Surface classification is a necessary first step in this process to separate facies which will have different albedo signatures and changes. Using the Landsat data archive, how much change in the interseasonal and interannual albedo of the accumulation and ablation areas can be seen on the chosen study area over the past 35+ years? What implications does this have for understanding cryospheric change in Svalbard and the Arctic? By grouping the investigations into these three sections, the thesis considers the in situ data and surface classes on their own, then fully quantifies the implications for glaciological applications of a wide range of multispectral imagers, and finishes with an application of multispectral classification of glacier surface types. 1.3 Thesis Structure This thesis is broken into nine chapters. Following this introduction which frames and poses the questions that this thesis aims to answer, the reader is brought through a review of the relevant literature related to glacier facies, multispectral remote sensing, glaciers surface classification, and in situ and remote reflectance measurements of snow and ice in Chapter 2. Chapter 3 introduces the data and study areas, and Chapter 4 discussed the methods used to analyze the data. Chapters 5, 6, and 7 are dedicated to the three groups of questions described above before being combined in a final discussion and conclusions which includes considerations for future research. 4 5Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 2. Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 2.1 Glacier Facies Glacier surfaces exhibit a range of zones – wet and dry, snowy and icy, clean and dirty. In order to understand better the changing conditions of a glacier’s surface, the area can be considered to be divided into a set of systematic, idealised facies that are characterised by a particular set of properties relating to the metamorphosis of the snow or ice surface; facies range from dry snow at colder, higher elevations through to melting ice near the glacier terminus (Benson 1960; Williams et al. 1991). Consider a glacier at the end of its melt season such that all possible facies are expressed (see Figure 2.1). At high elevations on the glacier is the “dry snow” facies – this is snow that has been deposited in the past year but is not melting. Descending the glacier, temperatures rise above the melting point, and the dry snow line is reached. Next follows the “percolation” facies where snowmelt permeates only part of that year’s snow accumulation. The “wet snow” facies is defined as the area where surface snowmelt permeates the entire snowpack. Some also consider a “slush zone” where the water content is so high that it overcomes the integrity of the snowpack. The lower limit of snow accumulation from a given year is known as the snow line. On some glaciers, water may percolate through the snow and seep out below the snow line, subsequently refreezing into superimposed ice (Brown et al. 1999). This is a means other than ice flow by which mass is transferred downglacier. Beneath the superimposed ice zone is glacier ice which began as accumulated snow, has subsequently been compressed, and has flowed downhill through the glacier. One facies not pictured in Figure 2.1 is firn: snow which has survived at least one ablation season (Cogley et al. 2011). Whether firn is exposed as a facies depends on the mass balance of a year relative to other years (see Figure 2.2). Figure 2.1: Cross section of a glacier showing idealised glacier facies. After Williams et al. (1991). Dry Snow Facies Percolation Facies Wet Snow Facies Glacier End of melt season surface Glacier Bed Accumulation AreaAblation Area Ice Facies Snow Line Glacier Terminus Wet Snow Line Slush Zone Equilibrium Line Superimposed Ice Zone Dry Snow Line Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 6 ice rn snow Snow Line (Equilibrium Line) Firn Line a) ice rn snow Snow Line (Equilibrium Line) b) At the end of the winter in cold regions, the entire glacier surface can be described by the “dry snow” facies, but as the melt season progresses, the facies succeed upglacier to higher elevations. Not all glaciers will display all facies. Sometimes a firn zone will not be visible, temperate glaciers will not retain a dry snow zone (e.g. Sigurðsson 1994), and many Antarctic glaciers do not have a wet snow zone. Although there are a wide range of glacier facies, they can be combined such that the glacier is divided into two larger regions: the accumulation zone and the ablation zone; the transition between these two areas is the line of net zero annual mass change known as the equilibrium line (Cuffey & Patterson 2010; Cogley et al. 2011). Above the equilibrium line is the area of net annual accumulation (the accumulation zone) including dry snow, percolation, wet snow, slush, and superimposed ice facies. Because superimposed ice (accumulation) can be located downglacier of firn (ablation), the equilibrium line would be thought of as abstract rather than physically observable in that instance. Below the equilibrium line the glacier surface is characterised by annual mass loss; this ablation area is made up of the firn zone and bare glacier ice. There are two other important facies transitions: the snow line and the firn line. First, the snow line (see Figure 2.1) is a transient feature which moves upglacier as melt progresses in the summer. At any given time it is a good indicator (although not a perfect representation) of the glacier’s equilibrium line (Dyurgerov 1996; Cuffey & Patterson 2010; Cogley et al. 2011). Second, when the snow line has progressed sufficiently upglacier, the transition between bare ice and firn is known as the firn line (see Figure 2.2). As has been described, each different configuration of facies is evidence of a different metamorphic history. Based upon this extensive terminology, it is an understatement to write that facies distributions vary across glaciers both within seasons and across years, and not all facies are necessarily present on all glaciers. In addition, beyond these zones which are considered ‘facies,’ further surface classes can be Figure 2.2: (a) Where mass balance is more negative than the previous year, the equilibrium line is above the firn line. (b) Where mass balance is less negative than the previous year, snow completely covers firn from previous years’ accumulation. After De Ruyter De Wildt et al. (2002). 7Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review identified in situ and remotely. For example, there are extensive areas of wind glaze and sastrugi in Antarctica (Orheim & Lucchitta 1987; Kuchiki et al. 2011; Scambos et al. 2012). The presence of snow algae imparts a reddish tinge to an evolving wet-snow facies (Takeuchi 2009), and dust will darken the snow surface as well (e.g. Painter 2011). Debris cover on the glacier can be considered another type of surface class (e.g. Shukla et al. 2009; Casey et al. 2012), as can volcanic ash deposited on glacier surfaces from a nearby eruption. In this thesis, ‘facies’ are considered to be the idealised zones of glacier surfaces which relate directly to accumulation and melt, while ‘surface classes’ are the zones which can be distinguished from the surface. Classification of glaciers surfaces is discussed in Section 2.6. Based on in situ reflectance measurements (see Section 3.5), this study will select a range of glacier surface classes identified in the field and investigate their identification in multispectral satellite imagery. 2.2 Linking Glacier Facies and Mass Balance As mentioned above, some glacier facies are associated with accumulation while others are indicative of ablation. Identification of accumulation versus ablation classes can therefore be used as a proxy for a glacier’s mass balance, often in combination with further data such as a digital elevation model (e.g. Braithwaite & Müller 1978; Dyurgerov et al. 2009; Rabatel et al. 2005). The equilibrium line is simplified spatially into the equilibrium line altitude (ELA) – the average altitude of net zero annual mass balance (e.g. Braithwaite 1984; M. König et al. 2004; Chinn et al. 2005; Shea et al. 2012). While easiest to think of the ELA as a physical concept, it can also be used in the abstract, for example on tropical glaciers where seasonal dynamics mean that annual maxima of both accumulation and ablation occur at the same time (Benn et al. 2005). Another way of simplifying the spatial considerations of the accumulation area and ablation area is the accumulation area ratio (AAR). This number is calculated by measuring the surface area of a glacier’s accumulation area and dividing by the glacier’s total surface area (e.g. Kulkarni 1992; Hock et al. 2007). The ELA has been shown to have a linear relationship with annual mass balance: ELA = ELA0 + kb (2.1) where k is a constant, b is the annual area-averaged mass balance, and ELA0 is the ELA when b = 0 (Braithwaite & Müller 1978). Equally, AAR = AAR0 + kb (2.2) where k is another constant, b is the annual area-averaged mass balance, and AAR0 is the AAR when b = 0. While ELA0 is highly dependent on local elevation, temperature, and precipitation conditions, a normal value for AAR0 would be ~0.67 (Braithwaite & Müller 1978). It is important to note that these relationships are only valid for land-terminating glaciers and do not take into account mass loss via other processes such as calving. In both cases, k is an empirical parameter. Usually, many years of in situ observations are required to calculate k for Equations 2.1 and 2.2, and most studies of this type require long-term mass balance records (e.g. Braithwaite 2002; Cogley et al. 2011). However, by monitoring the transient snow line and mass balance at different points during the melt season, a reasonably robust relationship between AAR and mass balance can be determined in just one melt season (Dyurgerov 1996; Hock et al. 2007). ELA and AAR provide a more immediate and nuanced view of glacier-climate interactions than Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 8 simply looking at the (admittedly easier to monitor) terminus of a glacier. Flow dynamics, bed properties, and complex time lags strongly influence the terminus. For example, looking at the AAR of retreating glaciers has the potential to predict whether a given glacier is likely to reach a new, smaller, steady state configuration or should be expected to disappear completely (Pelto 2010). Similarly, the ELA can be seen as a sensitive climate indicator which is calculable not only from observations on glaciers but also from regional temperature and precipitation data (Ohmura et al. 1992; Lie et al. 2003); this enhances the regional and prognostic capabilities of the ELA. Success of using glacier facies as a proxy for measuring glacier mass balance varies widely depending on location and year. Due to similarities in optical properties (Brown et al. 2005), persistent difficulty is had in identifying the equilibrium line / snow line and accumulation area rather than the more distinct glacier ice-firn line (e.g. Hall et al. 2000; Langley et al. 2008). Some argue that identifying the firn line altitude is actually a more meaningful parameter as it provides a temporally-averaged climate signal (Brown et al. 1999). However, it is the snow line rather than the firn line which provides a sensitive indicator of how a given year has influenced a glacier. Beyond ELA and AAR, glacier facies can be related to mass balance in other ways. Snow is bright (in much of the visible and near-infrared) and ice is darker, therefore as the melt season progresses the glacier as a whole gets darker overall - specifically in proportion to the relative contributions of different glacier facies. This effect is also enhanced by debris cover which is exposed in the ablation area as the seasons snow melts away. In this way, it is possible to monitor glacier albedo as a tool for monitoring glacier mass balance (Greuell & Oerlemans 2005a; Greuell et al. 2007; Dumont et al. 2012). 2.3 Glacier Facies and Energy Balance Shortwave radiation is crucial to the energy balance of a glacier. Glacier facies meaningfully contribute to this radiation balance and therefore to the surface energy balance of GIC. The energy balance of a glacier is the sum of all energy inputs and outputs, including shortwave radiation, longwave radiation, sensible heat flux, latent heat flux, geothermal heat flux, and rain heat flux (Hock 2005). Different glacier facies contribute in very different ways to the components of a glacier’s energy balance, and therefore, increased knowledge about glacier facies will lead to increased understanding of glacier surface processes and change (e.g. Wolken et al. 2009). A clear example of the interrelated nature of energy balance and glacier facies can be seen in the simple parameterization of the degree-day melt model where ice and snow have different degree day factors (e.g. Hock 2003). The most relevant component of a glacier’s surface energy balance to multispectral classification is the shortwave radiation balance because both deal with the same part of the electromagnetic spectrum. Albedo is defined as “the ratio of the reflected flux density to the incident flux density, usually referring either to the entire spectrum of solar radiation (broadband albedo) or just to the visible part of the spectrum” (see Section 2.7; Cogley et al. 2011). As will be seen clearly in later sections, the spectral reflectance of different glacier surfaces classes varies widely; because reflectance at individual wavelengths varies, so too will the broadband albedo. Information about the interannual and intra-annual evolution of glacier surfaces is also a key parameter in building energy balance models. Fuller consideration of glacier facies in glacier 9Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A reviewmelt modelling is gaining increasing traction within the glaciological community (Machguth et al. 2009; Dumont et al. 2010). This is true not just for GIC, but also for the larger ice sheets, where better classification and description of the unique properties of different facies improve melt model behaviour (e.g. van Angelen et al. 2012). 2.4 Snow, Ice, and Visible to Shortwave Infrared Radiation Remote sensing harnesses all wavelengths from the electromagnetic spectrum to gain information, but this thesis is primarily concerned with multispectral imagery, and particular reflected visible to infrared (VIR) wavelengths. Here, the VIR spectrum is considered to be limited to 350-2500 nm. Subsets include the visible (~390-700 nm), near infrared (NIR, ~700-1400 nm), visible-near infrared (VNIR, ~350-1400 nm), and the shortwave infrared (SWIR, ~1400-2500 nm). While snow is often considered to be a homogeneous substance with uniform properties, it is in fact a highly variable, evolving medium. From the moment it is deposited, snow crystals change shape, configuration, size, and water content, all of which in turn influence the reflective properties of the snow surface. In addition, snow surface spectral reflectance can be influenced by underlying snow, ice, and ground surfaces up to 20 cm water equivalent (w.e.) beneath the surface (Wiscombe & Warren 1980). Both single and multiple scattering contribute strongly to the spectral reflectance of a snowpack. Therefore, determining snow properties from spectral reflectance is quite complex (Gardner & Sharp 2010). The first consideration of snow optics is to understand the relationship of the incoming light, the snow surface, and the measurement instrument. At low solar elevation angles, snow is a specular reflector, while at angles above 10° fresh snow can be considered a diffuse reflector (Choudhury & Chang 1981). This simplification is valuable for some applications, but as the snow surface ages it no longer reflects equally in all directions. The bidirectional reflectance distribution function (BRDF; see section 2.7) describes reflectance as a function of both incoming light direction and view angle in order to interpret remote measurements for such anisotropic surfaces (Warren 1982). 2.4.1 Observations There is a strong history of measuring the spectral reflectance of snow and ice surfaces to build from, dating back to Hulburt (1928) and O’Brien and Munis (1975) who based their studies only on seasonal snow. Considering glaciers, Zeng et al. (1983) profiled some facies on a glacier in western China. Further studies measuring snow and ice spectral reflectance have been conducted on the summertime accumulation area of glaciers in Svalbard (Winther 1993), snow covered tundra in Svalbard (Gerland et al. 1999), Hintereisferner in Austria (Hendriks & Pellikka 2004), snow-covered and snow-free terrain in Greenland (Rott & Søgaard 1987), snow and ice cover in Antarctica (Grenfell et al. 1994; Casacchia et al. 2002), ablation facies on Peyto glacier (Cutler & Munro 1996), debris-covered glacier lithologies in Nepal (Casey et al. 2012), snow algae and seasonal glacier change (Takeuchi 2009), and to measure snow impurities (Kuchiki et al. 2009; Painter 2011). However, with respect to potential applications of surface classification, it is important to note that none of the available data sets spans the full range of the Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 10 VIR spectrum, at high spectral resolution, with a complete range of surface types present on an Arctic/ temperature glacier at the end of the melt season. Hendriks and Pellikka (2004) do observe a suite of glacier surface types but qualify these observations with the acknowledgement that wider sampling is necessary; this work builds on their study. A typical snow spectral response (see Figure 2.3) starts with very high reflectance in the visible, dropping off substantially towards the NIR. In other words, while snow is strongly reflective in the visible, it is almost “black” or strongly absorbing in parts of the NIR (Warren 1982). There is a peak near 1100 nm and a shoulder and quick drop-off around 1400 nm. There are triangular peaks at 1850 nm and 2250 nm; the former is asymmetrical while the latter is symmetrical. It is worth noting that the samples in Figure 2.3 are lab-prepared. Field spectra will be affected by processes such as impurity deposition, grain metamorphosis, and increased water content. This will be most notable in the SWIR and NIR, and the physical considerations will be addressed in the next section on modelling. In addition, Figure 2.3 shows water ice as having very low reflectance. Glacier ice will behave completely differently (see Section 3.5.4). Although the two are chemically identical, they are not so structurally; the same can be said of ice and snow, too. This highlights the importance of absorption versus scattering. Often highly reflective due to a high concentration of air bubbles, glacier ice has also been observed to be blue and even green in certain circumstances (Warren et al. 1993) Snow surface contamination and pollution also significantly influence spectral response; this makes intuitive sense as blue or white ice and snow begin to look grey and brown with soot or dust deposition. Kuchiki et al (2009) used a ground-based spectral radiometer to measure the mass concentrations of snow impurities (updated by Painter 2011). Decreased albedo as a result of increasing contaminant concentration is thus another way that older, more metamorphosed snow surfaces can be distinguished from fresher ones.     $"   "&" %'#% '  !% ' Figure 2.3: Snow and ice spectral reflectance from the ASTER Spectral Library 2.0 (Baldridge et al. 2009). 11 Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review It is worth noting that all of these measurements were taken as small-scale observations. Remote sensors would average such information across an entire (relatively large) pixel thus providing a more homogenised viewpoint. The impact of resolution on classification will be more fully addressed in Chapter 6. 2.4.2 Modelling Mathematical modelling of reflection by and radiative transfer through snow and ice yields significantly enhanced process-based understanding of the factors which will influence the spectral response of glacier surfaces. Initial modelling studies (Dunkle & Bevans 1956) were limited by lack of reflectance spectra data to constrain their models, especially in the infrared. Wiscombe and Warren (1980), Warren and Wiscombe (1980) and others afterwards (e.g. Gardner & Sharp 2010) much more satisfactorily modelled the reflectance of snow based on variables such as grain size, density, liquid water content, and black carbon content. Recent studies have also combined observation and modelling to show the role of chemical contaminants (e.g. black carbon, trace elements) in snow spectral reflectance (Carmagnola et al. 2012; Kokhanovsky 2013). Indeed, impurities such as soot, dust, and, most relevant to this study, volcanic ash, lower snow reflectance especially at wavelengths less than 900 nm; this effect is strongest in coarse grained snow and ice (Warren & Wiscombe 1980; Gardner & Sharp 2010; Kokhanovsky 2013), confirmed with experimental measurements (Painter et al. 2007). For deep snow, grain size is the dominant factor in determining snow surface albedo (see Figure 2.4; Wiscombe & Warren 1980; Nolin & Dozier 1993; König et al. 2001). A snowpack with a large grain size is characterised by a low albedo; conversely, a snowpack with a fine grain size is characterised by a high albedo. As snow ages, grain size grows and reflectance decreases (Nakamura et al. 2001). This Figure 2.4: Modelled reflectance of snow as a function of grain size (0 – 1 mm) and wavelength (0.4 – 2.5 μm). From Nolin & Dozier (1993). Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 12 physical basis for changing reflectance will aid in facies classification. This effect can also be mimicked by increasing the liquid water content of a snowpack. The water replaces air in voids and increases the apparent grain size, thereby resulting in a decreased albedo (Wiscombe & Warren 1980; Hall et al. 2009). Consider light propagation through snow and ice statistically. Light is scattered at the air-ice interface and absorbed while passing through ice. Increasing the effective grain radius effectively increases the relative likelihood that a photon will be absorbed rather than reflected out of the snowpack (Gardner & Sharp 2010). This understanding has led to modelling of quantitative relationships between spectral reflectance and snow grain size. Dozier et al. (1981) showed that NIR reflectance could be used to roughly estimate snow grain size, and this technology has progressed to the state where Tape et al. (2010) used NIR photography to record microscale variations in snowpack stratigraphy. Further recent applications of snow grain size measurement are numerous (Nolin & Dozier 1993; Aoki et al. 2000; Kokhanovsky & Zege 2004; Painter et al. 2007; Jin et al. 2008; Montpetit et al. 2012). Additionally, as can be seen in Figure 2.4, the response to a variable such as grain size can be amplified or diminished depending on the observed wavelength. It is important to consider this behaviour when designing an experiment. For example, a sensor is more sensitive to grain size and less sensitive to atmospheric aerosol interference or snow water equivalent if the correct channel is chosen (Dozier et al. 1981). In addition, the quality of a signal at a sensor must be considered; images of snow covered areas are often saturated in the visible spectra while retaining quantitative information in the infrared (Hall et al. 1988). An interesting case study in the variability of snow reflectance and complementary roles of observation and modelling is provided by the SWIR reflectance peaks. While on the one hand captured by many models and observations (e.g. Wiscombe & Warren 1980; Casacchia et al. 2001), on the other they are not present in many measured spectra (e.g. Hall et al. 1988; Boresjö Bronge & Bronge 1999). It has been suggested that liquid water content and snow metamorphosis will decrease reflectance in these bands while the presence of a dry snowpack and angular crystals near the surface (e.g. hoar frost) will contribute to the peaks (Casacchia et al. 2001). In most satellite-observed reflectances used for glacier facies identification, these SWIR bands show snow as dark, even darker than the surrounding geology, and hence they can be used for glacier delineation (Paul 2000; Albert 2002). 2.4.3 Spectral Response and Facies Identification As the preceding sections show, snow reflectance is heavily wavelength-dependent. In particular, the NIR has been seen as containing quantitative information about snow and ice surfaces. Although specific wavelengths are often determined by the available sensor, for grain size measurements Nolin & Dozier (1993) used a band centred at 1.03 μm and Li et al. (2001) used radiances at 0.96, 1.05, 1.24, and 1.73 μm; many studies agree with this chosen spectrum (Kokhanovsky & Zege 2004; Kokhanovsky & Schreier 2009). Glacier facies classification, too, has focused on the NIR to the exclusion of the visible, although snow studies have highlighted both ranges (Zeng et al. 1983). While there will be a fuller treatment of classification strategies in Section 2.6, it is worth mentioning here that both Sidjack and Wheate 13 Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review (1999) and Braun et al. (2007) cited some saturation in the visible and enhanced performance in the NIR as reasons for choosing linear combinations of input multispectral bands which contained large contributions from the NIR and minimal contributions from the visible. Based on these examples, it is natural to hypothesise that sensors with enhanced capabilities in the NIR (e.g. Sentinel-2 and ATM, see Section 2.5) will be able to classify glacier facies better than their counterparts. This belief will be put to the test in Chapter 6. 2.5 Multispectral Remote Sensing of Glaciers Remote sensing can broadly be defined as collecting information from an object without making direct physical contact with it but the definition is generally more restricted to sensors carried on an airborne or spaceborne platform. More specifically, multispectral remote sensing is a type of passive observation which measures reflectance of electromagnetic radiation from the Earth’s surface Table 2.1: Properties of popular and prominent multispectral sensors and their VIR bands. Sensor / Revisit Time Airborne Thematic Mapper (ATM) As ordered ASTER 16 days Landsat ETM+ (Landsat 7) 16 days Landsat MSS (Landsat 1-5) 18 days (1-3), 16 days (Landsat 4-5) Landsat TM (Landsat 4-5) 16 days Landsat OLI (LDCM a.k.a Landsat 8) 16 days MODIS Sub-daily Sentinel-2 5 days Band 1-10 1-3 4-9 1-5, 7 8 4-7 1-5, 7 1-7, 9 8 1-2 3-7 8-19, 26 1, 9, 10 2-4, 8 5-7, 8a, 11, 12 Spatial Resolution Sub-metre to metre-scale, depending on airborne platform flying height 15 m 30 m 30 m 15 m 60 m 30 m 30 m 15 m 250 m 500 m 1000 m 60 m 10 m 20 m Information from NASA 2002; NASA 2011; Drusch et al. 2012; Helder et al. 2012a; Helder et al. 2012b; Irons et al. 2012; http://arsf.nerc.ac.uk/instruments/atm.asp; http://asterweb.jpl.nasa.gov/; http://modis.gsfc.nasa.gov/; http:// landsat.gsfc.nasa.gov/about/mss.html; and http://landsat.gsfc.nasa.gov/about/tm.html. Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 14 Fi gu re 2 .5 : C om pa ris on o f t he sp ec tr al b an ds o f m aj or m ul tis pe ct ra l s en so rs . B ac kg ro un d sp ec tr um is th e s am e a s ‘ Fi ne S no w ’ i n Fi gu re 2 .3 .        # "      "   !   !                                           15 Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review simultaneously in multiple different wavebands at once. A typical multispectral scanner will have four to ten bands recording images in the visible – near infrared, shortwave infrared, and thermal infrared (Campbell 2006). Multispectral remote sensing takes advantage of coarse variations in spectral response. For example, while snow and cloud both reflect strongly in the visible, they respond quite differently in the infrared (Hall et al. 1995). Multispectral remote sensing images are some of the most prevalent, easily available, and versatile forms of data available for the Earth Observing glaciologist. Table 2.1 and Figure 2.5 describe the important properties of a range of popular and prominent multispectral imagers that are considered in this thesis. Multispectral data is often quite intuitive to use, as true-colour composites are easily generated. In addition, multispectral sensors are calibrated digital radiometers which generate at-sensor radiance and can be used to calculate at surface reflectance (see Section 4.3). Despite the fact that any multispectral sensor is easily thwarted by cloud cover, these data remain widely used for studying GIC. Glaciological uses of multispectral imagery include glacier maps and inventories (Paul 2000; Albert 2002; Kargel et al. 2005; Paul & Kääb 2005; Hendriks & Pellikka 2008), albedo calculation (Greuell et al. 2002; Knap, Reijmer, et al. 1999), identification of surface and basal crevasses (Luckman et al. 2012), feature tracking (Heid & Kääb 2012), interpolating digital elevation models (Pope et al. 2013), identifying ice sheet grounding lines (Bindschadler et al. 2011), and much more (Rees 2006; Pellikka & Rees 2010). Analysis techniques applied to glacial multispectral remote sensing include band thresholding (investigating snow properties, Gupta et al. 2005), band ratio combinations to reduce the effect of shadows (Hall et al. 1992) or enhance the effect of grain size variations (Boresjö Bronge & Bronge 1999), and spectral mixture analysis studies (e.g. snowline detection, Klein & Isacks 1999; sub- pixel variability in snowcover, Dozier & Painter 2004). Normalised differencing, and specifically the Normalised Difference Snow Index (NDSI), deserves a special mention. The NDSI is the ratio of the difference and sum of reflectance in the green and the SWIR (usually Landsat bands 2 and 5). Originally intended for distinguishing snow and cloud (Hall et al. 1995), the NDSI has been used widely to image snow cover, glacier extent, and even some glacier facies (e.g. Albert 2002; Gupta et al. 2005; Cea et al. 2007) 2.5.1 Choosing a Sensor There are a variety of factors which must be weighed in choosing an appropriate multispectral sensor; each separate investigation or task requires an imager which is fit for purpose. Major considerations include: • spatial resolution • spectral resolution / band wavelengths • radiometric resolution and range • temporal resolution / revisit time • data cost and ease of access • length of data archive / data availability • availability of pre-processed products Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 16 From the range of different options, it is clear why no one sensor will be optimal for all studies. For glacier facies classification, this thesis aims to understand what the benefits and limitations are for a range of sensors. With this goal in mind, sensors were chosen considering current popularity in the community and ease of data access while maintaining a diversity of spatial, spectral, radiometric, and temporal resolutions (see Table 2.1 and Figure 2.5). There are many, many other multispectral sensors which could just as easily have been considered, but the vast majority will have some analogue among the sensors considered here. 2.6 Classification and Glacier Facies Image calibration and processing (see Section 3.3) is an important first step in deriving meaning from multispectral remote sensing. The next step is classification, the process which takes quantitative information from every pixel and places each into one of a group of discrete categories. It is at the classification stage that data are not only manipulated but interpreted. For example, while a pixel has particular brightness values in bands 2, 3, and 5, classification interprets these and places a pixel into category A, B, or C (Campbell 2006). This section includes information both on remote sensing image classification in general, as well as how it has been applied to glacier surfaces. While all classifications are assessed on their accuracy, in the past computational efficiency was paramount. Although computational ease is still a factor, automation and transferability between scenes and data sets are becoming more important. 2.6.1 Classification Techniques There are a wide variety of classification techniques, the broadest categories of which are supervised and unsupervised classification. Supervised classification takes into account user input when directing the types of classes that are then output. While unsupervised techniques are based on an independent algorithm which uses only the distribution of data to determine the classes that it will produce (Campbell 2006). Unsupervised classification does not require prior knowledge of the study area, there is little associated human error, and classes are by definition spectrally unique. However, the user has little control over the identity of each class and spectral classes are not necessarily intuitive or meaningful. Supervised classification, on the other hand, requires users to “train” the algorithm by selecting examples of each class to be defined. This means that while supervised classifications are tailored to a site or image, by definition have applicable meaning, and can accept a wide range of spectral variability, users must also be skilled in the process and remain wary of errors imposed by incorrect or unrepresentative training data. The simplest classification technique is setting image thresholds. For example, any surface above a certain radiance or reflectance is termed snow. This can be done on an absolute scale or based on the image’s histogram, including a single band or some algebraic manipulation of multiple bands. Image thresholds can therefore require a range of levels of user intervention. However, because of their simplicity, thresholds can inherently be easily confounded and require careful application. Unsupervised classifications have had significant success in classifying glacier facies not only 17 Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review because they are easily reproducible but also because they are often able to exploit subtle features within data sets. While some statistically-based techniques show promise for automated identification of landcover (e.g random forests, Gislason et al. 2006; or PLDA, Brenning 2009), ISODATA (Iterative Self- Organizing Data Analysis; e.g. Nolin & Payne 2007; Wolken et al. 2009) and k-means classification (e.g. König et al. 2004; Barcaza et al. 2009) are the most widely and easily implemented clustering algorithms for glacier surface classification. In both, candidate clusters are initially selected (with or without user input) and their means are allowed to migrate in the spectral domain, optimizing the classification with each iteration; each pixel is re-compared to each cluster mean and nearest pairs are matched. Means are recalculated and the process is repeated until the pixel-to-cluster-means distances are minimised to within a certain percentage change or a set number of iterations is reached. Supervised classifications, on the other hand, embrace user input by taking into account intuitive training areas rather than selection between computer-defined classes. As mentioned above, both the strength and weakness of supervised classification is that the user directs the classes that are selected; this allows more control but is only as reliable as the training areas which are provided. Hard classifiers then place pixels into the “closest” class where “closest” can be defined in a variety of ways (e.g. spectral angle mapping and parallelepiped grouping), whereas fuzzy classifiers include a degree of membership that a pixel represents to a particular class. All classifications presented so far have been pixel-based. However, as humans naturally process images, it is intuitive to think not just about individual pixels but small groups or regions. This method of reducing image heterogeneity is known as image segmentation or object-based image analysis and can often serve to make the classifications more effective and intuitive (e.g. Baker et al. 2013). Object- based image analysis software has been expensive and limited, and therefore of limited use to the wider community, but it is possible that this will change. In a sense, different sized pixels can also been seen as different sized (square) objects. On the one hand, very high resolution imagery and point clouds can have very high success because pixels can match end-members well (e.g. Höfle et al. 2007). Conversely, lower spatial resolution can be an advantage as it increases pixel homogeneity (e.g. Hall 2008). The balance of these two factors is given significantly further consideration in Chapter 6. 2.6.2 Impact of Spatial Resolution on Multispectral Classification Coarse imagery often confounds land cover classification as a result of pixels representing a mixture of classes. In this case, increasing resolution would increase classification success. Resolution effects have been investigated in a range of different studies, including those looking at urban heat islands (Sobrino et al. 2012), coral reefs (Phinn et al. 2012), wetland vegetation (Michishita et al. 2012), disaster relief management (Battersby et al. 2012), and the impacts of natural gas drilling on forest land (Baker et al. 2013). The results of these comparisons have shown that agreement between spatial scales varied depending on the type of surface and the scale of inhomogeneity. In other words, although it is not immediately intuitive, increasing spatial resolution can decrease classification success. This happens in the case where coarser imagery serves to smooth out spatial inhomogeneity within classes. The question then is, what effects come in to play for glacier surface classification? To the author’s knowledge, no study has directly considered this question (see Section 6.2). Previous studies have used Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 18 higher resolution imagery (~10 and ~1 m) to assess the accuracy of glacier extent measurements using medium resolution (~30 m) imagery (Paul 2000; Paul et al. 2013). There is no significant difference in measured glacier area using imagery at 60 m resolution and finer; lower resolution imagery was not tested (Paul et al. 2002). For snow cover, 250 m resolution imagery from MODIS provides much better results than 500 m from the same sensor (Notarnicola et al. 2013a; Notarnicola et al. 2013b), but daily temporal resolution is likely not needed for glaciology and glacier surface topography is less complex than most land topography. In glacier albedo calculations, the scale of albedo variations is smaller than 30 m TM pixels (Reijmer et al. 1999), but because albedo can vary within facies this does not necessarily mean glacier facies vary on the same scales. Albedo variations within facies would then confound glacier surface classifications with spatial resolution finer than 30 m. 2.6.3 Multispectral Classification of Glacier Facies Throughout the literature, the whole range of different classification techniques has been applied to glacier surfaces. It should be noted that (automated) classification of glacier extent is considered to be a separate problem, one which has been largely solved, with the exception of debris-covered areas (Paul et al. 2013). This study therefore focuses on classification of different surface types once a glacier outline has been provided, highlighting the more subtle differences between facies rather than the quite significant differences between glacier and surrounding terrain. Many different techniques have been applied to multispectral data to identify glacier surface classes. Starting with measuring glacier surface reflectance and segmenting into different zones (Hall et al. 1988), simple multispectral techniques include thresholding of individual bands and conventional indices such as the Normalised Difference Snow Index (Albert 2002; Cea et al. 2007; Gupta et al. 2005; Williams et al. 1991). Albert (2002) showed that band ratios can be used to help identify glacier boundaries in the Andes and, Paul et al. (2002) did likewise in the Alps. Keshri et al. (2009) used band ratios to classify glacier surfaces in the Himalaya. In addition, spectral mixing has been applied to identify spectral end-members on glacier surfaces (Klein & Isacks 1999). Applications of multispectral classifications to glaciology include snowline identification (McFadden et al. 2011), employing unsupervised clustering for the same purpose (Aniya et al. 1996; Barcaza et al. 2009), classifying snow zonation with a similar method (De Angelis et al. 2007), and classifying snow, slush, and ice (Williams et al. 1991). Supervised classifications of multispectral data have been applied to glacier surfaces, for example, by using the green, near-infrared, and short-wavelength infrared to identify snow, ice, mixed ice and debris, valley, rock, and water in the Himalaya (Shukla et al. 2009), and in order to identify two types of ice and three types of snow across the entirety of Hielo Patagónico Sur (De Angelis et al. 2007). While Aniya et al. (1996) used a parallelepiped grouping to group facies, Heiskanen et al. (2002) used a maximum likelihood classification scheme. Albert (2002) tested multiple grouping techniques for supervised classifications identifying snow, ice, water, moist vegetation, and dry ground and found for that purpose that spectral angle mapping and optimised angle input rules performed better than parallelepiped grouping. However, this conclusion was based on a large range of classes (rock to snow) and should be reconsidered for application to more subtle transitions (ice, firn, etc.). Successful unsupervised classifications have been used for comparison to surface melt and glacier 19 Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review hydrology models (Braun et al. 2007; de Woul et al. 2006); classes must be nuanced for this task. De Woul et al. (2006) classified glaciers identified by three different masks (Landsat band 4 thresholding, NDSI, and normalised difference vegetation index). This resulted in 20-30 classes which were then manually reviewed and recoded into four classes of snow, firn, ice, and non-glacier/no-data. Similarly, Braun et al. (2007) identified firn, snow, and ice with an ISODATA classification clustering based on NDSI, a NIR/SWIR band ratio, and three principal components (see Section 4.4) which maximised NIR and minimised the visible. After initial processing of the Landsat data, 10-20 classes were returned and were subsequently merged by visual inspection of an experienced operator. In addition, final classes were compared to mass balance models of glaciers to good agreement, although further in situ corroboration is desired. These examples provide workable methods of identifying glacier facies using ISODATA classification and multispectral imagery and also suggest viable ways forward. It is these techniques that employ an ISODATA clustering algorithm which are used as a point of departure for this thesis. Beyond these two studies, most germane to this investigation are therefore other studies which applied clustering algorithms to multispectral data. Aniya et al. (1996) pioneered the use of ISODATA classification for glacier facies, employing Landsat TM bands 1, 4, and 5 to identify bedrock, accumulation facies, and ablation facies. Barcaza et al. (2009) used the same technique as well as a k-means unsupervised classification and corroborating aerial surveys to identify zones of shadow, bare ice, non-ice, slush, and wet snow from normalised difference images; they subsequently estimate the regional ELA by using a separately-derived digital elevation model. Aniya et al. (1996), on the other hand, employed ISODATA results to identify suitable training areas for a supervised classification, and De Angelis et al. (2007) did similarly for the Southern Patagonian Icefield, using an unsupervised classification of band ratios to suggest ideal training areas. 2.6.4 Alternatives to Multispectral Imagery for Classification Glaciological and remote sensing research literature is demonstrably well represented with studies that have used multispectral imagery to look at glacier surfaces. Nevertheless, this approach has its limitations - multispectral imagery can only be used under clear, daytime conditions. Beyond multispectral data, radar imagery and lidar data have been used to investigate glacier facies and can be complementary techniques. Lidar can be thought of as being akin to multispectral classification because it also looks at a reflected signal. However, the sensor is active rather than passive and uses just one wavelength which is specified by the instrument’s ranging laser. The strength of the return of the laser pulse can be processed to be consistent across swathes and then threshold can be applied to yield different classes (e.g. Lutz et al. 2003). Höfle et al (2007) use more of an object-based approach, taking into account both local roughness and reflected signal to classify ice, firn, snow, and crevasses. Despite these successes, lidar data can be expensive to collect and somewhat complex to process, and campaigns must be planned for specific targets. Radar, on the other hand, measures backscatter rather than shortwave radiance. It therefore depends on roughness and bulk and dielectric properties down to a certain penetration depth (Pellikka & Rees 2010). Synthetic aperture radar has been used to image dry and wet snow facies on the Greenland Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 20 Ice Sheet (Joshi et al. 1998), identify the firn line on Icelandic glaciers (Brown et al. 1999), and classify glacier surfaces and the snow line in Tibet (Huang et al. 2011; Li et al. 2012). Similar to multispectral techniques, there are often issues in distinguishing the finer gradations of facies like superimposed ice and the transient snow line (Brown et al. 1999; Demuth & Pietroniro 1999). Because dry snow is effectively transparent to microwave radiation, radar can be used in dry snow areas to get a better idea of facies despite a winter snow cover (König et al. 2002); ground penetrating radar can also be used to validate and compliment SAR facies classifications (Langley et al. 2008). While advantageous because it can penetrate clouds and measure during the polar night, radar is largely confounded in the wet-snow conditions which persist on many GIC. In addition, multispectral imagery is generally available at higher resolution and wider distribution at lower or no cost. Therefore, while recognizing the value of alternative measurements of glacier surfaces, this study will continue to focus on multispectral data. 2.7 Albedo Observation techniques based on shortwave radiation have the potential to allow researchers to obtain large amounts of information at high resolution from satellite or other remote platforms. With this in mind, the simplest optical variable of a surface is its albedo (α), or the proportion of incoming short-wave radiation that is reflected. This value is important for a surface’s energy balance: high-albedo surfaces (e.g. new snow, α ≈ 0.7-0.9) reflect the majority of incoming short-wave radiation, but low- albedo surface (e.g. glacier ice, α ≈ 0.2-0.4) absorb a significant proportion of incoming radiation (Cuffey & Patterson 2010). While fresh snow can sometimes be considered a Lambertian scatterer which reflects incident radiation equally in all directions (e.g. Scambos & Haran 2002), metamorphosed snow and ice are known to be anisotropic reflectors (e.g. Choudhury & Chang 1981; Kuhn 1985; Knap & Reijmer 1998). Therefore, more specific definitions of albedo are necessary. Reflectance spectra measured in this study, where incoming light is measured from all directions but reflected light is measured at a point, are hemispherical-directional reflectances (HDR, Schaepman-Strub et al. 2006). With assumptions about the surface’s BRDF, HDR can be extrapolated to bihemispherical reflectance (BHR, Schaepman-Strub et al. 2006), where incoming and outgoing radiation are both considered in all directions. For reference, the BHR corresponds to the “white-sky” albedo while directional-hemispherical reflectance (DHR) corresponds to “black-sky” albedo (Lucht et al. 2000). As discussed above, understanding a glacier’s albedo is important for both energy balance and classification, but changing reflectance has also been implicated in glacier melt-albedo feedbacks. That is, as it receives shortwave radiation, a glacier surface begins to melt. This, in turn, increases the snowpack’s effective grain size and therefore reduces the albedo. A further reduction in albedo causes more melting, thus initiating a positive feedback cycle which is not broken until the next fresh snowfall. Such a melt feedback can be further exacerbated near the snowline as the surface changes from snow to firn to ice. This melt-albedo feedback has been observed on the Greenland Ice Sheet (Box et al. 2012), and has been suggested as a possible mechanism for enhanced thinning of higher elevations of glaciers in Svalbard (James et al. 2012). Interestingly, a reverse, negative feedback mechanism caused by increased precipitation has 21 Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review been suggested as a mechanism for inhibiting melt of the interior Antarctic Ice Sheet (Picard et al. 2012). 2.7.1 Estimation of Albedo from Landsat Data Many different sensors have been used to calculate glacier albedo, for example MODIS (Greuell & Oerlemans 2004; Greuell & Oerlemans 2005b; Stroeve et al. 2006; Greuell et al. 2007; Dumont et al. 2012), AVHRR (Greuell & Oerlemans 2004; Greuell & Oerlemans 2005b), and even ground-based oblique photography (Dumont et al. 2011). However, because of the extensive and radiometrically consistent and calibrated Landsat archive (Markham & Helder 2012), this thesis will focus on Landsat. Field studies have been used to calculate a narrow to broadband (NTB) conversion from Landsat band reflectance (see Section 3.3) to albedo using Band 2 (green) and Band 4 (NIR). Originally done by Knap et al. (1999b) with data from Morteratsch Glacier, this was updated with a larger data set including Iceland, Greenland, and Antarctica (Greuell et al. 2002): α = 0.539rLandsat2 + 0.166rLandsat4(1 + rLandsat4) (2.3) where r is reflectance in the indicated multispectral band. With this NTB, the residual standard deviation is just 0.011 (Greuell et al. 2002). The parameterisation works better with lower albedos than higher albedos, although there is no clear reason presented. Multiple studies have shown that one unified parameterisation is valid for the range of snow and ice facies found on glaciers (Knap et al. 1999a; Greuell et al. 2002; Hendriks & Pellikka 2004). The largest difference between measured and modelled albedos is seen, as expected, in spatially heterogeneous areas (Knap et al. 1999a). Many studies have used Equation 2.3 for Landsat NTB albedo calculation. Further considerations include a topographic correction to account for the varying amounts of incident radiation (Knap et al. 1999a; Klok et al. 2003; Hendriks and Pellikka 2004), ensuring removal of shaded and saturated pixels, and including diffuse radiation with a radiative transfer model (Klok et al. 2003). However, after NTB calculation, the most common correction to improve albedo calculation is using a BRDF parameterisation to account for anisotropic reflection by glaciers surfaces. Equation 2.3 refers to hemispherical surfaces only and therefore a consideration of anisotropy does need to be included (Knap et al. 1999b; Knap & Reijmer 1998). Indeed, Knap et al. (1999a), who assumed Lambertian reflectance, invoke anisotropy as the largest likely source of error, estimating a potential error of 0.15. Similarly, Hendriks and Pellikka (2004) implicate anisotropic reflectance in the persistent underestimation of reflectance from satellite observations. Studies have shown that the higher the solar zenith angle (the angle between zenith and the sun), the less uniform the snow’s reflectance is (Choudhury & Chang 1981); this has been confirmed in both polar and alpine areas (Kuhn 1985; Winther 1993). More recent work has confirmed that glacier ice is a forward-scattering surface, and that the type of ice (i.e. wet, honeycombed, etc.) does not need to be known to be accurately accounted for (Knap & Reijmer 1998). In applying corrections, however, it is important not to employ any circular logic; reflectance should not be used to classify glacier surfaces to then apply a BRDF used to re-calculate reflectance (Greuell et al. 2002). Earlier NTB albedo calculations with Landsat assumed Lambertian reflectance with some success Chapter 2: Glacier Facies, Multispectral Remote Sensing, and Reflectance – A review 22 (e.g. Knap et al. 1999a). De Ruyter de Wildt et al. (2002; 2003) used BRDF measurements for snow from Koks (2001), but chose to use a Lambertian assumption on glacier ice due to the presence of volcanic ash and tephra. For clean ice, Greuell and de Ruyter de Wildt (1999) provide BRDF parameterisations for Landsat bands 2 and 4; all examples exhibit nadir-darkening and forward-scattering. The function covers a wide range of zenith angles (25-75°) and varies with albedo itself, thus requiring a numerical solution; a simplification is provided for low slope surfaces and near-nadir view angles. Klok et al. (2003) use BRDF corrections for both snow (Koks 2001) and ice (Greuell & De Ruyter de Wildt 1999) to investigate the temporal and spatial variation of the surface albedo of Morteratschgletscher, Switzerland, as derived from 12 Landsat images; their anisotropy corrections ranges up to 0.10, depending on wavelength, solar zenith angle and surface type. Although for MODIS rather than for Landsat, Dumont et al. (2012) employed multiple simulations and radiative transfer modelling to observe a significant improvement in glacier albedo calculation when using a BRDF correction for snow based upon the data presented in Dumont et al. (2010). A more in-depth consideration of the corrections which need to be accounted for in this study will be discussed and investigated in Chapter 7. 23 Chapter 3: Field Sites and Data 3. Field Sites and Data This chapter begins with background on the two field sites addressed in this thesis: Langökull, Iceland and the glaciers of the Brøgger Peninsula, Svalbard (see Figure 3.1). Remote sensing data sources (Landsat and Airborne Thematic Mapper imagery) are described and pre-processing steps are detailed. Collection and processing of in situ field spectra are explained, and the data for both Iceland and Svalbard are presented. The chapter concludes with a small data set available for validating glacier surface classification using Landsat imagery in Svalbard. 3.1 Langökull, Iceland Langjökull (Icelandic for ‘Long Glacier’) is Iceland’s second largest ice cap (~950 km2; Björnsson & Pálsson 2008) and oriented SW-NE in central western Iceland. With an equilibrium line altitude at ~1000 metres above sea level (m a.s.l.), the ice cap has six major outlet glaciers, two of which, Hagafellsjökul Vestari (West Hagafells Glacier) and Hagafellsjökul Eystri (East Hagafells Glacier) in the south, have Figure 3.1: Regional maps and false-colour Landsat ETM+ images of (a) the Brøgger Peninsula’s Austre Brøggerbreen (AB), Midtre Lovénbreen (ML), and Austre Lovénbreen (AL) and (b) Langjökull, Iceland. Landsat bands 4 (NIR), 3 (red), and 2 (green) are used for red, green, and blue respectively in the false-colour composite. 0 km 50 100 150 Longyearbyen Barents Sea78°N 79°N 80°N 77°N Austfonna Ny-Ålesund 0 50 100 150 km Langjökull Vatnajökull Mýrdalsjökull Eyjafjallajökull Drangajökull Hofsjökull 64°N 65°N 66°N Arctic Circle 24°W 22°W 20°W 18°W 16°W North Atlantic Ocean Reykjavík a) b) AB ML AL 0 5 10 15 km 0 1 2 3 kma) b) Chapter 3: Field Sites and Data 24 records of surging behaviour (Sigurðsson 1998; Björnsson et al. 2002). Surface slopes range from ~5° on outlet glaciers to ~1° near the summit, and the average slope is 3.4° (Pope et al. 2013). The surrounding topography generally consists of exposed volcanic bedrock interspersed with glaciofluvial outwash plains and moraines. Termini of Langjökull have been scientifically observed since 1933; as is typical Iceland-wide, Langjökull’s surge-type glaciers are dominated by their surge activities, while on a whole there has been a large retreat through much of the 20th century, with a slowing in the retreat beginning in the 1970s (Sigurðsson 1998). Both in situ GPS and interferometric synthetic aperture radar have been used to measure flow speeds of Langjökull, ranging from 0 m yr-1 at the ice divides up to a maximum speed of 75 m yr-1 (Palmer et al. 2009; Gourmelen et al. 2011). Volume change has been calculated from historical maps and digital elevation models back to 1937, and the University of Iceland has been running an in situ mass balance programme on Langjökull since 1996 (Pálsson et al. 2012; Jóhannesson et al. 2013). From 1996 to 2001, Langjökull lost 5.73 m w.e. or ~3% of its total mass (Björnsson et al. 2002), and from 1996 to 2009 the mass balance was -1.3 m w.e. yr-1 (Pálsson et al. 2012). The ice cap is predicted to completely disappear around 2140 (Björnsson & Pálsson 2008); this loss is correlated with low snow accumulation and high annual temperatures (Björnsson et al. 2002). In order to conserve its current state, Langjökull must experience increased winter precipitation and/or decreased annual temperatures (Flowers et al. 2007). Less studied than many other glaciers in Iceland, Langjökull was chosen for study because of available airborne multispectral data and ease of fieldwork access. 3.2 Brøgger Peninsula, Svalbard With proximity to the Ny-Ålesund research station, the glaciers of the Brøgger Peninsula (Brøggerhalvøya) in northwest Spitsbergen, Svalbard – especially Austre Brøggerbreen (AB; East Brøgger Glacier), Midtre Lovénbreen (ML; Middle Lovén Glacier), and Austre Lovénbreen (AL; East Lovén Glacier) – have been well studied over recent decades (see Figure 3.1). ML is an alpine-type polythermal glacier, less than 5 km long and ~5 km2 in area. Mass balance has been measured along a central stake line by the Norsk Polarinstitutt since 1967 (Rippin et al. 2003; Arnold et al. 2006; Kohler et al. 2007; Rees & Arnold 2007). AB is predominantly cold-based and has a central flowline of ~6 km and an area of ~10 km2. While ML has a main trunk with a couple of small tributaries, AB is significantly wider and has three separate and significant tributary basins (Björnsson et al. 1996; Barrand et al. 2010). AL is somewhat less studied than its two neighbours, ~4 km long and ~5 km2 in area, similar in size and shape to ML (Friedt et al. 2012). All three glaciers have a similar altitude range of ~50-650 m a.s.l. and flow generally northward toward Kongsfjord onto the coastal plain (Kohler et al. 2007; Barrand et al. 2010; Friedt et al. 2012). ML’s mass balance, surface energy balance, and topographic change have been studied throughout the 20th century (Arnold et al. 2006; Rees & Arnold 2006; Rees & Arnold 2007; Kohler et al. 2007; Rye et al. 2010). Barrand et al. (2010) studied the variability of elevation change of both ML and AB and compared geodetic and traditional mass balance measurements; their mass balance has also been reconstructed back to 1948 (Rasmussen & Kohler 2007). James et al. (2012) observed enhanced thinning 25 Chapter 3: Field Sites and Data of higher elevations on ML, AB, and other Svalbard glaciers. The volume change of AL has also been studied for the second half of the 20th century, and it has been the focus of hydrological studies during the International Polar Year (2006-2010) because of its confined catchment area (Friedt et al. 2012). Evolution of snow cover on AL has been monitored with oblique, ground-based time-lapse cameras (see Section 3.6; Laffly et al. 2012). The Brøgger Peninsula was chosen as a study area because of available multispectral imagery, extensive interest in the region by the research community, and ease of field access to the area. 3.3 Landsat Imagery Since the launch of Landsat 1 in July 1972, NASA’s Landsat missions have long been dominant multispectral imagers (Loveland & Dwyer 2012), and, in a boon to Earth Observing science, their data are now freely available online (Wulder et al. 2012). Many articles have been written on the multiple Landsat missions, but a good first reference is Loveland and Dwyer (2012). Landsat data are only increasing in value as the programme extends its long archive of radiometrically consistent data (Chander et al. 2009; Loveland & Dwyer 2012; Markham & Helder 2012). This study uses data beginning in 1976 with the Multispectral Scanning System (MSS; Helder et al. 2012b), as well as data from the Thematic Mapper (TM: Helder et al. 2012a) and Enhanced Thematic Mapper Plus (ETM+; NASA 2011). ETM+ from 2003 onwards is subject to data gaps from the Scan Line Corrector failure, resulting in omission of 22% of the image (NASA 2011), although radiometric performance is not affected (Markham et al. 2005). Field spectra are also used to simulate data collected with the Operational Land Imager (OLI; Irons et al. 2012; Wilson & Oreopoulos 2013) on Landsat 8. See Table 2.1 and Figure 2.5 for more information on the spatial and spectral capabilities of the Landsat sensors. As Chapter 2 described, Landsat data have been extensively used for glaciological studies including surface classification. More specific to reflectance and albedo, (near-)coincident reflectance measurements have been conducted in Greenland (Hall et al. 1990), sub-Arctic Canada and Alaska (Hall et al. 1992), Svalbard spring (Winther et al. 1999; Casacchia et al. 2001), Svalbard summer (Winther 1993), and in Antarctica (Boresjö Bronge & Bronge 1999). While coincident in situ measurements and georeferenced remote sensing images offer the ability to directly compare reflectance or albedo, they will still suffer from issues related to the footprint size of each measurement (i.e. 1 m field of view for in situ vs. 30 m pixel). This study is concerned with surface class reflectance and identifying similarities between spectral properties. Physical processes leading to the maturation of a glacier’s surface will be similar across melt seasons and across years, and comparisons are based on spectra rather than a pixel-based validation. Therefore, although there may be different proportions of surface classes between image capture and in situ measurement, using older imagery to obtain a complete glacier image is acceptable. This also allows for a larger variety of images to select from and therefore avoid cloud cover, non-ideal solar and viewing geometry, and data gaps. Chapter 3: Field Sites and Data 26 3.3.1 Obtaining Data Landsat ETM+ imagery for this study was obtained by search and download through NASA ECHO Reverb (http://reverb.echo.nasa.gov/) and the United States Geological Survey’s EarthExplorer service (http://earthexplorer.usgs.gov/). For preliminary classification and comparison with reflectance spectra (see Chapter 5), in order to see the full range of melt facies, images of Langjökull and Midtre Lovénbreen well into the melt season were selected. For Langjökull, an ETM+ image from 20 August 2000 (Path 219, Row 15) was chosen. For ML, AB, and AL, an ETM+ image from 26 July 2000 (Path 220, Row 3) was chosen. Images previous to the ETM+ scan line corrector failure are used to ensure complete scenes. It is acknowledged that a decade passed between satellite imaging and spectra measurements, but, as discussed above, a spectral comparison will still be valid and robust. Images for the timeseries of surface classification and albedo measurements of the Brøgger Peninsula glaciers (see Chapter 7) were selected from the entire Landsat archive. All cloud-free images were initially selected, but those in early spring and late autumn are omitted due to the very high solar zenith angle and concordant long shadows, which obscured much of the terrain. Table 3.1 details the date, sensor, and scene of all images that are used. 3.3.2 Pre-Processing Landsat images analysed in Chapter 5 were converted to HDR at nadir using dark-pixel correction (Rees 2006) and radiometric constants provided by the Landsat Science Data Users Handbook (NASA 2011). The radiometric constants provided in metadata files alongside Landsat imagery have been revisited in recent years to ensure that the entire Landsat record is radiometrically calibrated and consistent; this was ensured using pseudo-invariant ground control points, and radiance is reported to be accurate to within ~10% uncertainty (Markham & Helder 2012). Imagery used for albedo calculations were subjected to a more rigorous pre-processing workflow. The LEDAPS atmospheric correction processing scheme was investigated (Masek et al. 2006), but it included some assumptions which were not appropriate for application here. Instead, the LandCor processing chain (Zelazowski et al. 2011) was used to apply the 6S radiative transfer code for atmospheric correction and calculation of at-surface reflectance. LandCor is a MATLAB processing chain that prepares whole images for the pixel-based 6S script, then interprets the results back into images. This facilitates atmospheric correction for whole scenes, where geometry and atmospheric conditions can change substantially. The process includes specification of ozone, water vapour, and aerosol content within the atmosphere by using available measurements or parameterizations (e.g. from NASA’s TOMS, OMI, or MODIS instruments) and accounts for variations in view angle and zenith angle across the Landsat scene. Additionally, a BRDF must be specified, while solar and sensor geometry is determined from the Landsat metadata file. For this analysis, total water vapour and aerosol optical depth are derived from the “MODIS/ Aqua Aerosol, Cloud and Water Vapor Subset 5-Min L2 Swath 5 km and 10 km V5.1” product; it is used to create seasonal lookup tables for the region of interest. All three properties are scaled based on elevation and month. The digital elevation model (DEM) used for this purpose was a 1995 DEM from 27 Chapter 3: Field Sites and Data the Norsk Polarinstitutt, derived from six stereo-overlapping aerial photographs taken in August 1995 (Rippin et al. 2003; Kohler et al. 2007; Friedt et al. 2012). Where available, the analysis uses the ozone estimate measured on the day of the Landsat image, but when these data are unavailable a seasonal lookup table has been developed with the rest of the instrument’s archived data. All Landsat images were georeferenced to a common reference frame, and the DEM was georeferenced and resampled to match this frame. The above inputs were processed through LandCor in MATLAB, provided to the 6S atmospheric Table 3.1: Landsat images used for surface classification and albedo timeseries of glaciers on the Brøgger Peninsula, Svalbard (see Chapter 7 for use). Date 25/04/1976 28/04/1976 02/05/1976 03/05/1976 01/06/1976 03/06/1976 09/07/1976 09/07/1976 10/07/1976 16/09/1988 10/07/1999 21/06/2000 24/06/2000 26/07/2000 15/06/2001 30/06/2002 07/07/2002 11/07/2002 11/07/2002 08/07/2005 25/06/2006 19/07/2006 23/07/2006 26/07/2006 12/06/2007 09/06/2009 25/06/2009 20/06/2010 06/07/2010 20/08/2010 23/06/2011 25/07/2011 05/08/2011 17/08/2011 Sensor L2 MSS L2 MSS L2 MSS L2 MSS L2 MSS L2 MSS L2 MSS L2 MSS L2 MSS L4 TM L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L5 TM L5 TM L5 TM L7 ETM+ L5 TM L5 TM L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ L7 ETM+ Path 235 238 242 243 237 239 2 239 240 223 218 215 220 220 216 220 221 217 217 220 220 220 216 221 220 220 220 220 220 215 220 220 217 221 Row 3 3 2 2 3 3 3 3 2 2 3 4 3 3 4 3 3 3 4 3 3 3 4 3 3 3 3 3 3 4 3 3 3 3 SLC Off Off Off Off Off Off Off Off Off Off Off Chapter 3: Field Sites and Data 28 model that calculates the atmospheric contribution to measured radiance, and these results were then used by LandCor to calculate at-surface reflectance for each Landsat band. 3.4 Airborne Thematic Mapper Campaigns The Airborne Thematic Mapper (ATM) was designed to mimic many of Landsat 7 ETM+’s bands but with added spectral coverage and spatial resolution (30 vs. ~2 m, determined by flight height and processing; see Table 2.1 and Figure 2.5). In particular, the ATM’s near-infrared band (0.91-1.05 μm) is interesting because snow reflectance in this range is largely dependent on a snow surface’s physical characteristics (e.g. grain size, melting, metamorphism; see Section 2.4). The ATM, Landsat, and other multispectral sensors have much in common, and there is much that can be learned regarding their utility in glacial remote sensing with respect to their relative spectral and spatial capabilities. The UK Natural Environment Research Council (NERC) Airborne Research and Survey Facility (ARSF) flew campaigns over Langjökull, Iceland and Midtre Lovénbreen on 2 August 2007 and 9 Figure 3.2: False-colour ATM image of Midtre Lovénbreen collected on 9 August, 2003. ATM bands 4, 3, and 2 are used for red, green, and blue respectively. North 0 1 km 29 Chapter 3: Field Sites and Data August 2003 respectively. On these missions, the ATM was mounted inside the ARSF’s Dornier 228 aircraft. Simultaneously collected laser ranging data (Rees & Arnold 2007; Pope et al. 2013) were used to orthorectify the imagery. Azgcorr version 5.0.0, produced by Azimuth Systems UK and provided by the ARSF, was used to perform the orthorectification (Azimuth Systems 2005). The entirety of Langjökull was flown in 24 separate, overlapping strips. Inconsistencies were clearly visible across swaths of uniform reflectance; these result from the relatively high flying height over Langjökull, increased atmospheric interference, and differences in viewing geometry (e.g. Palubinskas et al. 2003). In addition, there were some inconsistencies in the georeferencing of the ATM imagery. Therefore, it was decided on Langjökull to use the ATM data only for visual inspection, to aid in interpreting Landsat imagery, and in planning fieldwork. Midtre Lovénbreen, on the other hand, fits almost entirely into one ATM swath (see Figure 3.2). Different flying height and clearer atmospheric conditions resulted in significantly reduced swath effects on the ATM data. ATM measurements are delivered as at-sensor radiance. No measurements of incoming radiation were available coincident with airborne data collection, and calibration to reflectance with pseudo-invariant off-glacier targets from Landsat imagery was attempted, but no nearby surface was found to be consistent in its reflectance. Therefore, ATM data were left as at-sensor radiance. These band values are used in Chapter 6 to investigate the impact of resolution on glacier surface classification. No surface anisotropy or slope corrections were implemented. The meteorological records at the nearby Ny-Ålesund research station also indicate that the glacier surface froze overnight and that frost deposition was likely. Angular crystals on the surface may increase surface reflectance (Casacchia et al. 2001). These combined effects may have some impact on glacier surface classification. However, as the ATM image is only compared to spatially degraded versions of itself (see Section 6.2), possible classification-confounding factors should not impact the results of this study. 3.5 In Situ Spectral Measurements Ground-truthing satellite data can be a tricky process because it is impossible to mimic the exact data being collected by the imager. Footprint sizes are different, atmospheric conditions are different, and band wavelengths are different. However, much of this can be overcome by judicious and extensive in situ measurements, atmospheric correction (see Section 3.3.2), and incorporation of spectral response functions for NTB conversions (see Section 4.1). For this thesis, HDRs (see Section 2.7) were collected during two field campaigns; as the review in Section 2.4 noted, to the author’s knowledge this is the first presentation of spectra from the full range of facies available on end-of-melt-season glaciers. In August 2010, data were collected on Midtre Lovénbreen. In August 2011, data were collected on the major western outlet of Langjökull. These two locations were chosen so as to include sampling of glaciers that had undergone different accumulation and melting histories (see Sections 3.1 and 3.2). Additional differences were introduced by the springtime eruption of Grímsvötn, which deposited windborne ash on Langjökull’s surface; such ash-deposition events occur in Iceland frequently. All spectra were measured with an ASD FieldSpec Pro on loan from the UK’s Natural Environment Chapter 3: Field Sites and Data 30 Research Council’s Field Spectroscopy Facility (FSF). Simply, the FieldSpec can be though of as a portable hyperspectral instrument. The FieldSpec has a spectral range of 350-2500 nm, a sampling interval of 1.4 nm (350-1000 nm) or 2 nm (1000-2500 nm), and a spectral resolution (FWHM) of 3 nm at 700 nm, 10 nm at 1400 nm, and 12 nm at 2100 nm; spectra are all resampled to 1 nm. The FieldSpec features one 512-element VNIR photodiode array (350-1000 nm) and two separate, thermoelectrically cooled, graded-index SWIR InGaAs photodiodes (1000-2500 nm). Spectra were collected following procedures suggested by the FSF (see Figure 3.3; MacArthur 2007c; MacArthur 2007b; MacArthur 2007a; Milton et al. 2009; Goetz 2012). All spectra were collected with a lens-less 18° foreoptic from a height of ~1 m yielding a footprint with a radius of ~32 cm; the lack of a lens in the foreoptic also removes the need for a scrambler to compensate for preferential focusing of the viewing area on particular sensor fibres (MacArthur et al. 2012). Flat and homogeneous surfaces were chosen for measurement. Surface radiance and radiance of a Spectralon reference panel were collected in quick succession in order to calculate surface reflectance, and light stability levels were monitored with a handheld light meter. The FieldSpec has two data collection modes: “White Reference” mode and “Raw DN” mode. In White Reference mode, the spectrum of the Spectralon panel is measured and stored in memory Table 3.2: Descriptions of the 7 surface classes measured on Midtre Lovénbreen and 16 surface classes measured on Langjökull. See Figure 3.4 for BHR (reflectance) spectra. The numbered classes on Langjökull are the result of classes that were observed to be similar in the field but had sufficiently distinct spectra. Glacier Midtre Lovénbreen Langjökull Surface Class Fine Snow Coarse Snow Slush ‘Dry’ Ice Dry semi-saturated ice Wet semi-saturated ice Saturated Ice New Drifted Snow 1, 2, 3, 4 Ashy Snow 1 & 2 Ashy Snow Pit Side Ashy Snow Pit Bottom White Ice 1 & 2 Wet Ashy Ice 1 & 2 Debris Ice Very Debris Ice Stream Ice Debris 1 & 2 Description Snow with grain size <2.5 mm Snow with grain size >2.5 mm Saturated snow of any grain size Glacier ice with a friable surface and minimal liquid water content Glacier ice with a visible but small amount of liquid water on the surface Glacier ice with a significant liquid water content on the surface Glacier ice which is completely saturated with water on the surface, either streaming or very shallow ponding Snow which fell during the period of data collection, redistributed by wind Snow from the preceding summer with a dusting of volcanic ash which had not ablated into large pits Snow from the preceding summer on the side of the large pits caused by differential melting resulting from volcanic ash deposition with a moderate ash concentration Snow at the bottom of ash-induced melt pits (up to 2m deep) with the highest ash concentration measured Glacier ice with a friable surface, minimal liquid water content, and no apparent ash Glacier ice with some liquid water content and apparent volcanic ash Glacier ice with some non-ash debris Glacier ice with the majority of the surface covered by non- ash debris Glacier ice overlain with very shallow flowing meltwater and accumulated debris 31 Chapter 3: Field Sites and Data for subsequent sample spectra such that the only data that are permanently recorded are reflectance. However, in Raw DN mode, the radiance of the Spectralon panel and the sample spectra are stored in independent files. Reflectance is then calculated in post-processing using tools provided by the FSF. All data were collected in Raw DN mode. Spectra from Midtre Lovénbreen were processed using the Excel worksheet provided by the NERC FSF (2013), and spectra from Langjökull were calculated with the then newly available MATLAB User Group Toolbox (Robinson & MacArthur 2011). See Section 4.6 for more information on the software used in this thesis. A variety of locations were selected on each study glacier for reflectance measurements. Based upon a preliminary investigative survey of the study area, multiple locations were selected for what were determined to be the area’s representative facies based on snow type, water content, ice condition, and debris content (see Table 3.2, Figure 3.5). Fresh snow can be a Lambertian (uniform) reflector (Scambos & Haran 2002), and while melting ice is known to be a somewhat forward-scattering surface (Knap & Reijmer 1998), for calculation of BHR from HDR (see Section 2.7) across the range of surface classes we assume Lambertian behaviour. For each location, to mitigate potential inhomogeneous reflectance, 5 samples at -90°, -45°, 0°, 45°, and 90° to the sun’s azimuth were averaged together (where at 0° the observer is pointing the sensor arm directly towards the sun; see Figure 3.3); these were duplicated on Midtre Lovénbreen. In addition, to increase signal strength and accuracy, each ‘individual’ sample collection was actually a rapid averaging of 50 individual spectra in a few seconds. After processing, spectra were grouped into distinguishable surface types based upon the mean and standard deviation Figure 3.3: Measurement setup for collecting surface reflectance spectra as seen on Langjökull. Note the FieldSpec in the backpack, controlling laptop on a belly board, sensor attached to the end of a pole, and Spectralon white-reference panel on the tripod in the foreground. Angles indicate the relative position of the operator for multiple measurements of the same sample. 0° 45° 90° -45° -90° θ ,Solar Azimuth S Chapter 3: Field Sites and Data 32 New Drifted Snow 1 Wet Semi-Saturated Ice ‘Dry’ Ice Coarse Snow New Drifted Snow 2 Saturated Ice Dry Semi-Saturated Ice Slush Fine Snow New Drifted Snow 4New Drifted Snow 3 Table 3.2 (cont.): Images of the surface classes defined in the earlier table, presented in the same order. Photos show people, tripods, shoes, and a fieldwork mascot to provide a sense of scale. 33 Chapter 3: Field Sites and Data White Ice 1 White Ice 2 Stream Ice Debris 2 Very Debris Ice Wet Ashy Ice 2 Ashy Snow Pit Bottom Ashy Snow 2Ashy Snow 1 Ashy Snow Pit Side Wet Ashy Ice 1 Debris Ice Stream Ice Debris 1 Chapter 3: Field Sites and Data 34 a) b) Figure 3.4: In situ surface reflectance spectra from (a) Midtre Lovénbreen and (b) Langjökull. The thickness of each line in (a) includes error bars based on the combination of many spectra in each class; error bars are spaced every 100 nm for (b). Data in the ranges 1350-1460 nm and 1790-1960 nm have been removed due to confounding water absorption in these wavelengths. Horizontal bars in (a) are an example of weighted averages of the continuous data to mimic bands of Landsat 7’s ETM+ (see Section 4.1). 35 Chapter 3: Field Sites and Data of each in situ data collection group; if spectra were within the standard deviation of another group across all wavelengths, the groups were merged. Similarly, if spectra fell outside their mutual standard deviations, groups were split: on Midtre Lovénbreen, snow surfaces were divided into coarse snow and fine snow; on Langjökull, surfaces which appeared to be the same in situ but whose spectra were not similar (as determined by standard deviation) are those which have been given numbers. 3.5.1 Midtre Lovénbreen data In situ surface reflectance spectra from Midtre Lovénbreen (see Figure 3.4a) behave in a very orderly manner as expected from model results as well as previously published snow and ice reflectance spectra (see Section 2.4; Wiscombe & Warren 1980; Zeng et al. 1983; Hall et al. 1987; Hendriks & Pellikka 2004; Winther 1993; Casacchia et al. 2001; Gardner & Sharp 2010). All spectra show a relatively high reflectance from 350 to ~800 nm due to primary and multiple scattering before darkening in the NIR and becoming dark in the SWIR due to high absorption by water content. All spectra also show a peak in reflectance at ~1100 nm as modelled and observed in many other studies. Surface classes are described in Table 3.2. Fine snow is the brightest, followed by coarser snow and then ice as defined by varying levels of surface water content. Slush is darker than some icy surfaces but brighter than fully saturated ice. As the snow ages and melts, this increases the effective grain size, both by grain growth itself as well as through liquid water content as noted in Section 2.4. 3.5.2 Langjökull Data Both recent snow accumulation and drifted ash from a nearby volcanic eruption (Grímsvötn) in the spring contributed to a very variable glacier surface when collecting in situ reflectance spectra from Langjökull (see Figure 3.4b). Impurities on the snow surface (such as ash) will contribute to a darkening in the spectral region where ice absorption is weak (< 900 nm) (< 900 nm, Gardner & Sharp 2010; Kokhanovsky 2013), especially for coarse-grained snow (Warren & Wiscombe 1980). If it were of interest, with some assumptions, a radiative transfer model (e.g. Kokhanovsky 2013) could be used to calculate the concentration of particulate matter on the snow surface. An ash-induced melt feedback resulted in pits, which melted to almost 2 m below the surrounding snow surface; very deep pits showed a much darker ash layer, the remnant of the Eyjafjallajökull eruption the previous year (2010). The uneven ash surface contributed to a larger number of surface classes than on Midtre Lovénbreen. In addition, the larger variability in spectra led to some difficulty in defining particular surfaces, creating in some cases more of a continuum; the result is higher levels of uncertainty for some classes. More important than the individual spectra is their crossing behaviour; physically, this is likely indicative of contributions from different impurities and water content, which impact reflectance at different wavelengths. While the Midtre Lovénbreen spectra stay parallel and sub-parallel and thus maintain the same relative brightness to each other at all wavelengths, this is not the case on Langjökull. Langjökull spectra intersect and cross each other. These changes in spectrum slope and concavity are its defining features. Nevertheless, the drifted snow classes, because they have few impurities and small grain sizes thanks to a smaller degree of grain metamorphosis, are still the brightest surface classes. Chapter 3: Field Sites and Data 36 100 m 200 m 300 m 400 m500 m 0 1 2 km 6 sites: Dry Semi-Saturated Ice Wet Semi-saturated Ice Saturated Ice ‘Dry’ Ice 6 sites: Dry Semi-Saturated Ice Wet Semi-saturated Ice 3 sites: ‘Dry’ Ice 6 sites: ‘Dry’ Ice 3 sites: Fine Snow 6 sites: Fine Snow 9 sites: Coarse Snow Fine Snow Slush 1300 m 1200 m 1100 m 1000 m 1 site: Ashy Snow 11 sites: Ashy Snow Pit Bottom Ashy Snow Pit Side New Drifted Snow Ashy Snow 9 sites: Ashy Snow Pit Bottom Ashy Snow Pit Side New Drifted Snow Ashy Snow 13 sites: Stream Ice Debris Very Debris Ice Wet Ashy Ice Debris Ice White Ice 0 1 2 3 4 km a) b) Figure 3.5: Locations of field spectra collected on (a) Midtre Lovénbreen and (b) Langjökull. Contours for (a) were derived from a 2005 LiDAR digital elevation model and for (b) from a 2004 SPOT-derived digital elevation model. 37 Chapter 3: Field Sites and Data One class, New Drifted Snow 1, is observed to be over 100% reflectance. Although seemingly impossible, there are two possible, physically based explanations for this observation. The first is that the irradiance values that correspond to the New Drifted Snow 1 spectra were taken at a slightly lower light condition; as light conditions were carefully monitored, this is unlikely. Therefore, it is almost certainly the second explanation: the FieldSpec measured some specular rather than purely diffuse reflectance. This would not occur when measuring the Lambertian Spectralon panel; thus observed radiance is larger than irradiance. 3.6 In Situ Facies Classification with Time-lapse Imagery On Austre Lovénbreen, observational work done by French colleagues has allowed for the unique opportunity to use ground-based measurements to evaluate the remote-sensing based classification performed in Chapter 7. A network of time-lapse cameras images the glacier three times per day, then the images are orthorectified and combined, and finally they are classified into accumulation or ablation areas. Evolution of snowcover from complete glacial cover to a minimum AAR and back to complete cover has been observed on Austre Lovénbreen (Laffly et al. 2012). A similar setup has been deployed by Huss et al. (2013) in the Swiss Alps. In addition, work has been done to move the data processing into the cloud to facilitate handling of the large amount of data (Ranisavljević et al. In Review). In this study, three days coincident with available satellite images have been processed for direct comparison with Landsat- based glacier surface classification (20 August 2010, 23 June 2011, and 17 August 2011, see Chapter 7). 38 39 Chapter 4: Methods 4. Methods This chapter contains information about the methods used in processing the data presented in Chapter 3 so that it can be combined and analysed in Chapters 5, 6, and 7. Topics include narrow-to- broadband (NTB) conversions from reflectance spectra to multispectral bandwidths and multispectral band reflectance into albedo. Also discussed are methods to reduce data dimensionality while retaining information, in particular principal component analysis (PCA) and various classification techniques. 4.1 Spectral Response Matching A NTB conversion is necessary to compare the nearly continuous FieldSpec data with the individual multispectral bands presented in Figure 2.5. This is done with a known relative spectral response function for each band (at 1 nm resolution, like FieldSpec data) obtained via NASA, ESA, and the NERC FSF; relative spectral response functions are only available for MODIS Bands 1 through 16, so higher bands are omitted in further analyses. Sentinel-2 and OLI are both approximated with a “top hat function” (i.e. uniform spectral response across each band) as no further data were available at the time of writing. The relative spectral response function describes the extent to which the multispectral sensor responds to light of a particular wavelength and is therefore multiplied by an in situ reflectance spectrum to produce a weighted average – in other words: individual, band-specific BHRs. Presented in equation form: rnb = (4.1) where rnb is the narrowband reflectance, r(λ) is the spectral reflectance, R(λ) is the relative spectral response, and λ is wavelength (Greuell & Oerlemans 2004). Integrals are presented from 0 to infinity as a mathematical convenience. Shortly outside the defined bandwidths presented by multispectral sensors, the spectral response drops to negligible values from one main peak; a notable exception in relative spectral response is Band 1 of the Airborne Thematic Mapper, which has a small second peak at 800 nm (see Figure 4.1). Thus, using the real rather than theoretical range of relative spectral response, the calculation of rnb becomes easily tractable. These rnb values can then be directly compared to the reflectance of (for example) classified ETM+ imagery (see Chapter 5). 4.2 Multispectral Sensor Radiometric Properties Beyond the NTB conversion itself, further calibration parameters from each multispectral band must be taken into account to fully simulate the measurements which would be taken by a multispectral sensor – had the user been holding an extremely portable version of a satellite rather than the FieldSpec’s 0 ∫ ∞ r(λ)R(λ)dλ 0 ∫ ∞ R(λ)dλ Chapter 4: Methods 40 photodiodes. Each sensor is limited by the minimum and maximum radiance values which it can measure (LMIN and LMAX respectively); from the digital number (DN) provided as data downloads, the at-sensor spectral radiance (Lλ) can then be calculated: Lλ = LMIN + (Q – QMIN) (4.2) where Q, QMIN, and QMAX are the DN as measured, minimum DN, and maximum DN respectively. All Q values are integers, and the range of Q defines the radiometric resolution of the sensor. Although field spectra (and some imagery, i.e. ATM) are already provided as at-sensor radiance, this step is important for understanding the linear interpolation inherent in recording remote sensing imagery. At-sensor spectral reflectance (ρλ) is then calculated from at-sensor radiance: ρλ = (4.3) where d is the Earth-Sun distance in astronomical units, ESUN,λ is the mean solar exoatmospheric solar irradiance, and θz is the solar zenith angle; the inclusion of π is the result of an assumption of Lambertian surface reflectance when estimating the influence of solar irradiance. Using Equation 4.3, the range of measureable reflectance values can be calculated by substituting LMIN and LMAX for Lλ:                       Figure 4.1: The relative spectral response functions for Bands 1-10 of the Airborne Thematic Mapper. Note Band 1 (heavy line), which has an anomalous second peak; all other bands have one well-defined peak. The grey bars are the bands as defined by the manufacturer. LMAX – LMIN QMAX – QMIN π Lλ d2 ESUN,λ cos(θz) 41 Chapter 4: Methods ρMIN = (4.4) ρMAX = (4.5) where LMIN and LMAX are properties of each multispectral band as well as the gain setting in which the sensor is configured. Low gain settings are used to increase the observable radiometric range. High gain settings, with a smaller radiometric range, are often used to observe darker regions. Both settings retain the same radiometric resolution, and so high gain images record finer gradations in brightness, albeit for a subset of brightness values. ESUN has a different value for each spectral band, but it is independent of gain setting. Just as radiance is linearly interpolated to a range of integers when recorded by a multispectral sensor, the field spectra must be rounded to a set of discrete values, thus quantizing the data. This is done using both the radiometric range (ρMIN to ρMAX) and radiometric resolution (QMIN to QMAX): For rnb < ρMIN, DNBAND = QMIN For ρMIN ≤ rnb ≤ ρMAX, DNBAND = || (rnb – ρMIN) + QMIN || (4.6) For rnb > ρMAX, DNBAND = QMAX where rnb is the narrowband reflectance calculated in Equation 4.1, ||x|| is notation for rounding x to the nearest integer, and DNBAND is the digital number which would be recorded by a multispectral sensor in place of the FieldSpec given the same signal input. DNBAND is easily converted back to a reflectance value, rBAND: rBAND = DNBAND (4.7) such that rBAND is reflectance for a particular sensor’s multispectral band emulated using field spectra. For this study, the radiometric calibration values for each sensor and setting are recorded in Table 4.1. QMIN, QMAX, LMIN, and LMAX are obtained from the data providers (ARSF 2010; NASA 2002; Drusch et al. 2010; Chander et al. 2009; NASA 2011; NASA 2013). ESUN is provided as a calibrated value for Landsat sensors, MODIS, and ASTER; for Sentinel-2 and ATM, other sensors ESUN is calculated using relative spectral response functions and the same solar spectrum used by the Landsat team for its calculation (see Equation 4.1; Thuillier et al. 2003). The Earth-Sun distance d is taken to be its mean value of 1. Solar zenith angles for both Midtre Lovénbreen and Langjökull are investigated using Landsat metadata files for a range of images. For Langjökull, θz = 60°, and for Midtre Lovénbreen, θz = 55°. Field spectra are not used to simulate Landsat-8 OLI data because calibration values were not available at the time of writing; the 12 bit radiometric resolution of the OLI is assessed by comparing the π LMAX d2 ESUN cos(θz) QMAX – QMIN ρMAX – ρMIN ρMAX – ρMIN QMAX – QMIN π LMIN d2 ESUN cos(θz) } Chapter 4: Methods 42 Sensor ATM ASTER Landat TM / ETM+ Landsat MSS MODIS Band (gain) 1 2 3 4 5 6 7 8 9 10 1 (normal) 1 (low 1) 1 (hi) 2 (normal) 2 (low 1) 2 (hi) 3 (normal) 3 (low 1) 3 (hi) 1 (L) 1 (H) 2 (L) 2 (H) 3 (L) 3 (H) 4 (L) 4 (H) 4 5 6 7 1 2 3 4 5 6 7 8 9 QMAX 65,535 65,535 65,535 65,535 65,535 65,535 65,535 65,535 65,535 65,535 254 254 254 254 254 254 254 254 254 255 255 255 255 255 255 255 255 255 255 255 255 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 LMAX 211.23 291.27 314.06 296.30 269.91 213.10 145.43 71.97 18.40 4.830 427 569 170.8 358 477 179 218 290 106.8 293.7 191.6 300.9 196.5 234.4 152.9 241.1 157.4 238 168.7 143.6 125.8 354.52 231.22 449.36 402.06 82.71 29.90 53.50 333.09 544.95 QMIN 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 LMIN 0.751 0.831 0.687 0.487 0.417 0.255 0.181 0.227 0.006 0 0 0 0 0 0 0 0 0 0 -6.2 -6.2 -6.4 -6.4 -5.0 -5.0 -5.1 -5.1 7.2 -0.6 -2.4 3.9 4.16 1.33 16.32 7.66 0.24 0.05 0.01 18.12 16.92 ESUN 1562 1935 1792 1671 1522 1326 1028 800.6 216.4 77.5 1847 1847 1847 1553 1553 1553 1118 1118 1118 1997 1997 1812 1812 1533 1533 1039 1039 1827 1569 1260 866.4 1577.9 971.4 2060.2 1839.2 454.6 239.7 98.9 1720.2 1875.9 43 Chapter 4: Methods effect of radiometric resolution on the full-spectrum data (see Section 6.1). To conclude, by performing “corrections” to the in situ spectral reflectance data, they can be used to simulate the radiometric and spectral capabilities of a suite of multispectral imagers while maintaining the basic radiometric properties as a control (see Chapter 6). 4.3 Principal Component Analysis The spectra measured in this study are highly multidimensional data, and multispectral data have a number of dimensions (in this case, bands) themselves. Therefore, some method is needed to reduce the dimensionality of the data for analysis – without losing important information – in order to understand what the most important input wavelengths are for glacier surface classification. A dataset may start with n dimensions, but the assumption is made that fewer than n dimensions are sufficient to capture the important information provided by the data. Information is considered to be ‘important’ in this context if it helps distinguish the input data points, whether pixels or spectra. In other words, some amount of autocorrelation (i.e. the reflectance at one wavelength is highly related to nearby wavelengths, effectively reducing the spectral resolution) is exploited to reduce the data’s dimensionality. Concisely put, PCA is a transformation that reduces the dimensionality of a dataset by reprojecting it into a new coordinate space (e.g. Boresjö Bronge & Bronge 1999; Sidjak & Wheate 1999). Each dataset is considered in a multidimensional space (whether 6 dimensions for ETM+ or 1000 dimensions for in Sensor MODIS Sentinel-2 Band (gain) 10 11 12 13 (lo) 13 (hi) 14 (lo) 14 (hi) 15 16 1 2 3 4 5 6 7 8 9 10 QMAX 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 4,095 LMAX 349.00 231.89 210.13 72.33 53.59 84.97 46.93 63.57 58.95 587.87 615.48 559.01 484.13 449.55 412.92 387.08 307.8 232.91 45.0 QMIN 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 LMIN 13.48 9.43 8.44 4.11 4.11 3.86 3.86 2.66 1.31 15.97 11.7 6.49 3.31 2.61 2.06 1.67 0.95 0.51 0.06 ESUN 1956.7 1857.2 1865.4 1517.4 1517.4 1474.9 1474.9 1278.8 955.8 1884.9 1969.3 1825.1 1545.2 1419.7 1289.9 1156.9 956.4 812.2 366.4 Table 4.1: Radiometric calibration values for a range of multispectral imagers. QMIN and QMAX are the minimum and maximum DNs which can be recorded by a sensor (i.e. radiometric resolution). LMIN and LMAX are the minimum and maximum radiance measureable by a sensor (i.e. radiometric range) and are reported in units of W m-2 sr-1 μm-1. ESUN is the mean solar exoatmospheric solar irradiance and is recorded in units of W m-2 μm-1. Chapter 4: Methods 44 situ spectra, in each case corresponding to the spectral resolution of the data), and a new coordinate space is defined. Thus, each principal component (PC) is a linear combination of the input data. Given n input parameters, n PCs will also be output. These PCs are ranked by their variance. The variance can be considered as a measure of the amount of important information contained within a given PC. The lower a PC, the higher the variance, and the more information it contains. For example, for an ETM+ image, the first principal component (PC1) describes most of the brightness information while PC6 will be largely noise. Crucially, PCs are uncorrelated. Put another way, reflectance of a pixel in Band 1 will most certainly be related to the reflectance in Band 2, while the value PC1 for that pixel will not be related to the value for PC2. Because PCs are not correlated, a variance threshold can be set, below which PCs can be ignored without losing valuable information. Thus, PCA can be used to reduce data dimensionality while retaining most of their information (gathering the data into fewer bands). As with multispectral bands, if field spectra are input, then PC1 will show the defining characteristics of the set while PC1000 will highlight noise. Once the primary PCs are identified, the factors used as coefficients for the linear combination of input components to create them can also be used to investigate and quantify the relative contribution of the input components, which is what is done here with the in situ spectra. For field spectra, the wavelengths with the largest magnitude coefficients are the most important for reflectance-based glacier surface classification. In this thesis, PCs are also used as input into unsupervised classifications. Andrews Plots (Spencer 2003), a method for graphic continuous multivariate data, were also considered as a method to either replace PCA or aid in the analysis of the linear combinations output by PCA. The use of Random Forest Classification (Breiman 2001) was also investigated for reducing data dimensionality; this method harnesses an ensemble of individual tree classifiers, each of which uses random inputs sampled independently and with the same distribution, all correlated with each other. Attempts at using these techniques, however, were not successful with the field spectra because of the very autocorrelation within the data that PCA exploits. This is just the kind of situation suited to PCA, therefore it is a valuable tool for analysing and fully utilizing the full range of data which overcomes autocorrelation without discarding complexity across the spectra. 4.4 Classification For comparison with in situ spectra in Chapter 5, ETM+ imagery was segmented into radiometric classes. Glaciers were roughly manually outlined, and within these polygons an unsupervised classification was performed (as noted in Section 2.6.3) on inputs of the Normalized Difference Snow Index, a ratio of Landsat bands 4 (NIR) and 5 (SWIR), and the second, third, and fourth principal components according to a scheme from de Woul et al. (2006) and Braun et al. (2007). These PCs were chosen by previous studies in order to minimize the visible bands (which sometimes saturate out) while maximizing the NIR and SWIR (which have lower reflectance for snow and ice). For each image, 10 separate classes (a manageable, arbitrary number) were identified with a clustering algorithm. Previous studies used ISODATA (Iterative Self-Organising DATa Analaysis; see examples in Section 2.6), an algorithm in which candidate clusters are initially (randomly) selected and their means are allowed to migrate in the spectral domain, optimizing the classification with each 45 Chapter 4: Methods iteration; each pixel is re-compared to each cluster mean, then assigned to the nearest cluster. Means are recalculated and the process is repeated until the pixel-to-centre distance means are minimized or a set number of iterations is reached. ISODATA algorithms are easily implemented in a variety of software environments. A similar clustering algorithm, also widely used, is k-means clustering (see examples in Section 2.6). Similar to ISODATA, k-means clustering is an iterative scheme which attempts to minimize the distance between points labelled to be in a cluster and the centre of that cluster. Using ISODATA, thresholds for merging and splitting classes can be used such that the output number of classes is not predefined, but k-means requires the number of output clusters to be pre-specified. A similar scheme, k-medoids, uses a point as the centre of a cluster rather than the average of a group. K-means clustering is easier to implement than ISODATA in some software environments which are not tailored to remote sensing imagery. In Chapter 5, classes are grouped using an ISODATA algorithm. Rather than combine classes (using either geographical or spectral knowledge; de Woul et al. 2006; Braun et al. 2007), each class is preserved and pixel values are averaged for radiometric comparison with in situ reflectance data. Chapters 6 and 7 also classify glacier surfaces using ATM imagery and Landsat (MSS, TM, and ETM+) data respectively. ISODATA is used in Chapter 6, but k-means clustering was used in place of ISODATA in Chapter 7 for ease of implementation. Instead of the inputs described in the first paragraph of this section, linear combinations presented in Chapters 5 and 6 are used. The 10 output classes were subsequently grouped into accumulation and ablation areas for statistical analysis (see Section 6.2), for comparison with time-lapse imagery (see Section 3.6), and used for calculation of albedo in glaciologically meaningful and repeatable areas (see Chapter 7). In addition to spatial classifications, field data are clustered based upon their spectral reflectance functions. To cluster the in situ data and derived values which simulate multispectral sensors, k-means clustering is used in Chapters 5 and 6. Linear combinations defined in each chapter are used as inputs to the classifications. 4.5 Statistical Classification Comparison This study requires that two main classification accuracy assessments be conducted. The first is a clustering analysis where the classification of field spectra is compared to knowledge from fieldwork (see Section 6.1). The second compares the results of classification of ATM imagery degraded to different spatial resolutions (see Section 6.2). Classification accuracy is often assessed by comparing results derived from remote sensing imagery and from ground-truth fieldwork; control point sampling and spatial distribution therefore become important considerations (e.g. Comber et al. 2012). In this study, the first comparison does not include a spatial component, and the second comparison uses other imagery rather than field data for comparison. Therefore, sampling and regional distribution are not factors that need to be addressed. Minimum mapping units (MMUs) are compared on a one-to-one basis, creating in the first instance a contingency matrix which shows the number of times (i.e. number of pixels or number of spectra) where two classifications agree or disagree. Matrix values are often converted to percentages Chapter 4: Methods 46 for ease of understanding. The contingency matrix provides information on errors of both omission and commission and is the basis for normalizing the counts and for calculating classification similarity statistics (Congalton 1991; Congalton & Green 1999; Foody 2002; Rees 2008). All spectra or pixels of a remote sensing image are used so there is no sampling bias. Pixels are chosen as the MMU (i.e rather than polygons, Stehman & Wickham 2011) because that is the level at which the classification is carried out. However, comparison of square pixels across different resolutions could be considered as using blocks rather than pixels for classification. MMUs are assumed to be homogeneous (Stehman & Wickham 2011). Choice of MMU can have significant impact on change detection, and more care needs to be taken coarsening high resolution imagery than downscaling low resolution imagery (Colditz et al. 2012); this will be kept in mind in Section 6.2 where comparison across different spatial resolutions is addressed. Creation of the contingency matrix also assumes that output classes are indeed identical. By virtue of being collected on the ground, this is known to be the case for field spectra. In addition, all ATM classification comparisons are evaluating radiometrically identical datasets (see Section 6.2), so this is a reasonable assumption. Pixel-based comparisons also assume that images are collected coincidentally in time and space; again, because all images are based on the same original ATM image, this assumption is easily satisfied. Because of the controlled nature of these comparisons, this study avoids many of the potential pitfalls of a more applied study of classification accuracy assessment (Colditz et al. 2012). From the information in the contingency matrix, the most basic statistic is “A,” or the overall accuracy agreement. This is the sum of MMUs for which the classifications agree divided by the total number of MMUs sampled. Put another way, A is “the trace (sum of the terms in the leading diagonal) of the contingency matrix normalized to the sum of its terms” (Rees 2008). However, A is subject to agreement by random chance and therefore overestimates classification accuracy. In response, Cohen’s Kappa (K, Cohen 1960) accounts for random effects within the classification comparison and remains an indexed value (i.e. perfect agreement result in K = 1). K = (4.8) where A (same as above) is the proportion of units for which the classifications agree, and A* is the proportion of units for which agreement is expected by chance. For a given contingency matrix, this can be calculated as: where n is the number of rows and columns, crc is the number of observations in row r and column c, and N is the total number of observations. In this way, K reduces the overestimation of classification success included in A. While some believe that K may go too far and underestimate error (Foody 2002), K is nevertheless widely accepted by the community as a metric for classification accuracy assessment (Congalton 1991; Monserud & A = ∑ cii (4.9) A* = ∑ ∑ cij ∑ cji (4.10) 1 N 1 N2 n i=1 n i=1 n j=1 n j=1 A - A* 1 - A* 47 Chapter 4: Methods Lower Bound 0.00 0.05 0.20 0.40 0.55 0.70 0.85 0.99 Upper Bound 0.05 0.20 0.40 0.55 0.70 0.85 0.99 1.00 Degree of Agreement No Very Poor Poor Fair Good Very Good Excellent Perfect Table 4.2: Descriptive terms for classification based upon Cohen’s K values, following Monserud and Leemans (1992). Leemans 1992; Rees 2008). K values can also be considered with qualitative descriptions rather than simply numerical values (see Table 4.2). It is worth noting that K is suitable for this study because the number of classes is identical in the pairs of cases which are compared; if this were not the case, another metric such as Cramer’s V may be appropriate (Rees 2008). 4.6 Software Tools As with any modern-day research, a very wide range of software tools is indispensible. Rather than attempting an exhaustive list, this section aims to acknowledge those tools most crucial to the production of this thesis. Where practicable, an effort was made to choose freely available software for analysis steps, although widespread use by the community or particular functionalities meant that proprietary software was often necessary. Imagery was investigated and analysed mainly in software packages ImageJ (http://rsb.info. nih.gov/ij/), MultiSpec (https://engineering.purdue.edu/~biehl/MultiSpec/), and ERDAS IMAGINE 9.3. Geospatial tasks were addressed in both ArcMap 9.3.1 and QuantumGIS (http://www.qgis.org/). Data analysis, especially that requiring increased automation and customisation was done in MATLAB R2010a (an open-source alternative would be Octave, http://www.gnu.org/software/octave/) and R (http://www.r-project.org/). Simple tasks were tackled in Microsoft Excel. Beyond data analysis, presentation is, of course, a key part of communicating results. Figures were created in Plot (http://plot.micw.eu/) as well as Adobe Illustrator CS5 and Adobe Photoshop CS5. The document was written in Microsoft Word, composed in Adobe InDesign CS5, and Zotero (http:// www.zotero.org/) was used to manage references. 48 49 Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 5. Using In Situ Spectra to Explore Landsat Classification of Glacier Facies The most effective glacier surfaces classification schemes in the literature using Landsat data have not been developed using full-spectrum reflectance (see Section 2.6). The bands of multispectral imaging sensors were designed according to atmospheric transparency and many other applications’ specifications (e.g. agriculture), not with reference to the spectral properties of snow and ice; it is not self-evident at all that multispectral sensors should be optimal for remote sensing of glaciers. Therefore, this study approaches the problem from a different angle. Using in situ reflectance data collected over the full range of surface classes present on a glacier at the end of the melt season (see Section 3.5), this study aims to answer these questions: What parts of the visible-infrared (VIR, 350-2500 nm) spectrum are best-suited to distinguishing glacier facies? Can the ETM+’s restricted spectral data discriminate glacier facies as effectively as continuous reflectance spectra? What do these results mean for robust glaciological applications of surface classification? The spectra and Landsat imagery which are used here have been described in earlier sections (see Chapter 3). Initially, a qualitative analysis of the spectra and comparison with satellite imagery classification is used to identify salient features in the data and to understand some of the physical meaning behind the spectra (Section 5.2). The qualitative conclusions are built upon with a quantitative treatment of the spectra, using principal component analysis (see Section 4.3) to understand what wavelengths are most important for distinguishing surface classes (see Section 5.3). A clustering analysis (k-means, see Section 4.4) is used to investigate the robustness of this technique for classification (see Section 5.4). The full spectra are able to simulate what satellite data would be (see Section 4.1), and this is used to develop the basis for a transferrable classification method which can be used on any glacier (see Section 5.5). This chapter finishes by considering this technique with other published studies in terms of ease of use, success, and real-world application. 5.1 Landsat ETM+ Classification ETM+ images from both the Svalbard glaciers and Langjökull are classified (for method, see Section 4.4); the classification and band reflectance for each class are displayed in Figure 5.1. In both locations we see bright (i.e. highly reflective) classes at higher elevation and darker classes at lower elevations. Both locations show classes which are easily interpreted as non-snow and non-ice areas. The Svalbard glaciers also exhibit shadowing near their valley heads. While we can say that the brightest classes are linked to accumulation and classes near glacier termini are related to ablation, there is little compelling evidence to draw an equilibrium line between any particular classes. Therefore, in both locations, while basic spectral and geographical information can be used to identify glacier facies, confident identification of intermediate facies is challenging. Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 50 Fi gu re 5 .1 : U ns up er vi se d cla ss ifi ca tio n of L an ds at E TM + im ag er y of g la ci er s o n (a ) S va lb ar d’s B rø gg er P en in su la a nd (c ) L an gj ök ul l I ce C ap , I ce la nd a nd th e re fle ct an ce o f e ac h cla ss (b , a nd d , r es pe ct iv ely ). Ba sic g eo gr ap hi ca l a nd sp ec tr al in fo rm at io n ca n be u se d to re co gn iz e so m e cla ss es , b ut id en tifi ca tio n of in te rm ed ia te fa ci es re m ai ns u nr eli ab le. W hi le S W IR w av ele ng th s a re u se fu l f or d efi ni ng g la ci at ed an d no n- gl ac ia te d ar ea s, ab so lu te b rig ht ne ss an d th e s lo pe o f s pe ct ra in th e V N IR re gi on ar e m or e u se fu l f or id en tif yi ng g la ci er fa ci es . 1 2 3 4 5 6 7 8 9 10 W et s em i-s at ur at ed Ic e S ha do w ed s no w M ix o f 4 a nd 9 /1 0 M el tin g Ic e ‘D ry ’ I ce S no w S no w S ha do w O ff- gl ac ie r O ff- gl ac ie r 1 2 3 4 5 6 7 8 9 10 Fr es h sn ow M el tin g Ic e O ff- gl ac ie r O ff- gl ac ie r W et a sh y ic e A sh y sn ow A sh y sn ow A sh y sn ow (p it bo tto m ) A sh y sn ow Fr es h sn ow a) b) c) d) C la ss I nt er pr et at io n C la ss I nt er pr et at io n 51 Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 5.2 Qualitative Spectral Analysis Through a qualitative comparison of ETM+ remote classifications and in situ spectra, we are able to identify not only what facies are present in the satellite imagery, but more importantly what information leads up to those conclusions; the qualitative analysis gives some motivation for the process of quantitative analysis. We do note that almost a decade has passed between image capture and spectra measurement. It is acknowledged that local conditions and processes do vary, and this will be borne in mind. Landsat imagery and spectra were collected at similar days of the year, and at both sites meteorological records do not indicate any anomalous melt or accumulation for years with Landsat or in situ measurements. In addition, for Langjökull, the 2000 Hekla eruption can be expected to have an influence similar to the 2011 Grímsvötn eruption. Nevertheless, as mentioned earlier, because we are comparing spectral measurements and not their spatial location, and because the same physical properties influence the glacier every year, these comparisons should remain valid. Starting with the Svalbard glaciers (see Figure 5.1a), the exclusive presence of class 8 proximal to known ridges leads us to call it shadow, and similarly the exclusively distal location of classes 9 and 10 (confirmed from field experience) indicates that they are rocky/off glacier; class 3 appears to be a mix of classes 9/10 and class 4. These examples show very well that bands in the SWIR (1450-2500 nm) can give information about on/off glacier, but do not distinguish on-glacier classes in this study of end-of-melt-season Arctic glaciers; other exceptions have been referenced in the literature, including classification of on-glacier debris lithologies (Casey et al. 2012) or classification of cold snow based on grain size algorithms (e.g. Bourdelles & Fily 1993). Both the geographic position and brightness of class 6 correspond very well with coarser snow. Classes 1 and 2 are have small numbers of pixels; due to its low perceived reflectance, class 2 appears to be mostly shadowed snow, but class 1 has another explanation. Class 1 matches very well with the measured spectrum of “wet semi-saturated ice” and has low reflectance but is associated with areas of surrounding higher or lower reflectance. We therefore conclude that class 1 indicates either crevassing in the accumulation area (as seen by Sidjak & Wheate 1999) or areas of lower water saturation in the ablation zone. Due to their lower reflectances, classes 4 and 7 match up well with wetter icy surfaces, which depending on location could indicate higher crevassing or melt at lower elevations. Class 5 is the only remaining undiscussed class and falls in between the bright snow and darker ice; thanks to correspondence in brightness with the “dry ice” we are able to identify that it is an ablation facies rather than an accumulation facies – a conclusion which in this case is not possible simply from visual interpretation of the imagery. It is worth noting that while good agreement was found between in situ spectra and satellite classes in the visible, the simple atmospheric correction and reflectance calculation does not appear as robust in the NIR (700-1400 nm). Nevertheless, comparison with in situ data allowed for identification and confirmation of classes which would not have been possible in their absence. As with the in situ spectra, the Langjökull qualitative comparison is again more complex than the Svalbard equivalent. While there is some disparity between the absolute reflectance of the in situ spectra and the satellite classes, the relative behaviour of classes is still consistent. Again, the on-glacier classes have reflectances which cross each other rather than remain in relative order (see Figure 5.1b). Appearing only in a small western area of the icecap, classes 1 and 10 are the brightest and smallest classes Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 52 – from this we deduce a localized and transient snow precipitation event. Class 2, on the other hand, is widely distributed across the entire periphery of the icecap and therefore is very likely melting ice. Due to their presence at the periphery of the icecap and inspection of the original imagery, classes 3 and 4 are concluded to be off-glacier (confirmed from field experience) and again bring us to the conclusion that SWIR is useful for distinguishing glaciated versus non-glaciated areas, but the VNIR is more important for distinguishing on-glacier surface types. Therefore, for the purposes of this study, contrary to Sidjak and Wheate (1999) and Braun et al. (2007), glacier classifications should emphasize the VNIR rather than attempt to avoid those wavelengths, despite possible saturation. The relative behaviour of classes 5 and 8 exhibit very similar behaviour to ‘wet ashy ice 2’ and ‘ashy snow pit bottom’ – with one class brighter in the visible, but dipping slightly underneath the other in the NIR. Based upon this comparison, we determine that class 5 is melting, ashy ice while class 8 is highly concentrated ash on a snow surface. Similarly, the relative behaviour of classes 6, 7, and 9 correspond to different types of ashy and metamorphosed snow surfaces. It appears that the ash (or other particulate matter) is largely what is causing the “crossing” behaviour within the VNIR. These qualitative comparisons serve to highlight that beyond absolute brightness within the VNIR, the slope and concavity within these wavelengths are important for distinguishing surface classes. To this end, it is likely that sensor-specific indices will prove to be the most useful in obtaining the most information from multispectral data, a similar conclusion to that reached by Casey et al. (2012). Figure 5.2: Principal component analysis of in situ spectra from Midtre Lovénbreen and Langjökull. Note the similarity of PCs between field sites. 53 Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 5.3 Full-Spectrum Principal Component Analysis To look at Landsat-appropriate indices (i.e. band combinations tailored to the spectral and radiometric properties of the Landsat TM and ETM+ sensor), and potentially those of other sensors, it is necessary to move from qualitative comparison to a quantitative analysis. As explained in Section 4.3, PCA allows coherent use of a highly multidimensional dataset. PCA of the full set of in situ spectra is taken for each field site; the coefficients of this PCA show the wavelengths that are the keys to clustering and identifying the different surface classes. Figure 5.2 shows the first, second and third PCs for the spectra from Midtre Lovénbreen and Langjökull, using a window of 350 to 1350 nm; because the qualitative comparison showed the importance of variation within the VNIR, only these wavelengths are used, so as to avoid “noise” in the result introduced by other wavelength ranges. In considering the PCs, because each PC is orthogonal, relative sign within PCs is meaningful, whereas across PCs it is not. PC1 is very similar between both field sites (representing 97.4% of the variance in the in situ spectra from Langjökull and 97.8% from Midtre Lovénbreen) – maintaining a constant sign without very much variation. In agreement with previous studies (Greuell et al. 2002; Knap et al. 1999b), this unsurprisingly indicates that the broadband albedo is the most important distinguishing factor between spectra. PC2 represents 2.2% of the variance in the Langjökull spectra and 1.6 % for Midtre Lovénbreen; also similar between Iceland and Svalbard, PC2 shows that the coefficients go from positive to negative, crossing at ~750 nm. This sign change indicates that of further importance is the difference in reflectance between the blue wavelengths (450-500 nm) and the NIR – in other words, the slope that was identified in qualitative comparisons of how quickly the reflectance drops in the higher wavelengths. In addition, the similarity of PC1 and PC2 between Iceland and Svalbard (r2 = 0.83 and 0.97 respectively) indicates that these conclusions should be largely transferrable across many glaciated areas; very similar principal components were also achieved when analysing Landsat imagery over glaciers by Sidjack and Wheate (1999) and the ETM+ images in this study as well, thus reinforcing this conclusion that the PCs are robust across glaciers and across time. PC3 contains only 0.4% of the variance for the spectra and also differs more significantly between Iceland and Svalbard (r2 = 0.31). Nevertheless, both sites show to varying degrees two sign changes which point to the importance of the relative brightness of mid-VNIR wavelengths (~650-1000 nm) compared to the longer and shorter wavelengths. 5.4 Spectrum Clustering Using Principal Component Analysis PCs were calculated using the suites of coefficients for the full spectra (see Figure 5.2 for coefficients). For Midtre Lovénbreen (see Figures 5.3a, b), while individual surface types are not fully distinguishable, it is very clear that wet classes, dry icy classes, and snow classes are quite distinct. While clustering algorithms were able to distinguish these 3 distinct groupings, there was no consistent or reliable demarcation within those classes. Thus, although more spectral information is available, we are not able to move beyond the results of Williams et al. (1991). We do not have data to test whether firn is successfully detected as in other studies (Braun et al. 2007; de Woul et al. 2006). A similar figure was presented by Zeng et al. (1983), but it was restricted to only subsets of the presented spectra and the grouping scheme used was not robustly repeatable. We note that through this grouping, as expected Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 54 Fi gu re 5 .3 : C lu ste rin g bi pl ot s o f i n sit u sp ec tr a fro m M id tre L ov én br ee n ca lc ul at ed u sin g fu ll- sp ec tr um in fo rm at io n: (a ) P C1 v s. PC 2 an d (b ) P C1 v s. PC 3 as sh ow n in F ig ur e 5. 2 an d ca lc ul at ed u sin g w ei gh te d av er ag es o f r es tr ic te d sp ec tr a w hi ch m im ic th e b an ds o f t he E TM +, (c ) L C1 v s. LC 2 an d (d ) L C1 v s. LC 3 as sh ow n in E qu at io ns 5 .7 -5 .9 . a) b) c) d) 55 Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies a) b) c) d) Fi gu re 5 .4 : C lu ste rin g bi pl ot s o f i n sit u sp ec tr a f ro m L an gj ök ul l c al cu la te d us in g fu ll- sp ec tr um in fo rm at io n: (a ) P C1 v s. PC 2 an d (b ) P C1 v s. PC 3 as sh ow n in F ig ur e 5 .2 an d ca lc ul at ed us in g w ei gh te d av er ag es o f r es tr ic te d sp ec tr a w hi ch m im ic th e b an ds o f t he E TM +, (c ) L C1 v s. LC 2 an d (d ) L C1 v s. LC 3 as sh ow n in E qu at io ns 5 .7 -5 .9 . Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 56 through the PCA eigenvalues, almost all the information necessary to distinguish classes is contained within the PC1 vs PC2 biplot as demonstrated with a merged k-means clustering of the PC1-PC2 biplot (see Table 5.1a) similar to the method of Landsat classification; PC3 is largely extraneous. The situation for Langjökull is again more complicated (see Figures 5.4a, b). It appears that in the case of a more complex, ashy surface, PC3 becomes more important in distinguishing surface types. While white ice is clearly visible in the PC1-PC2 plot, much like the drier classes in the previous case, the rest of the field is largely indistinguishable. Clustering algorithms confirm that the white ice classes are reliably distinct from the rest of the classes. In the PC1-PC2 plot, some drifted snow classes appear merged with ashy snow (of which there are also many types), and so this plot is not sufficient on its own. However, the biplot of PC1 and PC3 clearly distinguishes new snow from other surface types; this matches well with Hendriks and Pellikka’s (2004) simple threshold to identify fresh snow. Unlike the Svalbard case, PC3 on the Icelandic glacier is more extreme in its difference between the middle and outlying wavelengths, with a trough at lower wavelengths as well. This brings us to the conclusions that this type of third principal component can distinguish surface types of ashy or dust-covered glaciers from clean surface classes, although it is likely that a PC analysis is more complex than necessary for this purpose and that a simpler iterative or multi-component approach may be more appropriate for wider application. 5.5 Building Transferrable Linear Combinations for Classification With these analyses, we are able to identify which parts of the VIR spectrum are most important in distinguishing glacier surface types, however it does not remain clear whether the ETM+ is spectrally sufficient to accomplish this task. Therefore, in situ spectra are converted into ETM+-like spectra as described in Section 4.1, and then a PCA is performed as it was to the full spectra , resulting in the following linear combinations (LCs). Midtre Lovénbreen LC1 = 0.4662 B1 + 0.4512 B2 + 0.5112 B3 + 0.5383 B4 (5.1) Midtre Lovénbreen LC2 = 0.4933 B1 + 0.4200 B2 – 0.0446 B3 – 0.7604 B4 (5.2) Midtre Lovénbreen LC3 = 0.5257 B1 – 0.0741 B2 – 0.7737 B3 + 0.3456 B4 (5.3) Table 5.1: Contingency table of in situ spectra groupings as observed in the field (rows) from Midtre Lovénbreen as compared to (a) clustered classes of a merged k-means clustering of PC1 and PC2, full spectra and (b) clustered classes of a merged k-means clustering of LC1 and LC2, Landsat-band restricted spectra (columns). Values are percentage of total spectra measured during the Svalbard field campaign. a) Snow (in situ): Wet (in situ): Ice (in situ): b) Snow (in situ): Wet (in situ): Ice (in situ): Snow (classified): 41.0% 0% 0.5% Snow (classified): 41.5% 0% 1.2% Wet (classified): 0% 21.3% 0% Wet (classified): 0% 21.3% 1.2% Ice (classified): 1.7% 0% 35.5% Ice (classified): 1.2% 0% 33.6% 57 Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies Langjökull LC1 = 0.5298 B1 + 0.5234 B2 + 0.5072 B3 + 0.4337 B4 (5.4) Langjökull LC2 = 0.5712 B1 + 0.2162 B2 – 0.1557 B3 – 0.7764 B4 (5.5) Langjökull LC3 = 0.5667 B1 – 0.3556 B2 – 0.6001 B3 + 0.4384 B4 (5.6) Coefficients are then considered together and rounded to simple approximations of the PCs to give the following “ETM+ band” LCs. This has two benefits. One, it aids in a conceptual understanding of what each LC is emphasizing within the data. Two, it facilitates wider transferability of LCs by not being so specifically tailored to the field data. That is to say, LC1 is representative of the broadband albedo, LC2 emphasizes the difference between blue/green and red/infrared and LC3 highlights the difference between blue/infrared and green/red. While PCA is scene-dependant, these LCs are considered to be widely applicable because of the similarity of PCA results as discussed above. LC1 = B1 + B2 + B3 + B4 (5.7) LC2 = B1 + B2/2 – B3/2 – B4 (5.8) LC3 = B1 – B2/2 - B3 + B4/2 (5.9) Encouragingly, the clustering results from Midtre Lovénbreen for the restricted spectra are very similar to the full-spectrum calculations (see Figures 5.3c, d). Further confidence in classifications based on the ETM+ LCs is built because we work in a direction opposite of published studies (e.g. Boresjö Bronge & Bronge 1999); rather than interpret a satellite image and make inferences, the technique is based upon the full spectra and therefore radiometric qualities of the various surface classes are not in question. Again, wet classes, drier ice classes, and snowy classes are the three distinguishable features, although the snowy classes are slightly higher in LC2 as a result of higher reflectance at low wavelengths. The results of a contingency table of this clustering compared to field observations (see Table 5.1b) indicate that LC clustering is successful and very similar to the PC clustering. In this field site, virtually all information obtained with the full spectrum is reproduced using only ETM+ wavelengths, strongly indicating that this multispectral sensor is well suited to monitoring glacier surface classes in this region. For Langjökull, the restricted spectra still require all three LCs, but restricted LC2 is unable to resolve the spectra as well as with the full spectrum information (see Figures 5.4c, d). In particular, the drifted snow classes not only appear similar to the ashy snow classes, but also the white ice classes. Sidjak and Wheate (1999) noted that debris distinction was difficult, which they attributed to ill-defined spectral definition (i.e. poor spectral resolution). We have very well defined spectral information for debris and ash, although limited lithological knowledge (unlike Casey et al. 2012). While the LC1-LC3 biplot is still able to identify the new snow classes as well as the largely debris-ridden areas, which would be impossible given just LC1, it is not really possible to say that clean, icy classes are distinguishable. Therefore, in the case of the ashy and dusty glacier, the ETM+ wavelengths are insufficient compared to the full-spectrum surface reflectance values. Although this study started with a classification including PCA, NDSI, and band ratios after de Woul et al. (2006) and Braun et al. (2007), it now presents the data within a two dimensional feature space (or three, in some cases, as opposed to Sidjak and Wheate’s (1999) four PCs); LCs can then be easily Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies 58 clustered using a variety of algorithms. The LC biplots & clustering contingency tables indicate that for simple glaciers earlier methods are unnecessarily complicated, and the simpler clustering presented here provides an acceptable solution. The aims of glacier surface classifications include application to mass balance measurement, hydrological understanding, and assimilation into energy balance models. Because classes were chosen based upon the different surfaces observed in the field rather than an idealized structure, it is not immediately obvious what their use may be. Nevertheless, this study fits into the framework provided by the literature, providing classes similar to those of Williams et al. (1991) and particularly yielding classes which can be merged and interpreted as accumulation and ablation facies for application to mass balance studies following, for example, Rabatel et al. (2005). The range of surface classes is repeatable, can be characterized, and would be able to be assimilated into an energy balance scheme, too. 5.6 Section Summary In this study, full-spectrum in situ surface reflectance data from a range of end-of-melt-season glacier surface zones are combined with Landsat ETM+ imagery to inform, improve, and better understand satellite classifications of glaciers surface facies. By starting through investigation and analysis of these spectra, this study arrives at an incrementally more in-depth understanding of the limitations of Landsat ETM+ data to glacier surface classification. • Qualitative comparison of in situ data and unsupervised classifications lead to the conclusion that while VNIR compared to SWIR is important for identifying glacier versus non-glacier, contrary to published glacier facies classification techniques the VNIR wavelengths contain the most important information for distinguishing these glacier surface types. Not only is absolute brightness / broadband albedo important, but particular indices are crucial for this task. • PCA of the full-spectrum data shows that only two principal components are required to segment a “clean” glacier surface, such as Svalbard’s Midtre Lovénbreen, while more thought and a third PC are necessary for this task on ash-covered Langjökull. • Similarity of PCs from Iceland and Svalbard, as well as from separate ETM+ imagery, indicates that these are robust and transferrable across both time and field locations. • Full-spectrum data were used to simulate data obtained using Landsat’s ETM+, leading to a new group of band combinations. For Midtre Lovénbreen, the clustering with the restricted wavelengths performed as well as the full spectra in clustering surface facies, while this was not the case for Langjökull’s more complex surface. 59 Chapter 5: Using In Situ Spectra to Explore Landsat Classification of Glacier Facies • Just two or three transferrable linear combinations yields classes which can be reconciled with the traditional facies structure. For “clean” glaciers, while LC1 contains the primary information for delineating glacier facies, it is only through combination with LC2 that delineations can begin to be successfully drawn. A glacier with ash or other significant debris requires not only the inclusion of further LCs, but also more spectral information than a sensor like the ETM+ is able to provide, particularly within the VNIR. • The linear combinations presented here are transferrable across glaciers and time periods. Output classes are suited for a range of glaciological applications. It is confirmed that Landsat’s ETM+ has sufficient spectral resolution to effectively classify glacier surfaces as well as full spectra. Indeed, hyperspectral data do not give significantly more insight into classifying many glacier surfaces, but is crucial for ash and debris. Moving forwards, there are open questions about whether other multispectral sensors behave in a similar fashion. For further sensor- specific indices and LCs, the impact that pixel/footprint size may have on distinguishing surface classes and the role that radiometric accuracy may play (i.e. 8-bit ETM+ versus the 12-bit OLI) should be tested. 60 61 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 6. Using Field Spectroscopy and ATM Imagery to Assess the Role of the Spatial, Spectral, and Radiometric Properties of Multispectral Imagers in Glacier Surface Classification Multispectral sensors are suited for particular applications directly related to their temporal, spatial, spectral, and radiometric resolutions and ranges (see Section 2.5). While temporal resolution and length of archive determine what land cover changes are detectable and at what timescales, it is possible to use common data sets as experimental controls to understand the impact of the spatial, spectral, and radiometric properties of multispectral imagers on a given task. To quickly define a few terms as they are used here: spatial resolution is the pixel size; spectral resolution is the number of available sensor bands; temporal resolution is the revisit time to a particular site; radiometric resolution is the bit-rate of the recorded data; spectral range is the spread of available sensor bands; and radiometric range is defined by the highest and lowest radiance values that a sensor is configured to record. As with any remote sensing investigation, for glacier surface classification it is important to pick the best tool for the task. This chapter will build upon classification of linear combinations (LCs) of Landsat imagery in the previous chapter by providing similarly optimised LCs for a range of multispectral sensors. Through the use of in situ spectra (see Section 3.5) to investigate spectral and radiometric properties and ATM imagery (see Section 3.4) to investigate spatial resolution, this chapter aims to answer the question of what qualities define the best multispectral imagers for glacier surface classification. These experiments will begin with spectral and radiometric considerations, transition to an investigation of the impact of spatial resolution, and then combine the two to understand the advantages and limitations of a range of popular multispectral sensors to glacier surface classification. With these goals in mind, sensors were chosen to span a wide range of properties (spatial scales, spectral bands, and gain settings) and also chosen based on wide use and easy data access; in this chapter, the following are considered specifically: • Airborne Thematic Mapper (ATM) • Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) • Landsat’s Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) • Landsat Multispectral Scanner (MSS) • Moderate Resolution Imaging Spectroradiometer (MODIS), Bands 1-7 • Operational Land Imager (OLI) • Moderate Resolution Imaging Spectroradiometer (MODIS), Bands 8+ • Sentinel-2 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 62 See Table 2.1 and Figure 2.5 for details on the spectral and spatial properties of these sensors. Although many imagers could have been included, nearly all will find some analogue in those considered here. 6.1 Impact of Spectral and Radiometric Properties on Glacier Surface Classification The in situ field spectra from Midtre Lovénbreen and Langjökull were scaled and rounded as described in Sections 4.1 and 4.2. Sensor-specific linear band combinations were then used to reduce data dimensionality. These data were subsequently analysed and classified to elucidate the impact of imagers’ spectral and radiometric properties of the suite of imagers on glacier surface classification. 6.1.1 Linear Combinations In Chapter 5, in situ spectra were converted to satellite bandwidths using Landsat band relative response functions (see Section 4.1), and then spectra were clustered based upon simple linear combinations (LCs) of the band-simulating data. The pre-analysis process for doing this for other sensors is similar – the only difference is that the full-spectrum data were used to emulate different sensor-specific bands. Then, as with the Landsat-like data, a Principal Component Analysis (PCA; see Section 4.3) was run on the sensor-emulating visible and NIR bands. PCs were calculated separately for Langjökull (Lang) and Midtre Lovénbreen (ML) field spectra, as seen below: ML ATM LC1 = 0.3453 B1 + 0.3348 B2 + 0.3449 B3 + 0.3576 B4 + 0.3664 B5 + 0.3851 B6 + 0.3786 B7 + 0.3098 B8 (6.1) ML ATM LC2 = 0.1636 B1 + 0.4240 B2 + 0.3892 B3 + 0.2404 B4 + 0.05536 B5 – 0.2114 B6 –0.4703 B7 – 0.5683 B8 (6.2) ML ATM LC3 = 0.4019 B1 + 0.3771 B2 + 0.0602 B3 – 0.2428 B4 – 0.4462 B5 – 0.4601 B6 – 0.0119 B7 + 0.4719 B8 (6.3) Lang ATM LC1 = 0.3672 B1 + 0.3873 B2 + 0.3838 B3 + 0.3792 B4 + 0.3720 B5 + 0.3573 B6 + 0.3573 B7 + 0.3144 B8 (6.4) Lang ATM LC2 = 0.2323 B1 + 0.3986 B2 + 0.2383 B3 + 0.1149 B4 + 0.0063 B5 – 0.1794 B6 –0.4669 B7 – 0.6833B8 (6.5) Lang ATM LC3 = 0.5791 B1 + 0.3606 B2 – 0.1423 B3 – 0.3510 B4 – 0.3888 B5 – 0.3114 B6 – 0.0020 B7 + 0.3781 B8 (6.6) ML ASTER LC1 = 0.5385 B1 + 0.5751 B2 + 0.6158 B3 (6.7) ML ASTER LC2 = 0.6798 B1 + 0.1354 B2 – 0.7208 B3 (6.8) ML ASTER LC3 = 0.4979 B1 – 0.8068 B2 – 0.3181 B3 (6.9) Lang ASTER LC1 = 0.6102 B1 + 0.5919 B2 + 0.5266 B3 (6.10) Lang ASTER LC2 = 0.6072 B1 0.0775 B2 – 0.7908 B3 (6.11) 63 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification Lang ASTER LC3 = 0.5089 B1 – 0.8023 B2 – 0.3120 B3 (6.12) ML TM / ETM+ LC1 = 0.4662 B1 + 0.4512 B2 + 0.5112 B3 + 0.5383 B4 (6.13) ML TM / ETM+ LC2 = 0.4933 B1 + 0.4200 B2 – 0.0446 B3 – 0.7604 B4 (6.14) ML TM / ETM+ LC3 = 0.5257 B1 – 0.0741 B2 – 0.7737 B3 + 0.3456 B4 (6.15) Lang TM / ETM+ LC1 = 0.5298 B1 + 0.5234 B2 + 0.5072 B3 + 0.4337 B4 (6.16) Lang TM / ETM+ LC2 = 0.5712 B1 + 0.2162 B2 – 0.1557 B3 – 0.7764 B4 (6.17) Lang TM / ETM+ LC3 = 0.5667 B1 – 0.3556 B2 – 0.6001 B3 + 0.4384 B4 (6.18) ML MSS LC1 = 0.4768 B4 + 0.5075 B5 + 0.5454 B6 + 0.4666 B7 (6.19) ML MSS LC2 = 0.6532 B4 + 0.2862 B5 – 0.2920 B6 – 0.6373 B7 (6.20) ML MSS LC3 = 0.4497 B4 – 0.3256 B5 – 0.5909 B6 + 0.5854 B7 (6.21) Lang MSS LC1 = 0.5547 B4 + 0.5269 B5 + 0.5030 B6 + 0.3834 B7 (6.22) Lang MSS LC2 = 0.5269 B4 + 0.2274 B5 – 0.2255 – 0.3834 B7 (6.23) Lang MSS LC3 = 0.6058 B4 – 0.4504 B5 – 0.5039 B6 + 0.4197 B7 (6.24) ML MODIS 1-7 LC1 = 0.4625 B3 + 0.4757 B4 + 0.5049 B1 + 0.5380 B2 + 0.1242 B5 (6.25) ML MODIS 1-7 LC2 = 0.4518 B3 + 0.4106 B4 0.0331 B1 – 0.6953 B2 – 0.3779 B5 (6.26) ML MODIS 1-7 LC3 = 0.0189 B3 + 0.1004 B4 + 0.1060 B1 – 0.4121 B2 + 0.8992 B5 (6.27) Lang MODIS 1-7 LC1 = 0.5281 B3 + 0.5218 B4 + 0.5078 B1 + 0.4286 B2 + 0.0857 B5 (6.28) Lang MODIS 1-7 LC2 = 0.4495 B3 + 0.2094 B4 – 0.0468 B1 – 0.6353 B2 – 0.5902 B5 (6.29) Lang MODIS 1-7 LC3 = 0.6237 B3 – 0.1585 B4 – 0.4887 B1 – 0.1123 B2 + 0.5783 B5 (6.30) ML OLI LC1 = 0.0039 B1 + 0.0037 B2 + 0.0030 B3 + 0.0013 B4 + 0.0011 B5 (6.31) ML OLI LC2 = 0.4245 B1 + 0.4235 B2 + 0.4366 B3 – 0.4623 B4 – 0.4859 B5 (6.32) ML OLI LC3 = 0.3494 B1 + 0.3675 B2 – 0.2996 B3 – 0.0971 B4 + 0.8023B5 (6.33) Lang OLI LC1 = 0.4705 B1 + 0.4685 B2 + 0.4624 B3 + 0.4486 B4 + 0.3726 B5 (6.34) Lang OLI LC2 = 0.3523 B1 + 0.2706 B2 + 0.0836 B3 – 0.1291 B4 – 0.6085 B5 (6.35) Lang OLI LC3 = 0.5065 B1 + 0.2264 B2 – 0.2422 B3 – 0.4564 B4 + 0.1958 B5 (6.36) ML MODIS 8+ LC1 = 0.0422 B8 + 0.0416 B9 + 0.0409 B10 + 0.0406 B11 + 0.0404 B12 + 0.0392 B13 + 0.0391 B14 + 0.0424 B15 + 0.0416 B16 (6.37) ML MODIS 8+ LC2 = 0.2702 B8 + 0.2675 B9 + 0.2666 B10 + 0.2702 B11 + 0.2728 B12 – 0.2953 B13 – 0.3200 B14 – 0.3100 B15 – 0.2915 B16 (6.38) ML MODIS 8+ LC3 = 0.2944 B8 + 0.3144 B9 +0.3283 B10 – 0.3246 B11 – 0.3139 B12 – 0.0545 B13 – 0.0296 B14 – 0.1579 B15 – 0.3164 B16 (6.39) Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 64 Lang MODIS 8+ LC1 = 0.3213 B8 + 0.3204 B9 + 0.3191 B10 + 0.3172 B11 + 0.3163 B12 + 0.3043 B13 + 0.3029 B14 + 0.2919 B15 + 0.2546 B16 (6.40) Lang MODIS 8+ LC2 = 0.3648 B8 + 0.3209 B9 + 0.2581 B10 + 0.1871 B11 + 0.1597 B12 – 0.0204 B13 – 0.0322 B14 – 0.1479 B15 – 0.3125 B16 (6.41) Lang MODIS 8+ LC3 = 0.5098 B8 + 0.2989 B9 – 0.0632 B10 – 0.1396 B11 – 0.2056 B12 – 0.3635 B13 – 0.3594 B14 + 0.2608 B15 + 0.0083 B16 (6.42) ML Sentinel-2 LC1 = 0.2830 B1 + 0.3026 B2 + 0.3032 B3 + 0.3209 B4 + 0.3056 B5 + 0.3241 B6 + 0.3353 B7 + 0.3572 B8 + 0.3294 B9 + 0.2940 B10 (6.43) ML Sentinel-2 LC2 = 0.4209 B1 + 0.4679 B2 + 0.4237 B3 + 0.1603 B4 – 0.0027 B5 – 0.1137 B6 – 0.2447 B7 – 0.3152 B8 – 0.3110 B9 – 0.3602 B10 (6.44) ML Sentinel-2 LC3 = 0.4018 B1 + 0.2449 B2 – 0.0275 B3 – 0.4491 B4 – 0.4750 B5 – 0.2860 B6 – 0.0728 B7 + 0.0854 B8 + 0.1710 B9 + 0.4764 B10 (6.45) Lang Sentinel-2 LC1 = 0.3444 B1 + 0.3676 B2 + 0.3543 B3 + 0.3417 B4 + 0.3096 B5 + 0.3037 B6 + 0.2960 B7 + 0.3053 B8 + 0.2764 B9 + 0.2421 B10 (6.46) Lang Sentinel-2 LC2 = 0.4858 B1 + 0.4221 B2 + 0.2521 B3 + 0.0559 B4 – 0.0374 B5 – 0.1166 B6 – 0.2496 B7 – 0.3451 B8 – 0.3450 B9 – 0.4514 B10 (6.47) Lang Sentinel-2 LC3 = 0.5253 B1 + 0.1853 B2 – 0.2591 B3 – 0.4227 B4 – 0.3507 B5 – 0.2448 B6 – 0.0454 B7 + 0.0916 B8 + 0.1429 B9 + 0.4797 B10 (6.48) Coefficients were rounded and compared. As with full-spectrum PCs (see Figure 5.2) and Landsat LCs (see Section 5.5), the coefficients are very similar between field sites. A unified set of LCs was produced to simplify further analysis and make it transferrable across glaciers and time periods. Again, as with full-spectrum PCs and Landsat LCs, LC1 in each case is representative of VNIR albedo, LC2 emphasizes the difference in reflectance at blue / green wavelengths and red / NIR wavelengths, and LC3 highlights the difference in blue / NIR reflectance and green / red reflectance: ATM LC1 = B1 + B2 + B3 + B4 + B5 + B6 + B7 + B8 (6.49) ATM LC2 = B1/2 + B2 + B3 + B4/2 – B6/2 –B7 – B8 (6.50) ATM LC3 = B1 + B2 – B4 – B5 – B6 (6.51) ASTER LC1 = B1 + B2 + B3 (6.52) ASTER LC2 = B1 – B3 (6.53) ASTER LC3 = B1/2 – B2 – B3/2 (6.54) TM / ETM+ LC1 = B1 + B2 + B3 + B4 (6.55) TM / ETM+ LC2 = B1 + B2/2 – B3/2 – B4 (6.56) TM / ETM+ LC3 = B1 – B2/2 – B3 + B4/2 (6.57) 65 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification MSS LC1 = B4 + B5 + B6 + B7 (6.58) MSS LC2 = B4 + B5/2 – B6/2 – B7 (6.59) MSS LC3 = B4 – B5 – B6 + B7 (6.60) MODIS 1-7 LC1 = B3 + B4 + B1 + B2 + B5/2 (6.61) MODIS 1-7 LC2 = B3 + B4 – B2 – B5 (6.62) MODIS 1-7 LC3 = B3 – B4 – B1 – B2 + B5 (6.63) OLI LC1 = B1 + B2 + B3 + B4 + B5 (6.64) OLI LC2 = B1 + B2 – B4/2 – B5 (6.65) OLI LC3 = B1 + B2/2 – B3 – B4 + B5 (6.66) MODIS 8+ LC1 = B8 + B9 + B10 + B11 + B12 + B13 + B14 + B15 + B16 (6.67) MODIS 8+ LC2 = B8 + B9 + B10 + B11 + B12 – B13 – B14 – B15 – B16 (6.68) MODIS 8+ LC3 = B8 + B9 – B11 – B12 – B13 – B14 + B15 (6.69) Sentinel-2 LC1 = B1 + B2 + B3 + B4 + B5 + B6 + B7 + B8 + B9 + B10 (6.70) Sentinel-2 LC2 = B1 + B2 + B3 + B4/2 – B6/2 – B7/2 – B8 – B9 – B10 (6.71) Sentinel-2 LC3 = B1 + B2 – B4 – B5 – B6/2 + B8/2 + B9/2 + B10 (6.72) 6.1.2 Midtre Lovénbreen Clustering and Classification Analysis For this study, Midtre Lovénbreen provides an example of a largely “clean” and idealised glacier. LCs 1, 2, and 3 for each sensor (see Equations 6.49 to 6.72) and setting were input into a k-means classification (see Section 4.4); all cases are described in Table 6.1 and presented in Figure 6.1. Following earlier studies (De Woul et al. 2006; Braun et al. 2007), a larger number of clusters (here, seven) was defined and then manually merged by the user into three main classes. As presented in Section 5.4, these are defined as snow surfaces, ice surfaces, and wet surfaces; these three main classes are circled in Figure 6.1a. Classification success is assessed with known surface types observed during fieldwork, quantified with the Cohen’s Kappa (K, see Section 4.5), and also described with terminology from Monserud and Leemans (1992). The results presented in Table 6.1 are ranked in order of descending K. Overall, all sensors in all settings are largely able to classify the three clusters. Admittedly, this is the idealised case because there are no mixed pixels to contend with, but even the worst of the cases attains almost 80% success and is considered “very good.” The only setup which does not appear to group the wet classes as well as the other sensors is MODIS, using bands 8 and higher, in a low gain setting (Figure 6.1o); this is possibly due to truncation of contributions to the LCs at Band 16, thereby decreasing information from NIR wavelengths. Radiometric resolution on its own does not appear to be important in the “clean glacier” case on Midtre Lovénbreen, as can be seen by comparing Figures 6.1a, 6.1d, 6.1e, and 6.1c which are all clustered using full spectra but at unrestricted, 16 bit, 12 bit, and 8 bit resolution, respectively. However, when Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 66                                                                a) Figure 6.1: Biplots of the first and second linear combinations of in situ spectral data from Midtre Lovénbreen modified to mimic the spectral and radiometric properties of a range of multispectral imagers. See Table 6.1 for the details of all individual plots, (a) through (p). The ellipses in Figure 6.1a also show the three groups into which the data were classified.            b) c) d) e) f) ASTER, hi ASD Fieldspec ASD Fieldspec, 8 bit ASD Fieldspec, 12 bit ASD Fieldspec, 16 bit Landsat MSS 67 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification                                                                                                               g) i) k) m) o) h) j) l) n) p) ASTER, low 1 ASTER, normal TM / ETM+, LLLH MODIS 1-7 MODIS 8+, hi TM / ETM+, LLLL MODIS 8+, lo TM / ETM+, HHHH Sentinel-2 ATM Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 68 Sensor ASD FieldSpec ASD FieldSpec ASD FieldSpec ASD FieldSpec Sentinel-2 Landsat TM / ETM+ ASTER Landsat MSS ASTER MODIS 1-7 Landsat TM / ETM+ ATM Landsat TM / ETM+ ASTER MODIS 8+ MODIS 8+ Gain Settings - - - - - LLLH low 1 - normal - HHHH - LLLL hi hi lo Radiometric Resolution - 16 bit 12 bit 8 bit 12 bit 8 bit 8 bit 8 bit 8 bit 12 bit 8 bit 16 bit 8 bit 8 bit 12 bit 12 bit K 1 1 1 1 1 1 1 1 0.9081 0.6579 0.6559 0.5979 0.5841 0.5101 0.4879 0.3522 Monserud Agreement Perfect Perfect Perfect Perfect Perfect Perfect Perfect Perfect Excellent Good Good Good Good Fair Fair Poor Figure 6.2 a b c d e f g h i j k l m n o p Table 6.2: For Langjökull spectra: sensors, radiometric resolution, gain settings, K (Cohen’s Kappa) of clustering success, description of classification success (Monserud & Leemans 1992), and corresponding panel in Figure 6.2. Sensor ASTER ASD FieldSpec ASD FieldSpec ASD FieldSpec Landsat MSS ASTER ASD FieldSpec ATM ASTER Sentinel-2 Landsat TM / ETM+ Landsat TM / ETM+ MODIS 1-7 MODIS 8+ MODIS 8+ Landsat TM / ETM+ Gain Settings hi - - - - low 1 - - normal - LLLH HHHH - hi lo LLLL Radiometric Resolution 8 bit 8 bit 16 bit 12 bit 8 bit 8 bit - 16 bit 8 bit 12 bit 8 bit 8 bit 12 bit 12 bit 12 bit 8 bit K 0.9816 0.9742 0.9742 0.9742 0.9705 0.9669 0.9669 0.9631 0.9595 0.9484 0.9448 0.9267 0.9158 0.8938 0.8721 0.7977 Monserud Agreement Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Excellent Very Good Figure 6.1 b c d e f g a h i j k l m n o p Table 6.1: For Midtre Lovénbreen spectra: sensors, radiometric resolution, gain settings, K (Cohen’s Kappa) of clustering success, description of classification success (Monserud & Leemans 1992), and corresponding panel in Figure 6.1. 69 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification combined with a limited set of bands to use in LCs, the quantisation does begin to appear (e.g. Equations 6.52 and 6.53 used to produce Figures 6.1g and 6.1i, emulating ASTER). This example supports the importance of using LCs with a higher number of band combinations which better represent the full spectrum surface reflectance. Radiometer properties (i.e. radiance range and gain settings) do appear to play some role in classification accuracy, but it is hard to assess fully with such widespread success for the Midtre Lovénbreen spectra. Saturation is clearly visible in some Landsat settings; its distinctive signature is a linear alignment of spectra with high LC1 values for Landsat TM/ ETM+ HHHH and Landsat MSS in Figures 6.1l and 6.1f, respectively. Landsat TM/ETM+ LLLH replicates original spectra biplots better than Landsat LLLL (Figures 6.1k and 6.1p, respectively), which makes sense given the high to lower reflectance transition of glacier surfaces from visible to NIR wavelengths. Interestingly, ASTER in a high gain setting (Figure 6.1b) appears to have such success in classification because all of the bright, snowy surfaces were compressed into the same saturated point in the plot, thereby making the entire class almost entirely homogeneous. Ultimately, while varying slightly in performance, imager radiometric properties do not provide significant limitations or guidance in selecting the most appropriate multispectral imager for surface classification of clean glaciers. Therefore, more consideration should be given instead to the desired spatial resolution, revisit time, and data availability. 6.1.3 Langjökull Clustering and Classification Analysis Owing to the unique environmental setting of Iceland, Langjökull provides an example of a glacier with a more complex set and larger range of surface classes, furnished in this case by ash from the Grímsvötn eruption in spring 2011. For Langjökull spectra, a very similar analysis to Midtre Lovénbreen (in the preceding section) is performed. For all spectra, LCs 1, 2, and 3 for each sensor (see Equations 6.49 to 6.72) and setting were input into a k-means classification (see Section 4.4); all cases are described in Table 6.2 and presented in Figure 6.2. However, the results were merged into only two classes (clean ice and other; this is indicated by the circle in Figure 6.2a), and therefore started from six classes rather than the seven used for Midtre Lovénbreen. In addition, PC2 was scaled by a factor of 10 to aid classification success. Again, classification is assessed using knowledge from fieldwork and Cohen’s Kappa, which is described with terminology from Monserud and Leemans (1992). The results are presented in Table 6.2 ranked in order of descending K. While the “clean” glacier ice was easily classified across all sensors, and despite the apparently simpler task of dividing into two groups rather than Midtre Lovénbreen’s three, there is a much larger range of success between sensors and settings for the Langjökull spectra. For the ash-covered glacier, no imager emulation is fully able to represent the range of information contained in the full spectrum data. From the spread of data points, it appears this is largely due to lack of sensitivity in LC2, although saturation in LC1 also plays a role. Nevertheless, because of the simpler task (i.e. identifying clean ice), it is possible to achieve “perfect” classification success. Radiometric resolution, again, is not found to be important on its own, but the quantising effects are again seen in ASTER because of the smaller number of bands contributing to the LCs (see Figures 6.2g, Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 70                                                                                                                                                             Figure 6.2: Biplots of the first and second linear combinations of in situ spectral data from Langjökull modified to mimic the spectral and radiometric properties of a range of multispectral imagers. See Table 6.2 for the details of all individual plots, (a) through (p). The ellipse in Figure 6.2a highlights the clean ice spectra which were classified as distinct from the rest. a) b) c) d) e) f) ASD Fieldspec ASD Fieldspec, 16 bit ASD Fieldspec, 12 bit ASD Fieldspec, 8 bit Sentinel-2 TM / ETM+, LLLH 71 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification                                                                                                                g) i) k) m) o) h) j) l) n) p) ASTER, low 1 ASTER, normal TM / ETM+, HHHH ATM Landsat MSS MODIS 1-7 TM / ETM+, LLLL ASTER, hi MODIS 8+, hi MODIS 8+, lo Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 72 6.2i, and to a lesser extent 6.2n). Even without restricting spectral range and rounding for radiometric resolution only (no scaling, 16 bit, 12 bit, and 8 bit in Figure2 6.2a, 6.2b, 6.2c, and 6.2d, respectively), saturation does appear to be present in the data; this is because some spectra were measured at higher than 100% reflectance as a result of a specular component of reflectance (see Section 3.5.2). Indeed, radiometric resolution is overshadowed by other factors. For example, 12 bit MODIS sensors would be expected to perform well given their higher radiometric resolution compared to Landsat’s 8 bits. However, MODIS gain settings are tuned for darker land and ocean surfaces and so are less suitable to the task of glacier surface classification. As this demonstrates, radiometric range is more important than radiometric resolution. The importance of radiometric range and gain settings is also demonstrated by ASTER and Landsat. ASTER hi gain (Figure 6.2n), Landsat TM / ETM+ LLLL (Figure 6.2m), Landsat MSS (Figure 6.2h), and Landsat TM / ETM+ HHHH (Figure 6.2k) all show a linear feature influencing the higher ranges of both LC1 and LC2, the result of “sensor” saturation. This is more pronounced for ASTER than Landsat because more contributing “bands” are saturated, and LC2 shares more coefficients in common with LC1 for ASTER. Landsat MSS is actually very similar to Landsat TM / ETM+ HHHH, but Landsat LLLH performs better than other Landsat setups. For ASTER, as in the Midtre Lovénbreen case, the brightest classes (New Drifted Snow 1, New Drifted Snow 3, and White Ice 2) are compressed to a single point. The LC1-LC2 biplots for MODIS 8+ in both gain settings (Figure 6.2o and p) demonstrate an intriguing chevron shape. Various theories were considered for why this would occur, including lack of band representation in the NIR or perhaps particular placement of bands in higher and lower wavelengths causing anomalous effects in LC2. In order to find the real culprit, it is necessary to return to exactly what LC1 and LC2 are: MODIS 8+ LC1 = B8 + B9 + B10 + B11 + B12 + B13 + B14 + B15 + B16 (6.67) MODIS 8+ LC2 = B8 + B9 + B10 + B11 + B12 - B13 - B14 - B15 - B16 (6.68) Unlike for other bands, the magnitude of all coefficients of bands contributing to the LCs for MODIS 8+ are identical, except for a flip in sign for higher bands. Where reflectance in Bands 8 to 12 is much higher than 13 to 16, there is a positive linear pattern, and when the opposite is true there is a negative linear pattern. The linear pattern is more pronounced than it would be otherwise because snow and ice spectra have fairly uniform reflectance across the range of wavelengths observed by Bands 8 to 16. Although PCs are uncorrelated, rounding in the coefficients of the LCs caused an artefact in this case. Ultimately, as can be seen with the K rankings in Table 6.2, while a large number of sensors do a “perfect” job identifying clean ice on Langjökull, others perform very poorly. There is a pronounced division between the two, jumping from K = 0.9081 for ASTER normal gain down to 0.6579 for MODIS 1-7. It should be noted that these results hold only for a simple unsupervised classification; supervised or iterative approaches have the potential to yield more specific classes but would lose transferability and ease of implementation. For surface classification of “dirty” or ash-covered glaciers, these results indicate that sensors with “Good,” “Fair,” or “Poor” K values should be foregone in preference for the many alternative sensors which rank higher in performance. 73 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 6.1.4 Concluding Considerations on Spectral and Radiometric Properties Using full-spectrum in situ reflectance data to emulate the spectral and radiometric properties of a range of imagers and settings is a controlled experiment of sorts, removing uncertainty introduced by unknown or changing surface conditions. For both clean and dirty glacier surfaces, although radiometric resolution is largely insignificant, selecting the sensor / gain settings with the most appropriate radiometric range is very important. For the data presented here, Landsat TM / ETM+ LLLH, Sentinel-2, Landsat MSS, ASTER low 1, and ASTER normal perform the classification tests the best. The radiometric properties of the recently-launched OLI were not available at the time of writing, but by considering analogues for its spectral bands and radiometric resolution, it is possible to envision that it would yield results which are a cross between full-spectrum 12 bit results, Sentinel-2, and Landsat TM / ETM+ and would therefore be quite well suited to glacier surface classification, too. 6.2 Impact of Spatial Resolution on Glacier Surface Classification Beyond imagers’ spectral and radiometric properties, another key consideration in classification success is the spatial resolution of the imagery (see Section 2.6.2). To the author’s knowledge, no previous study has addressed this question. 6.2.1 Experimental Strategy and Considerations ATM imagery of most of Midtre Lovénbreen (see Section 3.4) was used in this experiment. LCs were calculated for the 2 m imagery (Equations 6.48 and 6.49), and the image is masked using a manual outline of the glacier. LCs were degraded using bilinear interpolation to 20, 30, 60 and 250 m pixels analogous to Sentinel-2, Landsat TM / ETM+ / OLI, Landsat MSS, and MODIS Bands 1-2, respectively; 500 m imagery is considered too coarse to resolve smaller mountain glaciers. LCs were then input into an ISODATA classification (10 classes, maximum 10 iterations, 95% convergence); output classes were merged by the user into meaningful glaciological classes (i.e. ablation and accumulation facies, see Sections 4.3 and 5.1). As spatial resolution is varied, the radiometric content of all images remains constant as a control. As alluded to in Section 2.6.2, the impact of resolution will depend on the fractal scale of the surfaces being considered. For this purpose, Midtre Lovénbreen is taken to be a representative glacier surface because there is no reason to believe otherwise. Most classification accuracy and quality assessments suffer from their inability to fulfil some assumptions, namely co-registration and random sampling (e.g. Comber et al. 2012); this experiment uses all pixels to assess accuracy, and because coarser imagery is created from finer imagery, co-registration is a not an issue. By definition, the detail available in the 2 m results will be blurred out (to the lower resolution), but it is unknown what beneficial or detrimental effects this may have on surface classification accuracy and what this will mean for individual sensors which provide imagery at a variety of spatial resolutions. In this case, only pixel-based classification is considered because pixel-based classifications have been Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 74 Class 1 2 3 4 5 6 7 8 9 10 2 m Pixels 124,425 196,744 57,694 67,947 107,323 158,711 205,371 179,783 116,938 72,040 Percent 9.7 15.3 4.5 5.3 8.3 12.3 16.0 13.0 9.1 5.6 Interpretation Shadow / Debris Shadow / Debris Ablation Ablation Ablation Ablation Ablation Ablation Accumulation Accumulation Total Area: Shadow/Debris: Ablation: Accumulation: 5.15 km2 25.0 % 60.4 % 14.7 % Table 6.3: Comparison of the results of the classifications presented in Figure 6.3 by considering relative area of each class as well as their aggregation into glaciologically meaningful groups of accumulation and ablation area. traditionally easier to implement with a range of software tools. Although object-based image analysis has been shown to have benefits for very-high resolution imagery (1 m), at any lower resolutions it does not produce statistically significantly different results from pixel-based classification (Baker et al. 2013). 6.2.2 Independent Resolution Assessment For each resolution, the ISODATA classification outputs 10 classes (see Figure 6.3). Classes 1 and 2 are mixes of shadow and thin debris cover, classes 3 through 8 are interpreted as ablation classes, and classes 9 and 10 are interpreted as accumulation classes. For reasons described in earlier sections concerning the potential presence of frost and the use of radiance rather than reflectance, from investigation of the visible imagery it appears this interpretation (i.e. merging of classes) may slightly overestimate the accumulation area. Nevertheless, this is deemed to be preferable to significant underestimation of the area of accumulation facies. Percentages of each class and the aggregated accumulation and ablation areas are presented in Table 6.3. Most classifications have very similar accumulation and ablation area measurements, with the exception of the 250 m pixel classification, which results in slightly less ablation and more accumulation, although the differences are below 5%. While these figures agree, that does not mean the classification results agree on the pixel level; to understand this, further analysis is necessary. Class 1 2 3 4 5 6 7 8 9 10 20 m Pixels 1,383 1,932 610 653 1,023 1,479 2,307 1,883 1,107 765 Percent 10.5 14.7 4.6 5.0 7.8 11.3 17.6 14.3 8.4 5.8 Interpretation Shadow / Debris Shadow / Debris Ablation Ablation Ablation Ablation Ablation Ablation Accumulation Accumulation Total Area: Shadow/Debris: Ablation: Accumulation: 5.26 km2 25.0 % 60.5 % 14.2 % 75 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 0 1 km a) 2 m c) 30 m e) 250 m b) 20 m d) 60 m Class 1 Class 2 Class 3 Class 4 Class 5 Class 6 Class 7 Class 8 Class 9 Class 10 Figure 6.3: Midtre Lovénbreen surface classification using ATM imagery at 2 m resolution (a) and the same image degraded to 20 m (b), 30 m (c), 60 m (d), and 250 m pixels (e). Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 76 Class 1 2 3 4 5 6 7 8 9 10 30 m Pixels 630 889 287 288 454 624 1,005 895 491 358 Percent 10.6 15.0 4.8 4.9 7.7 10.5 17.0 15.1 8.3 6.0 Interpretation Shadow / Debris Shadow / Debris Ablation Ablation Ablation Ablation Ablation Ablation Accumulation Accumulation Total Area: Shadow/Debris: Ablation: Accumulation: 5.33 km2 25.7 % 60.0 % 14.3 % Class 1 2 3 4 5 6 7 8 9 10 60 m Pixels 167 228 93 87 92 153 240 264 119 100 Percent 10.8 14.8 6.0 5.6 6.0 9.9 15.6 17.1 7.7 6.5 Interpretation Shadow / Debris Shadow / Debris Ablation Ablation Ablation Ablation Ablation Ablation Accumulation Accumulation Total Area: Shadow/Debris: Ablation: Accumulation: 5.55 km2 25.6 % 60.2 % 14.2 % Class 1 2 3 4 5 6 7 8 9 10 250 m Pixels 15 16 7 11 7 6 13 17 11 9 Percent 13.4 14.3 6.3 9.8 6.3 5.4 11.6 15.2 9.8 8.0 Interpretation Shadow / Debris Shadow / Debris Ablation Ablation Ablation Ablation Ablation Ablation Accumulation Accumulation Total Area: Shadow/Debris: Ablation: Accumulation: 7.0 km2 27.7 % 54.5 % 17.9 % Table 6.3 (cont.): Comparison of the results of the classifications presented in Figure 6.3 by considering relative area of each class as well as their aggregation into glaciologically meaningful groups of accumulation and ablation area. 77 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification Table 6.4: Comparison of the results of the classifications presented in Figure 6.3 by down-sampling results with majority filtering to lower resolutions. A is overall agreement and K is Cohen’s Kappa, both presented in Section 4.5; perfect agreement would lead to a value of 1 for both A and K. Starting Resolution 2 m 2 m 20 m 30 m 2 m 20 m 60 m 30 m 20 m 2 m A 0.6315 0.7217 0.5788 0.6414 0.5284 0.6126 0.2857 0.1923 0.1975 0.2405 K 0.5825 0.6837 0.5226 0.5927 0.4662 0.5605 0.1881 0.0878 0.1023 0.1425 Final Resolution 30 m 20 m 30 m 60 m 60 m 60 m 250 m 250 m 250 m 250 m A, grouped 0.9477 0.9459 0.9324 0.9084 0.9135 0.9090 0.7727 0.7612 0.6857 0.6912 K, grouped 0.9035 0.9017 0.8746 0.8360 0.8359 0.8326 0.5956 0.5681 0.4554 0.4322 A, ranked 1 2 3 6 4 5 7 8 10 9 K, ranked 1 2 3 4 5 6 7 8 9 10 6.2.3 Down-Sampled Resolution Assessment To begin to understand pixel-level agreement of glacier surface classification at various resolutions, a majority filter was used to down-sample high resolution images to lower resolutions (see Table 6.4). This was done using the algorithm built in to ArcGIS and forced to the same grid using nearest-neighbour interpolation in MultiSpec. Because class numbers are not indicative mathematically of their similarity or difference, mathematic convolution would not be meaningful. After filtering, the similarity of high resolution images and their down-sampled counterparts were assessed using direct pixel comparisons with both the overall accuracy agreement (A) and Cohen’s Kappa (see Section 4.5). As would be expected, similarity is the most meaningful between images with similarly sized pixels (see Table 6.4). However, the 2 and 30 m images and 2 and 20 m images are more similar than the 20 and 30 m images. While this could be the result of a resampling artefact, it does indicate that medium resolution imagery is doing a good job at reproducing the results obtained with high resolution imagery. The 2 m classification results are approximately as similar to the 60 m results as all of the medium resolution classifications are to each other, potentially indicating that 20 or 30 m imagery is more suited to glacier surface classification than 60 m resolution imagery. At the bottom, the low resolution imagery (250 m pixels) by a large step shows the lowest agreement with all other results; according to this result, MODIS imagery is not appropriate for glacier surface classification. The A and K rankings of the classification comparisons differ slightly in the relative position of the comparison of 60 m and higher resolution results, 250 m to 2 m and 20 m images. The inferences drawn above are consistently supported by both A and K, but divisions are more visible in the K values than in the A values. This lends some confidence to the conclusions, because K values should contain more signal and less false agreement than A values. Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 78 Starting Resolution 30 m 20 m 30 m 60 m 60 m 60 m 250 m 250 m 250 m 250 m A 0.7381 0.6393 0.5976 0.6383 0.5958 0.4977 0.2612 0.2342 0.2165 0.2145 K 0.7030 0.5918 0.5443 0.5902 0.5422 0.4317 0.1690 0.1400 0.1207 0.1190 Final Resolution 20 m 2 m 3 m 30 m 20 m 2 m 60 m 30 m 20 m 2 m A, grouped 0.9603 0.9585 0.9450 0.9338 0.9252 0.9065 0.6909 0.6833 0.6711 0.6647 K, grouped 0.9269 0.9233 0.8983 0.8778 0.8609 0.8253 0.4469 0.4223 0.3985 0.3883 A, ranked 1 2 3 4 5 6 7 8 9 10 K, ranked 1 2 3 4 5 6 7 8 9 10 Table 6.5: Comparison of the results of the classifications presented in Figure 6.3 by downscaling results to higher resolutions. A is overall accuracy agreement and K is Cohen’s Kappa, both described in Section 4.5; perfect agreement would lead to a value of 1 for both A and K. 6.2.4 Up-Sampled Resolution Assessment In addition to downscaling high resolution imagery, the low resolution imagery is resampled to high resolution by converting each pixel in the high resolution image to a host of small pixels of the same class (i.e. one 250 m pixel becomes 15,625 corresponding 2 m pixels). This was done using bilinear interpolation in ERDAS Imagine and forced to the same grid using nearest-neighbour interpolation in MultiSpec. Pixels in the original high resolution and the down-sampled classifications are compared and assessed using A and K (see Table 6.5). A and K ranks and relative values are in better agreement than in the previous section. For this round of comparisons, similarity in pixel size was the unambiguous driver of similarity in results. Again, there is a clear break between high (2 m) / medium (20, 30, or 60 m) resolution image results and any comparison to the 250 m resolution classification results. 6.2.5 Concluding Considerations on Spatial Resolution Each method used to compare glacier surface classification at different pixel sizes gives a slightly different impression of the importance of sensor spatial resolution. Figure 6.3 clearly shows the loss of detail associated with observation at lower resolutions, but the relative area of shadow, debris, accumulation and ablation facies is very similar among images of all classes. However, relative accuracy at different spatial scales is dependent on the scale of inhomogeneity within and between classes. For Midtre Lovénbreen, classes appear to be similarly behaved at high and medium resolutions, but the glacier is definitely more complex than 250 m pixels can capture. Although 15 m pixels (analogous to ASTER bands or fused ETM+ images) are not explicitly considered, in view of these results, such an analysis would appear to have been superfluous. Based upon this analysis, for glacier surface classification, high resolution imagery would indeed be desirable. It appears that even the highest resolution that MODIS is capable of providing (250 m) is insufficient for glacier surface classification. Medium resolution imagery is found to be adequate for the task, and 20 or 30 m imagery is preferable to 60 m imagery but not drastically. 79 Chapter 6: Impact of Spatial and Spectral Resolution on Glacier Surface Classification 6.3 Chapter Summary With the increasing availability of a range of multispectral data, it is important not only to know which are the best for glacier surface classification, but also to have transferable techniques to be able to take advantage of as many data sources as possible. For a range of multispectral sensors, the work in this chapter tested the impact of spatial, spectral, and radiometric sensor properties on glacier surface classification. Starting with in situ spectra, linear band combinations for all sensors were calculated for input into classification. It was found that: • The range and resolution of spectral bands did not limit classification results, • Radiometric resolution (i.e. 8 bit, 12 bit, or 16 bit data) was not a significant factor unless combined with a limited number of input bands, and • The most important consideration in choosing a multispectral sensor was radiometric range (i.e. gain settings). The most successful imagers have a low gain in the visible and higher gain in the near infrared. Using an airborne multispectral image, the impact of spatial resolution on classification was analysed using multiple methods. While the 2 m base imagery was taken to be the best representation of the glacier surface, it was also found that: • Course resolution imagery (250 m) is unsuitable for capturing glacier surface classes’ spatial complexity, and • Medium resolution imagery is found to be adequate for the task, but 20 or 30 m imagery is preferable to 60 m imagery. Based upon these findings, particular sensors (and settings) are recommended as well-suited to glacier surface classification. These are: • Landsat TM / ETM+ (LLLH), • Sentinel-2, • ASTER (low 1), • ASTER (normal), and • OLI (by comparison with analogous sensors). Landsat MSS is also found to be radiometrically well-suited for glacier surface classification, but its lower spatial resolution makes it a secondary selection. However, MSS has historical imagery whereas other sensors have more recent (or future) data ranges, and therefore these imagers could be used in conjunction with each other. This demonstrates, once again, that although priority can be given to sensor capabilities, temporal resolution and data availability will still remain important considerations. This chapter presents a healthy range of sensors whose data archives can be harnessed into the future. In the next chapter, this knowledge will be used to try to answer a glaciological question related to glacier surface properties. 80 81 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 7. Glacier Albedo on the Brøgger Peninsula, Svalbard Svalbard’s glaciers (totalling 36,000 km2) have been well studied and are considered to be early indicators of Arctic cryospheric change (e.g. Fleming et al. 1997). As mentioned in Chapters 1 and 2, James et al. (2012) recently showed that glaciers in Svalbard are thinning at an accelerating rate in their higher elevations, but without evidence to explain the process. This thinning could significantly increase the area of net ablation in an area which is highly sensitive to such changes (Hagen et al. 2003), with large consequences for future glacier retreat in the Arctic. One of the proposed mechanisms for higher altitude thinning is a melt-albedo feedback (see Section 2.7), which is best suited to investigation using multispectral satellite remote sensing. Inherent in understanding albedo change are questions of spatial variations in albedo. Previous studies have been challenged by difficulty in upscaling albedo measurements at a point to effectively compare with remote sensing data (e.g. Knap et al. 1999a; Stroeve et al. 2006; Box et al. 2012) because of the large variability of albedo on glacier surfaces, particular in the ablation area. In the context of the spectral work presented in the preceding chapters, this will be immediately obvious. Whole-glacier albedo has been used as a proxy for glacier mass balance (Greuell et al. 2007; Dumont et al. 2012; Kohler et al. 2012), while variability in albedo has been previously shown to be the least in a glacier’s highest and lowest reaches (Greuell et al. 2007). Because the hypothesis raised by James et al. (2012) deals specifically with the snow melt-albedo feedback, a necessary first step will be to classify the glacier into accumulation and ablation facies. The simple and transferrable techniques used earlier in the thesis (see Sections 4.4, 5.5, and 6.2) are well suited to this task, allowing albedo variation to be studied in thematic subsets. Addressing the hypothesis of melt-albedo feedback on Svalbard glaciers also requires investigation of temporal variability. While James et al. (2012) have surface elevation data dating back from 1961, finding albedo estimates back that far will be impossible. However, Landsat imagery provides the best solution available with easy data access dating back to 1976, medium spatial resolution (30 or 60 m pixels, depending on the sensor; see Section 2.5), and consistent and stable radiometric data to within 10% (Markham & Helder 2012). Using the Landsat data archive, this chapter aims to understand how seasonal and interannual albedo of the accumulation and ablation areas can be observed over the past 35+ years and what implications this has for glaciers in Svalbard and the wider Arctic. This chapter will begin by revisiting the processing steps applied to the Landsat imagery including considerations needed in albedo calculation and glacier surface classification, move on to analysis of seasonal variability in albedo, address interannual variability in albedo, discuss factors which inform interpretation of the data, and finally consider conclusions which can or cannot be drawn from the data presented. 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 82 7.1 Data Preparation The research area identified for this study is the Brøgger Peninsula (see Section 3.2), including two glaciers used in the James et al. (2012) study (Midtre Lovénbreen and Austre Brøggerbreen) and one glacier (Austre Lovénbreen) with in situ photographs of glacier surface classes (see Section 3.6, Laffly et al. 2012). Landsat images are listed in Table 3.1. 7.1.1 Atmospheric Correction of Landsat Imagery As described in Section 3.3.2, Landsat images were corrected with LandCor, an implementation of the 6S atmospheric radiative transfer model which includes corrections for atmospheric aerosols, ozone, water vapour, and some solar geometry effects. The LandCor routine (Zelazowski et al. 2011) includes a monthly mean total water vapour and aerosol optical depth as parameterised for each Landsat path and row. For the Svalbard site, all available MODIS data were investigated to create a scene-specific relationship between elevation and total water vapour or aerosol optical depth. However, the atmospheric load of both were found to be an order of magnitude lower than reported for the tropical site used by Zelazowski et al. (2011). In addition, the elevation range of the Svalbard site is much lower than the Zelazowski site. Both of these factors limit the ability to build a robust parameterisation. Therefore, the general parameterisation for total water vapour and aerosol optical depth in LandCor were determined to be sufficient. For the processed Landsat imagery, top-of-atmosphere reflectance and atmospherically-corrected reflectance were compared for all sensors (i.e. MSS, TM, and ETM+). Atmospheric correction was found to have an effect on the order of ~2% reflectance. Given the dry air and low aerosol composition of the atmosphere above Svalbard, such a small correction is not unexpected. 7.1.2 Albedo Calculation Equation 2.3 in Section 2.7.1 presented a method for estimating broadband albedo from Landsat data. Repeated below, it is: α = 0.539rLandsat2 + 0.166rLandsat4(1 + rLandsat4) (2.3) where α is broadband albedo and r is reflectance in the indicated multispectral band. With this NTB conversion, the residual standard deviation is just 0.011 (Greuell et al. 2002). Band 2TM/ETM+ and Band 4TM/ETM+ are valid for TM and ETM+; for MSS, Bands 4MSS and 6MSS were used. The wavelengths are largely coincident, but Band 6MSS edges a bit towards lower wavelengths than Band 4TM/ETM+; this would bias the MSS albedos slightly higher than TM and ETM+. This effect will be borne in mind when interpreting the results. Previous studies (see Section 2.7) have implemented a range of corrections to refine their albedo estimations. One widespread correction is that made for topography and the effect this will have due to variation in incident radiation (Knap et al. 1999a; Klok et al. 2003; Hendriks & Pellikka 2004). Hendriks and Pellikka (2004) estimate an approximately 10% variation with and without this topographic 83 7. Glacier Albedo on the Brøgger Peninsula, Svalbard correction for Hintereisferner, Austria. Dumont et al (2012) conducted a similar experiment, concluding that a topographic correction was more important for a Lambertian assumption than non-Lambertian (this will be considered below). All of these previous studies have attempted to accurately estimate broadband albedo. There is a crucial difference for the study here; this investigation is interested in the change in albedo, not absolute albedo itself. It will be important to keep this in mind when deciding what corrections are necessary to homogenise the data set. Multiple studies have included a topographic correction for albedo (Knap et al. 1999a; Hendriks & Pellikka 2004). When albedo is calculated with the above method, a horizontal surface is assumed. A glacier’s surface is often closer to horizontal than the surrounding area. Nevertheless, without this correction, a north-facing slope will appear darker than in reality (Knap et al. 1999a), for example. Topographic correction is dependent on the solar zenith angle, solar azimuth, surface slope, and the ratio of direct to diffuse radiation, although the latter can be parameterised for glacier and non-glacier surfaces. Despite glacier retreat, except at the very terminus the topography of a glacier’s surface will be largely determined by the basal topography; even there, averaged across an ablation area, the slope will remain largely stable. Thus, it is only changes in solar position that could meaningfully influence this study. Because Landsat is in a sun-synchronous orbit which has been consistent across many missions, this would have minimal impact at low latitudes. High latitudes, however, have larger seasonal variations in solar position, which could cause meaningful differences in topographic correction across scenes. Images used here do not extend into early spring or late autumn, minimizing potential effects. In addition, this would be an intra-annual rather than interannual bias which would be easily circumvented by comparing images taken at a similar time of year. Therefore, a topographic correction is not performed before albedo is calculated. A second, very important correction in estimating albedo from multispectral reflectance is based on the surface bidirectional reflectance distribution function (BRDF, see section 2.7.1), accounting for solar incidence angle and viewing angle from a non-Lambertian (non-uniform) scatterer. Indeed, Dumont et al. (2012) cite the difference in albedo between Lambertian and anisotropic reflectance to be similar to the scale of the topographic corrections above, on the order of 10%. Previous comparisons of field-spectra converted to broadband albedo and albedo estimations using albedo cited surface anisotropy as the dominant reason for persistent underestimation by satellite-derived values; while spectra were collected for that study under diffuse radiation, satellite images were obviously not (Hendriks & Pellikka 2004). However, it is worth noting that although Hendriks and Pellikka attribute the majority of the mismatch, on the scale of 10-20%, to be due to anisotropy, there was a two week window for the glacier surface to metamorphose and darken, which could contribute 0-10% of the darkening in itself (Greuell et al. 2007; Box et al. 2012; Kohler et al. 2012). Nevertheless, it has been well demonstrated than snow and ice surfaces exhibit forward-scattering and nadir darkening (e.g. Choudhury & Chang 1981; Kuhn 1985; Knap & Reijmer 1998; Dumont et al. 2010) and so any BRDF correction will increase calculated albedo values. Greuell and de Ruyter de Wildt (1999) provide a parameterisation for calculating this anisotropy correction for clean ice. Ice’s BRDF is actually fairly well-behaved near nadir; because glacier surfaces are largely low-slope and the change in 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 84 Landsat viewing angle is only 7.7° across the entire scene, a near-nadir simplification provided by Greuell and de Ruyter de Wildt (1999) is appropriate: f = rnb/αnb = 1+ (c4 + c5θz)αnb-m (7.1) where f is a proportionality factor relating narrowband reflectance and narrowband albedo; rnb is narrowband reflectance; αnb is narrowband albedo; c4, c5, and m are empirically derived constants; and θz is the solar zenith angle available in Landsat metadata. The solar azimuth is not needed because the near- nadir BRDF of ice is independent of azimuth. This correction was calculated for a range of reflectances using the range of solar zenith angles present in the Landsat imagery employed in this study (see Table 7.1). Table 7.1 shows that the BRDF correction factor decreases with increasing albedo, reaching values which agree with the ~10% correction cited earlier (Dumont et al. 2012). While there is some variation across scenes, it is small compared to the size of the correction. The impact this would have on the estimated albedo is perhaps best considered for two hypothetical cases. One, across the Landsat record, if the ice got darker since 1976, the change in a BRDF-corrected albedo would then be accentuated compared to a “Lambertian” albedo. Two, conversely, if the ablation area in fact brightened, the change would be diminished in a corrected albedo as compared with a non-corrected albedo. Nevertheless, although the magnitude of the change is different when a BRDF correction is implemented, the change itself is present in both cases. Therefore, again because the goal of this study is to investigate albedo change rather than absolute albedo, and because the difference in correction across all scenes is small, no BRDF correction is implemented. It should be noted that Equation 7.1 presents a parameterisation for ice, not snow. Some field data have been used to build a lookup-table from limited field data for snow reflectance (Dumont et al. 2010). The author was contacted, but it was impossible to obtain these data in a timely fashion. Still, it is expected that the conclusions above for ice will hold for snow, too, as they share similar reflective behaviour in many regards. In addition, when considering more than one surface when applying corrections, it is important not to employ any circular logic. Reflectance corrected with a BRDF should not be used to classify glacier surfaces to subsequently apply a facies-specific BRDF used to re-calculate reflectance (Greuell et al. 2002). This is another reason that a BRDF correction of reflectance and therefore albedo is not employed in this study. There are two days on which the study area was captured by multiple Landsat scenes: 9 July, 1976 and 11 July, 2002. These can be used to help assess the consistency of albedo observations. Although the coincident dates will have very similar solar geometry, they have difference viewing geometry within the scene. For the 1976 scenes, all albedos (averaged across facies or the entire glacier) are measured to be Table 7.1: Mean and standard deviation of BRDF correction factors calculated using Equation 7.1 and parameters for the suite of Landsat images used in this study. The range of reflectances was determined by the range for which empirical constants were determined to be valid. r f 0.3 1.43 ± 0.06 0.2 1.17 ± 0.09 0.5 1.22 ± 0.03 0.1 2.59 ± 0.19 0.4 1.30 ± 0.04 0.6 1.18 ± 0.03 85 7. Glacier Albedo on the Brøgger Peninsula, Svalbard within 2% of each other. In addition, an image was captured on 10 July, 1976 and there was a similar level of precision for individual glacier albedo measurements (although Austre Lovénbreen was not imaged on 10 July, 1976 due to cloud cover). 2002 also shows high levels of agreement between scenes; all albedo averages show variability of about 5%. Encouragingly, all near-coincident albedo values are within the ~10% level of error expected from earlier studies. 7.1.3 Classification Classification of glacier surfaces in conjunction with albedo calculations is an important step in being able to thematically analyse the albedo estimates. This will be especially true for identifying the snow-covered areas (as opposed to ice areas). If the hypothesis that melt-albedo feedback is causing glacier thinning in high elevations (James et al. 2012), then it is essential to isolate the snow facies; in this case, it would be expected to see darkening in accumulation areas and a constant albedo in ablation areas. As detailed earlier (see Sections 4.4, 5.5, and 6.2), the classification is undertaken using a k-means clustering of 3 linear combinations (LCs) of Landsat bands. LCs were the same used in Chapter 5. The k-means classification outputted 10 classes, which were then merged by the user into an interpretation of the broader accumulation area (snow classes) and ablation area (ice classes). Because of changing glacier extent from 1976-2012, as well as changing illumination conditions, masks were manually created for all non-shadowed areas of each glacier in each image. The success of this classification for the Brøgger Peninsula is assessed using time-lapse cameras which observe Austre Lovénbreen (See Section 3.6). The contingency table comparing camera-based and Landsat-based classification (see Table 7.2) confirms some success for the satellite classifications. A values (see Section 4.5) range from 0.7 up to 1.0. These A values are lower than those obtained from spectra-based classifications. However, this real- world classification does have to deal with mixed pixels, thereby increasing ambiguity of classes and (as expected) decreasing the classification accuracy. Figure 7.1 shows that the patterns of accumulation and ablation facies match fairly well, but not perfectly. No large bias in under or over estimation of classes can be seen (in this admittedly small sample); if any, there many be a small tendency to underestimate the accumulation area. Although not immediately ideal for mass balance measurements, this does mean that a) Accumulation (Landsat) Ablation (Landsat) b) Accumulation (Landsat) Ablation (Landsat) c) Accumulation (Landsat) Ablation (Landsat) Accumulation (in situ) 20% 15% Accumulation (in situ) 99% 0% Accumulation (in situ) 6% 10% Ablation (in situ) 15% 51% Ablation (in situ) 1% 0% Ablation (in situ) 2% 82% A: 0.71 A: 1.0 A: 0.88 Table 7.2: Contingency table comparing in situ photography-based and Landsat-based surface classifications for Austre Lovénbreen on 20 August, 2010 (a), 23 June, 2011 (b), and 17 August, 2011 (c). A is the overall classification agreement. 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 86 snow albedo measurements will not be biased towards darkening, which is encouraging for rigorously testing the melt-albedo hypothesis. However, this bias, if present, is small. Thus, the study moves forwards employing this accumulation area vs. ablation area classification scheme. 7.2 Seasonal Variability As explained above, all images are classified and albedo is calculated. As the glaciers evolve throughout each year, a strong seasonal signal should be present both in the TAAR (Transient Accumulation Area Ratio), as well as the albedo. However, over the 35 year satellite record, only 12 years are represented; among those, many years have just one or two images available. Therefore, images from all years are plotted together against the day of the year (DOY) to visualise seasonal change. For the TAAR (see Figure 7.2) there is a strong and consistent decrease from fully snow- covered to almost no accumulation area. It should be noted that these TAARs are the ratio of the area of accumulation facies and the area of the total glacier mask for each image; masks were manually traced to omit shadowed areas. Because shadowed areas are often near steep valley heads above the upper glacier, these TAARs will underestimate the real TAAR of the glacier. Nevertheless, the seasonal transition in TAAR is consistent across years. The increase in TAAR at the end of the summer / beginning of autumn observed in 1988 is not an artefact; meteorological records from Ny-Ålesund (available from the Norwegian Meteorological Institute) indicate that snow fell on 13 September, 1988 (DOY 256) as well as the day preceding image acquisition (DOYs 259, 260). All three glaciers behave similarly; while there Ablation Accumulation a) b) c) d) e) f ) 1 km 20 August, 2010 23 June, 2011 17 August, 2011 In Situ: Ablation Accumulation a) b) c) d) e) f ) 1 km 20 August, 2010 23 June, 2011 17 August, 2011 Landsat: Figure 7.1: Surface classifications of Austre Lovénbreen using oblique photography (a-c) and Landsat imagery (d- f) on 20 August, 2010 (a, d), 23 June, 2011 (b, e), and 17 August, 2011 (c, f). Stripes in the Landsat images are the result of the Landsat 7 SLC error. 87 7. Glacier Albedo on the Brøgger Peninsula, Svalbard                       "    % !#!#                      is some difference between the AARs of each glacier, they demonstrate no consistent relative behaviour. Similar to TAAR, albedo measurements are plotted for entire glaciers (see Figure 7.3a and b), accumulation facies (see Figure 7.3c and d), and ablation facies (see Figure 7.3e and f). Using MODIS data, with a significantly higher temporal resolution, earlier studies documented strong seasonal cycles, with albedo starting with fresh snow values of ~0.8, decreasing rapidly through the melt season to bare- ice values of ~0.5, and subsequently rebounding as melt decreased and accumulation recommenced (Stroeve et al. 2006; Box et al. 2012). In Svalbard in particular, this cycle was observed to take place between May (15 May ≈ DOY 135) and September (15 September ≈ DOY 258; Greuell et al. 2007). The seasonal darkening signal is seen strongly for the full-glacier observations. Based upon the shape of the season effect in earlier studies, it is fit with a 3rd order polynomial (α = -1.5340*10-7 DOY3 +7.8751*10-5 DOY2 – 0.0139DOY + 1.3727, r=0.77, see Figure 7.3a). Although it is a good empirical fit, the low albedo predicted early in the season may be an artefact of the data which limits extrapolation. This is used below (see Section 7.3) to remove the seasonal signal for interannual trend analysis. The seasonal cycle is the result of a combination of effects, including the widespread transition from snow to ice. This darkening signal can be used to track the ELA and is a proxy of mass balance in itself (Greuell et al. 2007; Dumont et al. 2012; Kohler et al. 2012). The brightening observed in early 1976 is likely the result of sensor variability or illumination effects, as there is no new snowfall in the meteorological record to explain it. However, albedos are measured to be consistently lower in the Landsat data than have been reported using MODIS data. Also evident is the apparently anomalously bright signal from 1999; meteorological records indicate that prior to image acquisition, it was ~2 °C and rainy at Ny- Figure 7.2: Estimated Transient Accumulation Area Ratio (TAAR) for Austre Lovénbreen, Midtre Lovénbreen, and Austre Brøggerbreen for years spanning 1976 to 2011 calculated using an unsupervised classification of non- shadowed areas of Landsat images. The strong seasonal transition from snow-covered to ice-covered glacier is visible across all years, as well as onset of autumn in 1988. Error bars on TAAR estimates are on the order of 5-10% (see Section 8.4.1). 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 88 Ålesund, so it is possible that some accumulation fell on the glaciers, in particular Austre Lovénbreen. In addition, seasonal darkening can also be seen for accumulation facies. This is expected because snow darkens as it melts, refreezes, compacts, ages, and metamorphoses. However, much of this evident trend is dominated by the observations from one year in particular, 2010. With this year removed, the darkening trend is somewhat muted. This is possibly because of a cap on the measureable albedo using Landsat (a result of sensor radiometric range), or the fact that most of the observations are concentrated in the middle of the melt season. Spring and autumn are quite underrepresented; although temperatures fell consistently below 0 °C and snow did fall prior to the 1988 autumn measurements, it           # !           !!! # "!                           "    % !#!#            "    % !#!#             "    % !#!#                 a) c) b) d)                      e) f) Figure 7.3: Albedo estimates for the study area displayed seasonally, broken up by glacier as well as facies. Different years are indicated by colour while different glaciers are indicated by symbol. (a) and (b) show albedo for entire glaciers, (c) and (d) show albedo for accumulation facies, and (e) and (f) show albedo for ablation facies. (a), (c), and (e) present averages of the data in (b), (d), and (f). As discussed above, from earlier studies, all measurements should be considered to have error bars on the scale of ~0.1. Whole Glacier: Accumulation Facies: Ablation Facies:            "    % !#!#             "    % !#!#            "    % !#!#  89 7. Glacier Albedo on the Brøgger Peninsula, Svalbard was a very shallow snowfall, so the underlying darker layers would not have been completely masked. Nevertheless, both the individual glacier measurements and the whole-image averages do document seasonal darkening of snow. The albedo of ablation facies, as expected, does not demonstrate a significant seasonal signal. Variability will be determined by differences in debris cover and the amount of liquid water present on the glacier surface. While some brightening would be expected in the autumn, this time period is not sufficiently resolved to see a signal. As with TAAR observations, for all of Figure 7.3, individual glaciers behave similar to each other, but there is no consistent ranking of the three glaciers’ albedos. 7.3 Interannual Observations The seasonally-displayed data (see Figure 7.3) are a starting point for investigation potential interannual variations in albedo which could contribute to a melt-albedo feedback. A consistent change in the timing or magnitude of seasonal effects would be indicative of a significant role played by the surface energy balance in the glaciers’ changing behaviour. However, no trend in recent darkening is observable in these data. Albedo measurements are plotted versus year (see Figure 7.4), rather than day of year to visualise potential long-term albedo. This shows the full length of the 35-year record, but also shows how sparse it is for much of the time. From the melt-albedo hypothesis, there is no expected long-term change in albedo for ablation facies (see Figure 7.4c), and none is observed. The hypothesis would be confirmed by a meaningful trend in albedo for accumulation facies (see Figure 7.4b). A trend in darkening in the last decade could potentially be interpreted from this time-series (in Figure 7.4a and b), but this is easily explained as the result of a high albedo observation on one day in 1999 and is not robust. Intra-annual changes in albedo overwhelm interannual variability. Therefore, an attempt was made to remove the seasonal effects using the 3rd order polynomial fit to the data in Figure 7.3a (see Section 7.2). This also removes the confounding and potentially spurious influence created by the spring measurement at the very beginning of the time series and the late-summer measurement at the very end. Only the whole glacier measurements were corrected because meaningful seasonal parameterisations could not be fit to the albedo values for accumulation and ablation facies. Nevertheless, considering Figure 7.5, there is no observable long-term trend in albedo. Thus, there is no evidence in the Landsat record which supports a melt-albedo feedback in high elevations of Svalbard glaciers. The whole-glacier albedo is related to mass balance, both summer balance and annual balance (Greuell et al. 2007). Brøgger peninsula glaciers have been shown to have had a mostly negative, but largely consistent, mass balance since 1968, although there is a slight trend to more negative mass balances since 2000 (Barrand et al. 2010). This trend could manifest itself in Figure 7.4a. However, as has been explained, any trend which may be interpreted would depend on just one observation from 1999 and therefore is not robust. In Figure 7.5, the seasonal effects are largely removed, so anomalously high and low values are easy to identify. However, for both Figures, the albedo measurements are not exclusively from the end of the melt season and so cannot be compared with annual mass balance measurements. 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 90              !#$ " !#$ " !&                                 Fi gu re 7 .4 : L on g- te rm al be do es tim at es fr om L an ds at im ag er y f or th re e s tu dy g la ci er s i n Sv al ba rd ; d iff er en t g la ci er s a re in di ca te d by th e s ym bo ls in th e l eg en d in (a ). M ea su re m en ts ar e av er ag ed o ve r t he w ho le g la ci er (a ), ac cu m ul at io n fa ci es (b ), or a bl at io n fa ci es (c ). A s d isc us se d ab ov e, fro m e ar lie r s tu di es , a ll m ea su re m en ts sh ou ld b e co ns id er ed to h av e er ro r b ar s on th e s ca le o f ~ 0. 1. a) b) c) 91 7. Glacier Albedo on the Brøgger Peninsula, Svalbard              !  !  #   Figure 7.5: Long-term albedo estim ates from Landsat im agery for three study glaciers in Svalbard. Th e data from Figure 7.4a have been corrected for seasonal effects using the 3rd order polynom ial fit in Figure 7.3a. A ll m easurem ents should be considered to have error bars of at least 0.1. 7. Glacier Albedo on the Brøgger Peninsula, Svalbard 92 7.4 Discussion While other satellite-derived albedo studies greatly benefitted from in situ observations to aid in interpretation (e.g. Knap et al. 1999a; Hendriks & Pellikka 2004; Stroeve et al. 2006; Box et al. 2012), that was unfeasible for this study. Attempts were made to obtain any such data from colleagues at the Norwegian Polar Institute, but without success. Even without in situ data, the albedo estimates in this study are consistently darker than those reported in earlier works. This can be explained as a combination of many factors, the largest of which is likely to be related to the reflective behaviour of snow and ice surfaces. As discussed above, no BRDF correction has been applied, therefore nadir darkening is expected, as documented by Hendriks and Pellikka (2004). Other factors which could skew albedo estimates include assumptions made in atmospheric correction (i.e. negligible adjacency effects, extrapolation of atmospheric variables), slope effects, and inaccuracies in the narrow-to-broadband calculation in Equation 2.3. These are considered in turn. Atmospheric correction has such a small impact (~2%, as described above) that although it improves the results of the study, it will nevertheless be overwhelmed by the magnitude of other corrections (~10%). Slope effects are also likely to contribute to the darkening of the albedo estimates calculated here. Again, however, because slope will be largely consistent, it will not impact the evaluation of seasonal or interannual trends. Finally, there is the narrow-to-broadband conversion. Equation 2.3 has been shown to work better for low albedos than high albedos, and the data here do not at any point stretch the upper bounds of the parameterisation. Although there was the possibility that MSS estimates might slightly overestimate albedo, this effect does not perturb interpretation of the time series. Nevertheless, there is no evidence in the time series which supports a melt-albedo feedback in high elevations of Svalbard glaciers. While it may be the case that this hypothesis is not valid, this is by no means proof against the melt-albedo hypothesis. The James et al. (2012) study reaches back to 1961; it is therefore possible that this Landsat record is simply not long enough to observe a change, or that the length of the record makes the observable darkening would be smaller and therefore harder to detect. Still, there would be an accelerating signal which is not observed here. Alternatively, it is possible that a signal is hidden within the error margins of the estimated albedo. No large implications are therefore able to be discussed for Svalbard because no hypothesis was definitely proven or disproven. It does open the possibility that dynamic changes, rather than albedo feedback, are responsible for thinning of glaciers in Svalbard. Thus, changes in accumulation and dynamic response must be better modelled to understand future changes to glaciers in Svalbard and the wider high Arctic. However, this is still largely speculative and opens many future research directions (see Section 8.5). 7.5 Section Summary Using the Landsat data archive, this chapter investigated the seasonal and interannual albedo variability of the accumulation and ablation areas for three glaciers on Svalbard’s Brøgger Peninsula from 1976 to 2011. Landsat data were converted to at-sensor reflectance using information contained in image metadata, to at-surface reflectance using an implementation of the 6S atmospheric radiative transfer model, and to glacier albedo using an empirically-determined parameterisation. At-surface reflectances 93 7. Glacier Albedo on the Brøgger Peninsula, Svalbard were also used to classify glacier accumulation and ablation areas. Classification from time-lapse cameras (for three days in 2010 and 2011) showed that classification was largely successful, and near-coincident image acquisition on two days showed the consistency of Landsat albedo estimates, despite omission of BRDF and slope-angle corrections. Albedo was considered for each facies and well as for the entire glacier, and for each glacier individually as well as for all three together. From the data: • TAAR consistently decreases through the spring and the transition to autumn is observed for one year. • A seasonal darkening signal is seen strongly for the full-glacier observations. • A seasonal parameterisation is used to remove this effect to aide in interannual, long-term trend analysis. • Seasonal darkening can also be seen for accumulation facies. • The albedo of ablation facies, as expected, does not demonstrate a seasonal signal. • All three individual glaciers behave similarly, but they demonstrate no consistent relative behaviour. • Although no in situ albedo measurements are available, nadir-darkening is attributed as the most likely cause that albedo measurements are lower than observed in earlier studies. • No long-term change in albedo for the whole glacier, ablation facies, or accumulation facies is observed. Intra-annual changes in albedo overwhelm interannual variability. Therefore, although the classification and albedo estimates are meaningful and robust, there is no evidence in the Landsat record which supports a melt-albedo feedback in high elevations of Svalbard glaciers. This could indicate that thinning of higher elevations on Svalbard glaciers is the result of dynamic thinning, thus drawing attention to interannual variability of accumulation rather than reflective properties. However, it is too soon to draw any meaningful conclusions, and further research is needed. 94 95 8. Discussion and Conclusions 8. Discussion and Conclusions The aims of this thesis were to conduct a focused analysis of glacier surface classification using in situ spectral reflectance data, to understand the implications for a range of multispectral imagers, and subsequently to apply these results to a case study. This research was motivated by the wide utility and importance of multispectral remote sensing of glacier surfaces. In addition, because multispectral imagers were not designed with glaciological parameters in mind, and because current techniques for glacier surface classification were not based on field data, questions were raised about the efficacy of current practices. This chapter summarizes and synthesises the discussions and conclusions of earlier chapters, thereby contextualising the investigations of this thesis with each other in addition to the wider literature. The results of previous chapters are first revisited (Section 8.1), and the original research objectives are then reconsidered (Section 8.2). This chapter also discusses cross-cutting considerations of this work such as transferability (Section 8.3) and downstream scientific impacts (Section 8.4) before moving on to future research directions (Section 8.5) and a final summary (Section 8.6). 8.1 Synthesis of Results 8.1.1 Using In Situ Spectra to Explore Landsat Classification of Glacier Facies Full-spectrum in situ surface reflectance data from a range of end-of-melt-season glacier surface zones on Midtre Lovénbreen (Svalbard) and Langjökull (Iceland) were combined with Landsat ETM+ imagery to inform, improve, and better understand satellite classifications of glacier surface facies. While many previous studies have attacked the problem of characterising the spectral properties of snow and ice, this is one of a very few to simultaneously view the entire range of classes present on a glacier surface. The unique perspective taken by Chapter 5 was an inverse approach to glacier surface classification, using in situ data to drive interpretation and analysis of satellite imagery. Highlights are bulleted below. • Comparison of in situ data and unsupervised classifications of Landsat imagery led to the conclusion that while VNIR (~350-1400 nm) compared to SWIR (~1400-2500 nm) is important for delineating glacier extent, variation within the VNIR is important for distinguishing these glacier surface types. • Based upon principal component analysis of the in situ spectra, a simple technique (using only two inputs for clean glaciers and three for ash-covered glaciers) for unsupervised classification of glacier surfaces using Landsat imagery was derived. • For clean glaciers, full-spectrum clustering was only slightly more precise than clustering limited to wavelengths which were measured by Landsat. For ash-covered glaciers, the hyperspectral data were important. 8. Discussion and Conclusions 96 • The linear combinations presented here for Landsat TM / ETM+ are transferrable across glaciers and time periods. Output classes are suited for a range of glaciological applications. These conclusions incrementally simplified and refined previous work, making meaningful glacier surface classification one step more accessible to glaciologists. In contrast to prior work, this section focused spectral investigations almost solely on the VNIR wavelengths to the exclusion of the SWIR. These conclusions opened questions about wider transferability across sensors and the capacity that each individual imager does or does not have. These questions were addressed in the subsequent chapter. 8.1.2 Impact of Spatial and Spectral Resolution on Glacier Surface Classification Increasing availability of multispectral data requires that researchers know what data types are best suited for their own research questions. For glacier surface classification, the radiometric, spectral, and spatial properties of a suite of popular sensors (ATM, ASTER, MSS, TM, ETM+, OLI, Sentinel-2, and MODIS) are investigated using data sets in common to perform controlled analyses. Linear combinations for all sensors were created based on principal component analysis of in situ spectra. Among these sensors, spectral resolution and range or radiometric resolution were not important on their own. The most important property which determined the suitability of a multispectral imager for glacier surface classification was its radiometric range. In particular, it was found to be beneficial to have a low gain in the visible and a higher gain in the NIR. Spatial resolution can, seemingly paradoxically, be either beneficial or detrimental to classification, depending on the fractal scale of the surface being classified. For glaciers, it was found that low resolution imagery (250 m pixels) is too coarse to represent the true complexity present on a glacier. However, medium resolution imagery (60 m, 30 m, or 20 m) did accurately represent the results derived from high resolution airborne imagery. Nevertheless, 30 m imagery was preferable to 60 m imagery. From this, it was inferred that inhomogeneity on glaciers is significant on a scale between ~60 and 250 m. Based upon these radiometric, spatial, and spectral requirements, the satellite imagers currently in use that are most suitable for glacier surface classification are Landsat TM / TM+ (with gain settings of LLLH) and ASTER (low 1 or normal gain). Both Sentinel-2 and the OLI on Landsat 8 are also expected to be similarly qualified. Nevertheless, because of the encouraging suitability of a range of sensors to glaciology, data availability will likely remain the most important consideration. 8.1.3 Glacier Albedo on the Brøgger Peninsula, Svalbard Using the Landsat data archive, the seasonal and interannual albedo variability of the accumulation and ablation areas for three glaciers on Svalbard’s Brøgger Peninsula from 1976 to 2011 were investigated. Landsat data were used to both estimate glacier albedo using an empirically-determined parameterisation as well as to classify glacier accumulation and ablation areas. Comparison with imagery from time-lapse cameras showed that Landsat classification was largely successful, and near-coincident image acquisition on two days showed the consistency of Landsat albedo estimates, despite omission of BRDF and slope- angle corrections. 97 8. Discussion and Conclusions Observations based on this time series included the following: • TAAR consistently decreased during the spring, and the transition to autumn was observed for one year. • A seasonal darkening signal was seen for both full-glacier observations and accumulation facies only. As expected, this was not visible in ablation facies. • All three individual glaciers behave similarly, but they demonstrate no consistent relative behaviour. • Although no in situ albedo measurements are available, nadir-darkening was attributed as the most likely cause for albedo measurements to be lower than observed in earlier studies. • No long-term change in albedo for the whole glacier (raw and corrected for seasonal effects), ablation facies, or accumulation facies was observed. Intra-annual changes in albedo overwhelmed interannual variability. Ultimately, although the data showed interesting information about the seasonal behaviour of these glaciers, there was no evidence in the Landsat record which supported a melt-albedo feedback in high elevations of Svalbard glaciers. 8.2 Achievement of Research Objectives This thesis provided motivation in Section 1.1 and enumerated several questions in Section 1.2, and it is against those metrics that the success of this investigation should be measured. In direct response to those questions, the achievements of the research thus far are stated below. Landsat classifications were qualitatively compared with field data to yield preliminary insights, the optimal wavelengths for investigations of glacier surface classification were successfully identified using principal component analysis (i.e. primarily VNIR broadband reflectance and secondarily the difference between blue/green and red/NIR reflectance), and this information was used to create simple and transferrable linear combinations of satellite bands which highlight these features. Limitations were identified for this framework, namely that it is suitable for clean glaciers rather than ashy glaciers. As discussed above, the spectral and radiometric properties of a range of sensors were successfully studied and their impact upon clustering was successfully quantified. Similarly, so was the impact of variable spatial resolution. The limitations thereof were discussed. Particular recommendations were also made for the most appropriate sensors for multispectral classification of glacier surfaces. However, the question of implications for glaciological applications were not considered in Chapter 6; they are addressed below (see Section 8.4). Next, the classification scheme was used for application to three glaciers in Svalbard to understand their albedo over the Landsat record, in particular to investigate the potential role of melt-albedo feedback in high-elevation thinning observed on Svalbard glaciers. While seasonal trends in glacier surface cover and albedo were successfully observed, no long-term trend in albedo was found. Thus, no evidence was found in the Landsat record to support melt-albedo feedback on Svalbard glaciers. This null result still successfully achieves the research aims, opening many future research questions (see Section 8.5). 8. Discussion and Conclusions 98 It is important to consider whether the successful answering of these research questions moves forward the wider research field which provided motivation for the study. Put another way, as a result of this thesis, do we know more about multispectral glacier surface classification? And do the results presented here progress the field? The answer to both questions is ‘yes.’ Although multispectral imagers were designed around physical limitations (i.e. atmospheric transparency) and particular tasks (e.g. agricultural monitoring), we now have evidence that these restrictions have not significantly impacted glacier surface classification. Indeed, this thesis does not need to provide suggestions for alternative bands because the current sensors perform satisfactorily. Thus, it is confirmed that we can have confidence in previous results and continue with similar techniques as well as those presented here. Transferability of results motivated some of this work and is addressed explicitly in Section 8.3. Indeed, methods developed in this thesis were successfully used to investigate processes on Svalbard glaciers. Although widespread conclusions eluded this case study, progress has been made and future research is already being planned. However, in assessing success of transferability, it is necessary to question an early assumption, namely that the simplest achievable classification scheme is the best approach suited to these tasks. The work presented in this thesis has demonstrated some success with unsupervised clustering (ISODATA and k-means), but because improvements in sensor technology are neither expected nor demonstrated to be necessary, future focus should move back to classification technique. Perhaps a rethinking, and study, of multi-step, hierarchical classification may be necessary to continue to improve glacier surface classification. In this way, improvement can be made in applications such as mass balance measurements and hydrological modelling, for example. Future research questions are addressed in Section 8.5. Ultimately, the investigations brought together and linked through this thesis are a wide consideration of multispectral imaging of glacier surfaces. From field data and theory through to the Landsat record and application, this dissertation contributes to the understanding of glacier facies, multispectral classification, and Arctic glacier behaviour. 8.3 Transferability of Results With any experiment or observation, the immediate questions are what generalisations can be made and to what extent conclusions can be extrapolated to a larger area. This is a topic relevant to all three results chapters presented earlier. Specifically, how representative are the field spectra of the glaciers they were measured on? How representative are the surfaces of Langjökull (spectrally) and Midtre Lovénbreen (both spectrally and spatially) of glaciers in general? Are the techniques presented here realistically able to be used in a wider context? And how representative are the three glaciers of the Brøgger Peninsula of the wider population of Svalbard glaciers? While it is beyond the scope of this study to definitively answer these questions, they must at least be considered. The selection of field spectra sampling locations was based upon the exploration of the field party. For Midtre Lovénbreen, the glacier is small and therefore it is highly unlikely that any major classes were omitted. Langjökull is much larger, and measurements were limited to a single outlet. Nevertheless, Landsat classification of this outlet indicates the presence of the full range of facies along the transect which was used, and therefore it is unlikely that any major classes were omitted there, either. 99 8. Discussion and Conclusions The question then turns to the relative proportion of each class as measured; to an extent, it is important to recognize that these relative proportions will have some impact on the statistics, in particular the proportions in any contingency table. In the case of Table 5.1 (considering full-spectrum vs. Landsat-simulated classification of clear glacier spectra), the vast majority of the data fall along the primary diagonal, which would not change. The off-diagonal components would change if the proportions of spectra were to change, but likely not significantly enough to impact any conclusions. So, classification success may be somewhat artificially high or low, depending on the real proportion of spectrally confusing classes. Nevertheless, there is no reason to expect bias from the samples that were collected, and the very high classification success minimises the effect any potential impact of sampling bias on the conclusions which were drawn. The ranking of sensors according to K values, however, could conceivably have been more impacted by different proportions of facies. For example, for Midtre Lovénbreen, inclusion of more ‘coarse snow’ and ‘dry ice’ spectra would likely have depressed all K values. Similarly, for Langjökull, including more spectra from the classes near the ‘white ice’ spectra could have had a similar impact. However, although the magnitude of the statistics would have changed, this would have impacted (beneficially or detrimentally) all simulated clustering analyses, and it is therefore unlikely that the conclusions thereof (based on relative rank) would change. This study based conclusions on all available data, choosing not to filter out spectra. Nevertheless, potential impacts of relative proportions of classes would be testable with further sampling campaigns. This study began with the hope of improving both the quality of multispectral classification of glacier surfaces as well as the efficiency of the techniques to do so. The outputs of the resulting classification scheme corroborated rather than replaced the results of earlier work. In other words, recommendations were not made which impacted the results of the classifications, simply the way in which they were achieved. So, direct comparison of classification results was not undertaken; indeed, there are no ground measurements with which to assess such a comparison. Instead, incremental improvement in efficiency is the contribution of this work. Future work will have to consider whether such an improvement is helpful for individual applications. It is crucial for this study that the spectra measured on Midtre Lovénbreen and Langjökull (and principal components thereof) are considered to be representative of other glaciers. Because spectra are a result of a combination of physical processes (snow formation, accumulation, compaction, melt, etc.), it is reasonable to expect that this will be the case. Indeed, the similarity of the first and second principal components of spectra between the two glaciers supports this interpretation. In addition, preliminary principal component analysis of satellite imagery from other glaciers at other times yield nearly identical principal component band combinations. Indeed, for a particular glacier, if the user wanted a highly customised band combination, it would be reasonable to use PC1 and PC2 of a site-specific PCA. Thus, because spectra from the two chosen field sites and satellite images from others agree upon the principal components, it is reasonable to assume that Midtre Lovénbreen and Langjökull are spectrally appropriate representative glaciers for study. Nevertheless, although the major classes of Midtre Lovénbreen and Langjökull will be spectrally representative of many glaciers, there are still surface features on other subsets of glaciers that are not considered but could and would impact transferability of results. The largest subset of glaciers will 8. Discussion and Conclusions 100 be those in colder climates, in particular those with dry snow facies, for example in Greenland and Antarctica. The metamorphism of the snow moves at a much slower pace, leading to brighter surfaces. There will also be areas of wind-scour, wind glaze, sastrugi, and blue ice zones (e.g. Kuchiki et al. 2011; Scambos et al. 2012; Dixon et al. 2013). While melt facies give information about mountain glacier mass balance, these facies indicate that other processes are at work (e.g. persistent wind patterns). It is entirely possible that the LCs presented here will be appropriate for classification of these facies; earlier work (Boresjö Bronge & Bronge 1999) classified snow and ice zonation in Antarctica (even including sea ice) using PCA as a guide. More work will be needed to confirm or deny this hypothesis, but that is beyond the scope of this study. Another subset of glaciers that should be considered is those which are covered by ash and debris. Ash cover and some debris cover was considered for Langjökull, and for a unique group of glaciers, volcanically-influenced facies will confound traditional facies classification while presenting other interesting information related to radiative forcing and implications for sudden mass balance change. The reflective properties of ashy glaciers will depend not only on concentration in the snowpack, but also the wetness, the ash lithology, the ash grain size and shape, and other site- and even eruption-specific factors. Ultimately, the classifications described here will not be applicable to all ashy glaciers; further field and modelling studies are necessary for a fuller understanding. Similarly, a range of site-specific factors will influence the understanding of glacier debris cover. However, debris cover-related facies are largely found in the ablation area of glaciers. Therefore, although they do not fit the classification LCs presented here (e.g. they will benefit from contribution of information in the SWIR), if a glacier outline is already available, their distribution will not severely influence the glaciological applications discussed here. Again, research is ongoing in this regard, especially with application to Himalayan glaciers (e.g. Casey et al. 2012; Racoviteanu & Williams 2012). Additionally, there are surface classes that don’t completely cover the snow or ice surface or totally mask it, but instead combine with it. Namely, dusty snow, black carbon, and snow algae. As discussed in previous work (e.g. Takeuchi 2009; Painter 2011; Painter et al. 2012; Kokhanovsky 2013), these surface types will have a seasonal impact on surface variability and energy balance. These surfaces will have a somewhat different spectral response from pure snow facies, but the severity of the change will be dependent on the concentration of the contaminants. Likewise, whether the LCs developed here will continue to be effective will depend on contaminant concentration. Because the soot/dust/algae will each be less variable than ash cover, it is much easier to use spectral information to gain quantitative information (e.g. radiative forcing, Painter et al. 2012) about the contamination. These contaminants are much more widespread than some of the other facies discussed above, and they will have influenced the spectral sampling considered in this study to a small extent; this effect is not expected to be significant because contaminants were not observed by eye. Further investigation and progress in quantifying their impact requires particular physical and chemical sampling strategies and remote sensing applications will likely require hyperspectral data rather than the multispectral data, albeit not exclusively. The analysis of the impact of spatial resolution on classification success (see Section 6.2) raises the question of whether Midtre Lovénbreen is also spatially representative of other glaciers. Processes which impact the spatial distribution and scale of facies include the accumulation distribution, wind and avalanche redistribution, melt patterns, and local slopes, to name a few. In these regards, there is nothing 101 8. Discussion and Conclusions that sets Midtre Lovénbreen apart as a special glacier. By contrast, factors such as crevasses or incised supraglacial streams are glacier-specific and may have some impact on the impact of resolution on surface classification. However, such features are small on Midtre Lovénbreen and would likely manifest themselves as an important difference between high- and mid-resolution images, rather than influencing the conclusion that low-resolution images are inappropriate for facies classification. Thus, while there is no reason to suspect that Midtre Lovénbreen would have a different spatial character than other glaciers, it is recognized as a limitation of this study. Investigation of a wider study area would confirm or confine this extrapolation, so it is suggested as a potential future research direction. Ease of implementation and transferability were invoked in the motivation of this research, so it is only fitting that they are assessed as it comes to a close. While linear combinations are provided for a range of multispectral sensors, are these indeed enough to ensure that the results are transferrable? They are a first step to ensure that the study is easy and fast to implement (no intermediate principal component analysis needed), and fewer components being input into a classification is always easier. However, although the basis of the classification is unsupervised, the classification did require manual glacier-tracing in this case. Wider inventories from, for example, the Randolph Glacier Inventory (Arendt et al. 2012) could be used to broaden the scope of the study area. As with previous studies (De Woul et al. 2006; Braun et al. 2007), there is necessary user-involvement in merging the classes, which is an area which requires further development. Considerations for an objective measure of merging could be considered here (e.g. relative Euclidean distance, F-test, etc.), but at some level this is merely providing more information for a subjective, expert-made decision. As in many remote sensing applications, such as glacier delineation (Paul et al. 2013), there is recognition that the best results will often require human intervention at some stage. While spectrally and spatially Midtre Lovénbreen can be considered to be a representative glacier, are Midtre Lovénbreen and its two neighbours actually representative of the behaviour of Svalbard glaciers, in particular with respect to climate, AAR, and albedo? Because of their proximity to the Ny- Ålesund research station, many, many glaciology studies have assumed that the behaviour of glaciers on the Brøgger Peninsula is applicable to much of the rest of Svalbard and the world. However, it is important to recognize that they are a biased sample in many regards, especially in terms of geography (western side of the archipelago), in aspect (all North-flowing), and in size (all fairly small). While their small size is the reason they are of interest to a Landsat study, their geography definitely has implications for their behaviour, resulting in unique insolation and precipitation. That acknowledged, the outcomes of this study were to observe seasonal behaviour successfully and not prove a (melt-albedo feedback) hypothesis. The main result of suggested research directions is as valid for the Brøgger Peninsula glaciers as elsewhere, although they will remain a very good subset of glaciers for case studies. 8.4 Potential for Glaciology Efficient and effective glacier surface classification has many implications for glaciological research. The most immediate use is application to widespread measurements of AAR (and therefore ELA) as mass balance proxies; this is given explicit treatment in Section 8.4.1. Some of the studies upon which the classification method developed in this thesis were used for 8. Discussion and Conclusions 102 validation of glacier melt and hydrology models (De Woul et al. 2006; Braun et al. 2007). The effective identification of wet facies on clean glaciers (see Chapter 5) predisposes this classification scheme to effective application to hydrological applications. Increased application to validate models studying glacier surface hydrology would therefore be appropriate. In addition to small mountain glaciers (e.g. Dahlke et al. 2012), there is increasing interest in water-saturated areas in Greenland, so this may be a promising research direction. Energy balance models may find some overlap with hydrological modelling, driving the liquid contributions to the glacier system. Van Angelen et al. (2012) demonstrated the importance of, for example, albedo parameterisations in regional mass balance models. Raised as a possibility by other studies (Machguth et al. 2009; Dumont et al. 2012), using classes derived from multispectral imagery could be a simplified method to drive particular parameterisations or provide validation for the models. This study provides another step in the direction of successfully applying such an assimilation or validation mechanism. In this consideration of downstream impacts, it is notable that there is no sudden paradigm shift. This is not unexpected. The full spectrum analysis conducted here demonstrated that the spectral information measured by many satellite sensors does not contribute meaningfully more than in situ spectra. As such, the methods presented will enable classification to be more incrementally more efficient and repeatable, but the results will not largely differ from those attainable previously. Other areas of research which will be impacted eventually include, as mentioned earlier, studies of climate variability, understanding of water resources, study of geomorphologic hazards, and investigation of high altitude and latitude biodiversity. These all tie back to better process-based understanding of glacier surface processes enabled by application of multispectral remote sensing. 8.4.1 Mass Balance Proxies A more quantitative treatment can be given to the impact this work will have on mass balance proxies. As described in Section 2.2, glacier surface classification can be used to calculate the AAR and combined with a DEM to determine ELA, both of which are proxies for glacier mass balance. Done either on an annual timescale at the end of the melt season (e.g. Rabatel et al. 2005; Rabatel et al. 2008; Shea et al. 2012), or theoretically ongoing throughout a melt season (Dyurgerov 1996), these are very helpful glaciological parameters obtainable from remote sensing. However, what uncertainty levels are there on these parameters? The linear regression between mass balance and AAR has its own uncertainties, but we can use the clustering done in this study to look at just the uncertainty of just the classification. AAR itself is the ratio of accumulation area (c) to the entire glacier area (AG): AAR = c / AG (8.1) Considering the effect of potential errors, this means that the lower bound on AAR (AARL) is determined by: AARL = (1 - EO)c / AG (8.2) 103 8. Discussion and Conclusions where EO is the error of omission. Similarly, the upper bound on AAR (AARU) is determined by: AARU = (1 + EC)c / AG (8.3) where EC is the error of commission. Therefore, the uncertainty in AAR measurements is a range based on the magnitude of omission and commission error percentages. These errors are considered using A (rather than K or A*) because the subject of interest is absolute success, not success relative to chance. As mentioned above, the number of spectra which were measured is not representative of the relative area of classes. Nevertheless, considering the percentage of omission and commission results for the range of band-simulated spectra (see Chapter 6.1) yielded TAAR “measurements” of ±5%. Different spatial resolutions (see Chapter 6.2) also contribute to uncertainty in AAR measurements. Taking the highest resolution of ATM image (2 m pixels) classification as the basis by which to assess lower resolution estimates of AAR, a larger proportion of commission than omission errors was observed. It is entirely possible that this is the case only for Midtre Lovénbreen or indeed just for the particular ATM images used in this study, so no wider conclusions are drawn from this point. Again however, AARs were found to differ in most cases by less than ±5% and at the extreme ±10%. However, these are just simulated comparisons without ground-truth. Figure 7.1 gives the chance, however, to compare TAAR from Landsat and in situ photography for Austre Lovénbreen. For 20 August, 2010, TAARIn Situ=0.35 and TAARLandsat=0.34; for 23 June, 2011, TAARIn Situ=0.99 and TAARLandsat=1.0; and for 17 August, 2011, TAARIn Situ=0.16 and TAARLandsat=0.08. From these comparisons, it appears that AAR is very much within the 10% error bound, and largely within the 5% envelope. The almost exact agreement of the 20 August and 23 June classifications may be coincidental, but it is nonetheless encouraging for the interpretation of Figure 7.2. There may be some magnitude-dependence to AAR error, as glacier surfaces will become more complex as AAR decreases. Nevertheless, this is just a snapshot of three dates, and there are likely many other cameras on other glaciers which could be used to analyse these numbers more in depth. 8.5 Future Research Excitement in research doesn’t stop with results. It is just the opposite. It propagates with wider implications and the many branches of future research directions. This section aims to consider some of the avenues opened by the research presented in this thesis and address where, given the time and support, the next steps would lead. Chapter 5 is one of the first studies to provide a full data set of spectra spanning the entire Landsat VIR spectrum, at high spectral resolution, with a complete range of surface types present on an Arctic/temperature glacier at the end of the melt season. Following on from this work, the author is already is discussions with colleagues at the University of Reading to archive the spectra within a larger spectral library, as well as with a colleague at the University of Southampton to use the spectra to help in a new project which aims to invert atmospheric aerosol properties from Landsat data. Indeed, because there is only one other study so similar to this that the author is aware of (Hendriks & Pellikka 2004), the logical next step is to continue to widen the available spectral library which describes the full range 8. Discussion and Conclusions 104 of facies present at the end of a melt season on a glacier’s surface. Such measurements would be able to further investigate and quantify the role of ash on glacier reflectance and classification (e.g. to refine the continuum observed on Langjökull), further test the transferability of the methods presented here, address the statistical considerations proposed in Section 8.4, and obtain spectra for classes which were not able to be measured here. In particular, the classifications considered here did not measure firn or superimposed ice, both of which would be very desirable for glaciological studies to be able to more reliably classify. With an interest in firn and superimposed ice, it may also be beneficial to combine multispectral classification with glacier facies classification based on synthetic aperture radar (SAR). SAR classification zones have been shown to be related to surface layer water content, snow depth, melt history, and facies extent from prior years (see Section 2.6.4). Some overlap will be had between the two classification techniques, as the multispectral groupings of clean glaciers also highlighted melt and water content, to an extent. Firn and superimposed ice only manifest at certain timeframes in the melt season and are difficult to visualize, therefore having a remote sensing technique which sees through clouds will result in a more reliable method. As with the framing of this multispectral study, many approaches have been taken to investigate SAR classification of glacier surfaces, but it still remains an active research area. Therefore, further study in this direction would have to be very carefully framed and would benefit from extensive field data to minimise speculative interpretations. In addition, further coincident measurements and ground-truth would facilitate the application of more complex classification schemes or potentially the application of object-based image analysis, either to multispectral or SAR data. Nevertheless, multispectral data will continue to have a very strong role in glacier surface classification. Use to delineate glacier extent, investigate the extent of glacier facies, measure mass balance proxies, and more will only continue to grow as more and better multispectral sensors’ data become available in the coming years. Teaming up with glacier modellers to integrate the facies and albedo measurements used here to improve local and regional mass balance or energy balance models would be a fruitful and exciting research direction, as well, thereby diversifying the use and impact of multispectral glacier surface classification. Chapter 6 used two data sets to investigate the wider application of glacier surface classification across many platforms. For the new sensors whose properties are not fully quantified (OLI on Landsat 8 and Sentinel-2), it will be very exciting to confirm the conclusions drawn here from the assumptions made here. Further principal component analysis of multispectral remote sensing imagery of glaciers would confirm the linear combinations of bands presented here for surface classification. In addition, the transferability of the conclusions about the necessity of medium-to-high spatial resolution has wide implications. It is very important from both pragmatic and efficient points of view to identify the lowest reasonable resolution of data to use for studying glacier surface properties. This is especially true if the techniques are going to be used across wide areas of long time series. Although the discussion above explained why it is not expected that the conclusions will change significantly, it would still be beneficial to confirm with ATM images from other glaciers, in other areas, and across larger areas as well. Beyond application of ATM images, coordinated campaigns allowing for intercomparison of coincident, multi- resolution data spanning UAVs, airborne platforms (e.g. with IceBridge or ARSF), and satellite images would lead to even more robust conclusions in the future. 105 8. Discussion and Conclusions Moving to the final section of the thesis, despite or perhaps because of the null result obtained investigating melt-albedo feedback on the Brøgger Peninsula glaciers, Chapter 7 opens potentially the most exciting future research questions. Although the Landsat record did not provide evidence for melt- albedo feedback on Svalbard glaciers, it does not necessarily mean that such a signal does not exist. Perhaps a wider study area is needed, although this is unlikely, because James et al. (2012) included Austre Brøggerbreen and Midtre Lovénbreen in their study. It was noted than interpretation was hampered by the sparse nature of the data. Therefore, although a shorter temporal record and at lower spatial resolution, investigation of albedo with sensors having higher temporal resolution (e.g. AVHRR, MODIS, and VIIRS) could prove to be beneficial. In addition to searching out a potential albedo signal which may not exist, it would be interesting to test the second hypothesis: dynamic thinning of Svalbard glaciers. In short, an increased load of accumulated of snow can cause an acceleration of glacier flow, thereby thinning the glacier. In other words, the surface elevation is a balance between melt/accumulation and redistribution of ice through glacier flow. Mass-budget approaches have been used in a variety of applications (e.g. Nuth et al. 2012; Fischer 2011). Elevation changes have already been calculated by James et al. (2012), and downscaled glacier melt models have been optimized to investigate the climatic (surface) mass balance of Svalbard glaciers (Rye et al. 2010). Therefore, the dynamic hypothesis for thinning in higher elevations of Svalbard glaciers is eminently testable. By definition, a thesis is meant to push forwards the boundary of knowledge in a particular field. This dissertation does so for multispectral classification of glacier surfaces, presenting new data, new methods, and a new application of Landsat data. These investigations yield many promising future research directions, including expanding glacier spectral libraries, further quantifying glacier surface classes, combining a range of remote sensing data sources to improve classification schemes and further test their transferability, and looking more in depth into the processes leading to enhanced thinning at high elevations of Svalbard glaciers. 8.6 Summary In final summary, this thesis aimed to conduct a focused analysis of glacier surface classification using in situ spectral reflectance data, to understand the implications for a range of multispectral imagers, and subsequently to apply these results to a case study. It achieved these research objectives directly by: • Presenting a new set of field spectra from a wide range of glacier facies, • Employing principal component analysis to develop linear combinations highlighting the wavelengths important for distinguishing glacier surface classes for Landsat and many other multispectral sensors, • Identifying the central role that radiometric range plays in wider application of multispectral classification of glacier facies, • Quantifying the impact of spatial resolution on glacier surface classification and thus recommending the use of medium and high resolution imagery for this purpose, 8. Discussion and Conclusions 106 • Describing, as a case study, seasonal cycles of facies transitions and albedo characteristics of three glaciers in Svalbard through the Landsat record (1976-2011), and • Failing to observe any long-term trend in decreasing albedo of accumulation facies of glaciers of the Brøgger Peninsula, Svalbard. In addition, themes which transcend the entire thesis were considered, including wider transferability and further applications within glaciology. The conclusions (both spectral and spatial) were widely determined to be transferrable, although further study would of course be necessary to confirm or deny this. The errors for application to estimating accumulation area ratio using the methods described here were quantified in particular. Additional future research directions were identified, including proposals which called upon integrating multiple data sets to refine surface classification schemes, studying their wider transferability, and test further hypotheses for the case study of accelerated thinning of high elevations of Svalbard glaciers. In sum, the work presented in this dissertation has contributed to the understanding of glaciological applications of multispectral remote sensing imagery, a field which will be sure to remain innovative and vital to glaciology for many years to come. 107 8. Discussion and Conclusions 108 109 9. References 9. References Albert, T.H., 2002. Evaluation of Remote Sensing Techniques for Ice-Area Classification Applied to the Tropical Quelccaya Ice Cap, Peru. Polar Geography, 26, pp.210–226. Anderson, J.B., 1999. Antarctic Marine Geology, New York: Cambridge University Press. Van Angelen, J.H. et al., 2012. Sensitivity of Greenland Ice Sheet surface mass balance to surface albedo parameterization: a study with a regional climate model. The Cryosphere, 6(5), pp.1175– 1186. De Angelis, H., Rau, F. & Skvarca, P., 2007. Snow zonation on Hielo Patagónico Sur, Southern Patagonia, derived from Landsat 5 TM data. Global and Planetary Change, 59(1-4), pp.149–158. Aniya, M. et al., 1996. The use of satellite and airborne imagery to invesntory outlet glaciers of the Southern Ratagonia Icefield, South America. Photogrammetric Engineering and Remote Sensing, 62(12), pp.1361–1369. Aoki, Teruo et al., 2000. Effects of snow physical parameters on spectral albedo and bidirectional reflectance of snow surface. Journal of Geophysical Research: Atmospheres, 105(D8), pp.10219–10236. Arendt, A. et al., 2012. Randolph Glacier Inventory [v2.0]: A Dataset of Global Glacier Outlines. Available at: http://www.glims.org/RGI/RGI_Tech_Report_V2.0.pdf [Accessed December 24, 2012]. Arnold, N.S. et al., 2006. Evaluating the potential of high-resolution airborne LiDAR data in glaciology. International Journal of Remote Sensing, 27(5-6), pp.1233–1251. ARSF, 2010. Airborne Thematic Mapper. Available at: http://arsf.nerc.ac.uk/instruments/atm. asp?cookieConsent=A [Accessed March 5, 2013]. Azimuth Systems, 2005. AZGCORR Users Manual. Available at: http://arsf.nerc.ac.uk/documents/ azgcorr_v5.pdf [Accessed February 25, 2013]. Bajracharya, S.R. & Mool, P., 2009. Glaciers, glacial lakes and glacial lake outburst floods in the Mount Everest region, Nepal. Annals of Glaciology, 50(53), pp.81–86. Baker, B.A. et al., 2013. Does spatial resolution matter? A multi-scale comparison of object-based and pixel-based methods for detecting change associated with gas well drilling operations. International Journal of Remote Sensing, 34(5), pp.1633–1651. Baldridge, A.M. et al., 2009. The ASTER Spectral Library Version 2.0. Remote Sensing of Environment, 113, pp.711–715. Baraer, M. et al., 2011. Glacier recession and water resources in Peru’s Cordillera Blanca. Journal of Glaciology, 58(207), pp.134–150. Barcaza, G. et al., 2009. Satellite-derived equilibrium lines in Northern Patagonia Icefield, Chile, and their implications to glacier variations. Arctic, Antarctic, and Alpine Research, 41(2), pp.174–182. Barrand, N.E., James, T.D. & Murray, T., 2010. Spatio-temporal variability in elevation changes of two high-Arctic valley glaciers. Journal of Glaciology, 56(199). Barry, R.G., 2011. The cryosphere – past, present, and future: a review of the frozen water resources of 9. References 110 the world. Polar Geography, 34(4), pp.219–227. Battersby, S.E., Hodgson, M.E. & Wang, J., 2012. Spatial Resolution Imagery Requirements for Identifying Structure Damage in a Hurricane Disaster: A Cognitive Approach. Photogrammetric Engineering and Remote Sensing, 78(6), pp.625–635. Benn, D.I. et al., 2005. Reconstruction of equilibrium-line altitudes for tropical and sub-tropical glaciers. Quaternary International, 138/139, pp.8–21. Benn, D.I. et al., 2012. Response of debris-covered glaciers in the Mount Everest region to recent warming, and implications for outburst flood hazards. Earth-Science Reviews, 114(1-2), pp.156–174. Benson, C.S., 1960. Stratigraphic studies in the snow and firn of the Greenland ice sheet. PhD. Pasadena, California: California Institute of Technology. Available at: http://etd.caltech.edu/etd/ available/etd-03232006-104828/ [Accessed February 25, 2009]. Bindschadler, R. et al., 2011. Getting around Antarctica: new high-resolution mappings of the grounded and freely-floating boundaries of the Antarctic ice sheet created for the International Polar Year. The Cryosphere, 5(3), pp.569–588. Björnsson, H. et al., 1996. The thermal regime of sub-polar glaciers mapped by multi-frequency radio- echo sounding. Journal of Glaciology, 42(140), pp.23–32. Björnsson, H. & Pálsson, F., 2008. Icelandic Glaciers. Jökull, 58, pp.365–386. Björnsson, H., Pálsson, F. & Haraldsson, H. H., 2002. Mass Balance of Vatnajökull (1991-2001) and Langjökull (1996-2001), Iceland. Jökull, 51, pp.75–78. Bolch, T. et al., 2012. The State and Fate of Himalayan Glaciers. Science, 336(6079), pp.310–314. Boresjö Bronge, L. & 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), p.225. Bourdelles, B. & Fily, M., 1993. Snow grain-size determination from Landsat imagery over Terre Adelie, Antarctica. Annals of Glaciology, 17, pp.86–92. Box, J.E. et al., 2012. Greenland ice sheet albedo feedback: thermodynamics and atmospheric drivers. The Cryosphere, 6(4), pp.821–839. Braithwaite, R.J., 1984. Can the mass balance of a glacier be estimated from its equilibrium line altitude? Journal of Glaciology, 30(106), pp.364–368. Braithwaite, R.J., 2002. Glacier mass balance: the first 50 years of international monitoring. Progress in Physical Geography, 26(1), pp.76–95. Braithwaite, R.J. & Müller, F., 1978. On the parameterization of glacier equilibrium line altitude. IAHS- AISH, 126, pp.263–271. Braun, M. et al., 2007. Comparison of remote sensing derived glacier facies maps with distributed mass balance modelling at Engabreen, northern Norway. Glacier Mass Balance Changes and Meltwater Discharge, IAHS 318, pp.126–134. Breiman, L., 2001. Random Forests. Machine Learning, 45(1), pp.5–32. Brenning, A., 2009. Benchmarking classifiers to optimally integrate terrain analysis and multispectral remote sensing in automatic rock glacier detection. Remote Sensing of Environment, 113(1), pp.239–247. Brown, I.A., Kirkbride, M.P. & Vaughan, R.A., 1999. Find the firn line! The suitability of ERS-1 and ERS- 111 9. References 2 SAR - data for the analysis of glacier facies on Icelandic icecaps. International Journal of Remote Sensing, 20(15), p.3217. Brown, I.A., Klingbjer, P. & Dean, A., 2005. Problems with the retrieval of glacier net surface balance from SAR imagery Julian A. Dowdeswell & Ian C Willis, eds. Annals of Glaciology, 42, pp.209–216. Campbell, J.B., 2006. Introduction to Remote Sensing, London: The Guilford Press. Carmagnola, C.M. et al., 2012. Snow spectral albedo at Summit, Greenland: comparison between in situ measurements and numerical simulations using measured physical and chemical properties of the snowpack. The Cryosphere Discussions, 6(6), pp.5119–5167. Casacchia, R. et al., 2002. Field reflectance of snow/ice covers at Terra Nova Bay, Antarctica. International Journal of Remote Sensing, 23(21), pp.4653–4667. Casacchia, R. et al., 2001. Radiometric investigation of different snow covers in Svalbard. Polar Research, 20(1), pp.13–22. Casey, K.A., Kääb, A. & Benn, D.I., 2012. Geochemical characterization of supraglacial debris via in situ and optical remote sensing methods: a case study in Khumbu Himalaya, Nepal. The Cryosphere, 6(1), pp.85–100. Cea, C., Cristobal, J. & Pons, X., 2007. An improved methodology to map snow cover by means of Landsat and MODIS imagery. IEEE Transactions on Geoscience and Remote Sensing, IGARSS 2007, pp.4217–20. Chander, G., Markham, B.L. & Helder, D.L., 2009. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sensing of Environment, 113(5), pp.893–903. Chinn, T.J., Heydenrych, C. & Salinger, M.J., 2005. Use of the ELA as a practical method of monitoring glacier response to climate in New Zealand’s Southern Alps. Journal of Glaciology, 51, pp.85–95. Choudhury, B.J. & Chang, A.T.C., 1981. On the angular variation of solar reflectance of snow. Journal of Geophysical Research, 86(C1), pp.465–472. Cogley, J.G. et al., 2011. Glossary of Glacier Mass Balance, Paris: UNESCO-IHP. Cohen, J., 1960. A Coefficient of Agreement for Nominal Scales. Educational and Psychological Measurement, 20(1), pp.37–46. Colditz, R.R. et al., 2012. Potential effects in multi-resolution post-classification change detection. International Journal of Remote Sensing, 33(20), pp.6426–6445. Comber, A. et al., 2012. Spatial analysis of remote sensing image classification accuracy. Remote Sensing of Environment, 127, pp.237–246. Congalton, R.G., 1991. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sensing of Environment, 37(1), pp.35–46. Congalton, R.G. & Green, K., 1999. Assessing the Accuracy of Remotely Sensed Data: Principles and Practices, Boca Raton, FL: Lewis Publishers. Cuffey, K.M. & Patterson, W.S.B., 2010. The Physics of Glaciers 4th ed., London: Academic Press. Cutler, P.M. & Munro, D.S., 1996. Visible and near-infrared reflectivity during the ablation period on Peyto Glacier, Alberta, Canada. Journal of Glaciology, 42(141), pp.333–340. 9. References 112 Dahlke, H.E. et al., 2012. Contrasting trends in floods for two sub-arctic catchments in northern Sweden – does glacier presence matter? Hydrol. Earth Syst. Sci., 16(7), pp.2123–2141. Demuth, M. & Pietroniro, A., 1999. Inferring Glacier Mass Balance Using RADARSAT: Results from Peyto Glacier, Canada. Geografiska Annaler. Series A, Physical Geography, 81(4), pp.521– 540. Dixon, D.A. et al., 2013. Variations in snow and firn chemistry along US ITASE traverses and the effect of surface glazing. The Cryosphere, 7(2), pp.515–535. Dozier, J. & Painter, T.H., 2004. Multispectral and hyperspectral remote sensing of alpine snow properties. Annual Review of Earth and Planetary Sciences, 32, pp.465–494. Dozier, J., Schneider, S.R. & McGinnis, D.F., 1981. Effect of grain size and snowpack water equivalence on visible and near-infrared satellite observations of snow. Water Resources Research, 17(4), pp.1213–1221. Drusch, M. et al., 2012. Sentinel-2: ESA’s Optical High-Resolution Mission for GMES Operational Services. Remote Sensing of Environment, 120, pp.25–36. Drusch, Matthias, Gascon, Ferran & Berger, M., 2010. GMES Sentinel-2 Mission Requirements Documentq (2.1). Dumont, M. et al., 2010. High-accuracy measurements of snow Bidirectional Reflectance Distribution Function at visible and NIR wavelengths – comparison with modelling results. Atmos. Chem. Phys., 10(5), pp.2507–2520. Dumont, M., Gardelle, J., et al., 2012. Linking glacier annual mass balance and glacier albedo retrieved from MODIS data. The Cryosphere, 6, pp.1527–1539. Dumont, M. et al., 2011. Monitoring spatial and temporal variations of surface albedo on Saint Sorlin Glacier (French Alps) using terrestrial photography. The Cryosphere, 5(3), pp.759–771. Dumont, M., Durand, Y., et al., 2012. Variational assimilation of albedo in a snowpack model and reconstruction of the spatial mass-balance distribution of an alpine glacier. Journal of Glaciology, 58(207), pp.151–164. Dunkle, R.V. & Bevans, J.T., 1956. An Approximate Analysis of the Solar Reflectance and Transmittance of a Snow Cover. Journal of Meteorology, 13, pp.212–216. Dyurgerov, M., 1996. Substitution of long term mass balance data by measurements of one summer. Gletscherkd. Glazial geol., 32, pp.177–184. Dyurgerov, M., Meier, M.F. & Bahr, D.B., 2009. A new index of glacier area change: a tool for glacier monitoring. Journal of Glaciology, 55(192), pp.710–716. Fischer, A., 2011. Comparison of direct and geodetic mass balances on a multi-annual time scale. The Cryosphere, 5(1), pp.107–124. Fleming, K.M., Dowdeswell, J. A. & Oerlemans, J., 1997. Modelling the mass balance of northwest Spitsbergen glaciers and responses to climate change I. M. Whillans, ed. Annals of Glaciology, Vol 24, 1997, 24, pp.203–210. Flowers, G.E. et al., 2007. Glacier fluctuation and inferred climatology of Langjökull ice cap through the Little Ice Age. Quaternary Science Reviews, 26(19-21), pp.2337–2353. Foody, G.M., 2002. Status of land cover classification accuracy assessment. Remote Sensing of Environment, 80(1), pp.185–201. 113 9. References Friedt, J.-M. et al., 2012. Assessing the relevance of digital elevation models to evaluate glacier mass balance: application to Austre Lovénbreen (Spitsbergen, 79N). Polar Record, 48(244), pp.2– 10. Gardner, A.S. & Sharp, M.J., 2010. A review of snow and ice albedo and the development of a new physically based broadband albedo parameterization. Journal of Geophysical Research-Earth Surface, 115, p.F01009. Gerland, S. et al., 1999. Physical and optical properties of snow covering Arctic tundra on Svalbard. Hydrological Processes, 13(1415), pp.2331–2343. Gislason, P.O., Benediktsson, J.A. & Sveinsson, J.R., 2006. Random Forests for land cover classification. Pattern Recognition Letters, 27(4), pp.294–300. Goetz, B., 2012. Making Accurate Field Spectral Reflectance Measurements. Gourmelen, N. et al., 2011. Ice velocity determined using conventional and multiple-aperture InSAR. Earth and Planetary Science Letters, 307(1-2), pp.156–160. Grenfell, T.C., Warren, S.G. & Mullen, P.C., 1994. Reflection of solar radiation by the Antarctic snow surface at ultraviolet, visible, and near-infrared wavelengths. Journal of Geophysical Research, 99(D9), pp.18,669–18,684. Greuell, W. et al., 2007. Assessment of interannual variations in the surface mass balance of 18 Svalbard glaciers from the Moderate Resolution Imaging Spectroradiometer/Terra albedo product. Journal of Geophysical Research, 112, p.D07105. Greuell, W. & Oerlemans, J., 2005a. Assessmenrt of the surface mass balance along the K-transect (Greenland ice sheet) from satellite-derived albedos. Annals of Glaciology, 42, pp.107–116. Greuell, W. & Oerlemans, J., 2004. Narrowband-to-broadband albedo conversion for glacier ice and snow: equations based on modeling and ranges of validity of the equations. Remote Sensing of Environment, 89(1), pp.95–105. Greuell, W. & Oerlemans, J., 2005b. Validation of AVHRR- and MODIS-derived albedos of snow and ice surfaces by means of helicopter measurements. Journal of Glaciology, 51(172), pp.37–48. Greuell, W., Reijmer, C.H. & Oerlemans, J., 2002. Narrowband-to-broadband albedo conversion for glacier ice and snow based on aircraft and near-surface measurements. Remote Sensing of Environment, 82, pp.48–63. Greuell, W. & De Ruyter de Wildt, M., 1999. Anisotropic Reflection by Melting Glacier Ice: Measurements and Parametrizations in Landsat TM Bands 2 and 4. Remote Sensing of Environment, 70(3), pp.265–277. Grinsted, A., 2013. An estimate of global glacier volume. The Cryosphere, 7(1), pp.141–151. Gupta, R.P., Haritashya, U.K. & Singh, P., 2005. Mapping dry/wet snow cover in the Indian Himalayas using IRS multispectral imagery. Remote Sensing of Environment, 97(4), pp.458–469. Hagen, J.O. et al., 2003. On the net mass balance of the glaciers and ice caps in Svalbard, Norwegian Arctic. Arctic Antarctic and Alpine Research, 35(2), pp.264–270. Hall, D.K. et al., 1987. Characterization of snow and ice reflectance zones on glaciers using landsat thematic mapper data. Annals of Glaciology, 9, pp.104–108. Hall, D.K. et al., 1990. Comparison of in situ and satellite-derived reflectances of Forbindels Glacier, Greenland. International Journal of Remote Sensing, 11(3), p.493. 9. References 114 Hall, D.K. et al., 2000. Evaluation of remote-sensing techniques to measure decadal-scale changes of Hofsjökull ice cap, Iceland. Journal of Glaciology, 46, pp.375–388. Hall, D.K. et al., 2009. Evaluation of surface and near-surface melt characteristics on the Greenland ice sheet using MODIS and QuikSCAT data. Journal of Geophysical Research, 114, p.F04006. Hall, D.K., Chang, A.T.. & Siddalingaiah, Honnappa, 1988. Reflectances of glaciers as calculated using Landsat-5 Thematic Mapper data. Remote Sensing of Environment, 25(3), pp.311–312, IN1, 313–321. Hall, D.K., Foster, J.L. & Chang, A.T.., 1992. Reflectance of snow measured in situ and from space in sub- Arctic areas in Canada and Alaska. Geoscience and Remote Sensing, IEEE Transactions on, 30(3), pp.634–637. 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(2), pp.127–140. Hall, J., 2008. The identification of different ice facies on the North Patagonian Icefield (Chile) from ETM+ and ASTER images using a rule-based classification. Journal of Maps Student Edition, pp.11–22. Harrison, S. et al., 2006. A glacial lake outburst flood associated with recent mountain glacier retreat, Patagonian Andes. The Holocene, 16(4), pp.611–620. Heid, T. & Kääb, A., 2012. Repeat optical satellite images reveal widespread and long term decrease in land-terminating glacier speeds. The Cryosphere, 6(2), pp.467–478. Heiskanen, J. et al., 2002. Assessment of Glaciological Parameters Using Landsat Satellite Data in Svartisen, Northern Norway. EARSeL eProceedings, 2, pp.34–42. Helder, D.L. et al., 2012a. Landsat 4 Thematic Mapper Calibration Update. IEEE Transactions on Geoscience and Remote Sensing, 50(6), pp.2400 –2408. Helder, D.L. et al., 2012b. Radiometric Calibration of the Landsat MSS Sensor Series. IEEE Transactions on Geoscience and Remote Sensing, 50(6), pp.2380 –2399. Hendriks, J. & Pellikka, P., 2004. Estimation of Surface Reflectances from Hintereisferner: Spectrometer Measurements and Satellite-Derived Reflectances. Zeitschrift für Gletscherkunde und Glazialgeologie, 38(2), pp.139–154. Hendriks, J. & Pellikka, P., 2008. Semi-automatic glacier delineation from Landsat imagery over Hintereisferner in the Austrian Alps. Zeitschrift für Gletscherkunde und Glazialgeologie, 41, pp.55–75. Hock, R., 2005. Glacier melt: A review of processes and their modelling. Progress in Physical Geography, 29(3), pp.362–391. Hock, R., 2003. Temperature index melt modelling in mountain areas. Journal of Hydrology, 282, pp.104– 115. Hock, R., Kootstra, D. & Reijmer, C.H., 2007. Deriving glacier mass balance from accumulation area ratio on Storglaciaren, Sweden. Glacier Mass Balance Changes and Meltwater Discharge, IAHS Pub 318, pp.163–170. Höfle, B. et al., 2007. Glacier Surface Segmentation Using Airborne Laser Scanning Point Cloud and Intensity Data. International Archives of Photogrammetry and Remote Sensing, 36(3), 115 9. References pp.195–200. Hopkinson, C. & Demuth, M., 2006. Using airborne lidar to assess the influence of glacier downwasting on water resources in the Canadian Rocky Mountains. Canadian Journal of Remote Sensing, 32(2), pp.212–222. Huang, L. et al., 2011. Classification and snow line detection for glacial areas using the polarimetric SAR image. Remote Sensing of Environment, 115(7), pp.1721–1732. Hubbard, B. et al., 2005. Impact of a rock avalanche on a moraine-dammed proglacial lake: Laguna Safuna Alta, Cordillera Blanca, Peru. Earth Surface Processes and Landforms, 30(10), pp.1251–1264. Hulburt, E.O., 1928. The Ultraviolet, Visible and Infrared Reflectivities of Snow, Sand, and Other Substances. Journal of the Optical Society of America, 17(1), p.23. Huss, M. et al., 2012. Conventional versus reference-surface mass balance. Journal of Glaciology, 38(208), pp.278–285. Huss, M. et al., 2013. Towards remote monitoring of sub-seasonal glacier mass balance. Annals of Glaciology, 54(63), pp.85–93. Huss, M. & Farinotti, D., 2012. Distributed ice thickness and volume of all glaciers around the globe. Journal of Geophysical Research, 117, p.F04010. IPCC, 2007. Summary for policymakers. In S. Solomon, ed. Climate change 2007: the physical science basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge: Cambridge University Press. Irons, J.R., Dwyer, J.L. & Barsi, J.A., 2012. The next Landsat satellite: The Landsat Data Continuity Mission. Remote Sensing of Environment, 122, pp.11–21. Jacobsen, D. et al., 2012. Biodiversity under threat in glacier-fed river systems. Nature Climate Change, 2(5), pp.361–364. James, T.D. et al., 2012. Observations of enhanced thinning in the upper reaches of Svalbard glaciers. The Cryosphere, 6(6), pp.1369–1381. Jin, Z. et al., 2008. Snow optical properties for different particle shapes with application to snow grain size retrieval and MODIS/CERES radiance comparison over Antarctica. Remote Sensing of Environment, 112(9), pp.3563–3581. Jóhannesson, T. et al., 2013. Ice-volume changes, bias estimation of mass-balance measurements and changes in subglacial lakes derived by lidar mapping of the surface of Icelandic glaciers. Annals of Glaciology, 54(63), pp.63–74. Joshi, M.D. et al., 1998. Classification of snow facies on the Greenland ice sheet using passive microwave and SAR imagery. In Geoscience and Remote Sensing Symposium Proceedings, 1998. IGARSS ’98. 1998 IEEE International. Geoscience and Remote Sensing Symposium Proceedings, 1998. IGARSS ’98. 1998 IEEE International. pp. 1852–1854 vol.4. Kargel, Jeffrey S. et al., 2005. Multispectral imaging contributions to global land ice measurements from space. Remote Sensing of Environment, 99(1-2), pp.187–219. Keshri, A.K., Shukla, A. & Gupta, R.P., 2009. ASTER ratio indices for supraglacial terrain mapping. International Journal of Remote Sensing, 30(2), pp.519–524. 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 9. References 116 Change, 22(1-4), pp.139–154. Klok, E.J., Greuell, W. & Oerlemans, J., 2003. Temporal and spatial variation of the surface albedo of Morteratschgletscher, Switzerland, as derived from12 Landsat images. Journal of Glaciology, 49(167), pp.491–502. Knap, W.H. et al., 1999a. Comparison of Landsat TM-derived and ground-based albedos of Haut Glacier d’Arolla, Switzerland. International Journal of Remote Sensing, 20(17), p.3293. Knap, W.H. et al. 1999b. Narrowband to broadband conversion of Landsat TM glacier albedos. International Journal of Remote Sensing, 20(10), pp.2091–2110. Knap, W.H. & Reijmer, C.H., 1998. Anisotropy of the Reflected Radiation Field Over Melting Glacier Ice: Measurements in Landsat TM Bands 2 and 4. Remote Sensing of Environment, 65(1), pp.93–104. Kohler, Jack et al., 2007. Acceleration in thinning rate on western Svalbard glaciers. Geophysical Research Letters, 34, p.L18502.Kohler, Jack et al., 2012. MODIS albedo used to detect ELA on Svalbard Glaciers. In AGU Fall Meeting. San Francisco, California, p. C13F–0689. Kokhanovsky, A.A., 2013. Spectral reflectance of solar light from dirty snow: a simple theoretical model and its validation. The Cryosphere Discussions, 7(1), pp.533–550. Kokhanovsky, A.A. & Schreier, M., 2009. The determination of snow specific surface area, albedo, and effective grain size using AATSR space-borne measurements. International Journal of Remote Sensing, 30(4), pp.919–933.Kokhanovsky, A.A. & Zege, E.P., 2004. Scattering Optics of Snow. Applied Optics, 43(7), pp.1589–1602.Koks, M., 2001. Anisotropic reflection of radiation by melting snow, Landsat TM bands 2 and 4. M.Sc. Universiteit Utrecht, Instituut voor Marien en Atmosferisch Onderzoek Utrecht.König, M. et al., 2002. Detection of superimposed ice on the glaciers Kongsvegen and midre Lovenbreen, Svalbard, using SAR satellite imagery. Annals of Glaciology, 34, pp.335–342. König, M. et al., 2004. Two methods for firn-area and mass-balance monitoring of Svalbard glacier with SAR satellite images. Journal of Glaciology, 50(168), pp.116–128.König, M., Winther, J.-G. & Isaksson, E., 2001. Measuring Snow and Glacier Ice Properties from Satellite. Reviews of Geophysics, 39(1), pp.1–27.Kuchiki, K. et al., 2009. Retrieval of snow physical parameters using a ground-based spectral radiometer. Applied Optics, 48(29), pp.5567–5582. Kuchiki, K. et al., 2011. Effect of sastrugi on snow bidirectional reflectance and its application to MODIS data. Journal of Geophysical Research, 116, p.15 PP. Kuhn, M., 1985. Bidirectional Reflectance of Polar and Alpine Snow Surfaces. Annals of Glaciology, 6, pp.164–167.Kulkarni, A.V., 1992. Mass balance of Himalayan glaciers using AAR and ELA methods. Journal of Glaciology, 38(128), pp.101–104. Laffly, D. et al., 2012. High temporal resolution monitoring of snow cover using oblique view ground-based pictures. Polar Record, 48(244), pp.11–16. 117 9. ReferencesLangley, K. et al., 2008. From glacier facies to SAR backscatter zones via GPR. IEEE Transactions on Geoscience and Remote Sensing, 46(9), pp.2506–2516.Li, W. et al., 2001. Snow grain size retrieved from near-infrared radiances at multiple wavelengths. Geophysical Research Letters, 28(9), p.PAGES 1699–1702.Li, Z. et al., 2012. Glacier Snow Line Detection on a Polarimetric SAR Image. IEEE Geoscience and Remote Sensing Letters, 9(4), pp.584 –588. Lie, Ø., Dahl, S.O. & Nesje, A., 2003. A theoretical approach to glacier equilibrium-line altitudes using meteorological data and glacier mass-balance records from southern Norway. Holocene, 13(3), pp.365–372. Loveland, T.R. & Dwyer, J.L., 2012. Landsat: Building a strong future. Remote Sensing of Environment, 122, pp.22–29.Lucht, W., Schaaf, C.B. & Strahler, A.H., 2000. An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Transactions on Geoscience and Remote Sensing, 38(2), pp.977 –998.Luckman, A. et al., 2012. Basal crevasses in Larsen C Ice Shelf and implications for their global abundance. The Cryosphere, 6(1), pp.113–123.Lutz, E., Geist, Th. & Stötter, J., 2003. Investigations of airborne laser scanning signal intensity on glacial surfaces – Utilizing comprehensive laser geometry modelling and orthophoto surface modelling (A case study: Svartisheibreen. In Proceedings, ISPRS Workshop on 3-D reconstruction from airborne laserscanner and INSAR data.MacArthur, A., 2007a. Field Guide for the ASF FieldSpec Pro - Radiance/Irradiance Measurements in Raw DN Mode. Available at: http://fsf.nerc.ac.uk/resources/guides/pdf_guides/ asd_guide_v2_RadIrrad.pdf [Accessed February 25, 2013]. MacArthur, A., 2007b. Field Guide for the ASF FieldSpec Pro - Raw DN Mode. Available at: http://fsf. nerc.ac.uk/resources/guides/pdf_guides/asd_guide_v2_dn.pdf [Accessed February 25, 2013]. MacArthur, A., 2007c. Field Guide for the ASF FieldSpec Pro - White Reference Mode. Available at: http://fsf.nerc.ac.uk/resources/guides/pdf_guides/asd_guide_v2_wr.pdf [Accessed February 25, 2013].MacArthur, A., MacLellan, C.J. & Malthus, T., 2012. The Fields of View and Directional Response Functions of Two Field Spectroradiometers. IEEE Transactions on Geoscience and Remote Sensing, 50(10), pp.3892 –3907.Machguth, H. et al., 2009. Calculating distributed glacier mass balance for the Swiss Alps from regional climate model output: A methodical description and interpretation of the results. Journal of Geophysical Research, 114(D19), p.D19106. Markham, B.L. et al., 2005. SLC-off Landsat-7 ETM+ reflective band radiometric calibration. In Earth Observing Systems X. Earth Observing Systems X. San Diego, CA, USA: SPIE, p. 58820D–11. Available at: http://link.aip.org/link/?PSI/5882/58820D/1 [Accessed February 16, 2009]. Markham, B.L. & Helder, D.L., 2012. Forty-year calibrated record of earth-reflected radiance from 9. References 118 Landsat: A review Keywords: Landsat, Radiometric Calibration, History. MSS, TM, ETM+. Remote Sensing of Environment, 122, pp.30–40. Masek, J.G. et al., 2006. A Landsat surface reflectance dataset for North America, 1990-2000. IEEE Geoscience and Remote Sensing Letters, 3(1), pp.68 – 72.McFadden, E.M., Ramage, J. & Rodbell, D.T., 2011. Landsat TM and ETM+ derived snowline altitudes in the Cordillera Huayhuash and Cordillera Raura, Peru, 1986–2005. The Cryosphere, 5(2), pp.419–430.Meier, M.F. et al., 2007. Glaciers dominate eustatic sea-level rise in the 21st century. Science, 317, pp.1064–1067.Michishita, R. et al., 2012. Bi-scale analysis of multitemporal land cover fractions for wetland vegetation mapping. ISPRS Journal of Photogrammetry and Remote Sensing, 72, pp.1–15. Milton, E.J. et al., 2009. Progress in field spectroscopy. Remote Sensing of Environment, 113(Supplement 1), pp.S92–S109.Monserud, R. & Leemans, R., 1992. Comparing Global Vegetation Maps with the Kappa-Statistic. Ecological Modelling, 62(4), pp.275–293. Montpetit, B. et al., 2012. New shortwave infrared albedo measurements for snow specific surface area retrieval. Journal of Glaciology, 58(211), pp.941–952. Nakamura, T. et al., 2001. Spectral reflectance of snow with a known grain-size distribution in successive metamorphism. Cold Regions Science and Technology, 32(1), pp.13–26. NASA, 2002. ASTER User Handbook, Version 2. Available at: http://asterweb.jpl.nasa.gov/ content/03_data/04_Documents/aster_user_guide_v2.pdf [Accessed March 5, 2013]. NASA, 2011. Landsat 7 Science Data Users Handbook. Available at: http://landsathandbook.gsfc. nasa.gov/ [Accessed July 1, 2011]. NASA, 2013. MODIS Calbiration Parameters. Available at: http://mcst.gsfc.nasa.gov/calibration/ parameters [Accessed March 5, 2013]. NERC FSF, 2013. Excel Post Processing Templates and User Guides. Available at: http://fsf.nerc. ac.uk/resources/post-processing/index.shtml [Accessed February 25, 2013].Nolin, A.W. & Dozier, J., 1993. Estimating snow grain size using AVIRIS data. Remote Sensing of Environment, 44(2-3), pp.231–238. Nolin, A.W. & Payne, M.C., 2007. Classification of glacier zones in western Greenland using albedo and surface roughness from the Multi-angle Imaging SpectroRadiometer (MISR). Remote Sensing of Environment, 107(1-2), pp.264–275. Notarnicola, C. et al., 2013a. Snow Cover Maps from MODIS Images at 250 m Resolution, Part 1: Algorithm Description. Remote Sensing, 5, pp.110–126. Notarnicola, C. et al., 2013b. Snow Cover Maps from MODIS Images at 250 m Resolution, Part 2: Validation. Remote Sensing, 5(4), pp.1568–1587. Nuth, C. et al., 2012. Estimating the long-term calving flux of Kronebreen, Svalbard, from geodetic elevation changes and mass-balance modelling. Journal of Glaciology, 58(207), pp.119–133. 119 9. References O’Brien, H.W. & Munis, R.H., 1975. Red and Near-Infrared Reflectance of Snow. CRREL Research Reports, 332.Oerlemans, J., 1994. Quantifying Global Warming from the Retreat of Glaciers. Science, 264(5156), pp.243–245. Ohmura, A., Kasser, P. & Funk, M., 1992. Climate at the equilibirum line of glaciers. Journal of Glaciology, 38(130), pp.397–411.Orheim, O. & Lucchitta, B., 1987. Snow and ice studies by thematic mapper and multispectral scanner landsat images. Annals of Glaciology, 9, pp.109–118.Painter, T.H., Bryant, A.C. & Skiles, S.M., 2012. Radiative forcing by light absorbing impurities in snow from MODIS surface reflectance data. Geophysical Research Letters, 39(17), p.L17502. Painter, T.H., 2011. Comment on Singh an others, “Hyperspectral analysis of snow reflectance to understand the effects of contamination and grain size”. Journal of Glaciology, 57(201), pp.183–185.Painter, T.H. et al., 2007. Contact spectroscopy for determination of stratigraphy of snow optical grain size. Journal of Glaciology, 53, pp.121–127.Palmer, S. et al., 2009. Ice velocity measurements of Langjokull, Iceland, from interferometric synthetic aperture radar (InSAR). Journal of Glaciology, 193(55), pp.834–838.Pálsson, F. et al., 2012. Mass and volume changes of Langjökull ice cap, Iceland, ⬚1890 to 2009, deduced from old maps, satellite images and in situ mass balance measurements. Jökull, 62, pp.81–95.Palubinskas, G., Muller, R. & Reinartz, P., 2003. Mosaicking of optical remote sensing imagery. In Geoscience and Remote Sensing Symposium, 2003. IGARSS ’03. Proceedings. 2003 IEEE International. Geoscience and Remote Sensing Symposium, 2003. IGARSS ’03. Proceedings. 2003 IEEE International. pp. 3955–3957 vol.6. Available at: 10.1109/ IGARSS.2003.1295326 [Accessed July 16, 2010].Paul, F. et al., 2002. Comparison of TM-Derived Glacier Areas with Higher Resolution Data Sets. EARSeL eProceedings, 2, pp.15–21.Paul, F., 2000. Evaluation of different methods for glacier mapping using Landsat TM. In EARSeL- SIG-Workshop Land Ice and Snow. Desden/FRG, pp. 238–245.Paul, F. et al., 2013. On the accuracy of glacier outlines derived from remote sensing data. Annals of Glaciology, 54(63), pp.171–182.Paul, F. & Kääb, A., 2005. Perspectives on the production of a glacier inventory from multispectral satellite data in Arctic Canada: Cumberland Peninsula, Baffin Island. Annals of Glaciology, 42, pp.59–66.Pellikka, P. & Rees, W.G., 2010. Remote Sensing of Glaciers, London: CRC Press.Pelto, M.S., 2010. Forecasting temperate alpine glacier survival from accumulation zone observations. The Cryosphere, 4(1), pp.67–75.Phinn, S.R., Roelfsema, C.M. & Mumby, P.J., 2012. Multi-scale, object-based image analysis for mapping geomorphic and ecological zones on coral reefs. International Journal of 9. References 120 Remote Sensing, 33(12), pp.3768–3797.Picard, G. et al., 2012. Inhibition of the positive snow-albedo feedback by precipitation in interior Antarctica. Nature Climate Change, 2(11), pp.795–798.Pope, A. et al., 2013. Combining airborne lidar and Landsat ETM+ data with photoclinometry to produce a digital elevation model for Langjökull, Iceland. International Journal of Remote Sensing, 34(4), pp.1005–1025. Rabatel, Antoine et al., 2008. 25 years (1981-2005) of equilibrium-line altitude and mass-balance reconstruction on Glacier Blanc, French Alps, using remote-sensing methods and meteorological data. Journal of Glaciology, 54, pp.307–314.Rabatel, Antoine, Dedieu, J.-P. & Vincent, C., 2005. Using remote-sensing data to determine equilibrium-line altitude and mass-balance time series: validation on three French glaciers, 19942002. Journal of Glaciology, 51, pp.539–546.Racoviteanu, A. & Williams, M.W., 2012. Decision Tree and Texture Analysis for Mapping Debris-Covered Glaciers in the Kangchenjunga Area, Eastern Himalaya. Remote Sensing, 4(12), pp.3078–3109. Radić, V. & Hock, R., 2010. Regional and global volumes of glaciers derived from statistical upscaling of glacier inventory data. Journal of Geophysical Research, 115, p.F01010. Radić, V. & Hock, R., 2011. Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise. Nature Geosci, 4(2), pp.91–94. Ranisavljević, E. et al., In Review. A dynamic and generic cloud computing model for environmental analysis using in situ data applied to glacier mass balance. Journal of Applied Earth Observation and Geoinformation.Rasmussen, L.A. & Kohler, Jack, 2007. Mass balance of three Svalbard glaciers reconstructed back to 1948. Polar Research, 26(2), pp.168–174.Rees, W.G., 2008. Comparing the spatial content of thematic maps. International Journal of Remote Sensing, 29(13), pp.3833–3844.Rees, W.G., 2006. Remote Sensing of Snow and Ice, Boca Raton: Taylor & Francis.Rees, W.G. & Arnold, N.S., 2007. Mass Balance and Dynamics of a Valley Glacier Measured by High-Resolution LiDAR. Polar Record, 43(04), pp.311–319. Rees, W.G. & Arnold, N.S., 2006. Scale-dependent roughness of a glacier surface: implications for radar backscatter and aerodynamic roughness modelling. Journal of Glaciology, 52(177), pp.214–222. Reijmer, C.H., Knap, W.H. & Oerlemans, J., 1999. The Surface Albedo of Vatnajokull Ice Cap, Iceland: A Comparison Between Satellite-Derived and Ground-Based Measurements. Boundary- Layer Meteorology, 92, pp.125–144. Reynolds, J.M., 1998. High-altitude glacial lake hazard assessment and mitigation: a Himalayan perspective. Geological Society, London, Engineering Geology Special Publications, 15(1), pp.25–34.Rippin, D. et al., 2003. Changes in geometry and subglacial drainage of Midre Lovenbreen, Svalbard, determined from digital elevation models. Earth Surface Processes and Landforms, 121 9. References28(3), pp.273–298. Robinson, I. & MacArthur, A., 2011. The FSF Post Processing Toolbox User Guide. Available at: http://fsf.nerc.ac.uk/user_group/user_group.shtml [Accessed May 30, 2011]. Roe, G., 2011. What do glaciers tell us about climate variability and climate change? Journal of Glaciology, 57(203), pp.567–577. Rott, H. & Søgaard, H., 1987. Spectral Reflectance of Snow-Covered and Snow-Free Terrain in Western Greenland. Zeitschrift für Gletscherkunde und Glazialgeologie, 23(2), pp.115–121.De Ruyter De Wildt, M.S., Oerlemans, J. & Björnsson, H., 2003. A calibrated mass balance model from Vatnajökull, Iceland. Jökull, 52, pp.1–20.De Ruyter De Wildt, M.S., Oerlemans, J. & Björnsson, H., 2002. A method for monitoring glacier mass balance using satellite albedo measurements: application to Vatnajökull, Iceland. Journal of Glaciology, 48(161), pp.267–278.Rye, C.J. et al., 2010. Modeling the surface mass balance of a high Arctic glacier using the ERA-40 reanalysis. Journal of Geophysical Research, 115(F2), p.F02014.Scambos, T.A. et al., 2012. Extent of low-accumulation “wind glaze” areas on the East Antarctic plateau: implications for continental ice mass balance. Journal of Glaciology, 58(210), pp.633–647.Scambos, T.A. & Haran, T., 2002. An image-enhanced DEM of the Greenland ice sheet J.-G. Winther & R. Solberg, eds. Annals of Glaciology, 34, pp.291–298. Schaepman-Strub, G. et al., 2006. Reflectance quantities in optical remote sensing - definitions and case studies. Remote Sensing of Environment, 103(1), pp.27–42.Shea, J. M., B. Menounos, R. D. Moore, and C. Tennant, 2013. An Approach to Derive Regional Snow Lines and Glacier Mass Change from MODIS Imagery, Western North America. The Cryosphere, 7(2), pp.667–680. doi:10.5194/tc-7-667-2013. Shukla, A., Gupta, R.P. & Arora, M., 2009. Estimation of debris cover and its temporal variation using optical satellite sensor data: a case study in Chenab basin, Himalaya. Journal of Glaciology, 55(191), pp.444–452. Sidjak, R.W. & Wheate, R.D., 1999. Glacier mapping of the Illecillewaet icefield, British Columbia, Canada, using Landsat TM and digital elevation data. International Journal of Remote Sensing, 20(2), pp.273–284. Sigurðsson, O., 1994. Comments on “Analysis of glacier faciers using satellite techniques” by Williams and others. Journal of Glaciology, 40(134), pp.202–203.Sigurðsson, O., 1998. Glacier variations in Iceland 1930-1995. Jökull, 45, pp.3–25.Sobrino, J.A. et al., 2012. Impact of spatial resolution and satellite overpass time on evaluation of the surface urban heat island effects. Remote Sensing of Environment, 117, pp.50–56.Spencer, N.H., 2003. Investigating Data with Andrews Plots. Social Science Computer Review, 21(2), pp.244–249. Stehman, S.V. & Wickham, J.D., 2011. Pixels, blocks of pixels, and polygons: Choosing a spatial unit for thematic accuracy assessment. Remote Sensing of Environment, 115(12), pp.3044– 9. References 122 3055.Stroeve, J.C., Box, J.E. & Haran, T., 2006. Evaluation of the MODIS (MOD10A1) daily snow albedo product over the Greenland ice sheet. Remote Sensing of Environment, 105, pp.155–171. Takeuchi, N., 2009. Temporal and spatial variations in spectral reflectance and characteristics of surface dust on Gulkana Glacier, Alaska Range. Journal of Glaciology, 55(192), pp.701–709.Tape, K. et al., 2010. Recording microscale variations in snowpack layering using near-infrared photography. Journal of Glaciology, 56(195), pp.75–80.Thuillier, G. et al., 2003. The Solar Spectral Irradiance from 200 to 2400 <i>nm</i> as Measured by the SOLSPEC Spectrometer from the Atlas and Eureca Missions. Solar Physics, 214(1), pp.1–22.Warren, S.G., 1982. Optical Properties of Snow. Reviews of Geophysics and Space Physics, 20(1), pp.67–89. Warren, S.G., Brandt, R.E. & Boime, R.D., 1993. Blue ice and green ice. Antarctic Journal, pp.255–256. Warren, S.G. & Wiscombe, W.J., 1980. A Model for the Spectra Albedo of Snow. II: Snow Containing Atmospheric Aerosols. Journal of the Atmospheric Sciences, 37, pp.2734–2747. Williams, R.S., Hall, D.K. & Benson, C.S., 1991. Analysis of glacier facies using satellite techniques. Journal of Glaciology, 37(125), pp.120–128. Wilson, M.J. & Oreopoulos, L., 2013. Enhancing a Simple MODIS Cloud Mask Algorithm for the Landsat Data Continuity Mission. IEEE Transactions on Geoscience and Remote Sensing, 51(2), pp.723–731. Winther, J.-G., 1993. Landsat TM derived and in situ summer reflectance of glaciers in Svalbard. Polar Research, 12(1), pp.37–55. Winther, J.-G. et al., 1999. Spectral reflectance of melting snow in a high Arctic watershed on Svalbard: some implications for optical satellite remote sensing studies. Hydrological Processes, 13, pp.2033–2049. Wiscombe, W.J. & Warren, S.G., 1980. A Model for the Spectral Albedo of Snow. I: Pure Snow. Journal of the Atmospheric Sciences, 37(12), pp.2712–2733. Wolken, G.J., Sharp, M.J. & Wang, L., 2009. Snow and ice facies variability and ice layer formation on Canadian Arctic ice caps, 1999–2005. Journal of Geophysical Research, 114, p.F03011. De Woul, M. et al., 2006. Firn layer impact on glacial runoff: a case study at Hofsjökull, Iceland. Hydrological Processes, 20, pp.2171–2185. Wulder, M.A. et al., 2012. Opening the archive: How free data has enabled the science and monitoring promise of Landsat. Remote Sensing of Environment, 122(0), pp.2–10. Zelazowski, P. et al., 2011. Reconciling satellite-derived atmospheric properties with fine-resolution land imagery: Insights for atmospheric correction. Journal of Geophysical Research, 116(D18), p.D18308. Zemp, M., Hoelzle, M. & Haeberli, W., 2009. Six decades of glacier mass-balance observations: a review of the worldwide monitoring network. Annals of Glaciology, 50, pp.101–111. Zeng, Q. et al., 1983. A study of spectral reflection characteristics for snow, ice and water in the north of China. Hydrological Applications of Remote Sensing and Remote Data Transmissions, 145, pp.451–462. 123 Fjord the Reindeer 124