DISTRIBUTION AND EVOLUTION OF ICE APRONS IN A CHANGING CLIMATE IN THE MONT-BLANC MASSIF (WESTERN EUROPEAN ALPS)

Ice Apron (IA) is a poorly studied ice feature, commonly existing in all the world’s major mountain regions. This study aims to map the locations of the IAs in the Mont Blanc massif (MBM), making use of the very high-resolution optical satellite images from 2001, 2012 and 2019. 423 IAs were identified and accurately delineated in the MBM on the images from 2019, and their topographic characteristics were studied. We generated our own Digital Elevation Model (DEM) at 4 m resolution since the freely available products predominantly suffer from significant inconsistencies, especially in steep mountain areas. Results show that most IAs exist at elevations above the regional Equilibrium Line Altitude (ELA), on steep slopes, on concave surfaces, on northern and southern aspects and on the most rugged terrains. They are also commonly associated with steep slope glaciers as 85% of them occur on these glaciers’ headwalls. A comparison between 2001 and 2019 shows that IAs have lost around 29% of their area over a period of 18 years. This is significant and the rate of area loss is very alarming in comparison with the larger glacier bodies. We also studied the effect of topographic parameters on the area loss. We found that topographic factors like slope, aspect, curvature, elevation and Terrain Ruggedness Index (TRI) strongly influence the rate of area loss of IAs.


INTRODUCTION
The impacts of climate change on mountain environments are of great concern to the scientific community. Like other mountain terrains in the world, the Mont-Blanc massif (MBM) has experienced a general trend of glacier retreat since the end of the Little Ice Age (LIA) despite small re-advances culminating in 1890, the 1920s and 1980s (Bauder et al., 2007). In the MBM, the loss of glacier surface area as recorded from 1985 to 2007 is more than 7% (109.5 km 2 in 1985 to 102.4 km 2 in 2007) in relation to the ELA (Rabatel et al., 2013). Similarly, the recorded loss of glacier surface area was 24% of the total area from the end of the LIA to 2008 (Gardent et al., 2014). The reported loss of ice thickness is also noteworthy. For example, the loss of ice thickness at the front of the Mer de Glace glacier (1400 m a.s.l) from 1979 to 2003 was up to 60 m (Berthier, 2007), while that of Argentière glacier (1500 m a.s.l.) was 80 m from 1994 to 2013. Similarly, the surface of Talèfre glacier lowered by 5-10 m between an altitude of 3000 -3500 m a.s.l from 1979 to 2003 (Berthier, 2007). The surface of Géant glacier also lowered by 20 m between 1992 and 2012 at 3550 m a.s.l. (Ravanel et al., 2013). The glacial retreat and shrinkage also concur with the rise of the Equilibrium Line Altitude (ELA) by 170 m between 1984 and 2010 in the western Alps (Rabatel et al., 2013). However, all the above examples focus on primarily large glacial systems like cirque and valley glaciers. Researchers have overlooked small ice features like Ice Aprons (IAs) until now; hence very few studies have been dedicated to their existence and consequent disappearance. One of the primary reasons for the disregard towards IAs could be attributed to their small size; hence their disappearance, as an effect of global warming, does not drastically affect the water resource availability. However, the presence of IAs is considered as a necessary condition for * Corresponding author climbing many mountaineering routes (Mourey et al., 2019). Carbon 14 dating analysis of an ice core extracted from an IA on the Triangle du Tacul showed the age of ice to be around 2640 ± 130 years (Guillet et al., in review). Considering the old age of the ice preserved in the IAs, their shrinkage and eventual disappearance in the next decades should also be regarded as a loss of an important glacial heritage. Guillet and Ravanel, 2020 defined IAs as very small (typically smaller than 0.1 km 2 in extent) ice body of irregular outline, lying on slopes >40 • , regardless of whether they are thick enough to deform under their own weight. The same study also showed that IAs have been losing area since the end of the Little Ice Age (LIA ≈ 1855) as a result of the changing climate (Guillet and Ravanel, 2020). However, the impact of topographic factors along with climatic variables on the area loss of IAs has not yet been understood. Since, the topographic factors for a bedrock are non-dynamic and remain constant over time, their impacts can be modelled accurately. Keeping this into perspective, this article primarily aims to define the spatial distribution and characteristics of IAs in the MBM in order to understand their evolution. The second aim is to assess the impacts of global warming and topographic factors like the slope, aspect, elevation, curvature, topographic ruggedness index (TRI) on the area lost by IAs over the past two decades.

