Antoine Guillemot

Agnès Helmstetter

Éric Larose

Laurent Baillet

Stéphane Garambois

Raphaël Mayoraz

Reynald Delaloye

Seismic monitoring in the Gugla rock glacier (Switzerland): ambient noise correlation, microseismicity and modeling

Keywords: Rock glacier, environmental seismology, seismic noise, coda waves, micro-seismicity

A network of seismometers has been installed on the Gugla rock glacier since October 2015 to estimate seismic velocity changes and detect micro-seismicity. These two processes are related to mechanical and structural variations occurring within the rock glacier. Seismic monitoring thus allows a better understanding of the dynamics of rock glaciers throughout the year. We observed seasonal variations in seismic wave velocity and micro-seismic activity over the three years of the study. In the first part of our analysis, we used ambient noise correlations to compute daily changes of surface wave velocity. In winter, seismic wave velocities were higher, probably due to refreezing of the active permafrost layer and cooling of the uppermost permafrost layers, leading to increased overall rigidity of the medium. This assumption was verified using a seismic model of wave propagation that estimates the depth of P-and S-wave velocity changes from 0 down to 10 m. During melting periods, both a sudden velocity decrease and a decorrelation of the seismic responses were observed. These effects can probably be explained by the increased water content of the active layer. In the second part of our study, we focused on detecting micro-seismic signals generated in the vicinity of the rock glacier. This seismic activity (micro-quakes and rockfalls) also exhibits seasonal variations, with a maximum in spring and summer, which correlate with acceleration of the rock glacier's displacement rate. In addition, we observed short bursts of micro-seismicity, both during snowfall and during rapid melting periods. All these observations could be included in systems monitoring permafrost and slope destabilization, to complement existing records.

Introduction

Rock glaciers are the most prominent features in alpine permafrost [START_REF] Schoeneich | Velocity Changes of Rock Glaciers and Induced Hazards[END_REF]. They appear as lobate or tongue-shaped bodies of frozen debris with interstitial ice, a varying amount of interlaced or isolated ice lenses, or a core of massive ice [START_REF] Arenson | Borehole deformation measurements and internal structure of some rock glaciers in Switzerland[END_REF] [START_REF] Whalley | Rheological models of active rock glaciers: Evaluation, critique and a possible test[END_REF] [START_REF] Jansen | Rock glacier dynamics: Stick-slip motion coupled to hydrology[END_REF]. These geomorphological features are the result of permafrost creep, which corresponds to the steadystate deformation of the ice or the ice-rock mixture [START_REF] Haeberli | Permafrost creep and rock glacier dynamics[END_REF] [START_REF] Jones | Rock glaciers and mountain hydrology: A review[END_REF]. On a large temporal scale, rock glaciers can be used as a climate indicator as their spatial distribution correlates with the permafrost limit [START_REF] Jansen | Rock glacier dynamics: Stick-slip motion coupled to hydrology[END_REF].

In many high mountain regions, warming of perennially frozen ground significantly influences slope stability due to accelerated movement in warmer conditions [START_REF] Kääb | On the response of rockglacier creep to surface temperature increase[END_REF][START_REF] Staub | Using Near-Surface Ground Temperature Data to Derive Snow Insulation and Melt Indices for Mountain Permafrost Applications[END_REF]. As a result, rock glaciers may locally play an important role in rock mass transport and gravitational processes [START_REF] Jansen | Rock glacier dynamics: Stick-slip motion coupled to hydrology[END_REF].

Active rock glaciers are creeping landforms composed of debris supersaturated with ice, moving downslope at a rate of a few cm/yr to several m/yr [START_REF] Schoeneich | Velocity Changes of Rock Glaciers and Induced Hazards[END_REF]. When their dynamics are altered (rapid acceleration, development of specific morphological features, including cracks), active rock glaciers can be considered to be destabilized [START_REF] Marcer | Evaluating the destabilization susceptibility of active rock glaciers in the French Alps[END_REF]. Destabilization and catastrophic collapses of rock glaciers are facilitated by increasing temperatures of permafrost. Their thermal response to climate has been studied over significant periods [START_REF] Bodin | Two decades of responses (1986-2006) to climate by the Laurichard rock glacier, French Alps[END_REF], and shows a high sensitivity when ground temperatures stays near 0 °C, as recorded over at least the last two or three decades [START_REF] Kääb | On the response of rockglacier creep to surface temperature increase[END_REF]. A destabilization event can be preceded by high horizontal velocities, high front advance rates or crevasse-like cracks opening [START_REF] Roer | Observations and considerations on destabilizing active rock glaciers in the European Alps[END_REF].

According to [START_REF] Delaloye | Rapidly moving rock glaciers in Mattertal[END_REF], the main causes of destabilization are:

-Morphological forcing, (steep slope angle, higher shear stress, convex bedrock topography, local thinning of the rock glacier); -Mechanical forcing, due to external factors (overloading, landslides, rock avalanche and rockfalls); -Thermal forcing, affecting the rheological properties of ice (ductility), increasing pore water pressure that decreases friction at depth, where the creeping process develops (shear horizon).

An acceleration of the rock glacier's movement can increase the rate of mass wasting events at its front, as well as altering the front geometry, depending on the material available and the topography. Although large ruptures of rock glaciers are not known to occur, they may collapse, thus threatening large areas downslope, particularly when the rock glacier front overlooks steep slopes or torrential catchment areas. Secondary processes mobilizing the debris released -such as torrential and debris flows-can further increase the hazard related to rock glaciers.

An example of a huge landslide of frozen debris including part of a rock glacier is provided by the collapse in August 2006 of the frontal part of the Bérard rock glacier, in the southern French Alps. Some 500 000 m 3 of material slid downslope one month after the first signs of destabilization were observed. The collapse was probably favored by the slope convexity in the area of the scarp, the nature of the fine debris composing the majority of the rock glacier, and the softness of the bedrock. The final collapse may have been triggered by the hydro-niveo-meteorological conditions during the preceding hours or weeks [START_REF] Bodin | The 2006 Collapse of the Bérard Rock Glacier (Southern French Alps)[END_REF].

With these issues in mind, continuous meteorological surveys combined with GPS geodetical recordings are used to link climate evolution to long-term topographic reactions. These data then contribute to early warning systems for destabilization in the short-term. Although rock glacier deformation is distributed over the whole thickness of the rock glacier, surface displacement is currently the easiest variable to characterize [START_REF] Wirz | Short-term velocity variations at three rock glaciers and their relationship with meteorological conditions[END_REF][START_REF] Buchli | Kinematic investigations on the Furggwanghorn Rock Glacier, Switzerland[END_REF]. In addition, obtaining information about the deeper structure and deformations of a rock glacier is highly challenging. Geotechnical data collected from field campaigns are very valuable but require a lot of work and resources to collect and treat. Indeed, it is difficult to drill boreholes and install inclinometers into a rock glacier to monitor local temperatures, water content and deformation [START_REF] Bodin | Multi-Annual Kinematics of an Active Rock Glacier Quantified from Very High-Resolution DEMs: An Application-Case in the French Alps[END_REF]. Also, most geophysical investigations (georadar, seismic reflection tomography, electric resistivity tomography) require onerous field resources to provide valuable results, and cover only one period of observation. Even when several geophysical methods are combined to provide complementary data, interpretation can be very difficult [START_REF] Kneisel | Advances in geophysical methods for permafrost investigations[END_REF]. In particular, the highly heterogeneous nature of rock glaciers (voids, ice, rock debris) tends to attenuate seismic and radar waves, masking the bulk structures [START_REF] Maurer | Geophysical imaging of alpine rock glaciers[END_REF]. However, passive seismic monitoring using relatively low frequency (2-20 Hz) ambient noise and micro-seismicity suffers less from scattering and attenuation in these conditions than active high frequency seismics. Thus, passive monitoring can be used to detect changes in the rock glacier properties at depth and to characterize the damage and deformation processes occurring within the medium, in complement to other investigations (kinematics, geophysics).

