Lucie Esteoulle

Francois Daubeuf

Mayeul Collot
email: mayeul.collot@unistra.fr

Stéphanie Riché

Thierry Durroux

David Brasse

Patrice Marchand

Iuliia A Karpenko
email: i.karpenko@unistra.fr

Andrey S Klymchenko
email: andrey.klymchenko@unistra.fr

Dominique Bonnet
email: dominique.bonnet@unistra.fr

François Daubeuf

A near-infrared fluorogenic dimer enables background-free imaging of endogenous GPCRs in living mice

Keywords:

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

G protein-coupled receptors (GPCRs) are the largest family of the transmembrane receptors in humans. GPCRs are involved in virtually all aspects of human physiology in health and disease. 1 Not surprisingly, GPCRs are molecular targets of more than 30% of currently drugs on the market. 2 Each member of the GPCR family has a unique tissue-dependent expression and localization pattern which is crucial for playing its physiological and physiopathological roles. Thereby, to correlate a level of GPCR expression to a disease, it is crucial to access their spatial distribution at the cell but also the organismal levels.

Regarding the whole-organism molecular imaging techniques, fluorescence-based contrast agents emit non-ionizing radiation and have longer shelf lifetime comparing with radioisotopebased probes. Moreover, fluorescence properties of organic dyes can be modulated in a wide range by chemical modifications. [3][4][5] Although fluorescent probes are readily used in studying GPCRs in living cells, [6][7][8] the whole animal fluorescence imaging of receptors is still in its infancy. Few reports have been focused on the imaging of transgenic GPCRs in living mice. For instance, Ma et al. 9 imaged the α1-AR receptor in a xenografts model using a far-red dye-ligand conjugate. More recently, Alcobia et al. used the β2-adrenergic receptor fused to the bioluminescent reporter NanoLuc enabling the detection of the receptor-ligand binding by BRET. 10 However, to the best of our knowledge, no example of fluorescence imaging of endogenous GPCRs has been reported in mice, mostly due to the low expression level of many endogenous GPCRs and the lack of appropriate fluorescent probes.

Ideally, a fluorescent probe for the in vivo imaging of endogenous GPCRs should meet the following requirements: 1) high affinity and selectivity for its target; 2) absorption and emission in the near-infrared (NIR) region to minimize the light scattering in tissues and to enhance tissue penetration; 11,12 3) a fluorogenic character to "turn on" its fluorescence after binding to the target receptor to ensure a high signal-to-noise ratio. Fluorogenic dyes have been successfully used for background-free detection and imaging of various analytes, [13][14][15] and have been developed by us for the detection of ligand-GPCR binding in living cells. 16,17 Recently, we have reported the concept of fluorogenic squaraine dimers with environmentsensitive folding which allowed for the visualization of GPCRs in living cells in no-wash conditions with a high signal-to-noise ratio. 18 In aqueous medium, the formation of the dimer of Haggregate type resulted in complete fluorescence quenching of the probe. In contrast, once bound to the receptor, the fluorophores were exposed to a hydrophobic environment of the biomembrane, which led to dissociation of the dimer and recovery of fluorescence. Although the squaraine dimer is a powerful tool for receptor labelling in living cells, it displays absorption and emission in the far-red region which is not optimal for the in vivo imaging. This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins

Please do not adjust margins

Here, based on the dimer concept, we developed NIR fluorogenic probe of unprecedented brightness allowing for the first time the background-free detection of an endogenous GPCR, the oxytocin receptor (OTR), in living mice. The OTR is known to be involved in the modulation of complex social behavior such as social recognition, attachment, empathy, trust, [START_REF] Kendrick | Behavioral Pharmacology of Neuropeptides: Oxytocin[END_REF] and is proposed as a potential therapeutic target for the treatment of the autistic spectrum disorders. [START_REF] Yamasue | Behavioral Pharmacology of Neuropeptides: Oxytocin[END_REF] In mice, the OTR is highly expressed in the uterus during pregnancy and in the mammary glands during late pregnancy and lactation.

21 However, direct in vivo optical imaging of this receptor and GPCRs in general remains a challenge so far.

Results and discussion

The key element in the design of the NIR fluorogenic dimer probe for the OTR is the choice of the fluorophore. In addition to operate in the NIR window (700-950 nm), the fluorophore should be bright, photostable and sufficiently water-soluble.

For this purpose we designed and synthesized an original cyanine derivative 2 (Cy5.5) decorated with a PEG8 chain (Scheme 1 and ESI), which is intended to compensate for the lipophilic character of the dye and to avoid non-specific interactions.

16,18,22

The construct bearing two fluorophore copies and the OTR ligand was built using PEG linker 1, presenting two amino groups and an alkyne moiety (Scheme 1). The alkyne was introduced by performing the synthesis of 1 on SPOrT resin, developed previously in our group, 23 using a Fmoc/tBu approach and HBTU/HOBt activation. First, two PEG3 units were assembled using Fmoc-NH-PEG3-COOH amino acid, which was obtained as described by Soriano, et al. 24 Then, Fmoc-L-Lys(Fmoc)-OH amino acid was introduced to create the bifurcation. Finally, the chain was elongated with another PEG3 unit, and the terminal Fmoc protecting groups were removed. The resulting dimeric chain 1 was cleaved from the solid phase in a mixture of TFA/H 2 O/TIS 95/2.5/2.5 (v/v) and isolated by semi-preparative HPLC. The introduction of the fluorophore 2, performed in solution using PyBOP in situ activation, yielded the NIR dimer dCy5.5-PEG. In a similar way, by skipping the introduction of Fmoc-L-Lys(Fmoc)-OH, the monomeric PEG chain was synthesized and coupled to the fluorophore 2 to yield mCy5.5-PEG (Scheme 1).

To characterize the fluorogenicity resulting from the dimerization of the NIR cyanine, the absorption and fluorescence properties of the dimer dCy5.5-PEG and the monomer mCy5.5-PEG were evaluated in solvents of different polarities. Both the dimer and the monomer were highly fluorescent in organic solvents, with fluorescence quantum yields (QY) ranging from 26 to 59% (Table 1) and the fluorescence maxima situating around 710 nm. However, in contrast to the monomer mCy5.5-PEG which was fluorescent in water (QY = 22%), the fluorescence in aqueous medium of the dimer dCy5.5-PEG was almost negligible (QY = 0.4%).

To confirm the formation of the intramolecular dimer we compared the absorption spectra of two dyes. The monomer presented similar absorption spectra in water and in MeOH (Fig. 1A) with the absorption maxima around 680 nm.

Table 1. Photophysical properties of mCy5.5-PEG and dCy5.5-PEG. Solvent mCy5.5-PEG dCy5.5-PEG λabs, nm [a] λem, nm [b] QY [c] λabs, nm [a] λem, nm [b] QY [c] Although the absorption spectrum of the dimer in MeOH was identical to that of the monomer having the maximum at 682 nm, its absorption spectrum in water presented a blue-shifted maximum at 629 nm and a long-wavelength shoulder. This new band can be assigned to the non-fluorescent intramolecular dimer of the H-aggregate type, which is highly favored in aqueous medium. 18 Indeed, the intramolecular Haggregate in dCy5.5-PEG quickly disappeared upon the addition of MeOH to water, which resulted in the shift of the absorption maximum to 682 nm (Fig. 1B) and the recovery of the fluorescence (Fig. 1C). As a consequence, dCy5.5-PEG presented excellent fluorogenic properties, with up to 140-fold higher QY in organic solvents than in water (Fig. 1D). For comparison, the monomer mCy5.5-PEG is characterized by . Then, given its strong QY (56% in DMF), dCy5.5-PEG is, to the best of our knowledge, one of the brightest fluorogenic NIR dyes reported to date. Being encouraged by strong fluorogenicity of the Cy5.5 dimer, we envisaged its grafting to a ligand targeting OTR (Fig. 2). Recently, we demonstrated that carbetocin (CBT), a peptidic ligand for OTR, could be modified with bulky fluorophores with or without PEG chains at position 8 while preserving the binding properties for its receptor. 16 In Please do not adjust margins Please do not adjust margins fluorescent in aqueous solution, creating a strong background (Fig. 3). In sharp contrast, the background of the image with dCy5.5-PEG-CBT remained completely dark, probably because in solution the dimeric probe existed in the form of the nonfluorescent H-aggregate.

The ultimate challenge was to visualize the endogenous OTR on a lactating mouse model using the new probe. dCy5.5-PEG-CBT (7.5 nmol) was injected by a tail intravenous (i.v.) route in lactating Swiss mice 11 days after delivery. After 30 min, mice were imaged in a small animal living imaging system. As shown in Fig. 4A, strong fluorescence in mammary glands was detected, with practically negligible off-target signal, except for liver, the organ expected to accumulate the injected dyes. Representative images of at least 3 biological replicates.

