**1. Introduction**

Mountain glaciers are visible indicators of climate change, providing some of the clearest evidence of global warming, perturbation in the atmospheric flow pattern and consequent precipitation regime. Due to their relatively small size and high mass turnover rates, these glacial bodies demonstrate extremely sensitive to temperature and precipitation modifications [1,2]. For the European continent, Alpine glaciers are a key component of local and regional hydrogeological cycles. The globally observed retreat of these ice masses over the last century has significant impacts on the surrounding environment, affecting both physical and socioeconomic systems [3–5]. The loss of freshwater storage and its consequent release into seas and oceans, modifications in resource availability for water consumption, irrigation and power generation, glacier-related natural hazards (e.g., landslides, avalanches and outburst of glacial lakes, causing floods and debris flows) are between the potential consequences of the general shrinkage trend.

Variations in glacier volume and mass are therefore primary targets of investigation for the understanding of ongoing modifications and the forecast of possible future scenarios. Remote sensing techniques (e.g., geodetic surveys, satellite images, GNSS or LiDAR surveys) can help to monitor these fluctuations [6–9] from accurate time-lapse reconstructions of the topographic surface of the glacier. Comprehensive knowledge of the ice bottom depth and morphology is however needed to provide a total volume estimate and reliable mass balance evaluations [10–12]. Ice thickness characterization and monitoring are also of primary importance for the modeling of future glacier dynamics [13], hydrological projections [14], glacier-related natural hazards [15] and ice core analyses [16].

Traditional glaciological measurements, i.e., local probe deployment for snow and firn/ice thickness variations or density measurements, are usually limited to the shallowest part of the glacier and do not extensively cover wide areas of investigations. Despite their relevance, ice bottom depth and morphology are often poorly known, mainly due to inadequate characterization methods and logistical issues. Geophysical methods can overcome these limitations, mitigating the operational efforts and enlarging the depth of investigation and data density of traditional techniques.

In the last decades, the ground-penetrating radar (GPR) technique has assumed an increasingly important role in the glacial exploration for the imaging of subsurface conditions [17]. The non-magnetic and low-conductive snow/ice column is indeed generally favorable for the efficient propagation of electromagnetic (EM) pulses, enabling a successful mapping of bottom morphology. In addition, GPR can also help in the evaluation of the glacier bottom properties, the search for internal water floods or underground channels and the recognition of internal glacial structures (as snow layering and crevasses detection). Several GPR applications are reported in the literature on Arctic and Antarctic ice caps and high-elevation icefields of the Tibetan Plateau [18–21]. Growing applications are also recorded in the Alpine context. Del Longo et al. [22] acquired GPR profiles on a Dolomitic glacier to identify the ice bottom depth. Binder et al. [23] investigated two glaciers in the Austrian Alps, for determination of total ice volume and ice-thickness distribution. Merz et al. [24] acquired a dense grid of helicopter-borne GPR, combined with ground-based seismic and geoelectric profiles, on a rock glacier of western Swiss Alps to reconstruct the 3D bedrock topography. Airborne GPR surveying generally overcomes problems in accessibility of the steeper glacier sectors and wave scattering on large boulders in close vicinity of the antenna, which are particularly disturbing in case of debris-cover and rock glaciers. Dossi et al. [25] performed quantitative 3D GPR analysis to estimate the total volume and water content of a glacier in Eastern Alps.

Techniques complementary to GPR surveys have also been applied. Maurer and Hauck [26] reviewed possible additional methods and proposed the joint interpretation of GPR, geoelectric and seismic surveys for analyzing the subsurface conditions of two rock glaciers in the Eastern Swiss Alps, with a special focus on the distribution of ice and water, the occurrence of shear horizons and the bedrock topography. However, geoelectric and reflection seismic methods require long arrays to reach considerable investigation depths and high-resolution imaging of the glacier subsurface, thus complicating field operations. Single-station passive seismic measurements may offer an effective and logistically affordable method to mitigate the problems of multi-channel active surveys, thanks to the use of portable and compact broadband seismometers and no need of energetic active sources. Picotti et al. [27] presented a pioneer study on the application of single-station passive measurements

to determine ice thickness on five Alpine glaciers. Data were processed with the horizontal-to-vertical spectral ratio (HVSR, [28,29]) technique and validated using GPR, geoelectric and active seismic profiling methods. The HVSR method is based on the computation of the ratio between the spectra of the horizontal and vertical components of a triaxial seismic station, to retrieve the fundamental resonance frequency of a site [30,31]. Generally, the resonance originates from the trapping of seismic waves in sites with sufficiently high S-wave impedance contrast (e.g., soft sediments overlying a stiff bedrock). For these sites, a clear peak is found in the HVSR curve. The frequency of this peak is related to the depth of the interface between the two materials [32]. The origin of the HVSR peak is still debated and generally related to a superposition of Rayleigh-, Love- and/or shear-wave resonances [33].

In a glacial environment, the impedance contrast between ice and bedrock is expected to generate resonance phenomena and similar effects on the HVSR curves. Comparing HVSR data with other geophysical imaging results, Picotti et al. [27] showed that the resonance frequency depicted in the HVSR curves is correlated with the ice thickness at the site, in a wide range of ice bottom depths, from few tens of meters to more than 800 m. The reliability of the method mainly depends on the coupling of the sensor at the glacier surface and on the basal impedance contrast. However, beside the intrinsic limitation of the 1D approximation, very few Alpine glaciers lay directly on stiff bedrock for all their extent, while a basal layer of subglacial deposits and debris of significant thickness is generally present below the ablation zone [34], thus attenuating the vertical impedance contrast.