Seismic monitoring provides information relating to:

-The evolution of the rock glacier's bulk elastic properties, by estimating changes in seismic surface wave velocity computed from ambient noise correlation processing;

-Damage and fracturing processes, as well as mass wasting events, by detecting seismic signals emerging from the noise that are generated by these processes.

As a result of mechanical waves propagating in a 3D continuous medium, seismic wave velocities depend on elastic features (mostly rigidity and density), and thus are sensitive to changes in water pressure, freezing, and damage. Thanks to recent progress in methodology and computational technology, seismic data can now track very minute changes in the subsurface, induced by the environment (see for instance the growing interest in environmental seismology [START_REF] Larose | Environmental seismology: What can we learn on earth surface processes with ambient noise[END_REF]). The approaches used consist in placing natural (passive) or artificial (active) seismic sources to study alterations to surface wave propagation related to external environmental changes (thermal and hydro-meteorological evolution, erosion processes).

As a passive seismic method, Coda Wave Interferometry has been successfully used to monitor tiny changes in several subsurface contexts, including volcanoes [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia --2006[END_REF], landslides [START_REF] Mainsant | Ambient seismic noise monitoring of a clay landslide: Toward failure prediction: SEISMIC NOISE MONITORING OF A LANDSLIDE[END_REF], and permafrost [START_REF] James | Improved moving window cross-spectral analysis for resolving large temporal seismic velocity changes in permafrost[END_REF]. This technique is based on ambient noise correlation and consists in comparing correlograms over time (to a reference averaged over a long or stable period, for example).

Coda waves arrive later than ballistic waves, and mostly result in partitioned surface and body waves, which have been diffused by scattering at depth [START_REF] Obermann | Lapse-time-dependent coda-wave depth sensitivity to local velocity perturbations in 3-D heterogeneous elastic media[END_REF]. The basic idea behind this method is that the ambient noise cross-correlations between two sensors (correlograms) represent the impulse response between them as if one was a source (or a fraction of this impulse response that occurs to be stable from one day to another). Comparing the coda part of correlograms thus makes it possible to analyze changes to seismic waveforms that have extensively surveyed the medium by following a complex path. Repeated measurements make this method compatible with probing the properties of the soil over time. Since sensors are fixed and ambient noise sources (due to natural phenomena or human activity) are statistically stable over time and space [START_REF] Hadziioannou | Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments[END_REF], alterations to noise cross-correlations are only due to elastic changes in the subsurface surveyed. These cross-correlation changes may be due to variations in surface wave velocity, but can also be induced by structural changes affecting the waveforms. Using the stretching method (described further below, in Methods), we can separately quantify both effects [START_REF] Larose | Environmental seismology: What can we learn on earth surface processes with ambient noise[END_REF].

In addition, seismic monitoring of landslides or glaciers can detect, locate and characterize seismic activity induced by several phenomena, such as rock falls, debris flows, crack propagation, basal sliding, or meltwater runoff [START_REF] Helmstetter | Seismic monitoring of Séchilienne rockslide (French Alps): Analysis of seismic signals and their correlation with rainfalls[END_REF], Helmstetter et al. 2015a,b, Gimbert et al. 2016). Properties of the seismic events (frequency, duration, location, temporal distribution, etc.) can be used to discriminate the different mechanisms of deformation (creeping, stick-slip, fracture) [START_REF] Provost | Towards a standard typology of endogenous landslide seismic sources[END_REF]. The temporal evolution of seismic activity may be related to meteorological variations and displacement rate, and can thus be used in studies to characterize rock glacier dynamics.

This study aims at further our understanding of how environmental forces drive rock glacier dynamics. Our long-term goal is to define seismic observables predicting mass wasting events or destabilization, and to implement them in early warning systems.

Context and instrumentation of the Gugla rock glacier

In the Valais Alps (southwestern Switzerland), a large number of active rock glaciers have been regularly investigated over the last decade ( [START_REF] Delaloye | Rapidly moving rock glaciers in Mattertal[END_REF][START_REF] Merz | Multidisciplinary geophysical investigations over an alpine rock glacierMultigeophysics over a rock glacier[END_REF][START_REF] Wirz | Short-term velocity variations at three rock glaciers and their relationship with meteorological conditions[END_REF][START_REF] Buchli | Kinematic investigations on the Furggwanghorn Rock Glacier, Switzerland[END_REF]).

Among them, the Gugla (also called Gugla-Breithorn or Gugla-Bielzug) rock glacier is located above the Bielzug torrent, on the slope between the Breithorn and Gugla summits. Its tongue-shaped morphology covers about 130 m in width, 600 m in length, and is up to 40 m thick in its terminal part. Most of the rock glacier is inclined at less than 20°, although its terminal part (roughly the last 100-150 m) is steeper. From GPS and InSAR data, annual surface velocity was seen to increase continuously from 2.8 m/yr in August 2007 to more than 5 m/yr in November 2010, particularly for the southern terminal part [START_REF] Delaloye | Rapidly moving rock glaciers in Mattertal[END_REF]) [START_REF] Wirz | Short-term velocity variations at three rock glaciers and their relationship with meteorological conditions[END_REF]. Since 2010, surface velocities have been about 5 m/yr at the front, with a peak of activity in the southern part culminating, in 2013, at a displacement of more than 15 m/yr. This increase in velocity has also propagated to the rooting zone (from 0.6 m/yr in 2008 to 2 m/yr in 2018). About 150 m upstream of the front (where the topography changes) a fractured area has developed with visible transverse scars. This destabilization could be the result of a combined effect of the topography, of the landform geometry (larger and thicker below the topographical step) and of the increasing permafrost temperatures over the last two decades [START_REF] Delaloye | Rapidly moving rock glaciers in Mattertal[END_REF].

The surrounding rock walls are mainly composed of highly fractured paragneiss, from which rocks regularly fall and supply the rock glacier. The steep convex slope dominates a torrential gully, allowing efficient sediment connectivity with a 0.64 km 2 catchment basin. Yearly from 2012 to 2018 (except 2017), one or more debris flows were triggered from an area immediately downslope of the rock glacier front and reached the main valley close to the village of Herbriggen (St Niklaus, VS). These events threaten railways and roads between Visp and Zermatt (Figure 1a). Interestingly, the Europaweg hiking trail passes directly through the frontal zone of rock glacier, and was closed to the public from 2019. Since the risk of runout onto the village remains current, embankments and a reservoir have been erected to mitigate and manage material carried by torrential flows. The volumes involved ranged from 500 to more than 5000 m 3 per event [START_REF] Kummert | Erosion and sediment transfer processes at the front of rapidly moving rock glaciers: Systematic observations with automatic cameras in the western Swiss Alps[END_REF]. The total mass wasting that supplies the rock glacier has been roughly estimated at between 5000 and 10 000 m 3 /yr, based on the volume of unstable rock zones and the evolution of the eroded front [START_REF] Delaloye | Blockgletscher und Hangrutschungen in Permafrostgebieten. Projekt "Mattertal[END_REF]. All recent debris flows occurred following intense snow melt in the catchment area upstream of the rock glacier, or following long-lasting or repeated rainfall. Their magnitude also depended on the debris available beneath the rock glacier's front.

