Mountain glaciers analyzed by remote sensing

Elise Colin
16 min readMar 1, 2022

Like last year, 31 first-year students at CentraleSupélec worked for a week on remote sensing of glaciers. This year, they focused on the Guil watershed in the French Alps, the Karakoram range which span the border of India, Pakistan and China and Southeastern Alaska in the United States. They studied the type of snow cover, surface glacier velocities and glacier surface textures for the detection of surges and the study of intraseasonal changes.

Students : Guillem Carcanade, Alban Wales, Antoine Manyeres, Jérémie Sabot, Tommy Champon, Gaspard Borderon, Yann Bruneel, Xavier Fiat, Lauriane Mousset, Marie Saubusse, Benjamin Vissa, Kim Leconte, Rémi Forasetto, Maggy Zhou, Elliot Hernandez, Alexandre Dessagne, Paul Lardillier, Racim Menasria, Théo Cavina, Gauthier Debuisschert, Louise Mui, Etienne Vaquier, Henri Coquinot, Adrien Gleize, Matiss Péchoultres, Nicolas Barbarin, Flore Bergerot, Pauline Paquereau, Martin Mollicone, Paul Tabbara, Ilian Kara-Mostefa

This project was proposed by Laurane Charrier, a PhD student from LISTIC, University Savoie Mont-Blanc and Onera, assisted by Flora Weissgerber, Nathan Letheule and Elise Colin at Onera.

Summary of the work

First, a seasonal version of the “REACTIV” code has been developed (REACTIV_season), it seems to be very useful for the detection of surges in mountainous regions.

Second, since more and more glacier velocity estimation data are available online (ITS_LIVE, RETREAT, etc), we have shown that these data are useful to provide a global overview of glacier dynamics thanks to user-friendly websites.

Third, the evaluation of snow cover in optical images is temporally restricted by cloud cover. In radar, it is very interesting but the main track of future improvements will consist in separating vegetation areas from bare ground areas.

This study highlights the complementary between the analysis of glacier surface velocity and radiometric changes as well as the complementary between optical and SAR images.

Methods

Ice velocities

There are several methods for calculating glacier displacements from two images taken on two dates. We were able to implement one of them, developed at Onera: GeFolki. This method allows us to calculate the optical flow, or deformation field, existing between two images, in each pixel.