Noteworthy, no toxicity was observed even after seven days post-injection of the probe. To demonstrate the specific labelling of OTR, dCy5.5-PEG-CBT was injected in the presence of a 60-fold excess of non-fluorescent CBT (Fig. 4B). In that case, only the liver of mice was fluorescent, leaving the mammary glands non-labelled. The absence of mammary gland labelling was also observed in naïve mice (Fig. 4C), which is not expected to overexpress oxytocin GPCR in the glands region. Finally, the administration of the monomeric probe mCy5.5-PEG-CBT resulted in a strong off-target fluorescence, which can be seen in the image (Fig. 4D) using the equivalent intensity scale (maximum value of this scale is 20-fold larger than the minimum value). These results highlight the advantage of using a fluorogenic dimer probe to increase the signal-to-noise ratio for the in vivo imaging.

Conclusions

In conclusion, based on the concept of fluorogenic dimers with environment-sensitive folding, we developed bright fluorogenic NIR probe allowing the specific, background-free and unprecedented imaging of endogenous OTR in living mice. We anticipate that this concept can be readily implemented for other GPCRs and other NIR dyes, which therefore opens up fascinating perspectives of non-invasive and non-ionizing fluorescence cartography of GPCRs in living animals.

Experimental section Chemical Synthesis.

Lys(N3)-CBT. The synthesis was performed on Fmoc Rink Amide AM resin (0.21 mmol, loading 0.7 mmol/g, 300 mg). The cleavage of Fmoc protecting groups was performed in 20% piperidine in DMF (5 mL; 2 times for 15 min). Fmoc-protected amino acids were coupled in DMF (5 mL) for 45 min using HBTU (3.8 equiv.) and HOBt (4 equiv.) with DIEA (12 equiv.) as activating agents, except the introduction of Fmoc-Cys(Mmt)-OH (5 equiv.) which was carried out using HATU (4.9 equiv.) with tetramethylpiperidine (10 equiv.) in DMF (5 mL) for 45 min. 4-Bromobutyric acid (5 equiv.) was introduced using DIC (5 equiv.) and HOBt (5 equiv.) in DMF (5 mL) for 24 hours. To remove the cysteine Mmt protecting group the peptide was treated with TFA/TIS/DCM 1/5/94 (v/v/v; 12 mL; 7 times for 2 min). The removal of Mmt was monitored by analytical RP-HPLC. The intramolecular cyclisation was performed in 1.4 M NH 3 in MeOH/THF 1/4 (v/v, 5 mL) for 4 hours at room temperature. The peptide was cleaved from the resin by TFA/H 2 O/TIS 95/2.5/5.5 (v/v/v; 15 mL) treatment for 3 hours at room temperature. The filtrate was added dropwise to 120 mL of cold Et 2 O, centrifuged for 5 min at 3000 rpm at 4 °C. The solvent was removed, the solid was washed once with cold Et 2 O, which was then removed by centrifugation for 5 min at 3000 rpm at 4 °C and decantation. The crude peptide was dried and purified by semi-preparative RP-HPLC using a linear gradient (10% to 60% in 30 min) of solvent B in solvent A, affording Lys(N3)-CBT (60 mg, 28%) as a white solid. Dimeric PEG chain 1. The synthesis was performed on SPOrT resin (0.12 mmol, loading 0.6 mmol/g, 200 mg). The cleavage of Fmoc protecting groups was performed in 20% piperidine in DMF (2 mL; 2 times for 20 min). Fmoc-NH-PEG3-COOH (2 equiv.) was introduced in DMF (2 mL) for 45 min using HBTU (1.9 equiv.) and HOBt (2 equiv.) with DIEA (6 equiv.) as activating agents. Fmoc-Lys(Fmoc)-OH (4 equiv.) was introduced in DMF (2 mL) for 45 min using HBTU (3.8 equiv.), HOBt (4 equiv.) and DIEA (12 equiv.).

Scheme 1 .

 1 Scheme 1. Synthesis of the NIR fluorogenic dimer dCy5.5-PEG and the monomer mCy5.5-PEG. Conditions: a) piperidine/DMF, b) Fmoc-NH-PEG3-COOH, HBTU, HOBt, DIPEA, DMF, c) Fmoc-L-Lys(Fmoc)-OH, HBTU, HOBt, DIPEA, DMF, d) TFA/TIS/H2O, e) 2, PyBOP, DIPEA, DMF.

Fig. 1 .

 1 Fig. 1. Spectroscopic properties of monomeric and dimeric Cy5.5 derivatives: (A) Absorption spectra of mCy5.5-PEG in MeOH and water. Absorption (B) and fluorescence (C) spectra of dCy5.5-PEG in water-MeOH mixtures. D) Ratio of QY in organic solvents to that in water for mCy5.5-PEG and dCy5.5-PEG.

Fig. 2 .

 2 Fig. 2. Structures of the conjugates of the NIR fluorogenic dimer and the corresponding monomer with the carbetocin ligand. only < 2.7-fold difference in QY between water and organic solvents. As the absorption spectrum of dCy5.5-PEG in the

Fig. 3

 3 Fig. 3 Confocal microscopy studies of dCy5.5-PEG-CBT and mCy5.5-PEG-CBT on living HEK 293 cells expressing OTR-GFP fusion under no-wash conditions. Cells were incubated with the ligands for 5 min at room temperature prior to the imaging.

Fig. 4 .

 4 Fig. 4. In vivo images of lactating (A, B, D) or naïve (C) mice injected i.v. with 7.5 nmol of dCy5.5-PEG-CBT (A and C), 7.5 nmol of dCy5.5-PEG-CBT and 450 nmol of CBT (B) or 7.5 nmol of mCy5.5-PEG-CBT (D) 30 min prior to the imaging.

  t R =11.04 min (>95% purity [220.8 nm]); HRMS (ESI) calc. for C 45 H 68 N 14 NaO 12 S ([M+Na] + ): 1051.4760; found: 1051.4776.

-PEG-CBT.

  CuSO 4 (1 equiv., 2.04 μmol, 20.4 μL of 0.1 M aqueous solution), sodium ascorbate (1.2 equiv., 2.44 μmol, 24.4 μL of 0.1 M aqueous solution) and TBTA (1.2 equiv., 2.44 μmol, 24.4 μL of 0.1 M DMF solution) were pre-activated during 20 min at room temperature in a total volume of water/DMF 2/8 (v/v) of 70 µL. Lys(N3)-CBT (1.2 equiv., 2.44 μmol, 2.52 mg) and dCy5.5-PEG (1 equiv., 2.04 μmol, 6.6 mg) were added to the mixture followed by 700 µL of water/DMF 2/8 (v/v). The reaction mixture was stirred for 3 hours at 37 °C. The crude product was purified by semi-preparative RP-HPLC using a linear gradient (20% to 70% in 30 min) of solvent B in solvent A, to obtain the desired product(3.5 mg, 40%) as a blue solid. t R = 14.2 min (>95% purity [220.8 nm]); HRMS (ESI) calc. for C 209 H 309 N 30 O 48 S ([M+H] CuSO 4 (1 equiv., 3.25 μmol, 32.5 μL of 0.1 M aqueous solution), sodium ascorbate (1.2 equiv., 3.9 μmol, 39 μL of 0.1 M aqueous solution) and TBTA (1.2 equiv., 3.9 μmol, 39 μL of 0.1 M DMF solution) were pre-activated during 20 min at room temperature in a total volume of water/DMF 2/8 (v/v) of 100 µL. Lys(N3)-CBT (1.2 equiv., 3.9 μmol, 4.01 mg) and mCy5.5-PEG (1 eq., 3.25 μmol, 6 mg) were added to the mixture followed by 1100 µL of water/DMF 2/8 (v/v). The reaction mixture was stirred for 3 hours at 37 °C. The crude product was purified by semi-preparative RP-HPLC using a linear gradient (20% to 70% in 40 min) of solvent B in solvent A, to obtain the desired product (4.1 mg, 44%) as a blue solid. t R = 12.8 min (>95% purity [220.8 nm]); HRMS (ESI) calc. for C 137 H 207 N 24 O 34 S ([M+2H] 3+ /3): 921,4978; found: 921.4959. This journal is © The Royal Society of Chemistry 20xx Please do not adjust margins Please do not adjust margins

			Please do not adjust margins
	min). Fmoc-NH-PEG3-COOH (2 equiv.) was introduced in DMF
	(0.5 mL) for 45 min using HBTU (1.9 equiv.) and HOBt (2 equiv.)
	with DIEA (6 equiv.) as activating agents. The monomeric chain
	was cleaved from the resin by TFA/H 2 O/TIS 95/2.5/2.5 (v/v/v)
	treatment for 3 hours at room temperature. The filtrate was
	precipitated with cold Et 2 O, centrifuged for 5 min at 3000 rpm
	at 4 °C and the solvent was removed by decantation. The
	residue was washed with cold Et 2 O, centrifuged one more time
	and the solvent was removed by decantation. The crude
	peptide was dried and purified by semi-preparative RP-HPLC
	using a linear gradient (0% to 30% in 40 min) of solvent B in
	solvent A, to obtained the monomeric PEG chain 3 (13.3 mg,
	24%) as a brown oil. t R = 5.48 min (>95% purity [220.8 nm]);
	MS (ESI): calc. for C 36 H 66 N 8 O 13 ([M+2H]	2+ /2) 409.24; found
	409.24.
	dCy5.5-PEG. Pegylated cyanine 2 (1.9 equiv., 13 mg, 0.012
	mmol) and the dimeric PEG chain 1 (1 equiv., 9 mg, 0.006
	mmol) were solubilized in 262 µL of dry DMF. PyBOP (2 equiv.,
	6.66 mg, 0.012 mmol) and DIEA (6 equiv., 6.36 μL, 0.0385
	mmol) were added to the mixture. The reaction mixture was
	stirred for 1 hour at room temperature. The crude product was
	purified by semi-preparative RP-HPLC using a linear gradient
	(15% to 60% in 30 min) of solvent B in solvent A, to obtain
	dCy5.5-PEG (16 mg, 77%) as a blue solid. t R = 14.71 min (>95%
	purity [220.8 nm]); MS (ESI): calc. for C 164 H 241 N 16 O 36
	([M+H]	3+ /3) 1003.58; found 1003.58.
	mCy5.5-PEG. cyanine 2 (1.2 equiv., 9.59 mg, 0.009
	mmol) and the monomeric PEG chain 3 (1 equiv., 7 mg, 0.008
	mmol) were solubilized in 308 µL of dry DMF. PyBOP (1.2
	equiv., 4.7 mg, 0.009 mmol) and DIEA (6 equiv., 7.46 μL, 0.045
	mmol) were added and the reaction mixture was stirred for 1
	hour at room temperature. The crude product was purified by
	semi-preparative RP-HPLC using a linear gradient (15% to 60%
	in 30 min) of solvent B in solvent A, to obtain the desired
	product (8.1 mg, 58%) as a blue solid. t R = 13.15 min (>95%
	purity [220.8 nm]); MS (ESI): calc. for C 92 H 138 N 10 O 22 ([M+H]	2+ /2)
	867.50; found 867.50.
	dCy5.5		The dimeric chain was
			cleaved from the resin by TFA/H 2 O/TIS 95/2.5/2.5 (v/v/v)
			treatment for 3 hours at room temperature. The filtrate was
			precipitated with cold Et 2 O, centrifuged for 5 min at 3000 rpm
			at 4 °C and the solvent was removed by decantation. The
			residue was washed with cold Et 2 O, centrifuged one more time
			and the solvent was removed by decantation. The crude
			peptide was dried and purified by semi-preparative RP-HPLC
			using a linear gradient (0% to 30% in 40 min) of solvent B in
			solvent A to obtain the dimeric PEG chain 1 (20 mg, 23%) as a
			brown oil. t R = 5.51 min (>95% purity [220.8 nm]); MS (ESI):
			calc. for C 52 H 96 N 12 O 18 ([M+2H]	2+ /2) 588.35; found 588.35.
			Monomeric PEG chain 3. The monomeric chain was
			synthesized on a SPOrT resin (0.060 mmol, loading 0.6
			mmol/g, 100 mg). The cleavage of Fmoc protecting groups was
			performed in 20% piperidine in DMF (0.5 mL; 2 times for 20
			Please do not adjust margins

3+ /3): 1346,4127; found: 1346.4097. mCy5.5-PEG-CBT.

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 3

