The Potential of Earth Observation for the Analysis of Cold Region Land Surface Dynamics in Europe—A Review

Abstract

Cold regions affect global, regional and local climate; oftentimes they are relevant for water supply, host valuable ecosystems, and support human livelihood. They are thus eminently important for human society. In the context of ongoing climate change, monitoring and understanding cold region land surface dynamics is essential for environmental scientists, stakeholders and decision makers. However, the definition of cold regions remains inexplicit, and no up-to-date cold region maps or overarching spatial analyses exist. For example, Europe has densely populated cold regions, but hardly an article exists that provides a solid overview of Earth Observation (EO) based applications assessing cold region land surface dynamics in Europe. With this review article we aim at closing this gap by providing an overview of EO-based techniques for cold region observation in Europe, focusing on the dynamics of glaciers and snow. We present a novel spatial delineation of cold regions for Europe before analyzing the benefits and limitations of different EO sensor types and data processing methods for EO based cold region research. Furthermore, we identify research gaps and discuss challenges for future studies.

1. Introduction

Worldwide, cold regions are undergoing rapid change. Amongst those changes, one of the most obvious is the near-global glacier retreat over the last century [1,2,3]. Besides, changes to ice, snow, and frozen ground (e.g., sea ice shrinking, earlier snowmelt onset, and permafrost thawing) provide obvious evidence of a changing climate. Since cold regions are sensitive to global warming, they can provide indications and increase comprehension of climate fluctuations. Thus, obtaining continuous, timely and reliable information about cold region land surface dynamics (e.g., geographical extent and characteristics) is an essential task for supporting management activities as well as policy- and decision-making processes. Many Earth Observation (EO) based research groups such as the Group on Earth Observations Cold Regions Initiative (GEOCRI) have been established. They have carried out projects around the globe and produced geospatial information products on subjects relevant for cold regions. They include the Interactive Multisensor Snow and Ice Mapping System (IMS) provided by the National Oceanic and Atmospheric Administration/National Environmental Satellite, Data, and Information Service (NOAA/NESDIS) [4,5], the Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover products (MOD10/MYD10) [6], the Global Land Ice Measurements from Space (GLMS) project [7,8], the global inventory of glacier outlines from the Randolph Glacier Inventory (RGI) [9], the European Space Agency’s (ESA) GlobSnow product coordinated by the Finnish Meteorological Institute (FMI) [10,11], the GlobPermafrost initiative (2016–2019) launched by ESA [12], the Global Record of Daily Landscape Freeze/Thaw Status (FT-ESDR) from the National Aeronautics and Space Administration’s (NASA) Making Earth System Data Records for Use in Research Environments (MEaSUREs) program [13], the German Aerospace Center/German Remote Sensing Data Center’s (DLR/DFD) Global SnowPack data set [14], the Global Boreal Forest Mapping (GBFM) project led by the Earth Observation Research Center (EORC) of the Japan Aerospace Exploration Agency (JAXA) [15], the HiMAC (High Mountain and Cold Regions Using big Earth Data) service of the Digital Belt and Road Initiative (DBAR) [16], and the “Special Report on Climate Change and Oceans and the Cryosphere” to be published in the future by the Intergovernmental Panel on Climate Change (IPCC).

It is noteworthy that the impacts of climate change differ significantly among regions [17]. Northern Europe is one of the northernmost populated areas, where many settlements have more than 100,000 inhabitants. Because of the presence of ice and snow, high latitude regions and high-mountain areas like the European Alps are also one of the most climate sensitive regions in Europe [18]. In recent years, Europe as a whole has experienced great climate variability and extremes [19,20,21], such as continuously increasing average temperatures particularly in northern Europe [17,22], more frequent high-temperature extremes and less frequent low-temperature extremes [22,23,24], increasing annual precipitation in northern Europe and decreasing precipitation in some areas of southern Europe [24], and increasing mean sea level, excluding the northern Baltic Sea [24,25,26,27]. Thereby, significant changes in land, freshwater, and marine ecosystems have occurred in Europe [28]. For example, advancement in the mean onset dates of spring events (e.g., bud burst, breaking hibernation, flowering, migration, snowmelt) has affected much of Europe’s flora and fauna [29,30,31,32]. Annual maximum water discharge from melting has increased in parts of northwestern Europe [33,34], and the water quality of the River Meuse (western Europe) has increasingly deteriorated due to more frequent and intense droughts induced by climate change [35]. In such a context, evidence-based adaptation strategies are urgently required. To contribute to the development of such strategies, cold region land surface dynamics must be monitored in Europe not only to determine local vulnerability, but also to provide a mesoscale assessment of climate change.

1.1. Definitions of Cold Region

“Cold region” is a term that is often used by environmental scientists, engineers, physicists and other scientists. It was first used by the German scientist Wladimir Köppen (1846–1940), who spoke of a “Kalter Gürtel” (“cold belt”) in 1884 [36]. Cold regions are defined either descriptively or quantitatively. A typical descriptive definition of cold regions based on geographical extent was provided by GEOCRI. According to GEOCRI, cold regions include polar regions and high mountain areas, namely, the Arctic, Antarctic, high-latitude oceans, the Himalaya Third Pole, and mountainous cold areas [37]. At the regional scale of Europe, cold regions were defined by Heal [38] as the European Arctic and Alpine areas. In China, Li and Liu [39] identified the northern area of the Qinling Mountains and the Huaihe River as cold regions.

Air temperature can be used to delineate a cold region quantitatively. Köppen [36] characterized cold regions as areas with 1–4 moderate months (i.e., mean temperature between 10 °C and 20 °C) and 11–8 cold months (i.e., mean temperature less than 10 °C). Chen et al. [40] mapped cold regions in China based on the mean temperature of the coldest month (below 3 °C), fewer than 5 months with mean temperature above 10 °C, and mean annual temperature below 5 °C. In addition to air temperature, Yang et al. [41] located the cold regions in China via thresholding 10 different climate factors. To identify the southern boundary of cold regions, a 30-cm frost penetration depth occurring once in 10 years is a generally accepted criterion [42]. To date, the most comprehensive and widely-acknowledged definition was provided by the Cold Regions Research and Engineering Laboratory (CRREL) [43], which located the cold regions in the Northern Hemisphere according to air temperature (0 °C and −18 °C isotherms), snow depth (30 and 61 cm isolines), ice cover (100 and 180 annual mean unnavigable days), and frozen ground extent (permafrost and 30 cm frost penetration).

1.2. Characterstics of Cold Regions and Their Significance

Cold regions include a variety of land-cover types (e.g., snow, glacier, permafrost, boreal forest, tundra) and living organisms existing in a cryospheric environment. Frozen water (in its various forms) is one of the cryosphere’s essential characteristics. It includes snow, glaciers and ice caps, ice shelves and ice sheets, sea ice, lake and river ice, frozen ground and permafrost. These elements play a pivotal role in Earth’s climate system by altering the surface energy budget, the water cycle, primary productivity, surface gas exchange and sea level [3]. In general, cold region characteristics are not only passive indicators reflecting climate change, but also significant and long-lasting factors influencing physical, biological and social systems throughout vast parts of Earth’s surface [3].

Snow is the most widespread cryospheric component during winter, with a maximum extent of about 45.2 million km² on the land surface in the Northern Hemisphere [44]. Snow cover influences Earth’s energy balance, river discharge, and the hydropower potential; it can be an environmental hazard and is a factor in tourism [45,46,47]. Frozen ground is the second most extensive component in cold regions. Its active layer of periodic melting plays a key role in hydrological and geomorphological regimes, which can therefore have a drastic impact on engineering projects [42]. Permafrost thaw has led to severe consequences such as carbon/methane emission [48,49], which may further exacerbate global warming. Other cryosphere components are ice sheets, ice caps and glaciers, which are primary reservoirs affecting not only water supply but also sea level [44,50,51]. Biospheric components also play an important role in cold regions. The boreal forest is the largest terrestrial biome [52]. Together with the tundra, boreal forests are essential habitats for diverse flora and fauna. The significance of boreal forest and tundra is linked to their biodiversity [53], atmospheric pollution uptake [54,55], organic nitrogen mineralization [52], effects on regional and global climate [56], and their economic value [57,58,59] for countries such as Russia, Canada, and Sweden.

1.3. Objectives and Scope of This Paper

An overview of EO-based research on terrestrial cold region land surface dynamics in Europe is given below, with a special focus on glacier and snow dynamics. In this context we present a novel spatial delineation of cold regions for Europe before analyzing the benefits and limitations of different EO sensor types and data processing methods for EO based cold region research. Furthermore, we identify research gaps and discuss challenges for future study.

So far, there is no up-to-date cold region map for Europe. Therefore, we undertook a quantitative analysis of 30 years (1986–2015) of reanalysis data to delineate cold region boundary dynamics in Europe. In addition, we explored the potential of EO for mapping cold regions (Section 2). Thereafter, we looked at relevant characteristics in the delineated cold regions in Europe. A particular focus was on glaciers (Section 3.1) and snow (Section 3.2), which are the most widely distributed cold region components in Europe. Study areas of the reviewed articles included single countries, continents, and the Northern Hemisphere. Although Greenland belongs to Europe (Denmark) politically, geographically it is actually part of the North American continent. In this review paper focusing on the cold region land surface dynamics of Europe, a large number of publications on Greenland have accordingly been excluded. We then discuss the current status of EO monitoring applications in cold regions of Europe as well as the corresponding advantages and disadvantages of different types of EO data. We close with a discussion including the identification of current research gaps and needs (Section 4).

2. The Potential of Earth Observation for Cold Region Delineation

We presented selected definitions of cold regions in Section 1.1, which differ according to the objectives of the research. Engineers usually prefer a definition according to frost penetration depth due to its significance for construction work. For environmental scientists, air-temperature-based definitions are more useful. However, a comprehensive (i.e., multiple parameter based) delineation integrating air temperature, snow depth, and frost penetration depth may be beneficial as it provides additional information on atmosphere-land interaction.

Moreover, cold regions have highly temporal dynamics in extent and condition, which deserves particular attention in the context of ongoing global climate change. For the purpose of cold region extent delineation, temporally continuous data of high to medium resolution is preferred, and must be analyzed in a quantitative way. For more than 40 years, EO sensors have collected a tremendous amount of data at different spatiotemporal resolution. This makes EO a capable candidate for cold region mapping and monitoring. Currently, no cold region map of Europe based on an analysis of EO data collected over a long term exists, since precisely extracting aforementioned parameters (i.e., air temperature, snow depth and frozen ground cover) solely from EO data remains challenging. Also, the lack of an unambiguous definition of cold region hinders the mapping procedure. Thus, one of main objective of this study is to explore the usefulness of EO data for delineating European cold regions in a quantitative way.

2.1. Multi-Temporal Cold Region Mapping

