Guillaume Copie

Fabrizio Cleri
email: fabrizio.cleri@univ-lille1.fr

Younes Makoudi

Christophe Krzeminski

Maxime Berthe

Frédéric Cherioux

Frank Palmino

Bruno Grandidier

Surface-induced optimal packing of two-dimensional molecular networks

Keywords: numbers: 68, 43, Hn 07, 79, Cz 02, 70, Ns

High-density packing in organic crystals is usually associated with an increase of the coordination between molecules. Such a concept is not necessarily extended to two-dimensional molecular networks self-assembled on a solid surface, for which we demonstrate the key role of the surface in inducing the optimal packing. By a combination of scanning tunneling microscopy experiments and multi-scale computer simulations, we study the phase transition between two polymorphs. We find that, contrary to intuition, the structure with the lowest packing fraction corresponds to the highest molecular coordination number, due to the competition between surface and intermolecular forces. Having the lowest free energy, this structure spreads out as the most stable polymorph over a wide range of molecular concentrations.

Organic crystals in two or three dimensions (2D, 3D) can display polymorphism because of the competition among relatively weak forces between molecules. [1] Achieving controlled and reproducible crystal structures is of particular interest since polymorph crystals might have different physical properties. Variations in color, solubility, wettability, hardness have a strong impact on potential applications. On the other hand, avoiding undesired crystallization is a major concern in manufacturing thin films, drugs or organic dyes. [1][2][3][4] Among the thermodynamic forces affecting the formation of the crystal phase (temperature, solvent, additives and impurities, template), the activity (or concentration) of free molecules is one of the primary factors. For example, concentration was found essential to account for the preferential growth of α-glycine 3D crystals in comparison with other known polymorphs; [5] concentration was also the key parameter in the self-assembly of 2D molecular crystals with different geometries. [6,7] Coupled to the direct observation of molecular patterns by scanning tunnelling microscopy (STM), manipulating the concentration has proved to be an efficient route to selectively synthesize the polymorph of choice. [8,9] In this letter, we focus on the role of the surface in driving structural transitions of 2D molecular networks, self-assembled on a Si surface, as a function of molecule concentration. We used a combination of computer simulations and properly designed benchmark experiments to highlight the key role of the surface. Notably, all the above studies are consistent with the common wisdom that molecular packing should increase with molecule concentration. [START_REF] Kitaigorodskii | Organic Chemical Crystallography[END_REF] As a result, upon increasing concentration one should observe structures with the highest possible packing density, since higher coordination number implies higher contact density amongst the molecules, as verified for 3D organic crystals. [START_REF] Brock | [END_REF][12][13][14] However, such a straightforward consideration might be insufficient for 2D organic lattices assembled on a surface, because of the broken symmetry and the consequent role played by the molecule-surface interactions.