Acknowledgements

This work was supported by the Agence Nationale de la Recherche (FluoroPEP, ANR-16-CE18-0030), ERC Consolidator grant BrightSens 648528, the IdEx Unistra, the LabEx MEDALIS (ANR-10-LABX-0034), the Centre National de la Recherche Scientifique, the University of Strasbourg. L.E. was supported by a fellowship from the Ministère de l'Education Nationale, de l'Enseignement Supérieur et de la Recherche. We are grateful to Dr. Dominique Bagnard (INSERM 1119 -BMNTS Lab) for giving us access to the NightOwl luminograph, Christel Valencia and Dr. Pascal Villa (Plate-forme de Chimie Biologique Intégrative de Strasbourg, UMS3286) for prodiving us with the HEK293 Cells, Dr. Delphine Garnier and Dr. Estefania Oliva for mass spectrometry (PACSI platform GDS3670) and Dr. Dmytro Dziuba for fruitful suggestions and proofreading the manuscript.

Notes and references

Conflicts of interest

There are no conflicts to declare.


Ludovic Moreau
email: ludovic.moreau@univ-grenoble-alpes.fr

Pierre Boué

Agathe Serripierri

Jérôme Weiss

Daniel Hollis

Ildut Pondaven

Stéphane Garambois

Éric Larose

Agnès Helmstetter

Laurent Stehly

Benjamin Vial

Agnés Helmstetter

Gregor Hillers

Olivier Gilbert

Sea Ice Thickness and Elastic Properties From the Analysis of Multimodal Guided Wave Propagation Measured With a Passive Seismic Array

published or not. The documents may come

Introduction

The effects of global climate changes are strongest in the Arctic, which is currently under-going the most intense warming of any region on Earth, for example, IPCC (2014). A spectacular signature of these changes is the accelerated decline of the sea ice cover, at a rate faster than that forecasted by climate models, both in terms of ice extent [START_REF] Stroeve | Arctic sea ice decline: Faster than forecast[END_REF] or average thickness [START_REF] Rampal | IPCC climate models do not capture Arctic sea ice drift acceleration: Consequences in terms of projected sea ice thinning and decline[END_REF]. Significant multi-decadal trends are also observed for mechanically related variables, such as an acceleration of sea ice drift and deformation [START_REF] Rampal | Arctic sea ice velocity field: General circulation and turbulent-like fluctuations[END_REF] and a mechanical weakening at the Arctic basin scale [START_REF] Gimbert | Recent mechanical weakening of the Arctic sea ice cover as revealed from larger inertial oscillations[END_REF].

Arctic sea ice has an influence on climate at the global scale, because of the key role it plays in several atmosphere-ice-oceans interactions, such as the thermohaline circulation in icy waters, or the Arctic oscillation [START_REF] Chylek | Arctic air temperature change amplification and the Atlantic Multidecadal Oscillation[END_REF][START_REF] Uttal | Surface heat budget of the Arctic Ocean[END_REF]. Positive feedbacks between sea ice decline, modification of the albedo of the Arctic ocean, and the warming of arctic waters, likely influences this evolution and may lead to a sea-ice-free Arctic in summer within few decades, an emblematic example of a climatic tipping point [START_REF] Steffen | Trajectories of the Earth System in the Anthropocene[END_REF][START_REF] Stroeve | Trends in Arctic sea ice extent from CMIP5, CMIP3 and observations[END_REF]. Hence, understanding this sea ice decline in conjunction with global warming is one of the main challenges of modern climatology. In this matter, a finer description of the dynamic and thermodynamic processes involved in sea ice models is needed. However, accurate parameterization of these models require more complete data sets of sea ice properties: salinity, temperature, mechanical properties, and thickness.

The use of seismic methods to study sea ice has been considered for more than 60 years. Sea ice acts as a waveguide for seismic waves. As such, it supports dispersive and multimodal wave propagation, which can be used to estimate its mechanical properties and thickness. Pioneering works in this matter are the experiments by [START_REF] Anderson | Preliminary results and review of sea ice elasticity and related studies[END_REF], who inferred the elastic constants of the ice from seismic waves velocities, and by [START_REF] Hunkins | Seismic studies of sea ice[END_REF], who demonstrated the influence of temperature on the elastic constants of the ice, and used the dispersion of the flexural wave to estimate its thickness. [START_REF] Yang | Experimental characterization of elastic waves in a floating ice sheet[END_REF] measured broadband multimodal propagation, using a linear array of triaxial geophones and hammer blows. [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF] inferred the elastic properties and thickness of sea ice by inverting seismic data obtained with a similar experimental setup.

Table 1 compares sea ice properties measured from several sets of seismic data. Young's modulus extends between 3.1 and 8.5 GPa, and these should not be considered as the extremes values. These differences are, for the most part, essentially related to ice composition. Thinner and younger ice tends to be more porous and with a larger brine content [START_REF] Weeks | The growth, structure and properties of sea ice[END_REF], resulting in lower values of Young's modulus [START_REF] Mellor | Mechanical behavior of sea ice[END_REF].

Due to the hostile environment and heavy logistics in polar regions, such studies have remained rare. With the need of exhaustive and accurate field data to meet the issue of climate change, and the rapid technological and methodological advances of the last decade, there has been a reconsideration of seismic methods for polar seismology [START_REF] Kanao | Recent advance in polar seismology: Global impact of the International Polar Year[END_REF]. In particular, passive seismology has addressed difficulties related to the need of active sources. Based on continuous recordings of the ambient seismic field, it is possible to turn seismic stations into virtual sources via basic signal processing [START_REF] Sabra | Extracting time-domain Green's function estimates from ambient seismic noise[END_REF][START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]. This was recently applied for passive seismic studies of sea ice. [START_REF] Marsan | Sea-ice thickness measurement based on the dispersion of ice swell[END_REF] inferred Arctic sea ice thickness from measurements of the flexural wave dispersion. [START_REF] Sutherland | Observations of wave dispersion and attenuation in landfast ice[END_REF] measured the dispersion and attenuation of the flexural wave in landfast ice. [START_REF] Marsan | Characterizing horizontally-polarized shear and infragravity vibrational modes in the Arctic sea ice cover using correlation methods[END_REF] observed the propagation of both the flexural and shear-horizontal waves from ambient seismic noise recordings in the Arctic.