Study Area
The MBM is located in the north-western Alps between France, Switzerland, and Italy. It has an area of about 550 km 2 , of which the glaciers cover about 145 km 2 (about 26 %) ( Figure. 1). It is the highest and most glacierized massif in the French Alps (Gardent et al., 2014). The highest point in the MBM is the summit of the Mont Blanc at 4809 m a.s.l. There are also around a dozen peaks with elevations ≥ 4000 m a.s.l. Because of the considerable high elevation of the MBM, there are about 100 glaciers, including 12 large glaciers ( ≥ 5 km 2 ) bordered by steep sidewalls. The largest glacier in the MBM is the Mer de Glace which is 12 km long, with an overall surface area of 30 km 2 . The MBM shows an asymmetry across the range in its geomorphometry. Out of twelve, six of its largest glaciers are located on its North-Western (NW) aspect where slopes are gentle, and glaciers are well fed by the westerly winds and melting is reduced by the protection of the shaded northern aspect. As compared to this, the South-Eastern (SE) aspect is characterized by smaller glaciers and, in general, by steeper slopes bounded by very high sub-vertical rock walls. This asymmetry is also evident by the difference in climatic conditions observed on the two sides of the MBM. For example, the Mean Annual Air Temperature (MAAT) recorded in Chamonix at an altitude of 1,044 m a.s.l is +7.2 • C, while that in Courmayeur (1,223 m a.s.l) is + 10.4 • C.

DATASETS
The datasets used in the study are summarized in table 1. The following sections provide a detailed description of the each type of data.

Data
Resolution (m) Acquisition