Two deep-seated landslides that affect the site are worth noting. First, the so-called Breithorn landslidelocated above and to the right of the rock glacier, including the entire south-west facing slope -which is moving downwards at a rate of 0.4 m/yr [START_REF] Delaloye | Rapidly moving rock glaciers in Mattertal[END_REF]. The second deep-seated landslide, called Längenschnee, affects the left side of the site, including the south part of the front (CREALP 2015). Its displacement rate (from 1m/yr to 8m/yr) is about the same order of magnitude than the rock glacier's rate, so that it may still affect the cliffs surrounding the left side of Gugla rock glacier, and thus contribute to gravitational processes.

Eight boreholes and one geophysical campaign for seismic refraction profiles were conducted on the site in 2014 (CREALP 2015(CREALP , 2016) ) to investigate the composition and thickness of surface layers, and locate possible shear stress horizons. The data showed a layered structure of the rock glacier (Fig. 1c) with a 5 m thick macro-porous active layer composed of gneiss blocks. Deeper layers were found to be permanently at or below 0 °C throughout the year, and thus correspond to permafrost. An intermediate layer between 5 m and around 10 m deep, hereafter "ice-rich permafrost", is composed of coarse frozen materials embedded in an icy matrix containing ice lenses. Then, a layer composed of fine frozen sandy materials, termed "icepoor permafrost" hereafter, extends to the bedrock (located at about 20 m below the surface, with high variability). Shear horizon depths are highly variable (from 2 m to 23 m), but mainly located in the latter layer, and thus the rock glacier moves thanks to an overall creeping dynamic, rather than sliding on the bedrock.

A meteorological station (named GUG2) continuously records snow height, air surface temperature, and precipitation (Fig. 2). It is located 350 m to the North of the rock glacier. Surface displacements have been monitored by various GNSS sensors installed on the rock glacier along the decade. We selected one of these sensors (GUG1, see Fig. 2(a), red curve), set up next to our area of study since February 2012. In addition, three webcams provide hourly images showing different parts of the rock glacier front.

With the aim of investigating thermal behavior of the rock glacier at depth, from Autumn 2014 to Summer 2017, two thermistor chains (0.25 °C uncertainty, 0.065 °C resolution) continuously recorded temperatures at different depths (Fig. 3). Over this period, the results show that the annual temperature increase affects zones to a depth of around 10 m, whereas no thermal variations occur at deeper levels. The active layer thickness was estimated at 4.5 m (+/-20%), while deeper layers are permanently close to 0° C or less (CREALP 2016). From these data, we extracted some additional information relating to:

annual alterations to the active layer, by tracking the zero-curtain effect. During phase transition of water/ice, the temperature remains at 0 °C due to latent heat release (freezing) or consumption (thawing). In our case, we observed such zero-curtain effects from 1 m until 4 m, indicating a gradual freezing of the active layer during winter, and a complete thawing during summer (Fig.

3(b));

thermal diffusivity of the active layer, by evaluating the time-lag between air temperature at the surface and temperature variations at depth. If latent heat effects due to water phase changes and convection processes are neglected, the temperature is only determined by conduction, and its evolution at depth can be fitted by sinusoidal curves with a 1-year period. Thus, we can estimate the time-lag between changes to surface temperature and effects on ground temperature at different depths, and these data can be directly related to the thermal diffusivity coefficient by a theoretical relation. Finally, we obtained several values of the order of 1.2 ± 0.1 mm 2 /s, which confirms previously published rough estimations [START_REF] Hartlieb | Thermo-physical properties of selected hard rocks and their relation to microwave-assisted comminution[END_REF].

The seismological network monitors the lower part of the rock glacier, in the range from 2570 to 2720 m asl. It is composed of five seismic sensors (Sercel L22 geophones with a resonance frequency of 2 Hz) that have been recording data on the site since October 2015. Two of them (C2 and C4) are set up on the glacier's longitudinal axis, whereas the others are placed on the two stable sides. Thus, sensors C1-C3 are around 100 m apart, and form a cross-section in the middle part of the glacier, whereas sensors C4-C5 form a half cross-section 100 m downstream (Fig. 1(b)). All the seismic sensors have been set up on the surface of the rock glacier, and connected together to the digitizer (Nanometrics Centaur, sampling rate 200 Hz) with wires insulated by shealth. They are cemented on the top of relatively large, stable and flat boulders, ensuring a good signal coupling. A little tube in plastic shelters each of them, and permits to shield off any influence of rain, wind and snow.

Methods

-Ambient noise correlation

Since the last decade, ambient noise correlation techniques have allowed seismic wave propagation data to be extracted from the correlation of continuous noise recordings from two seismometers [START_REF] Campillo | Long-Range Correlations in the Diffuse Seismic Coda[END_REF]. The cross-correlation obtained is a reconstruction of the impulse response of the medium between the two sensors, and includes ballistic waves together with reflected ones [START_REF] Larose | Correlation of random wavefields: An interdisciplinary review[END_REF]) and the subsequent coda.

Raw seismic data was collected hourly on the Gugla rock glacier site. Pre-processing consisted of pruning and preparing these signals for cross-correlation, and comprised the following steps: (a) zero-mean fixing and detrending; (b) whitening (Fourier transform normalization between 1 and 30 Hz, with a Hanning window apodization) as a spectral normalization to minimize the influence of persistent monochromatic sources and to broaden the frequency range of measurements [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broadband surface wave dispersion measurements[END_REF]); (c) clipping (highamplitude removal by setting a maximum threshold equal to three-times the standard deviation of each trace) as a temporal normalization to eliminate spurious influence due to quakes and instrumental irregularities [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broadband surface wave dispersion measurements[END_REF]).

Subsequently, cross-correlation of signals from the two seismometers can be processed. For this study, only the data from sensors C1 and C2 were processed, because both provided good quality records (no gaps in time), and their locations allowed enough of the rock glacier's area to be investigated (Fig. 1). The results were normalized relative to the total energy (computed as the geometric mean of the squared signal traces), and finally time-averaged to produce daily correlograms. This normalization and averaging improved both the Signal-to-Noise ratio (SNR) and the stability of the signal obtained.