The present paper follows up on our work initiated at the laboratory scale to infer the thickness and mechanical properties of the ice from the observation of guided wave modes [START_REF] Moreau | Monitoring ice thickness and elastic properties from the measurement of leaky guided waves: A laboratory experiment[END_REF]. In a layer of ice floating on water, the wavefield is composed of modes that are of the same nature as those that propagate in a stress-free plate, that is, the Lamb wave modes and the guided shear-horizontal (SH) wave modes. Lamb waves are decomposed in two families known as antisymmetric modes, A n , and symmetric modes, S n [START_REF] Lamb | On waves in an elastic plate[END_REF]. We distinguish between the fundamental modes (n = 0) that propagate at all frequencies and the higher-order modes (n = 1, 2, 3, ...) that become propagative only at frequencies above their cutoff frequency. Antisymmetric modes produce an axial displacement component that is antisymmetric and a normal displacement component that is symmetric, with respect to the middle-plane of the plate. For symmetric modes, this is the opposite. As such, antisymmetric modes generate flexural motion, and symmetric modes generate traction-compression motion.

The solid-liquid interface, however, induces a change in the boundary conditions. This does not modify the propagation of the guided SH modes (because there is no coupling with shear motion in water), but it breaks the symmetry of the problem for the Lamb modes. As a consequence, the elastic wavefield in a free ice layer is different from that in a floating ice layer. The main differences are the existence of an additional mode, similar to the Scholte mode that propagates at the liquid-solid interface between two semi-infinite media, as well as a modification of the dispersion branches of the guided waves. For these reasons, the wavefield in a floating ice layer is referred to as quasi-symmetric (QS n ), quasi-antisymmetric (QA n ), and quasi-Scholte (QS) modes.

Below the cutoff frequency of the higher-order guided modes, the elastic wavefield in the ice is theoretically distributed on four propagating modes: QS, QA 0 , QS 0 , and the fundamental shear-horizontally polarized mode (SH 0 ). The energy distribution of these modes on the displacement components depends on the product between the frequency of the exciting wavefield and the thickness of the ice. For example, when this product remains under 50 Hz•m, most of the energy on the axial and transverse components comes from the QS 0 and SH 0 modes, while on the vertical component of the displacement the energy comes almost exclusively from the QS mode, because of the strong energy leakage of the QA 0 mode in water. In such conditions, the wavefield can be approximated with only three propagating waves: the flexural wave, the axial wave, and the shear-horizontal wave [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF]. These are essentially the asymptotic behavior of the guided modes, when the displacement field across the ice thickness is considered linear for the QS mode and constant for the QS 0 mode, as shown in Figure 1a. This approximation is valid only for frequency-thickness values under 50 Hz•m. For higher values, the approximation no longer holds, and the dispersion branch of the QS mode starts to deviate from that of the flexural wave (Figure 1b). If frequency-thickness is increased further, the energy on the vertical component of the displacement becomes progressively dominated by the QA 0 mode, the QS 0 mode becomes dispersive, and higher-order modes become propagative [START_REF] Moreau | Monitoring ice thickness and elastic properties from the measurement of leaky guided waves: A laboratory experiment[END_REF].

As far as seismology on sea ice is concerned, when active sources are used it is possible to measure wave propagation at frequency-thickness values exceeding 1000 Hz•m. There have also been reports of higher-order modes propagating at lower frequency-thicknesses [START_REF] Yang | Experimental characterization of elastic waves in a floating ice sheet[END_REF], although the origin of these modes could not be explained, because the maximum frequency-thickness was around 100 Hz•m values. For passive seismology (i.e., seismology without active sources), the majority of reports are for environments characterized by frequency-thickness below 0.5 Hz•m [START_REF] Marsan | Sea-ice thickness measurement based on the dispersion of ice swell[END_REF][START_REF] Sutherland | Observations of wave dispersion and attenuation in landfast ice[END_REF], where wave excitation comes mainly from wind and swell. In this paper, we demonstrate that multimodal and broadband propagation of guided waves can be extracted far beyond this value from passive recordings and used to monitor sea ice properties. To this end, an array of 247 geophones was installed in the Van Mijen fjord (Svalbard) to record the ambient seismic field for 27 days. Recordings exhibit wave propagation at frequencies up to the sampling limit of 500 Hz, with a level of energy showing that waves may propagate at even higher frequencies.

Using array processing methods, we analyze several aspects of this rich data set to demonstrate that passive seismic recordings on sea ice can produce results at least as good as active seismic experiments.

In section 2, we introduce the array, the seismic acquisition, and the ground penetrating radar profiles. In section 3, we discuss the properties of the seismic wavefield and the dispersion curves of sea ice are calculated on 1 and 5 March 2019 from active sources, icequakes, and seismic noise. They are subsequently inverted for sea ice thickness, Young's modulus, and Poisson's ratio, yielding values that are consistent with in situ measurements and observations, as well as with previous reports on elastic properties for first-year ice.

Material and Methods

The aim of the deployment is to demonstrate that broadband propagation of guided waves can be extracted from passive seismic data measured in the ice and used to estimate ice thickness and its mechanical properties. A total of 247 FairFieldNodal Zland geophones (flat frequency response down to the cutoff frequency of 10 Hz for 1C and 5 Hz for 3C) were deployed on the ice of Van Mijen fjord in Vallunden Lake (Figure 2a). This part of the fjord, located 2 km south-east of Sveagruva, is surrounded by a moraine (Figure 2b). Because of its connection to the fjord by a canal of width ∼10 m, Vallunden Lake is subject to tidal flows and forcing.

It is a quiet zone that is well suited for scientific experiments on sea ice. Each node is equipped with a spike for better coupling with the ground. However, this is not ideal in the ice. Instead, the nodes were installed directly in the ice without their spike. To maximize the coupling, a milling tool was specifically designed to drill the ice at the diameter of the nodes. They were installed in the holes at about half their height and covered back with snow to insulate them for preserving their battery life. At the time of the deployment, the internal temperature of several nodes was measured, before and after covering them with snow, showing an increase from -21 to -16 • C .

Array Deployment

Recording the three components of the displacement is technically not necessary to measure the QS and QS 0 modes. However, the fundamental shear-horizontal (SH 0 ) mode can be measured only on the horizontal channels. Moreover, the QS 0 mode is easier to measure on the horizontal channels, because it produces a weak vertical displacement at frequency-thickness values below 1000 Hz•m. Therefore, in practice the use of 3C geophones helps with data interpretation. However, due to limited geophones availability, the deployment contains both 1C and 3C stations. We set the sampling frequency of the 1C and 3C stations to 500 and 1000 Hz, respectively. A consistent sampling rate of 1000 Hz was prohibitive because of the shorter battery life of the 1C geophones.

The geophones were divided into five zones, including one large, two-dimensional array, and four smaller linear arrays. The distribution of the sensors is presented in details in Figures 3a and3b, with blue and red dots indicating 1C and 3C geophones. The main array is a large squared area with sides of length 48 m (Figure 3b). The other four arrays consist of four stations with an aperture of 12 m. It is noteworthy that, during deployment, a 1C geophone was replaced with a spare 3C sensor because it was deficient (see station 133-517 in Figure 3b). Moreover, stations 113-509 and 125-542 (marked with a cross in Figure 3b) were deployed but failed to record, which leaves the main array with 231 operating sensors. The purpose of the main array is to measure guided waves propagating in all directions. The linear arrays act as distant virtual sources, for the calculation of the noise correlation functions [START_REF] Sabra | Extracting time-domain Green's function estimates from ambient seismic noise[END_REF][START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF].

Based on an ice thickness of 50 cm, the spatial sampling of the 1C geophones in the main array (4 m) is sufficient to measure the slowest mode, QS, and the fastest mode, QS 0 , in the 1-20 Hz and 30-250 Hz frequency ranges, respectively. With the denser spatial sampling of the 3C geophones (1 m), these two modes can be measured in extended frequency ranges: 1-150 Hz for QS and 30-500 Hz for QS 0 . This high spatial density is not necessary for monitoring ice thickness and mechanical properties, because the guided wave inverse problem is very well constrained. The quality of the inversion depends essentially on the number of modes that can be measured and the bandwidth on which they are measured. The more modes and the larger this bandwidth, the better constrained the inversion. The SH 0 mode is not dispersive, and the QS 0 mode becomes dispersive only at much higher frequency-thickness values (above 1000 Hz•m). If the QS, QS 0 , and SH 0 modes can all be measured, then a well-constrained inversion can be expected even in narrower frequencies bands, for example, in the 1-5 Hz bandwidth. In that case, a spatial sampling of 10-15 m could be sufficient. Of course, for sea ice with a thickness x times greater than that at Vallunden Lake, the spatial sampling can be reduced by the same factor x. One of the objectives of this experiment is to underpin research on the trade-off between accuracy and array density and should help setting new benchmarks for future inversions of sea ice seismic data.

