Supplement of Brief communication: Detection of glacier surge activity using cloud computing of Sentinel-1 radar data

For studying the flow of glaciers and their response to climate change it is important to detect glacier surges. Here, we compute within Google Earth Engine the normalized differences between winter maxima of Sentinel-1 C-band radar backscatter image stacks over subsequent years. We arrive at a global map of annual backscatter changes, which are for glaciers in most cases related to changed crevassing associated with surge-type activity. For our demonstration period 2018–2019 we detected 69 surging glaciers, with many of them not classified so far as surge type. Comparison with glacier surface velocities 5 shows that we reliably find known surge activities. Our method can support operational monitoring of glacier surges, and some other special events such as large rock and snow avalanches.


Global results
Here we present in more detail the global results of our detection of glaciers with surge activity. Fig. S1 shows the global coverage of the Sentinel-1 NDI images for the period 2018-2019 as calculated from the available data in Google Earth Engine (GEE). To be able to calculate such a difference image we need radar backscatter imagery in both winters (defined as January -March for the Northern Hemisphere and June -August for the Southern Hemisphere), 2018 as well as in winter 2019. The VH polarized imagery in GEE covers a large fraction of the global land area. But northeast Canada, Greenland, the eastern edge of Svalbard, Franz Josef Land, Novaya Zemlya, and eastern parts of Severnya Zemlya are not covered by this in either ascending or descending path. In order to cover these glacierized area, we need to select the HV polarized images that gives us almost complete coverage of the Canadian Arctic and most of Greenland's local glaciers in descending path. There are no images over Antarctica. In ascending orbit we can complement with parts of Greenland and Svalbard. Both with VH and HV polarization, the descending orbit gives most coverage. So we use descending orbit as default and use ascending orbit images to complement for the areas not covered by descending orbit, and a few cases where the mountain topography makes ascending more suitable than descending orbit (see also Section 2). Combination of ascending and descending paths leaves some glacierized areas, such as Franz Josef Land and smaller areas in Greenland and Canadian Arctic, not covered for the 2018-2019 data. Where possible, we use NDI images from other years to fill data gaps.
In total we find surge activity for 69 glaciers and glacier tributaries in the study period 2018-2019. Table S1 provides a list of the surge-type glaciers we find, including the surge classification given in the RGI [RGI Consortium, 2017] based on Sevestre and Benn [2015]. Fig. S2 shows the distribution of these glaciers on a world map. The majority of the surges we have detected are found in Alaska (18 glaciers with surge activity), Svalbard (14) and High-Mountain Asia (13). In both the table and the figure, the glaciers are grouped based on whether the surge activity detection is based on increased backscatter (31 cases), decreased backscatter (29), or a combination of increased and decreased radar backscatter (9). In the table we have also included the 18 glaciers that show a change in backscatter but for which we found it unclear if this was caused by surge activity.
Of the 69 glaciers on which we have detected surge activity, 3 are classified as possible surge type (class 1) in the RGI, 5 are classified as probable surge type (class 2), 16 are classified as observed surging (class 3), 10 are classified as no evidence of surge (class 0), and 35 have no surge-type classification (class 9) in the RGI. This means that almost two-thirds of the glaciers for which we detect surge activity in our study period 2018-2019 have not been identified as surge-type glaciers previously. There are large regional differences in the classification of the surge-type glaciers we detect. Out of the 14 surging glaciers we found in Svalbard, 12 are classified in RGI as having either direct or indirect evidence of surge activity, but 16 of the 18 glaciers we find surging in Alaska do not have a surge classification in the RGI. 6 out of the 10 glaciers that we find to have surge activity while the RGI classification states no evidence of such activity are located in RGI region South Asia (West). Figure S1: Coverage (in grey) of NDI images based on Sentinel-1 radar images as available in Google Earth Engine for the periods used in this study, 1st January 2018 -1st April 2018 and 1st January 2019 -1st April 2019: a) VH polarization, descending orbit b) VH polarization, ascending orbit c) HV polarization, descending orbit d) HV polarization, ascending orbit. Figure S2: Global overview of the location of glaciers with surge activity found in this study. Glaciers for which the surge activity is derived from a backscatter increase between winter 2018 and winter 2019 are shown in red, a decrease in backscatter in blue, and glaciers which show both types of changes, an increase as well as parts with a decrease in backscatter, are shown in cyan. Table S1: GLIMS ID, RGI ID, name, region, latitude, longitude, and surge type as given in the RGI, for the glaciers where we have detected surge activity over the year 2018-2019 from Sentinel-1 radar images. Sometimes a local inventory code is given as a name, if both inventory code and a name are given in the GLIMS data set [GLIMS andNSIDC, 2015 (updated 2018)], we have only displayed the name. If no name is given, we have put a "-". Latitude and Longitude are the coordinates of the glacier centre point as given in GLIMS. Surge type classifies the surge activity from Sevestre and Benn [2015] as given in the RGI: 0 no evidence of surge; 1 surge possible; 2 surge probable; 3 observed surge; 9 not assigned. based on NDI image from years 2017-2018 * called "Hassanabad Glacier I" in GLIMS and RGI ** called "Hassanabad Glacier II" in GLIMS and RGI *** name in RGI is given as "Tupungato Sur/Tunuyan" **** a tributary to Walsh glacier ***** RGI name is SU5X14309227 Garmo † used as example in methods section and shown in Fig 2 Descending and ascending path The Sentinel-1 constellation is in a sun-synchronous orbit with and inclination of 98.18°and data are available from both ascending and descending path. As mentioned in the previous Section, we use mainly radar images from descending paths in the global study. The main reason to use the descending path is the better coverage of these images in GEE. In general we find no difference in the results between using images from ascending and descending paths where both are available, although in some particular occasions the two can provide complementing information. In Fig. S3 we show a comparison between the normalized difference index (NDI) images from descending and those from ascending paths. The first example in Fig. S3 (a-c) is the NDI image for the same area in the Pamirs as used in Fig. 1. As the line of sight is different, the valley walls that give a noisy signal and the ones that are in the radar shadow are different for the ascending and descending NDI images. But the NDI for the glacierized areas are very similar, such that the detection of surge activity is the same for the NDI derived from ascending path as it is for the NDI image from descending path. The second example ( Fig. S3 d-f) shows the NDI for both descending and ascending path over Shispher glacier and the Sentinel-2 image for the same area in the Karakoram, Pakistan. Shispher Glacier is known to have surged recently [Rashid et al., 2020]. The surging glacier tongue blocked the river from a tributary valley which led to an ice-dammed lake that threatened the downstream village of Hassanabad and the Karakoram Highway. Shispher glacier is surrounded by extremely steep terrain. The NDI images shown in Fig. S3 show a lot of noise around the highest peaks. On the descending NDI image shown in Fig. S3, the whole mountain flank to the west of the glacier tongue is very bright as well, while this part is not covered by ice. This mountain flank is under SAR foreshortening and layover that appears stretched out through the orthorectification process. Some artificial differences in this process create strong differences in the backscatter stack maxima. This example highlights that differences in radar backscatter brightness can be caused by acquisition and processing effects in addition to natural ones. In this specific case, the ascending geometry creates much less such problems and the surge can be detected much more clearly. The Shispher Glacier example also demonstrates how glacier lake changes can be detected using our method. Changes in the ice-dammed lake to the west of the glacier tongue are clearly visible as backscatter changes in the descending path (for the ascending path the lake lies in the radar shadow), based on the fact that smooth water surfaces create particularly strong forward-scattering, and thus reduced backscatter (red arrow in Fig. S3d). The noise pattern in the difference images from the ascending paths is clearly different from that of the descending paths. As the radar signal comes from the opposite direction, the effects of overlay, foreshortening and shading differ. This is most prominent in the example of Shispher Glacier, due to the extreme topography of the Karakorum where it is located. In the particular case of Shispher Glacier, the signal from the glacier surge is much better discernible from the acquisition and processing effects in the image taken in ascending path. However, for the vast majority of the detected surges, as in the three other examples shown here, the signal from the glacier surface is almost identical for both paths and the surge detection is not influenced by the choice of path azimuth. This also indicates that the surge detection is not dependent on the angle between direction of glacier flow and radar satellite azimuth, likely due to an over large parts chaotic nature of surge-related crevasses, without directional preference.  Fig. 1. The second to fourth example show the NDI images of Shispher Glacier, the glacier with GLIMS ID G077483E35705N, both in the Karakoram Range, and Kluane Glacier, Alaska, respectively. The red arrow in d points at the decrease in radar backscatter due to proglacial lake changes.