At increasing supersaturation, molecules start to cluster into structureless droplets on the surface. The loss of entropy, however, is soon outweighed by intermolecular and surface adhesion enthalpy, leading to the nucleation of ordered phases. [15] During this self-assembly process, the interaction with the substrate can contrast the intermolecular attraction, leading to the competition and coexistence of different ordered phases at a given value of activity (the gas pressure, translating into a given number n of molecules adsorbed in a finite surface patch). Since steric occupation competes with maximization of the coordination, the coordination of a low-density (LD) molecular network may be higher than for any other higher-density (HD) lattices, and such a LD network may become the most stable phase. In the following we give a startling example of lattices formed by planar molecules deposited on a boron-doped silicon surface for which, contrary to simple intuition, the predominant phase at high surface coverage is indeed a low-density and high-coordination one. 1,3,5-tri(4"-cyano-4,4'-biphenyl) benzene molecules (TCNBB), consisting of three identical arms each made of two phenyl rings, attached to a central phenyl at the 1, 3 and 5 carbons, and terminated by a CN group, were synthesized by a Suzuki cross-coupling reaction between adequate precursors. [16] The TCNBB molecules were subsequently evaporated at 270 • C in ultra high vacuum onto a clean and defect-free, boron-doped Si(111)

( √ 3 × √ 3)R30
• surface, indicated as Si:B(111) in the following. During the evaporation, the Si crystal was maintained at room temperature. Then the surface was imaged by STM at 77 K. The deposition of TCNBB led to two different networks (Figure 1). In the center of the STM image, a low density (porous) phase is observed, whereas on the outer margins of the image, molecules appear to form a denser, compact phase with a prominent feature of parallel rows. Individual molecules are usually resolved by probing the filled electronic states.

At such bias voltage, the ternary symmetry of the molecules clearly appears (Fig. 1(b),(c)).

The LD phase displays a honeycomb structure, each molecule connecting three adjacent pores, defined by the end-group of six molecules. In the high-density (HD) phase, molecules are arranged in alternating parallel rows. Note that all three arms of each molecule interact with one arm of a neighbour molecule in the LD, whereas this interaction only involves one arm per molecule in the HD phase.

When the surface is fully covered with TCNBB, large-scale STM images at random locations invariably show the LD network, Fig. 1(d), with only small patches of the HD coexisting in some images. By further increasing the surface coverage rate, we did not observe any additional phase transitions in the first monolayer. These experimental results indicate that the LD network is the most stable phase up to the highest coverage, at odds with recent studies of 2D self-assembled crystals. [8,9] In order to rationalize these observations, we performed multi-scale computer simulations by means of the theoretical framework we recently developed. [15,17] Firstly, densityfunctional theory and metadynamics (MtD) provide a best guess for the lowest-energy ad-sorption sites of single molecules and dimers on the surface. Then, molecular dynamics (MD) simulations using empirical force fields are performed with fully 2D-periodic molecular networks. The resulting stable structures are compared to the experimentally observed networks, finding an excellent agreement. In the case where several phases coexist, further MD simulations also allow to identify the kinetic pathway between the phases. In both the MtD and MD simulations, the MM3 force field is adopted, [18] including hydrogen bonding, and the improved description of Van der Waals forces already assessed in our previous study. [17] In practice, the infinitely-ranged VdW forces are rapidly decaying and are summed up to a cut off radius of 7 Å, large enough to avoid the need for a correction. Then, with the guidance of STM images, we built a LD and a HD network, using the dimer from MtD as the basic unit. The coverage rate is n LD = 0.33 and n HD = 0.37 molecules/nm 2 for the infinite (fully periodic in 2D) lattices. MD simulations were carried out on both the LD and HD networks. For both cases, the comparison of the calculated network structure against the experimental STM images is extremely good, as shown in In a simple interpretation, molecule-surface interactions may be taken proportional to n, while molecule-molecule interactions are rather proportional to the number of contacts, n 2 .

Then, an energy function of the type E(n) = αn + βn 2 should exist for any sublattice of given symmetry, at least for islands of not too large size compared to the range of molecular forces. If another phase competes with the first one, with different symmetry, coordination β ′ and surface interaction α ′ , the energy difference at zero temperature is ∆E(n) = (α -

α ′ )n + (β -β ′ )n 2 . At a critical activity n C = (α ′ -α)/(β -β ′
), a phase transition may be observed, depending on the value of the coefficients α, α ′ , β and β ′ . [15] Since the energy differences found in the present MD simulations (i.e., the coefficients (α -α ′ ) and (β -β ′ ))

have opposite sign, a phase transition can occur. We can therefore guess the following kinetic phase diagram, as a function of n: from n=0 up to n C denser patches of the HD phase form, initially stabilized by the predominant molecule-surface interaction; such local structures tend to disappear at increasing n > n C , when the patches turn into the LD phase, because the intermolecular interaction quickly gains over the surface term; this phase persists up to n LD , when further molecules cannot be added to the first monolayer without disrupting the symmetry.

Remarkably, the strongest intermolecular interactions correspond to the LD packing.

Such a sequence in 2D crystallisation is reminiscent of Ostwald's "rule of stages", [21,22] stating that the first crystal phase formed would be the least stable polymorph, followed by phases of increasing stability. Here, phase transformation kinetics provide a basis for this empirical rule. Note that turning a HD patch into a LD one actually increases the fraction of occupied surface. This also implies a slight increase in the fraction of broken bonds along the free perimeter of the patch, which however contribute only a rapidly vanishing n 1/2 surface tension term.

Since in the experiments some small patches of HD phase are sometimes observed to coexist with the predominant LD phase, we checked the possible existence of a kinetic barrier between the two structures, by using an elastic-bands MD method. Therefore, patches of molecules locally arranged in the HD network could be kinetically trapped, and would not spontaneously transform into the (thermodynamically favored) LD phase, even for n > n C , until a threshold temperature is reached Finally, we used KMC simulations [20] to characterize the finite-temperature kinetics and long-time stability of the LD phase, in the two limiting cases of (i) homogeneous nucleation from randomly adsorbed molecules, and (ii) heterogeneous nucleation from an isolated patch of HD phase.

To study the stability of the LD phase in the case (i) of homogeneous nucleation, we ran several KMC simulations at constant-{NVT}, with molecules randomly placed on the surface, at different activities and effective temperatures k B T /∆E ≈ 1. Up to a coverage of about 0.9 monolayer (the maximum occupancy for the LD phase, since n LD /n HD =0.9) KMC simulations always show the formation of a complete LD network, starting from homogeneously dispersed molecules. Upon increasing n the LD network cannot grow further, and nuclei of HD appear, however without growing at any temperature. These first KMC simulations are consistent with our previous results, which identify the LD phase as the thermodynamically stable state.

The case of heterogeneous nucleation (ii) was also simulated by KMC at constant-{µVT}, by including a model process of continuous deposition of molecules from the gas phase. We start the KMC simulation (Figure 3) with an isolated nucleus of HD network mimicking a n < n C situation, in which free molecules condensing from the gas in disordered droplets on the surface initially go into the closest energy minimum, i.e. the HD (see the above discussion about Ostwald's rule). Every N timesteps τ (or effective rate R = 1/N) we add one molecule on the surface at a random position around the border of the HD patch. Between two subsequent depositions, all the molecules can move either by jumps to a neighboring site, or turn on the site in steps of ±60 • . For R ≥ 2 × 10 -4 τ -1 , the HD patch grows, and new HD nuclei can also be observed for the fastest rate R = 2 × 10 -3 τ -1 . Conversely, upon slowing down R the HD phase starts to disappear; at the same time, the LD phase is nucleating and starts to grow. The effective surface diffusion coefficient is found to increase linearly with 1/R. Clearly, as R slows down, molecules diffuse away from the HD island and nucleate the LD patches. Such findings are qualitatively consistent with the experimental results at variable activity,[20] thus lending support to our kinetic explanation about the persistence of patches of the HD within the LD phase.

In conclusion, by combining experimental STM images with multi-scale computer simulations, ranging from DFT, to molecular dynamics and metadynamics, to kinetic Monte Carlo, we studied the competition between a low-density and high-density molecular network adsorbed on a passivated Si:B(111) surface. In agreement with the experiments, the low-density phase is the most stable one, thanks to a subtle balance between intermolecular and surface force contributions. Upon increasing the activity, the gain in molecular coordination switches the balance in favor of the former. Notably, an ordered lattice coupling lower density and higher coordination is quite a rare occurrence in 3D extended systems.

Our results suggest that, thanks to the surface interactions, this could be a much more common finding in 2D ordered systems. represents the number of steps between two successive landings. For the faster deposition rates (X=500) the HD network can grow (rectangular symmetry, grey shading). When slowing down the deposition rate, the HD network starts to disappear and the LD network starts to nucleate (X=5-10,000, hexagonal symmetry, red shading). At even slower rates (X ≥ 50,000), the LD network takes over the entire network.

[ 19 ]

 19 Finally, kinetic Monte Carlo (KMC) simulations on a frozen Si:B(111) surface, with energy parameters derived from MD, give insight into the long-time growth kinetics leading to the preferential self-assembly of the LD network.We performed several MtD simulations [20] by placing a molecule at random sites on the surface and relaxing the system energy by steepest-descent method. The minimum-energy adsorption site coincides with the Si adatoms, as shown in Figure2(a). Competing minima at nearly the same energy exist on the hollow sites (centers of three adatoms), separated from the adatom site by a potential barrier of about 5 kcal/mol. Next, we studied the free energy landscape when a second molecule is added on the surface. The first molecule (on the left in Fig.2(b)) is initially placed on a Si adatom; the second molecule is randomly placed within the interaction radius of the first one. By averaging over several MtD runs, the dimer configuration shown in Fig.2(b) is obtained, with minor variations from one run to the next. In this configuration, the extra intermolecular interactions force both molecules to shift, so that their centers occupy two adjacent hollow sites.Molecular flexibility plays an important role in stabilizing the observed networks, by allowing to optimize the molecular interactions. In Fig.2(a) the single molecule is flat on the surface, while addition of a second molecule (Fig.2(b)) involves a twist of the phenyl rings. By looking at the free energy differences, the hollow site in Fig.2(a) is at -5.5 kcal/mol, while it increases to -8.1 kcal/mol once the dimer is formed. Clearly, this difference comes entirely from the molecule-molecule interaction, including about 35% from H-bond formation and 65% from the VdW interaction, optimized by the twisting of adjacent arms.

Figure 2 (

 2 Figure 2(c) and 2(d). After quenching down to T=0 K, the two ordered structures differ in energy by 2.68 ± 0.23 kcal/mol. However, contrary to common wisdom, it is the LD network that yields the lowest energy. In fact, while the number of H-bonds per molecule is the same in the two lattices, the intermolecular interaction is larger for the LD network, in which all the three arms of a molecule are paired with another molecule, thereby saturating the intermolecular interactions. By breaking down the energy contributions, we find the VdW intermolecular energy and the electrostatic energy to be lower for the LD network by -5.1 ± 0.3 kcal/mol and -3.1 kcal/mol respectively, yielding an overall energy difference ∆E M M (LD-HD)= -8.2 ± 0.6 kcal/mol. Notably, the largest part of the total intermolecular energy is due to VdW forces, the average energy per molecule from H-bonds being only about 0.2 kcal/mol, equal for both networks. Conversely, by looking at the molecule-surface interactions, we find a lower energy for the HD lattice, equal to ∆E M S (LD-HD) = +5.3±3.6 kcal/mol Since at such a low temperature the total internal energy E(NV T )=E M M +E M S practically coincides with the free energy G(µV T ) at constant volume and chemical potential, the LD phase has indeed the lowest free energy, in agreement with the experimental STM images.

  [20] We find that both the intermolecular (MM) energy and the molecule-surface (MS) energy display a migration barrier in going from the HD to the LD, of E M = 1.5 kcal/mol for the MM and E M = 4 kcal/mol for the MS interaction, both quite larger than room-temperature fluctuations.[20]

FIGURESFIG. 1 .

 1 FIGURES


Agnès Helmstetter

Luc Moreau

Barbara Nicolas

Pierre Comon

Michel Gay

Intermediate-depth icequakes and harmonic tremor in an Alpine glacier (Glacier d'Argentière, France) : evidence for hydraulic fracturing?

We detected several thousand deep englacial icequakes on Glacier d'Argentière (Mont-Blanc Massif) between 30 March and 3 May 2012. These events have been classified in eight clusters. Inside each cluster, the waveforms are similar for P-waves and S-waves, although the time delay between the P-waves and S-waves vary by up to 0.03 s, indicating an extended source area. Although these events were recorded by a single accelerometer, they were roughly located using a polarization analysis. The deepest events were located at a depth of 130 m, 60 m above the ice/bed interface. The clusters are separated in space. The largest cluster extends over about 100 m. For this cluster, the strike of the rupture plane is nearly parallel to the direction of the open crevasses, and the dip angle is 56 degrees. Deep icequakes occur in bursts of activity that last for a few hours and are separated by quiet periods. Many events occurred on 28 and 29 April 2012, during the warmest days, when snowmelting was likely important. The distributions of interevent times and peak amplitudes obey power laws, as also observed for earthquakes, but with larger exponents. The polarity of the P-waves for all of the events is consistent with tensile faulting. Finally, between 25 April and 3 May, we observed a gliding harmonic tremor with a fundamental resonance frequency that varied between 30 Hz and 38 Hz, with additional higher frequency harmonics. During this time we also observed shallow hybrid events with high-frequency onsets and a monochromatic coda. These events might be produced by the propagation of fractures and the subsequent flow of water into the fracture. The strongest resonance was observed just after a strong burst of deep icequakes, and during an unusually warm period when the snow height decreased by 60 cm in one week. The resonance frequency shows a succession of several sharp decreases and phases of progressive increases. One of the strongest negative steps of the resonance frequency on 28 April coincides with a burst of deep icequakes. These events appear to be associated with the propagation of fractures, which can explain the decrease in the resonance frequency. Finally, we observed an acceleration of glacier flow on 29 April, suggesting that meltwater had reached the ice/bed interface. These observations suggest that deep icequakes are due to hydraulic fracturing, and that they can be used to track fluid flow inside glaciers.

Introduction

Alpine glaciers produce a wide variety of seismic signals. The most frequent seismic events are shallow icequakes that are associated with fracture opening or propagation [START_REF] Neave | Icequakes on the Athabasca Glacier[END_REF], Mikesell, et al.(2012)] and to the collapse of large ice blocks [O' Neel et al.(2007), [START_REF] Roux | Microseismic activity within a sérac zone in an Alpine glacier (glacier d'Argentière, Mont Blanc, France)[END_REF]]. Basal icequakes are much less frequent. These can be produced by the opening or closing of tensile faults near the ice/bed interface due to changes in the basal water pressure [START_REF] Walter | Evidence for near-horizontal tensile faulting at the base of Gornergletscher, a swiss Alpine glacier[END_REF], [START_REF] Walter | Deep icequakes: What happens at the base of Alpine glaciers?[END_REF]], or by basal slip [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]]. Seismic events have also been detected at intermediate depths that are below the usual depth of open crevasses, but are above the ice/bed interface [START_REF] Deichmann | Evidence for deep icequakes in an Alpine glacier[END_REF], [START_REF] Walter | Basal icequakes during changing subglacial water pressures beneath Gornergletscher, Switzerland[END_REF], [START_REF] Walter | Moment tensor inversion of icequakes on Gornergletscher, Switzerland[END_REF], [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF]]. These events might be due to hydraulic fracturing. Given sufficient water supply, a water-filled fracture can propagate to the base of a glacier when the overburden stress at the fracture tip is offset by the density contrast between the ice and water [ Van der Veen(1998), Van der Veen(2007), [START_REF] Carmichael | Seismic multiplet response triggered by melt at Blood Falls, Taylor Glacier, Antarctica[END_REF]]. Under certain circumstances, bottom crevasses might also propagate upward [ Van der Veen(1998)].