To achieve the aim to assess cold region land surface dynamics, a “reference map” is initially needed. Most of the criteria mentioned in Section 1.2 [36,40,41,42,43] emphasize intra-annual variation, but neglect the inter-annual changes. The ERA-Interim Archive from the European Centre for Medium-Range Weather Forecasts (ECMWF) [60] is a continuously updated global reanalysis of data set collected since 1979. In consideration of CRREL’s definition [43], we used the ERA-Interim air temperature and the snow depth data (based on land stations, ships, and drifting buoys observations) span from 1986 to 2015 to produce a reference map of the extent of the European cold region. Afterwards, we compared the result with mean Snow Cover Duration (SCD) during 2000–2014 derived from Global SnowPack data set [14] to explore the potential of EO in cold region delineation.

We did not consider the criterion “ice cover” for two reasons: one is the lack of spatiotemporally continuous data; the other is its questionable usefulness, also noted by CRREL [43]. Accordingly, we use the following criteria:

(1)

Air temperature below 0 °C observed half of the time during the coldest month of the year;

(2)

Maximum observed snow depth on the ground no less than 30 cm;

(3)

30 cm frost penetration into the ground.

It has to be mentioned that the required snow depth data should be actual accumulation but not water equivalent according to CRREL [43]. In the ERA-Interim Archive the unit of snow depth is meters of water equivalent. Thus, a conversion from snow depth in water equivalent to actual accumulation was conducted empirically according to the equations proposed by Chang et al. [61]. Moreover, the extent of 30 cm frost penetration is approximated by calculating 100 degree days of freezing temperature (based on 0 °C) as originally proposed by CRREL [43]. Next, to combine the results from the three above-mentioned criteria, we introduce a stability metric (S_CR $\in$ [0,1]), that is, a normalized index of the votes from m binary classifier (criterion) over n years that represents the frequency of a certain area being classified as a cold region over the past n year:

$$S_{CR}= \frac{1}{mn} {\displaystyle \sum }_{i=1}^m{\displaystyle \sum }_{j=1}^n{V}_{i,j}$$

(1)

where V_i,j denotes the binary result from the ith classifier at jth year. In this article, m = 3 (i.e., three considered criteria) and n = 30 (i.e., a 30 year time span). The result is shown in Figure 1. Different results were acquired based on the criterion selected for determining cold region extent. Apart from Southern Europe, most of Europe was classified as a cold region according to air temperature thresholding. Based on 100 freezing degree days, only northernmost Europe was considered to be a cold region. In comparison, the geographical patterns are visible either based on snow depth or the stability product. The cold region stability patterns were highly consistent with most descriptive definitions (e.g., [37,38]). Geographically, the Carpathian Mountains, European Alps, Finland, Iceland, Norway (including the Svalbard Archipelago) and Sweden are covered.

2.2. Potential of Snow Cover Duration as an Indicator for Cold Region Mapping

Due to the difficulty of precisely retrieving the aforementioned parameters (i.e., air temperature and snow depth) only using EO data, alternative indicators are needed. Given the high climatic sensitivity as well as mature existing mapping methods, temporal snow cover parameters derived from EO data hold great potential for delineating a changing cold region boundary. Furthermore, an EO-derived cold region boundary would improve the spatial resolution (e.g., 500 m using MODIS data) resulted from the reanalysis data (i.e., 80 km from ERA-Interim). Among the temporal snow cover parameters, we explored the suitability of 15-year mean SCD derived from Global SnowPack data set [14] (Figure 2) from DLR/DFD. The DLR/DFD Global SnowPack is a cloud-free time series of daily snow cover processed from the operational 500 m MODIS daily snow cover products MOD10A1 and MYD10A1 [6], and has been validated using European Climate Assessment & Dataset (ECA&D) station data in Europe (details see [14]). To ensure the consistency of spatial resolution, we downscaled the spatial resolution of mean SCD data from 500 m to 80 km (the spatial resolution of the ERA-Interim data set and the cold region stability result). We performed a correlation analysis between mean SCD (2000–2014) and cold region stability. The results indicate high correlation (R² = 0.734 with a p-value < 0.001).

3. The Potential of EO for the Analysis of Cold Region Characteristics

In cold regions of Europe, long-term monitoring of environmental and social indicators is required [17]. Conventional in-situ measurements are time-consuming, cost intensive, logistically demanding, and challenging, particularly in cold regions with low accessibility [62]. Hence, ground data are often spatiotemporally discontinuous, i.e., not up-to-date, sparsely distributed, and they do not cover more than a decade [63,64]. To address such problems, satellite based EO data provide a promising alternative to in-situ measurements. EO sensors have collected a tremendous amount of data at different temporal resolutions. Thanks to the free and open data policies of archives for long-term missions, such as the Advanced Very High Resolution Radiometer (AVHRR), MODIS, Landsat, Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), and the novel European Sentinel fleet, large volumes of satellite data are now accessible at no charge [65,66]. With satellite remote sensing we can continuously observe and monitor Earth’s surface and ensure objective, repeated, synoptic, consistent and large-scale information for land surface dynamics in cold regions.

3.1. Glaciers

The history of systematic large-scale glacier monitoring reaches back to 1894, when the International Glacier Commission was established in Zurich, Switzerland, at the 6th International Geological Congress. So far, abundant long-term records of glaciological observation exist for Europe, including photography, paintings and other graphics. Over time, the records allow glaciologists to trace back glacier changes to the 19th century [67,68]. Furthermore, these long-term glaciological records are providing valuable information for evaluating the response of glaciers to climate change [69].

Nowadays, EO is the backbone of glacier monitoring in a continuous and large-scale manner (Figure 3). Altogether, 109 EO-based glacier publications were reviewed, whose study area(s) are located in the cold regions of Europe. The dates of these publications span a period from 1993 to 2017. In Figure 4, basic information on study objective(s), study area(s) and exploited sensor(s) is summarized. A study may cover multiple objectives, study areas and/or EO data types. Therefore, the presented statistics are inherently “multi-counted”. As shown in Figure 4, glacier motion has aroused the greatest interest among glaciologists and geoscientists, being covered in almost one-third of the reviewed glacier related articles. It is followed by interest in glacier mass balance (24%) and glacier delineation (22%) studies. However, fewer publications (11%) concerning glacier volumetric change studies were found. The spatial distribution of the study areas show a strong bias towards Norway including Svalbard, followed by the European Alps. Glacier inventories are normally built for 5–10 years intervals and therefore yield the smallest number of reviewed articles. More than half of the studies focus on Norway (including Svalbard). Some 30% of the studies were carried out in the European Alps, and therein the French Alps and the Swiss Alps were much more frequently studied. Notably, Iceland also received great attention, mainly due to the existence there of Europe’s second largest ice cap (in area), Vatnajökull, after the Austfonna ice cap on Svalbard. From a sensor perspective, optical EO data were used in over 60% of the reviewed glacier-related articles. Landsat fleet data are the most favoured for analyses. Synthetic Aperture Radar (SAR) data, especially from the European Remote Sensing-1/2 (ERS) archive, are also frequently used. The value of ERS-1/2 data is significantly pronounced in SAR interferometry thanks to the three-day repeat cycle during the “Ice Phase” and 1-day temporal baseline of ERS-1/2 tandem operations. Despite the small quantity of solely Light Detection And Ranging (LiDAR)-based studies, the value of LiDAR for Digital Elevation Model (DEM) calibration and verification should not be ignored.

3.1.1. Glacier Delineation

Glacier extent is an obvious parameter describing glacier dynamics. According to the guidelines provided by GLIMS [70], one should circumscribe the glacier bodies of ice and snow as well as all debris-covered parts at the end of the melt season to produce glacier outlines. To date, various techniques exist in the field of glacier delineation using remote sensing data, with optical imagery playing a major role. The most straightforward method of glacier delineation is manual digitization. For better visual interpretation, the background image is often displayed as a false color composite to enhance image contrast, such as with ASTER bands 4, 3, 2 [71], TM5, 4, 3 or TM4, 3, 2 [72,73,74], ETM+ 742 [75]. Several (semi-)automated techniques have been developed to extract glacier bodies of snow and ice from EO data. Among them, taking advantage of the significant spectral contrast between the Short Wave Infra-Red (SWIR) and the visible band, the band ratio method (thresholding ratio image, e.g., [76,77,78,79]) is the most widely acknowledged and utilized algorithm. Taking Landsat TM band configuration as an example, commonly selected band ratios (hereafter in this section we exemplify them by Landsat TM bands) are: (1) TM3/TM5 [80,81,82,83], (2) TM4/TM5 [76,77,78], (3) Normalized Difference Snow Index (NDSI) [84]:

$$\mathrm{NDSI}=\frac{\mathrm{TM}2-\mathrm{TM}5}{\mathrm{TM}2+\mathrm{TM}5}$$

(2)

In the case of TM3/TM5, an additional threshold for the blue band (i.e., TM1) is required to further improve the mapping result in shadowcast areas [85]. Lastly, manual correction and median filtering are recommended to improve the results [79]. Further information with regard to the threshold sensitivity, different band combinations and spatial resolutions can be found in Paul et al. [86]. By far, this algorithm has become the backbone for delineating glacier outlines in many glacier inventories (e.g., [9,80,81,82,83,87,88]). Other classification algorithms have also been applied, but scarcely for glacier delineation: Serandrei-Barbero et al. [89] applied fuzzy contextual classification to map the glaciers in the Italian Alps; Robson et al. [90] implemented Object-Based Image Analysis (OBIA) to classify more than 100 glaciers in the Austrian Alps for three time steps semi-automatically.

For the sake of evaluating glacier outlines obtained from the above-mentioned methods, a precursory comprehensive comparison among manual digitization, Digital Number (DN) based band-ratio thresholding (e.g., TM3/TM5 and TM4/TM5), a supervised classification (Maximum Likelihood) and a unsupervised classification (ISODATA) was carried out by Paul [91]. This study was furthered by later researches [74,79,85,86,92] in which atmospheric effects, threshold sensitivity and spatial resolution as well as a comparison between simple ratio and NDSI were investigated. The key conclusions from these studies are:

• Generally, DN based simple ratio methods had the best results (fast, robust and strict to the pixel boundary). Digitization and the band-ratio method are consistent in clean glaciers;
• Simple ratios are relatively insensitive to threshold selection. Yet the additional TM1 threshold is of higher threshold sensitivity;
• Atmospheric correction for TM2 is important when NDSI is used;
• Manual correction is recommended to correct miss-classifications caused by debris cover, proglacial lakes, shadowed areas, and clouds.
• Using thresholds of the panchromatic band (e.g., OLI8) divided by SWIR could improve the spatial resolution, but this increases the workload in manual correction.