Since Rayleigh wave sensitivity is depth-dependent, seismic data were analyzed at different frequency bands to provide information on various depths within the medium investigated. In general, the probed depth decreases as the frequency increases. Therefore, the relative velocity changes (dV/V) and correlation coefficients measured produce slightly different patterns depending on the frequency filtering applied. For this study, seismic results were examined at frequency ranges between 1.5 and 14 Hz with a constant frequency bandwidth 𝑓 ()* -𝑓 (,-= 4 𝐻𝑧. This range was chosen to achieve a compromise between depth resolution (small bandwidth) and signal-to-noise (SNR) ratio in the correlograms (large bandwidth).

To compute relative wave velocity changes, a filtering step must be applied to increase the SNR ratio in the daily correlograms:

-For low frequencies (below 10 Hz), a Wiener filter was applied with K=L=8, which corresponds to an 8-day window and a 0.04-s time-lag [START_REF] Moreau | Improving ambient noise correlation functions with an SVD-based Wiener filter[END_REF]. The Wiener filter aims to minimize noise in a random process through a statistically based maximization of the coherence in several occurrences of the same signal corrupted by an additive and stationary noise.

-For high frequencies (above 10 Hz), raw correlograms are very noisy. As recommended by Moreau et al. [START_REF] Moreau | Improving ambient noise correlation functions with an SVD-based Wiener filter[END_REF], Singular Value Decomposition (SVD) and Wiener filters were combined.

SVD is commonly used to separate signal-to-noise subspaces, by retaining only the dominant singular values in the correlogram matrix. The choice between which singular vectors to retain or reject comes with a tradeoff between noise reduction and signal loss [START_REF] Moreau | Improving ambient noise correlation functions with an SVD-based Wiener filter[END_REF]. Here, the selection parameters were heuristically set to keep the 30 first singular values, and a Wiener filter was applied with K=10; L=7, corresponding to a 7-day window and a 0.05-s time-lag.

The time-window used to measure dV/V should start after the direct Rayleigh wave arrivals and should stop when the coda can no longer be distinguished from noise. Based on the distance between the two seismometers (around 100 m), the mean direct wave velocity (500 m.s -1 ) (Geo2X 2014), and after plotting the mean correlogram signal waveform, the time-window was set as t1=0.3 s to t2=0.8 s for both the causal and acausal parts of the correlogram (Fig. 4).

If the material undergoes a homogeneous velocity change 𝑉 → 𝑉 + 𝑑𝑉, the time axis for the last record 𝜑 7 (𝑡) will be stretched by a factor 𝜀 defined as follows:

𝜀 = - 𝑑𝑉 𝑉 = 𝑑𝑡 𝑡 (1)
The stretching method then corrects this effect by computing the correlation between the stretched signal 𝜑 7 (𝑡(1 + 𝜀)) and a reference signal 𝜑 < (𝑡) for several 𝜀 , values.

Without aliasing, the most accurate estimation of relative velocity change maximizes the correlation between these two signals, as explained below [START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia --2006[END_REF][START_REF] Hadziioannou | Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments[END_REF]:

𝐶𝐶 , = 𝐶𝐶(𝜀 , ) = ∫ 𝜑 7 (𝑡(1 + 𝜀 , )) 𝜑 < (𝑡)𝑑𝑡 ∫ ?𝜑 7 (𝑡(1 + 𝜀 , )) 7 𝜑 < (𝑡) 7 𝑑𝑡 (2) 𝐶𝐶 @ 𝑑𝑉 𝑉 A = max E F 𝐶𝐶 , (3)
The reference period is the period over which daily correlograms are averaged for comparison with the stretched signal for the day of interest. At each frequency band surveyed (between 1.5 and 14 Hz) the most stable 60-day period of the survey was taken as the reference. This corresponded to a frozen period from

February to March 2016.

The stretching approach is preferred to the cross-spectral moving window [START_REF] Clarke | Assessment of resolution and accuracy of the Moving Window Cross Spectral technique for monitoring crustal temporal variations using ambient seismic noise[END_REF]) because it provides direct access to the correlation value CC (the maximum correlation coefficient between the reference and stretched signals), together with the daily relative velocity change dV/V (the corresponding dilation value). CC is a measure of the similarity of the correlation response with respect to the reference, and thus reveals structural changes in the medium investigated that may be caused by several factors (water infiltration, rock deformation, cracks opening/closing, etc.). This value provides a direct estimate of decorrelation (Kd, discussed further below), which can be defined as:

𝐾𝑑 = 1 -𝐶𝐶 (
𝑅𝑀𝑆 (𝜀) = √1 -𝐶𝐶 7 2𝐶𝐶 Q 6S T 7 𝑇 𝜔 O 7 (𝑡 7 U -𝑡 < U )
(5)

This formula applies as long as the two correlated waveforms satisfy some assumptions (stationarity, noise-like and Gaussian, similar spectra, long time-window relative to the inverse of the central frequency, same amplitude), which was the case here.

-Detection, classification and localization of seismic events

The method described by [START_REF] Helmstetter | Seismic monitoring of Séchilienne rockslide (French Alps): Analysis of seismic signals and their correlation with rainfalls[END_REF] was used to extract seismic events from the continuous seismic records. Events were detected on the signal spectrogram, then averaged over all sensors (C1-C5). An event was defined as a SNR exceeding 3 in the 2-20 Hz frequency range.

Many different processes generate seismic signals. They may be natural (earthquakes, avalanches, debris flows, rockfalls, fracture opening, basal slip, wind, storms, etc.) or anthropogenic (trains, helicopters, hikers, etc.). For this study, the focus was on seismic events produced by the internal deformation of the rock glacier (basal slip, fracture propagation) and by erosion at the front of the rock glacier (debris detachment).

Therefore, two classes of seismic events were defined:

-"quakes", characterized by a short duration (less than 10 s), a low average frequency (less than 10 Hz) and a strong coherence between traces (greater than 0.5);

-"rockfalls", with a longer duration (longer than 10 s), a broad frequency range with a mean frequency of around 10 Hz, and larger amplitudes on sensors C4 and C5 located closer to the frontal zone. Signals were composed of successive peaks probably generated by successive rock impacts, similar to rockfall seismic signals detected in previous studies (e.g., [START_REF] Helmstetter | Seismic monitoring of Séchilienne rockslide (French Alps): Analysis of seismic signals and their correlation with rainfalls[END_REF].

Signals were automatically classified based on their duration, average frequency and inter-trace correlation.

Only limited time intervals (2015/10/1-2016/8/1) and the largest events (peak amplitude larger than 10 µm/s) were manually checked. The seismic signals, durations and average frequencies were very similar to other signals recorded on landslides at short distances (for a review, see [START_REF] Provost | Towards a standard typology of endogenous landslide seismic sources[END_REF]).

Figure 5 shows examples of seismic signals for each class of events. For quake events, the signals were very similar for all sensors, with only small time delays and variations in amplitude between sensors. In contrast, the rockfall signals lasted much longer, covered a broader frequency range, and generated more variability between sensors.

The mass-wasting event that occurred on June 6 th 2016 was one of the largest recorded. It can be seen on the webcam images taken a few hours before and after the event (Fig. 6). Comparison of the two images reveals that several mass-wasting events occurred during the night in the frontal zone of the rock glacier.

Seismic signals recorded on landslide or rock glaciers are often difficult to locate because P-and S-waves are difficult to separate, and because signals are emergent and arrival times picking lacks accuracy. The beam-forming method is a good alternative to manual picking when using a dense array of sensors [START_REF] Lacroix | Location of Seismic Signals Associated with Microearthquakes and Rockfalls on the Sechilienne Landslide, French Alps[END_REF]). This method maximizes the average correlation between sensors after shifting each trace to account for the source-sensor travel time. This travel time can be estimated as a function of the epicenter of the source and the apparent seismic wave velocity. A small apparent velocity indicates a shallow source and/or acoustic waves. In contrast, a large apparent velocity indicates a deep source, such as basal slip, deep fractures, or an external source (e.g. earthquake, distant avalanche or landslide). The number of sensors used for this study was too small to obtain an accurate source location; nevertheless, it was possible to estimate the azimuth of the source and the apparent velocity, providing useful information to identify the source process.

Results and discussion

1 -Ambient noise correlation 1.1 Temporal changes of dV/V

At very low frequencies ([1.5-5.5 Hz])

At the lowest frequencies (between 1.5 and 5.5 Hz), the seismic velocity determined from crosscorrelations between C1 and C2 showed moderate amplitude (+/-1%) changes without seasonal pattern (Fig. 7). This result suggests that the medium surveyed at these frequencies is located deeper than the zone affected by environmental seasonal fluctuations, such as temperature. The large error bars are probably a consequence of the low level of seismic signal (due to the lack of low frequency sources).

At low frequencies ([4-8 Hz] and [7-11 Hz])

As a general overview, 𝑑𝑉/𝑉 estimates for the 4-8 Hz and 7-11 Hz frequency ranges showed roughly the same pattern (Fig. 8). These 𝑑𝑉/𝑉 values are well-fitted by a sinusoidal curve with a 1-year period, with a higher seismic velocity in winter than in summer. In relation to the reference period (February to March 2016), all relative velocity changes oscillated between around +/-1.5% over the seasons, but the duration and intensity of increases or decreases varied depending on the year.

The slight decreasing linear trend observed in both cases corresponds to around 7. 10 XY relative velocity variations per year and may be due to the relative displacement of the sensors with the rock glacier.

Alternatively, it may indicate changes to the internal structure of the glacier over the study period. If we consider the annual displacement of the sensor located on the rock glacier (maximum 2 m along the longitudinal axis, since the sensor on the rock glacier side is fixed over the transverse axis), the maximal expected strain is 8. 10 XY . This value can be directly related to seismic velocity changes between the two sensors 𝜀 = -(Δ𝑣 𝑣 ⁄ ), and is consistent with the order of magnitude of the linear trend. Based on this result, the trend is probably not due to interannual structural changes. Moreover, this trend would not be very significant as we have only three years of measurement data. However, the trend should be closely monitored in the coming years, to identify higher velocity changes that could be effectively linked to interannual permafrost degradation. For the moment, these time series can be approximated in the first order by a sum of a linear trend and a 1-year periodic sinusoidal modulation. Similar variations are also observed for atmospheric thermal forcing and surface displacement rate time series. In the 4-8 Hz and 7-11 Hz frequency ranges, the 𝑑𝑉/𝑉 and air temperature time series recorded showed a significant anticorrelation (R~-0.4), with a lag time of 25 days. If we assume thermal forcing to be the main driver of such variations, this lag time correlates with thermal wave propagation with depth. This wave propagation is modulated by several effects inside the active layer, such as air convection and latent heat loss due to water phase transitions. These effects make the link between thermal indicators and depth-dependent seismic velocities more difficult to interpret than in a homogeneous medium, where thermal conduction dominates.

At high frequencies ([10-14 Hz])

At high frequencies, dV/V variations were less seasonally modulated than at lower frequencies. During all three melting periods covered by our data (red boxes on Fig. 9), we observed a sudden velocity decrease (around -1.5%), together with a simultaneous CC drop to below 0.6, corresponding to a decorrelation event (Fig. 9). These episodes were generally followed by rapid dV/V increases, that mostly occurred during wet summer periods. During late summer (September), we observed another apparently chaotic dV/V decrease (especially in 2017), which lasted until the first snowfalls in early winter.

At very high frequencies (above 14 Hz)

Above 14 Hz, filtered correlograms showed huge variations during summer periods, namely some peaks of velocity change of about 15% (not presented here). These peaks occurred only during summer, which is the period when the resonance frequency for the whole medium was around 20 Hz. These extreme values may be due to resonance phenomena at these frequencies [START_REF] Colombero | Integration of ambient seismic noise monitoring, displacement and meteorological measurements to infer the temperature-controlled long-term evolution of a complex prone-tofall cliff[END_REF], and are therefore difficult to compare with results at other frequency bands. The monitoring of resonance frequency and surface wave interpretation at frequencies above 14 Hz are thus out of the scope of this article, and such issues will be addressed in future studies.

Modelling

Seasonal changes to rigidity

At low frequencies (between around 4 Hz and 10 Hz), the clear seasonal pattern of velocity variations leads us to assume that dV/V is mainly sensitive to global variations in rigidity, which are due to freezing of the porous medium in the shallow part of the rock glacier.

To test this assumption, and to better constrain the physical processes, we build a seismic velocity model of the Gugla rock glacier (Fig. 10), that would be considered as the reference in summer. For this, we considered the values of Vp and Vs as a function of depth that were measured in summer. We then used the observed seasonal fluctuations of surface wave velocity changes (dV/V amplitudes) at different frequencies in order to estimate seasonal fluctuations of Vp and Vs. We could thus locate and quantify the seasonal changes in seismic velocity for P and S waves., that would account for the rigidity changes due to freezing processes causing dV/V variations.

At the first step of seismic modeling of the Gugla-Breithorn rock glacier, we collected seismic wave velocities (𝑉 ] (Fig. 10 (b)) and 𝑉 ^ (Fig. 10 (c)) for P-and S-waves, respectively) from geophysical prospections carried out in summer months, such as compressional and shear wave tomographies obtained by seismic refraction methods. These tomographies were obtained for one transversal profile of the rock glacier, and we selected values from a unique point on the profile, corresponding to the midpoint between sensors C1 and C2. Despite no quantitative density data from boreholes, we selected usual values from drilling campaigns in such frozen debris [START_REF] Haeberli | Core drilling through rock glacier permafrost[END_REF][START_REF] Wagner | Creep of alpine permafrost, investigated on the Murtel rock glacier[END_REF]) : the rock glacier density was assumed constant at 2000 kg/m 3 , whereas the bedrock density was fixed to 2800 kg/m 3 (Fig. 10 (d)).