Finally, another class of seismic signals detected on glaciers is produced by the resonance in subglacial water-filled cracks and conduits. This class includes low-frequency events [START_REF] St | Hydraulic transients: A seismic source in volcanoes and glaciers[END_REF], Métaxian, et al.(2003)], hybrid events [START_REF] West | Glacier microseismicity[END_REF]], and tremors [START_REF] Winberry | Seismic observations of transient subglacial water-flow beneath MacAyeal Ice Stream, West Antarctica[END_REF], [START_REF] Heeszel | Humming glaciers[END_REF], [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF]]. Hybrid events are characterized by high-frequency onset followed by lower-frequency coda, similar to signals detected on volcanoes [START_REF] West | Glacier microseismicity[END_REF], [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF]]. They are interpreted as due to the opening of englacial fractures and the subsequent flow of water into the fractures. Subglacial water flow can also produce harmonic seismic tremors that can last for several hours [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF], [START_REF] Heeszel | Humming glaciers[END_REF]]. However, harmonic tremor can also be produced by the repeated failure of small subglacial asperities [START_REF] Winberry | Nucleation and seismic tremor associated with the glacial earthquakes of Whillans Ice Stream, Antarctica[END_REF]].

In this report, we have studied intermediate-depth icequakes detected at Glacier d'Argentière (Mont-Blanc massif) in the spring of 2012. Most seismic signals are likely associated with shallow sources, as either serac falls or crevassing [START_REF] Neave | Icequakes on the Athabasca Glacier[END_REF], [START_REF] Deichmann | Evidence for deep icequakes in an Alpine glacier[END_REF], [START_REF] Roux | Microseismic activity within a sérac zone in an Alpine glacier (glacier d'Argentière, Mont Blanc, France)[END_REF], [START_REF] Walter | Basal icequakes during changing subglacial water pressures beneath Gornergletscher, Switzerland[END_REF]]. We have also identified two classes of deep icequakes that are characterized by higher frequencies, shorter durations, and distinct P-waves and S-waves. In the first part of this two-part submission [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]], we have analyzed the first class of deep events, and we suggest that these events were associated with repeated ruptures of asperities that were located at the ice/bed interface. Here, we focus on the second class of deep events, which are located at intermediate depths that are deeper than the depth of open crevasses (about 20 m) but are above the bedrock. We describe the distribution of these events in time, space and magnitude, and discuss the possible mechanisms that are responsible for these events. We also analyze the spectral content of the seismic noise, which shows a resonance at a frequency of about 35 Hz. Several methods used in this study (e.g., detection and location of seismic sources) were also used for the analysis of the basal icequakes. More details on these methods are given in [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]].

Study area and instrumentation

Glacier d'Argentière is a temperate glacier in the Mont Blanc massif in the French Alps. This glacier has a length of about 10 km, from its top part at an altitude of 3400 m a.s.l., to its tongue at an altitude of 1600 m a.s.l.. We installed an accelerometer in the ablation zone, at an elevation of 2380 m, about 600 m upstream from Lognan icefall. At this location, the glacier moves on average by 18 cm per day, and the maximum ice thickness is approximately 200 m [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on Glacier d'Argentière[END_REF]]. Figure 1 shows an aerial photograph of the glacier, with the location of the accelerometer. The seismometer was in operation during five time windows that covered 44 days between February and June 2012: 17-25 February, 12-17 March, 30 March to 12 April, 25 April to 3 May, and 5-19 June. The gaps in the monitoring were due to difficult glacier access, insufficient power supply and limited data storage capacity. The accelerometer was a MEMS-based Sercel DSU3 that recorded in continuous mode at a sampling rate of 1000 Hz. Due to the high frequency content of the deep icequake signals, the data was undersampled. In future studies, the sampling rate should be increased to 2000 Hz, to avoid potential aliasing effects and to improve location accuracy. The sensor was installed in winter by digging 2 meters into the snow layer down to the snow/ice interface. On 9 May, there was still 2.5 m of snow above the sensor. Therefore, we are confident that the orientation of the sensor remained correct throughout the time interval of 30 March to 3 May that was analyzed in this study.