Since mapping shadowcast glacierized areas is the most critical point by far [85], TM3/TM5 with TM1 correction or NDSI is therefore often preferred. However, the local environmental settings concerning vegetation and/or water conditions must always be considered to select the most proper ratio [74].

Presently, automated identification of debris cover remains a well-recognized challenge particularly for optical data [81,92,93], mainly due to its spectral similarity to the surrounding terrain. Even though debris cover can be properly digitized manually in most cases, it would be extremely cumbersome on a large scale. To tackle the problem, some EO-based techniques have been developed mainly based on the thermal band(s). A pioneer study was undertaken by Lougeay [94], who explored the potential of thermal remote sensing to identify different moraine types and glacier ice based on thermal contrast. The research was verified by Taschner and Ranzi [95], Ranzi et al. [96] and Mihalcea et al. [97] in the Italian Alps using Landsat TM/ETM+ and ASTER images. These studies revealed that the strongest thermal contrast appears around midday. Also, several empirical methods using the relationship between debris thickness and temperature have been developed (e.g., [96,97]). An obvious drawback of these methods is the lack of spatiotemporal transferability, since field measurement and recalibration are imperative. As alternatives, Paul et al. [98] developed a semi-automatic multisource (i.e., multispectral data and DEM) method to map the supraglacial debris, and compared it to multispectral data based Artificial Neural Network (ANN) classification. They found this multisource approach to be superior to ANN classification solely based on multispectral data. Furthermore, Foster et al. [99] exploited a physically based model to predict the thickness of supra-glacier debris cover from ASTER thermal band data. However, results in regions where debris is thicker than 0.5 m remain error-prone. Unlike optical imagery, SAR images are rarely used for glacier delineation. Since spaceborne SAR sensors are usually of single frequency and radar penetrates dry snow and ice, this limits their suitability for correctly delineating glaciers [100,101]. Nevertheless, an Interferometric Synthetic Aperture Radar (InSAR) coherence image is still useful for mapping glacier extent using low coherence (due to glacier motion and glacier surface condition changes) as an indication of glacier existence [90,102]. The concept has been realized by a four-step decision tree introduced by Atwood et al. [102], including coherence indicators, surface slope indicators, morphological operations and patch size analysis.

3.1.2. Glacier Motion

Motion is one of the most important characteristics of glaciers. Glacier retreat/advance can be easily derived by differentiating bi-(/multi-)temporal glacier outlines. For instance, Holobâcă [103] identified the retreat of the Mount Elbrus glacier system using bi-temporal change detection. The glacier outlines were derived with a thresholded ratio image from two Landsat-TM images (1985 and 2007) acquired at the end of the ablation season. The results revealed an overall glacier retreat rate of 0.4% per year. Such change detection provides an intuitive areal change, which is essential for evaluating the intensity of the glacier-wide advance/recession. To extract the velocity field and the displacement within the glaciated area(s), more advanced techniques are needed. So far, there are two leading methods for velocity field extraction, image offset-tracking and D-InSAR (Differential InSAR).

Image offset-tracking can be applied to both optical and SAR data. A 2-D velocity field can be derived by determining the movement of identifiable supra-glacier objects (e.g., crevasses, moraines) or speckles, using a window (also known as a search chip/patch) moving with an n-pixel step to find the peak cross-correlation. In the optical field, this is usually termed “image matching”. For instance, Kääb et al. [104] presented the work of extracting annual and summer surface velocity fields based on 15-m resolution ASTER 3N band and Landsat ETM+ pan band imagery during four time periods between 1999 and 2002 in Kronebreen, Svalbard. The authors pointed out the high complementarity of “image matching” to existing D-InSAR products. Berthier et al. [105] measured the surface displacement of glaciers in the Mont Blanc area by matching successive SPOT5 image pairs. The results were validated against the Differential Global Positioning System (DGPS) data. It showed an uncertainty on the order of one-fifth of the pixel size over a sensor revisit time, which is comparable to the ERS D-InSAR result. Heid and Kääb [106] noted that the attitude-variable subpixel noise inherent in whisk-broom systems (e.g., Landsat TM/ETM+) reduces the accuracy compared to ASTER images. Recently, Kääb et al. [107] applied the algorithm to Sentinel-2A data with a particular focus on radiometric and geometric performance. The authors identified bands 4, 8 and 11 of Sentinel-2 as seemingly the most suitable for glacier horizontal motion derivation based on the image matching method. Meanwhile, 12-bit radiometric resolution improves the performance in shadowcast areas and avoids saturation in snow/ice in comparison to conventional 8-bit sensors. A detailed comparison of different image matching techniques can be found in [74,108]. Also, to measure glacier motion in Europe, speckle/feature and coherence tracking are widely applied to TerraSAR-X (X-band) data for the Swiss Alps [109] and Svalbard [110], ERS-1/2 SAR (C-band) in Svalbard [111], Radarsat-2 (C-band) in Svalbard [112], and JERS-1 (L-band) in Svalbard [113] using amplitude/phase patterns and coherence. When intensity is used, it is commonly called a cross-correlation optimization procedure; coherence tracking refers to the fringe visibility algorithm or coherence optimization procedure. The patch size of intensity depends on the coherence. If coherence is preserved, the speckle patterns are correlated. Small patch size is therefore efficient for achieving high accuracy [111]. Otherwise, a very large image-path size is required, which is a disadvantage for accuracy and computational efficiency [111].

D-InSAR is an alternative to offset-tracking. It includes: a two-scene technique (DEM elimination technique), a three-scene technique (three-pass technique) and a four-scene technique (four-pass technique) [114]. D-InSAR is amongst the most accurate techniques [109,111]. Whereas, it is restricted to a short temporal baseline in order to avoid the loss of coherence. A significant drawback of D-InSAR is that it only depicts motion in the slant-range direction. Only in a special case, when dual-azimuth image pairs (ascending and descending passes) are available and coherence is retained, might D-InSAR conquer this problem. Otherwise, an azimuth displacement component derived from coherence tracking is often combined with the D-InSAR result to build the motion both in slant-range and in the azimuth direction. Strozzi et al. [115] rendered the surface motion of Unteraargletscher (Swiss Alps) in the line-of-sight (LOS) direction by applying D-InSAR to two pairs of ERS-1/2 imagery with a one-day temporal baseline. They elucidated the influenced of the LOS: displacement perpendicular to the LOS is not detectable by D-InSAR. Murray et al. [116] employed dual-azimuth D-InSAR and intensity offset-tracking of ERS imagery to explore the surge mechanism in northern Svalbard. When a 3-D velocity field is desired, the flow direction is needed. This can be obtained by assuming that the glacier ice bodies flow parallel to the surface [117]. It has to be mentioned that the flow direction assumption is error-prone. In fact, flow tends to be upward in the ablation zone and downward in the accumulation zone [118]. Alternatively, especially when a DEM is unavailable, the surface is assumed to be horizontal. In this case, a 2-D velocity field is enough.

3.1.3. Glacier Elevation and Volumetric Change

In terms of measuring the response of glaciers to climate change, glacier volumetric change is a more suitable indicator than glacier extent variation since the glacier extent appears to shrink with time delay compared to the glacier’s meltdown in height [119]. The volumetric change can be calculated by multiplying the average glaciated area by the glacier elevation change within a certain time span [120]. To calculate glacier elevation changes, spaceborne EO provides two main types of elevation products, DEM (interferometric DEM and stereoscopic DEM) and height information from satellite altimetry measurements. These (quasi-)global elevation products are often treated as baseline topography for mapping glacial elevation change. There are three global elevation data sets in common use:

(1)

The Shuttle Radar Topography Mission (SRTM) C-band DEM covers 60°N to 57°S with spatial resolution levels of 1 and 3 arc sec; the vertical accuracy of SRTM DEM is of ±16 m at 90% confidence [121].

(2)

The ASTER Global Digital Elevation Model (GDEM, updated to Version 2 on 17 October 2011), is the most representative stereoscopic DEM. It has spatial coverage from 83°N to 83°S (~99% of Earth’s landmass) with 30 m resolution in 1° × 1° tiles. The vertical accuracy of ASTER GDEM is specified as 17 m at the 95% confidence level [122].

(3)

The Ice, Cloud, and land Elevation Satellite—Geoscience Laser Altimeter System (ICESat-GLAS), which is not only an elevation data input for DEM differencing, but also important as auxiliary data for calibrating and validating interferometric DEM and stereoscopic DEM.

Digital Elevation Model (DEM) subtraction is the most straightforward way to calculate glacier elevation change. The DEM(s) can be derived from photogrammetry using aerial photography, airborne SAR interferometry, airborne LiDAR, and topographic maps. Regarding the EO data, optical sensors of stereo geometry (e.g., ASTER, SPOT5, IKONOS, Quickbird, WorldView, Pléiades) have enduring popularity. Of these optical sensors, ASTER is of particular interest thanks to its free accessibility (unlimited public access since 1 April 2016), extraction software based on automated algorithms (e.g., SilcAst, PCI Geomatics), comparably short revisit time interval, long-term operation (16-day revisit time from 2000 to present), and relatively high resolution at near-global coverage (15 m nadir and backward for band 3). Before the subtraction, reprojection, co-registration and resampling are imperative for heterogeneous DEMs. Paul et al. [74] suggested downscaling the finer resolution image using algorithms more complicated than nearest neighbour to avoid artefacts and horizontal disalignment within a sub-pixel. Moreover, when the pixel size difference is at least two times, a block average filter is recommended. Regarding co-registration, Nuth and Kääb [123] developed a three-step workflow for DEM co-registration as well as bias correction. The algorithm was tested in Svalbard and New Zealand using STRM, ICESat, ASTER GDEM and automatically generated ASTER and SPOT5-HRS DEMs. Meanwhile, the authors found significant elevation-related biases in some stereoscopic DEMs due to the satellite acquisition geometry. Berthier et al. [124] evaluated the performance of sub-meter stereo imagery from Pléiades for glacier DEM generation, and found that the 12-bit radiometric resolution makes Pléiades superior to ASTER and SPOT (both are 8-bit encoding). This is because higher radiometric resolution improves the contrast over snow/ice and reduces the signal saturation risks.

