Florent Gimbert

Ugo Nanni

Philippe Roux

A Helmstetter

S Garambois

A Lecointre

A Walpersdorf

B Jourdain

M Langlais

O Laarman

F Lindner

A Sergeant

C Vincent

F Walter

The RESOLVE project: a multi-physics experiment with a temporary dense seismic

published or not. The documents may come

stick slip and englacial fracturing sources as well as novel constraints on the heterogeneous nature of the noise field generated by subglacial water flow and on the link between crevasse properties and englacial seismic velocities. We outline in which ways further work using this dataset could help tackle key remaining questions in the field.

INTRODUCTION

The deployment of large and dense seismic arrays is becoming increasingly common in various geophysical contexts thanks to new technological developments of autonomous wireless seismographs and increased computational power. Spatially dense arrays enhance the characterization of high frequency (>1 Hz) body and surface waves propagating in the subsurface, such as for example in near-surface fault systems exhibiting hundreds to few tens of meters long structures (e.g., the Newport-Inglewood Fault, see [START_REF] Lin | High-resolution 3D shallow crustal structure in Long Beach, California: Application of ambient noise tomography on a dense seismic array[END_REF], and the San Jacinto Fault, see [START_REF] Roux | A methodological approach towards high-resolution surface wave imaging of the San Jacinto Fault Zone using ambient-noise recordings at a spatially dense array[END_REF]). The improved resolution provided by dense arrays increases the completeness of impulsive seismic event catalogs [START_REF] Vandemeulebrouck | The plumbing of Old Faithful Geyser revealed by hydrothermal tremor[END_REF], thus allowing source spatio-temporal dynamics and subsurface structure to be studied in greater details [START_REF] Meng | Detection of small earthquakes with dense array data: example from the San Jacinto fault zone, southern California[END_REF][START_REF] Chmiel | High-sensitivity microseismic monitoring: Automatic detection and localization of subsurface noise sources using matched-field processing and dense patch arraysHigh-sensitivity microseismic monitoring[END_REF]. Dense arrays also help to detect other sources of radiation (e.g. tremor and anthropogenic sources) compared to what is possible with single stations or regional networks [START_REF] Inbal | Localized seismic deformation in the upper mantle revealed by dense seismic arrays[END_REF][START_REF] Li | High-resolution seismic event detection using local similarity for Large-N arrays[END_REF][START_REF] Meng | Detection of small earthquakes with dense array data: example from the San Jacinto fault zone, southern California[END_REF]. Despite its strong potential, dense array monitoring is still limited in the study of glaciers, although it could be used to address some of the key open questions in the field of glaciology.