A permanent meteorological station is installed near the seismic sensor (Figure 1, circle). This station provides temperature, but has no measure of precipitation. The meteorological data are available from the GLACIOCLIM database [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on Glacier d'Argentière[END_REF]]. We also used data from two stations from Météo-France. One of these stations is installed 4 km north of the network (Le Tour village) at an elevation of 1500 m a.s.l., and this was used to obtain the precipitation data. The other station is located in the Aiguilles Rouges massif, at approximately the same elevation (2365 m a.s.l.) as our network and at about 7 km to the west; this station provided the snow height measurements.

Under Glacier d'Argentière, the water-collection tunnels constructed by Emosson S.A. allow access to several natural and permanent subglacial cavities beneath the Lognan ice fall, at an elevation of 2180 m a.s.l. The glacier basal sliding velocity has been measured since 1970 in one of these cavities, using a cavitometer [START_REF] Vivian | Subglacial cavitation phenomena under the glacier d'Argentière, Mont Blanc[END_REF], [START_REF] Bocquet | Measurement of ice movement in subglacial cavities: a new cavitometer beneath the Glacier d'Argentière (Mt Blanc, France)[END_REF], Moreau(1999)]. This instrument is made of a wheel that is kept in contact with the ice by a spring that acts on an articulated arm.

Detection and classification

We applied an automatic detection method based on the spectrogram of the seismic signal to obtain a catalog of seismic events [START_REF] Helmstetter | Seismic monitoring of Séchilienne Rockslide (French Alps) : analysis of seismic signals and their correlation with rainfalls[END_REF]]. This method provides a list of events with their characteristics (e.g., peak amplitude, duration, dominant frequency). We then looked manually at all of the events to classify them as shallow icequakes, serac fall, or deep icequakes. Examples of waveforms for these different types of events are shown in Figure 2. The signals in this Figure 2 were converted into ground velocity to allow for comparisons with signals presented in previous studies. Unlike shallow icequakes, the deep events are characterized by very short durations (< 0.1 s), an impulsive onset, distinct P-waves and S-waves, and the absence of surface waves [START_REF] Deichmann | Evidence for deep icequakes in an Alpine glacier[END_REF], [START_REF] Walter | Basal icequakes during changing subglacial water pressures beneath Gornergletscher, Switzerland[END_REF]]. The deep events were further classified as either basal or intermediate-depth icequakes, based on the source depth. The basal events were analyzed in [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]], while the present study focuses on the intermediate-depth events. These intermediate-depth events were located at least 60 m above the ice/bed interface, but generally below the surface crevassing zone. These events are referred to as "deep icequakes" in the following.

By looking at the seismograms of these deep icequakes, they could be classified into several clusters, with very similar waveforms within each cluster. The time delay between the P-waves and S-waves varies by up to 0.03 s within one cluster, corresponding to a change in the source-sensor distance of 87 m. To detect events that belonged to each cluster, we performed cross-correlation of a template event with the full continuous database. The template event was first chosen randomly by picking one event in our list of deep events. As the time delay between the P-waves and S-waves varies significantly within each cluster, we separately cross-correlated the signals for the P-waves and S-waves and allowed for a variation of 50% in the time delay between the P-waves and S-waves relative to the template event. We fixed a minimum correlation of 0.6 to classify signals into clusters. We then inspected visually all of the signals, and rejected those events that were significantly different by eye from the template event. We sometimes found waveforms that were different from the template, but were nonetheless typical of deep icequakes. In this case, these waveforms were associated with a new cluster, and were used as new template events. When an event was detected by several templates, we chose the template with the highest correlation. We then updated the template signal by stacking all of the events in each cluster. Before stacking the signals, we aligned the signals in time, and normalized each signal according to its peak amplitude. We cut each signal into two parts between the P-wave and S-wave arrival times. Each part of the signal for the P-waves and S-waves was aligned separately, by cross-correlation with the template signal.

We identified eight clusters of deep icequakes, which are labeled as A-H, and these contained between 220 and 5578 events. The waveforms of the template event of each cluster are shown in Figure 3. For each cluster and each component, the gray shaded areas in Figure 3 represent the standard deviations of the signals for all of events in the cluster (after normalization and alignment). All of the signals show impulsive P-waves and S-waves that were separated by a few hundredths of a second. The polarity and the relative amplitude of the P-waves and S-waves for the three components vary from one type of events to another. Types A and B differ only by the presence of multiple P-wave arrivals for the second type, while the waveforms for cluster A are simpler. Similarly, the waveforms for clusters C and G are similar to clusters D and H, respectively, but the signals for clusters C and G are simpler. The events of clusters B, D and H occurred about one month later than those for clusters A, C and G. For all types of events, the polarity of the P-wave arrival on the vertical component is upward. This is compatible with tensile failure, as suggested by [START_REF] Walter | Moment tensor inversion of icequakes on Gornergletscher, Switzerland[END_REF]] from a moment tensor inversion of deep icequakes on Gornergletscher.

Most of the events that were identified occurred between 30 March and 3 May, 2012, with only nine deep icequakes detected earlier. We thus focused on this period in this study. Note that we might have missed events of different types. We did not look systematically at all of the events. We only identified eight different types of deep icequakes, and we searched for similar events in the full continuous database by cross-correlation with a template.

Location of seismic sources

We can roughly locate the source of these icequakes even with a single sensor. First, the distance between the source and the sensor can be estimated from the time delay between the P-waves and the S-waves

d = t S -t P 1/V S -1/V P , (1)
where t P and t S are the arrival times of the P-waves and S-waves, respectively. We assumed a uniform seismic-wave velocity of 3600 m/s for the P-waves, and 1610 m/s for the S-waves. These values were measured by [Vallon(1967)] for Glacier du Tacul, which is a temperate glacier in the Mont Blanc massif that is at about the same altitude as our seismometer. To estimate more precisely the time delay t S -t P between the P-waves and S-waves, we fit each signal of a given cluster using the following signal, as in [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]]:

y t (t) = y t,P (t -∆t P ) + y t,S (t -∆t S ) (2)
where y t,p and y t,s are the parts of the template event that correspond to the P-waves and S-waves, respectively, and ∆t P and ∆t S are the time delays, respectively. For each signal, we searched for the values of ∆t P and ∆t S that maximize the correlation between this signal and the template event y t . The inversion started with a grid search with maximum time shifts of 1 ms for the P-waves and 20 ms for the S-waves. Then the best values of ∆t P and ∆t S were chosen as the starting values of a simplex optimization [START_REF] Nelder | A simplex method for function minimization[END_REF]]. We obtained a relative accuracy for the time delays on the order of 0.1 ms, about one tenth of the sampling rate. However, the absolute accuracy was much larger, as it was on the order of 1 ms. In addition, the direction of the source can be estimated by analysis of the polarization of the particle motion [Vidale(1986)]. We selected the time window associated with the P-waves, from the picked arrival time t P to t P + 0.006 s. We then applied a singular value decomposition to the analytic signal, without integration or filtering. The azimuth α of the source and the apparent angle of incidence φ were estimated from the direction of the first eigenvector. The apparent angle of incidence of the P-waves is, however, different from the direction of propagation, because it is affected by reflections at the free surface. We thus used the method of [START_REF] Neuberg | Effects of volcano topography on seismic broad-band waveforms[END_REF]] to correct for this effect. Assuming a planar free surface, the true angle of incidence φ c can be estimated from the apparent angle φ by:

φ c = sin -1 V P V S sin φ 2 . (3)
Before applying this correction, we subtracted the surface slope along direction α from the apparent incidence angle. The surface slope around the seismometer reached 8 degrees along the glacier flow direction. Equation (3) has no real solution for φ > φ * with φ * = 2 sin -1 V S V P = 53 degrees. The source depth z is given by z = d cos(φ c ) and it varied between 0 and about 130 m. For incidence angles φ ≈ φ * , small errors on φ yield large uncertainties on φ c , and therefore on the depth. We found 13% of events with an estimated incidence angle φ > φ * , which should not be observed due to free surface effects. For these events, we could not estimate φ c nor the source depth. Values of φ > φ * might be due to uncertainties in the estimated apparent incidence angle or to diffracted waves on the fractures. An example of waveforms and particle motion is shown in Figure 4. The polarization of the P-waves and S-waves is not perfectly linear, as would be expected for body waves. This might indicate heterogeneity, surface effects, or the influence of reflected waves at the ice/bed interface.

The polarization analysis and the correction (3) were validated with numerical simulations, earthquakes with known location, and basal icequakes detected by two 3-component stations [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]]. In this study, we performed additional tests on synthetic seismograms to determine the influence of undersampling. We generated 1000 synthetic seismograms with uniform distributions of source depth between 20 m and 150 m, epicentral distances between 0 m and 100 m, assuming tensile faulting on planes with dip between 0 and 90 degrees, and strike between 0 degrees and 360 degrees. We assumed a point source with a rise time of 1 ms (i.e., duration of displacement) and an attenuation factor Q = 20 for both P-waves and S-waves [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]]. We computed the ground acceleration with a sample rate of 3000 Hz, and then undersampled the signal down to 1000 Hz. We applied the same polarization analysis as for the real signals of deep icequakes, using a sample length of 6 ms. We also assumed a standard error of 1 ms on t P and t P . We obtained a standard error of 12 degrees for the azimuth and 2 degrees for the incidence angle. The median absolute errors are much smaller, because the error distributions are not Gaussian, with median errors of 0.4 degrees for the azimuth and 0.2 degrees for the incidence angle. Surprisingly, the vertical location error (standard error, 3.2 m) is smaller than the horizontal error (standard error, 4.3 m). The vertical location error is larger for the shallowest events, because the incidence angle is close to φ * , and small variations of the incidence angle when φ ≈ φ * yield large variations of the corrected incidence angle φ c estimated from Equation (3).

The location uncertainty is, however, much larger for real signals, and is difficult to evaluate due to poor signal-to-noise ratio, coarse sample rate, and heterogeneity of the seismic wave velocities. One way to evaluate the accuracy of the location is to analyze the distance between successive events. Events that occur very close in time are likely to also be very close in space. We selected all of the events with an amplitude A > 200 µm/s 2 and a correlation with the template C > 0.7. For each event, we measured the distance between this event and the next event of the same cluster, and obtained a median distance of 15 m. However, this method overestimates the location accuracy, because successive events are likely in a similar location, with a similar source mechanism and source-receiver path. [START_REF] Helmstetter | Basal icequakes recorded beneath an alpine glacier (Glacier d'Argentière, Mont Blanc, France) : evidence for stickslip motion? in press[END_REF]] applied the same location method to basal icequakes, and found good agreement between the icequake depths and the glacier thickness, with a mean difference of 2 m and a standard error of 18 m. This further suggests that the location accuracy (one standard error) is on the order of 18 m, although the depth of the shallow events is more uncertain when φ ≈ φ * .

A map of the icequakes is shown in Figure 5, for all of the events with A >300 µm/s 2 , C ≥ 0.7 and φ < φ * . Each cluster has a distinct location and extends over several tens of meters. Most events are deeper than the maximum depth of open crevasses (about 20 m below the surface [Paterson(1994)]), but shallower than the ice/bed interface. Clusters A, C and G, which occurred before 13 April, were much shallower than the other clusters that were only active after 25 April. The average depth was 30 m before 25 April, and 70 m after this date. The deepest events at ≈ 130 m occurred in the central part of the glacier, where the ice thickness was approximately ≈ 190 m [START_REF] Vincent | Glacier thickening and decay analysis from 50 years of glaciological observations performed on Glacier d'Argentière[END_REF]]. This implies that the fractures that generated these icequakes did not reach the glacier bottom, or that the fracture propagation below 130 m depth was aseismic.

The depth of the icequakes for cluster B varied with time. Figure 6 shows a complex pattern with successive increases and decreases in the source depth with time. Before 25 April, the depth of the icequakes ranged between 20 m and 50 m, with a mean of 30 m. When the temperature rose after 25 April, up to 9 degrees C (Figure 7c), the average depth increased up to 70 m, and the maximum icequake depth reached 130 m. In the morning of 28 April, the depth increased gradually from 60 m to 85 m in 1.4 hour, at a mean rate of 18 m/hr. During the last days of the acquisition after 30 April, the temperature decreased below freezing and the mean icequake depth decreased. These observations suggest that higher temperatures increased meltwater production, which raised the water level in fractures and promoted propagation of the fractures at depth.

There is likely a smooth transition between the shallow and deep seismic activities. For the shallowest events detected, the P-waves and S-waves overlapped. Therefore, closer events cannot be detected in the first step of our detection procedure, when separately correlating the continuous signal with the P-waves and S-waves of the template signal.

To analyze the spatial structure of the largest cluster, B, we used principal component analysis of the three-dimensional cloud of points [START_REF] Michelini | Application of the principal parameters method to the 1983 Coalinga, California, aftershock sequence[END_REF]]. We estimated the plane that best matched cluster B from the first two eigenvectors. We then computed the strike direction and the dip angle of this plane. We found a strike of 67 degrees, which is nearly parallel to the direction of the open crevasses near this cluster, a dip angle of 56 degrees, and a planarity coefficient of 0.7. This coefficient is defined by 1 -2λ 3 /(λ 1 + λ2), where λ 1 > λ 2 > λ 3 are the three eigenvalues, and it ranges between 0 (no structure in space) and 1 (all points belong to the same plane).

The polarity of the P-waves for all of the clusters and all of the components are consistent with tensile failure (positive in the direction away from the source), as also observed by [START_REF] Carmichael | Seismic multiplet response triggered by melt at Blood Falls, Taylor Glacier, Antarctica[END_REF]]. On the other hand, shear failures produce mixed first-motion polarities. A single-station and single-source polarity measurements thus cannot exclude shear sources. However, as all of the eight clusters with different source azimuths have compressive P-wave polarities, this argues against a shear source mechanism.