Similarly, users can also use InSAR to generate interferometric DEMs. However, radar penetration of the snow and ice makes interferometric DEM subtraction more complicated for glacial studies [123,125,126]. Nevertheless, SAR still provides valuable elevation change information, especially when the glacier surface experienced ice-water phase change (e.g., precipitation and melting). Magnússon et al. [127] explored the potential of InSAR for detecting regional glacier uplift and subsidence. They identified the localized anomalies in the vertical ice motion of Vatnajökull using four InSAR scenes (each one is over 24 h) derived from ERS-1/-2 data. The anomalies were later explained by the subglacial hydrology. Even though ICESat-GLAS has the best vertical accuracy (0.15 m [128]) and the best consistency [123] of these elevation products, it has two major drawbacks for glacier elevation change measurement: approximately 70 m diameter footprint with 170 m ground spacing along track, and shift of the footprint centers [129]. To tackle the problem, two leading methods exist for glacier elevation change estimation based on altimetry data: the repeat-track method and the cross-over method [130,131,132]. To date, ICESat data are mainly utilized for the Antarctica and Greenland ice sheets. In Europe, Moholdt et al. [131] tested two ICESat repeat-track methods for deriving glacier elevation changes in Svalbard and validated the results against cross-over points and the SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies (SPIRIT) DEM. The results confirmed the validity of repeat-track methods for mapping short-term Arctic glacier elevation changes. The authors also suggested the application of a seasonal data filter before the plane fitting for the purpose of mean annual elevation change rate calculation. For the first time, elevation and volumetric change in the European Alps (the Grosser Aletschgletscher in the Swiss Alps) were derived from ICESat in combination with external DEMs by Kropáček et al. [133]. Due to rugged terrain and limited ICESat footprints in the mid-latitude mountain glaciers, the repeat-track method and the cross-over method were found to be inappropriate. Moreover, ICESat observations are proven to be inaccurate for areas with >10° slope [134]. Those observations (>10° slope) are hence replaced by the mean value of the two adjacent records. To obtain the elevation changes, the authors fitted a linear function to ICESat-DEM differences along ICESat tracks. The authors emphasized the use of high-quality DEM for accurate results retrieval. Nuth et al. [135] and Kääb [129] compared the performance among ICESat, DEMs, and topographic maps to derive elevation change in Svalbard. Both studies shed light on the value of ICESat for investigating glacier elevation/volumetric changes at a regional scale.

3.1.4. Glacier Mass Balance

Simple changes in glacier surface area or length are not adequate to represent the variation in glacier mass balance at a decadal scale [136]. The term “glacier mass balance” refers to the variation of the mass of (part of) the glacier over a certain time span [137]. Conventionally, glacier mass balance is directly measured via snowpits and stakes. Europe is one of the areas with the richest and longest glacier mass balance records, including Storglaciären (Sweden), Limmern, and Plattalva glaciers (Switzerland), Sarennes Glacier (France) and Storbreen (Norway) [62,138,139]. However, such glaciological measurement has an obvious drawback, the laboriousness of long-term and large-scale monitoring. EO-based mass balance retrieval has thus gained increasing attention in recent years.

The most widely-applied EO-based method is the “geodetic method” [140]. Hannesdóttir et al. [120] employed this method to reconstruct the mass balance of glaciers in the southeast part of the Vatnajökull ice cap for a seven-year span from ~1890 to 2010 based on multiple sources including glacial geomorphological features, historical photographs, maps, aerial images, DGPS measurements, a LiDAR survey and several optical satellite images (i.e., SPOT5, MODIS, Landsat). The glacier-wide geodetic mass balance, B_G, can be expressed as:

$$B_G=\overline{\rho }\times \frac{∆H}{∆t}$$

(3)

where $\overline{\rho }$ denotes glacier-wide average density, which is either obtained from field measurement or oftentimes assumed to be 900 m·kg⁻³ [141,142]. The areal-mean thickness/elevation change ∆H is averaged through a certain time-span ∆t [141]. The uncertainties of this method have been well-documented [141,143,144]. Apart from bias in ∆H (see Section 3.1.3) when converting elevation changes into water equivalent volume changes, the simplified assumption of ice density $\overline{\rho }$ is error-prone especially in compressible firnpack [142,143,144]. In addition, the time-span ∆t is significantly uncertain with regard to the acquisition dates of the EO data [141]. Nevertheless, Pope et al. [145] compared mass balance of the whole Langjökull ice cap derived from a glaciological survey, a geodetic method (using a DGPS interpolated DEM, a SPOT5 stereoscopic DEM and a LiDAR DEM), and a surface mass balance model during 1997–2007. The three methods led to good and consistent results. Zemp et al. [146] reported a RMSE = 0.23 m water equivalent (m. w. e.) of the annual discrepancy between glaciological measurement and the geodetic method for 12 European glaciers over 50 recording periods.

To construct mass balance time series, the geodetic method is usually hampered by the scarcity of long-term consistent annual elevation data, so alternatives are desired. Equilibrium Line Altitude (ELA) is an appropriate indicator for glacier mass balance [147]. ELA refers to the mean altitude of the equilibrium line which divides the accumulation zone and the ablation zone and therefore has zero climatic mass balance [137]. Based on ELA variability, Rabatel et al. [136] reconstructed glacier-wide mass balance time series for 30 glaciers in the French Alps from 1983 to 2014. The ELA-based mass balance (B_E) can be obtained via a two-step calculation following Equations (4) and (5) [63,148]. Firstly, the theoretical ELA under the steady state (glacier-wide mass balance = 0), ELA_theo, over the period under consideration is expressed as:

$$EL{A}_{theo}=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}EL{A}_t+\frac{\overline{B}}{\frac{\partial b}{\partial z}}$$

(4)

where ELA_t is the ELA at the tth year. $\overline{B}$ represents the average glacier-wide mass balance over the studied n years, which can be determined either by using field observation or a geodetic method (Equation (3)). $\partial \mathrm{b}/\partial \mathrm{z}$ is the mass balance gradient (activity coefficient) across the (Equilibrium Line) EL, and is fixed at 0.78 m. w. e. per 100 m for the Alps [63]. In the second step, ELA-based glacier-wide mass balance (B_E) is computed as:

$$B_E=\left(EL{A}_{theo}-EL{A}_t\right)\frac{\partial b}{\partial z}$$

(5)

In this step, the variation of ELA changes is transferred into mass balance variation by the mass balance gradient. The mass balance gradient varies significantly according to glacier locality (climate and therefore accumulation/ablation characteristics) [149]. However, when the specific mass balance gradient is unknown, an average value near ELA for a given period and climatic characteristics does not weaken the outcome [63]. More crucially, the EO acquisition date must be restricted to the end of the hydrological year, during which the glacier must be high enough to preserve the accumulation zone [63,136]. The method was validated by field measurement data and yielded a Mean Absolute Error (MAE) between 0.24 and 0.51 m. w. e. [63].

As an alternative to ELA, albedo from satellite measurements is useful for depicting glacier mass balance in a time-series manner. Dumont et al. [150] explored the relationship between measured annual mass balance and MODIS-derived glacier albedo. They observed strong correlation (R² > 0.9) between the albedo and mass balance and argued that albedo is more informative than ELA with regard to glacier surface. The “satellite albedo method” was developed and applied for Vatnajökull (Iceland) using NOAA-AVHRR imagery [151] and Svalbard using MODIS data [64]. The satellite-derived mass balance (B_S) can be obtained by:

$$B_s=\frac{1}{L_f}{{\displaystyle \int }}_{t_1}^{t_2}max\left\{I_0{\tau }_{atm}\left(1-\overline{{\alpha }_{sat}}\right)+Q_0,0\right\}{d}_t$$

(6)

where I₀ denotes the daily average incoming radiation at the top of the atmosphere, τ_atm is the atmospheric transmissivity, $\overline{{\alpha }_{sat}}$ represents the mean satellite-derived albedo over the glacier surface, and Q₀ is the sum of the net long-wave radiative flux and the turbulent fluxes. τ_atm and Q₀ can be derived from in-situ measurements [64]. The entire melt energy in the ablation season (t₁ to t₂) is calculated by integrating at a one-day time step between t₁ and t₂. Lastly, the amount of melt is converted by dividing the total melt energy by the latent heat of fusion, 1/L_f. The method’s disadvantage is its incapacity to retrieve absolute mass balance. Interannual variability in mass balance is more informative than is an absolute value with regards to the climate sensitivity of glaciers [64]. One aspect of the glacier mass balance is still missing, which is the spatial distribution. Dumont et al. [152] assimilated MODIS-derived albedo and snowpack models to estimate the spatially distributed mass balance of Glacier de Saint-Sorlin in the French Alps. The results were compared with Automatic Weather Station (AWS) data and showed nonbias (RMSE < 0.5 m. w. e.).

In addition to the aforementioned algorithm for glacier mass balance calculation, other methods also exist, such as the Accumulation Area Ratio (AAR) method [153] and the “flux gate” method [144] (only for tidewater glaciers). However, they were seldom applied in Europe or were restricted to a specific region (e.g., Svalbard) for the purpose of quantitative mass balance evaluation. Hence, for detailed information, we suggest the reader refer to the individual authors.

3.1.5. Glacier Inventory

The ambitious concept of building a global glacier inventory was first proposed in the 1950s; it was later known as the World Glacier Inventory (WGI) under the direction of F. Müller [9]. Cogley [154] extended the WGI to the WGI-XF with records for more than 130,000 glaciers. After half a century of effort it still only has a global coverage of around 48%, and it includes only point data, which omits a vast amount of spatial information. Nowadays, the integration of EO in glaciological studies, the further development of EO-based methodologies, and free access to some EO data sets greatly facilitates the generation of such a global inventory [9,78,79,155]. On 22 February 2012 the first complete global glacier inventory (excluding Greenland and the Antarctic)—Randolph Glacier Inventory (RGI, version 1.0)—was released. It is EO data based (including data from the Landsat fleet, ASTER, SPOT5 HRS, SRTM DEM, ASTER GDEM, etc.) and was developed within a short time (1–2 years) to fulfill the urgent need of the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR5) [9]. At a regional scale, glacier information has been well-documented in Europe for over a century. Meanwhile, a large number of glacier inventories exists based on either conventional or EO data.

Principally, glacier inventories include two kinds of information, glacier extent and glacier-related parameters (e.g., identification code, coordinates, total surface area, length, elevation, mean aspect and slope). For the purpose of glacier delineation in regional glacier inventories, orthophotographs and historic maps are a valuable resource because they contain information gathered before the availability of EO data). The primarily EO-based glacier inventories for Europe (excluding Greenland) are summarized in

Table 1. References for existing glacier inventories derived from satellite imagery for European glaciers (excluding Greenland).

Location	Year	Sensor	DEM	Delineation Technology 1	References
Globe (RGI)	1999 or later	Landsat fleet, etc.	SRTM, SPOT DEM, ASTER GDEM2	T + D	[9,156]
European Alps	2003	Landsat fleet	STRM	T	[81]
Austrian Alps	1985, 1992	Landsat fleet	—	T	[78]
	1985, 2003, 2013	Landsat fleet, ALOS PALSAR	SRTM, Local DEM	S	[90]
French Alps	1967/71, 1985/86, 2003, 2006/09	Landsat fleet	Local DEM	S	[157]
Italian Alps	1987, 2009	Landsat fleet	Local DEM	T	[87]
Swiss Alps	1985	Landsat fleet, SPOT	Local DEM	S/U/T/M	[79,158]
	1985, 1992, 1998/99	Landsat fleet	Local DEM	T	[88]