By using Geopsy software 1 , a dispersion curve was computed from this input model (𝑉 ] , 𝑉 ^ and density along depth), corresponding to the theoretical Rayleigh wave phase velocities for the fundamental mode, at each frequency, over the summer period (Fig. 11(a), red curve).

(http://www.geopsy.org)

To check its validity, this dispersion curve was then inverted to obtain a seismic model (P-and S-wave velocity profiles) for summer. The inverted parameters are 𝑉 ] and 𝑉 ^ throughout the rock glacier depth, whereas both rock glacier density and bedrock parameters are kept fixed for the inversion.

As expected, the results of this inversion converge toward the input model, considered as the reference model (Fig. 11(b) and 11(c), red curves), which is quite similar (8% and 20% of uncertainties, for Vp and Vs respectively) to the geophysical data collected during summer. Thanks to this preliminary result, the accuracy of dispersion curve inversion method has been validated.

The observed dV/V maximum amplitudes _ `a a b ()*

were set as the amplitude of the fitted sinusoidal curve of seasonal variations, at each frequency band. These values could then be used to estimate the dispersion curve in winter, by applying the following formula:

@ 𝑑𝑉 𝑉 A ()* = 𝑉 c,-def -𝑉 ^g((ef 𝑉 ^g((ef (6)
This expression directly provides the analytic relation between the surface wave velocities in winter and in summer:

𝑉 c,-def = 𝑉 ^g((ef @1 + @ 𝑑𝑉 𝑉 A As expected, the results of the winter model showed an increase in compressional and shear wave velocities within the active and ice-rich permafrost layers, and allowed us to quantify these seasonal rigidity changes.

More precisely, P-wave velocity increased by 50% within the active layer and by 20% within the ice-rich permafrost layer (from the surface to a depth of around 10 m), whereas S-wave velocity increased by 40% within the active layer and by 25% within the ice-rich permafrost layer. These increases are in accordance with literature data for compressional wave velocity changes in gneiss following freezing [START_REF] Draebing | P-wave velocity changes in freezing hard lowporosity rocks: a laboratory-based time-average model[END_REF]. Considering the higher porosity in our case, together with field measurements from geophysical methods on other permafrost environment [START_REF] Kneisel | Advances in geophysical methods for permafrost investigations[END_REF], this evolution may be related to the partial freezing of the water-filled pores within the two layers (from the surface to around 10 m depth), that changes the ice-water ratio profile and consequently the overall rigidity.

However, these results of inversion must be considered with caution: although the S-wave model found for the winter appears robust, the corresponding P-wave model is not well constrained, due to the low sensitivity of Rayleigh waves to P-wave velocity. Hence the results presented here are one possible model in accordance with the expected physical processes responsible for the seasonal changes measured.

Thermal data

The following results were determined from thermal data acquired in a borehole set up from 2016 to 2018 near the seismic sensors (named F2, see location in Fig. 1). These data correspond to daily ground temperatures from the surface to 10.5 m depth, and show indications of freezing and temperature changes at depth. Thus, the active layer from the surface to 4.5 m depth was completely thawed in the summer, with temperatures considerably exceeding 0 °C, whereas deeper layers (from 4.5 m to 10.5 m depth) remained permanently below 0 °C, in line with the definition of permafrost (Fig. 3). The data also showed a zerocurtain effect from 1 m to 4.5 m depth, where the temperature remained at freezing-point during the solidliquid transition. Based on the depth at which it was observed, this zero-curtain effect indicates partial or complete freezing of the active layer at least. Deeper within the ice-rich permafrost layer, seasonal temperature variations were smaller but not negligible (between -5 °C in winter and -0.5 °C in summer), revealing an effect of thermal forcing as well. Indeed, the seasonal thermal wave propagates down to around 10 m depth, its amplitude decreasing with depth, in line with conductive models. This below freezing-point temperature cycling is responsible for the seismic velocity changes measured at these depths.

The seismological modeling was supported by these thermal data: both approaches converged toward a seasonal freezing cycle within the active and ice-rich permafrost layers down to around 10 m.

Role of water infiltration

At high frequencies ([10-14 Hz] and above), data presented in Figure 9 show that relative velocity changes (around -1.5%) often combine with simultaneous decorrelation of seismic responses (correlation coefficient CC drops from nearly 1 to less than 0.6) during melting periods. These observations suggest that dV/V and decorrelation at high frequencies are mainly driven by density and incompressibility changes due to water infiltration within the active layer. At these frequencies, seismic waves illuminate a shallower medium than at lower frequencies. In the shallower part of the Gugla rock glacier (a few meters deep, mainly composed of the active layer), heterogeneities, porosity and water infiltration are significant. Consequently, dV/V at these high frequencies is more sensitive to variations in liquid saturation, as a result of induced density changes, than to annual freeze-thawing of the medium as a whole.

In addition to changing noise sources, decorrelation (CC drops) may be caused by structural changes linked to scattering [START_REF] Planès | Decorrelation and phase-shift of coda waves induced by local changes: multiple scattering approach and numerical validation[END_REF]) that could be induced by cracks opening, rock movements, or water infiltration (water content increase) in the medium surveyed [START_REF] Théry | Tracking fluids in multiple scattering and highly porous materials : toward applications in non-destructive testing and seismic monitoring[END_REF]. As a result, fluid injection into a highly heterogeneous porous medium can lead to decorrelation events (measured through the variable 𝐾 `= 1 -𝐶𝐶, which is well correlated with water content). Coda Wave Decorrelation is thus an alternative measurement to track fluid levels inside rock glaciers.