Digital Elevation Model
Digital Elevation Models (DEM) files are provided in several formats, but the elevation data is usually grouped in regular grids (rasters) or triangulated irregular networks (TIN). With the advent of remote sensing techniques, DEMs are now predominantly generated using remote sensing techniques. Some of the standard remote sensing techniques are photogrammetry (Coveney and Roberts, 2017), spaceborne and airborne Interferometric Synthetic Aperture Radar (InSAR) and Light Detection And Ranging (LiDAR). The benefits of the techniques mentioned above are that a large spatial area can be mapped in considerably reduced time and fewer people. Also, the cost of generating DEMs through these methods is comparatively much lower than other techniques. An important point to note before promptly using a DEM file for any application is that the final product delivered to the user is usually a result of heavy modelling and pre-processing steps typically comprising many errors (systematic and random) uncertainties. The accuracy of a DEM is measured as vertical and horizontal accuracy. Vertical accuracy is simply the difference between modelled height and the actual elevation on the ground calculated from a standard reference. Horizontal accuracy relates to the positional shift or bias in the X and Y directions between features (Ground Control Points) that can be accurately identified in the DEM and on a map. Depending on the application, even minor elevation errors can considerably affect the accuracy and usefulness of the results (Holmes et al., 2000). Many topographic products like slope, aspect, curvature etc, that form the basis for advanced topographical studies, also suffer from inconsistencies of the original DEM products (Januchowski et al., 2010). Assessment of the nature and extent of these errors is generally very complex and tedious and not readily available to spatial data users. (Wechsler, 2007). This leads to many uncertainties in the elevation data that is hard to comprehend. Commonly, when a DEM product is produced locally by a user, information about its accuracy and precision is well documented and understood. This is often not the case with global products, which are pre-processed, and any information about errors introduced during processing is rarely readily available. Considering the theme of our work and the small dimensions of the objects we are working on, it was of paramount importance to have a high spatial resolution DEM with reasonable accuracy to avoid errors due to DEM inaccuracies. Keeping this in mind, we decided to put additional efforts into building a new DEM instead of relying on a global product that is freely available.
As part of the KALIDEOS -Alpes project, we acquired multiannual stereoscopic sub-meter optical images from the Pleiades constellation over the MBM. For each acquisition, a 4 m DEM was computed from a pair of stereo images of 0.5 m resolution from 2019 (25/08/2019) using the Ames Stereo Pipeline (ASP) (Shean et al., 2016). ASP uses rational polynomial coefficients (RPCs) provided with the Pleaides images, eliminating the requirement of a large number of high accuracy GCPs. We used the same parameters as (Marti et al., 2016) for our processing. The point cloud derived from ASP stereo routines was gridded into a single 4m DEM (point2dem function into ASP). Then ortho-images (12-bit dynamic range) were generated for both stereo acquisitions using the 4 m DEM generated by ASP with the DEM gaps filled with ALOS world 3D DEM available freely to ensure sharp and gap-free images. The robustness of the Pleiades data processing based only on RPCs has been tested previously by  and (Berthier et al., 2014), and it is considered sufficient to co-register the images and DEMs relatively using automatic co-registration methods.
In this sense, we applied an automatic co-registration method given (Nuth and Kääb, 2011), using a LiDAR 2 m DEM as a (reference DEM), to co-register the Pleiades acquisition DEM (source DEM) ( Figure. 2a.). We used the RGI (Randolf Glacier Inventory) glacier contours, forest extensions and manually delineated polygons to mask non-stable areas with an opensource script (https://github.com/dshean/demcoreg) (Shean et al., 2016) (Figure. 2b.). Then, the ortho-images were shifted (translation-only) using corresponding co-registration values (shift values for x, y and z in meters) detailed in figure. 2f. Figure. 2c. shows the difference in elevation between the reference DEM and the source DEM before co-registration. Figure. 2d. shows the results after the co-registration process. A strict visual inspection was also applied to improve the accuracy of the matching.
(a) (f) Figure 2. High resolution DEM processing chain for Pleiades stereo images Figure. 2e. shows the elevation difference (error in coregistration) considering only the stable areas. The distribution of errors can be visualized by a histogram of the sampled errors, where the number of errors (frequency) within certain predefined intervals is plotted (Höhle and Höhle, 2009). Figure. 2f. shows the histogram of the errors ∆h (elevation difference between the reference and source DEM) in meters for the stable areas. The accuracy estimate before and after the co-registration is shown by the normalized median absolute deviation (nmad) and the median value calculated together. The co-registration process was stopped when an acceptable value for nmad and the median was achieved. At the end of the processing we gen-erated a precisely co-registered high resolution DEM which was used to compute various topographic indices like the slope, aspect, curvature, elevation and Topographic Ruggedness Index (TRI).

High-resolution typological map of the Mont Blanc massif
Geomorphological or typological maps are among the most valuable tools for understanding the Earth's physical nature since they provide a detailed description of different landforms (Dramis et al., 2011). These maps include accurate information about the spatial properties of landforms like the shape, size, orientation, slope, relief, etc. In many studies, landforms are distinctly categorized into different types based on their origin, the process of their deposition or formation (morphogenesis), effects of the bedrock lithology/ structural control, the rate of dynamic or genetic processes (morphodynamics), the relative or absolute age (geochronology) and the type of bedrock and near-surface deposits. Many global products or inventories that claim to demarcate the glacier boundaries accurately are available to download. The most well-known of these are the GLIMS (Global Land Ice Measurements from Space) and the RGI. Although the glacier outlines from global inventories are reasonably accurate for global analysis and big glaciers, their accuracy on a local scale and smaller glacier bodies is questionable. These inventories also do not provide a sub-classification of the different types of glaciers. Considering the above factors and looking at the importance of having a highly accurate glacier inventory, we decided to prepare our database using very high-resolution satellite and aerial photographs available for the MBM. Unlike a general glacier outline, we divided the glaciers into subclasses, which are familiar to geomorphologists and glaciologists. at 20 cms resolution and the Pleiades stereo images at 50 cms resolution for areas not covered by the orthophotos. The final glacier inventory is the only inventory available at such a high resolution and with a sub-classification of glacier types and other features. For this study only the glaciers on steep slopes were considered. We mapped and divided the study area into the following classes: cirque glaciers, slope glaciers, slope glaciers CGS (connected to other glacial systems), slope glaciers HG (with hanging front) ( Figure.3).

Optical images
Since the work relies predominantly on the utilization of very high resolution optical images, the most important aspect of this study was to collect images from various sources to cover the entire study period. The datasets mentioned in table 1 according to their dates of acquisition were used for accurately delineating the boundaries of IAs for the years 2001, 2012 and 2019. For an accurate comparison between images it is imperative that all the optical images are precisely co-registered. Since our main data source for 2012 and 2019 comes from the Pleiades constellation, the two sets of images were already co-registered and no further processing was required for a comparison between these two sets of images. There was however a need to co-register the orthophotos of 2001 with the Pleiades images of 2012 and 2019. For this, we performed an automatic image to the image registration process in ENVI 5.0 (Exelis Visual Information Solutions). The process includes locating and matching several feature points called tie points in the base image and a warped image selected for registration. Here, we used the Pleiades image of 2019 as a base image and the orthophotos from 2001 as a warped image. We used the general cross-correlation algorithm, which works well for registering images with similar modality (e.g., registering an optical image with another optical image). The normalized cross-correlation between the patch in the base image and the patch in the warped image is computed as the matching score. The minimum matching score threshold we set was 0.7, and the tie points with a matching score less than this value were considered outliers and removed. Since we are working with a set of orthophotographs, the transformation model between the base image and the warp image fits an RST transform (rotation, scaling and translation) model. Maximum Allowable Error Per Tie Point was set at 10 pixels initially to get a dense network of tie-points. The search window size was set at 255, while the matching window size was set at 61. Our very high resolution image pairs resulted in a tie point network of 38 points with an overall average RMSE of 5 m. Fine co-registration procedures were performed by deleting tie-points with very high RMSE and manually selecting GCPs by identifying features on the base image and then marking accurately the same feature on the warped image. Features like road crossings, permanent settlements, bridges, train stations, peaks of high mountains and refuge huts in the study region make for good GCP points (Förstner and Wrobel, 2016). A collection of 21 manual GCP points spread throughout the study area were added, and a final accuracy stated in terms of the RMSE was achieved to be 0.22 m. This is less than half the pixel resolution of the images in question, and hence the results of the co-registration can be considered to be highly precise (Han and Oh, 2018). The misalignments or co-registration accuracies suffer at the outer edges of the image and in larger areas without GCPs. We verified the misalignment at the image borders to be no more than 2 pixels (or 1 m).

Location and morphometric characteristics of the Ice Aprons in the Mont Blanc massif
A total of 423 ice aprons were identified in the MBM by visual interpretation and manually digitized in a GIS environment.
Since the size of these glacier bodies is extremely small, we used a combination of very high spatial resolution images like the orthophotos available from Geoportail IGN (0.2 m resolution), panchromatic images from SPOT 6 (1.5 m resolution), Pleiades stereo images (0.5 m resolution) and high-resolution images available from Google Earth. Problems associated with the lack of coverage, cloud cover, illumination, and shadow that make visual interpretation difficult and the presence of seasonal snow cover were overcome by not relying on a single dataset. The seasonal variation of snow cover can lead to an overestimation of the IA area. Hence, it is crucial to select datasets at the end of the summer period (preferably late August or early September). The primary dataset we relied on was the Pleiades images from late August 2019, as these were the latest highresolution dataset available, plus the only dataset we had from the end of the ablation period. The total area occupied by IAs in 2019 in the MBM is 4.214 km 2 , which is 0.76 % of the area of the massif. Analysis of the location of IAs with mean altitude shows that only 49 IAs out of 423 (12%) exist below the regional Equilibrium Line Altitude (ELA) ( Figure. 4a.). We fixed the ELA line very conservatively at 3100 m a.s.l based on the calculations of (Rabatel et al., 2013), where they suggested the regional ELA to be located at 3035 ± 170 m for the western Alps. The heatwaves of 2015 and 2017 would have likely resulted in the ELA rising higher than the 3100 m threshold we have considered, but we did not account for this as an accurate measurement for this change is not available. With no recharge source, the IAs below the ELA are incredibly susceptible to melting due to global warming. Also, the rise in ELA will result in an increasing number of IAs becoming susceptible to melting in the future. However, most of the IAs show mean altitude values between 3300 -3500 m a.s.l ( Figure. 4a.). An analysis with the aspect shows that most IAs exist on mean northern aspects (N, NW, NE), while the eastern and western aspects are the least favourable ( Figure. 4b.). There is also a considerable number of IAs existing on mean southern aspects. The trend is similar to the glaciers in MBM, with the most significant glaciers found on the NW aspect, while smaller glaciers are found on the SE aspect. Another crucial factor to characterize the location of IAs is the slope. Since we have already restricted our criteria for selecting IAs to steep slopes (>40 • ) We observed that most IAs in MBM exist on mean slopes between 53 − 65 • (Figure. 4c.). In comparison, very steep slopes (>70 • ) and more minor steep slopes (between 40 − 50 • ) are less favoured. Also, a majority of IAs lie on concave surfaces, as compared to convex surfaces. (Figure. 4d.) To understand the topographic locations of IAs in the MBM, we made a final comparison with the Topographic Ruggedness Index (TRI). TRI was calculated by the methodology given by (Sappington et al., 2007). We divided our study region into four subclasses of TRI based on the degree of ruggedness. Figure. 4e. shows that most IAs exist in TRI classes of medium and high ruggedness compared to other classes.

Evolution of Ice Aprons
Using the inventory of glaciers we prepared from very high resolution orthophotos from 2015, and the locations of IAs identified previously, we tried to understand the evolution and relation of the IAs with respect to the other glacier systems. We  observed that around 85 % of the total IAs mapped in 2019 exist above a steep slope glacier or on the head walls of a cirque glacier. These IAs were part of large glacial systems in the past, but are separated today to form a new glacial entity by the presence of a bergschrund. Bergschrund is defined commonly as a deep crack near or at the top of a glacier, separating moving ice (glacier below) from ice that is not moving (IA above) ( Figure. 5 a.). There are many possible reasons for the development of a bergschrund, the primary of which is the difference in flow velocities, both tangential to the surface and vertical from the upper part (which has lower flow velocities) to the lower part (higher flow velocities). One of the contributing factors for the disparity in the flow velocities leading to the formation of a bergschrund is the increased melting due to global warming effects. Thus, the upper part of the glacier is separated from the main glacier body and thus forms an IA. Hence, we see an increasing number of IAs of this type in every mountain environments of the world. Figure 5a. explains the relation of IAs with respect to the steep slope glaciers and Figure 5b. shows few examples from the MBM.

Effects of topographic factors on relative area loss of Ice Aprons
Similar to the methodology applied to identify IAs in images from 2019, IAs were also identified and accurately delineated for the years 2001 and 2012 using a combination of very highresolution optical images mentioned in the smaller size of IAs makes them more affected by the global warming than large glaciers. Also, topography and local climate conditions effects are far more pronounced for IAs than for other glaciers. Results presented in figure 6 show that IAs present on western aspects (west and south-west), low altitudes, convex curvatures, high TRI values, and steep slopes (except for slopes >70 • ) values have lost more area since 2001 as compared to other factors. Hence, IAs present in these topographic conditions seem to be most impacted by global warming. The results are consistent for the analysis of the results from 2012 and 2019, as can be seen by the different colour bars in figure 6. In comparison, IAs located on north facing slopes, at high altitudes, low slope angles (except slopes>70 • ), concave curvatures, and low TRI values look more protected. The percentage of the relative area lost by IAs in these topographic conditions is the least. This shows the significant influence the topography has along with the climatic variables on the area loss of IAs in the MBM.

CONCLUSIONS
The first part of the analysis we performed for this study dealt with accurately identifying the locations and the area occupied by the IAs in 2019. We managed to accurately delineate 423 IAs in the MBM using a combination of very high-resolution optical satellite images. The total area of IAs in the MBM in 2019 is 4.21 km 2 . Most of the IAs occur on the N and S aspects, at elevations between 3300 and 3500 m a.s.l, on concave plan curvature terrains, on steep slopes between 53 and 65 degrees and generally medium to high ruggedness defined regions. This, in short, explains that IAs generally exist in the most complex topographic settings in the mountain regions. The relation of IAs with steep slope glaciers was also well documented by our research. An overwhelming majority of the IAs occur in relation to an existing steep slope glacier and above a bergschrund. This indicates a common link in their evolution genesis or history. Although, from this study, it is hard to pass a generalized comment of the evolution of IAs, one reasonable possibility was discussed along with the formation of a bergschrund. Further research into this hypothesis needs to be carried out to ascertain the development of IAs. However, this study provides initial evidence, which can lead to a future research direction on this topic. On comparing the area of IAs from 2001, 2012 and 2019, it was evident that the effects of climate change on the area loss of IAs is drastic. The rate of area loss of IAs is significantly higher compared to more prominent geomorphic features like the cirque and valley glaciers. It is also evident that the effects of the topographic parameters like slope, aspect, curvature, elevation, and TRI are more pronounced in the IAs. The small size of IAs make them more susceptible to changes. From 2001 to 2019, IAs have lost around 29% of their area in a time span of 18 years. This is alarming because, at this rate, we will lose most of the IAs that existed in 2001 before the end of this century. As previously stated, the loss of ice volume from an IA and their eventual disappearance by the next decades can also be regarded as a loss of an important glacial heritage. In the future, more interest and research focus needs to be diverted towards these small but critical glaciological features.