Iceland	1890–2010	Landsat fleet, MODIS, SPOT	Local DEM	D	[120]
Norway	LIA-2003	Landsat fleet	—	D	[73]
	1999	Landsat fleet	Local DEM	T	[83]
	1966, 1997, 2003, 2006	Landsat fleet	Local DEM	T	[82]
	1947/85, 1988/97, 1999–2006	Landsat fleet	Local DEM	T	[80]
Svalbard	1930s–2010	ASTER, Landsat fleet, SPOT	SPOT DEM, ASTER GDEM2	D	[159]
	2000–2006	ASTER, Landsat fleet,	ASTER generated	T+D	[71]

¹ Delineation Technology for EO data: T = Thresholded Ratio Image; S = Supervised Classification; U = Unsupervised Classification; D = Digitization. LIA: Little Ice Age.

. The Landsat fleet is by far the most often used backbone of glacier inventory generation. The Landsat fleet has several merits compared to ASTER or SPOT, which are also popular sensors in glaciological studies (see Figure 3): large swath width, free accessibility, long-term records, georeferencing and orthorectification [80]. The methods were summarized in Section 3.1.2 above. As seen in Table 1, ratio image thresholding and manual digitization are the most frequently applied delineation methods. Regarding the updating interval, 5 to 10 years is suggested [72] in order to reveal highly variable changes as well as the spatial patterns of changes [83]. Topographical information is the most often extracted glacier parameter, and can be automatically calculated from existing freely available DEMs. The suitability of national and global DEMs with regard to glacier-specific topographic parameter extraction was investigated by Frey and Paul [126] for a sample of 1786 Swiss glaciers. The study concluded that the first choice is local, high quality DEMs. When such DEMs are not available, SRTM DEMs are preferred to ASTER GDEMs. Notwithstanding better results obtained from SRTM DEMs, both DEMs are useful for calculating topographic parameters in glacier inventories. In terms of the calculation, Paul et al. [155] recommended the standardized calculation of relevant glacier parameters (9 basic and 7 auxiliary parameters). Apart from topographical attributes, mass balance information is one of the most essential variables for glacier studies. Mass balance information is seldom included. Recently, some inventories have started to include or update mass balance information. For example, Hannesdóttir et al. [120] included the average geodetic mass balance for the outlet glaciers of southeast Vatnajökull. The latest RGI (Version 5.0, released 20 July 2015) linked some of the RGI glaciers to mass balance data from the WGMS mass balance database [156].

Relying on EO and other aforementioned methods, there are several general findings on the contextual, geoscientific analysis of European glaciers:

• Glaciers in Europe are shrinking, a few have completely disappeared between the “Little Ice Age” and the present [73]. The average glacier area loss is approximately −0.4%·a⁻¹~−2%·a⁻¹ for land-terminating glaciers since the 1970s [80,81,157,160], and the reduction rate is increasing in some of the observed glaciers [157]. Over the past three decades, the reduction rate is comparably small (ca. −0.23%·a⁻¹) in the whole Svalbard Archipelago [159] and below −0.2% in the Svaritisen and the Jostedalsbreen regions of Norway [82,83].
• In Svalbard, the velocity of surging glaciers ranges from 1.9 m·d⁻¹ to 5 m·d⁻¹ during the past 20 years [112,116,161,162,163], for example, at Fridtjovbreen, Kronebreen, Kongsbreen, Monacobreen, Nathorstbreen, and Perseibreen. For land-terminating glaciers, Berthier et al. [105] identified the highest velocity as 500 m·a⁻¹ on the icefalls of glaciers in the Mont Blanc region in the year 2003.
• Generally, the glaciers in Europe are thinning at the rate of −0.1 m·a⁻¹ to −6 m·a⁻¹ at the beginning of the 21st century [124,129,131,133,135,164], which is especially accentuated in the European Alps [62].
• European glaciers are experiencing mass loss [71,120,136]. In terms of the geodetic mass balance, the highest mass loss (approximately −1.5 m. w. e.·a⁻¹) was observed by Hannesdóttir et al. [120] in the outlet glaciers of southeast Vatnajökull (Iceland) from 2002 to 2010. The average glacier-wide mass loss in the French Alps is between −0.39 m. w. e.·a⁻¹ and −0.93 m. w. e.·a⁻¹ between 1979 and 2011 [136], and around −0.02 m. w. e.·a⁻¹ to −0.41 m. w. e.·a⁻¹ in Kongsvegen and Kronebreen (Svalbard, Norway) for the periods 1966–1990/95 and 1990/95–2007 [144].

For detailed information, the reader should refer to the individual authors or the review paper on European mountain cryosphere dynamics written by Beniston et al. [69].

3.2. Snow

Quantitative analyses of snow are vital for a variety of society-relevant research and planning applications, including water resources management as well as understanding and predicting likely future climate changes [69,165]. For the cold regions of Europe, a total of 102 publications were found that used EO data to observe snow cover. These reviewed articles were published in the time frame from 1980 to 2016. The foci of the investigated papers can be categorized into five topics: snow cover extent (34%), snow cover fraction (17%), snow characterization (25%), snowmelt (15%) and EO snow products (9%) (Figure 5). In addition to the methodic study objectives, more information concerning study areas and applied sensors are illustrated in Figure 6. More than 40% of the reviewed publications focused on the European Alps. Amongst the European Alps, the Italian Alps gained the most attention (17%), almost twice as much as the second most investigated Alpine area (the Swiss Alps, 9%). Apart from the European Alps, Finland and Norway are the two other hot spots which are considered in some 15% and 20% of the found articles, respectively. Optical sensors (Landsat and MODIS) are the most commonly applied. Active and passive microwave (PM) sensors were employed in 25% of the studies.

3.2.1. Snow Cover Area

Snow Cover Area (SCA) is a basic parameter of snow observation. In the optical EO domain, band-ratio thresholding is the leading algorithm for mapping SCA. It utilizes the significant signal contrast between the visible band and the SWIR band to distinguish snow cover from ambient surroundings. Among the existing threshold-based algorithms, SNOMAP [166,167] is the most quintessential example. The originally proposed NDSI ≥ 0.4 [166,167] has been universally applied to identify snow covered area. Alternatively, the Normalized Difference Snow and Ice Index (NDSII) that uses the red band instead of the green band in the calculation can be employed for sensors without a green band, such as SPOT-VEGETATION. The results derived from NDSII have proven to be very comparable to the SCA derived using NDSI [168]. Based on the NDSII algorithm, Xiao et al. [169] mapped SCA using SPOT-VEGETATION data from 1998 to 2001 in the pan-Arctic (north of 45°N). Whereas, using solely NDSI/NDSII may erroneously include very dark targets (e.g., black spruce forest, water bodies). To avoid such misclassification, further threshold(s) were introduced in addition to NDSI ≥ 0.4. For example, Winther and Hall [170] included $\rho$b4 (reflectance at band 4) ≥ 0.11 for Landsat TM imagery, Xiao et al. [169] used $\rho$b3 ≥ 0.11 for SPOT-VEGETATION images, Foster et al. [171] additionally restricted $\rho$ b2 ≥ 0.1 and $\rho$ b4 ≥ 0.1 to MODIS data, etc. A well-recognized problem, identifying snow within forested regions, should also be noted. Usually, this problem is addressed by supplementarily computing the Normalized Difference Vegetation Index (NDVI). Afterwards, a NDSI—NDVI threshold field can be employed to identify snow in forest pixels according to Klein et al. [172]. Moreover, to further eliminate contaminated snow pixels (e.g., cloud, coastal sand, aerosol influenced pixels), a thermal mask (283 K) was suggested by Romanov et al. [173].

In addition to the NDSI algorithms, several other methods were also developed and applied in Europe. For instance, (1) maximum likelihood classification was employed by Turpin et al. [174] in Arctic Sweden and the Swiss Alps using Landsat TM imagery; (2) decision trees were used by Malcher et al. [175] to classify snow, forest snow, fractional snow, and clouds in the Austrian Alps using MODIS and MERIS data; (3) Hüsler et al. [176] developed a test-score aggregation technique for snow detection using historical NOAA-AVHRR data over the European Alps; (4) for geostationary satellite data like Meteosat Second Generation-Spinning Enhanced Visible and Infrared Imager (MSG-SEVIRI), Siljamo and Hyvärinen [177] introduced a multi-rule thresholding based snow cover algorithm and demonstrated its application in Europe; (5) Pepe et al. [178] used crisp and fuzzy soft classifiers on MODIS and ASTER images to extract SCA in the Austrian Alps and the central Italian Alps. (6) Metsämäki et al. [179] developed a linear interpolation based method for SCA estimation in Finland using NOAA-AVHRR data. (7) Notarnicola et al. [180] proposed a new method (using a so-called EURAC algorithm) for mapping near-real-time snow extent over Central Europe (42°~51°N, 5°~30°E) at 250 m resolution, based on the 250 m resolution MODIS bands and NDVI values.

Due to the inherent challenges of SAR sensors—dry snow penetration and single-frequency configuration—mapping SCA using SAR data independently is difficult. Yet PM sensors with multi-frequency channels hold potential for mapping snow cover utilizing the inverse correlation between surface emission and radar frequency. Thereupon, to map the snow cover of the Northern Hemisphere, Grody [181] and Grody and Basist [182] developed a decision tree based snow identification algorithm using the Defense Meteorological Satellite Program-Special Sensor Microwave Imager (DMSP-SSM/I) PM data. Kongoli et al. [183] adapted the algorithm [181,182] to the Advanced Microwave Sounding Unit (AMSU)-A and AMSU-B and compared the result with the IMS snow cover product. They found that the IMS product was superior in identifying early season snow cover and separating snow-free land. The drawbacks of PM with regard to SCA include restriction to dry snow, difficulty in determining thin snow (<5 cm) [171], very coarse resolution (e.g., SSM/I with 25 km resolution), and attenuation of forest [184]. Contemporaneously, an optical-PM-scatterometer blended algorithm (i.e., AFWA NASA Snow Algorithm, ANSA) was developed by the US Air Force Weather Agency (AFWA) and NASA jointly [171,185,186]. Therein, three EO data/products are used, namely MODIS, the Advanced Microwave Scanning Radiometer (AMSR)-E, and the Quick SCATterometer (QuikSCAT). Moreover, the authors demonstrated the possibility of mapping dry, incipient melting and actively melting snow using AMSR-E (37 GHz), AMSR-E (19 GHz) and QuikSCAT (13.4 GHz), respectively [171]. The algorithm has been validated by Casey et al. [187] using in-situ data in Finland, which resulted in good consistency.