The orientation of the fractures at depth can also be observed from boreholes. [Vallon(1981)] drilled a borehole into Glacier d'Argentière in 1980, at an elevation close to our network (a few hundred meters above Lognan ice fall), but closer to the side, at about 200 m from the right side of the glacier. He observed two families of fracture planes with a dip angle that decreased with depth, from about 80 degrees at 40 m down to about 0 degrees at the ice/bed interface at a depth of 240 m. The first family of fractures corresponded at the surface to marginal crevasses with an angle of about 45 degrees from the margin. These fractures had a right-lateral motion that was consistent with the glacier motion being slower near the sides than at the center of the glacier. [START_REF] Fountain | Fractures as the main pathways of water flow in temperate glaciers[END_REF]] also observed fractures at depth in boreholes for a small Swedish glacier. They suggested that steeply dipping fractures with openings on the order of several centimeters provide the main pathways of water flow in temperate glaciers. These features can be observed from near the surface down to a few meters above the glacier bottom. Their strikes are generally parallel to the direction of the surface crevasses.

Distribution of events in time and size

Figure 7 illustrates the temporal variations of the seismic activity, meteorological parameters, snow height, and basal displacement. The seismic activity was very intermittent, with quiet periods that were separated by bursts of activity with several hundreds of events per hour. Each burst usually involved only one cluster, and lasted for a few hours. In contrast, the rate of shallow icequakes (Figure 7b, gray line) was much less sporadic than deep events. Shallow icequakes were more frequent after 27 April, although the fluctuations of the rate of shallow events were not clearly related to the rate of deep events. During bursts of deep icequakes, there were more deep events than shallow ones. However, when looking at the whole duration of the study, there were many more shallow events than deep icequakes.

During the first acquisition period, which was between 30 March and 12 April, the temperature was beneath freezing most of the time (Figure 7c), so that meltwater production was very small. The highest peak of activity occurred on 4 April at noon (14:00 local time), with 507 events in one hour. Before 12 April, there is no apparent correlation between the icequakes and the precipitation. Strong precipitation occurred between 21 and 25 April, with a total of 40 mm of water that yielded 57 cm of snow at the meteorological station of Aiguilles Rouges (elevation of 2365 m, about 7 km to the west of the study area). Then between 25 and 27 April, the mean daily temperatures rose from -7 degrees C up to 9 degrees C. During this warming episode, the snow height decreased by 60 cm (snowmelt and/or compaction). This strong increase in the temperature was followed on 28 April by a strong peak of the rate of deep icequakes, and by acceleration of the glacier flow. During the following days, both the temperature and the rate of deep icequakes decreased on average. Two small bursts of deep icequakes on the evenings of 30 April and 1 May coincided with precipitation.

[ [START_REF] Carmichael | Seismic multiplet response triggered by melt at Blood Falls, Taylor Glacier, Antarctica[END_REF]] detected icequakes at Taylor Glacier, Antarctica, which were clearly asso-ciated with episodes of surface melt. Surface melt events produced an increase in the background rate of icequakes. These icequakes were located near surface cracks, but their depths were not constrained.

The slow variation of icequake activity observed in this study is more similar to our observation of the rate of shallow icequakes than with intermittent deep icequakes.

The rate of deep icequakes shows only weak daily fluctuations (Figure 8). There is a strong peak of activity around 14:00 local time, although this is due to a single burst of activity on 4 April. This peak disappears when removing the events for this day, although there are still slightly more events during the day than at night. The daily distribution does not change when the minimum amplitude of the events is increased from 300 µm/s 2 up to 1000 µm/s 2 , which suggests that it is not affected by daily changes in noise level and completeness threshold.

To characterize the temporal clustering of the deep events, we estimated the probability density function of the interevent times, as shown in figure 9a. For events with peak amplitude A > 100 µm/s 2 , this distribution was fitted by a power-law of exponent 1.66 for the interevent times longer than 10 s. When increasing the amplitude threshold to A > 1000 µm/s 2 , we found an exponent 1.48 for the interevent times longer than 50 s. The cut-off for times shorter than 0.1 s is likely a detection bias, as events separated by a time shorter than the signal duration are difficult to separate. The roll-off for times longer than 10 5 s (about one day) is due to the limited duration of our catalog. In the time interval of 0.1 s to 10 s, there is an intermediate power-law regime with a smaller exponent, of ≈ 0.3. This change of exponent for the interevent times of about 10 s is likely real and not due to a detection problem, because this time is much longer than the signal duration and because this crossover occurs at a longer time ≈ 50 s for the largest events with A > 1000 µm/s 2 . This distribution of interevent times is similar to that obtained for earthquakes, but with an exponent 1.6 larger than the value ≈ 1 that is commonly reported for earthquakes [START_REF] Bak | Unified scaling law for earthquakes[END_REF]]. For earthquakes, the powerlaw distribution of the interevent times was explained by earthquake interactions. Each earthquake can trigger aftershocks with a rate that decays as a power law of time since the triggering event, and increases exponentially with the magnitude of the triggering event [START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF]]. However, this model does not explain the temporal clustering of deep icequakes. We have estimated the average rate of deep icequakes as a function of the time since a deep icequake and of the amplitude of this icequake, as was done for earthquakes in [START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF]]. Unlike for earthquakes, we found that the average rate of activity after an event does not depend on its amplitude. Clustering in time might therefore arise from short-term fluctuations of an external forcing mechanism, which might be water pressure, rather than interactions between events.

The deep icequakes show a wide range of peak amplitudes, ranging from 97 µm/s 2 up to 13950 µm/s 2 . After the integration of the signals, the peak ground velocity ranged between 62 nm/s and 11474 nm/s. The distribution of the peak amplitude is shown in Figure 9. This can be fitted by a power law with an exponent of 3.5 for amplitudes between 1000 µm/s 2 and 10000 µm/s 2 . Previous studies also found a power-law distribution of icequake amplitudes, but with a smaller exponent ≈ 1 [START_REF] Roux | Microseismic activity within a sérac zone in an Alpine glacier (glacier d'Argentière, Mont Blanc, France)[END_REF], [START_REF] Faillettaz | Icequakes coupled with surface displacements for predicting glacier break-off[END_REF], Dalban Canassy et al.(2013)]. These data sets were likely dominated by shallow events, which might have a different magnitude distribution. As in the time domain, the distribution of the peak amplitudes for icequakes is similar to that observed for earthquakes, although the exponent is much larger for icequakes than earthquakes [START_REF] Bak | Unified scaling law for earthquakes[END_REF]]. [Weiss(1997)] showed that a power law distribution of peak amplitudes corresponds to an exponential distribution of magnitudes (i.e., the Gutenberg-Richter law) with the same exponent b, even if the amplitudes are not corrected for attenuation with distance (geometric and anelastic attenuation). For earthquakes, the exponent b of the Gutenberg Richter law is larger for induced earthquakes than for tectonic events [Grünthal(2013)]. For geothermal projects, where earthquakes are induced by hydraulic fracturing, the b value is 1.94, compared with b = 1.16 for tectonic earthquakes. The large b value found for our database of deep icequakes might thus be related to the triggering mechanism, which is suspected to be hydraulic fracturing. However, the failure mechanism is different for induced earthquakes (i.e., mostly shear failure) and for deep icequakes (i.e., tensile failure), and the b value might depend on the failure mechanism. Indeed, [Schorlemmer et al.(1997)] suggested that the b value for earthquakes varies systematically for different styles of faulting, with higher b values corresponding to lower stress.

Seismic noise and resonance frequency

We observed a resonance of the seismic signal during part of our acquisition. This resonance was seen both in the seismic noise and in the coda of shallow and deep icequakes. An example of a hybrid icequake, with high-frequency onset followed by lower frequency monochromatic coda is shown in Figure 10. Between 25 April and 3 May, 45% of the shallow events had a peak frequency between 30 Hz and 40 Hz, compared with only 14% before 13 April. Similar hybrid events were detected by [START_REF] West | Glacier microseismicity[END_REF]] in Alaska. These events were interpreted as brittle fractures that were driven by fluids, or hydraulic fracturing, followed by resonances caused by the rush of water into the newly opened space. The spectrogram of the continuous seismic signal is shown in Figure 11a for the time interval of 30 March to 3 May. The signal was first clipped at a maximum amplitude of 100 µ m/s 2 , to weaken the influence of seismic events. We computed the power spectrum of the vertical component of the signal for nonoverlapping time windows of 1 min. We then took the geometric average over 10 successive windows. During the first acquisition period, between 30 March and 12 April, there is no apparent resonance. Resonance at a frequency of about 35 Hz appeared during the gap between 12 and 25 April, and lasted until the end of the second acquisition on 3 May 2012. The resonance frequency f r is shown in Figure 11b. For each time window of 10 min, we normalized the power spectrum by the average spectrum, and we computed the maximum value of the normalized spectrum in the frequency range of 20 Hz to 50 Hz. The maximum value of the normalized spectrum is defined as the amplification and this is shown in Figure 11b. The main resonance frequency varied with time between 30 Hz and 38 Hz. Several other weaker resonances are also visible at higher frequencies, and these show the same temporal variations as the dominant frequency. The resonance at about 250 Hz on 28 and 29 April is also visible in the seismograms of deep icequakes that occurred during this time period. Figure 3 shows that clusters B, D, E, F and H, which were active during the harmonic tremor, had longer coda for the P-waves, with a frequency of about 250 Hz. In contrast, clusters A, C and G were only detected before 12 April, and had much simpler waveforms.