In addition, Coda Wave Interferometry (dV/V measurements) simultaneously shows sensitivity to meltwater infiltration through the drop in seismic velocity noted in early summer months. At the end of winter, the snow cover on the rock glacier surface starts to melt and to percolate through the surface. We hypothesize that this process induces active-layer melting together with pore-filling by liquid water, altering the medium's density and incompressibility, and thus decreasing its velocity. Depending on the water residence time and preferential paths onto the porous materials, this melt-water may progressively drain downward, resulting in a dV/V re-increase. Snow melting from above the watershed together with heavy rainfall can occur throughout the following summer period, and thus promote refilling of the hydrogeological system.

We checked this water tracking hypothesis by modeling the rock glacier as a porous medium, and quantifying the influence of the moisture profile on seismic velocities through a poroelastic framework. In this view, the Biot-Gassmann theory describes the dependence of P-and S-wave velocity on fluid saturation changes, and can thus mimic the influence of moisture profile changes on Rayleigh waves [START_REF] Voisin | Seismic noise monitoring of the water table in a deep-seated, slow-moving landslide[END_REF].

Most poroelastic materials can be modeled as an effective medium with a single fluid phase and a single solid phase, described by homogenized parameters at the mesoscale [START_REF] Pride | Hydrogeophysics: Water Science and technology library[END_REF]. With this in mind, we modeled the entire rock glacier composed of the three layers indicated above (active, ice-rich and ice-poor permafrost layers, see Fig. 10 According to this model, the most important poroelastic attributes were the effective porosity 𝜙, together with solid parameters of the porous skeleton, that determine the effective mechanical moduli 𝐾 if (compressibility) and 𝐺 if (shear). A set of these parameters were adjusted to remain consistent with P-and S-wave velocity values for the reference model, obtained by seismic refraction campaigns (Geo2X 2014 p.

2). Porosity and water saturation levels were estimated from borehole data. The effective porosity of the active layer 𝜙 was fixed based on a consistent evolution with lithologic composition throughout the glacier's depth (Fig. 10(e)). Inside the active layer, we assumed a general decrease in porosity due to compaction, from a high value of 𝜙 = 0.35 for the two shallowest layers to 𝜙 = 0.15, as commonly used for permafrost. The porosity value for the basal layer was assumed to be slightly higher (𝜙 = 0.3), due to its structure consisting of a dense pore concentration within fine frozen materials.

Finally, we considered a water saturation profile derived from borehole data, obtained in July 2014 by the CREALP (CREALP 2016).

Given this set of poroelastic parameters for each sublayer of the rock glacier, following [START_REF] Voisin | Seismic noise monitoring of the water table in a deep-seated, slow-moving landslide[END_REF],

we computed the variations in Vp and Vs as a function of water saturation from 0 to 1 using the analytic expressions for both velocities (equations 9.19 and 9.15 in [START_REF] Pride | Hydrogeophysics: Water Science and technology library[END_REF] for Vp and Vs, respectively). The general result indicated that Vp increases with high saturation level values, whereas Vs decreases slightly.

These velocity variations are then the prerequisite result when modeling the effect of water filling on the Rayleigh wave velocity. We next modified the moisture depth (defined as the interface below which water saturation is equal to 1, and equal to the summer value above) in 0.5-m steps corresponding to each sublayer.

Whereas the moisture profile is undoubtedly more complex and smoother than a unit step function, the lack of moisture data means we have no choice but to approximate it in this simple manner. Assuming the linearity of seismic velocity changes with respect to depth, we investigated water filling for the three layers separately, in order to determine how the water influenced each of them.

By modifying the moisture depth (and thus the saturation profile), we obtained new seismic models of Vp and Vs (Fig. 12(a)). The respective dispersion curves were then computed using Geopsy software, and compared to the reference curve (corresponding to the reference model presented above, in Fig. 10). The differences between dispersion curves are directly related to Rayleigh wave velocity changes, and can be considered as theoretical dV/V induced by water filling, with respect to depth and frequency. Since the influence of water saturation in deep permafrost layers (below 5 m) is negligible (theoretical dV/V below 0.2 %) and not detectable from other effects, the final results are shown only for moisture changes in the active layer (0 to 4.5 m depth, Fig. 12(b)), where this effect is more significant and detectable by ambient noise correlation.

From these results, we noted the drop in dV/V (around -1.5%) at high frequencies (10 Hz and above), which resulted in the complete filling of the active layer, and matches the observed dV/V values during wet periods.

Using this poroelastic model, we can thus conclude that most of the seismic velocity changes observed at high frequencies can be explained by water infiltration due to snow-or ice-melt or precipitation. Assuming that these seismic velocity changes are only due to poro-elastic effects, the water infiltration is mostly located within the active layer, but water can also percolate to deeper layers (ice-rich and ice-poor permafrost layers), depending on preferential paths, crack openings and residence time for the water in the whole hydrological system. This quantification may be refined downwards, probably as a result of combined processes of percolation, drainage and partial refreezing of melting water in the active layer, which help buffer the influence of fresh water passing through the rock glacier.

-Seismic events

We detected 24 552 rockfall signals and 24 861 quakes between 2015/10/1 and 2019/1/1 with a peak ground motion velocity larger than 1 µm/s. The temporal evolution of the seismic activity is shown in Figure 13. Seasonal variations in the rate of events were more pronounced for rockfalls than for quakes, with a peak of activity in late spring and summer. The daily rate of quakes correlated well (𝑅~0.7) with the displacement rate (Fig. 2(a), red curve): an increase in the number of quakes was almost simultaneously followed by an increase in surface displacement. However, no clear evidence of such correlation occurred between the daily rate of rockfalls and other seasonal forcing (displacement rate nor temperature). The very high peaks of rockfall activity (more than 150 events per day) often occur during melting periods, and could thus be explained by the rebalancing of the mass at the front of the rock glacier in early summer. Indeed, during the winter months the freezing process limits erosion (very few rockfall signals), whereas the rate of displacement (Fig. 2(a), red curve) remains non-negligible. Hence, a mass excess appears at the front of the rock glacier, that is promptly eroded upon thawing, causing these observed peaks in the rate of mass-wasting events.

In addition, there are bursts of quake activity that last for a few days. Some of these bursts were associated with snowfall (indicated by green boxes in Fig. 13), whereas others correspond to rapid melting episodes (red boxes in Fig. 13). However, increased quake activity was not observed with every episode of snowfall.

The specific conditions involved in triggering quakes are further discussed below.

The increase in quake activity in summer can be explained by the increase in displacement rate and in pore water pressure. Indeed, crevasse propagation is promoted by an increase in water pressure, which allows

crevasses to remain open (van der Veen 1998).

The increase in quake activity during snowfall was more unexpected. Since most bursts of activity occurred during periods of very cold weather, they were not linked to melt or an increase in water pressure. Therefore, the only possible explanation for the observed effect is that the weight of the snow layer increased the stress on the underlying structure. The increase in stress may have been sufficient to induce a transition between stable (aseismic) and unstable (stick-slip) regimes, as expected from the rate-and-state friction law [START_REF] Dieterich | Modeling of rock friction 1[END_REF][START_REF] Scholz | Earthquakes and friction laws[END_REF]. Similar observations of snowfall-triggered seismicity have been reported at Mount-Rainier volcano [START_REF] Allstadt | Swarms of repeating stick-slip icequakes triggered by snow loading at Mount Rainier volcano[END_REF], where they were likely located at the base of the glacier covering the volcano. A more detailed analysis of these bursts will require further work, but preliminary results have been presented elsewhere [START_REF] Helmstetter | Repeating quakes triggered by snow-falls at Gugla rock-glacier: transition between stable slip and stickslip[END_REF].

Daily fluctuations in the rate of seismic events were also observed, as shown in Figure 14. Rockfalls occurred more often in the afternoon than at other times of day and correlated well with temperature.

Previous studies indicated that rockfall activity increases shortly after rainfall [START_REF] Helmstetter | Seismic monitoring of Séchilienne rockslide (French Alps): Analysis of seismic signals and their correlation with rainfalls[END_REF][START_REF] Delonca | Statistical correlation between meteorological and rockfall databases[END_REF]). In the case of the Gugla rock glacier, we suggest that melt-water production increases with temperature, thus enhancing mass-wasting activity due to the same processes by which rainfall can trigger rockfalls. This process could involve fracture opening due to increased water pressure, but the shallow layers of the rock glacier are likely too permeable to allow for an increase in water pressure.