Monitoring the temporal SCA signatures from a time series is also crucial for snow cover status evaluation. For this purpose, metrics like Snow Cover Duration (SCD) Snow Cover Start (SCS) and Snow Cover Melt (SCM) [188] are the most utilized. Based on these metrics, Dietz et al. [32] assessed the temporal snow cover characteristics in Europe between 2000 and 2011 using MODIS-Terra/Aqua daily snow cover products (MOD10A1 & MYD10A1). The authors also pointed out that the usability of optical data is strongly limited by the relatively frequent cloud obstruction and polar darkness in Europe during winter. To fill these gaps, the following algorithms were developed especially for MODIS images and used often in combination: (1) a Terra-Aqua combination, (2) a conservative backwards/forwards temporal filter, (3) a regional snowline-based correction and (4) a seasonal filter [32,188,189,190,191]. These techniques were integrated and improved, for example, by Da Ronco and De Michele [192], to bridge cloud cover gaps in MODIS daily snow cover products for the Po River Basin.

3.2.2. Snow Cover Fraction

“Mixed pixels” is a well-known inherent problem of remotely sensed imagery; especially in medium and coarse imagery. In this regard, fractional snow cover (FSC) may provide superior information to a binary classifier. Until now, fractional snow cover can be retrieved using (semi-)empirical methods, Spectral Mixture Analysis (SMA) and reflectance modelling.

The Norwegian Linear Reflectance-to-Snow-Cover Algorithm (NLR) and the MODIS-based NDSI-FRA (FRAction of snow-covered area) algorithm are two prevalent (semi-)empirical methods to retrieve FSC, which were introduced by Andersen [193] and Salomonson and Appel [194], respectively. The NLR method was established assuming a linear relationship between the pixel signal and FSC [193]. Initially, the algorithm was developed for estimating FSC in Norway based on NOAA-AVHRR data. In fact, it can be regarded as a 2-EM (End-Member) SMA model [195]. Later, Zhu and Woodcock [196] modified the NLR algorithm to build the snow detection module of Tmask (multiTemporal mask, an automated cloud, cloud shadow, and snow detection method) for Landsat data. Using NDSI values, the NDSI-FRA algorithm has been applied to generate MOD10A1 (Terra) and MYD10A1 (Aqua) FSC products [197]. It was originally developed to estimate FSC using MODIS-Terra imagery [194]. Subsequently, Salomonson and Appel [198] expanded the algorithm to data available from the Aqua platform as well. Due to the non-function of MODIS Aqua band 6, the intercept and slope of the linear regression line were determined by using the number of SNOMAP-classified Landsat ETM+ snow pixels within a MODIS 500 m grid as a response variable; these are −0.01 (−0.64) and 1.45 (1.91) for Terra (Aqua).

Spectral unmixing has been employed for sub-pixel FSC retrieval since the 1990s, e.g., [199,200,201,202]. Foppa et al. [203] rendered a linear spectral mixture analysis (LSMA) to NOAA-AVHRR data for the whole European Alps. The authors suggested that the presented processing chain is appropriate for operational and Near Real Time (NRT) applications since it is simple, fast, objective, and reproducible. Later, Foppa et al. [204] validated the results quantitatively over the European Alps based on higher resolution satellite imagery (ASTER). Although their results showed good agreement (R² = 0.78 and MAE = 10.44%), a significant bias towards low FSC was observed in LSMA-based results. Also in the European Alps, Veganzones et al. [205] applied a much advanced unmixing technique to MODIS data. The algorithm includes two parts: a (partially) Constrained Least Squares Unmixing (CLSU) algorithm and an Endmember Induction Algorithm (EIS). Meanwhile, SMA is also a leading technique for solving the “snow in forest” issue. Vikhamar and Solberg [206] developed a constraint LSMA method, SnowFrac, to estimate FSC in the forested areas of Jotunheimen (Norway) using Landsat TM data. The SnowFrac method has three submodels (BirchMod, ShadMod and DiffusMod) integrated in a major unmixing model named SnowFor, simulating the pixel-wise spectral variation due to the presence of snow, trees and snow-free ground. Afterwards, Vikhamar and Solberg [207] applied the SnowFrac method to MODIS-Terra data and validated it using higher resolution images (Landsat ETM+ with 30 m spatial resolution). The validation data set was calculated by firstly classifying 30 m Landsat ETM+ images into binary snow cover maps. Thereafter, the authors aggregated the 30 m classification results into 500 m resolution images at the optimized co-registration position. The accuracy exceeded 80% with R² values of 0.85 (year 2000) and 0.95 (year 2001). Such a validation technique has the drawback of underestimating at low FSC and overestimating at higher FSC [208].

In terms of optical reflectance model-based methods, SCAmod is the most widely used. It was developed by the Finnish Environment Institute (SYKE) initially for estimating SCA in Finland from NOAA-AVHRR data [209,210,211]. The resulting SCA is actually the mean FSC [210], which is a function of observed reflectance ($\rho$_λ,obs(SCA)), snow ($\rho$_λ,snow), ground ($\rho$_λ,ground) and forest ($\rho$_λ,forest) reflectance at wavelength λ under the effective transmissivity τ_λ (Equation (7)).

$$SC{A}_{mod}=\frac{\frac{1}{{\tau }_{\lambda }^2}·{\rho }_{\lambda , obs}\left(SCA\right)+\left(1-\frac{1}{{\tau }_{\lambda }^2}\right)·{\rho }_{\lambda , forest}-{\rho }_{\lambda , ground}}{{\rho }_{\lambda , snow}-{\rho }_{\lambda , ground}}$$

(7)

The model is effective and transferable to any optical EO data. Moreover, it does not need any auxiliary land-cover information apart from a water mask, which can be derived from EO data simultaneously [210]. Salminen et al. [212] found that the SCAmod method performed better when using MODIS and MERIS data than did NOAA-AVHRR images at the end of the snowmelt season. Salminen et al. [212] also emphasized the role of atmospheric correction in retrieving more precise snow information. The accuracy (MAE) of the SCAmod method ranges from 5% to 12% [212]. For a larger scale, Metsämäki et al. [208] presented the application of SCAmod over northern and eastern Europe using MODIS-Terra and ENVISAT-AATSR data. Together with the MOD10_L2 fraction snow product, the authors further validated the SCAmod results against in-situ data. They have concluded that SCAmod (RMSE = 0.11) performs significantly better than does the MODIS FSC product (RMSE = 0.20).

3.2.3. Snow Characterization

To characterize a snowpack, Snow Depth (SD), Snow Grain Size (SGS), Snow Water Equivalent (SWE), and Snow Impurity (SI) are the most often-used parameters. Of these parameters, optical data are useful for retrieving SGS and SI, since snow reflectance depends on the impurity concentration in the visible wavelength range and on SGS in the NIR range [213,214]. Fily et al. [215] discovered an empirical linear relationship between Landsat TM reflectance and SGS. The study was carried out in the French Alps using either a single band reflectance or band-ratios as exploratory variables. The authors suggested the use of a band at 1.2 µm because of the highly pronounced SGS effect and low influence due to atmosphere and SI at this wavelength. Painter et al. [201] developed a Multiple Endmember SMA (MESMA) based model, MEMSCAG (Multiple End Member Snow Covered Area and Grain size), to model FSC and SGS using hyperspectral data. Later, Painter et al. [216] adapted this method to multispectral data with a new model called MODSCAG (MODIS Snow Covered Area and Grain size). The validation against field measurements revealed a slight overestimation (MAE = 51 µm) of SGS from MODSCAG [216]. In order to map light-absorbing snow impurities (e.g., black carbon, mineral dust, and volcanic ash), Di Mauro et al. [217] introduced a Snow Darkening Index (SDI):

$$\mathrm{SDI}=\frac{\mathrm{RED}-\mathrm{GREEN}}{\mathrm{RED}+\mathrm{GREEN}}$$

(8)

The index has been successfully applied to Landsat OLI imagery to monitor mineral dust on snow in the European Alps [217].

So far, PM remains the backbone in SWE retrieval. Hallikainen and Jolma [218] compared the correlation coefficients between Brightness Temperature (BT) functions and SWE. Generally, BT-difference functions perform better than single BT or BT-ratio functions. Later, Hallikainen and Jolma [219] performed a correlation analysis for 17 BT metrics and SWE map-derived values in Finland. Their conclusion suggested that BT_10V − BT_37V and BT_18V − BT_37V are the most suitable (i.e., have the highest correlation coefficients) for estimating SWE in small areas. Given the advantage of these high correlations, there are several empirical methods for SD and SWE calculation that are summarized in

Table 2. Empirical equations for Snow Water Equivalent (SWE) and Snow Depth (SD) calculation.

Equation	Annotation	Reference
$SD=-\frac{GT-0.085}{0.036}$	$GT=\frac{B{T}_{37H}-B{T}_{18H}}{37-18}$	[220]
$SWE=0.27·SD$		[220]
$∆T=F\left(SWE\right)-F\left(SWE=0\right)$	$F=B{T}_{18H}-B{T}_{37H}$	[221]
$SD=1.59·\left(B{T}_{18H}-B{T}_{37H}\right)$		[61]
$SWE=a_0+a_1·\left[\left(B{T}_{19V}-B{T}_{37V}\right)+\left(B{T}_{19V}-B{T}_{19H}\right)\right]$	α0 and α1 are empirical coefficients	[222]
$SWE=c·\left(B{T}_{18H}-B{T}_{37H}\right)$	c is the slope of linear fit 4.8 mm·K−1	[223]
$SD=1.59·\left(B{T}_{19H}-6\right)-\left(B{T}_{37H}-1\right)$		[223]
$SW{E}_{non-forested}=a_0-a_1·\left(B{T}_{36.5V}-B{T}_{18.7V}\right)$	α0 and α1 are empirical coefficients	[224]
$SW{E}_{forested}=a_0-a_1·\left[\left(B{T}_{365V}-B{T}_{18.7V}\right)-B{T}_{18.7V}-B{T}_{18.7H}\right]$		[224]
$SWE=a_0+a_1·\left(B{T}_{23V}-B{T}_{89V}\right)$	α0, α1, α2, α3 are empirical coefficients	[183]
$SWE=a_2+a_3·\left(B{T}_{23V}-B{T}_{31V}\right)$		[183]

.