In this study, we acquired both GPR profiles and single-station passive seismic measurements on the terminal lobes of the debris-covered Belvedere Glacier (Macugnaga, NW Italian Alps), with the aim of ice thickness estimation and bottom morphology reconstruction. Very poor and conflicting reference information on ice thickness is indeed available for the glacier, from previous geophysical investigation attempts [35,36]. GPR data acquired in different seasons and with different antennas were processed and manual picking of the bottom surface was performed on the radargrams showing the clearest subsurface imaging. Retrieved ice depths were then spatially interpolated to reconstruct the 3D bottom morphology and compared with previous geophysical results [35,36]. Some single-station passive seismic tests were performed in the same area imagined with GPR profiles. Seismic recordings were processed with the HVSR method and interpreted in the light of the GPR results, for a combined evaluation attempt of the ice bottom morphology and properties.

#### **2. Study Site**

Belvedere Glacier is a debris-covered glacier located NE of the highest peaks of Monte Rosa Massif, in NW Italian Alps (Figure 1). Thanks to the debris cover and the favorable solar exposure, its frontal sectors reach considerably low altitudes and end approximately 2 km W of Macugnaga (VB). Belvedere Glacier represents the terminus of four higher-elevation glaciers (Nordend, Monte Rosa, Signal and Northern Locce Glaciers). Measurements of the Italian Glaciological Committee (CGI-CNR, http://www.glaciologia.it/i-ghiacciai-italiani) carried out in 2006 indicated a surface area of 1.46 km2, a maximum length of 3091 m, spanning in elevation from 2397 m to 1770 m above sea level (a.s.l.), with an average slope of 8◦. The terminal portion of the glacier is bilobate. Both lobes are currently exhibiting a visible retreat. Nowadays, the bigger N lobe has an average length of 650 m, reaching a minimum elevation of about 1810 m a.s.l. (i.e., 40 m higher than measurements of 2006), while the S lobe has a length of 350 m, reaching an elevation of 1840 m a.s.l. at its front. The two lobes are separated by a median morainic relief, hosting the chair lift station and the Belvedere Mountain Hut (Figure 1).

**Figure 1.** (**a**) Geographical location and (**b**) aerial view of the north-eastern glaciers of the Monte Rosa Massif. The main glacial masses are delimited by light blue areas (modified from the Italian Glaciological Committee (CGI-CNR), http://www.glaciologia.it/i-ghiacciai-italiani).

Despite several glaciological studies were carried out on site since the beginning of the 20th century [37,38], the first attempt towards ice thickness characterization is reported by De Visintini [35]. P-wave reflection seismic measurements were indeed carried out on the glacier in 1957 for ice bottom contouring. Explosive sources were adopted, and the active shots were recorded with a 12-channel analog seismograph. Seven local areas of the glacier were surveyed with seismic profiling, from the vicinity of the Zamboni Zappa Mountain Hut (Figure 1), down to the confluence between the two lobes. A global contour map of ice bottom depth was achieved from the interpolation of these measurements, even if some glaciers parts (e.g., the two terminal lobes) were not investigated during this survey.

A later study by VAW-ETH [36] was performed using the GPR technique, with the aim of completing ice bottom characterization in unexplored compartments. A low-frequency GPR instrumentation (USGS Monopulse-radar, with variable central frequency in the range 1–5 MHz) was adopted, with 40 sparse measurements located along 9 transverse profiles covering almost all the glacier length. For each measurement, the time delays between the transmitted EM pulse and the received echoes were recorded. The highest-amplitudes echoes were referred to ice bottom reflections, their depth was estimated from the time delay (considering a constant velocity of the EM pulse in ice) and an interpretation of the bottom morphology was retrieved from data interpolation along the 9 cross-sections.

In the upper part of the glacier (W and SW of Zamboni Zappa Mountain Hut) both radar and seismic results seemed to locate the bedrock at depths of 200–250 m. Progressing northwards, major discrepancies in ice bottom determination arose between the two techniques. Particularly, radar bottom reflectors were found to be located approximately 100 m higher than the seismic reflectors in the central part of the glacier. The two techniques, sensitive to different physical changes in the subsurface, seemed therefore to identify different interfaces. This difference in ice bottom location was interpreted as the result of the presence of a layer of subglacial deposits with relevant thickness (up to

more than 100 m in the central part of the glacier) between ice and bedrock. Glacial deposits revealed transparent to the seismic reflection profiling, but were identified by the radar investigation, probably due to the significant contrast in electrical properties between these sediments and ice.

In a more recent work, Diolaiuti et al. [39] estimated the changes in ice volume between 1957 and 1991 and measured the debris thickness above the whole glacier area. The debris cover was indeed found to slow down the ablation rate of Belvedere Glacier, with respect to other similar glaciers of the same region lacking this peculiarity. A contour map of ice thickness was also presented. However, this map was based just on the results of reflection seismic [35] and did not consider the possible lower ice thicknesses, due to the presence of thick subglacial deposits, as highlighted by the later radar study [36].

In the last decades, geophysical prospection experienced greatly advanced in instrumentation, acquisition and processing techniques. Digital multi-channel recordings and high-sampling big data storage capacities allow for continuous radar profiling with respect to the study of VAW-ETH [36]. Passive seismic acquisitions require lighter instrumentation and no active sources with respect to the measurements of De Visintini [35]. At the light of these advances, new geophysical prospections at the site can therefore help in understanding the discrepancies between the previous works and provide new and more extensive information on the glacier bottom morphology and properties. Despite these considerations, the presence of blocks and debris with variable thickness and lateral distribution at the surface, the rough topography and heterogeneities inside the investigated glacial mass deeply complicate data acquisition and potentially affect the quality of the final results.