For the analysis of glacier velocities, we also used existing databases: ITS LIVE [gardner et al., 2018] (https://its-live.jpl.nasa.gov/) , that provides a comprehensive and temporally dense multi-sensor record of land ice velocity, and RETREAT [Friedl et al., 2021b] (http://retreat.geographie.uni-erlangen.de/search) which provide glacier surface velocities from Sentinel-1 images.

Study of SAR radiometric changes: Leclerc and REACTIV methods

Two complementary methods implemented on the Google Earth Engine platform were used: the Normalized Difference Index (NDI) from [Leclercq et al., 2021] and the Reactiv algorithm [Colin Koeniguer et al., 2019]. Both techniques aim to detect the increase in radar backscatter related to the opening of crevasses on the glacier surface, the roughness becoming more important. The opening of crevasses, linked to the increase in glacier surface velocity, is an indirect indicator of glacier surge.

[Leclercq et al., 2021] compute a normalized difference between aggregated images of maximum backscatter over the winter periods of two different years. If the backscatter increases between two years, the index will be negative meaning that the glacier is likely at the end of a surge. On the contrary, if this index is positif, the glacier may be at the beginning of a surge.

The REACTIV method [Colin Koeniguer et al., 2019] allows to visualize the evolution of the glacier with more than 3 images by associating the hue channel to the corresponding acceleration periods. The method considers a whole time series, and it looks if the profiles correspond to the realizations of the same speckle over time, potentially including a deterministic component. This detection is done using the coefficient of temporal variation. The saturation, null for the theoretical average value of this coefficient, increases progressively with it as one moves away from a homogeneous temporal population. Thus, in the images produced, the colored glaciers are potential glaciers in surge.

Following the project week, a group of students initiated a seasonal version of REACTIV: in this one, only images restricted to the winter seasons are taken into account, to get rid of the natural seasonal variations. This version will be called REACTIV_season.

Study of snow cover: Nagler Method, Lievens Method and NDSI

The Normalized Difference Snow Index (NDSI) snow cover is an index that is related to the presence of snow in a pixel. It is computed by the ratio of the difference in VIS and SWIR reflectance from multispectral optical images. A pixel with a positive NDSI is considered to have some snow present whereas a pixel with a negative NDSI is a snow free land surface [Riggs et al., 2016].

In SAR images, two kinds of snow have very different behaviors: wet snow and dry snow. Dry snow is composed of solid and gaseous molecules whereas wet snow is also partly composed of liquid water. This distinction is important because the presence of wet snow reveals a snowfall in the last few days or snowmelt. Using polarimetry, we can exploit the difference of scattering of the wave between wet and dry snow: the wave propagates into dry snow and its polarization changes, while it is reflected on wet snow. From dual pol VH and VV radar data, the [Nagler et al., 2016] method is a method that allows the detection of the presence of wet snow. The [Lievens et al., 2019] method is a method that allows us to calculate the evolution of a dry snow depth over time, from a reference image, using the Snow Index (SI, in dB).

The Guil basin

The Guil basin in the Hautes-Alpes is a mountainous area where the peaks can reach 3000 m. The Guil, a tributary of the Durance, flows into the artificial lake of Serre-Ponçon equipped with a hydroelectric dam, and is fed by the melting snow in the region. It is thus important to be able to calculate the snow cover of the basin and to estimate the volume of snow at any time of the year.

We therefore applied the Nagler method [Nagler et al., 2016] , to detect wet snow. A new collection of binary images that reflects the evolution of melt over the year was created. Then, we calculated the dry snow depth using the Lievens method [Lievens et al., 2019]. This method is based on an estimation of the volume backscatter, which can come either from the dry snow or from the vegetation in which the radar waves penetrate. We could thus observe that if the vegetation is not separated from the snow by external data, the accumulated snow height over the basin does increase in autumn but does not decrease in May/June when melting is observed.

In order to obtain a snow mask, on which to apply the [Lievens et al., 2019] method, we computed snow masks thanks to the thresholding of the NDSI index computed from the Sentinel-2 data, and other criteria proposed from MODIS [Hall et al., 2021]. The main limitation of this measurement is the cloud cover, so we had to restrict the dates of exploitable acquisitions. Nevertheless, we observed an initial gradual increase in snow cover from October 2018 to February 2019 until reaching 80% of the basin surface, before observing the slow disappearance of snow. To better understand the limits of optical images, and to detect recurrent observation or interference windows, we also investigated some statistics on cloud cover.

Karakoram

The Karakoram range, situated between China, India and Pakistan is covered by 28 to 50 percent of ice. This region is special since the glacier mass change was slightly positive from 2000 to 2016 [Brun et al., 2018], a phenomenon called the Karakoram anomaly. Moreover, 45.3 % of the Karakoram glaciers correspond to surge-type glaciers [Guillet et al., 2021] . A surge-type glacier is a glacier which experiences an increase in flow velocity by a factor of 10 to 100 over a short period of time ranging from a few months to a few years.

These surge events represent a risk for the populations of the valleys when ice-dammed lakes formed downstream of the glaciers overflow. Such an event was observed in August 10, 2018 for the Kyagar glacier for example. The corresponding flood discharge was 1570 m/s. However, the control mechanisms of the surges are not completely understood, the techniques of identification and study of the surges are still interesting to develop. The aim of the study is to work on methods of identification of surge-type glaciers over a large scale based on radar and optical radiometric changes compared with estimates of glacier velocities.

Identification of surges through the analysis of radiometric changes over time

From the image obtained by the method of Leclerc, we have developed a binary mask to quickly identify areas of strong backscatter corresponding to glaciers in surge. Our aim is to set up an automatic detection of these areas where the current method still requires a visual work of identification on the image obtained.

Mask obtained by the Leclerc Method to identify rapid backscattering changes
The REACTIV_season algorithm, applied here on the Karakoram chain, allows to quickly identify particular change phenomena, and the years in which the main change is recorded.

The Reactiv method has proven to be more sensitive and has allowed us to identify a greater number of glaciers on the surface. It is also faster since it analyzes the evolutions over four years at the same time. On the other hand, it does not allow to differentiate between acceleration and deceleration and therefore remains complementary to the Leclercq method. These two methods, which can be automated, allow a rapid and large-scale identification of surge phenomena in a mountain range.

Then, we completed our analyses locally with optical images. We created timelapses over several years from EO-Browser false color compositions (B08-B04-B02) to facilitate the identification of ice. These timelapses confirm or not the considered glacier surges . In addition, a pointing carried out on each glacier using Tracker makes it possible to obtain the profile of speed and position of the terminus of the glacier and thus to determine the position in time and the intensity of the surge. We were able to attest the surges in particular of the Yanatsugat, South Rimo and Khurdopin glaciers.

The combination of these methods thus allows a fast and global overview of surge-type glaciers at the scale of a mountain range. However, it remains limited by the quantity of radar and optical images (especially before 2015), by the meteorological conditions and by the important noise of some SAR images which can complicate the automatic detection by the algorithms that we have started to develop.

Glacier Speed

Velocities computed over the Kyagar glacier.

Velocities were mainly studied over the Kyagar glacier. The Gefolki optical flow algorithm was applied to obtain velocity maps of the glacier between each consecutive month throughout 2015. Once the maps were obtained, we first checked the precision of our results by comparing the velocity in the chronological and anti-chronological order. The median absolute difference between chronological and anti-chronological velocities is on average 0.0007 m/d in 2015.

Visual comparisons as well as graphs characterizing the difference in velocity reveal that the GeFolki results and data from Retreat or [Round et al.,2017] are roughly similar, although GeFolki finds slightly lower values. The highest differences are observed at the glacier borders and over feature-less regions, characterized by a flat slope and a permanent snow cover.

With every method, we observe a zone on the bottom right end side of our study area showing very high velocity magnitudes. To determine the reasons for this, we first designed an HSV map of the velocity directions over the year. The displacement direction was opposite to the one expected according to the topographic slopes. Georeferencing our velocity images and projecting them onto the glacier allowed us to compare our results with a topographic profile. This analysis confirms that the observed velocity zone is likely an artifact. This analysis highlights the interest to use velocity direction and the topography to identify artifacts in velocity maps.

Analysis of Retreat and ITS_LIVE data

Surface velocities for Rimo glacier on the left and Miar glacier on the right. The datasets are from RETREAT on the top and ITS_LIVE at the bottom.

The study of the Kyagar Glacier on the Retreat platform has made it possible to highlight the period of surge that the glacier underwent until June 2015 by studying the evolution of the longitudinal monthly velocity profile over the years. However, because Sentinel-1 was launched in April 2014, the start date of the surge remains unknown. Using Landsat-8 scene pair velocities from the ITS_LIVE platform is complementary: the period of surge appears clearly, it extends from the beginning of 2014 to early 2016.

By repeating this analysis on several Karakoram glaciers, we identified glacial surges in recent years.. In order to detect them over a large scale, we have set up an identification algorithm which takes in input velocity profiles from the Retreat platform. The detection is based on a threshold of variation in longitudinal velocity averaged over three years, which detects the presence or not of a surge, and if so, its start and end dates. Although imperfect due to the singularity of the characteristics of each glacier, this algorithm represents a first approach in the processing of Retreat data. In addition, we compared the evolution of the velocities obtained on Retreat and ITS_LIVE over the same point on the Kyagar glacier: the trends are globally similar, and the surge appears clearly. The surrounding glaciers highlight several typical profiles: the Miar glacier, whose velocity profile is quasi-sinusoidal (no surge), and the Rimo glacier, with a significant peak (marked surge). We note that ITS_LIVE also allows us to estimate the position of slope breaks on a glacier, and to quickly determine the areas where velocities are high thanks to annual velocity mosaics.

Snow Cover Analysis

Velocity magnitude over a point on the Rimo glacier based on ITS_LIVE dataset

First, an analysis of radar images was performed on an area (Kangri glacier) close to a river and without any vegetation according to optical NDVI studies. We quantified the surface area of wet snow and the volume of dry snow as a function of time. We then correlated these data with the mean temperature fluctuations over the year in the Karakoram area, in order to observe seasonal variations.

A second area near the Kyagar glacier raises our attention. We observed a drop in the radar signal, indicating the appearance of wet snow, certainly due to melting, and then the reappearance of dry snow, over a few images and therefore over only a few weeks. During the melting episode, we could observe the propagation of the decrease of the radar signal from the bottom of the glaciers to the summits. We hypothesized a strong heat peak to explain this melting. But the absence of meteorological data and of exploitable optical images which could have allowed us to observe, possibly, the formation of lakes resulting from the melting, did not allow us to confirm this hypothesis. The study of this area over several years has confirmed that this melt is a cyclic episode that repeats itself every year.

We then carried out analyses on accumulation zones upstream of the Kyagar and Rimo glaciers with the aim of illustrating a correlation between the velocity of the glaciers studied previously, and the quantity of wet snow. The graphs obtained show a particularly large amount of wet snow during the summers of 2017 and 2018 (i.e. before the start of the glacier surge) in the case of the Rimo glacier.

This observation could be in favor of [Jiang et al 2021] assumption that the Rimo surge of 2019 is hydraulically controlled.

Alaska, Malaspina Glacier

REACTIV_by_season applied to Alaska, reveals strong variations on the Agassiz Lake
REACTIV_season applied to Malaspina Glacier

Alaska glaciers have lost 66.7 Gt/yr between 2000 and 2019, which corresponds to 25% of the world glacier mass loss over this period [Hugonnet et al., 2021]. Moreover, this region is characterized by a high number of glacier surges. Hence, Alaskan glacier temporal changes are crucial to study. Here, we focus our attention on Southeastern Alaska and particularly on Malaspina, Seward and Valarie glaciers. Malaspina glacier is the world’s largest piedmont glacier covering roughly 2200 km² of the coastal foreland. It is partially fed by Seward glacier which is a surge type glacier. The two form the complex Malaspina-Seward. Valarie glacier is a tributary of the Hubbard glacier, the world’s largest non-polar tidewater glacier.

A first qualitative analysis of the glacier and its main evolutions was carried out thanks to optical images in summer, which allowed us to observe the displacement, crevasses and to highlight the glacier surface flow with the lines of sediments of Malaspina glacier.

Leclerc’s method was applied to visualize changes between the years 2018 and 2019. It was used to highlight, around the glacier lobe, the shifts in sediment lines and flow beds.

In order to visualize the temporal evolutions, we used the REACTIV code. To study these evolutions between the different winter periods from year to year, we derived a version of the algorithm, REACTIV_season, allowing us to select a fixed season through the different years of the studied period. We have adapted the code to select and compare several consecutive winters, from winter 2017–2018 to winter 2020–2021, without considering summers. Our analyses focused around Lake Agassiz, the site of changes that caught our attention. The code also highlighted the movement of sediment over several years, marked by lines of different colors for each winter.

Lake Agassiz shows colored areas of change for colors associated with radiometric maxima reached in winter 2019–2020. We hypothesized that these areas could correspond to frozen surfaces covered with dry snow, corresponding to strong backscatter, while they would be covered with wet snow for other winters. These assumptions were confirmed by the snowpack study, where a dry snow was detected on this lake only during the winter of 2019–2020.

In parallel, the study of the snowpack and its evolution over five years has been the subject of a particular study, both using the NSDI index in optical, and using methods dedicated to SAR images for dry snow and wet snow. Seasonality effects are well detected: dry snow depth is greater in winter than in summer. Moreover, a linear regression performed on the pooled data showed a decrease in snow depth from year to year, which could be a sign of global warming.

On the other hand, a dynamic study of the glaciers was performed. Since we already have velocity maps using the ITS-LIVE method (NASA), we first tried to obtain similar maps by applying the GeFolki optical flow algorithm to pairs of LandSat8 images spaced 16 days apart, with 15m pixels, for two periods: end of April and end of August. Technically, the main difficulty lies in the use of different reference frames for the ITS-LIVE images and the LandSat images on which we calculated the flow. After passing in the reference frame and residual correction of biases between projections, we were able to generate the differences between the images obtained. Both methods lead to similar results with a relative uncertainty on the values estimated at 2.7%.

Finally, the study of radar backscatter changes of the Valerie glacier over four winters, has highlighted an increase in VH polarization backscatter between winter 2019–2020 and winter 2020–2021. In addition, the glacier velocity obtained using ITS-LIVE and RETREAT data increases strongly over the same period. The Malaspina glacier is also accelerating: its maximum speed is 1800 meters per year between 2014 and 2020 against more than 3000 meters per year in 2021. This finding suggests that a surge upstream of the two glaciers could have occurred and propagated; probably in the Mt. Augusta area. An analysis of the glacier velocity in this area and a map of the evolution of the snowpack were therefore carried out. However, these analyses have shown a stability for several years and have therefore contradicted the hypothesis of a correlation between the two phenomena. This point would need further investigations.

Bibliography

[Lievens et al., 2019] Lievens, H., Demuzere, M., Marshall, H. P., Reichle, R. H., Brucker, L., Brangers, I., … & De Lannoy, G. J. (2019). Snow depth variability in the Northern Hemisphere mountains observed from space. Nature communications, 10(1), 1–12.

[Nagler et al., 2016] Nagler, T., Rott, H., Ripper, E., Bippus, G., & Hetzenecker, M. (2016). Advancements for snowmelt monitoring by means of Sentinel-1 SAR. Remote Sensing, 8(4), 348.

[Hall et al., 2021] Hall, D. K., Riggs, G. A., Salomonson, V. V., Barton, J. S., Casey, K., Chien, J. Y. L., … & Tait, A. B. (2001). Algorithm theoretical basis document (ATBD) for the MODIS snow and sea ice-mapping algorithms. Nasa Gsfc, 45.

[Tsai et al., 2019] Tsai, Y. L. S., Dietz, A., Oppelt, N., & Kuenzer, C. (2019). Remote sensing of snow cover using spaceborne SAR: A review. Remote Sensing, 11(12), 1456.

[Charrier et al., 2020] Charrier, L., Godet, P., Rambour, C., Weissgerber, F., Erdmann, S., & Koeniguer, E. C. (2020, September). Analysis of dense coregistration methods applied to optical and SAR time-series for ice flow estimations. In 2020 IEEE Radar Conference (RadarConf20) (pp. 1–6). IEEE.

[Colin Koeniguer et al., 2019] Koeniguer, E. C., & Nicolas, J. M. (2019). Change detection in SAR time-series based on the coefficient of variation. arXiv preprint arXiv:1904.11335.

[Leclercq et al., 2021] Leclercq, P. W., Kääb, A., & Altena, B. (2021). Brief Communication: Detection of glacier surge activity using cloud computing of Sentinel-1 radar data. The Cryosphere, 15(10), 4901–4907.

[ITSLIVE] https://github.com/nasa-jpl/its_live

https://its-live.jpl.nasa.gov/

“Velocity data generated using auto-RIFT (Gardner et al., 2018) and provided by the NASA MEaSUREs ITS_LIVE project (Gardner et al., 2019).”

[Gardner et al., 2018] Gardner, A. S., Fahnestock, M. A., Agram, P. S., Scambos, T., Nilsson, J., Paolo, F. S., … & Dehecq, A. (2018, December). ITS_LIVE: A new NASA MEaSUReS initiative to track the movement of the world’s ice. In AGU Fall Meeting Abstracts (Vol. 2018, pp. C14A-02B).

[Riggs et al., 2016] Riggs, G. A., Hall, D. K., & Román, M. O. (2016). NASA S-NPP VIIRS Snow Products Collection 1 User Guide.

[Friedl at al., 2021a] Friedl, P., Seehaus, T., & Braun, M. (2021). Global time series and temporal mosaics of glacier surface velocities derived from Sentinel-1 data. Earth System Science Data, 13(10), 4653–4675.

[Friedl et al., 2021b] Friedl, P., Seehaus, T., & Braun, M. (2021, April). RETREAT: A new freely available data set of Sentinel-1 glacier velocities in regions outside the polar ice sheets. In EGU General Assembly Conference Abstracts (pp. EGU21–2740).

[Jiang et al., 2021] Jiang, Z., Wu, K., Liu, S., Wang, X., Zhang, Y., Tahir, A. A., & Long, S. (2021). Surging dynamics of South Rimo Glacier, Eastern Karakoram. Environmental Research Letters, 16(11), 114044.

[Brun et al., 2017] Brun, F., Berthier, E., Wagnon, P., Kääb, A., & Treichler, D. (2017). A spatially resolved estimate of High Mountain Asia glacier mass balances from 2000 to 2016. Nature geoscience, 10(9), 668–673.

[Guillet et al., 2021] Guillet, G., King, O., Lv, M., Ghuffar, S., Benn, D., Quincey, D., & Bolch, T. (2021). A regionally resolved inventory of High Mountain Asia surge-type glaciers, derived from a multi-factor remote sensing approach. The Cryosphere Discussions, 1–31.

[Hugonnet et al., 2021] Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., … & Kääb, A. (2021). Accelerated global glacier mass loss in the early twenty-first century. Nature, 592(7856), 726–731.

--

--

Elise Colin

Researcher with broad experience in Signal and Image processing, focusing on big data and IA aspects for Earth observation images, and medical images.