Artificial Neural Network (ANN) and modelling are alternative techniques for SWE and SD retrieval using PM data. Santi et al. [225] presented an ANN-based algorithm (HydroAlgo) for retrieving SD from AMSR-E data. Santi et al. [226] applied this algorithm to the eastern part of the Italian Alps between 2002 and 2011. The validation indicated better consistency between the HydroAlgo-derived and the measured SD using AMSR-E descending scenes (RMSE = 14.184 cm) than ascending-orbit images (RMSE = 24.652 cm). In terms of the modelling, Pulliainen et al. [227] introduced the Helsinki University of Technology (HUT) snow microwave emission model for SWE retrieval. Essentially, this is a model-based auto-inversion method. Pulliainen and Hallikainen [228] tested the model in the Kemijoki River Basin (Finland) using four-winter SSM/I data sets. Also, the study pointed out that the results were weather-condition (air temperature and precipitation) dependent; therefore, additional information (e.g., the near-surface air temperature) could improve SWE retrieval performance. Tedesco et al. [229] compared the performance of an ANN-based technique, an empirical algorithm (Spectral Polarization Difference method, SPD [222]) and the iterative inversion of the HUT snow emission model in retrieving SWE and SD. The authors found that the best results were obtained from the ANN method trained with experimental data. Satisfactory results were also received using the other two methods when SWE is less than 17 cm or excluding springtime measurements. For mesoscale studies, assimilation of spaceborne PM data and ground-based measurements is indispensable for improving SWE estimation. Pulliainen [230] presented a technique for assimilating PM observations and field measurements to map SWE and SD in boreal and sub-arctic zones. According to the validation, the assimilation technique yielded more accurate SWE and SD prediction than simple field observation interpolation. Later, based on Pulliainen’s [230] assimilation method, Takala et al. [231] produced a 30-year daily SWE product using SMMR, SSM/I and AMSR-E data acquired between 1979 and 2009 for the Northern Hemisphere.

In terms of active microwave data, Hallikainen et al. [232] explored the feasibility of integrating SSM/I and QuickScat data to calculate SWE over 21 test areas in Finland. The results emphasized that a PM-scatterometer-combined SWE retrieval method could improve results derived solely from a PM data set. Recently, in Phase-A of the Cold Regions Hydrology high-resolution Observatory (CoReH₂O) mission supported by ESA, several SWE retrieval algorithms utilizing dual-frequency SAR measurement have been developed. For example, Montomoli et al. [233] assessed the performance of the SWE retrieval algorithm developed for CoReH₂O mission in boreal forests of northern Finland. They found that the satisfactory SWE results can only be obtained in sparse forested (<30%) areas. Cui et al. [234] developed a bi-continuous-VRT model to estimate SWE from X- and Ku-band SAR data based on absorption loss. Meanwhile, the authors discovered good relationships between the single scattering albedo and snow optical thickness at X- and Ku-band [234].

3.2.4. Snowmelt

Cold regions cover a large proportion of snowmelt-dominated river basins (>45°N and S, with some exceptions) over the globe. These regions are particularly sensitive to ongoing climate change (precipitation and temperature) with regard to river discharge. Earlier runoff in the spring or winter and decreased flow during summer and autumn have been observed/projected under the circumstance of increasing temperature [45]. In this regard, long-term and NRT snowmelt information is essential for hydropower and hazard management. With EO it is possible to not only retrieve annual cyclic snow cover parameters (i.e., SCD, SCS, SCM; see Section 3.2.1), but also to monitor the snowmelt processes.

Conventionally, wet snow can be mapped using multi-temporal SAR observations. In principle, multi-temporal SAR based algorithms utilize the significant seasonal variation of snow surface backscatter due to changes of snow liquid water content and variable surface roughness. Typically, the backscattering coefficient decreases as snow wetness increases. Accordingly, Koskinen et al. [235] and Pulliainen et al. [236] introduced the Helsinki University of Technology (TKK) SCA method. The wet snow cover area (SCAws) can be calculated based on observed backscattering (${\sigma }_{obs}^{°}$), backscattering of reference surface (${\sigma }_{ref}^{°}$, dry snow or bare ground), and backscattering of wet snow (${\sigma }_{ws}^{°}$) according to Equation (9).

$$SC{A}_{ws}= \frac{{\sigma }_{ref}^{°}-{\sigma }_{obs}^{°}}{{\sigma }_{ref}^{°}-{\sigma }_{ws}^{°}}$$

(9)

The algorithm was validated by Luojus et al. [237] using model simulations from the Watershed Simulation and Forecasting System (WSFS); meanwhile the authors have calibrated the model using ground truth data. The results demonstrated the better performance of the TKK SCA method in forested areas compare to open areas. Slightly different, Nagler and Rott [238] developed a change detection algorithm for wet-snow detection in the Ötztal and Leutasch regions (Austrian Alps) using repeat-pass SAR backscattering. The author formulated a ratio (${\sigma }_{ws}^{°}/{\sigma }_{ref}^{°})$ between the backscattering coefficient of wet snow and a reference image, pixel by pixel. In the original paper, the authors suggested a threshold of −3 dB for the ratio image to separate wet snow from peripheral land-cover in Radarsat and ERS SAR data. Recently, Nagler et al. [239] updated and verified the backscattering ratio based algorithm for mapping snowmelt areas in Iceland and the European Alps using Sentinel-1 data. The validation revealed good consistency between Sentinel-1 snow maps and Landsat-based snow maps during summer (reference data). A slight trend of underestimation of snow extent was found in Sentinel-1 snow maps.

To determine the date of snowmelt onset, the temporal variations in backscattering and/or brightness temperature can be employed. For instance, Rawlins et al. [240] derived the timing of snowmelt in the pan-Arctic during the spring of 2000 using daily QuikSCAT-SeaWinds backscatter time series. The authors identified the primary snowmelt date as the day (during 60 to 182 Julian day) on which the difference between the daily backscattering and the previous 5-day mean reached a maximum. For multi-temporal snowmelt retrieval, Rotschky et al. [241] analyzed the spatiotemporal pattern of snowmelt-refreeze cycles in Svalbard using 9-year (from 2000 to 2008) QuikSCAT time series. For this purpose, the authors determined when the local backscatter dropped below a threshold that was calculated by multiplying mean winter (January–March) backscatter by an empirical factor c. Combing active (QuikSCAT) and passive (SMM/I) microwave satellite data, a study of snowmelt onset detection between 2000 and 2009 was conducted by Wang et al. [242] using previously published algorithms for the pan-Arctic . The study showed a positive trend (24.5 m·d⁻¹) between the mean snowmelt onset date and elevation for approximately 58% of the terrestrial Arctic.

Worldwide, to simulate and predict streamflow in snowmelt-dominated drainage basins, the Snowmelt Runoff Model (SRM) [243] is the most acknowledged model. Martinec [244] initially developed the SRM for small European river basins, and it has subsequently been successfully applied worldwide [243]. Swamy and Brivio [165] explored the potential of optical EO data (Landsat data set) together with ground meteorological and hydrological information to improve snowmelt runoff modelling for the Italian Alps. They conducted a linear regression analysis revealing a high correlation (r = 0.94) between simulated and measured discharges. Nagler et al. [245] presented a practical data assimilation scheme for short-term runoff forecasts in Alpine basins exemplified by the Ötztal drainage basin (Austrian Alps). The scheme assimilates the SCA retrieved from EO data and meteorological observations into a semi-distributed runoff model. During the model test runs, the authors found that the resultant runoff simulations were more influenced by meteorological-forecast-induced errors than was the satellite-derived snow cover information.

3.2.5. Operational Snow Products

Since 1966, weekly binary snow cover charts of the Northern Hemisphere have been provided by NOAA/NESDIS based on manual inspection by trained meteorologists using optical satellite imagery [223,246,247]. After the successful development of IMS in 1997, snow extent maps can be generated at daily intervals in a more efficient and accurate manner [4,5]. Such products are usually only suitable for studies at a global scale due to their medium resolution (1 km introduced in 2014). In December 1999, MODIS-Terra was launched. One year later, the MODIS-Terra 500 m daily snow cover product (MOD10A1) became available through the National Snow and Ice Data Center (NSIDC). The product is generated using the SNOMAP algorithm [166,167]. During the Globsnow-1/2 project [10], ESA also developed its SCA and SWE products [248] based on ERS-2/ATSR-2 (1995–2003) and ENVISAT-AATSR (2002–2012) data. For more snow product information and data access, we suggest that readers explore what is offered by the National Snow and Ice Data Center (NSIDC) [249]. Given the existence of various snow-related products, intercomparison and validation are essential. Toward this end, ESA launched the SnowPEX (Satellite Snow Product Intercomparison and Evaluation Exercise) project [250] within the Quality Assurance framework for Earth Observation (QA4EO) to intercompare and validate the major satellite snow products for the Northern Hemisphere and the globe. It is expected that the results will be published in the near future.

At present, snow products are usually derived from monotype EO data (e.g., optical or PM data). However, a practical, data-independent assimilation scheme is desirable as it holds the potential to combine the advantages of different types of EO sensors while avoiding data-inherent disadvantages (e.g., cloud obstruction, snow/ice penetration, coarse resolution). For such purpose, Foster et al. [171] presented a 25-km blended global snow product (ANSA) including SCA, FSC, SWE, onset of snowmelt, and snowmelt area. It is a combination of AMSR-E, MODIS and QuikSCAT data products. In Europe, Solberg et al. [251] introduced a confidence index (also known as a quality index). The confidence index is a function of time and single-sensor quality indicators for each pixel. A single-sensor quality indicator therein is a function of acquisition geometry, retrieval algorithm reliability, etc. Ultimately, the “best pixel” can be identified by the maximum confidence index among the multi-type sensor products pixel-wise. To provide snow monitoring service in Scandinavia and the European Alps, the concept has been applied in the EO-HYDRO project within ESA’s Earth Observation Market Development (EOMD) program [252]. Later, Solberg et al. [253] presented and summarized the confidence index based multi-sensor algorithm for snow cover and snowmelt monitoring for Norway and Sweden using MODIS-Terra and ENVISAT -ASAR data time series.

Major findings with regards to EO-based snow studies in Europe were mainly pursued by Dietz et al., Di Mauro et al., Rotschky et al., Wang et al., Foster et al., Dedieu et al., Hüsler et al., and Dye [32,217,241,242,254,255,256,257]. It is agreed by most authors that:

• The pattern of snow cover duration is strongly latitude and altitude dependent. Iceland and Norway have the longest snow cover duration.
• Average snow cover duration is shortening (5 to 6 days/decade between 1972 and 2000 [257]) in Europe, which is particularly accentuated at low elevations (700–900 m.a.s.l.). Also, the snowmelt onset date is shifting towards early spring at the rate of 3 to 5 days/decade, projected by Dye [257].
• An increasing amount of impurities are contained in European snowpack, frequently in spring, which increased the surface radiation of snowpack and may result in earlier snowmelt.

To date, direct monitoring of permafrost from space remains challenging. Since permafrost covers only a limited area in Europe, the number of publications is limited. Based on our literature research, EO-based research on lake and river ice in Europe is limited as well. As to cold region-related-hazard research using EO, we found some recently published review papers, such as [258,259]. For these reasons, the above-mentioned topics are excluded from this article. Given their uniqueness and spatial distribution, boreal forest, tundra and living organisms will not be studied here. We focused instead on glacier and snow, as well as their interaction with the surrounding environment.