Glaciers exhibit damage zones created by surface and/or basal crevasses [START_REF] Walter | Using glacier seismicity for phase velocity measurements and Green's function retrieval[END_REF][START_REF] Lindner | Crevasse-induced Rayleigh-wave azimuthal anisotropy on Glacier de la Plaine Morte, Switzerland[END_REF][START_REF] Zhan | Seismic Noise Interferometry Reveals Transverse Drainage Configuration Beneath the Surging Bering Glacier[END_REF][START_REF] Sergeant | On the Green's function emergence from interferometry of seismic wave fields generated in high-melt glaciers: implications for passive imaging and monitoring[END_REF] as well as complex three dimensional structures associated with firn/debris layers, bedrock topography or englacial water conduits [START_REF] Cuffey | The Physics of Glaciers[END_REF]. These features are known to undergo large spatial and temporal changes and to play an important role in glacier dynamics and thermo-dynamics [START_REF] Durand | Impact of bedrock description on modeling ice sheet dynamics[END_REF][START_REF] Scherler | Spatially variable response of Himalayan glaciers to climate change affected by debris cover, 3[END_REF][START_REF] Gilbert | The influence of water percolation through crevasses on the thermal regime of a Himalayan mountain glacier[END_REF]. However, conventional geophysical techniques such as radar sounders capable of resolving englacial structural features (e.g., [START_REF] Evans | Glacier Depth-Sounding from the Air, 5039[END_REF][START_REF] Church | Detecting and characterising an englacial conduit network within a temperate Swiss glacier using active seismic, ground penetrating radar and borehole analysis[END_REF] are not suited for evaluating detailed spatial and temporal changes as well as their effects on the overall glacier behavior. Instead, active and passive surveys using dense seismic arrays may enable accurately monitoring these changes (e.g. using englacial seismic velocities, see [START_REF] Lindner | Crevasse-induced Rayleigh-wave azimuthal anisotropy on Glacier de la Plaine Morte, Switzerland[END_REF] and [START_REF] Preiswerk | Effects of geometry on the seismic wavefield of Alpine glaciers[END_REF] and thus yield unprecedent constraints on glacier structure and temporal evolution.

Glaciers and ice sheets generate a large variety of seismic signals, from impulsive transients to emerging and sustained tremors [START_REF] Podolskiy | Cryoseismology[END_REF][START_REF] Aster | Glacial seismology[END_REF].

Impulsive arrivals from basal stick-slip events have been observed in numerous glaciological contexts [START_REF] Weaver | Seismic Evidence for Discrete Glacier Motion at the Rock-Ice Interface[END_REF][START_REF] Allstadt | Swarms of repeating stick-slip icequakes triggered by snow loading at Mount Rainier volcano[END_REF]Helmstetter, Nicolas, et al., 2015;[START_REF] Lipovsky | Tremor during ice-stream stick slip[END_REF][START_REF] Lipovsky | Glacier sliding, seismicity and sediment entrainment[END_REF][START_REF] Walter | Distributed acoustic sensing of microseismic sources and wave propagation in glaciated terrain, 1[END_REF] and may yield crucial information on glacier basal motion, which exerts a primary control on glacier and icesheet dynamics and the associated eustatic sea level rise [START_REF] Ritz | Potential sea-level rise from Antarctic ice-sheet instability constrained by observations[END_REF][START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF]. The mechanisms giving rise to stick-slip sliding as well as its effect on largescale ice flow, however, remain poorly constrained [START_REF] Lipovsky | Glacier sliding, seismicity and sediment entrainment[END_REF]. Dense array monitoring could allow improving the detection of stick-slip events, yield more accurate inversions of locations and show if and how fast stick-slip asperities migrate.

Impulsive events from englacial fracturing are also commonly observed [START_REF] Neave | Icequakes on the Athabasca Glacier[END_REF][START_REF] Roux | Observation of surface seismic activity changes of an Alpine glacier during a glacier-dammed lake outburst[END_REF][START_REF] Mikesell | Monitoring glacier surface seismicity in time and space using Rayleigh waves[END_REF][START_REF] Podolskiy | Nocturnal Thermal Fracturing of a Himalayan Debris-Covered Glacier Revealed by Ambient Seismic Noise[END_REF][START_REF] Garcia | Joint geodetic and seismic analysis of surface crevassing near a seasonal glacier-dammed lake at Gornergletscher, Switzerland[END_REF], and may help elucidate the role of crevassing in iceberg calving, the disintegration of ice-shelves and the occurrence of serac falls [START_REF] Faillettaz | Evidence of log-periodic oscillations and increasing icequake activity during the breaking-off of large ice masses[END_REF][START_REF] Krug | Combining damage and fracture mechanics to model calving[END_REF][START_REF] Lipovsky | Ice Shelf Rift Propagation and the Mechanics of Wave-Induced Fracture[END_REF]. The improved detection and resolution provided by dense arrays could provide novel constraints on crevasse depth and rupture propagation rates, which are needed to test models (van der [START_REF] Van Der Veen | Fracture mechanics approach to penetration of surface crevasses on glaciers[END_REF][START_REF] Weiss | Subcritical crack propagation as a mechanism of crevasse formation and iceberg calving[END_REF][START_REF] Tsai | A model for turbulent hydraulic fracture and application to crack propagation at glacier beds[END_REF] and thus better understand ice sheet integrity.

Recent seismic investigations have also reported widespread emergent and sustained tremor signals generated by resonances in moulins or water-filled crevasses [START_REF] Helmstetter | Intermediate-depth icequakes and harmonic tremor in an Alpine glacier (Glacier d'Argentière, France): Evidence for hydraulic fracturing?[END_REF][START_REF] Roeoesli | Seismic moulin tremor[END_REF][START_REF] Aso | Seismologically Observed Spatiotemporal Drainage Activity at Moulins[END_REF], subglacial water flow [START_REF] Bartholomaus | Subglacial discharge at tidewater glaciers revealed by seismic tremor[END_REF][START_REF] Eibl | Seismic ground vibrations give advanced early-warning of subglacial floods, 1[END_REF][START_REF] Lindner | Glaciohydraulic seismic tremors on an Alpine glacier[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF] and subglacial sediment transport [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF]. The possibility to calculate physical characteristics of subglacial water flow as well as of subglacial sediment transport from seismic tremor observations [START_REF] Tsai | A physical model for seismic noise generation from sediment transport in rivers[END_REF][START_REF] Gimbert | A physical model for seismic noise generation by turbulent flow in rivers[END_REF][START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Bakker | Field Application and Validation of a Seismic Bedload Transport Model[END_REF] is particularly appealing. These processes play an important role in ice sliding speeds [START_REF] Zwally | Surface Melt[END_REF][START_REF] Schoof | Ice-sheet acceleration driven by melt supply variability[END_REF][START_REF] Tedstone | Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming[END_REF] and bedrock erosion [START_REF] Beaud | Efficacy of bedrock erosion by subglacial water flow[END_REF], and yet hydraulic measurements at the ice bed are notoriously difficult, with traditional approaches such as borehole techniques providing point measurements, only (e.g., [START_REF] Iken | Combined measurements of Subglacial Water Pressure and Surface Velocity of Findelengletscher, Switzerland: Conclusions about Drainage System and Sliding Mechanism[END_REF][START_REF] Iken | Mechanisms of fast flow in Jakobshavns Isbrae, West Greenland: Part I. Measurements of temperature and water level in deep boreholes[END_REF]. Theoretical links between discharge, pressure regime, sediment transport rates and geometry of the subglacial drainage system and seismogenic hydraulic noise sources remain to be more fully tested and dense seismic arrays may provide the necessary spatial extent and resolution for doing so.

Properly evaluating the knowledge gain that dense seismic arrays may provide to address the above-mentioned challenges requires (i) monitoring a glacier that gathers the processes of interest, (ii) deploying instrumentation that covers scales and durations over which significant changes operate, and (iii) acquiring complementary observations to test the seismicallyderived findings and incorporate these into a wider glaciological context. Here we present data and first analysis from a 98 sensor array deployed over 35 days during early spring 2018 on an Alpine Glacier, the Argentière Glacier in the French Alps (Fig. 1). We also provide and analyze key complementary observations from Ground Penetrating Radar (GPR), drone imagery, Global Navigation Satellite System (GNSS) positioning and in-situ instrumentation of basal glacier sliding velocities and subglacial water discharge. We argue that the selected glacier, the time period of investigation as well as the completeness of the present dataset satisfies all of the three above-mentioned conditions. Through application of spectral analysis, template matching, matched-field processing and eikonal wave tomography, we demonstrate that use of the present dataset enhances spatial resolution of basal stick slip activity and near surface crevassing. We further provide novel constraints on the degree of heterogeneity of the seismic noise field generated by subglacial water flow and the variations of englacial seismic velocities. We finally outline how future work using this dataset could help overcome classical observational limitations and address key challenges in the field.

EXPERIMENT DESIGN FIELD SITE

The Argentière Glacier is located in the Mont Blanc Massif (French Alps, 45°55' N, 6°57'E, Fig. 1a) and is the second largest French glacier. It is about 10 km long, covers an area of about 12 km 2 , extends from an altitude of 1700 m above sea level (a.s.l.) up to about 3600 m a.s.l. with an equilibrium line altitude at about 2900 m a.s.l. [START_REF] Vincent | A Nonlinear Statistical Model for Extracting a Climatic Signal From Glacier Mass Balance Measurements[END_REF]. Over the past three decades, glacier total mass balance has been negative (0.7 meter water equivalent loss per year on average over the period 1976-2019), glacier snout has been retreating by several hundreds of meters (815 m retreat between 1990 and 2019), glacier surface elevation has decreased by several tens of meters and glacier surface velocities have decreased by about a factor of two in the ablation zone [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on Glacier d'Argentière, Mont Blanc area, France[END_REF]. The upper part of the glacier is constricted in a typical U-shaped narrow valley where ice sits on granite. The lower part of the glacier is characterized by a sharper incised, V-shaped valley where ice sits on metamorphic rocks [START_REF] Vallon | Contribution à l'étude de la Mer de Glace -Alpes françaises[END_REF][START_REF] Hantz | Waterways, Ice Permeability at Depth, and Water Pressures at Glacier D'Argentière, French Alps[END_REF][START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on Glacier d'Argentière, Mont Blanc area, France[END_REF]. The glacier generally exhibits temperate bed conditions [START_REF] Vivian | Subglacial Cavitation Phenomena Under the Glacier D'Argentière, Mont Blanc, France[END_REF], i.e. basal ice temperature is at the pressure melting point and water flow occurs at the interface as a result of year-round basal melt and summer surface melt [START_REF] Cuffey | The Physics of Glaciers[END_REF].

The monitored site is located in the lower part of the glacier (about 2 km from the glacier front) and at about 2400 m a.s.l. (Fig. 1b). In this area the surface slope is gentle (1-2%) and crevasses are restricted to an area of about 200 m from the glacier sides. The glacier flows at a yearly average velocity of about 60 m.yr -1 in its center, about half of which is due to sliding at the ice-bed interface and the other half to internal ice deformation [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF]. Internal ice deformation likely occurs primarily through ice creep except near the glacier sides where crevasses are large and potentially deep, such that englacial fracturing could also play a role. A strong seasonality is observed in glacier dynamics, with summer (typically May to September) velocities equal to about 1.5 times winter velocities [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF]. This behavior is known to result from melt water input lubricating the ice-bed interface and enhancing basal sliding [START_REF] Lliboutry | Une théorie du frottement du glacier sur son lit[END_REF][START_REF] Lliboutry | General theory of subglacial cavitation and sliding of temperate glaciers[END_REF][START_REF] Cuffey | The Physics of Glaciers[END_REF].

Previous seismic observations at this site report various seismogenic sources associated with surface and intermediate depth crevassing [START_REF] Helmstetter | Intermediate-depth icequakes and harmonic tremor in an Alpine glacier (Glacier d'Argentière, France): Evidence for hydraulic fracturing?[END_REF], basal stick-slip (Helmstetter, Nicolas, et al., 2015), subglacial water flow [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF], and serac instabilities in the glacier front [START_REF] Roux | Microseismic activity within a serac zone in an alpine glacier (Glacier d'Argentière, Mont Blanc, France)[END_REF].

SEISMIC INSTRUMENTATION AND GEOPHYSICAL CHARATERISATION OF GLACIER STRUCTURE AND DYNAMICS

Seismic instrumentation

Sensors of the dense seismic array (red dots in Fig. 1b) are Fairfield ZLand 3-component nodal seismographs with a sampling frequency of 500 Hz (hereafter referred to as nodes). These sensors have a low-corner cut-off frequency of 5Hz, a sensitivity of 76.7 V.m -1 .s -1 (see [START_REF] Ringler | Laboratory Tests of Three Z-Land Fairfield Nodal 5-Hz, Three-Component Sensors[END_REF] for a detailed laboratory analysis of sensor response characteristics) and a typical power autonomy of about 35 days. We deployed the nodes on 24 April 2018 when the glacier was entirely covered by a snow layer of about 3 m thick. We placed the sensors about 40 m apart from each other in the along-flow direction and about 50 m apart in the across-flow direction in order to enable subwavelength analysis in the 4-50 Hz frequency range of interest.

Given that snow melt occurs at an average rate of 2-4 cm.day -1 water equivalent at this location [START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF], we decided to bury the nodes into snow about 30 cm below the surface to ensure the sensors be (i) well coupled to their surroundings and maintained levelled over a week-long time period until snow melt uncovered them and (ii) shallow enough for the GNSS signal to pass through the snow layer and ensure proper reception for time synchronization. Given the limited melt that occurred over the first half of the monitoring period, we had to re-bury the sensors only once over the monitored period, on 11 May 2018. This strategy ensured that little data was lost due to melt-out-induced tilt.

We supplemented the seismic array with one three-component borehole seismic station placed at 5 m below the ice-surface (orange dot in Fig. 1b). This Geobit-C100 sensor was connected to a Geobit-SRi32L digitizer and provides higher sensitivity (1500 V.m -1 .s -1 ), higher frequency sampling (1000 Hz) and a lower low-corner cut-off frequency (0.1 Hz) compared to the nodes. This seismic station is the same as the one used for the two-year long seismic study of [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF].

Recovery of surface and bed digital elevation models from structure from motion surveys and ground penetrating radar

We construct a digital surface elevation model based on a drone geodetic survey that we conducted on September 5, 2018 when the glacier surface was snow free and crevasses could be identified. We used a senseFly eBee+ Unmanned Aerial Vehicle and acquired a total of 720 photos using the onboard senseFly S.O.D.A. camera (20 Mpx RGB sensor with 28 mm focal lens). We generate a digital elevation model of 10-cm resolution (see white contours in Fig. 1b) using differential Global Positioning System (GPS) measured ground control points (see green stars in Fig. 2a) and the Structure for Motion algorithm implemented in the software package Agisoft Metashape Professional version 1.5.2. A detailed description of the processing steps can be found in [START_REF] Brun | Quantifying volume loss from ice cliffs on debris-covered glaciers using high-resolution terrestrial and aerial photogrammetry[END_REF] and [START_REF] Kraaijenbrink | Object-based analysis of unmanned aerial vehicle imagery to map and characterise surface features on a debris-covered glacier[END_REF].

We calculate a crevasse map (black dots in Fig. 2a) based on the surface digital elevation model, which has been shown to be more reliable and precise than using optical/radar images [START_REF] Foroutan | Automatic mapping and geomorphometry extraction technique for crevasses in geodetic mass-balance calculations at Haig Glacier, Canadian Rockies[END_REF]. We first apply a 2D highpass filter with a low cut-off wavelength of 10 m and then define any location with elevation lower than -50 cm as being part of crevasses.

Finally, we apply a 2D median filter with a 1 by 1 m kernel in order to remove artifacts from boulders and moraines.

To establish a digital elevation model of the glacier bed we primarily use Ground Penetrating Radar (GPR) data acquired using a system of two transmitting and receiving 4.2 MHz antennas connected to a time triggered acquisition developed especially for glacial applications by the Canadian company Blue System Integration Ltd [START_REF] Mingo | An integrated lightweight ice-penetrating radar system[END_REF]. The GPR signal processing consists of correcting for source time excitation. We use both dynamic corrections to reproduce a zero-incidence acquisition from data acquired with a 20 m offset between source and receiver (Normal Moveout correction) and static corrections to highlight elevation variations along a profile. We do so using a constant wave velocity of 0.168 m.ns -1 that is typical for ice [START_REF] Garambois | Ground-penetrating radar and surface nuclear magnetic resonance monitoring of an englacial water-filled cavity in the polythermal glacier of Tête Rousse[END_REF]. We then apply a [1-15 MHz] Butterworth band-pass filter followed by a squared time gain amplification to the signal in order to increase signal-tonoise ratio. We show an illustration of the processed GPR data in Fig. 2b, where the direct airwave first arrival is followed by a large reflectivity V-shape pattern reaching 3000 ns around the center of the profile. This latter profile corresponds to the ice/bedrock interface, although its apparent shape is biased by reflections off the closest point on the ice-bed interface rather than off the bed portion directly below the instrument. We correct for this bias by applying a frequency-wavenumber Stolt migration technique [START_REF] Stolt | Migration by Fourier transform[END_REF] and convert time into distance using the constant wave velocity of 0.168 m.ns -1 . We note that prior to migration we add null traces (i.e. with null amplitudes) in places where harsh glacier surface conditions (mainly crevasses) prevented us from acquiring data. As illustrated in Fig. 2c the migration process is effective in correcting the artefacts due to the geometrical variation of the interface along the profile, which now appears smooth and continuous. We then pick the ice-bed reflection (yellow line in Fig. 2c) over all GPR profiles, such that a three-dimensional bed DEM can be reconstructed.

We reconstruct a three-dimensional bed DEM over a larger area than that covered by GPR surveys by incorporating additional constraints like glacier edge elevation as measured from drone imagery (blue area in Fig. 2a) and bed elevations obtained through rock drilling to the ice-bed interface from bedrock excavated tunnels located further down-glacier (purple area in Fig. 2a). Furthermore, we interpolated all data using a kriging method onto a 10 by 10 m grid. From different first onset pickings we estimate that the recovered depth uncertainty is of about 5 m below the seismic array and likely on the order of a few tens of meters outside of the array where observations are sparser.

In Fig. 2a we show the ice thickness map (using 25-m bin contours) as reconstructed based on subtracting the bed DEM from the surface DEM. The glacier bed exhibits a gently dipping valley, with a maximum ice thickness of about 255 m at the center of the seismic array. Glacier thickness decreases relatively sharply on the glacier margins where surface crevasses are observed. We also observe that bed elevation significantly increases down glacier, which results in a decrease by more than 150 m in glacier thickness (Fig. 2d). Beyond these generic characteristics we identify two interesting reflectivity features in the migrated GPR images (yellow ellipses in Fig. 2c) that correspond to localized scattering observed near the surface and a large reflectivity pattern observed just above the deepest portion of the interface. The near surface scattering feature could be caused by deep crevasses, and the deeper feature could be caused by englacial and/or subglacial water conduits as recently proposed by [START_REF] Church | Detecting and characterising an englacial conduit network within a temperate Swiss glacier using active seismic, ground penetrating radar and borehole analysis[END_REF], who made similar GPR observations in a temperate glacier and were able to bolster their interpretation with in-situ borehole observations.

Meteorological and water discharge characteristics

We use air temperature and precipitation measurements obtained at a 0.5 h time step with an automatic weather station maintained by the French glacier-monitoring program GLACIOCLIM (Les GLACIers un Observatoire du CLIMat; https://glacioclim.osug.fr/), which is located on the moraine next to the glacier at 2400 m a.s.l. (green diamond in Fig. 1b).

Precipitation is measured with an OTTPluvio weighing rain gauge. Subglacial water discharge is monitored at a 15 min time step in tunnels excavated into bedrock by the Emossons hydraulic power company, which are located 600 m downstream of the array center (at 2173 m a.s.l.) near the glacier ice fall (see blue star in Fig. 1b).

Temperature generally increases over the instrumented period, from a multi-daily average of about 0° C at the beginning of the measurement period to about 5 °C at the end (Fig. 3a). This drives the general increase in water discharge, which varies from few tenths of m 3 .s -1 to several m 3 .s -1 over the period. Episodic rain events also occur during the instrumented period, but have little to no effect on subglacial discharge likely as a result of the snow cover acting as a water storage buffer [START_REF] Fountain | Water flow through temperate glaciers[END_REF].

Glacier dynamics instrumentation and general features

We evaluate changes in glacier dynamics over the instrumented period by means of two observational methods. The first one is unique to the present site, and consists of basal sliding velocity measurements made continuously in the down glacier serac fall area (see red star in Fig. 1b) by means of a bicycle wheel placed directly in contact with the basal ice at the extremity of an excavated tunnel [START_REF] Vivian | Subglacial Cavitation Phenomena Under the Glacier D'Argentière, Mont Blanc, France[END_REF][START_REF] Moreau | Explications et synthèse des variations de l'hydrographie sous-glaciaire du glacier d'Argentière, Mont-Blanc, grâce aux mesures de l'écoulement du glacier sur son lit rocheux de 1970 à 1998[END_REF]. The wheel is coupled with a potentiometer that retrieves its rotation rate, which is then recorded digitally and converted back to a sliding velocity at a 1-s sampling time with a displacement increment resolution of 0.07 mm. The second type of measurements consists of 4 glacier surface and 1 reference bedrock GNSS stations (yellow stars in Fig. 1b) of type Leica GR25 acquiring the GNSS signals every second. This temporary array is supplemented by a permanent ARGR GNSS station from the RESIF-RENAG network (http://renag.resif.fr) on the bedrock close to the glacier 3 km uphill (yellow star in Fig. 1a). The GNSS antennas on the glacier are installed on 8-m long aluminum masts anchored 4-m deep in the ice and thus emerging about a meter above the snow surface at the beginning of the measurement period. The temporary station placed next to the glacier side provides a reference for validating kinematic GNSS processing approaches, evaluating station positions from every single set of GNSS signal recordings (i.e. every second, as opposed to static processing, which cumulates GNSS signals over a much longer time). We conduct such kinematic processing using the TRACK software ( [START_REF] Herring | Introduction to GAMIT/GLOBK, Department of Earth[END_REF], http://geoweb.mit.edu/gg/docs.php). Our processing chain includes the on-line tool SARI (https://alvarosg.shinyapps.io/sari/) for the removal of outliers that arise from low satellite coverage in the glacier valley and to perform a de-trend and re-trend analysis to estimate and correct for offsets due to manual antenna mast shortening as snow melt progresses. We also correct for multi-path effects induced by GNSS signal reflections from the ground, although we find that those are attenuated by the combination of GPS and GLONASS signals thanks to their different sidereal periods (~24 h for GPS and ~8 days for GLONASS). We finally calculate position time series at a 30-s time step sufficient to capture glacier dynamics and subsequently evaluate three-dimensional velocities by the linear trends of the position components. The horizontal velocity is calculated as 𝑣 ℎ = √𝑣 𝑁 2 + 𝑣 𝐸 2 where 𝑣 𝑁 and 𝑣 𝐸 are the North and East components, respectively.

To facilitate comparison of basal sliding and surface velocity here we smooth both timeseries at a 36-hr timescale (Fig. 3b), since daily down to sub-daily fluctuations in basal sliding velocities are largely affected by unconstrained variations in the local ice roughness in contact with the wheel, as for example, when an entrained rock passes over the wheel. Although basal sliding velocity is to be lower than surface velocity, here both quantities have similar absolute values because the sliding velocity is measured at a place where the glacier is much steeper (25% slope as opposed to 1-2%) and thus driving stress is much larger than at the GNSS locations. We observe an increase in basal sliding velocity from 4.5 mm/h to more than 6 mm/h at the very beginning of the monitored period. Such an acceleration is commonly observed in spring on alpine glaciers [START_REF] Iken | Combined measurements of Subglacial Water Pressure and Surface Velocity of Findelengletscher, Switzerland: Conclusions about Drainage System and Sliding Mechanism[END_REF][START_REF] Mair | Influence of subglacial drainage system evolution on glacier surface motion: Haut Glacier d'Arolla, Switzerland[END_REF][START_REF] Vincent | Sliding velocity fluctuations and subglacial hydrology over the last two decades on Argentière glacier, Mont Blanc area[END_REF] and is known to correspond to water pressurization caused by an increase in water input at the bed due to surface melt water supply, which causes the reduction of friction and thus the enhancement of sliding [START_REF] Lliboutry | Une théorie du frottement du glacier sur son lit[END_REF][START_REF] Lliboutry | General theory of subglacial cavitation and sliding of temperate glaciers[END_REF][START_REF] Iken | The Effect of the Subglacial Water Pressure on the Sliding Velocity of a Glacier in an Idealized Numerical Model[END_REF][START_REF] Schoof | The effect of cavitation on glacier sliding[END_REF][START_REF] Gagliardini | Finite-element modeling of subglacial cavities and related friction law[END_REF]. This acceleration is not seen in the GNSS observations, which could be due to the glacier seasonal acceleration occurring earlier at this location. We also observe one major acceleration event in the location of the dense seismic array occurring between 4

May and 8 May 2018 likely due to the large concomitant increase in water discharge (see blue line in Fig. 3a) causing basal water pressurization [START_REF] Cuffey | The Physics of Glaciers[END_REF].

FIRST RESULTS

SEISMIC NOISE CHARACTERISTICS

We investigate the spatial and temporal variability of seismic power 𝑃 (in dB) across a wide range of frequencies by applying Welch's method [START_REF] Welch | The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms[END_REF] over 4 seconds-long vertical ground motion timeseries (with 50 % overlap) prior to averaging power (in the decibel space) over 15 minutes-long time windows. This two-step strategy allows limiting the influence of impulsive events (which are studied in more details in the next sections) on the seismic power while enhancing that of the continuous background noise [START_REF] Bartholomaus | Subglacial discharge at tidewater glaciers revealed by seismic tremor[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF]. In Fig. 4 we present 1-100 Hz spectrograms (i.e. seismic power at any given frequency and time) over the first half of the instrumented period (15 April to 14 May 2018)

together with timeseries of 2-20 Hz frequency median seismic power at five different stations of the array, four stations located on the four array sides and one located in the array center (see node numbers in Fig. 1b and Fig. S1 for spectrograms across all stations and over the entire frequency range and experimental period). Time periods when sensors tilted beyond their specifications (and thus were no longer deemed functional) as a result of snow melt causing them to be no longer buried are manifested by drastically reduced seismic power across the whole frequency range (see node 6 from 8 May to 11 May 2018). Fortunately, sensor tilt only occurred at a small number of seismic stations (11 out of 98) and during a restricted time duration (less than 2 days on average, see Fig. S1). We also observe that spectrograms do not undergo significant change from prior to after sensor reinstallation on 11 May 2018. This suggests that these are not significantly affected by potential changes in sensor coupling to snow, which is pleasant given that a previous study found that coupling can strongly affect nodes recorded signals [START_REF] Farrell | Persistent Noise Signal in the FairfieldNodal Three-Component 5-Hz Geophones[END_REF].

All stations generally experience similar multi-day (e.g. four days' average, see black lines in Fig. 4) variations in seismic power that are highly correlated with multi-day discharge variations, although seismic power precedes discharge variations by about a day or two. The likely reason is that seismic power is controlled by the hydraulic pressure gradient, which is highest during periods of rising discharge [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF][START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF]. Although shorter term (e.g. diurnal) variations in seismic power are also similar across stations when discharge is low (from 24 April to 28 April and from 30 April to 4 May) and anthropogenic noise dominates [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF], the picture is different at higher discharges when seismic power is dominated by subglacial water flow [START_REF] Nanni | Quantification of seasonal and diurnal dynamics of subglacial channels using seismic observations on an Alpine glacier[END_REF]. On 29 April and from 5 May to 10 May seismic power exhibits pronounced (up to 10 dB) and broad frequency (1-100 Hz)

short time scale (sub-diurnal to diurnal) variations that are particularly marked at certain stations (e.g. node 6 (Fig. 4a), node 44 (Fig. 4c) and node 50 (Fig. 4d)) and not at others (e.g.

node 38 (Fig. 4b) and 95 (Fig. 4e)). We also observe that at certain stations seismic power appears to be continuously or intermittently enhanced within narrow frequency bands. For instance node 38 systematically presents higher seismic power above 20 Hz. These discrepancies suggest that measurements of ground motion amplitude are sensitive to heterogenous and intermittent subglacial water flow, although certain features discussed here could be due to extraneous noise sources associated with sensor coupling [START_REF] Farrell | Persistent Noise Signal in the FairfieldNodal Three-Component 5-Hz Geophones[END_REF], to localized sources other than subglacial fluid flow or to site effects.

DETECTING AND LOCATING BASAL STICK SLIP IMPULSIVE EVENTS USING TEMPLATE

MATCHING

We detect high-frequency (>50 Hz) basal stick-slip events using template matching. This follows a two-step analysis as in Helmstetter et al. (2015). We first build a catalog of events through applying a short-term-average over long-term-average (STA/LTA) detection method [START_REF] Allen | Automatic earthquake recognition and timing from single traces[END_REF] to the continuous high-pass filtered signal (>20 Hz) using a STA time window of 0.1 s and a LTA time window of 1 s. We identify an event when the STA/LTA ratio exceeds a factor of 2. We then manually select all events with short duration (<0.2 s) and high average frequency (>50 Hz) and define groups of events referred to as clusters when their correlation with each other exceeds 0.8. For each cluster, we compute the average waveform to define its "template" signal (using a time window of 0.25 s). We take the sum of seismograms normalized by their peak amplitude and weighted by the square correlation between each event and the template event, iterating this procedure several times until convergence. We visually check that events present distinct P and S wave arrivals and use a polarization analysis to ensure that they are not associated with surface waves (Fig. 5a). We then use the template matching filter method [START_REF] Gibbons | The detection of low magnitude seismic events using array-based waveform correlation[END_REF] to further detect smaller amplitude events not picked with the above strategy but exhibiting a correlation higher than 0.5 with the template signal.

We first conduct the analysis using the borehole station ARG, which has a higher sensor sensitivity, signal-to-noise and sampling rate compared to the nodes. We identify 31 active clusters during the dense array experiment. Interestingly, these clusters constitute a large part of the 46 clusters identified on a much longer period (from December 2017 to June 2018, using the borehole sensor which ran almost continuously, see Fig. 6). Although the amplitude of these signals varies strongly through time (from 1 ⋅ 10 -7 m.s -1 to 4 ⋅ 10 -6 m.s -1 Fig. 6a), waveform characteristics remain strikingly similar (Fig. 6b and6c). All 46 identified clusters exhibit similar characteristics to that shown in Fig. 6, and their activity does not appear to be temporally correlated with each other, nor with external drivers related to meteorology, hydrological or glacier dynamics.

We also apply the template matching algorithm on a subset of 10 nodes covering the whole study area, using the same 31 clusters as previously identified using station ARG. At each time step we compute the correlation coefficient between each template signal and the continuous signal averaged over all nodes and all components (using a time window of 0.35 s instead of 0.25 s in order to match signal duration at all selected nodes). The large number of sensors allows us to lower the correlation threshold from 0.5 to 0.2 while reducing the number of false detections. Indeed, events belonging to different clusters with very different locations can end up being correlated above the detection threshold when using one station, while using several stations this scenario is much more unlikely. We detect 79% more events using the nodes compared to using the station ARG. The newly detected events are mostly smaller amplitude events. Most (83%) of the events detected using ARG are also detected using the nodes. This

shows that increasing the number of sensors allows detecting more events and reducing the number of false detections despite signal-to-noise ratio and sampling rate not being optimum.

We determine the position of the 31 identified clusters by first manually picking on each node the P and S arrival times associated with the event in each cluster that is associated with the largest correlation with the template event, and then inverting for the location of each event and the associated P and S wave velocities assuming velocities are homogeneous and identical for all events. We only consider first arrivals that are usually geometrically predicted to be direct (as opposed to refracted) waves for most sensors and most events. Assuming a simple 1D velocity model with 𝑉 𝑃 =3620 m/s in the ice and 𝑉 𝑃 =4300 m/s in the bedrock (see Fig 10b)

and a glacier thickness of 200 m, the direct wave is faster for epicentral distances shorter than 306 m. Moreover, even when the refracted wave is faster, it is usually less impulsive and has a smaller amplitude than the direct wave. We estimate P and S wave velocities using a grid search inversion with a step of 10 m.s -1 and the Nonlinloc software [START_REF] Lomax | Probabilistic Earthquake Location in 3D and Layered Models[END_REF] to locate clusters. We assume a standard error of arrival times of 2 ⋅ 10 -3 s for P waves, 4 ⋅ 10 -3 s for S waves and of 3.5 ⋅ 10 -3 s for calculated travel times. We can see in Fig. 5a that the picked arrival times (black circles) are consistent with the computed travel times (green lines). The root-mean-square error for this event is 2.4 ms, which corresponds to about one sample (2 ms).

We show the locations of basal icequakes versus depth in an average transverse section in Fig. 7a and on a two-dimensional map in Fig. 8. They are mainly located in the down-glacier part of the array and in the central part of the glacier or near the right side, while there is no event observed towards the left side. Icequake depths range between 80 m and 285 m, and are in good agreement with the bedrock topography estimated from the radar profiles. Uncertainty on absolute source depth is on the order of 10 m (see errorbars in Fig. 7a), and the estimated seismic wave velocities of 𝑉 𝑃 =3620 m.s -1 and 𝑉 𝑆 =1830 m.s -1 (Fig. 7b) are in good agreement with velocities measured on other alpine glaciers [START_REF] Podolskiy | Cryoseismology[END_REF]. 𝑉 𝑆 is much better constrained by the data compared to 𝑉 𝑃 (Fig. 7b).

SYSTEMATIC LOCATION OF EVENTS USING MATCHED-FIELD PROCESSING

Contrary to in the previous section where a priori constraints on waveform characteristics and wave velocity are used to target basal stick-slip events, we next test location of a wide range of seismic events generated by impulsive or emergent sources with no a priori knowledge on waveform characteristics and minimal a priori knowledge on medium properties. The rationale is that the limited a-priori knowledge for source identification is balanced by the high spatial and temporal resolution provided by array processing techniques, which may provide spatial or temporal characteristics facilitating source identification [START_REF] Vandemeulebrouck | The plumbing of Old Faithful Geyser revealed by hydrothermal tremor[END_REF][START_REF] Chmiel | High-sensitivity microseismic monitoring: Automatic detection and localization of subsurface noise sources using matched-field processing and dense patch arraysHigh-sensitivity microseismic monitoring[END_REF].

We conduct Matched-Field Processing (MFP), which consists of recursively matching a synthetic field of phase delays between sensors with that obtained from observations using the Fourier transform of time-windowed data. We obtain the synthetic field from a source model with a frequency-domain Green's function that depends on 4 parameters, which are the source spatial coordinates 𝑥, 𝑦 and 𝑧 and the medium phase speed 𝑐. MFP output is normalized to range from 0 to 1 with higher values corresponding to better matches between modelled and observed signal phases, and therefore a higher confidence in true source location. Here we use a spatially homogeneous velocity field within the glacier, which the advantage of a fast-analytical computation, although it also results in a higher degree of ambiguity between 𝑧 and 𝑐. Contrary to classical beamforming techniques in which a planar wave front is often assumed, our MFP approach considers spherical waves and allows locating sources closer to and within the array. To build a large catalog of events, we apply MFP over short time windows of 1-s with 0.5-s overlap, across 16 frequency bands of ±2 Hz width equally spaced from 5 to 20 Hz and over the entire study period. Calculating source locations over such a large number of windows requires minimizing computational cost. We do so by using a minimization algorithm that relies on the downhill simplex search method (Nelder-Mead optimization) of [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] and [START_REF] Lagarias | Convergence Properties of the Nelder--Mead Simplex Method in Low Dimensions[END_REF] instead of using a multidimensional grid search approach. As the exploration of the solution space is characterized by a certain level of randomness, we maximize the likelihood that our minimization technique finds a global minima and thus the dominant source over the considered time window through (i) starting the optimized algorithm from a set of 29 points located at a depth of 250 m inside and near the array (see black crosses in Fig. 9d) with a starting velocity 𝑐=1800 m.s -1 and (ii)

taking the highest MFP output out of the 29 inversions found after convergence.

In Fig. 9b,c we present two examples of events located inside and outside the array and associated with a high MFP output of 0.92. The half-size of the focal spot in the MFP output field gives a measure of the location uncertainty [START_REF] Rost | Array Seismology: Methods and Applications[END_REF], which is about 10 m for events located inside the array and can increase up to 40 m when for events up to 100 m away from the array edges. Gathering all sources over one continuous day of record, we find that the associated MFP output distribution exhibits a heavy tail towards high values (red area in Fig. 9a for an example at 13 Hz). Such a heavy tail is not obtained for a random field, in which case MFP output exhibits a distribution shifted towards almost one order of magnitude lower values. This suggests that most identified sources correspond to real and detectable seismic events. Well resolved seismic events with MFP outputs higher than 0.8 are located near the surface and delineate crevasse geometries, such that they likely correspond to englacial fracturing (red dots in Fig. 8). Few (less than one percent) of these events are however located outside of the glacier and likely correspond to rock falls. Typical waveforms associated with englacial fracturing events are dominated by surface waves arrivals (Fig. 5b), although P waves arrivals as well as arrivals showing hyperbolic moveout (black arrows in Fig. 5b) are also distinguishable. Although P waves arrivals associated with surface crevassing events are not commonly observed [START_REF] Walter | Moment Tensor Inversions of Icequakes on Gornergletscher, SwitzerlandMoment Tensor Inversions of Icequakes on Gornergletscher, Switzerland[END_REF][START_REF] Helmstetter | Intermediate-depth icequakes and harmonic tremor in an Alpine glacier (Glacier d'Argentière, France): Evidence for hydraulic fracturing?[END_REF], their observation here may result from improved detection thanks to the dense seismic array.

Arrivals showing hyperbolic moveout likely correspond to reflected waves at the glacier/bedrock interface.

USING EVENT CATALOGS FOR STRUCTURE INVERSION

Dense-array techniques for seismic imaging often involve interferometry analysis on continuous seismic noise. Such techniques however require an equipartitioned wavefield inherited directly from homogenously distributed noise sources and/or indirectly from sufficiently strong scattering [START_REF] Lobkis | On the emergence of the Green's function in the correlations of a diffuse field[END_REF]Fichtner et al., 2019). These conditions strongly limit the applicability of such techniques on glaciers where sources are often localized and waves in ice are weakly scattered [START_REF] Sergeant | On the Green's function emergence from interferometry of seismic wave fields generated in high-melt glaciers: implications for passive imaging and monitoring[END_REF]. An alternative way is to use localized and short-lived sources with known positions [START_REF] Walter | Using glacier seismicity for phase velocity measurements and Green's function retrieval[END_REF] as those previously identified using our systematic MFP technique, which are numerous and more evenly distributed in space (Fig. 8).

We consider the catalog of sources associated with MFP outputs larger than 0.6, located near the surface (z<10m) and close to the array (within a radius of 400 m from the array center).

With these criteria our catalog includes about 10 6 sources gathered over the 35 days of continuous recordings. In order to further demonstrate that our MFP calculations yield reliable velocities (i.e. the ambiguity between 𝑧 and 𝑐 is limited for these sources), we use the velocities given from our MFP calculation (which for shallow events recover the dominant surface waves) to construct dispersion curves, as opposed to classical f-k analysis [START_REF] Capon | High-resolution frequency-wavenumber spectrum analysis[END_REF]. We infer surface wave phase velocity at each frequency between 3.5 Hz and 25 Hz by fitting a Gaussian function to the probability density distributions of velocities in each frequency bin, and taking the center of the Gaussian function as the most representative velocity in that frequency bin (see Fig. 10a (inset) for an example at 13 Hz). We note that the presently constructed dispersion curve is similar to the one that would be obtained using a classical f-k analysis (not shown). We find that surface wave velocity increases gently from 1560 m.s -1 to 1630 m.s -1 as frequency decreases from 25 Hz down to 7 Hz, and then increases sharply up to 2300 m.s -1 as frequency decreases down to 3.5 Hz. These observations can be reproduced using a three-layer one-dimensional elastic model (using the Geopsy package, [START_REF] Wathelet | Geopsy: A User-Friendly Open-Source Tool Set for Ambient Vibration Processing[END_REF]) that incorporates a gentle velocity increase (from 1670 to 1720 m.s -1

for 𝑉 𝑠 ) at 40 m depth and a drastic velocity increase (from 1720 to 2800 m.s -1 for 𝑉 𝑠 ) located between 200 and 220 m depth (Fig. 10b). These values were obtained by trial and error tests.

The slightly slower velocities and density within the first 40-m deep layer may be due to surface crevasses, and are consistent with surface events being associated with smaller P wave velocities than those associated with stick-slip events at the ice/bedrock interface (Fig. 5). The 200-to 220-m deep drastic discontinuity results from the ice/bedrock interface, consistent with the radar-derived average glacier thickness beneath the seismic network (Fig. 2a).

We go one step further and perform two-dimensional surface wave inversions from eikonal wave tomography [START_REF] Roux | The San Andreas Fault revisited through seismicnoise and surface-wave tomography[END_REF][START_REF] Lin | High-resolution 3D shallow crustal structure in Long Beach, California: Application of ambient noise tomography on a dense seismic array[END_REF][START_REF] Mordret | Nearsurface study at the Valhall oil field from ambient noise surface wave tomography[END_REF]. We first extract ~200,000 Rayleigh wave travel times using the best (associated with MFP outputs larger than 0.9) seismic events and then perform a simple linear inversion for the slowness (starting from a homogeneous initial model with a phase velocity of 1580 m.s -1 , see Fig. 10a) assuming straight rays as propagation paths and an a-priori error covariance matrix that decreases exponentially with distance over 10 m. The weight of the spatial smoothing is chosen at the maximum curvature of the standard trade-off analysis (L-curve) based on the misfit value [START_REF] Hansen | The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems[END_REF], and the inversion produces a residual variance reduction of ~98% relative to the arrival times for the homogeneous model. In Fig. 8 we show the Rayleigh wave phase velocity maps obtained as a result of the travel-time inversion on a regular horizontal grid with steps of 5 m and using 13-Hz Rayleigh waves, which have largest sensitivity between 20 m and 60 m depth (Fig. 10c) according to kernel sensitivity computations performed on the three layer elastic model (Fig. 10b) using the code of [START_REF] Herrmann | Computer Programs in Seismology: An Evolving Tool for Instruction and Research[END_REF]. We observe that locations with higher crevasse density are generally associated with lower phase velocities, as observed in the left glacier side and in the down glacier part of the array. This observation is however not systematic, since high velocities are also observed in the right glacier side and in the up glacier part of the array where crevasses are also present. This could be explained by shallower crevasses or by crevasse orientations, which affect different wave propagation directions in these regions. This latter potential source of bias could be investigated by explicitly accounting for anisotropy in the tomography inversion scheme [START_REF] Mordret | Nearsurface study at the Valhall oil field from ambient noise surface wave tomography[END_REF].

DISCUSSION

INTERPRETING SPATIAL AND TEMPORAL VARIATIONS IN GROUND MOTION AMPLITUDES

Although our seismic array observations generally exhibit spatially homogenous multi-day changes in seismic power, there exist specific times when changes in seismic power are spatially heterogeneous (Fig. 4). A surprising observation is that these heterogeneous changes are observed down to the lowest frequencies (3 to 10 Hz) associated with wavelengths larger than the inter-station spacing, such that the observed spatial heterogeneity is unlikely solely caused by wave attenuation. It remains to be investigated as to which processes mainly cause the observed spatial variability in seismic power. Punctual sources identified from the MFP analysis could be used to investigate the respective control of wave attenuation, wave scattering and site effects on amplitude field heterogeneity and its potential dependency on site attributes like crevasse density, glacier thickness or snow layer thickness. Full waveform modelling combined with wave polarity analysis could also be conducted in order to further understand how wave focusing in the near field domain as well as source heterogeneity and directivity may cause heterogenous amplitude wavefields. Incorporating these constraints into an improved model describing the control of both source and wave propagation physics on the seismic wave amplitude field [START_REF] Gimbert | Subseasonal changes observed in subglacial channel pressure, size, and sediment transport[END_REF] could allow using our dense array observations to infer the spatial variability in subglacial water flow parameters such as subglacial channel size and pressure.

PHYSICS OF STICK SLIP EVENTS

We demonstrate that dense array observations provide enhanced resolution on stick-slip motion. Applying template matching on an array of sensor as opposed to on a single station enables detecting many more events within clusters. Using the whole array for location inversions also allows significantly reducing location uncertainties. Future studies may focus on applying template matching across all sensors of the array in order to detect and locate more events within clusters and potentially more clusters. With our present analysis we find that events are all located in the down-glacier part of the array and in the central part of the glacier or near the right side, while there is no event observed towards the left side (Fig. 7 and8). This provides further observational support that specific bed conditions (e.g. water pressure, bed shear stress, bed roughness, bed topography, carried sediments) are necessary for these events to occur [START_REF] Zoet | The effects of entrained debris on the basal sliding stability of a glacier[END_REF][START_REF] Lipovsky | Glacier sliding, seismicity and sediment entrainment[END_REF]. Further insights into the physics controlling the spatio-temporal dynamics of these events could be gained by performing relative event location within each cluster using double-differences [START_REF] Waldhauser | A Double-Difference Earthquake Location Algorithm: Method and Application to the Northern Hayward Fault, California[END_REF] instead of simply inferring single cluster locations as presently done.

These improvements could allow identifying whether or not stick-slip asperities migrate. crevasses, englacial conduits) acting as scatterers. One could also combine MFP with eigenspectral decomposition to reveal weaker noise sources that would otherwise be hidden within the background noise [START_REF] Seydoux | Detecting seismic activity with a covariance matrix analysis of data recorded on seismic arrays[END_REF]. Additional constraints for seismic imaging may also be provided through identifying specific events generating indirect arrivals of particular interest for structural analysis, such as in bed-refracted waves shown in Fig. 5

Supplemental Material (Main Page, Tables , andFigures)

Click here to access/download;Supplemental Material (Main Page, Tables, and Fig. S1: Spectrograms calculated at all stations over the whole period. The lower panel shows a map view of the node number nomenclature. The reader should refer to Fig. 1 for absolute positioning of the dense array. The color scale represents seismic power, as calculated in Fig. 4, and is identic across all panels.

FIGURES Fig. 1 :

 1 FIGURES

Fig. 9 :

 9 Fig. 9: (a) Probability distribution functions (PDF) of MFP outputs obtained at 13 Hz when applying MFP on one day (1 May 2018) of real data (red) and on a numerically-generated random wavefield (green). The bottom panels (b) and (c) show the phase fields observed over a 1-s time window at 13 Hz for two selected events. Locations obtained from MFP using our minimization process are shown by the pink arrow/green crosses, while the contour lines show 0.1 and 0.8 MFP outputs iso-contours calculated by applying a grid search over the glacier

  USING MFP TO RETRIEVE SOURCES AND STRUCTURAL PROPERTIESSystematic MFP analysis with adequate parametrization opens a route to continuous, automatic, and statistics-based monitoring of glaciers. A wide diversity of seismic sources may be identified and studied separately with this technique by scanning through different values of MFP outputs. High MFP outputs may be used to study the dynamics of crevasse propagation with particularly high spatio-temporal resolution. Such observations may allow to better understand the underlying mechanisms associated with crack propagation, in particular through providing an opportunity to better bridge the gap between laboratory and theoretical material physics of crack propagation (van der[START_REF] Van Der Veen | Fracture mechanics approach to penetration of surface crevasses on glaciers[END_REF][START_REF] Weiss | Subcritical crack propagation as a mechanism of crevasse formation and iceberg calving[END_REF] and crevasse propagation under realistic glacier conditions in which water is expected to play an important role (van der[START_REF] Veen | Fracture propagation as means of rapidly transferring surface meltwater to the base of glaciers[END_REF]. Lower MFP outputs may be used to locate spatially distributed

sources generating coherent signals over only a limited spatial extent. These distributed sources may include tremor sources (e.g. water flow) or various glacier features (e.g.

  Wehrlé and B. Urruty for their help in the field, as well as Martin, Fabien and Christophe for mountain guiding the group.Induced Acceleration of GreenlandIce-Sheet Flow, Science 297, no. 5579, 218-222, doi: 10.1126/science.1072708.

	ground penetrating radar, drone imagery, GNSS positioning and in-situ instrumentation of Palanstjin, A.
	basal glacier sliding velocities and subglacial water flow discharge. We show that a wide range
	of glacier sources and structure characteristics can be extracted through multiple seismic
	processing techniques such as spectral analysis, template matching, matched-field processing
	and eikonal wave tomography. Future studies focusing more specifically on each aspect of the
	herein presented observations may yield novel quantitative insights into spatio-temporal
	changes in glacier dynamics and structure.
	(black
	arrows).
	SUMMARY
	We present a dense seismic array experiment made of 98 3-component seismic stations
	continuously recording during 35 days in early spring 2018 on the Argentière Glacier, French
	Alps. The seismic dataset is supplemented by complementary observations obtained from

ACKNOWLEDGEMENTS

This work has been supported by a grant from Labex OSUG (Investissements d'avenir -ANR10 LABX56). IGE and IsTerre laboratories are part of Labex OSUG (ANR10 LABX56). Complementary funding sources have also been provided for instrumentation by the French "GLACIOCLIM (Les GLACIers comme Observatoire du CLIMat)" organization and by l'Agence Nationale de la recherche through the SAUSSURE (, ANR-18-CE01-0015) and SEISMORIV (ANR-17-CE01-0008) projects. We thank C. Aubert, A. Colombi, L. Moreau, L. Ott, I. Pondaven, B. Vial, L. Mercier, O. Coutant, L. Baillet, M. Lott, E. LeMeur, L. Piard, S. Escalle, V. Rameseyer, A.

DATA AND RESOURCES

Raw seismic data can be found at: Roux, P., Gimbert, F., & RESIF. (2021). Dense nodal seismic array temporary experiment on Alpine Glacier of Argentière (RESIF-SISMOB) [Data set]. RESIF -Réseau Sismologique et géodésique Français. https://doi.org/10.15778/RESIF.ZO2018 (see also link http://seismology.resif.fr/#NetworkConsultPlace:ZO%5B2018-01-01T00:00:00_2018-12-31T23:59:59%5D). Processed data used in this paper can be found at: https://doi.org/10.5281/zenodo.3701519 for meteorological, subglacial water flow discharge and glacier sliding speed data https://doi.org/10.5281/zenodo.3971815 for bed thickness, surface elevation, nodes positions, crevasses positions, surface velocity, noise PSDs, event occurrences and locations derived from template matching for stick-slip events and MFP for englacial fracturing events https://doi.org/10.5281/zenodo.3556552 for drone orthophotos