Field

Ground Penetrating Radar Survey

On 28 February 2019, two linear GPR profiles were acquired over a distance of 250 m with 800 MHz shielded antenna of Mala Geosciences: a north-south profile between stations 501-104 and 301-104 and an east-west profile between stations 401-104 and 201-104. Both GPR profiles were processed with a 450-1100 MHz band-pass filter and late arrivals amplification with a time squared gain function. Figure 4 shows the resulting processed images for both profiles, which are dominated by the air wave followed by a large reflection, whose first onset varies between 2 and 3 ns along the profile. This could represent the base of the snow cover. This large event is followed by other reflections with the same polarity, which might be multiples. An additional waveform, which arrives around 10 ns, exhibits slight variations with an opposite reflectivity. This could indicate the reflection at the sea ice/sea water interface but with these GPR data, such interpretation remains speculative.

The velocity of electromagnetic waves in snow strongly depends on its properties [START_REF] Liu | Measurement of dielectric permittivity and thickness of snow and ice on a brackish lagoon using GPR[END_REF] and ranges from 16.8 cm/ns (highly dense compacted snow with a very weak porosity) to 25.5 cm/ns (very porous fresh snow). In the absence of any common middle point acquisition, we use an average velocity of 21.1 cm/ns. From this estimate, the snow thickness varies between 21.1 cm (2 ns) and 32.2 cm (3 ns), which is consistent with observations in the field, despite a large uncertainty of ±20%.

In sea ice, the velocity may vary between 10.6 and 15 cm/ns, depending on its brine content and temperature. Assuming an average velocity of 12.8 cm/ns, the second reflection would correspond to an ice thickness around 52 cm with a large uncertainty of ±17%. This range is consistent with the local drilling information (62 cm) and with estimates from seismic methods (54 cm, see section 3.5). Here, the interpretation of the GPR images appears speculative and uncertain. As noted by [START_REF] Liu | Measurement of dielectric permittivity and thickness of snow and ice on a brackish lagoon using GPR[END_REF] this problem appears when high-frequency antennas are used, because the high conductivity of sea ice drastically decreases the penetration of a GPR signal. Consequently, if GPR is able to detect the snow cover with a good resolution, the detection of sea ice bottom remains uncertain and strongly depends on the saline properties of sea ice, which is not the case with seismic methods. The vibrating sources and artificial noise will be useful for future in-depth analyses of signals with correlation-based processing methods. For example, snowmobiles are a relatively high-frequency seismic noise source, up to ∼200 Hz. Such acquisitions are therefore of interest for investigating the convergence of the noise correlation function.

Vibrating Sources, Impulsive Sources, and Artificial Noise Sources

Note that, in this section, we are introducing the full seismic experiment so that it can be referred to in our future works on this rich data set. However, in the following, only the impulsive sources to the east of the main array are considered, that is, those located at stations 401-101 to 401-104 (see Figure 3a). The purpose here is twofold. First, we are demonstrating that the dispersion curves obtained with passive data are comparable to those obtained with active sources. Second, we are demonstrating that the noise correlation function converges quickly toward the Green's function of the ice layer, and by extension that it can be used for monitoring ice properties with a temporal resolution that is less than a day.

Results

Noise Properties: Power Spectral Density

Figure 5 shows the spectrogram of one full day of recording on 5 March 2019, from station 125-525 at the center of the main array. The spectrogram contains thousands of impulsive and broadband events. Part of these events are probably due to the specific conditions of the ice at the Lake, where it is directly in contact with the shores and influenced by the mechanical stress induced by tidal effects. For example, tidal cracks were observed along the shoreline. Some cracks were frozen or some partially frozen, depending on the tidal phase. Water in and outflow through the connecting channel had a maximum speed of 3 m/s and probably contributed to the triggering of icequakes too. However, because they occur at any time and on every day, it is fair to assume that there is also a large amount that is not connected to the specific situation at the lake. For the same reason, it is also fair to assume that many icequakes also happen in freely floating sea ice, where plate inner deformation due to relative drift, swell, and thermal expansion of the ice, are known to induce cracking.

There is also a more consistent trend at lower frequencies between 1 and 25 Hz, where signal exhibits stronger energy levels. The origin of this seismic noise is a mix between wind and water current, under the ice, in the canal. Note the presence of more energetic signals between 8:30 a.m. and 5:30 p.m. They were generated by research field works that took place about 300 m north-east of the main array. They lasted 12 days, between 4 and 16 March 2019. Since these noise sources are not natural and very directional, they are not of interest for this study and will be disregarded in the following. Therefore, next we focus on data recorded in the absence of human activity.

Noise Properties: Wavefield Isotropy

In order to identify the azimuthal distribution of the wavefield around the array, we apply classical beamforming to the data. This processing aims at investigating a seismic wavefield in terms of its directionality and phase slowness, based on its spatial coherency across an array [START_REF] Rost | Array seismology: Methods and applications[END_REF]. When applied to a full day of data, it provides insights in the origin of the seismic noise. When applied to a few seconds of data, this allows the identification of impulsive microseismic events as well as their direction of propagation. In practice, beamforming consist of delaying the waveforms with the phase slowness of a plane wave propagating through the array with a given azimuthal direction, and summing the resulting signals. The slowness and angle that return the maximum output are closest to those of the wavefield recorded in the time window.

Because beamforming relies on the spatial coherency of the wavefield, it is sensitive to spatial sampling. As such, it is restricted to frequencies where geophones spacing allows a sampling of the wavelength that satisfies Nyquist's criterion to prevent aliasing. Spatial sampling in the main array is limited by the 1C stations, separated with a 4-m spacing. We are therefore restricted to frequencies where the wavelength is at least 8 m, that is, under 16 Hz for the QS mode.

We applied this processing to 12 hr of noise recording without nearby human activity on 3 March 2019 between 0h00 and 12h00. Figure 6 shows the beamforming obtained from the vertical displacement in four frequency bands: 3-5, 7-9, 11-13, and 15-17 Hz. This displacement component is strongly dominated by the QS mode at such frequencies; hence, the beamformer is mainly associated with this mode. Consequently, the aliasing limit mentioned previously is visible in the 15-17 Hz, where beamformer overlapping starts to occur near a slowness of 8 s/km.

The circular amplitude of the beamformer around the array center indicates that noise sources are distributed in all azimuthal directions. The intensity is slightly lower to the south of the array. This is likely the consequence of the array position relative to the edges of the lake, which are closer to the north, east, and west than they are to the south. Hence, waves emitted from the southern part of the lake may have a lower amplitude due to more energy leakage in water. The isotropic ambient wavefield allows to look for specific microseismic events to act as sources that originate from a desired direction, in view of computing the dispersion curves of sea ice.

Typical Microseismic Events

In this section, we focus on the first 8 hr of recording on 5 March 2019. The corresponding signal is shown in Figure 7 for station 125-525. It contains hundreds of events with amplitudes that vary by up to two orders of magnitude, as shown in the between 1:30 and 2:30 a.m. By scanning this signal with 4-s-long moving time window and applying beamforming to the data, it was possible to identify events with significant amplitude that originated from the east. Three of such events are shown in color in Figure 7. The figure also shows the beamforming associated with the third event. It was computed in the 15-17, 11-13, 7-9, and 3-5 Hz frequency bands, velocity of 139, 118, 93, and 60 m/s, respectively, which reveals the strong dispersion of the signals. Figure 8 shows a zoom on event 3, as well as the distribution of the vertical displacement component in the main array. These snapshots of the wavefield were extracted at the times marked with a red cross in event 3. For better visualization, data were first interpolated on a denser grid with a 1-m interval. This is a clear visualization of seismic wave dispersion in sea ice, with high frequencies arriving first, which is typical of the QS and QA 0 modes at low frequencies-thickness values. It is noteworthy that for values of the order of 1 Hz•m, the displacement is dominated by the QS mode on the vertical component, the QS 0 mode on the radial component, and the SH 0 mode on the horizontal component. Hence, the QA 0 mode contributes to the wavefield with a negligible effect, and the high dispersion observed in the beamforming is that of the QS mode, which is used next to estimate ice thickness.

Dispersion Curves of Sea Ice

In this section we calculate the frequency-wavenumber spectrum of the waveforms propagating along line 125 of the main array, where the denser spatial sampling allows measurements without aliasing up to 150 Hz for the three modes. This analysis is performed on 1 and 5 March 2019 to recover the dispersion curves of sea ice, based on the wavefield generated by active sources and icequakes, or the wavefield reconstructed from the seismic noise via correlation methods. We outline specific array processing that takes advantage of the multiplicity of sources for enhanced calculation of the dispersion curves.

Array Processing for Calculating Dispersion Curves

A classical way of computing the dispersion curves of a wavefield is to apply a Fourier transform to the time and space dimensions of the waveforms recorded along a line of geophones, which yields the frequency-wavenumber spectrum. If several sources are available, the spectra can be averaged over the number of sources. However, signals with a low signal-to-noise ratio (SNR) may result in frequency-wavenumber spectra of poor quality. A better way of taking advantage of several sources consists in combining them, following the processing described in [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF]. The reader is referred to [START_REF] Minonzio | Guided wave phase velocity measurement using multi-emitter and multi-receiver arrays in the axial transmission configuration[END_REF] and [START_REF] Moreau | Monitoring ice thickness and elastic properties from the measurement of leaky guided waves: A laboratory experiment[END_REF] for more details regarding its implementation and application to the present problem. Only the main steps of this processing are outlined here.