Other observed processes
As mentioned in the Discussion section, surge activity is not the only natural process that can cause temporal difference in radar backscatter brightness on glaciers. Although glacier surges are the main focus of the paper, we would like to stress that the method we propose here to detect glacier surges also has potential in other fields of geoscience. Further exploration of the method in terms of different stack statistics and other differencing periods than the yearly interval we have used could further improve the detection of a range of other processes. We have observed difference in backscatter brightness in accumulation areas where we suspect changing firn conditions to be the cause, changes in glacial lakes, change in sea-ice conditions, landslides, a glacier detachment, and snow avalanche activity. Figure S4 shows four examples of these processes.
The first image of the figure (Fig. S4a) shows the 2018-2019 normalized difference image where we observe increased backscatter from the surface of a nameless glacier in the Alaska Range, Alaska (GLIMS ID G213541E63316N) that is similar to the signal of a glacier surge. Studying the Sentinel-2 optical images of summer 2017 (Fig. S4b) and summer 2019 (Fig. S4c) it becomes clear that this increase in radar backscatter is caused by the deposits of a large landslide, rather than by surge activity. The landslide is first seen on the Sentinel-2 image from the 7th of July 2018, and must have happened between the Sentinel-2 acquisitions on the 6th and the 7th of July 2018. The landslide deposits have a higher roughness and thus backscatter than the underlying glacier surface, which leads to an increase in the backscatter brightness from winter 2018, before the event, to winter 2019.
The second example (Fig. S4 d-f) shows both a landslide and snow avalanches on Tlikakila Glacier Fork, Alaska Range. Again, the landslide causes an increase in radar backscatter brightness, and must have occurred between acquisition of the last radar images of winter 2018 and the acquisition of the summer 2018 Sentinel-2 composite image shown in Fig. S4f. Note that there is a band of no increase in backscatter where the landslide covers the middle moraine on the glacier, i.e. did not change roughness and backscatter. In addition to the landslide, we observed several small patches of decreased radar backscatter. These correspond to snow patches that are visible in summer 2018 ( Fig. S4f) but not in summer 2017 (Fig. S4e) or summer 2019 (not shown) and we believe to be the remainders of snow avalanche deposits. These patches of decreased radar backscatter can thus be explained by avalanche activity during the winter of 2018 that increased the stack maximum backscatter in this period on the avalanche deposits. During summer, the avalanche deposits melted and during the following winter there was less or no avalanche activity at these locations, such that the stack maximum backscatter in this period was lower than it was the winter before. In our global analysis we find a large number of snow avalanche signs, over glaciers and around.
The third example (Fig. S4 g-i) shows the observation of the glacier detachment of Sedongpu Glacier in eastern Tibet with GLIMS ID G094940E29811N that occurred in October 2018 [Kääb et al., 2021]. On the NDI image we see clear changes in radar backscatter on the lower part of the glacier and in the main valley where deposits of the glacier detachment blocked the Yarlung Tsangpo river. Comparison of two composite Sentinel-2 images, one from autumn (October -December, a season with less clouds than summer) 2017 and the second from autumn 2018, show that in autumn 2018 the lower part of the glacier has collapsed and that the river in the valley beneath has been filled with ice and rock debris.
A last example (Fig. S4 j-l) shows a decrease in radar backscatter between winter 2018 and winter 2019 over the firn area of Langjökull, Iceland. This decrease in radar backscatter is prominent on all Icelandic glaciers for this year, and changes in backscatter in glacier firn areas can be found in all regions. We believe these changes have to do with firn processes such as changes in the melt water content, formation of superimposed ice or ice lenses in the snow pack. Figure S4: Examples of features other than surges visible in the NDI images. Here we show a landslide on a glacier with GLIMS ID G213541E63316N, Alaska Range (a-c); snow avalanches and another land slide on Tlikakila Glacier Fork, Alaska Range (d-f); the detachment of Sedongpu Glacier, Eastern Tibet (g-i); and decreased backscatter in the firn area of Langjökull, Iceland (j-l). Shown are the 2018-2019 NDI images (a,d,g,j), summer 2017 Sentinel-2 "before" images (b,e,k), autumn (october-december) 2017 Sentinel-2 "before" image (h), summer 2019 Sentinel-2 "after" image (c,l), summer 2018 Sentinel-2 "after" image (f), and autumn 2018 Sentinel-2 "after" image (i). The Sentinel-2 images have a cloud mask in black.