4. Discussion

In the previous sections, we have provided an overview of EO-based studies on cold region land surface dynamics in Europe. It is obvious that numerous challenges for EO-based derivation of cold region land surface dynamics in Europe still exist. They also define the pathway for future analyses, based on existing methodological and geoscientific research gaps.

4.1. Challenges with Respect to Spatiotemporal Scale and Study Area Settings

Often, EO-based derivation of cold region land surface dynamics in Europe is undertaken at very different scales. Glacier studies are mostly presented at an individual glacier scale, and only a limited number of studies focus on a regional scale (e.g., the European Alps, the Svalbard Archipelago). The time interval of EO-based glacier studies varies from monthly to decadal, while most of the studies are multi-temporal observations at an annual interval. Snow cover studies on the other hand can be found even on the hemispherical scale with much higher temporal resolution (e.g., daily). This is due to the fact that high spatial resolution is not so critical when analyzing large scale snow cover characteristics, while snow cover is highly variable in time and therefore requires high temporal resolution. The cold regions in Europe are situated mostly in high latitude areas near the polar circle and in mountainous areas. Typical problems regarding high latitude areas are the occurrence of polar darkness and frequent cloud cover, which has a great influence on optical data availability. On the other hand, the total number of acquired observations increases owing to the increased overlap of individual paths of polar orbiting satellites at high attitudes. Also, some areas in high latitudes are simply not covered by EO data. In mountainous areas, the complex and steep terrain is problematic for SAR image processing.

4.2. Challenges with Respect to Data Availability, Costs and Suitability of Different EO Sensor Types

The selection of EO data is not only based on the extent of the study area, required accuracy and specified time span, but also on the suitability and availability of data as well as their cost. Very high-resolution optical and SAR data are generally quite expensive. Optical imagery of medium to coarse resolution (e.g., Landsat, ASTER, MODIS, AVHRR) is available for free and has been frequently used in cold region dynamic analyses. Among these sensors, geostationary sensor provides the highest temporal resolution, which could provide NRT information with regards to the general cold region land surface dynamics. In terms of SAR and LiDAR data, their availability is limited and data gaps often exist (e.g., between the termination of ENVISAT-ASAR and the launch of Sentinel-1). At a first glance, the day-night all weather operation of SAR sensors, makes them particularly intriguing for studying cold regions in Europe. However, ice/snow penetration, the complex topography and related geometric distortion effects in SAR data (layover, foreshortening, SAR shadows etc.) make this data not very suitable for the analyses in complex terrain. Generally, optical data have strong capabilities for retrieving almost all glacier- and snow-related parameters for cold region monitoring. Especially now with the Sentinel fleet in orbit (combination of Landsat with Sentinel-2A/B) the very high revisit time makes these data suitable for monitoring highly dynamic snow cover changes.

The benefits and limitations of optical and SAR data are further summarized in

Table 3. Benefits and limitations of optical imagery for cold region monitoring.

Aspect	Benefits	Limitation
	Good data availability and accessibility	No usable data under cloud cover and polar darkness
General	Good global spatiotemporal coverage	Relatively few spectral bands to distinguish snow/ice
		Atmospheric absorption and scattering
Glacier delineation	Fast, robust, semi-automated algorithm	Problematic identification in cases of debris cover
Glacier elevation change	Automated algorithm	High optical contrast is needed
	No snow/ice penetration	Limited number and accessibility of stereo-geometric sensors
Glacier mass balance	Usefulness for building time series	SLA does not always represent ELA
		Albedo only provides relative mass balance
Glacier motion	Long-term historic record	Error in image orthorectification processing and whisk-broom acquisition
Snow Cover Area	Fast, robust, semi-automated algorithm	Confusion between snow and dark targets (e.g., black spruce forest, water bodies)
Snow Cover Fraction	Fast, robust, straightforward algorithm	Proper endmember selection is required
		Limited number of bands in multispectral data
Snow Grain Size & Impurity	High sensitivity	Only for the upper layer of the snowpack
Snow Water Equivalent & Depth	N/A	N/A
Snowmelt	High temporal resolution	Cannot distinguish between wet snow and dry snow

and

Table 4. Benefits and limitations of Synthetic Aperture Radar (SAR) imagery for cold region monitoring.

Aspect	Benefits	Limitation
General	Day-night all-weather operation	Snow/ice penetration
		Complex format and big data size.
		Propagation effects in ionosphere
Glacier delineation	Debris cover is less problematic	Unable to detect glaciers based on single frequency amplitude/phase image
Glacier elevation change	Very accurate	Coherence requirement (InSAR & D-InSAR)
		Penetration correction is necessary
Glacier mass balance	N/A	N/A
Glacier motion	Most accurate (D-InSAR) and very	Need for retained coherence (Coherence tracking and D-InSAR)
	accurate for offset-tracking	Only slant-range motion (D-InSAR)
		Requires large image patches for incoherent intensity tracking
Snow Cover Area	N/A	Only able to detect wet snow
Snow Cover Fraction	N/A	N/A
Snow Grain Size & Impurity	N/A	N/A
Snow Water Equivalent & Depth	Physical model based method	Mostly single frequency sensors
Snowmelt	Can detect wet snow	N/A

for optical and SAR sensors, respectively. The only spaceborne LiDAR sensor, ICESat-GLAS, provided data during 2003–2009. It has a relatively large footprint (70 m) and ground spacing (170 m) along track, which is problematic for monitoring small glaciers. Nevertheless, the value of ICESat-GLAS observations as a calibration and validation data set should be emphasized. In general, PM data are less influenced by the atmosphere (apart from precipitation clouds); they offer high temporal resolution at global coverage. However, the spatial resolution of PM data is coarse, which makes this data type less suitable for regional studies. In terms of snow monitoring, it remains the major data set used in retrieving SWE. There are difficulties in using PM data to detect: (1) wet snow cover, and (2) thin snow cover (<5 cm snow) in heavily forested areas. In addition, PM data is strongly influenced by different snowpack physical properties. Underutilization of the recently available sensors is recognized. For example, there are only a few studies using Sentinel-2 for snow monitoring in Europe. For monitoring glaciers in Europe, researchers so far rarely used Sentinel-1, Sentinel-2 and TerraSAR-X data. Also, existing studies relying on these sensors often focus on methodological transferability and neglect some key features of the new sensor, such as the abundant bands on Sentinel-2. However, given that these sensors are only in orbit for a relatively short time, it is expected that studies in this field will be published soon.

4.3. Challenges with Respect to Method Applicability

Several different methods have been applied for the EO-based derivation of cold region land surface dynamics in Europe at different scales. For each topic, some methods have been widely applied, e.g., ratio-image thresholding for glacier delineation and snow cover area mapping, InSAR and offset-tracking for glacier motion monitoring, DEM subtraction for estimating glacier elevation/volumetric change, the “geodetic method” for glacier mass balance estimation, spectral unmixing and modelling for snow cover fraction derivation, empirical relationship for snow water equivalent retrieval, and snow runoff models for snowmelt prediction. Apart from the empirical algorithms developed for snow characterization (limited by dependency on field data), the other algorithms have been widely applied and validated in different regions at different scales. However, the aforementioned methods are usually based on a specific satellite sensor only. Their application may suffer from some sensor-inherent drawbacks. Currently, fully transferable multi-sensor-based methods, which combine the advantages of several sensors, or which have proven their easy transferability in space and time are still relatively rare.

Because existing semi-automatic algorithms for glacier delineation are mainly based on optical data, there are significant drawbacks in mapping glacier debris cover, as automatic algorithms cannot distinguish between this debris and the surrounding terrain. The required post manual correction obstructs algorithm automation. Combining novel InSAR coherence based algorithms with such optical-based algorithms thus hold a great potential for a fully automatized glacier delineation procedure. In terms of snow cover area monitoring, the most well recognized problem is snow cover detection in forested areas, which has a significant influence on the accuracy of snow cover product. Furthermore, the lack of a large-scale reference/validation data set complicates the accuracy assessment for snow cover products. Another challenge for snow-related studies is SWE and SD retrieval, particularly with regard to accuracy and spatial resolution.

5. Conclusions and Outlook

The cold regions in Europe belong to the most climate sensitive areas of the continent. Concern about and interest in cold region land surface dynamics in Europe have risen considerably in recent years and are expected to increase in the future. Earth Observation (EO) provides a unique data source for objective, repeated, synoptic, consistent, long-term, and large-scale analysis of cold region land surface dynamics. Especially EO-derived Snow Cover Duration (SCD) is a useful indicator for delineating changing cold region boundaries in a quantitative way. We reviewed over 250 SCI journal papers from the past 30 years for this area. An immense number of studies have been carried out in order to develop accurate and transferable methods for the retrieval of cold region land surface dynamics. Amongst these studies, Norway (including the Svalbard Archipelago) and the European Alps are most frequently investigated. Optical sensors are the backbone for this task owing to their suitability for cold region monitoring. Also, SAR data are increasingly employed for EO based analyses. In this domain, ratio-image thresholding, image offset-tracking, and InSAR are the most frequently applied methods.

Based on EO, there is an overall agreement and understanding that in Europe: (1) glaciers are shrinking; (2) glacier elevation is lowering and glaciers are experiencing mass loss; (3) glacier surging is especially accentuated within the Svalbard Archipelago; (4) snow cover duration is shortening; and (5) the snowmelt onset date is shifting towards earlier times in spring.

However, numerous challenges remain in analyzing cold regions in Europe, such as:

• Delineating glaciers in a fast and fully automated manner is difficult;
• Correction for snow/ice penetration of SAR data with regard to snow/glacier monitoring is desirable;
• Mapping snow in forested areas, and precisely retrieving Snow Water Equivalent (SWE) and Snow Depth (SD) determination at high to medium resolution remain challenging;
• Further studies on synergetic multi-sensor EO data application in cold regions is desired;
• Key features of recently available sensors (e.g., Sentinel fleet) are so far underexplored;
• Suitable validation data sets for assessing snow cover and glacier extent are scarce;
• The pattern of glacier dynamics at a meso-/continental and long-term scale is unclear;
• In some areas, regional responses of snow distribution to ongoing climate change remain uncertain;

Future studies should not only focus on the descriptive and quantitative retrieval of cold region land surface dynamics in Europe, but also be directed towards a better understanding of these dynamics, especially as a response to ongoing climatic changes. It is expected that great benefits arise from the novel European Sentinel fleet, with Sentinel-1 A/B, Sentinel-2 A/B, and Sentinel 3A already in orbit. In combination with 30 years of Landsat data, open ENVISAT archives, and additional expected missions such as NASA’s ICESat-2, envisaged for launch in 2018, or Germany’s envisaged Tandem-L, cold region studies for Europe as well as globally will be strengthened.