1. The matrix of transmit-receive signals has three dimensions: sources, receivers, and time. The first step is the application of the Fourier transform to the temporal dimension of this matrix. 2. At each frequency, the resulting Fourier-domain matrix is sliced into two-dimensional transmit-receive matrices. These matrices are then decomposed into singular values, with the underlying idea that the different levels of modal energy are distributed onto a set of orthogonal singular vectors. This allows a heuristic separation of the noise and signal subspaces, in a classical way for singular value-based filters. 3. The last step consists of defining test vectors that are representative of the wave propagation problem in the receivers basis, for example, plane waves with a given wavenumber. The projection of test vectors onto the singular vectors of the receivers basis leads to a scalar product that is maximized when the wavenumber in the test vector matches that of the waves in the measured wavefield.

This processing significantly enhances the identification of the dispersion branches in the frequency-wavenumber space, for two reasons: (i) the separation of the noise and signal subspaces and (ii) the projection of test vectors is not weighted by singular values. Yet the energy information is only contained in singular values, not in singular vectors, which all have a norm equal to unity. Hence, in theory all modes stand out with the same spectrum intensity, despite their different relative energy in the wavefield.

Dispersion Curves of Sea Ice from Passive and Active Seismic Events

The first data set used for computing the dispersion curves is the matrix of signals recorded during the 5 s following the active impulsive sources at stations 401-101 to 401-104, on 1 March 2019. The second data set is the matrix of signals recorded on 5 March 2019 during the 5 s of propagation of the three microseismic events shown in Figure 7. They exhibit energy at frequencies from 0 Hz up to at least the sampling limit of 500 Hz. These events all originate from a direction that almost coincides with that of line 125, the difference being less than 7 • . Hence, from the apparent wavenumber measured along the stations in row 125 during the propagation of these events, the wavenumber of the guided modes propagating in the actual east-west direction can be estimated with less than 1% error. A first attempt at calculating the dispersion branches was made by applying a classical frequency-wavenumber Fourier transform to the data sets, and averaging the 10.1029/2019JC015709 spectra over the number of sources. Only the QS appears clearly in the resulting diagrams, in a frequency band limited between 5 and 35 Hz. The QS 0 mode is barely visible, and the SH 0 mode is not detected (see Figures B1-a and B1-b). This is essentially due to the fact that most of the energy goes to the QS mode at these frequency-thickness values, and to energy leakage in water.

Application of the processing described in section 3.4.1 to the three components of the displacement significantly enhances the dispersion curves. The resulting diagrams for the active and passive data sets are shown in Figures 9a and9b.

These dispersion curves can be converted into modal velocities, using the following three-step process.

1. For each frequency, extract the couples (f m , k m ) with maximum intensity in the dispersion curves. This is performed for the three modes: m = Q S , Q S 0 , SH 0 . 2. Convert the frequency-wavenumber couples, (f m , k m ), into frequency-velocity couples, (f m , c m ), with c m = 2𝜋f m ∕k m . 3. For each mode, find the best fit of all (f m , c m ) couples. S 0 and SH 0 are not dispersive at the frequencies of interest; hence, we used linear interpolation. For the QS mode, however, we used spline interpolation to account for the nonlinear pattern.

It is also possible to calculate the standard deviation of the absolute error of the interpolations. The extracted velocities and the histogram of the interpolation error are shown in Figure C1. In the following we use these standard deviations as a measure of uncertainty for the velocities. They are also converted into the standard deviation of the estimated Young's moduli, Poisson's ratios, and thicknesses.

From the active data set recorded on 1 March 2019, we measure nondispersive propagation of the SH 0 and QS 0 modes at a velocity of 1,260 ± 15 and 2,200 ± 84 m/s, respectively, while the QS mode propagates with a velocity varying between 58 ± 1 m/s at 4 Hz and 277 ± 1 m/s at 60 Hz. From the passive data set recorded on 5 March 2019, we measure nondispersive propagation of the SH 0 and QS 0 modes at a velocity of 1,360 ± 48 m/s and 2,330 ± 63 m/s, respectively, while the QS mode propagates with a velocity varying between 60 ± 2 m/s at 4 Hz and 287 ± 2 m/s at 60 Hz.

These velocities correspond to an apparent wavefield propagating under the array. Therefore, any variations of ice thickness or elastic properties along the array are averaged in the frequency-wave number diagram. It is possible to account for thickness variations [START_REF] Moreau | Measuring the wavenumber of guided modes in waveguides with linearly varying thickness[END_REF] and potentially for variations of elastic properties as well, using more advanced signal processing, but this is out of the scope of this paper and is left for future investigations. The variations of guided modes velocities between 1 and 5 March 2019 demonstrate the high sensitivity of the dispersion curves to ice thickness and elastic constants. This property of the guided wave modes are in favor of a very well-constrained inverse problem, which is actually already used extensively a the ultrasonic scale, in the fields of nondestructive testing and medical ultrasound to estimate bone mechanical resistance and thickness.

Dispersion Curves of Sea Ice From the Noise Correlation Function

Because it might not always be possible to isolate suitable icequakes for estimating sea ice properties, an alternative is to use the so-called noise correlation function (NCF). Historically, ambient seismic noise was used to qualify the contribution of all unwanted sources that corrupt the quality of conventional seismic recordings. For example, microseismic events or human activity generate broadband seismic noise, while wind or swell generate low-frequency noise. In passive seismology, the NCF is calculated by correlating the ambient seismic noise (or ambient seismic field) recorded between station pairs. It can be shown that it converges toward the impulse response, or Green's function, of the medium [START_REF] Sabra | Extracting time-domain Green's function estimates from ambient seismic noise[END_REF][START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]. Because the NCF integrates the contribution of all noise sources, impulsive or more stationary, it naturally results in more energetic signals with a broader frequency content than those from isolated events. The aim of this section is to demonstrate that one can achieve results from the NCF that are at least as good as those obtained from icequakes or active sources.

As shown in section 3.2, the ambient seismic field is isotropic at the scale of a few hours. This is ideal for the convergence of the NCF. The impulse response of the ice is known on 1 March 2019 from the impulsive active sources (jumps). Hence, for comparison we use the seismic noise recorded on the same day, from 0h00 to 12h00, at the stations of the linear array to the east (stations 401-101 to 401-104) and those of line 125 of the main array. Stations of the linear array are used as virtual sources, and those of line 125 are used as the receivers. Recordings were first split in 30-min time windows, and spectral whitening was applied in the 1-150 Hz frequency band. The cross-correlations were then calculated between each virtual source and all the receivers, resulting in 24 cross-correlations for each stations pair. To obtain the final NCF, the filter introduced in Moreau, Stehly et al. ( 2017) was applied to each set of cross-correlations, which were finally stacked. The purpose of this filter is to minimize the incoherent information in the singular vectors of the matrix of pre-stack correlations with the Wiener filter, for improved convergence of the NCF. An example of NCF from the vertical displacement component is shown in Figure 10b when the virtual source 401-104 is used.

For comparison, the Green's function between the same stations pairs was estimated by correlating signals recorded in the 5 s following the impulsive active sources, located exactly at the position of the virtual sources used for the NCF. This is shown in Figure 10a for the vertical displacement from the impulsive source at station 401-104. Note that the signals shown in Figures 10a and10b are filtered in the 1-10 Hz frequency band, for visualization purposes. The causal and a-causal part of the NCF are symmetric, a typical feature that is expected when the source distribution is isotropic. There is a very good agreement between the causal part of the NCF from virtual source 401-104 and the response from the active impulsive source at the same station, indicating convergence of the NCF. By reciprocity, the a-causal parts of the NCF converge toward the Green's function propagating in the opposite direction. Note the strong dispersion of the waveforms, typical of the QS mode. This was expected because the vertical displacement is dominated by this mode in this frequency band.

Based on these observations, and also on the fact that noise sources are broadband, with energy up to at least 500 Hz (Figure 5), it is expected that the dispersion branches can be measured from the NCF in a larger frequency band, because this processing integrates signal over a much longer duration than the icequakes. For example, Figure 10c shows the dispersion curves obtained from the NCF of the three displacement components. In terms of wavenumber values (or velocities), they are in prefect agreement with those calculated on the same day along line 125 of the main array, using the impulsive sources (Figure 9a). This proves that the velocities of the modes calculated using the NCF are the same as those inferred using the active impulsive sources. An interesting feature is that the dispersion curve of the QS mode vanishes at frequencies beyond 60 Hz when using the active sources, but it can be measured on a larger frequency band, up to 150 Hz when using the NCF. This indicates that the inverse problem may be better constrained with the NCF.

Inversion of the Dispersion Curves for Sea Ice Thickness, Young's Modulus, and Poisson's Ratio

