Spatio-Temporal Distribution of Supra-Glacial Ponds and Ice Cliffs on Verde Glacier, Chile

Known for their important role in locally enhancing surface melt, supraglacial ponds and ice cliffs are common features on debris-covered glaciers. We use high resolution satellite imagery to describe pond-cliff systems and surface velocity on Verde debris-covered glacier, Monte Tronador, and Southern Chile. Ponds and ice cliffs represent up to 0.4 and 2.7% of the glacier debris-covered area, respectively. Through the analyzed period and the available data, we found a seasonality in the number of detected ponds, with larger number of ponds at the beginning of the ablation season and less at the end of it. Using feature tracking, we determined glacier surface velocity, finding values up to 55 m/yr on the upper part of the debris-covered area, and decreasing almost to stagnation in the terminus. We found that larger ponds develop in glacier zones of low velocity, while zones of high velocity only contain smaller features. Meanwhile, ice cliffs appeared to be less controlled by surface velocity and gradient. Persistent ice cliffs were detected between 2009 and 2019 and backwasting up to 24 m/yr was measured, highlighting significant local glacier wastage.

Supraglacial ponds are common features on DCGs, where they have been observed for decades (Iwata et al., 1980;Kirkbride, 1993). Studies about supraglacial ponds on DCGs are mostly centered in High Mountain Asia (e.g., Sakai and Fujita, 2010;Miles et al., 2018;Watson et al., 2018a;Chand and Watanabe, 2019), where they can locally represent significant portion of the DCGs. Seasonal and interannual variations in ponds count and area have been observed (Steiner et al., 2019), as well as rapid filling and draining (Miles et al., 2017b), associated with englacial hydrological networks (Miles et al., 2017a;Miles et al., 2019;Watson et al., 2018b). Ponds also regulate DCGs runoff, as they act as buffer reservoirs (Irvine-Fynn et al., 2017). Exposed ice cliffs often appear at the ponds margins and can help identify the emplacement of drained features (Miles et al., 2017a).
Ponds and ice cliffs are known for their important role in locally enhancing surface melt (Sakai et al., 2000;Benn et al., 2012;Salerno et al., 2012;Maurer et al., 2016), as they absorb and transfer atmospheric energy into the glacier ice (Sakai et al., 2002;Steiner et al., 2015;Brun et al., 2018;Miles et al., 2018). This induces the DCGs to lose as much mass as debris-free glaciers despite the debris layer, which is known as the "debris-cover anomaly" Brun et al., 2018;Bisset et al., 2020). Based on energy balance modeling and satellite-based pond distribution, Miles et al. (2018) estimate that supraglacial ponds contribute to 12 ± 2% of the total ice mass loss in Langtang catchment, Nepal. Thompson et al. (2016) and Reid and Brock (2014) found that ice cliff covering only 5 and 1.3% of glacier surface contribute to 40 and 7.4% of the ablation on DCGs, respectively, in the Himalaya and Alps. Similarly, Brun et al. (2018) found that ice cliffs have a net ablation rate 3.1 ± 0.6 times higher than the average glacier tongue surface. Also, ponds are likely to coalesce and form large glacial lakes, enhancing glacier calving (Chikita et al., 1998;Röhl, 2008;Sakai and Fujita, 2010) and increasing glacial lake hazards (Quincey et al., 2007).
Mapping and inventorying supraglacial ponds and ice cliffs is therefore important due to their role as ablation hotspot and their link with the englacial hydrological network. Understanding of their magnitude and spatio-temporal variation will become increasingly significant on glaciers with negative mass balances (Deline, 2005;Stokes et al., 2007;Benn et al., 2012;Thakuri et al., 2014;Tielidze et al., 2020). In this context of glacier shrinking, mountain ranges are thought to show some transition from debris-free to debris-covered and rock glaciers, with an increase of their relative importance on the freshwater availability, quality and timing Knight et al., 2019).
The present study depicts the first assessment of supraglacial ponds and ice cliffs distribution on a DCG in Chile. Our objectives are to: 1) Document interannual changes in ponds and ice cliffs distribution using Google Earth historical imagery. 2) Analyse the inventory by determining horizontal velocity through feature tracking and extracting surface gradient from digital elevation model (DEM). 3) Detect persistent ice cliffs and measure the associate backwasting.

