The Google Earth platform freely gives access to optical imagery of high spatial resolution, from SPOT and DigitalGlobe (e.g., QuickBird, Worldview-1 and 2, and IKONOS), and orthorectified based on the DEM from the Shuttle Radar Topography Mission–SRTM (Schmid et al., 2015). Moreover, Google Earth provides an historical catalogue with sufficient temporal resolution to enable investigation of interannual processes. Google Earth has thus been widely used as a main or supporting tool when performing cryosphere-related inventories, mainly of rock glaciers (Rangecroft et al., 2014; Schmid et al., 2015; Nagai et al., 2016; Charbonneau and Smith, 2018; Jones et al., 2018; Pandey, 2019), but also debris-free (Tielidze et al., 2020) and debris-covered glaciers (Alifu et al., 2016a,b), as well as glacial lakes (Wilson et al., 2018).

In the present study, the debris-covered tongue of Verde glacier is analysed with eight high resolution scenes (<0.5 m), whose acquisition dates range from November 2009 to October 2019 (Table 2). Although the image distribution is uneven over the year, the dataset contains three images from end of summer (March 2012, March 2016, and March 2018), one from end of winter (September 2015), as well as four from spring and start of summer (November 2009, December 2012, November 2013, and October 2019). The scenes present clear atmospheric conditions, mostly snow-free, and show limited areas of shadow. Error in horizontal positioning between images will be assessed by marking non-moving features, such as outcrops and trees, outside of the glacier area. No image processing was performed for interpretation.

The inventory was conducted following the method employed by Bhushan et al. (2018), who analysed high-resolution imagery from Google Earth to identify supraglacial ponds and ice cliffs for 10 DCGs in the Zanskar Basin of Western Himalaya. By-eye ponds and ice cliffs identification and manual digitisation were realised directly in the Google Earth platform through the built-in geographic information system. The ponds and ice cliffs polygons were then exported into the QGIS geographic information system software where their distribution and geometric properties were assessed. Regarding the ice cliffs, the inventory records the planimetric area. Error in area determination was estimated by multiplying feature perimeter by the spatial resolution (Casassa et al., 2014), here rounded at 0.5 m.

The inventory allowed the tracking of persistent ice cliffs over the study period to assess retreat, which is defined here, as the observed distance between the same ice cliff between the initial and the final satellite scenes, divided by the time interval. In order to remove the displacement due to glacier flow and to only assess backwasting due to ablation, the observed ice cliff retreat was corrected for the local mean surface velocity (Steiner et al., 2019). Discriminating between the observed retreat and the assessed backwasting is particularly important for ice cliff retreat occurring parallel to flow direction, in which case the difference is maximum.

Horizontal surface velocity was determined following the method described by Immerzeel et al. (2014) and Wigmore and Mark (2017), who manually tracked topographic features on high-resolution photographs acquired using unmanned aerial vehicles. In the present study, this approach was adapted by tracking evenly spaced and clearly distinguishable boulders, at the glacier surface using the Google Earth historical imagery (Figure 2). For each identified boulder, a vector was drawn using the Google Earth path tool, which each node corresponded to the boulder position at the respective imagery date. These vectors were exported into QGIS to calculate horizontal displacement for each mapped boulder, which was then converted to surface velocity using time delay between images. Finally, velocity was interpolated to the tongue surface, using spline and kriging interpolation methods, both techniques having been applied on feature tracking on DCGs (Immerzeel et al., 2014; Wigmore and Mark, 2017). Interpolations were performed using the respective QGIS built-in tools. Error in distance calculation was rounded at the pixel size (0.5 m). Error due to image displacement was assessed by tracking non-moving features, such as trees and outcrops, outside of the glacier area.

Surface gradient was extracted from the SRTM DEM. However, this is made difficult due to the topography of DCGs being generally very coarse (Nicholson and Benn, 2012), and the SRTM showing random speckle noises (Stevenson et al., 2010), also referred as coherent multiplicative noises, that inherently exist in these products. Into QGIS, a low pass filter was thus initially applied, which is a common method to smooth SRTM DEMs (Wendi et al., 2016). Surface gradients were then grouped into four classes (0–2, 2–6, 6–10, and >10), as proposed by Reynolds (2000), and applied by Quincey et al. (2007) and Miles et al. (2017b) to analyse ponds and ice cliffs distribution. Root Mean Square Error (RMSE) was applied to compute the uncertainty of the two interpolation methods in comparison with the observations.