From the approximation of the seismic wavefield in a floating ice layer by [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF], it is straightforward to determine Young's modulus and Poisson's ratio based on the modes velocity:

𝜈 = 1 -2 ( c QSH 0 c QS 0 ) 2
, and

E = 𝜌 × ( c QS 0 ) 2 ( 1 -𝜈 2 ) . ( 1
)
However, it should be noted that the use of equation ( 1) is possible only because the approximation of the QS 0 by the axial wave remains valid at our frequencies of interest, that is, between 0 and 150 Hz (Figure 1b).

10.1029/2019JC015709

We find values E = 3.9 ± 0.2 GPa and Poisson's ratio 𝜈 = 0.34 ± 0.02 on 1 March 2019 and E = 4.4 ± 0.2 GPa and Poisson's ratio 𝜈 = 0.32 ± 0.02 on 5 March 2019.

To estimate sea ice thickness, the dispersion curve of the QS mode is used. While the approximation in [START_REF] Stein | Inversion of pack ice elastic wave data to obtain ice physical properties[END_REF] holds up to a frequency-thickness of 300 Hz•m for the QS 0 mode, it holds only up to 50 Hz•m for the QS mode (Figure 1b). However, on 1 March 2019, sea ice thickness was greater than 0.5 m. Given that we measure the propagation of the QS mode up to 150 Hz (Figure 10c), another model has to be used. For meaningful results we use finite element (FE) modeling instead of the flexural wave approximation.

To this end, we modeled the seismic wavefield propagating in the sagittal plane of an ice layer floating on an infinite water column. FE simulations were made with the commercially available software COMSOL Multiphysics. For water, we used a density 𝜌 w = 1010 kg/m 3 and a speed of sound c = 1410 m/s. For sea ice, we used a density of 910 kg/m 3 [START_REF] Timco | A review of sea ice density[END_REF], and the elastics constants E and 𝜈 found previously on each day. Thickness was varied between 0.4 and 0.7 m with a 1-cm increment. The seismic source was defined as a vertical force distributed on a 20-cm-long zone at the free surface of the ice. We used three different positions of the source, all inline with an array of 44 receivers. They were located 92, 96, and 100 m away from the first receiver and sent a two-cycle tone burst with central frequency 75 Hz. The processing of the synthetic data was identical with that of the field data. However, because the simulation is two-dimensional, it was applied only to the vertical and longitudinal components of the wavefield.

Note the different relative amplitudes of the modes in the spectra from the synthetic and field data. This is essentially due to the fact that the displacement field used to measure the QS 0 mode in the synthetic data was taken at the center of the ice layer, where its energy is strongest. In the field data, we have only access to data measured at the surface of the ice, where the QS mode is dominant. In the spectra calculated from the icequakes, however, the difference in relative amplitude is not so significant. This is because the source mechanism is not the same when jumping on the ice or when cracking occurs. The former generates essentially a vertical motion, and thus, energy goes mainly to the QS mode, while the latter generates motion in the whole thickness and energy is distributed more evenly between the modes.

The dispersion curves calculated on 1 March 2019 from the impulsive sources or the NCF are in very good agreement with those for a simulated ice layer of thickness h = 54±3 cm, Young's Modulus E = 3.9 GPa, and Poisson's ratio 𝜈 = 0.34 (Figure 9c). The dispersion curves calculated on 5 March 2019 from the icequakes correspond to a simulated ice layer of thickness h = 58 ± 3 cm, Young's Modulus E = 4.4 GPa, and Poisson's ratio 𝜈 = 0.32 (Figure 9c).

These parameters are consistent with field observations and with values generally found in literature for sea ice:

-Thickness estimation on 1 March 2019 is consistent with that from the GPR acquisition along the east-west direction on the same day (52 cm). Drilling through the ice on 1 March 2019 revealed a slightly larger thickness of 62 cm. However, this does not mean that our thickness estimation from the dispersion curves is wrong, for two reasons. First, drilling was performed 20 m to the south of line 125 (see exact position in Figure 3b), and the spatial variability of the ice thickness was most likely around 15-20%. Second, as explained in section 3.4.2, the dispersion curves correspond to an average velocity along the propagation of the wavefield; hence, our thickness estimations correspond to an average value along the stations of line 125. -Poisson's ratio is typical for sea ice [START_REF] Mellor | Mechanical behavior of sea ice[END_REF][START_REF] Schulson | Creep and fracture of ice[END_REF]. -Young's moduli are consistent with values reported in [START_REF] Hunkins | Seismic studies of sea ice[END_REF] and [START_REF] Anderson | Preliminary results and review of sea ice elasticity and related studies[END_REF] for first-year ice, which is more porous and with a large brine content.

Discussion

The value of Young's modulus for sea ice is about half that of fresh-water ice [START_REF] Gammon | Elastic constants of artificial and natural ice samples by Brillouin spectroscopy[END_REF]. However, in the Van Mijen fjord, the salinity in March is generally around 6‰ [START_REF] Hoyland | Ice thickness, growth and salinity in Van Mijenfjorden Svalbard, Norway[END_REF] and the brine-liquid volume was shown to vary between 95‰ in 2010 and 30‰ in 2018 [START_REF] Karulin | Features of determining the ice flexural strength and the elastic modulus based on floating cantilever beam tests[END_REF]. At Vallunden Lake, the ice contained much brine and was quite porous, with water drawn to the surface from capillarity effect. Young's modulus is known to decrease when porosity and brine content increase [START_REF] Langleben | Elastic parameters of sea ice, Ice and snow; properties, processes, and applications[END_REF][START_REF] Mellor | Mechanical behavior of sea ice[END_REF]. The ice was also very brittle in the top 30 cm. Moreover, the ice in the field had a temperature of about -1 • C in the bottom 10 cm. For ice temperatures above -3 • C, Young's modulus 10.1029/2019JC015709 is significantly reduced [START_REF] Hunkins | Seismic studies of sea ice[END_REF]. These observation simply a gradient of mechanical properties through the ice layer. The snow layer covering the ice surface may also have an influence on the apparent Young's modulus [START_REF] Schulson | Creep and fracture of ice[END_REF]. Our onsite radar acquisitions show that the thickness of the snow was between 20 and 35 cm along row 125 (Figure 4b).

The combined effects of ice porosity, brine, and the snow cover affect the speed of the guided modes, which propagate with a reduced velocity in comparison with fresh-water ice. For example, the dispersion curves measured in [START_REF] Moreau | Monitoring ice thickness and elastic properties from the measurement of leaky guided waves: A laboratory experiment[END_REF] for a floating layer of pure ice with thickness of 1 cm, Young's modulus E = 9.4 GPa, and Poisson's ratio 𝜈 = 0.33 indicated a velocity of 3,215 m/s for the QS 0 mode. The velocity of the QS mode was 430 m/s for a frequency of 3,700 Hz that is for a frequency-thickness of 37 Hz•m.

In sea ice at Vallunden Lake, for the same frequency-thickness (i.e., at a frequency of 60 Hz and thickness of 62 cm), the QS mode propagates with a reduced velocity of 287 m/s. This indicates that the effective elastic properties estimated in sea ice with seismic methods are necessarily less than those of fresh-water ice measured at the ultrasonic scale, where mesoscale heterogeneities have no influence.

The differences between the elastic properties estimated on 1 and 5 March 2019 remain within the error bars. This is expected, since the ice should not undergo significant changes of properties within only 4 days. The estimations of the ice thickness remains within the error bars too, but we believe that the 4-cm increase in thickness is not related to uncertainties. Between 1 and 26 March 2019, drillings indicate a thickness increase of about 20 cm, which is not surprising given the cold weather conditions in the field, with an average temperature of about -25 • C.

Overall, the use of passive seismic data seems very promising for monitoring the thickness and mechanical properties of sea ice. However, the monitoring over the 27-day period is left for future works, since investigating NCF requires significant additional analyses that are beyond the scope of this paper. Future works will also investigate inversions based on a forward model that accounts for ice porosity, gradient of mechanical properties, and a snow layer. This should pave the way toward more systematic applications of passive methods for monitoring sea ice. The interest lies not only in the context of climate change at the global scale but also at a more local scale for monitoring the ice condition in places where fast ice is used for traveling.

Concluding Remarks

This paper presented a proof of concept to use passive seismic data for estimating sea ice thickness and properties. We show preliminary results from seismic data acquired on sea ice with an array of geophones deployed in the Van Mijen fjord near Sveagruva (Svalbard). The array recorded continuously the ambient seismic field for 27 consecutive days between 28 February and 26 March 2019, resulting in a sea ice seismic data set of unprecedented quality, both in terms of spatial and temporal sampling and also in terms of the acquisition duration. Based on passive acquisitions and appropriate signal processing, we show that broadband propagation of the seismic waves guided in sea ice can be extracted with the same level of accuracy compared to data from active seismic experiments. The dispersion curves of three fundamental guided modes are measured either with an active source or with passive seismic data, in the 1-150 Hz frequency range. The dispersion curves obtained on 1 March 2019 with active sources or the NCF show very good agreement with synthetic data in a layer of floating ice with thickness h = 0.54 m, Young's modulus E = 3.9 GPa, and Poisson's ratio 𝜈 = 0.34. Those obtained on 5 March 2019 from icequakes show very good agreement with synthetic data in layer of floating ice with a thickness of h = 0.58 m, Youngs modulus E = 4.4 GPa, and Poisson's ratio 𝜈 = 0.32. These values are consistent with field observations and reports in the literature. A major perspective of this work is the monitoring of these parameters for the whole duration of the experiment, for example, by producing a daily three-dimensional map of the ice properties. To this end, instead of using just one linear array, the four linear arrays will be used as virtual sources to calculate the NCF and recover wave propagation in all directions.