REGIONAL CONTEXT AND AREA OF INTEREST
The Southern Andes cover more than 4,500 km from northernmost Chile to the southern tip of South America in Tierra del Fuego. The extensive latitudinal range from subtropical to subantarctic domains (20°-55°S), steep elevation gradients, and the north-south orientation perpendicular to the prevalent atmospheric circulation cause significantly different climatic conditions, and consequently, a great variety of ice mass types along the Southern Andes (Lliboutry, 1998). North of 35 S, knows as the Dry Andes, encompasses a high (3,500-6,900 m a.s.l.) mountain range with arid (<500 mm/yr) condition to the north and modest annual precipitation amounts (1,000 mm/yr) to the south (Viale et al., 2019). Due to the topographic relief, geological context, and climatic setting, this region hosts the largest debriscovered glacier area of the Southern Andes (Barcaza et al., 2017;Zalazar et al., 2020). This region also contains the largest concentration of rock glaciers in the Andes (Barcaza et al., 2017), and includes many complex units that start as clean-ice glaciers at high elevations, gradually turn into debris-covered glaciers further down, and finally end as rock glaciers at the lowest sectors Kinnard, 2015, Monnier andKinnard, 2017;Zalazar et al., 2020). Recent inventories of the Argentina side of this region found that between 60 and 70% of the glacierised area (1,200 km 2 ; excluding rock glaciers) are partially (10%) to totally (˃90%) debris-covered .
South of the 35°S, in the Wet Andes, where the elevation of most peaks usually does not exceed 4,000 m a.s.l., the precipitation amounts increase considerably, exceeding 2,000 mm/yr (Viale et al., 2019). Although, at the southern part of this regions the topographic and climatological conditions allow the development of numerous and extensive glacierized areas, which constitute the largest glacierized surface in South America. In the north part, also known as the North Patagonian Andes (35°to 45°S), glaciers are smaller than those located further south. Although glaciers along this region are mainly clean ice or debris free (98% of the glaciated area), debriscovered glaciers can still be found due to local conditions such as rock-fall and stagnation. This is the case of Verde and other valley glaciers at Monte Tronador, where rock-falls and avalanches below a massive bedrock cliff present between 1,700 and 1,400 m a.s.l., allow the concentration of debris over the glaciers tongues (Ruiz et al., 2017).
The seasonal variation in the North Patagonian Andes is driven by the location and intensity of the southern hemisphere westerlies (Garreaud et al., 2009), which bring abundant precipitations between April and September (Aravena and Luckman, 2009). At these latitudes, orographic effect induces an increase of the annual mean precipitation from the Pacific coast to the western slopes of Chile, where it reaches more than 3,000 mm (Viale and Garreaud, 2015). Ruiz et al. (2015) measured more than 3,000 mm w.e. between May and September 2013 on a stake located close to the ice divide between Alerce and Castaño Overa glaciers, and close to the equilibrium-line altitude (ELA), which lies at 2000 m a.s.l. (Carrasco et al., 2005;Condom et al., 2007).
Verde glacier (8.15 km 2 , 41.21°S, 71.91°W, Table 1) lies on the southern flank of Monte Tronador (3,475 m a.s.l.), an extinct stratovolcano located in the North Patagonian Andes on the Chile-Argentina border ( Figure 1). The glacier is divided into three distinct parts: the accumulation zone on the upper slope of the Monte Tronador between 2,500 and 3,475 m a.s.l. an icefall between 1,400 and 2,500 m a.s.l. and a debris-covered tongue down to 980 m a.s.l. The debris-covered tongue is about 3.05 km 2 , which makes Verde glacier one of the most extended DCGs in Chile (Barcaza et al., 2017). Verde glacier has not shown significant retreat or advance over the last decades (Reinthaler et al., 2019), and its terminus is still in contact with the Little Ice Age moraines (Ruiz et al., 2015), which contrasts with the fast retreating behaviour of the glaciers in the region (Paul and Mölg, 2014). Between 2000 and 2012, the Verde glacier has a neutral mass balance (−0.08 ± 0.09 m w.e./yr), which contrast with the negative mass balance of Manso (−0.5 ± 0.1°m w.e./yr) and Casa Pangue (−0.29 ± 0.1 m w.e./yr), the other debris-covered valley glaciers at Monte Tronador ( Figure 1). Meanwhile, the most considerable ice thickness changes of the Manso and Casa Pangue glaciers are concentrated in their lower debris-covered tongues. The Verde glacier's lower part does not show significant elevation change between 2000 and 2012 (mean elevation change of −0.6 ± 0.5 m; Ruiz et al., 2017). By applying cross-correlation to Pléiades satellite images, Ruiz et al. (2015) measured surface displacements over Monte Tronador glaciers between March and June 2012. They found that these glaciers follow a radial flow pattern. At Verde glacier, maximum surface speeds of <390 m/yr were estimated on the steep icefall area. Meanwhile, the lower reaches of the debris-covered tongues of Verde and Casa Pangue glaciers are almost stagnant. They found that low-elevation debris-covered glacier tongues show increasing velocities at the beginning of the accumulation season, probably in response to an increase in water input to the subglacial system from winter rainfall events at low elevations and a decrease in meltwater production at higher elevations.
Recently, Zorzut et al. (2020) calculated the ice thickness distribution of Monte Tronador glaciers and found that the gently sloped debris-covered tongue of the Verde glacier is one of the thicker parts of the whole glaciers of Monte Tronador, with a maximum estimated thickness of around 180 ± 60 m and a total volume of 0.59 ± 0.2 km 3 .