A quick drop in f r from 33 Hz down to 30 Hz occurred on 28 April at 14:00 UTC (16:00 local time), during the largest burst of the deep icequakes. The amplification increased up to a factor of 20 until midnight on 28 April, and then it decreased. A first peak of amplification was observed on 27 April, and this coincided with a small burst of deep and shallow icequakes. The strongest resonance occurred during a relatively quiet period that followed the main burst of deep icequakes. After this time, the resonance frequency and the amplification fluctuated, and the rate of shallow and deep icequakes remained higher than before 25 April. However, there was no obvious link between the resonance frequency and amplification and the rate of deep icequakes after 29 April.

Other studies have also observed seismic resonance due to fluid flow inside glaciers [START_REF] St | Hydraulic transients: A seismic source in volcanoes and glaciers[END_REF]Qamar(1979), Métaxian, et al.(2003), [START_REF] West | Glacier microseismicity[END_REF]), Heeszel et al.(2014), [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF]]. [START_REF] St | Hydraulic transients: A seismic source in volcanoes and glaciers[END_REF]] attributed the seismic source of monochromatic low-frequency icequakes to hydraulic transients that were generated by abrupt changes in water flow through subglacial conduits. Another mechanism for harmonic or hybrid signals is the resonance of fluid-filled cavities [Chouet(1986), Métaxian, et al.(2003), [START_REF] Winberry | Seismic observations of transient subglacial water-flow beneath MacAyeal Ice Stream, West Antarctica[END_REF], [START_REF] Heeszel | Humming glaciers[END_REF]]. Assuming Stoneley guided waves in a fluid-filled fracture, the resonance frequency decreases with increasing crack length and increases with crack thickness [Korneev(2008), [START_REF] Tary | Interpretation of resonance frequencies recorded during hydraulic fracturing treatments[END_REF]]. A resonance frequency of 20 Hz can be explained by a crack length of 5 m and a thickness of 2 cm [START_REF] Heeszel | Humming glaciers[END_REF]]. However, the resonance frequency alone cannot constrain the fracture geometry, because there is a trade-off between fracture width and fracture length. Observed temporal changes in f r might thus be explained by the propagation of fractures or by their progressive opening or closing. It is important to note that the relevant length that controls the resonance frequency in this model is not the full length of the fracture, but only the part that is filled with water. In our case, there are likely important temporal variations of the water level inside the fractures if they are connected with the surface. Therefore, our observed changes in the resonance frequency probably highlight temporal changes in the water level, rather than changes in the length or width of fractures. Similar observations were reported by [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF]], who observed harmonic tremor induced by a moulin, with temporal fluctuations of the resonance frequency correlated with temporal fluctuations of the water level measured inside the moulin.

Finally, gliding harmonic tremor can also be produced by the repeated failure of small subglacial asperities. [START_REF] Winberry | Nucleation and seismic tremor associated with the glacial earthquakes of Whillans Ice Stream, Antarctica[END_REF]] recorded repeating basal icequakes and a gliding harmonic tremor during the rupture propagation of a glacial earthquake underneath Whillans Ice Stream in Antarctica. The emergent tremor was composed of discrete slip events, although these events were not clearly discernible. The temporal variation of the frequency was correlated with the surface speed of the ice stream. We do not think that this mechanism can explain our observations, as in our case the resonance was also visible in the coda of the shallow icequakes (Figure 10).

Our results are similar to the study of [START_REF] Heeszel | Humming glaciers[END_REF]], who observed a gliding harmonic tremor and high frequency deep basal icequakes during a glacier lake outburst at Gornergletscher. The harmonic tremor had a slowly decaying resonance frequency of about 20 Hz with several higher frequency harmonics and several sudden step changes. The resonance frequency of deep icequakes progressively varied between 100 Hz and 200 Hz. The P-wave polarity for these events suggests that they were due to crack closing. [START_REF] Heeszel | Humming glaciers[END_REF]] interpreted the increase in the resonance frequency with time as due to closing fractures. Similar observations of gliding harmonic tremor have also been observed for volcanoes [Konstantinou and Schlindwein(2008)] and during hydraulic fracturing experiments in the hydrocarbon industry [START_REF] Tary | Interpretation of resonance frequencies recorded during hydraulic fracturing treatments[END_REF]], due to the resonance of fluid-filled cavities.

Conclusions

We detected deep icequakes at Glacier d'Argentière at depths down to 130 m, which was shallower than the ice/bed interface. Most of the events were deeper than the maximum depth of open crevasses. Although such events have been detected in previous studies [START_REF] Deichmann | Evidence for deep icequakes in an Alpine glacier[END_REF], [START_REF] Walter | Basal icequakes during changing subglacial water pressures beneath Gornergletscher, Switzerland[END_REF], [START_REF] Walter | Moment tensor inversion of icequakes on Gornergletscher, Switzerland[END_REF], [START_REF] Röösli | Sustained seismic tremors and icequakes detected in the ablation zone of the Greenland ice sheet[END_REF]], this is the first time that so many events have been detected and located. These events are clustered in time and space. The distribution of the interevent times and the amplitudes can both be fitted by power laws, although with exponents that are larger than commonly observed for earthquakes [START_REF] Bak | Unified scaling law for earthquakes[END_REF]]. We suggest hydraulic fracturing as a potential source mechanism for these events for several reasons. First, the polarities of the P-waves for all of the events are consistent with tensile failure. Fractures cannot propagate within the ice below about 30 m, and they remain open unless they are filled with water [Van der Veen(1998), Van der Veen(2007)]. Second, the largest burst of deep icequakes coincided with strong warming and rapid decrease in snow height, as well as with transient acceleration of glacier motion. These observations suggest that meltwater had infiltrated down to the ice/bed interface. Moreover, the presence of hybrid events and of a gliding harmonic tremor, which are analogous to volcanic tremor and to seismic resonance observed during industrial hydraulic-fracturing experiments, suggests the existence of water-filled cavities. Finally, the large exponent b of the amplitude distribution is reminiscent of that observed for earthquakes induced by geothermal activity and triggered by hydraulic fracturing [Grünthal(2013)].