Over the 2009–2019 period, a total of 151 ponds and 940 ice cliffs were identified (Table 2 and Figure 3). An example of mapped pond and ice cliff are shown on Figure 4. Both the ponds and the ice cliffs cover most of the tongue area, with highest density in the middle part. It is important to note that most of the ice cliffs located in the upper part of the tongue (Figure 3B) are associated with supraglacial channels.

The most extensive pond coverage was observed on March 30, 2012 (13,517 ± 649 m²), which represents 0.44 ± 0.02% of the debris-covered area. Regarding the maximum ice cliff cover (82,364 ± 9,355 m²), it occurs on November 25, 2009 and represents 2.70 ± 0.30% of the glacier tongue. Most of the ice cliffs (866 of 940) were not in contact with ponds at the moment of their mapping. High interannual variability is observed in both ponds and ice cliffs coverage (Figure 5).

A total of 191 boulders were tracked at the glacier surface, allowing velocity measurements over most of the debris-covered tongue (Figure 6A). Some boulders show displacements not concordant with the flow line direction between two dates, which is likely due to topography effects and roll overs. Ten control points were tracked around the glacier tongue (Figure 6A), allowing to determine a maximum total displacement error of 0.3 m/yr along the flowline.

Both spline and kriging interpolations methods (Figures 6B,C) show statistically similar results (Table 3). Highest velocities are found in the upper part of the glacier tongue, just below the steep icefall, with a regular decrease towards the glacier terminus, which is almost stagnant. This is consistent with the surface slope becoming gradually shallower towards the snout (Figure 6D). Low velocities are also found on the margins of the glacier. Our results are concordant with findings from Ruiz et al. (2015) which shown that higher velocities are associated with the ice fall and found almost stagnant conditions in the lower reach and lateral margins of the debris-covered tongue (Figure 6D). Otherwise, some discrepancies appear in the upper limits of the covered-tongue, likely associated with a lack of mapped boulders in this area of the glacier and a greatest uncertainty while applying interpolation techniques.

Eleven ice cliffs were identified to persist over at least two scenes of the period of interest. Retreat varies from a few meters to a maximum distance of 201 m between November 2009 and October 2, 2019 for a ∼200 m wide ice cliff (Figure 7B). After correction for local surface velocity, the highest backwasting rate was assessed to 24 m/yr for a south-facing ice cliff (third arrow from the north on Figure 7), while the measured retreat was about 9 m/yr.

Although supraglacial ponds and ice cliffs are observed across the whole elevation range, their distribution is uneven over the tongue (Figure 3). Four main clusters can be identified, from North to South, in particular when describing the ice cliffs distribution (Figure 3B). The first cluster counts few ponds and two channels of numerous ice cliffs parallel to the flowline. The formation of these ice cliffs seems to be associated with supraglacial channels visible on the upper part of the debris-covered tongue. The second and the third clusters show the largest populations of ponds and ice cliffs. A clear gap, showing no features, is visible between the first and the second clusters, which we attributed to a huge quantity of debris deposited on the tongue by a landslide or rock avalanche originated on the western flanks of the valley (Figure 1). A smallest gap is observed between the second and the third cluster, likely associated with a former landslide or rock avalanche, forming visible arched ridges on the glacier surface. These thick deposits of debris appear to be inhibiting the formation of supraglacial ponds and ice cliffs at the glacier surface. Finally, the fourth cluster, located near the front, shows smaller density of ponds and ice cliffs, whereas reduced surface velocity and gentle slopes would be adequate for their formation (see Figure 3). We attributed this small amount of feature to a significant debris thickness, highlighted by the presence of trees and vegetation on the glacier surface (Figure 1). The lowest ponds and ice cliffs numbers in the areas with thicker debris layers is concordant with results from Steiner et al. (2019) in the Langtang catchment (Himalaya), who observed few features near the front, where the debris layer is the thickest. The absence of ponds and cliffs in areas of thick debris is likely due to suppression of ablation, as the insulation threshold might be exceeded.

All mapped ponds take place on the tongue with gradient below 10 (Figure 8). Above this threshold, all the available meltwater is able to drain away, as previously reported by (Reynolds, 2000) on DCGs in Bhutan. Specifically, 14, 74 and 12% of the ponds were observed on slopes <2°, between 2–6 and between 6–10, respectively. These distributions are concordant with observations made on DCGs in Himalaya (Liu et al., 2015; Chand and Watanabe, 2019).