Alternatively, unconsolidated blocks of rocks may be swept along by melt-water flowing down the surface of the rock glacier.

In contrast to rock fall, quakes were more frequent during the night than during the day. This pattern may be an artefact of the daily changes in seismic noise, as anthropogenic noise, for example, is higher during the day, thus decreasing the detection threshold. Indeed, when the minimum amplitude was increased from 1 to 5 µm/s, the daily fluctuations disappeared for quakes, but were still very clear for rockfalls.

The seismic events detected were difficult to accurately localize due to the small number of sensors and because signals are emergent, with no clear P-or S-waves. We selected all quakes and rockfall signals detected by at least five sensors between October 2015 and October 2018 for which the amplitude exceeded 1 µm/s. We applied the beamforming method developed by [START_REF] Lacroix | Location of Seismic Signals Associated with Microearthquakes and Rockfalls on the Sechilienne Landslide, French Alps[END_REF] to estimate the epicenter and the apparent velocity. Because it is impossible to estimate both the source depth and the seismic wave velocity, we assumed that the sources were located at the surface and considered a uniform seismic wave velocity. Signals were filtered between 3 and 20 Hz. We then applied the beam-forming method over a 4-s window centered on the peak amplitude. To remove erroneous locations, we further selected events with an average correlation between channels (after shifting signals in time to account for the travel time) exceeding 0.5. The set of events located included 6113 rockfalls and 12 080 quakes (Fig.

15(a)

). Many points were located along the boundary of the grid. The reason for this positioning is because events that are too far from the network cannot be accurately located using the beamforming method, and only the source azimuth can be well constrained (Fig. 15(b)). As expected, many rockfall seismic signals are located near the front of the rock glacier. This implies that most seismic signals classified as "rockfalls" are not real rockfalls (i.e., rocks falling from a cliff) but rather rock debris detachments from the front of the rock glacier. A few rockfalls are also located in the southern part of the rock glacier, near the front. This part of the rock glacier is surrounded by cliffs from which rocks occasionally fall, as recorded on the webcams. In addition, these cliffs are affected by the Längenschnee landslide, which may be partly responsible for the instability in this area. Quakes were mainly located East (upstream) of the seismic network, near the southern boundary of the rock glacier, an area where many open fractures are visible. The apparent velocity was larger for quakes (median 1805 m/s) than for rockfalls (1401 m/s) (Fig. 15(c)). This difference suggests that quakes are located at depth. The inclinometric data revealed the presence of several shear zones at depth between 3 and 23 m (CREALP 2015). These zones could produce some of the quake events. Quakes detected during snow-falls have quasi-periodic occurrence times and a narrow distribution of amplitudes, a pattern typical of stick-slip basal icequakes (Helmstetter et al., 2015b;[START_REF] Helmstetter | Repeating quakes triggered by snow-falls at Gugla rock-glacier: transition between stable slip and stickslip[END_REF]. Apart from these bursts, quakes occur randomly in time with a broad distribution of amplitudes, a behavior typical of icequakes associated with crevasse propagation (Helmstetter et al., 2015a). Because the seismic signals are emergent (smooth onset of the first waveforms), it was impossible to estimate the source process (slip or fracture opening) from the polarity of the P-wave arrival.

Conclusion

As a first passive seismic monitoring of a rock glacier in the field, this study confirm that ambient noise correlations could improve our understanding of permafrost degradation. In the long-term, ambient noise The seismic network was also used to detect seismic events generated by the rock glacier ("rockfalls" and "quakes"). The results revealed seasonal variations in the rate of events, correlated with annual changes in displacement velocity and in dV/V. In addition, rapid bursts of micro-seismicity occurred during strong melting episodes in the spring, simultaneous with rapid drops in dV/V. Additional work will be needed to better detect and characterize the source of these events.

Seismic monitoring could be used to detect mass wasting events and structural changes at depth in real time, with a view to creating an alarm system. Quantification of eroded volume from seismic signals could be useful to predict secondary processes, such as debris flows due to the accumulation of debris in the frontal zone. The template matching method could improve the detection and classification of seismic events (Helmstetter et al. 2015b). To accurately locate the source and identify the source process (fracture opening or basal slip?), temporary experiments using a dense network should be performed. Theoretical dV/V values predicted from dispersion curve differences between the reference curve (computed from the seismic model described in Figure 10 during a dry period) and curves computed for "wet models" with the changes to moisture depth indicated above. Over these frequencies, the more the active layer is filled with water, the slower the Rayleigh wave velocity, with the highest sensitivity recorded at high frequencies (10-14Hz)

Figure 13: Rate of seismic events for quakes (black) and rockfalls (red). Green boxes depict snowfalls, whereas red boxes indicate intense melting over the three years of data, associated with simultaneous bursts of rockfall activity. Blue arrows highlight peaks of quake activity associated with snowfall.

Figure 14: Daily temperature variations (black) and fluctuations in the rate of rockfalls (red) and quakes (blue). Circles and solid lines indicate the rate for events with a peak amplitude larger than 1 µm/s, whereas dashed lines with dots correspond to larger events with A>5 µm/s. The black curve represents the average hourly temperature in June (between 2015 and 2018), when rockfalls were most frequent.

4 )

 4 Uncertainties on the estimated dV/V values have to be addressed. With this in mind,[START_REF] Weaver | On the precision of noise correlation interferometry[END_REF] established a theoretical formulation of the statistical fluctuations of the stretching factor 𝜀 due to variations of the ambient noise field.[START_REF] Weaver | On the precision of noise correlation interferometry[END_REF] computed 𝜀 by comparing two signals that only differ from a statistical noise, that models the influence of any spurious measurement effects, like a change in source distribution. Due to this noise, an apparent stretching statistically occurs (𝜀 ≠ 0), and its root mean square (𝑅𝑀𝑆(𝜀)) allows to identify statistically significant variations of 𝜀 due to wave velocity change in the medium. This value has been confirmed by laboratory experiments and real seismic dataset, and depends on the correlation coefficient value 𝐶𝐶, the start and the end times, 𝑡 < and 𝑡 7 for the time-window processed, the inverse of the frequency bandwidth 𝑇 = 1 (𝑓 ()* -𝑓 (,-) ⁄ , and the central pulsation 𝜔 O , as below[START_REF] Weaver | On the precision of noise correlation interferometry[END_REF]):

  This modified dispersion curve takes velocity changes into account at each central frequency of the filtered band, and is interpolated between these points (Fig.11(a), blue curves). Inversion of this dispersion curve finally converges toward a new seismic model (a "winter model"), which provides a potential distribution of seismic velocities in winter throughout the glacier's depth (Fig. 11(b) and 11(c), blue curves).

  (a)), re-sampled as 0.5 m thick sublayers. A single P-wave velocity, S-wave velocity, porosity and saturation were attributed to each homogeneous sublayer based on observations and estimations (see Fig. 10(b) to (e)).

  correlation will provide new measurements to probe the changes occurring in rock glaciers. Coda Wave Interferometry allows accurate monitoring of rigidity and density within the whole medium, whereas Coda Wave Decorrelation provides an indicator of structural changes in surface layers, mainly due to snow melting and precipitations. Freeze-thawing effects are revealed by seasonal variations in seismic wave velocities, and are located between the surface and around 10 m depth. During melting periods, decorrelation and seismic velocity changes are indicative of water infiltration through the active layer (0 to 4.5 m). After further statistical analysis on available data and definition of threshold values above the observed fluctuations, these observables could potentially be used as precursor signals to predict overall destabilization threatening the downstream inhabited valley. Thermal and mechanical modeling of the Gugla rock glacier will further exploit the recorded seismic data.

FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: Meteorological data from the Gugla-Breithorn site during the seismic monitoring period. (a) Air temperature (black) and surface displacement rate (red) computed from GPS measurements. (b) Depth of snow cover (black) and the corresponding daily melting index (red), computed based on the decrease in snow depth over the preceding three days. Green areas indicate

Figure 4 :

 4 Figure 4: Wiener-filtered daily correlogram between sensors 1 and 2 during the survey period. Color bar indicates the normalized values of the daily correlation. Red boxes delimit causal and anti-causal parts of the correlogram used as input for the stretching method.

Figure 5 :

 5 Figure 5: Seismic signals (a,b) and spectrograms (c,d) for a rockfall (a,c) and a quake event (b,d). Seismograms high-pass filtered above 1 Hz are plotted for sensors 1-5. The peak amplitude A in µm/s is indicated for each sensor. The spectrogram (in log scale) is shown for sensor 1.

Figure 6 :

 6 Figure6: Webcam images taken on 2016/6/6 at 21:00 (a) and 2016/6/7 at 6:00 (b), a few hours before and after collapse of two regressive erosion sites (red circles). These events were detected by the seismic network (see seismic signals in Figure5).

Figure 7 :Figure 8 :

 78 Figure 7: Daily variations in surface wave velocity (dV/V in black, with error bars in gray) obtained from ambient noise correlation, normalized relative to a fixed reference period (February / March 2016) in the 1.5-5.5 Hz frequency band. For each dV/V value, the corresponding correlation coefficient is represented by the red curve. Green boxes indicate major snowfall whereas red boxes indicate the dates of significant melting periods. Blue lines highlight the dates of intense rainfall.

Figure 9 :Figure 10 :

 910 Figure 9: Same as Figure 7 over the frequency range 10-14 Hz.

Figure 11 :Figure 12 :

 1112 Figure 11: (a) Red: Reference dispersion curve from the reference seismic model in summer periods. Blue: interpolated dispersion curve from the summer data, integrating seasonal changes in amplitudes as a function of frequency of measurement.This new dispersion curve is thus related to winter periods. (b) Results of dispersion curve inversions, converging toward a model with the minimum number of misfit values. The figure shows the P-wave velocity profile (Vp) corresponding to the summer (red curve) and winter (blue curve) periods, as the results of inversion of the corresponding dispersion curves. (c)Same results for the S-wave velocity (Vs). The winter model shows slightly higher velocities than the summer model, with differences located in active and permafrost layers for Vp, and even deeper for Vs.

Figure 15 :

 15 Figure 15: (a) map of rockfalls (red dots) and quakes (blue dots). The green line delimits the contour of the rock glacier. The gray lines indicate the topography (one line every 10 m). Black triangles correspond to seismic sensors; (b) rose diagram of source azimuth relative to the center of the network (x=0, y=0) for rockfalls (red) and quakes (blue); (c) plot of apparent seismic velocity for rockfalls (red) and quakes (blue).

Acknowledgments

The GPS data, time-lapse photographs and the meteorological station are managed by ETH Zürich (J. Beutel, PermaSense project, http://data.permasense.ch) and Fribourg University (R. Delaloye, Geomorphology Research Group). Some valuable information about the geophysical campaign, data from boreholes and their interpretation were shared with permission from CREALP. We are particularly grateful to Benjamin Vial and Mickaël Langlais (ISTerre), Guillaume Favre-Bulle (CREALP) and Ludwig Haas (Wallis canton) and the geological department of Wallis for their invaluable assistance with fieldwork, site maintenance and seismic data retrieval. This work is supported by the OSUG@2020 Labex, the VOR-UGA program, the CNRS-INSU program, and the ANR LabCom GEO3iLab.