Comparison with surface velocity measurements
We compare the detection of surge activity from changing in Sentinel-1 radar backscatter with glacier surface velocity measurements of two regions, Svalbard and Alaska. The goal of this comparison is to determine whether the surges we detected can be confirmed by corresponding changes in observed glacier surface velocity. In Alaska we control for false positives by comparing the velocity fields of glaciers we have detected surge activity from the NDI images only. For Svalbard, we also check the velocity fields of all glaciers in the region for surges other than the ones we identified from radar backscatter differences that we may have missed with our method, i.e. in Svalbard we control for both false positives and false negatives. Figure S5 shows own Sentinel-1 derived 12-day glacier velocity maps for Januaries 2018 and 2019, and December 2019 [see e.g. Strozzi et al., 2017, for details on the method; data can be downloaded from https://www.nve.no/hydrology/glaciers/copernicus-glacier-service/glacier-velocity/], as well as the velocity difference between January 2018 and 2019 and the locations of surge activities from our backscatter-difference method. The glacier velocity maps show a significant amount of noise in the measured 12-day glacier velocity in the accumulation area of the glaciers, especially in the January 2018 (Fig. S5a) and January 2019 (Fig. S5b) and therefore also in the velocity difference map (Fig.  S5d). The increasing/decreasing surge activity of Monacobreen, Sonklarbreen, Negribreen, Arnesenbreen, Strongbreen, Recherchebreen, Svalisbreen, Wahlenbergbreen, Penckbreen and Tunabreen are cleary reflected in the velocity changes, mostly even at the scale depicted here. For Pulkovbreen and Ganskijbreen, the change in velocities is hard to measure and recognize due to the small size of the glaciers. The optical images (not shown), though, clearly show a surge of Ganskijbreen, and a weaker or partial surge for Pulkovbreen. The latter two glaciers are examples where the backscatter-based method makes surge detection easier compared to velocity change detection. To conclude, there is a January 2018 to January 2019 velocity increase in a north-eastern outlet glacier of Austfonna which does not show a signal in the NDI images. This is the only glacier instability in the velocity maps that is not picked up by the NDI images, but it is an open question whether this acceleration is related to surge activity as in the optical imagery there are no signs of a surge.

Alaska glacier surface velocity
In Alaska we use annual glacier surface velocity provided by NASA MeaSUREs ITS_LIVE project [Gardner et al., 2020]. ITS_LIVE generated surface velocity from optical satellite images using auto-RIFT [Gardner et al., 2018], and we use the annual velocity for the years 2015, 2016 and 2017 with a resolution of 120 m. We have to exclude the velocity field for 2018, as currently the data for this year have limited quality [Gardner pers. comm.]. This is due to a low number of scenes processed for this year so far. Still we expect that looking at the surface velocity in the years 2015 -2017 will provide an indication of surges that we detect with radar images from the period winter 2018 -winter 2019. If a peak in surge activity is reached in the years 2015-2017, we should be able to detect a decrease in radar backscatter in 2018-2019. In addition, as shown in the example of Negribreen, Svalbard, increased radar backscatter can occur even after the peak in surface velocity.  surge. For Reid Glacier surface velocities are just slightly higher in 2017 than they are in the two proceeding years, such that there is no evidence for a glacier surge in the surface velocity data of this glacier.
For the glaciers we observe a decrease in the radar backscatter, the interpretation of the velocity fields, shown in Fig. S7, is less straightforward. For only two glaciers (Walsh Glacier and Fisher Glacier), we find a strong decrease in surface velocity over the 2015 -2017 period such that these glaciers can be interpreted as being at the end of a surge. Five glaciers (a tributary of Walsh Glacier, Turner Glacier, Hubbard Glacier, La Perouse Glacier, and Le Comte Glacier) exhibit an increase in glacier surface velocity in 2017. As we have no reliable data for the surface velocity in 2018, the year used for the radar backscatter analysis, we can not exclude that in this year the velocity and surface roughness dropped significantly again. That would be needed to explain a decrease in radar backscatter. In that case the surface roughness has to respond faster than we observed in our example of Negribreen. To conclude, Tkope Glacier has no marked change in glacier surface velocity.
For the 2 glaciers that showed both increased and decreased radar backscatter, the surface velocity is shown in Fig. S8. Margerie Glacier has much higher velocity in 2017 than it has in 2015 and 2016, which indicates a surge event like for the other glaciers we discuss above. For Klutlan Glacier, the ITS_LIVE data give a high velocity in 2015, after which the glacier flow is lower in 2016 and increases again in 2017. The velocity in 2017 is partly higher than in 2015 but more unevenly distributed over the glacier. We believe the lower surface velocity in 2016 could be due to a incorrect processing of the fast flow of this glacier, as Altena et al. [2019] find a continuous propagation of the surge bulge. The passage of the surge bulge can explain why we find both areas with increased radar backscatter as well as areas with decreased backscatter on the surface of Klutlan glacier.
ITS_LIVE velocity data also provide a formal uncertainty to the velocity fields. The uncertainty in the velocity for the glaciers we discuss here is shown in Figures S9-S11. Uncertainty in surface velocity is higher in the first two years (2015, 2016) than it is in 2017, and in general the uncertainty is much lower than the observed velocity for the glaciers of our interest, especially for the velocities in 2017, such that it is unlikely that the marked increases and decreases in surface velocity are due to measurement errors.
Based on the glacier surface velocity shown in Fig. S6 -S8, we can conclude that for 16 of the 18 glaciers the ITS_LIVE velocities support the detection of glacier surge activity from the normalized difference image of Sentinel-1 radar data, as the surface velocity has a significant change in the years 2015 to 2017. Two glaciers do not have a significant change in surface velocity during the years 2015 to 2017 that can support our surge activity observations from Sentinel-1. For both Reid Glacier and Tkope Glacier, the radar signal observed in 2018-2019 has disappeared in the 2019-2020 normalized difference image and the glacier extent seems nearly unchanged over the period 2017 -2019. Therefore, we have most likely misinterpreted an increase and a reduction in backscatter, respectively, for surge activity while in fact it was caused by another process, most likely a process in the firn pack. To conclude, we should mention that if we look at Fig S6, we see in the surface velocity plots for the glacier with GLIMS ID G216046E62056N there are two glaciers that show a similar significant increase in surface velocity: glacier G216046E62056N in the centre of the plots and Copper Glacier (GLIMS ID G216225E62044N) to the east of it. However, we do not see a signal in the 2018-2019 normalised difference image for Copper Glacier, while we do see an increase in radar backscatter for glacier G216046E62056N.