The ice cliffs distribution shows a similar pattern with less than 1% observed on areas with gradient above 10, while 10, 73 and 16% of the mapped cliffs are located on slopes <2°, between 2–6 and between 6–10, respectively. These similar distributions of ponds and ice cliffs with respect to gradient could signify their formations are associated processes. Indeed, despite the fact that only 70 of 907 ice cliffs were in contact with a pond at the moment of their mapping, we observed that most of the backwasting ice cliffs (Table 4) are in contact with a pond at the initial stages of their formation (Figures 7B,C). Afterwards, the cliffs lose contact as they retreat and the ponds are either drained or filled with debris.

The distribution of ponds with respect to their surface area is strongly controlled by both surface gradient and surface velocity (Figure 8). Large ponds (>1,000 m²) only develop on glacier surface of relatively low velocity (<17 m/yr) and low gradient (<6). This results in large ponds being constrained to the lower-half and the margins of the tongue (Figure 3A).

Regarding the distribution of ice cliffs with respect to their area, no clear influence of surface velocity nor surface gradient was observed (Figure 9). Nevertheless, it appears that large ice cliffs only develop where surface velocity is low (Figure 9A).

The observed high variability in ponds and ice cliff distribution and their interannual evolution seems to highlight a very dynamic behaviour of Verde glacier, which is concordant with observed high surface velocity. These observations tends to demonstrate that the supraglacial ponds do connect and disconnect with the englacial conduits as crevasses open and close, respectively.

No long-term trend in ponds and ice cliffs count or covered area was detected between 2009 and 2019. On the other hand, and based on the available data, we observed a seasonal signal in the number of supraglacial ponds and in the number and area of ice cliff mapped through the year (Table 2; Figure 5). We found a maximum of supraglacial ponds and ice cliffs during spring and at the beginning of summer (November 2009, December 2012, November 2013, October 2019), and a minimum of ponds and ice cliffs, at the end of summer (March 2012, March 2016, March 2018), and at the end of winter (September 2015). We suggest the number of ponds and number and area of ice-cliff might be related to the evolution of the hydrological network through the year, as previously proposed in other studies (Reynolds, 2000; Sakai et al., 2000). At the beginning of the ablation period (spring and beginning of summer), the hydrological network is not well connected yet, so meltwater accumulated into depressions and form supgraglacial ponds. This early pond formation initiates the exposure of associated ice cliffs. During the ablation season, the crevasses and englacial conduits open and enhance connectivity. Similar cycles have been observed in the Eastern Himalaya and Karakoram, with an increase in ponds number at the beginning of spring (April) and a sharp decrease in June-July (Narama et al., 2017). Authors attributed the increase to inflow of meltwater from snow and ice, and the later decrease to a enhanced connectivity to the englacial drainage network. The earlier part of the cycle was also observed in High Mountain Asia during the pre-monsoon season, and related to the observed reduction of the snow cover (Miles et al., 2017b; Chand and Watanabe, 2019).

Finally, during the beginning of the accumulation season (winter), there is no meltwater at the surface to fill the ponds. This small amount of ponds and their reduced cover in winter are concordant with results from the Everest region (Chand and Watanabe, 2019).

The smallest amount of ice cliffs at the end of summer might also be related to their fast-vanishing behaviour, as few of them have been observed to persist year to year which could also explains that no seasonality was observed regarding the ponds area.

Most of the persistent and backwasting ice cliffs were found to be south-facing (Table 4), which is consistent with results from Sakai et al. (2002) and Buri and Pellicciotti (2018), who found that persistent ice cliffs are mostly north-facing in the northern hemisphere. In all cases, persistent ice cliffs appear to be pole-facing. In contrast, the sun-facing ice cliffs do not have time to extent spatially, and do not survive the ablation season, as they receive the highest amount of radiations. They are thought to not contribute significantly to annual mass balance (Buri and Pellicciotti, 2018). According to Sakai et al. (2002), the difference in stability between the pole-facing and the sun-facing ice cliffs is caused by a difference in incoming radiations. Pole-facing cliffs receive mostly long-wave radiation from the surrounding debris cover, which reaches the lower portion more than the upper portion. This results in the ice cliff being steep enough to not be covered by debris and persist over time. On the other hand, sun-facing cliffs are exposed to direct and more intense short-wave radiation, which affect mainly their upper portion. These cliffs have thus shallower slope and are more likely to get covered by debris.

The growth and demise of a persistent south-facing ice cliff on glacier Verde is illustrated on Figure 7C. This south-facing ice cliff initially emerges as a concave feature in contact with a pond. The ice cliff then grows spatially and reach a maximal extent, before it tails off as a convex feature. In the later stages, there is no apparent pond adjacent to the persistent ice cliffs.