Google Earth Imagery
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.

Ponds and Ice Cliffs Inventory
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 builtin 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.

Glacier Characteristics
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 highresolution 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 , 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.

Ponds and Ice Cliffs Inventory
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 2 ), which represents 0.44 ± 0.02% of the debris-covered area. Regarding the maximum ice cliff cover (82,364 ± 9,355 m 2 ), 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).

Surface Velocity
A total of 191 boulders were tracked at the glacier surface, allowing velocity measurements over most of the debriscovered 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. 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.

Ice Cliff Backwasting
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 FIGURE 2 | Tracking of a boulder located at the surface of the Verde glacier. Each dated square represent the same boulder on distinct cropped images from Google Earth. Red arrow on Figure 6A indicates the path of the present boulder. Coordinates in meters, universal transverse mercator (UTM) projection, zone 19, World Geodetic System (WGS84) datum.
Frontiers in Earth Science | www.frontiersin.org June 2021 | Volume 9 | Article 681071 from the north on Figure 7), while the measured retreat was about 9 m/yr.

Ponds and Ice Cliff Distribution and Formation
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    Ruiz et al. (2015) and the elevation along the same profile, respectively. 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 2 ) 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 lowerhalf 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.

Seasonality
No long-term trend in ponds and ice cliffs count or covered area was detected between 2009 and 2019. On the other hand, and . 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.

Characteristics of Backwasting Ice Cliffs
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 polefacing. 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 . 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.

CONCLUSION
We have compiled a complete inventory of supraglacial ponds and ice cliffs of the debris-covered tongue of Verde glacier, in the Chilean Northern Patagonian Andes. The inventory was based on Google Earth historical imagery covering the 2009-2019 period. Coverage ranges from 0.02 to 0.44% for ponds, and from 0.51 to 2.70% for ice cliffs.
Supraglacial ponds and ice cliffs are present at every elevation range of the debris-covered tongue. Ponds distribution is controlled by surface characteristics, with large features only occurring on areas of reduced velocity and low gradient. Ice cliffs distribution is less clear and appears to be slightly controlled by surface velocity. Debris thickness appears to also play a role in the presence or not of ponds and cliffs, as areas of thicker debris covers exhibit none of the two features. Based on available data, supraglacial ponds count and icecliff count and area exhibit a seasonal signal, with a maximum during the beginning of the ablation season and a minimum at the end of the ablation season and through the accumulation seasons. We attribute this seasonality to the opening and closure of the hydrological network through the year. Finally, we identified eleven persistent ice cliffs, most of them south-facing (or pole-facing). Ice cliff backwasting up to 200 m was observed, which highlights that ablation processes take place on the glacier surface, contrasting with the stagnant terminus.
Further researches need to be conducted to address the role of supraglacial ponds and ice cliffs on glacier mass balance, as well as their impact on the hydrological regime of DCGs. A continuous in-situ monitoring of ponds and ice cliff evolution could give insights into their interannual evolution, as well as the frequency and velocity of drainage and filling.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
TL gathered and prepared all data, performed all processing and calculations, and made all figures and tables. TL and LR contributed to the discussion of results, and shared the writing of the paper.