A strong increase in the temperature occurred between 25 and 28 April and produced a decrease in the snow height of 60 cm, which was probably mainly due to snowmelting. Meltwater might have infiltrated into open fractures and promoted their propagation. These fractures were progressively filled with water and generated the resonant seismic noise and the hybrid signals. A large number of deep icequakes occurred on 28 April at greater depths than before 12 April. The acceleration of the glacier on 29 April further suggests that the meltwater had reached the ice/bed interface, although the clusters of deep icequakes detected do not reach the ice/bed interface. The harmonic tremor intensified on the evening of 28 April, and weakened in the following days. The resonance frequency showed several decreasing steps that were sometimes associated with bursts of deep icequakes. These steps might be due to crevasse propagation, which decreases the resonance frequency of water-filled cavities. The resonance frequency also showed a progressive increase between 32 Hz on 24 April and 38 Hz on 3 May. Fluctuations of the water level can modify the length of the resonating fluid-filled cavity, and can also explain variations of the resonance frequency. As several factors influence the resonance frequency (e.g., water level, fracture width and length), it is difficult to interpret our observations without additional measurements.

[ [START_REF] Walter | Deep icequakes: What happens at the base of Alpine glaciers?[END_REF]] detected deep basal icequakes in Switzerland with the seismic moment tensor corresponding to either fracture opening or fracture collapse, with both event types occurring at the same location. In this study, we did not detect any deep icequakes with dilatational first motion ("down"), that would indicate fracture closing. We searched for crack-collapse events by reversing the polarity of the template signal for cluster B and looking for similar events in the full continuous seismic signal; however, no events were found. This suggests that in our study, while fracture opening produces seismic signals, fracture closing is due to creep and is aseismic.

Many deep icequakes occurred before 12 April, when the temperature was below freezing most of the time, and when no resonance was observed. These events were thus unlikely to have been triggered by hydraulic fracturing. These events were as deep as 60 m, which was deeper than the maximum depth of the open crevasses, given by the depth where the stress intensity factor is equal to the fracture toughness of the ice [ Van der Veen(2007)]. Water-free fractures can propagate below that depth, but will then progressively close by creep, due to the ice overburden pressure. This study focussed on the time period between 30 March and 3 May 2012, when most of the deep icequakes occurred. Several weeks of seismic data are also available for February and March 2012, although very few deep icequakes were detected during that time. No deep event was detected between June 5 and 19, although the meltwater production was much higher in June than in April. An interpretation here is that the water pressure was much higher early in the spring. The water pressure decreases in summer, when an efficient drainage system has developed [START_REF] Fountain | Water flow through temperate glaciers[END_REF]] and it might not be sufficient for hydraulic fracturing. To conclude, further studies are required with a larger number of seismic sensors, and a higher sampling rate, to better constrain the location of the deep englacial icequakes and their source mechanism.

Figure 2 :

 2 Figure 1: (a) Orthophotograph of the glacier taken in 2003 (from GLACIOCIM) (a) and (b) zoom around the seismometer. The seismometer is shown by a triangle and the GLACIOCLIM meteorological station by a circle. The blue diamond in (b) marks the location of the cavitometer in a subglacial cavity. The grid origin is 45.964 degrees N and 6.973 degrees E.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure3: Waveforms of template events for all of the clusters. Dotted lines indicate the arrival times of the P-waves and S-waves. Gray patches show the variability of the signal for all of the events of each cluster (one standard deviation of signals aligned on the P-waves and S-waves, and normalized according to the peak amplitude). The number of events in each cluster and the range of peak acceleration A are indicated in each plot.

Figure 7 :

 7 Figure 7: (a) Time series of the signal peak amplitudes for each cluster. (b) Number of shallow and deep events per hour with peak amplitude A > 300 µm/s 2 . Note that the largest peak on 4 April is not to scale (507 events per hour). (c) Hourly precipitation (from the Le Tour Météo-France station) and temperature (from the GLACIOCLIM database). (d) Snow height (from Aiguilles Rouges Météo-France station) and daily basal displacement measured by the cavitometer beneath Lognan ice fall. Gray shaded areas (a,b), when no data were available.

Figure 8 :

 8 Figure8: Diurnal distribution of the icequakes for different amplitude thresholds: A > 300 µm/s 2 (blue, red) or A > 1000 µm/s 2 (black). Times are local (Greenwich Mean Time +2 h). The strong peak at 14:00 disappears when events that occurred on 4 April were excluded (red).

Figure 9 :Figure 10

 910 Figure 9: (a) Probability density function of the interevent times for different amplitude thresholds. (b) Probability density function of the peak amplitudes. Dashed lines show the power-law fits.

Figure 11 :

 11 Figure 11: (a) Spectrogram of seismic noise (in a logarithmic color-scale). The thick vertical line indicates a gap in the monitoring between 12 and 25 April, 2012. The white rectangle indicates the time and frequency interval considered in (b). (b) Temporal variation of the resonance frequency f r and amplification factor (amplitude of the normalized spectra for f = f r ) after the gap. (c) Hourly rate of shallow (solid line) and deep (histogram) icequakes. Circles show daily basal motion. (d) Temperature (from GLACIOCLIM database) and hourly precipitation (from Le Tour Météo-France station). (e) Snow height at Aiguilles Rouges Météo-France station.

  Figure 6: Icequakes depths (dots) and median values for a sliding window of 10 events (red), for all icequakes of cluster B with peak amplitude A > 300 µm/s 2 and correlation C ≥ 0.7. The blue line indicates the hour of the day (local time), and each vertical dashed line indicates a new day.

	(a) ) amplitude (µm/s 2	10 3 4 10			A B C shallow E F G D H
								no data
		03/30 10 2		04/04		04/09	04/14	04/19	04/24		04/29	05/04
	(b)	400
					A	E
	number of events /hr	0 10 20 30 40 50 60 70 80 90 03/30 icequake depth (m) 100 200 300 0 100 110 120 130 03/30 0	0 04/25	04/04 04/28 04/04	500 B C D shallow F G H icequakes 04/09 1000 04/09	1500 04/14 icequake index 2000 04/19 no data 04/14 04/19	2500 04/24 04/29 04/24	04/30	3000 04/29 05/02 04/29	3500 6 8 10 12 14 16 20 22 24 18 4 05/03 2 0	05/04 hour of day (GMT+2) 05/04

Acknowledgments. We thank the Editor, the Associate Editor and two anonymous reviewers for carefully reading the manuscript and providing many constructive suggestions. We thank Fabian Walter, Francesca Raimondi, Olivier Michel, Claudia Röösli, Pierre Dalban, Christian Vincent, Daniel Amorese, Jean-Baptiste Tary and Michel Bouchon for useful discussions. The meteorological data were obtained from Météo-France and from the GLACIOCLIM database. GLACIOCLIM is supported by INSU and OSUG. We thank Laurent Ott, Olivier Michel and others for participating in the installation and maintenance of the seismic network. This work was supported by PEPS CNRS and by a grant from Labex OSUG@2020 (ANR10 LABX56). P. Comon is supported by the ERC AdG 2013-320594 DECODA. ISTerre and SigmaPhy team are part of Labex OSUG@2020 (ANR10 LABX56). The seismic data used in this study can be downloaded from ftp://ftp.osug.fr:/pub/DATA_JGR_ARGENTIERE.