We are confident that the approach introduced in this paper can be complementary to the already existing methods for monitoring sea ice properties. Applications are for both landfast ice and ice floes. The origin of mechanical stress in ice floes is mainly due to interior deformation resulting from relative drift and thermal expansion of the ice, while at Vallunden Lake it is dominated by tidal forcing. Hence, there might not be as many icequakes in ice floes as those recorded in this experiment, but not as many are necessary. Moreover,

Figure 1 .

 1 Figure 1. (a) Displacement field generated by the QS mode (top) and the QS 0 mode (bottom) in a floating ice layer of thickness h = 50 cm at a frequency of 100 Hz, with Young's modulus E = 4.5 GPa, Poisson's ratio 𝜈 = 0.33, and density 𝜌 = 910 km/m 3 . (b) Wavenumber-thickness versus frequency-thickness dispersion curves of the ice layer. Black dots show the theoretical values and red lines the approximation of the QS and QS 0 modes by the flexural and axial waves, which is valid when the frequency-thickness product is less than 50 Hz•m.

Figure 2 .

 2 Figure 2. (a) Location of the seismic array in the Van Mijen fjord near Sveagruva (Svalbard), with (b) a zoom around the array area in Vallunden Lake, a part of the fjord that is surrounded by a moraine and connected to the fjord by a channel. The gray scale show land which altitude is less than 25 m. All land above 25 m is shown in white to emphasize the shore line.

Figure 3 .

 3 Figure 3. (a) The 247 stations of the array, including the main central array, and the four linear arrays to the north, east, south, and west. Arrows indicate the positions of the two GPR profiles. (b) Aerial view of the main array, with station numbering and a photo of one geophone lines. Geophone station numbers are given as intersection numbers of east-west and north-south profiles. The two crosses indicate stations that were installed but failed to record: 125-542 and 113-509. Station 133-517 was originally a 1C station but was replaced by a spare 3C station due to technical problems. Red circles are for 1C stations and blue circles for 3C stations. The large arrowheads indicate the positions of ice thickness measurements.

Figure 4 .

 4 Figure 4. Processed ground penetration radar profiles acquired on 28 February 2019 along the (a) north-south direction between stations 501-104 and 301-104 and (b) east-west direction between stations 401-104 and 201-104 (see Figure 3a for exact profiles locations).

  A series of active seismic experiments was conducted on 1 and 26 March 2019, at the times and positions indicated in Table A1, including vibrating and impulsive sources. Vibrating sources are 10, 15, and 60 s sweeps with a shaking device at frequencies varying between 1 and 450 Hz. Impulsive sources are person weight (80 kg) drops (jumps) from 1 m height onto the ice. Moreover, artificial seismic noise was created by driving five snowmobiles around the array for several minutes. The first driving sequence was on 1 March 2019 between 14:25 and 14:30, and the second on 26 March 2019, between 15:49 and 15:54 (UTC).

Figure 5 .

 5 Figure 5. Spectrogram of the signal recorded at the center of the main array by station 125-525 on 5 March 2019. Hundreds of broadband events are recorded every hour. The increase of activity between 8:30 and 17:30 corresponds to field work on the lake.

Figure 6 .

 6 Figure 6. Beamforming of one full day of ambient seismic noise recorded on 3 March 2019. The distribution of noise sources is isotropic, with a slightly lower amplitude to the south. Note the aliasing limit in the 15-17 Hz frequency band, where the QS mode has a wavelength of about 8 m, causing a duplication of the beamformer that is very close to overlapping itself near the slowness of 8 s/km.

Figure 7 .Figure 8 .

 78 Figure 7. Top: 8 hr of recording from the vertical component of station 125.525 on 5 March 2019 with a zoom between 1:30 and 2:30 a.m. Colored signal indicates events 1, 2, and 3 used to compute the dispersion curves of the ice (see Figure9). Bottom: Beamforming of the data from event 3 in the 15-17, 11-13, 7-9, and 3-5 Hz frequency bands.

Figure 9 .

 9 Figure 9. (a, b) Dispersion curves of the ice at Vallunden Lake obtained using the three displacement components measured along row 125, from (a) the four impulsive sources to the east of the main array on 1 March 2019 and (b) the propagation of the three events isolated on 5 March 2019 (Figure 7). (c) Dispersion curves obtained from the two displacement components of a two-dimensional finite element simulation of seismic waves propagating in a layer of floating ice with a thickness of h = 54 cm, Young's modulus E = 3.9 GPa, Poisson's ratio 𝜈 = 0.34, and density 𝜌 = 910 kg/m 3 . (d) Same as (c) with h = 58 cm, E = 4.4 GPa, Poisson's ratio 𝜈 = 0.32, and density 𝜌 = 910 kg/m 3 .

Figure 10 .

 10 Figure 10. (a, b) Correlograms between station 401-104 and the stations of line 125 in the array, filtered in the 1-10 Hz frequency band for better visualization purpose. (a) correlations from the ambient wavefield recorded for 12 hr on 1 March 2019. (b) Correlations from signals recorded during the 5 s following the impulsive source at station 401-104 on 1 March 2019. (c) Dispersion curves of the ice at Vallunden Lake along line 125, obtained from the NCF of the three displacement components with the stations of the east linear array as virtual sources.

Table 1

 1 Sea IceProperties and Thickness From Various Seismic Studies Between 1958 and 2019 SH 0 velocity (m/s) QS 0 velocity (m/s)

	E (GPa)	𝜈	Thickness (m)

  March 2019, 10 m east of station 141-549. A thickness of 62 cm was measured. It was also measured on 26 March next to stations 125-511 and 126-525, giving a thickness of 70 and 73 cm. Around position N77.8845 • , E16.776 • , about 170 m north of station 101-549, several measurements by another group of researchers gave ice thicknesses that varied between 75 and 80 cm between 4 and 16 March 2019. This indicates substantial spatial and temporal thickness variations of the ice cover at Vallunden Lake.

Acquisitions 2.2.1. Ice Drillings

Sea ice thickness was measured in the fjord on 1

Acknowledgments

Journal of Geophysical Research: Oceans 10.1029/2019JC015709 reverberation in ice floes may play a important role in this matter, since each reflection at the edge acts as a secondary noise source, which is desirable for the convergence of the NCF.

Key aspects for a transfer of this concept to systematic monitoring of sea ice with the ambient seismic wavefield are battery life and the determination of optimal array size and density for best spatial coverage and accuracy. The former is a technological issue. The geophones used in this experiment allow continuous recordings for 40 days at a sampling rate of 500 Hz. The next generation of geophones will be able to record continuously for several months and transmit data via a 4G connection. Moreover, continuous recording of the wavefield is not necessary to extract the dispersion curves from the NCF, because it converges after only a few hours of ambient field records are processed. The latter will be tackled by studying the sensitivity of the inverse problem to measurement errors, which are directly linked to spatial sampling and frequency bandwidth.

Measuring broadband propagation of the three guided modes from passive seismic data is an important step toward the assessment of a more optimal array design in terms of spatial sampling. For example, a linear array of 30 geophones with 10-to 15-m spacing might allow the extraction of the dispersion curves of the three modes on a bandwidth that is sufficient for estimating sea ice properties. This remains to be confirmed with future investigations, by focusing on a systematic search of suitable icequakes and the calculation of the NCF, as these may play a complementary role to achieve this goal. ). The full data set (total size = 3.2 To) is available upon request to the corresponding author. A subset can be accessed instantly via a direct download link. To request this link, please send an e-mail (resif-dc@ujf-grenoble.fr) with your affiliation and IP address. This subset contains the continuous recordings of all stations on 1, 3, and 5 March 2019, from which the results presented in this paper were obtained (with 3C stations data downsampled at 500 Hz). ISTerre is part of Labex OSUG@2020. This research was funded by the Agence Nationale de la Recherche (ANR, France) and by the Institut Polaire Français Paul-Emile Victor (IPEV). Most of the computations presented here were performed using the CIMENT infrastructure (https:// ciment.ujf-grenoble.fr), which is supported by the Rhone-Alpes region (Grant CPER07_13 CIRA: http://www. ci-ra.org) and France-Grille (http:// www.france-grilles.fr). We thank Pr. Aleksey Marchenko from the University Centre in Svalbard for the valuable advices and discussions. G. Hillers acknowledges support funding from the University of Helsinki Faculty of Science. All authors also want to thank the reviewers, the Editor, Dr. L. Padman, and Dr. G. Taylor for their comments that helped to improve the manuscript.