Thomas Paviet-Salomon
email: thomas.paviet-salomon@ensta-bretagne.org

Member, IEEE Julien Bonnel
email: jbonnel@whoi.edu

Clément Dorffer
email: clement.dorffer@ensta-bretagne.fr

IEEE Barbara Nicolas Member

Thierry Chonavel
email: thierry.chonavel@imt-atlantique.fr

Dag Tollefsen
email: dag.tollefsen@ffi.no

David P Knobles
email: dpknobles@kphysics.org

Preston S Wilson
email: pswilson@mail.utexas.edu

Angélique Drémeau
email: angelique.dremeau@ensta-bretagne.fr

Estimation of frequency-wavenumber diagrams using a physics-based grid-free compressed sensing method

Keywords: Underwater acoustics, compressed sensing, normal mode propagation, Seabed Characterization Experiment, SBCEX17

Shallow water propagation can be described using modal theory. For low frequency sources, propagated signals are composed of a few dispersive modes, each of them propagating with its own frequencydependent wavenumber. Modal estimation, and particularly wavenumber estimation, is of great interest in seabed characterization but classically requires a large and dense horizontal line array (HLA). The compressed sensing (CS) paradigm, which allows one to reduce the number of sensors, has been used to overcome this limitation. However, CS performance is directly linked to the discrete basis used in the process and is known to degrade with basis mismatch. To mitigate this issue, the current paper proposes a physics-based grid-free approach to perform wavenumber estimation using a HLA with a limited number of sensors and a single broadband source. The proposed method has three main features: it starts with a speed correction to prevent wavenumber aliasing (using water sound speed at the array location), it then embeds physical prior (the modal dispersion relation) at the core of the CS framework, followed by a CS grid-free approach. The performance of the method is quantified on simulated data using the Jaccard's distance. The method is then applied successfully on experimental data from the 2017 Seabed Characterization Experiment.

I. INTRODUCTION

In shallow-water environments (water depth D < 200 m), the acoustic propagation of a low-frequency (f < 100 Hz) signal is conveniently described by normal mode theory. Under this assumption, the signal is modeled by the sum of a small number of propagating modes (less than ten in this work). The modes are fully characterized by their frequency-dependent wavenumbers and amplitudes. These two modal features have been extensively used as input data for inverse problems, such as source localization [START_REF] Wilson | Matched mode localization[END_REF], [START_REF] Gall | Performance analysis of single-receiver matched-mode localization[END_REF] and/or environmental estimation [START_REF] Frisk | Shallow water waveguide characterization using the Hankel transform[END_REF], [START_REF] Ballard | Inversion for range-dependent water column sound speed profiles on the New Jersey shelf using a linearized perturbative method[END_REF]. As a result, modal estimation is an important topic for ocean acoustics. The present paper focuses on modal wavenumber estimation.

When range aperture is available, wavenumber estimation is similar to spectral analysis in the spatial dimension. Assuming that the signal is collected on a long horizontal line array (HLA) with a source in the endfire direction, the easiest way to estimate the wavenumber spectrum is to compute a spatial Fourier transform (SFT) in the horizontal dimension. Further, if the source is broadband, one can also compute a second Fourier transform in the temporal dimension (TFT) to obtain a frequency-wavenumber (f-k) diagram, which fully characterizes modal dispersion [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF], [START_REF] Nicolas | Geoacoustical parameters estimation with impulsive and boat-noise sources[END_REF].

The main drawback of this simple procedure is that it suffers from the traditional limitations of the Fourier-based spectral analysis. In our context, a large horizontal aperture and a large number of sensors (i.e. a long and dense HLA) are required to properly resolve the modal wavenumbers. To mitigate this issue, the SFT can be replaced by more advanced spectral estimation methods. Pioneer studies successfully replaced the SFT with an auto-regressive (AR) estimator [START_REF] Becker | Evaluation of an autoregressive spectral estimator for modal wave number estimation in range-dependent shallow water waveguides[END_REF], [START_REF] Ballard | Geoacoustic inversion for the New Jersey shelf: 3-D sediment model[END_REF]. Since the propagation in shallow-water environments is described by a small number of modes, considering sparse models constitutes an interesting alternative. These models are in particular involved in the compressed sensing (CS) paradigm [START_REF] Candès | An introduction to compressive sampling[END_REF] which, as it will be emphasized in the following, is of particular interest here, but has been more broadly exploited in underwater acoustics (e.g. [START_REF] Chapman | Deconvolution of marine seismic data using the l1 norm[END_REF], [START_REF] Gerstoft | Single and multiple snapshot compressive beamforming[END_REF]). CS has first been used to estimate modal wavenumber in ocean acoustics by Le Courtois et al. [START_REF] Courtois | Compressed sensing for wideband wavenumber tracking in dispersive shallow water[END_REF]. Interestingly, CS has also been used to estimate modal wavenumber in other scientific fields, e.g. structural health monitoring using Lamb waves [START_REF] Harley | Dispersion curve recovery with orthogonal matching pursuit[END_REF], [START_REF] Sabeti | Sparse wavenumber recovery and prediction of anisotropic guided waves in composites: a comparative study[END_REF]. A special issue of the Journal of The Acoustical Society of America is dedicated to CS in acoustics, and the editorial introduction [START_REF] Gerstoft | Introduction to special issue on compressive sensing in acoustics[END_REF] includes a comprehensive review of CS.

Traditionally, CS algorithms expand the signal on the elements of a finite basis, i.e. a discrete grid.

When the components of the signal do not match the grid, basis mismatch occurs, which degrades the CS performances [START_REF] Duarte | Spectral compressive sensing[END_REF]. Off-grid CS approaches have been proposed to mitigate this issue [START_REF] Knudson | Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit[END_REF], [START_REF] Denoyelle | Theoretical and Numerical Analysis of Super-Resolution Without Grid[END_REF].

They have been notably applied to underwater acoustic problems, such as plane wave beamforming [START_REF] Xenaki | Grid-free compressive beamforming[END_REF] or modal estimation using a vertical line array (VLA) and a source at multiple ranges [START_REF] Park | Grid-free compressive mode extraction[END_REF]. Later, [START_REF] Denoyelle | Theoretical and Numerical Analysis of Super-Resolution Without Grid[END_REF] DRAFT Paviet-Salomon et al. [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF] proposed to use grid-free CS to perform modal wavenumber estimation using a single broadband source and a small number of sensors. As in Ref. [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF], they further embedded a physical assumption, the dispersion relation, within the CS framework to relate wavenumbers from one (temporal) frequency to the next. The method was successfully validated on data simulated in a Pekeris waveguide [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF].

In the present article, we extend the work of Paviet-Salomon et al. [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF]. The method is improved by implementing a frequency-dependent shift of the wavenumber spectra. This processing, which effectively creates an equivalent baseband version of the spatial signal thanks to a speed correction, is traditionally used for seismic applications [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF], [START_REF] Nicolas | Geoacoustical parameters estimation with impulsive and boat-noise sources[END_REF]. Here, it prevents wavenumber aliasing and thus allows wavenumber estimation over a relatively wide frequency band. On the other hand it modifies the dispersion relation that is used at the core of the CS framework.

In the present work, the performances of the proposed method is benchmarked on realistic simulations, and shown to be superior to the state of the art. The method is also applied on experimental marine data collected during the 2017 Seabed Characterization Experiment (SBCEX17) [START_REF] Wilson | Guest editorial an overview of the seabed characterization experiment[END_REF]. This article considers a combustive sound source (CSS) signal recorded on a HLA with an aperture of 1 km. The proposed method allows the estimation of the wavenumbers of the four first modes from 10 Hz to 100 Hz using as few as 10 sensors. These experimental results, obtained on SBCEX17 data collected on the New England Mud Patch, illustrate the method's robustness to complex environments with vertical stratification. Indeed, the New England Mud Patch features a complex layered seafloor. It notably has a layer of mud whose upper part is believed to be slower than water [START_REF] Belcourt | Depth-dependent geoacoustic inferences with dispersion at the new england mud patch via reflection coefficient inversion[END_REF], [START_REF] Ballard | In situ measurements of compressional wave speed during gravity coring operations in the new england mud patch[END_REF], which clearly impacts the propagating modes [START_REF] Bonnel | Trans-dimensional geoacoustic inversion using prior information on range-dependent seabed layering[END_REF].

The remainder of the paper is organized as follows. Section II briefly introduces modal propagation, and focuses on the dispersion relation which will be at the core of the proposed method. Section III reviews the CS framework. It then describes on-grid and off-grid methods. Section IV presents the proposed wavenumber estimation procedure. First, Sec. IV-A reviews speed correction, a process that prevents wavenumber aliasing. Then, Sec. IV-B details the physics-based grid-free wavenumber tracking algorithm. Lastly, Sec IV-C introduces the Jaccard distance, a suitable metric to assess the performances of the method. Applications are presented in Sec. V, which covers a simulated benchmark and comparison with the state-of-the-art, as well as experimental results obtained using data from the 2017 Seabed Characterization Experiment. Section VI summarizes and concludes the article.

May 27, 2021 DRAFT

II. ACOUSTIC PROPAGATION IN DISPERSIVE SHALLOW WATER ENVIRONMENTS

A. Received signal

This paper focuses on low-frequency (f < 100 Hz) acoustic propagation in shallow water (D < 200 m).

In this context, as stated in Sec. I, the propagation is described by the modal theory. Consider a frequencydomain signal s(f ) emitted by a source at depth at z s . The received pressure on a sensor located at the distance r and depth z can be written as [27,Chap 5]:

y(f, r) Q s(f ) √ r M (f ) m=1 A m (f, z s , z)e -jrkrm(f ) + n(f, r), (1)
where Q is a constant factor, M (f ) is the number of propagating modes at frequency f , k rm (f ) is the horizontal wavenumber of the m th mode, and A m (f, z s , z) is its amplitude. The quantity n(f, r) stands for the TFT of the measurement noise. Modal attenuation is usually included in Eq. ( 1) as the imaginary part of k rm . Here, the associated term e -r [krm(f )] is included in A m , so that k rm can be considered a real number.

We consider a HLA and a source in the endfire direction (i.e. the source is aligned with the array).

The geometrical attenuation factor 1/ √ r in Eq. ( 1) can be compensated for if the source position is known. As a result, the only significant range-dependence in Eq. ( 1) is within the phase of the modes.

Note that this is still a fair assumption even if the source position is unknown, as long as it is in the endfire direction. Indeed, most of the range variability in Eq. ( 1) is driven by the complex exponential e -jrkrm(f ) , and the 1/ √ r can be considered as constant along the HLA, provided that the source/array distance is large enough. As a result y(f, k), the SFT of y(f, r), provides a direct measurement of the wavenumber spectrum (i.e. with peaks at k rm (f )). For a broadband source, y(f, k) is called a f-k diagram [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF]. Figure 1

B. Dispersion relation

In any environment, the horizontal wavenumbers k rm are linked to their vertical counterparts k zm by the dispersion relation. At a given frequency f , this relation is

2πf c 2 = k rm (f ) 2 + k zm (f ) 2 , (2)
May 27, 2021 DRAFT where c is the speed of the sound. Note that in theory, both k zm (f ) and c are depth dependent. This is ignored here, as was done in previous studies [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF], [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF], [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF] that use the dispersion relation to guide modal estimation. Indeed, the potential impact of this assumption is small enough that it can be embedded into noise and/or other uncertainties.

As suggested by Le Courtois and Bonnel [START_REF] Courtois | Compressed sensing for wideband wavenumber tracking in dispersive shallow water[END_REF], one can discretize the frequency axis (with f = ν∆ f , ν ∈ N and ∆ f the size of a frequency bin) and relate the wavenumbers attached to two successive frequency indices using

k rm [ν + 1] 2 = k rm [ν] 2 + (2ν + 1) 2π∆ f c 2 + [ν], (3)
where

k rm [ν] = k rm (ν∆ f ) and [ν] = k zm [ν + 1] 2 -k zm [ν] 2 .
In shallow-water environments, the vertical wavenumbers k zm weakly depend on the frequency [27,Chap. 5]. As a result, the quantity is smaller than the other terms of the equation and can be neglected.

This hypothesis has successfully been used in previous studies that took advantage of Eq. ( 3) at the core of modal estimation scheme [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF], [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF], [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF]. In the present paper, Eq. (3) will be used as a physical prior to enhance wavenumber recovery using an off-grid CS algorithm.

III. COMPRESSED SENSING

A. Sparse representation

Using the discretized framework presented in Sec. II-B, we further assume that the received signal has been measured on a HLA with L elements. The received signal y(f, r) is now denoted using y ν , a May 27, 2021 DRAFT column vector of size L. Assuming that the horizontal wavenumber space is discretized into a grid of size N , Eq. ( 1) can be expressed as

y ν = Da ν + n ν , (4)
where D ∈ C L×N is a dictionary made up of (spatial) Fourier component, a ν ∈ C N is the wavenumber spectrum at the (temporal) frequency ν (i.e. the transposed version of the ν th line of the f-k diagram of interest) and n ν ∈ C L is the additive noise along the HLA. The (l, n) th element of D is defined as

d nl = e -jrlκn
, where r l is the distance between the l th sensor and the source, and κ n is the n th element of the wavenumber search grid. Note that Eq. ( 4) requires knowledge about the HLA configuration (i.e. sensor spacing), but does not require the source/array distance to be known, provided that the source is at the endfire direction. In this case, r l can be referenced to the first sensor of the HLA (i.e.

r 1 = 0). The resulting solution is similar to the one obtained using the absolute source/array distance, up to a multiplicative phase shift. Note that the n th column of D will be noted d κn hereinafter, i.e.

d κn = [1, . . . , e -jrLκn ] T .
The discrete wavenumber spectrum a ν can be estimated through the SFT of y ν , which is equivalent to a least-square estimation process. However, since the number of propagating modes is small, most of the elements of a ν are null. As a result, the use of sparse recovery (SR) is well adapted to estimate a ν [START_REF] Courtois | Compressed sensing for wideband wavenumber tracking in dispersive shallow water[END_REF], [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF]. Within this framework, D is seen as an overcomplete dictionary (i.e. L N ) while a ν contains many zero entries. The corresponding SR problem can be expressed as

âν = argmin aν y ν -Da ν 2 2 , subject to a ν 0 ≤ M ν , (5)
with ||a ν || 0 the l 0 pseudonorm of a ν which simply counts the numbers of non-zero entries of a ν and M ν the maximum number of modes expected to propagate at frequency ν.

B. On-grid CS algorithm

Eq. ( 5) is a combinatorial problem, and many heuristic methods have been developed to solve it.

They can roughly be divided into three families. A first group of methods replaces the l 0 -norm by a l p -norm (with 0 < p ≤ 1). This leads to a relaxed problem which can be solved efficiently by standard optimization procedures. A well-known approach is the Basis Pursuit (BP) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. A second group of methods are greedy algorithms which build up the sparse vector a ν from successive greedy decisions.

One of the most popular versions of such algorithms is Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. The last group of methods are Bayesian algorithms that express the sparse representation as the solution of a Bayesian inference problem. Different Bayesian methods consider different prior models, estimation problems and statistical tools to solve them [START_REF] Olshausen | Sparse coding with an overcomplete basis set: A strategy employed by v1?[END_REF], [START_REF] Fevotte | A bayesian approach for blind separation of sparse sources[END_REF].

Various on-grid CS methods have been used to estimate wavenumbers in ocean acoustics, including an iterative reweighted least squares method [START_REF] Courtois | Compressed sensing for wideband wavenumber tracking in dispersive shallow water[END_REF] and two Bayesian approaches [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF], [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF]. Interestingly, the two Bayesian approaches [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF], [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF] use the dispersion relation Eq. ( 3) to connect horizontal wavenumbers from one frequency to the next. However, Ref. [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF] assumes that the number of modes is constant over the bandwidth of interest and that the [ν] term in Eq. ( 3) can be fully neglected. These two assumptions restrict the modal recovery to a limited bandwith (in Ref. [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF], the provided examples were limited to a 20 Hz bandwidth, which kept the number of modes constant). On the other hand, Ref. [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF] also uses

Eq. ( 3) with [ν] = 0 to predict wavenumber at frequency ν + 1, but it introduces freedom around the predictions to account for [ν] = 0. The associated drawback is that the method includes several parameters that must be tuned by hand.

In the present article, we consider an off-grid CS method that also relies on the dispersion relation Eq. ( 3). We will later see that the proposed method enables tracking the wavenumbers over a relatively wide bandwidth, while the number of modes does not need to be known a priori. Before that, traditional off-grid CS is briefly presented.

C. Off-grid CS algorithm

On-grid CS methods have known limitations that notably occur when the non-zero coefficients of a ν do not match the grid points. This issue, called basis-mismatch [START_REF] Chi | Sensitivity to basis mismatch in compressed sensing[END_REF], [START_REF] Duarte | Spectral compressive sensing[END_REF], can be mitigated by using a very fine grid. However, this raises questions about the coherence of the grid (which in turns impacts CS performance) and also may induce numerical instabilities.

These concerns led to the development of grid-free setting methods. Practically, a grid-free version of the relaxed version of Eq. ( 5) can be obtained by replacing the l 1 norm (only valid in a finite dimensional setting) by the total variation norm [START_REF] Candès | Towards a mathematical theory of super-resolution[END_REF]. Although the underlying theory is complex, offthe-shelf toolboxes such as the CVX software (ConVeX: a library for convex optimization) [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF] are available to solve the problem. This approach was notably used to estimate wavenumbers (and modal depth functions) considering a VLA and a source at multiple ranges [START_REF] Park | Grid-free compressive mode extraction[END_REF]. However, the main drawback of this procedure is its high computational cost.

Simultaneously, a continuous version of OMP has been proposed in Ref. [START_REF] Knudson | Inferring sparse representations of continuous signals with continuous orthogonal matching pursuit[END_REF] while a continuous version of BP has been proposed in Ref. [START_REF] Ekanadham | Recovery of sparse translation-invariant signals with continuous basis pursuit[END_REF]. The approach is based on the idea of interpolating the dictionary components in between existing grid points. However, an expensive computational cost is still associated with those methods. Rather, we will use a traditional OMP method, coupled with gradient descent performed after each OMP step [START_REF] Chen | Sparse off-grid doa estimation method with unknown mutual coupling effect[END_REF]. The main advantage of this method is to obtain a continuous behavior while being simple, it barely increases the computational cost of the considered SR method. It May 27, 2021 DRAFT was previously used in our context (estimation of modal wavenumbers with a HLA) on data simulated in a Pekeris waveguide [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF].

IV. WAVENUMBER ESTIMATION

This section describes the proposed method to estimate modal wavenumbers. A modified dispersion relation that will be used within the CS framework is first presented.

A. Dispersion relation after speed correction

In any waveguide, the modal wavenumbers are known to be contained within the interval

[ 2πf cmax , 2πf cmin ],
where c max and c min are the maximal and minimal sound speed in the environment [27,Chap. 5]. At any frequency, the wavenumber support of the signal is limited to a relatively narrow (wavenumber) band:

the spatial signal is said to be a bandpass signal. A small sensor spacing is thus required to prevent wavenumber aliasing. An example of an aliased f-k diagram is presented in Fig. 2a. The simulated environment is the same as in Fig. 1, except that the frequency now goes up to 200 Hz. One can see that above 60 Hz, wavenumbers are aliased. This is because for f > 60 Hz, the maximal wavenumber value 2πf cmin is above the (spatial) Nyquist frequency k s = 1 2∆r , with ∆r the (spatial) sampling rate (i.e.

the sensor spacing).

Aliasing can be prevented by reducing the spacing between sensors, but this is not always possible at sea. Since the wavenumber spectrum is a bandpass signal, aliasing can also be prevented by translating the modal spectrum toward smaller wavenumber (i.e. baseband processing). Mathematically, this is done by time-shifting the received data using c min sound speed, which in the frequency domain corresponds to

ỹ(f, r) = e jr 2πf c min y(f, r), (6)
with ỹ(f, r) the shifted signal. The wavenumber support of ỹ becomes [ 2πf cmax -2πf cmin , 0], which drastically reduces the Nyquist frequency. This process is common in geophysics [START_REF] Yilmaz | Seismic data analysis: Processing, inversion, and interpretation of seismic data[END_REF], [START_REF] Freire | Application of singular value decomposition to vertical seismic profiling[END_REF], [START_REF] Duncan | Median filter behaviour with seismic data 1[END_REF] and has also found applications in ocean acoustics [START_REF] Nicolas | Geoacoustical parameters estimation with impulsive and boat-noise sources[END_REF], [START_REF] Nardin | 1-200 Hz wave propagation in shallow water[END_REF]. It will be called speed correction hereafter. Speed correction is applied on the simulated Pekeris example, and the result is presented in Fig. 2b. One can see that wavenumbers can now be recovered without aliasing up to 200 Hz.

Although the speed correction process is simple, it modifies the dispersion relation. The original wavenumbers k rm [ν] and their shifted versions krm [ν] are linked through with γ = 2π∆f cmin . Using Eq. ( 7) into Eq. ( 3), the discretized dispersion relation becomes

k rm [ν] = krm [ν] + γν, (7)
( krm [ν + 1] + γ(ν + 1)) 2 = ( krm [ν] + γν) 2 + (2ν + 1)γ 2 + [ν]. (8)
Assuming that [ν] = 0, the wavenumber at frequency ν + 1 can be predicted using

kpred rm [ν + 1] = -γ(ν + 1) + ( krm [ν] + γν) 2 + (2ν + 1)γ 2 . ( 9
)
Note that in the proposed method, the (incorrect) assumption [ν] = 0 will be mitigated by looking for wavenumbers in a small interval centered around the predicted value kpred rm [ν + 1].

B. Tracking modes

As stated in the introduction, the tracking method proposed in this article is largely similar to the one in Ref. [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF]. The only difference is that we now work after speed correction. As a result, the physical prior that is embedded into the grid-free CS algorithm is now Eq. ( 8), instead of Eq. ( 3). This enables tracking wavenumbers in configurations where sensor spacing is large and creates wavenumber aliasing.

In such situation, the method from Ref. [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF] cannot be used. Note that if sensor spacing is small enough May 27, 2021 DRAFT to prevent aliasing without speed correction, then the proposed method gives result similar to the method in Ref. [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF]. The proposed method is summarized below.

The proposed procedure starts at ν = 1 with a traditional OMP step, evaluating the most correlated atom of the dictionary D with the signal y ν . This step is then completed by a detection operation, which compares the resulting correlation to a given threshold T 0 . If selected, a mode is considered to propagate and the dispersion relation ( 8) is used to predict the interval in which the wavenumber of this mode is likely to be at the next frequency. The current estimate, as well as the predicted wavenumber at the next frequency, are then refined by a gradient-descent step, respectively based on the selected atom and within the predicted interval. For the predicted wavenumber of mode m at the next frequency ν + 1, a detection threshold T m,ν+1 is again applied on its estimated amplitude to prevent the propagation of false alarms.

All these steps are then repeated at the next frequency in order to detect possible new modes and to propagate those already predicted to higher frequencies.

The overall procedure is summarized in the Algorithm 1. An important parameter is ξ in Eq. ( 11): it gives the algorithm freedom to look for krm [ν + 1] around the predicted value kpred rm [ν + 1]. In a given context, the value of ξ is empirically determined by simulations to maximize the method's performances.

Algorithm 1 Physics-based grid-free Orthogonal Matching Pursuit for data with speed correction.

0. Initialization : ∀ν ∈{1, .., F }, rν =yν , Sν =∅, I 0,ν =∅, M =0.

For ν = 1 : F 1. Find new propagating wavenumber

κn = argmax κn | rν , dκn |, ( 10
)
where dκn is the n-th column of D∪ m∈{0,...,M } Im,ν made up of Fourier atoms not in ∪ m∈{0,...,M } Im,ν .

If | rν , dκn | ≥ T 0 -Set M = M + 1.
-Apply gradient-descent algorithm to refine previous estimate and get krm

[ν]. Set Sν = Sν ∪ krm[ν].
-Compute corresponding coefficients and update residual : rν = yν -D Sν âν,Sν with âν,Sν = D + Sν yν , with D + Sν the Moore-Penrose pseudo-inverse matrix of D Sν , dictionary made up of Fourier atoms specified by Sν .

If M > 1, propagate existing modes

For m = 1 : M -1 -Apply gradient-descent on interval Im,ν and get krm

[ν]. If | rν , d krm[ν] | ≥ Tm,ν , set Sν = Sν ∪ krm[ν].
-Compute corresponding coefficients and update residual : rν = yν -D Sν âν,Sν with âν,Sν = D + Sν yν .

Predict propagating intervals for next frequency

For all m ∈ {1, . . . , M }, define

I m,ν+1 = kpred rm [ν + 1] -ξ; kpred rm [ν + 1] + ξ . (11)
Detection thresholds may depend on false alarm probabilities. More particularly, considering a Gaussian complex circular noise assumption with variance σ 2 , we can define T 0 as

T 0 = σ -2 log β 0 (12)
May 27, 2021 DRAFT where β 0 is a given false alarm probability considered for the detection of a new mode. We adopt a similar definition for T m,ν , but we add further physical considerations. If a mode m is detected first at frequency ν c,m with threshold T 0 , we know that the mode will still be present at higher frequencies, suggesting thus to decrease the threshold for increasing frequencies. This can be done by considering:

T m,ν = T 0 - ν-νc,m i=1 1 -T∞ T0 2 -T∞ T0 i T 0 , (13)
where T ∞ = σ √ -2 log β ∞ stands for a desired limit threshold depending on an asymptotic false alarm probability β ∞ under the same Gaussian complex circular noise assumption. More details on the adopted strategy with regard to these thresholds can be found in [START_REF] Paviet-Salomon | Dispersive grid-free orthogonal matching pursuit for modal estimation in ocean acoustics[END_REF].

C. Performance of the algorithm : the Jaccard distance

An interesting question arises in quantifying the performance of the proposed method. Since the mode number is not known a priori, it is important to account for mismatch between the estimated number of modes, the true number of modes, and the fact that this number may vary with frequency. As a result, the problem can be seen as the detection and tracking of an unknown number of targets. A relevant metric in this context is the Jaccard distance [START_REF] Denoyelle | Theoretical and Numerical Analysis of Super-Resolution Without Grid[END_REF], [START_REF] Jaccard | Étude comparative de la distribution florale dans une portion des alpes et des jura[END_REF], [START_REF] Gower | Similarity, dissimilarity, and distance, measures of[END_REF].

The Jaccard distance D J is built using traditional detection theory features: the true positive (TP), false negative (FN) and false positive (FP) rates. Here, the quantities TP, FN and FP are evaluated at every frequency of interest using an acceptance radius R J . To explain the definition of TP, FN and FP, let us assume that we have a theoretical wavenumber k theo rm [ν] and an estimated wavenumber k r [ν]:

• if [k theo rm [ν] -R J , k theo rm [ν] + R J ] contains a single estimated mode k r [ν], then k r [ν] is a TP; • if [k theo rm [ν]-R J , k theo rm [ν]+R J ] contains several estimated modes, then the closest k r [ν] from k theo rm [ν]
is a TP while the others are FP (it is implicitely assumed that R J is smaller than half of the smallest distance between all theoretical modes, i.e. the acceptance radii do not intercept);

• if there is no estimated mode k r [ν] such that |k theo rm [ν] -k r [ν]| ≤ R J then k theo rm [ν] is a FN.
Mathematically, the Jaccard's index J I is first defined

J I = TP TP + FP + FN ( 14
)
and leads to the definition of the Jaccard distance

D J = 1 -J I ( 15
)
which is such that 0 < D J < 1. Note that D J = 0 means that TP = 1 and FP = FN = 0. On the other hand D J = 1 means that T P = 0, while D J 0 means that FP and/or FN are much larger than TP.

May 27, 2021 DRAFT The relevance of the Jaccard distance is illustrated using data simulated in our Pekeris waveguide example. For the sake of simplicity, wavenumbers are estimated using a classical OMP method. OMP is an iterative algorithm, selecting a new wavenumber at each iteration. Two stopping criteria are classically considered for this algorithm: it can stop either when the estimated reconstruction error drops below a given threshold or when the number of non-zero coefficients achieves a given number. This latter criterion is very convenient in our context, because the number of iterations is equivalent to the number of modes to be recovered. Although practical, this requires accurate prior information about the number of modes at each frequency. On the other hand, the stopping criterion based on the energy of the residual is less demanding in terms of prior information. This criterion is similar to what is used in the proposed method.

It will thus be used to illustrate the interest of the Jaccard distance.

Estimation performance is illustrated in Fig. 3. The TP rate increases when the number of sensors increases (because more data makes the estimation easier), and increases when the frequency decreases (because modal density decreases which also facilitates the estimation). On the other hand, the FP rate is high at very low frequencies (because the number of modes is very low, and thus a single false positive drastically increases the FP rate). The FP rate also tends to increase with the number of sensors. This is due to the chosen stopping criterion which implies an increasing number of iterations in OMP, and thus increases the likeliness to make false positive mistakes.

May 27, 2021 DRAFT Overall, the Jaccard distance D J is consistent with the behavior of TP and FP. As an example, D J is high (estimation performance is poor) when TP is small (e.g. for small number of sensors) and/or when FP is high (e.g. at very low frequencies). Also, when the sensor number is large enough, D J is relatively constant, because TP and FP tend to compensate for each other. The Jaccard distance is thus a good metric to summarize global estimation performance and will be used in the following.

V. APPLICATION

In this section, the performance of the method is assessed on realistic simulations. The method is also applied on experimental marine data collected during the 2017 Seabed Characterization Experiment (SBCEX17). The considered applications focuses on low-frequency data (f < 100 Hz). In our context, this makes the number of modes relatively small (5 or less), which allows for modal estimation with a small number of sensors. This choice is consistent (and allows comparison) with previous modal estimation studies based on sparsity [START_REF] Courtois | Compressed sensing for wideband wavenumber tracking in dispersive shallow water[END_REF], [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF], [START_REF] Niu | Block sparse bayesian learning for broadband mode extraction in shallow water from a vertical array[END_REF].

This section starts with a description of SBCEX17, whose context will be used to define the simulated environment.

A. SBCEX17

The SBCEX17 was dedicated to the understanding of sound propagation in fine-grained sediment. It was a multi-year, multi-institutional and multi-disciplinary effort which took place on the New England Mud Patch (NEMP). The NEMP is located about 100 km south of Cape Cod. The area is characterized by a relatively flat bathymetry (water depth D ∼ 70-75 m) and a thick upper sediment layer of mud.

A preliminary environmental survey was conducted in 2015. It notably included an intensive coring effort, as well as a high-resolution seismic survey. The main experiment took place in March-April 2017.

It involved three research vessels, and many acoustic sources and receivers were deployed, covering frequencies from ∼ 10 Hz to 10 kHz. A previous special issue of the IEEE Journal of Oceanic Engineering is dedicated to SBCEX17. Its editorial introduction [START_REF] Wilson | Guest editorial an overview of the seabed characterization experiment[END_REF] provides a succint overview of SBCEX17.

B. Experimental context

The SBCEX data was collected on a long HLA (aperture ∼ 1 km). The considered source is a

Combustive Sound Source deployed at 18:38 UTC on March 18 2017 at CSS station 29, located at (40.4983N, 70.5842W). The source signal is a powerful broadband impulse followed by several secondary impulses called bubble pulses [START_REF] Mcneese | An impulsive source with variable output and stable bandwidth for underwater acoustic experiments[END_REF]. Note that here, the specific source waveform does not matter as long as the source signal to noise ratio (SNR) is good enough in the frequency band of interest. The spacing (from 0.8 m to 72 m), with a symmetry around the array center and a higher density in the middle of the array. The specific source/array configuration has been chosen so that the source is in the endfire direction, and the distance between the source and the HLA is about 10 km. Note that in this context, the range variability along the HLA induced by the 1/ √ r term in Eq. ( 1) is negligible and can be ignored.

C. Simulation framework

First, the performance of the method is assessed using simulated data that mimics the experimental context. The environment is modeled using results from an inversion study performed using the HLA data [START_REF] Tollefsen | Ship-of-opportunity noise inversions for geoacoustic profiles of a layered mud-sand seabed[END_REF]. The environment is modeled as follows:

• water column: depth D = 65 m, sound speed gradient from 1468 m/s at the top of the layer to 1469 m/s at the bottom;

• sediment layer: thickness h = 5 m, sound speed c sed = 1500 m/s, density ρ sed = 1.6 g/cm 3 , attenuation α sed = 0.1 dB/λ;

• basement: c bas = 1700 m/s, density ρ bas = 2.0 g/cm 3 , attenuation α sed = 0.2 dB/λ.

Acoustic propagation is simulated in this environment over 0-100 Hz using the normal mode code ORCA [START_REF] Westwood | A normal mode model for acousto-elastic ocean environments[END_REF] with a frequency bin size of XX Hz, which led to simulation of YY frequencies. The environmental impulse response along the array is simulated in the range-frequency domain for a source at r = 10 km in the endfire direction. The array configuration mimics the experimental geometry. A bidimensional (2D) Gaussian white noise is added to the impulse response in the range-frequency domain, and the signal to noise ratio (SNR) is evaluated as the power of the (2D) range-frequency impulse response divided by the power of the (2D) range-frequency noise. Speed correction (see Sec. IV-A) is then applied with c min = 1468 m/s. Finally, f-k diagrams are computed using four different methods:

1) OMP: a traditional OMP algorithm;

2) SoBaP: a soft Bayesian Pursuit method that uses the dispersion relation to relate wavenumbers from one frequency to the next [START_REF] Drémeau | Reconstruction of dispersion curves in the frequency-wavenumber domain using compressed sensing on a random array[END_REF];

3) COMP: a grid-free OMP algorithm that uses gradient descent (COMP: Continuous OMP); 4) CPOMP: the method proposed in this article (CPOMP: Continuous and Physics-based OMP, i.e.

which uses the dispersion relation).

May 27, 2021 DRAFT

The methods OMP, SoBaP and COMP have been chosen as representative of the state-of-the-art. Note that OMP and COMP have been implemented here using the mode number as a stopping criterion (i.e. at frequency f , each algorithm looks for M (f ) wavenumbers). This gives them a strong advantage over the proposed method, which automatically determines the mode number. For CPOMP, we set β 0 = 0.001, β ∞ = 0.5 and ξ = 10 -3 .

Beside, all the methods considered here rely on the construction of the dictionary D (see Sec.III-A).

To construct this dictionary, one needs to sample the wavenumber axis with N elements. This is done by choosing a grid size and a maximum wavenumber value of interest. The grid size is defined as one fifth of the smallest value between 2 propagating modes over the frequency band of interest. Further, the maximum wavenumber value is computed considering the Nyquist theorem applied to the maximum space between two consecutive sensor (here ∆ k = 5.10 -4 rad/m and κ n = -0.09 rad/m).

D. Performance study

A first set of simulations studies the impact of the number of sensors on the f-k diagram estimation for a given SNR (12 dB). For a given number of sensors, sensors are randomly selected, except for the two sensors at the extremities of the HLA which are always selected (these two sensors are on both ends of the HLA so that total aperture of the array is always the same). The obtained performance is shown in Fig. 4 for 2 different frequency bands: 0-50 Hz and 50-100 Hz. Clearly, performance increases with number of sensors for all the methods, and the proposed method (CPOMP) outperforms the state of the art. Performance gain is minimal in the 0-50 Hz band, because most methods (OMP, COMP, CPOMP) provide satisfactory results. However, at higher frequencies where more modes are propagating, the gains of the proposed method becomes evident, for all number of sensors. As an example, with 40 sensors, CPOMP has D J < 0.1. This corresponds to an excellent f-k diagram recovery, as is illustrated at the end of this subsection. Interestingly, SoBaP which is supposed to be the most advanced of the state-of-the-art methods is the one with the worst performance. This is because SoBaP has been developed for data without speed correction. Although the SoBaP dispersion relation has been modified here to take into account the corrected dispersion relation, this is clearly not enough to correctly track wavenumbers. It is likely possible to further modify SoBaP to regain performance, but this is beyond the present scope.

A second set of simulations studies the impact of SNR for a given number of sensors [START_REF] Fevotte | A bayesian approach for blind separation of sparse sources[END_REF], and for The Jaccard distance D J study fully characterizes the performance of the four methods. However, it is difficult to relate a specific D J value to an actual wavenumber tracking result. To mitigate this issue, an example of wavenumber tracking with 40 sensors and SNR = 12 dB is illustrated in Fig. 6. It is clear that in this configuration, CPOMP performs really well (D J = 0.04). OMP (D J = 0.13) and COMP (D J = 0.13) also perform well, although they make a substantial number of false alarms. On the other hand, SoBaP (D J = 0.71) suffers from false alarms and missed detections, which explains its poor D J value.

E. Experimental results

The proposed method is now applied to the experimental CSS data described in Sec. V-B. Speed correction is first applied on the data using c min = 1459 m/s, a value empirically determined to timealign the data along the array. Wavenumber estimation is then performed using CPOMP and COMP.

The objective here is to exclusively compare the proposed method (CPOMP) to the best method of the state-of-the-art (COMP).

Wavenumbers are estimated using 10, 40 and 64 sensors. Experimental results are presented in Fig. 7.

Since no perfect ground truth is available to assess the experimental performance, the CPOMP 64-sensor May 27, 2021 DRAFT estimation is used as a reference (black dots in Fig. 7) to visually evaluate the results. One can see that CPOMP is consistently better than COMP. With 40 sensors, CPOMP gives results that are comparable to the 64-sensor reference. This is consistent with the simulated study which shows that CPOMP performance is roughly constant between 40 and 64 sensors (see Fig. 4). Interestingly, CPOMP also provides a good wavenumber estimation using only 10 sensors, although the last mode is missed. On the other hand, the COMP estimation with 10 sensors suffers from the small number of sensors, which leads to a consequent number of false alarms. The article presents a physics-based grid-free CS method to estimate modal wavenumber using a broadband source and a HLA with a small number of sensors. The method is based on three ideas. First of all, the method benefits from speed correction, so that a range-frequency signal can be conveniently sampled without wavenumber aliasing. Second, the method uses a grid-free framework to mitigate known CS drawbacks associated with basis mismatch. Last but not least, the method embeds physical information within the CS framework: it uses the dispersion relation (compensated for speed correction) to relate wavenumbers from one frequency to the other.

May 27, 2021 DRAFT

The proposed method has been benchmarked on simulations using the Jaccard distance as a performance metric. It was demonstrated that the method outperforms the state-of-the-art. Further, the method is experimentally validated on experimental data collected during SBCEX17. It notably allows a good estimation of modal wavenumbers from 0 to 100 Hz using as few as 10 sensors.

Although the proposed method works with a small number of sensors, it still requires a large horizontal aperture. As a result, a potential application for the method is modal estimation using synthetic aperture, as obtained using a fixed receiver and a moving source. This context is particularly appealing for geoacoustic inversion, where estimated modal wavenumbers can be used as an input to estimate the seafloor geoacoustic properties [START_REF] Ballard | Geoacoustic inversion for the New Jersey shelf: 3-D sediment model[END_REF], [START_REF] Frisk | Modal mapping experiment and geoacoustic inversion using sonobuoys[END_REF]. Following the CS paradigm, the proposed method can be used to collect fewer samples, but better ones. This has practical consequences for ocean acoustics, since the production of man-made source signal underwater is now considered as pollution [START_REF] Duarte | The soundscape of the anthropocene ocean[END_REF]. Reducing the number of samples (which is equivalent to the number of source signals in a synthetic aperture context) is an important perspective to reduce the noise footprint of ocean acoustic experiments.

  presents a simulated f-k diagram. The simulation has been run for a Pekeris waveguide representative of a classical shallow water scenario. For consistency, the environmental parameters are the same as those used in Refs. [21], [22]: water depth D = 130 m, sound speed c water = 1500 m/s and density ρ water = 1 g/cm 3 ; basement sound speed c bas = 2000 m/s and density ρ bas = 2 g/cm 3 . The f-k diagram was obtained using the FT of data simulated on a long and dense HLA with 240 regularly spaced sensors, with sensor spacing ∆r = 25 m (the array length is nearly 6 km).

Fig. 1 .

 1 Fig. 1. Simulated f-k diagram in a Pekeris waveguide.

Fig. 2 .

 2 Fig. 2. a) : Frequency-wavevumber diagram in the same configuration as in Fig. 1, except that frequency goes to 200 Hz: wavenumber are aliased for f > 60 Hz. b) : f-k diagram in the same configuration, but after speed correction: there is no aliasing.

Fig. 3 .

 3 Fig. 3. Performance of the estimation of modal wavenumbers in a simulated Pekeris waveguide using OMP: TP (left), FP (middle) and DJ (right).

  a 64 element HLA with a horizontal aperture of 1016 m deployed by the Norwegian Defence Research Establishment (FFI). The array orientation was roughly West-East, and the positions of the hydrophones at the HLA extremities are (40.4984N, 70.4677W) and (40.4983N, 70.4557W). The hydrophone spacing along the array is not uniform. The array elements have approximately logarithmic

SNR from 0

 0 Fig. 4. Jaccard's distance as a function of the number of sensors for a SNR of 12 dB.

Fig. 5 .

 5 Fig. 5. Jaccard's distance as a function of the SNR -using 32 sensors

Fig. 6 .

 6 Fig. 6. Wavenumber estimation on simulated data with 40 sensors and a SNR of 12 dB. Black: true wavenumbers ; blue: OMP, mallow: SoBaP, green: COMP, red: CPOMP.

Fig. 7 .

 7 Fig. 7. (black) Recovery of CPOMP with 64 sensors, (green) recovery with COMP, (red) recovery with CPOMP

[START_REF] Denoyelle | Theoretical and Numerical Analysis of Super-Resolution Without Grid[END_REF] DRAFT

ACKNOWLEDGMENT This work was supported by the Delegation Generale de L'Armement (DGA), the Office of Naval Research (ONR), the ONR Global, the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-16-IDEX-0005) operated by the French National Research Agency (ANR). May 27, 2021 DRAFT


Camille Delvigne
email: camille.delvigne@uclouvain.be

Abel Guihou

Jan A Schuessler

Paul Savage

Franck Poitrasson

Sebastian Fischer

Jade E Hatton

Katharine R Hendry

Germain Bayon

Emmanuel Ponzevera

Bastian Georg

Alisson Akerman

Oleg S Pokrovsky

Jean-Dominique Meunier

Pierre Deschamps

Isabelle Basile-Doelsch

Silicon isotope analyses of soil and plant reference materials: an inter-comparison of seven laboratories

Keywords: silicon isotopes, reference materials, soil, plant, inter-comparison of measurements Isotopes du silicium, matériaux

de référence, sol, plante, intercomparaison de mesures

Main message

This inter-laboratory comparison exercise establishes the δ 30 Si and δ 29 Si values for four new soil reference materials and one plant reference material, for future use as secondary Si isotope reference materials.

Introduction

Silicon (Si), the second most abundant element in the Earth's crust, is a ubiquitous rockforming element and is considered as a key component of global biogeochemical cycles, due to links with carbon cycling and other nutrient systems. Chemical weathering of silicate minerals controls atmospheric CO 2 concentration and thus regulates climate over geological timescales (e.g., [START_REF] Berner | The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years[END_REF]. Silicate weathering also releases dissolved Si into the hydrosphere, and thereby ultimately controls the biological pump, because Si is an essential nutrient in marine and freshwater ecosystems [START_REF] Nelson | Production and dissolution of biogenic silica in the ocean: Revised global estimates, comparison with regional data and relationship to biogenic sedimentation[END_REF][START_REF] Ragueneau | Si and C interactions in the world ocean: Importance of ecological processes and implications for the role of diatoms in the biological pump[END_REF].

Although the oceanic Si cycle has been a key research focus for decades, the terrestrial Si cycle has since attracted more scientific attention after the important and complex role of vegetation was recognized [START_REF] Conley | Terrestrial ecosystems and the global biogeochemical silica cycle[END_REF][START_REF] Street-Perrott | Biogenic silica: a neglected component of the coupled global continental biogeochemical cycles of carbon and silicon[END_REF][START_REF] Struyf | Emerging understanding of the ecosystem silica filter[END_REF]. Processes occurring at the land-ocean interface and anthropogenic impacts are now actively studied using Si isotopic variations or "fingerprints" [START_REF] Engström | Temporal isotopic variations of dissolved silicon in a pristine boreal river[END_REF][START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF][START_REF] Hughes | The riverine silicon isotope composition of the Amazon Basin[END_REF][START_REF] Mangalaa | Silicon cycle in Indian estuaries and its control by biogeochemical and anthropogenic processes[END_REF][START_REF] Vandevenne | Landscape cultivation alters δ³⁰Si signature in terrestrial ecosystems[END_REF].

Silicon has three stable isotopes: 28 Si, 29 Si, 30 Si with average relative isotope abundances of 92.23%, 4.67% and 3.10%, respectively [START_REF] Barnes | Absolute Isotopic Abundance Ratios and the Atomic Weight of a Reference Sample of Silicon[END_REF][START_REF] Bièvre | Table of the isotopic compositions of the elements[END_REF].

Pioneering Si isotope investigations in natural systems showed the potential to provide better understanding of the biogeochemical cycle of Si in both modern and past environments (see [START_REF] Basile-Doelsch | Si stable isotopes in the Earth's surface: A review[END_REF], for an early review). However, studies remained limited because of sometimes hazardous, but often delicate and time-consuming analytical methods. The advent of MC-ICP-MS instruments, combined with improved sample preparation techniques (e.g., [START_REF] Georg | New sample preparation techniques for the determination of Si isotopic compositions using MC-ICPMS[END_REF], have subsequently allowed Si isotope measurements to flourish and led to great advances in our understanding of the global Si cycle (see [START_REF] Poitrasson | Silicon Isotope Geochemistry[END_REF], for a more recent review).

For analyses, external matrix-matched reference materials are required to allow interlaboratory comparisons and to check for data quality for everyday analytical workparticularly in systems where standard-sample bracketing is the traditional way to measure isotope ratios. To date, only three inter-laboratory comparisons have been published for Si isotopes [START_REF] Reynolds | An inter-laboratory comparison of Si isotope reference materials[END_REF][START_REF] Hendry | Silicon isotopes in Antarctic sponges: An interlaboratory comparison[END_REF][START_REF] Grasse | GEOTRACES inter-calibration of the stable silicon isotope composition of dissolved silicic acid in seawater[END_REF]. The first one focused on three pure SiO 2 materials: Diatomite (δ 30 Si= +1.26 ± 0.20 ‰; from hereon, all reproducibilities will be given as 2 standard deviations unless otherwise stated), IRMM-018 (δ 30 Si= -1.65 ± 0.22 ‰) and Big Batch (δ 30 Si= -10.48 ± 0.54 ‰; [START_REF] Reynolds | An inter-laboratory comparison of Si isotope reference materials[END_REF]. The Diatomite reference material is still commonly included in analytical sequences to check for accuracy and precision; however, the last two are less frequently utilised because they are either highly fractionated (Big Batch) or potentially inhomogeneous (IRMM-018; [START_REF] Reynolds | An inter-laboratory comparison of Si isotope reference materials[END_REF]. For a better matching between δ 30 Si values of external standards and samples, and to account for any microstructural effect, a second inter-laboratory comparison reported on a Southern Ocean sponge material as an additional reference material for sponge spicule analyses [START_REF] Hendry | Silicon isotopes in Antarctic sponges: An interlaboratory comparison[END_REF]. A third inter-calibration study was motivated by the need to reliably compare the increasing number of Si isotope datasets of natural waters. External standards such as Diatomite (which is provided as a powdered, near-pure SiO 2 sample) cannot account for any biases introduced by the specific sample preparation and matrices involved when dealing with natural waters. These samples may contain Si with concentration below 20 µmol l -1 , which may require pre-concentrations steps or the processing of large amounts of sample, with potentially high salinity. With this in mind, [START_REF] Grasse | GEOTRACES inter-calibration of the stable silicon isotope composition of dissolved silicic acid in seawater[END_REF] characterised two seawater samples (ALOHA 300, [Si]= 9 µmol l -1 , δ 30 Si= +1.66 ± 0.13 ‰; ALOHA 1000, [Si]= 113 µmol l -1 , δ 30 Si= +1.25 ± 0.06 ‰;). The analyses of these "ALOHA" standards is now recommended when laboratories are reporting on seawater Si isotope data. In addition to this collection of reference materials, there are more than a dozen silicate rock reference materials that have enough extant published data to be suitable as external standards for silicate rock analyses (e.g., [START_REF] Abraham | δ 30 Si and δ 29 Si determinations on USGS BHVO-1 and BHVO-2 reference materials with a new configuration on a Nu plasma multi-collector ICP-MS[END_REF][START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF][START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF][START_REF] Savage | High temperature silicon isotope geochemistry[END_REF]. Among these reference materials, the Hawaiian basalt BHVO-2 (provided by the USGS) is the most commonly analysed. A valuable compilation of δ 30 Si values suggest that the reference value for BHVO-2 is δ 30 Si= -0.28 ± 0.03 ‰ [START_REF] Savage | High temperature silicon isotope geochemistry[END_REF] and more data are added continuously (see the GeoReM database).

Among the possible matrix effects on Si isotope measurements, the presence of anions such as sulfate [START_REF] Van Den Boorn | Sulfur-induced offsets in MC-ICP-MS silicon-isotope measurements[END_REF] or dissolved organic carbon (DOC) [START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF] has been reported to produce mass bias changes during MC-ICP-MS analysis. The issue remain debated as matrix effects induced by sulfate on Si isotopes were not subsequently confirmed [START_REF] Georg | New sample preparation techniques for the determination of Si isotopic compositions using MC-ICPMS[END_REF][START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF][START_REF] De Souza | Southern Ocean control of silicon stable isotope distribution in the deep Atlantic Ocean[END_REF] and could potentially depend on plasma conditions of the MC-ICP-MS [START_REF] Fietzke | Experimental evaluation of elemental behavior during LA-ICPMS: influences of plasma conditions and limits of plasma robustness[END_REF][START_REF] Yu | The impact of MC-ICP-MS plasma conditions on the accuracy and precision of stable isotope measurements evaluated for barium isotopes[END_REF]. A recent study provides constraints on matrix anion threshold concentrations and the extent of resulting biases, but the role of instrument type and different sample introduction systems (wet and dry plasma) has not been fully addressed [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF]. However, as far as we are aware, there are no Si isotope data for reference materials with an organic-rich matrix. Because such a matrix is prone to induce δ 30 Si bias during MC-ICP-MS analysis, a lack of suitable reference materials is an issue of major importance for the analysis of soil or plant material.

The aim of this paper is to provide the community with a set of easily available, wellcharacterized reference materials containing organic material as part of the matrix for use as secondary reference materials in Si isotope analysis. We propose four soil reference materials (i.e., GBW-07401, GBW-07404, GBW-07407, TILL-1), representative of solid samples with a complex matrix with various Si-bearing phases (i.e., amorphous and/or crystalline silicate minerals, refractory minerals like iron oxy-hydroxides) and organic matter.

In addition, the plant reference material ERM-CD281 (rye grass) was selected to represent the chemical composition of vegetation. It contains both organic matter and amorphous silica, which represents the form of Si generally accumulated in grass and then soils upon plant decay. This approach is an essential step in supporting δ 30 Si data comparison among and within laboratories over time.

Materials and methods

Investigated materials

This study investigates four soil reference materials GBW-07401, GBW-07404, GBW-07407, TILL-1 and one plant reference material ERM-CD281. The GBW soils "set", developed by the Institute of Geophysical and Geochemical Exploration (IGGE) of the Chinese Academy of Geological Sciences (CAGS) is not commercially available anymore. Because of this, an aliquot (~g) of the batches analysed in this study can be obtained upon request to Isabelle Basile-Doelsch at CEREGE. Reference material TILL-1 is provided by the CCRMP (Canadian Certified Reference Materials Project) and ERM-CD281 is provided by the EC-JRC (European Commission -Joint Research Center), and both are readily available to the community.

The reference material GBW-07401 is described as a dark brown podzolitic soil collected in a lead-zinc ore area, Heilongjiang, China; GBW-07404, as a limy-red soil obtained from Guangxi, China, and; GBW-07407, as a lateritic soil developed on basalt collected from Leizhou Peninsula, Guangdong, China (Wang et al., 2013). The reference material TILL-1 combines horizons B and C of a soil collected 25 km North-West of Lanark, Ontario, Canada [START_REF] Lynch | Provisional elemental values for four new geochemical soil and till reference materials, TILL-1, TILL-2, TILL-3 and TILL-4[END_REF]. Among those soil reference materials, SiO 2 content ranges from 32.67 to 62.60 % m/m, organic carbon (C org )content varies from 0.62 to 1.80 % m/m and sulfur content varies from 5 to 310 µg g -1 (more details in Table 1). The plant reference material ERM-CD281 is a rye grass, harvested in the United Kingdom in 1983, containing 0.13 % m/m Si (C org unknown) and 34 µg g -1 S (Table 1).

Sample preparation and silicon isotope measurements

The reference materials were analysed for their δ 29 Si and δ 30 Si values in seven laboratories: University of Bristol Isotope Group, United Kingdom (BIG); Centre Européen de Recherche et d'Enseignement en Géosciences de l'Environnement, France (CEREGE); Géosciences Environnement Toulouse, France (GET); German Research Centre for Geosciences, Germany (GFZ); Institut Français de Recherche pour l'Exploitation de la Mer, France (IFREMER); St Andrews Isotope Geochemistry, University of St Andrews, United Kingdom (STAiG); and Water Quality Centre -Trent University, Canada (WQC). The reference materials, provided as homogeneous powders by suppliers, were split and spread among the different laboratories via mail. Laboratories each received an aliquot of about 1 g of GBW soil reference materials, 0.2 g of TILL-1 and 3 g of ERM-CD281. Each laboratory conducted sample preparation and Si isotope measurements using their own protocols as detailed in Table 2. Sample preparation protocols of all laboratories are adapted from [START_REF] Georg | New sample preparation techniques for the determination of Si isotopic compositions using MC-ICPMS[END_REF] and hence are similar and consist of an alkaline fusion step followed by purification on a cation exchange resin. In the following, we highlight the differences between the sample preparation and measurement protocols. More details on the procedures used in each participating laboratory are given in previous publications or in Supplementary Materials. For BIG see [START_REF] Hatton | Investigation of subglacial weathering under the Greenland Ice Sheet using silicon isotopes[END_REF]; for CEREGE see Supplementary Materials; for GET see [START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF]; for GFZ see [START_REF] Schuessler | Testing the limits of micro-scale analyses of Si stable isotopes by femtosecond laser ablation multicollector inductively coupled plasma mass spectrometry with application to rock weathering[END_REF] and [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF]; for IFREMER see [START_REF] Bayon | The silicon isotopic composition of fine-grained river sediments and its relation to climate and lithology[END_REF]; for STAiG see [START_REF] Savage | Silicon isotopic variation in enstatite meteorites: Clues to their origin and Earth forming material[END_REF]; for WQC see [START_REF] Georg | New sample preparation techniques for the determination of Si isotopic compositions using MC-ICPMS[END_REF].

One major difference between protocols relates to the use of calcination to mineralize organic matter. CEREGE, GET and WQC calcinated the samples before fusion while BIG, GFZ and IFREMER did not; STAiG processed replicates both with and without calcination (see Supplementary Materials for detailed calcination procedures). All laboratories performed NaOH fusion in silver crucibles except for IFREMER who used NaOH and Na 2 O 2 in glassy carbon crucibles [START_REF] Bayon | The silicon isotopic composition of fine-grained river sediments and its relation to climate and lithology[END_REF]. Finally, in order to purify the samples, all laboratories used column chromatography with cation exchange resin (either AG50W-X12 or AG50W-X8, see Table 2) based on the procedure developed by [START_REF] Georg | New sample preparation techniques for the determination of Si isotopic compositions using MC-ICPMS[END_REF]. After alkali fusion, samples are normally taken up in a weakly acidic solution; many laboratories perform a measurement of solution Si concentration here (prior to Si purification), firstly to check for quantitative yield from the fusion step, and also to calculate the required amount of solution for column chemistry.

After column chemistry, the purified Si solutions were analysed by ICP-OES, (MC-)ICP-MS or by spectrophotometer using the molybdenum blue colorimetric method to i) verify quantitative Si column recovery (i.e., complete recovery within analytical reproducibilities of the concentration measurements of 100 ± 5%) and ii) check the purity of the Si sample solutions, with a particular focus on confirming the lack of magnesium (Mg) within the sample in the laboratories using Mg doping. Purity was better than 95% for Si, except for the plant ERM CD-281, where elevated levels of phosphorous and sulfur (mass ratios PO 4 /Si up to 5 and SO 4 /Si up to 6), respectively were still present. Complete removal of these anions was achieved neither by high temperature volatilisation during NaOH fusion at 750°C for 10 minutes (without prior calcination; GFZ) or calcination (CEREGE, no SO 4 data), nor by the Si column chemistry. However, as discussed below, all impurities were still within the tolerance limits found to not cause any bias in Si isotope measurements using the MC-ICP-MS setup (instrument, sample introduction system, analyte matrix, Mg doping,…) used by the participating laboratories in this study.

Silicon isotope measurements were performed by MC-ICP-MS, either Neptune Plus MC-ICP-MS or Neptune MC-ICP-MS (both from Thermo Fisher Scientific, Germany) (Table 2).

The analyses were performed in medium resolution (m/m ~ 4500, 5 and 95% peak side definition, see [START_REF] Weyer | High precision Fe isotope measurements with mass resolution MC-ICP-MS[END_REF]; except for GET which utilised high resolution (m/m ~7-10000)) using between 0.5 -3.6 µg Si per analysis, which resulted in a 3 -20 V ion beam on 28 Si. Four laboratories carried out analyses in 'wet plasma mode', three in 'semi-dry plasma mode' using a desolvating nebulizer system (ESI Apex-HF, ESI Apex Q). The main advantage of a desolvation unit is to reduce polyatomic interferences from oxygen, nitrogen, carbon and hydrogen and to improve sensitivity. However, use of a desolvation unit increases matrix effects and results in a less stable signal. Subtraction of the blank signal (measured on the pure acid matrix solution) from all sample Si intensities was performed before subsequent data processing by all laboratories (BIG, GFZ, CEREGE, GET, STAiG). For all laboratories, full procedural blank Si intensities did not exceed 1% of the sample Si signal (Table 2).

All measurements were performed using the standard-sample-bracketing technique using NBS 28 (also named NIST SRM 8546) as the bracketing standard, except for IFREMER where an in-house standard previously calibrated against the NBS 28 was used. With the exception of STAiG and WQC laboratories, the purified Si solutions were also doped with Mg (Si:Mg ratio of 1 to 3, depending on the laboratories; Table 2) prior to mass spectrometric analysis for an additional online mass bias drift correction using an exponential mass bias equation [START_REF] Cardinal | Isotopic composition of silicon measured by multicollector plasma source mass spectrometry in dry plasma mode[END_REF][START_REF] Engström | Chromatographic Purification for the Determination of Dissolved Silicon Isotopic Compositions in Natural Waters by High-Resolution Multicollector Inductively Coupled Plasma Mass Spectrometry[END_REF][START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF][START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF], (Table 2). When Mg-doping was utilised, Mg isotopes were analysed along Si isotopes using dynamic mode (magnet jump alternating between Si and Mg isotopes) for each sample doped with Mg. Magnesium doping is also useful to eliminate variations due to sample matrix effects or insufficiently stable laboratory temperature [START_REF] Cardinal | Isotopic composition of silicon measured by multicollector plasma source mass spectrometry in dry plasma mode[END_REF][START_REF] Engström | Chromatographic Purification for the Determination of Dissolved Silicon Isotopic Compositions in Natural Waters by High-Resolution Multicollector Inductively Coupled Plasma Mass Spectrometry[END_REF][START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF][START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF]. Additionally, BIG samples were doped with 50 µl 0.01 mol l -1 sulfuric acid per ml of sample to reduce anionic matrix effects. Silicon isotope compositions are reported relative to the NBS 28 in the delta notation as δ 30 Si and δ 29 Si:

δ 30, 29 Si (‰) = (R sample /R NBS 28 )-1
where R is the ratio 30 Si/ 28 Si or 29 Si/ 28 Si in a sample, and R NBS 28 the mean value for reference material NBS 28 calculated using the mean of the measured Si isotope ratio immediately measured before and after the unknown sample (i.e., standard-sample bracketing).

The well-characterized Diatomite and/or BHVO-2 reference materials, for which large datasets are available in the literature, were measured routinely for quality control as unknown samples during the runs by all laboratories to check for accuracy and precision. These data are provided in Table 2, and are well in agreement with the published values compiled in Table S1 (i.e., mean δ 30 Si = -0.28 ± 0.08 ‰ (U, k = 2, N = 61 where N is the number of data published)).

At least three aliquots of each reference material used in this study were processed through the whole sample preparation, hereafter referred to as full procedure replicates, except for GET (N=1), WQC (N=1) and the additional uncalcinated aliquot of GBW-07401 and GBW-07407 processed by STAiG. Each full procedure replicate was then analysed between 3 and 14 times either on different days or on the same day but not consecutively depending on laboratories (see Table 2 for details). Each mass-spectrometric analysis comprises between 20 to 36 signal integrations of 4 to 8.4 s, which is a total 80 to 210 s of Si signal acquisition per analysis. For each laboratory, replicate measurements carried out for each reference material full procedure replicate were averaged using the arithmetic mean and its associated 2 standard deviation (2SD; data are reported in Table S2).

Results and discussion

Data quality of participating laboratories

An essential measurement control is given by the three-isotope plot (Fig. 1). All data (from the reference materials full procedural replicates reported by each laboratory) fall on a straight line with a slope of about δ 29 Si ~ 0.5103 x δ 30 Si (r 2 > 99 %). The linear correlation encompasses the theoretical slopes calculated for both equilibrium (0.5178) and kinetic (0.5092) mass-dependent fractionation processes [START_REF] Engström | Silicon isotopic composition of boreal forest vegetation in Northern Sweden[END_REF], suggesting no major influence of isobaric interferences or matrix effects. In addition, the δ 30 Si and δ 29 Si values of the Diatomite and BHVO-2 reference materials measured by all laboratories agree well with compiled values [START_REF] Reynolds | An inter-laboratory comparison of Si isotope reference materials[END_REF][START_REF] Savage | High temperature silicon isotope geochemistry[END_REF]; our compilation in Table S1). The absolute differences between published values and those produced by participating laboratories are smaller than 0.06 ‰ and 0.03 ‰ for δ 30 Si and δ 29 Si, respectively, which are within the 2SD reported by each laboratory (Table 2).

WQC reported δ 30 Si calculated from measured δ 29 Si assuming an equilibrium massdependent fractionation. Comparing WQC δ 30 Si results with values reported by other laboratories (Fig 2.) shows no notable difference between their calculated δ 30 Si and other laboratories measured δ 30 Si, within their respective 2SD. This approach is therefore valid for the set of reference materials analysed in this study. However, this approach should be used with caution since it does not allow a data quality check for mass-dependent fractionation using the three-isotope plot (Fig 1 .).

Individual results of soil reference materials

The δ 30 Si and δ 29 Si values obtained for each reference material (representing full procedural replicates) by participating laboratories are illustrated in Fig. 2 and provided in Table S2. The Si isotope composition of the five reference materials reported are indistinguishable within their 2SD, both between the full procedural replicates for a given laboratory, and also between laboratories. The δ 30 Si data range between -0.30 ‰ < δ 30 Si < -0.20 ‰ for GBW-07401; -0.90 ‰ < δ 30 Si < -0.67 ‰ for GBW-07404; -2.05 ‰ < δ 30 Si < -1.71 ‰ for GBW-07407; -0.21 ‰ < δ 30 Si < -0.12 ‰ for TILL-1. This range of δ 30 Si values spans the typical δ 30 Si spectrum of soil materials ranging from -2.7 ‰ to +0.1 ‰ [START_REF] Poitrasson | Silicon Isotope Geochemistry[END_REF]. The most negative Si isotope compositions are measured in the reference materials with the lowest SiO 2 /Al 2 O 3 mass ratio (GBW-07407), i.e., the soil which has undergone the most desilification (Fig. 3). This is consistent with the general weathering trend, wherein the Si isotopic composition of secondary clay minerals formed during chemical weathering is lower than that of the primary silicate material (e.g., [START_REF] Opfergelt | Variations of δ 30 Si and Ge/Si with weathering and biogenic input in tropical basaltic ash soils under monoculture[END_REF][START_REF] Georg | Silicon fluxes and isotope composition of direct groundwater discharge into the Bay of Bengal and the effect on the global ocean silicon isotope budget[END_REF]Savage et al., 2013). Furthermore, the formation of more desilicified 1:1 clays is accompanied by a larger (negative) Si isotopic fractionation than the formation 2:1 clays [START_REF] Georg | Silicon fluxes and isotope composition of direct groundwater discharge into the Bay of Bengal and the effect on the global ocean silicon isotope budget[END_REF]Opfegelt et al., 2012).

Individual results of the plant reference material ERM-CD281

Of the laboratories involved in this study, two did not analyse the ERM-CD281 rye grass reference material because they had not previously processed highly organic-rich material (i.e., plant matter) and methodological development was beyond the scope of this study.

The δ 30 Si and δ 29 Si obtained for the plant reference material full procedure replicates by each laboratory are illustrated in Fig. 2 and provided in Table S2. The Si isotope compositions reported are indistinguishable within the 2SD both between the full procedural replicates for a given laboratory and between laboratories. The δ 30 Si data range between -0.36 ‰ < δ 30 Si < -0.26 ‰. In addition to the data acquired during this study, another published result for ERM-CD281 is available (δ 30 Si = -0.28 ± 0.08 ‰, 2SD, n = 3; [START_REF] Frick | Silicon isotope fractionation and uptake dynamics of three crop plants: laboratory studies with transient silicon concentrations[END_REF]). These published analyses were performed at GFZ using the same protocol in the same laboratory as reported for GFZ in this study. The Si isotope values for ERM-CD281 fall within the large range of δ 30 Si values reported in the literature for biogenic silica in plants (-2.3 to +6.1 ‰;[START_REF] Opfergelt | Silicon isotopes and continental weathering processes: Assessing controls on Si transfer to the ocean[END_REF]. The range of Si isotope compositions is large because the δ 30 Si values of biogenic silica in plants depend on several factors such as soil parent material lithology, soil weathering degree and on complex isotope fractionation processes during Si uptake and deposition [START_REF] Opfergelt | Variations of δ 30 Si and Ge/Si with weathering and biogenic input in tropical basaltic ash soils under monoculture[END_REF].

Influence of sulfur and DOC contents

The importance of chemical purification in removing so-called matrix effects has already been described for several isotopes systems measured by MC-ICP-MS (see for example Galy et al., 2001, for Mg isotopes). For Si isotopes, significant effects induced by the presence of sulfur in rock samples were first reported by [START_REF] Van Den Boorn | Sulfur-induced offsets in MC-ICP-MS silicon-isotope measurements[END_REF], i.e., a Si isotope bias of up to 1.4 ‰ for a SO 4 /Si mass ratio higher than 0.02. Subsequent studies have discussed the possibility of sulfur-induced bias on Si isotope measurements [START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF][START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF][START_REF] Chemtob | Silicon isotope systematics of acidic weathering of fresh basalts, Kilauea Volcano, Hawai'i[END_REF], and highlighted the potential for several parameters to influence the magnitude of the bias, notably in the dry plasma mode [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF]. Thus, care must be taken when comparing studies and drawing conclusions is not that straightforward. In particular, comparing sulfur-induced bias between S-doped solutions and S-rich rocks is likely irrelevant. However, in this study, sulfur-induced matrix effects are less likely to be significant, because the 720°C fusion step of solid samples used by all laboratories will at least partially remove the sulfur, with full removal achieved by heating at 1350°C ( [START_REF] Van Den Boorn | Sulfur-induced offsets in MC-ICP-MS silicon-isotope measurements[END_REF]. Despite different analytical settings, this can partly explain the large decrease from +1.4 to +0.3 ‰ (SO 4 /Si mass ratio ~0.3) between [START_REF] Van Den Boorn | Sulfur-induced offsets in MC-ICP-MS silicon-isotope measurements[END_REF] and [START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF] datasets as well as the absence of Si isotope bias measured on rock samples with SO 4 /Si mass ratios ≤0.14 by [START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF]. To allow comparison between studies, it would be useful to provide SO 4 /Si ratios of the Si analyte solutions in addition to the usual initial sample SO 4 /Si ratios before processing.

Restricting comparison to S-doped pure Si solutions (usually doped NBS 28), it appears that several parameters may reduce the isotopic bias. [START_REF] Van Den Boorn | Sulfur-induced offsets in MC-ICP-MS silicon-isotope measurements[END_REF] noticed that HCl 0.1 mol l -1 matrix reduces the bias by a factor about 1.5 compared to a HNO 3 matrix. Also, [START_REF] Chemtob | Silicon isotope systematics of acidic weathering of fresh basalts, Kilauea Volcano, Hawai'i[END_REF] noticed no significant bias due to sulfur for SO 4 /Si mass ratios <1 with a 0.002% HF and 2% HNO 3 matrix despite very similar analytical parameters (i.e., the use of a Neptune, a CETAC desolvator and no Mg doping). It appears that the overall sample aliquot matrix plays a role in the magnitude of the S-induced bias in addition to subtle differences in instrumental set-up. Similarly, [START_REF] Georg | New sample preparation techniques for the determination of Si isotopic compositions using MC-ICPMS[END_REF] noticed no deviation from expected values for SO 4 /Si mass ratio up to 48 using a Nu instrument, a CETAC dessolvator in a 0.1 mol l -1 HCl matrix without Mg doping. Similarly, [START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF] and [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF] investigated the effect with Mg doping. For SO 4 /Si ratios < 1, the Mg doping had no effect on the magnitude of the S-induced bias [START_REF] Zambardi | Precise Determination of Silicon Isotopes in Silicate Rock Reference Materials by MC-ICP-MS[END_REF][START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF]. Moreover, [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF] found that the addition of Mg to the Si analyte solution resulted in higher tolerance to impurities during MC-ICP-MS analyses when an Apex desolvator without a membrane was used. Remaining anionic impurities (PO 4 3-, SO 4 2-, and NO 3 -) had no effect on the Si isotope mass bias (with accuracy within ± 0.14 ‰) if the mass ratio of those anions to Si is lower than 6 [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF]. To overcome Sinduced bias, different techniques were suggested: removal of sulfur by heating the samples up to 1350°C (Van den [START_REF] Van Den Boorn | Sulfur-induced offsets in MC-ICP-MS silicon-isotope measurements[END_REF] or doping both samples and bracketing standard with significant amounts of sulfuric acid to overwhelm any variation in the sample [START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF].

In this study, no S-induced bias on Si isotopes measurements is expected for soil reference materials as SO 4 /Si mass ratios are below 0.01 (SO 4 /Si mass ratio ~0.003 for GBW-07401; 0.002 for GBW-07404; 0.005 for GBW-07407 and 0.002 for TILL-1). In contrast, the plant reference material ERM-CD281 has a SO 4 /Si mass ratio about 7.8 that could potentially be problematic for Si isotopes measurements. This was verified in Si analyte solutions where elevated levels of phosphorous and sulfur (mass ratios PO 4 /Si up to 5 and SO 4 /Si up to 6), respectively, were still present after NaOH fusion (with or without prior calcination) and Si column purification. Complete removal of these anions was achieved neither by high temperature volatilisation nor by the Si column chemistry. However, all impurities were still within the limits required to avoid any bias for all different analytical configurations used by different laboratories (prior calcination or not, wet plasma or Apex desolvation, Mg doping or no Mg doping, HCl or HNO 3 matrix). This is consistent with the study of [START_REF] Oelze | Mass bias stabilization by Mg doping for Si stable isotope analysis by MC-ICP-MS[END_REF] where no effect on the Si isotope mass bias is observed if the mass ratio of anions to Si is lower than 6 with an Mg correction (1 without an Mg correction). The Si isotope results of ERM CD-281 measured by laboratories using a regular spray chamber (wet plasma) agree well with those measured using an Apex sample introduction system. Similarly, the prior calcination of the sample has a limited efficiency (Fig. 2) and does not seem mandatory for SO 4 removal in such samples. In all cases, the chemical preparation and analytical settings seem to be robust enough to avoid matrix effects caused by residual impurities in the sample during Si isotope measurements. [START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF] also showed that measurable concentrations of dissolved organic carbon (DOC) in freshwater samples can affect their Si isotopic determination. Given the large spectrum of naturally occurring DOC forms, DOC-induced bias is suspected to be variable from one DOC form to another even before considering combined effects. This would make correcting for the presence of DOC in samples unworkable, due to the nearimpossibility of matching the matrix in the bracketing standards. Therefore, removal of organic matter during sample preparation seems to be the best option to ensure reliable analysis of organic carbon-rich samples.

The variability of the SiO 2 and organic carbon contents of the soils reference materials used in this study allow for evaluation of the influence of varying matrix to analyte ratios with or without calcination steps. When comparing the results with or without organic matter mineralization, no notable offset is observed, regardless of the considered reference material (Fig. 2). In addition, STAiG laboratory calcinated their reference materials but also processed one aliquot each of GBW-07401 and GBW-07407 without calcination (Fig. 2). The results for both calcinated and uncalcinated aliquots of these two reference materials are identical within their 2SD (Fig. 2 and Table S2). In addition, at GFZ the removal of organic carbon by the NaOH fusion method (750°C, 10 minutes) -without prior calcination -was investigated. Effective organic carbon removal was found for all soil and plant samples processed at GFZ during this study using semi-quantitative ICP-OES analyses (carbon 193.027 nm emission line) by evaluation of C remaining in the post-fusion and purified Si solutions prior to Si isotope analyses. All solutions showed C-counts indistinguishable from pure acid used for dilution for the analyses, whereas significant C intensities were recorded in sample aliquots that were not treated by alkali fusion and subsequent Si column purification. It thus seems that the presence of organic carbon has no impact on the accuracy or reproducibility of Si isotopes measurements if the high temperature alkaline fusion digestion is used (with or without prior calcination) as the remaining amount of organic matter is insufficient to induce Si isotopic bias. Several factors could explain this. The amount of organic matter in uncalcinated aliquots was insufficient to induce any bias. Indeed, C org /Si mass ratios of soil reference materials were below 0.1 while [START_REF] Hughes | Controlling the mass bias introduced by anionic and organic matrices in silicon isotopic measurements by MC-ICP-MS[END_REF] observed a shift in freshwater samples from DOC/Si > 0.2. However, the Si isotope measurement of the high C org /Si plant reference material ERM-CD281 analysed at IFREMER and GFZ -without prior calcination step -are not affected. In this case, it could be that the organic carbon in this sample is not refractory and is oxidised during the high-temperature fusion step as suggested by C analyses at GFZ. Finally, it is also possible that the DOC form present in soils or plants does not induce bias, at least not on the Thermo Neptune-family of MC-ICP-MS instruments.

More tests would be required to determine the critical DOC/Si ratio and if it is DOC-form sensitive.

3.5. Assigning δ 30 Si and δ 29 Si to the reference materials analysed in this study

Results from all seven laboratories are in excellent agreement for all reference materials (Fig. 2), despite the use of different sample preparations (calcination versus no calcination;

NaOH fusion versus NaOH and Na 2 O 2 fusion) and different analytical settings (Mg doping vs no doping, wet plasma sample introduction vs semi-dry plasma; see Table 2). Therefore, the results reported by the laboratories for each reference material full procedural replicates (with the exception of the δ 30 Si reported by WQC) were weighted and combined together to calculate averages and their associated combined uncertainties (U, k = 2). Final isotopic characterisation of the five reference materials analysed here are given in Table 3. Following the EURACHEM/CITAC Guide CG4, the combined uncertainties were calculated by error propagation, assuming independent variables, of the standard deviations obtained on full procedure replicates for each reference material, also including the standard deviation of the averages of full procedure replicates for each reference material. The combined uncertainties were multiplied by a coverage factor of k = 2 to obtain the expanded uncertainties for each reference material. The uncertainties (U, k = 2) range from 0.04 ‰ to 0.12 ‰ for δ 29 Si and from 0.06 ‰ to 0.17 ‰ for the δ 30 Si. The larger magnitude of U for δ 30 Si compared to those for δ 29 Si correspond with the lower natural abundance of 30 Si compared with 29 Si and the concomitant deterioration in counting statistics of the detectors.

In addition, the uncertainties of these results (Table 3) are comparable with standard deviation obtained on multiple measurements of the reference materials Diatomite and BHVO-2 (Table 2) used by the laboratories to check for accuracy. Overall, the reference materials homogeneity is valid at the scale of the minimum sample mass used in this study, i.e., >3 mg for the soils reference materials, and at least 50 mg for rye grass ERM CD281. The soils RM GBW-07401, GBW-07404, GBW-07407 and TILL-1 offer a wide range of δ 29 Si and δ 30 Si values comparable to and therefore useful for future Si isotope environmental studies.

The uncertainties obtained in this inter-laboratory comparison are more than one order of magnitude lower than the range of δ 29 Si and δ 30 Si measured in natural samples: -5.7 to +6.1 ‰ [START_REF] Opfergelt | Silicon isotopes and continental weathering processes: Assessing controls on Si transfer to the ocean[END_REF] and are therefore sufficient to allow the robust interpretation of small variations in the δ 29 Si and δ 30 Si.

Conclusion

In this study, the δ 29 Si and δ 30 Si values of four soil and one plant reference materials were determined and compared between laboratories. Different sample preparation procedures and analytical settings provide consistent results for both δ 29 Si and δ 30 Si values for all reference materials. Therefore, all protocols used by the participating laboratories for Si isotope analysis provide comparable results. Also, no notable isotopic bias was observed between calcinated and uncalcinated reference materials with different C org /Si ratios.

However, caution still stands for analyses of organic carbon-rich samples with a C org /Si ratio above the studied reference materials, where biases in δ 30 Si measurements could potentially occur. This is particularly the case when sample preparation does not involve calcination, leading to potential matrix effects during MC-ICP-MS analyses that could be DOC-form dependent. The data for soil and plant reference materials from this study can be reliably used by other laboratories for routine quality control of Si isotope ratio measurements during environmental and geochemical studies. Moreover, the wide range of δ 29 Si and δ 30 Si covered by these soil reference materials provided here allows analysts to select the reference materials that are best suited for their study.

Acknowledgements

At CEREGE, this project was initiated by an EIL funding from the Société Française des Isotopes Stables (SFIS) and also supported by the Agence Nationale pour la Recherche (ANR, France) through the project EQUIPEX ASTER-CEREGE and the project BIOSiSOL (ANR-14-CE01-002). C.D. is funded by the "Fonds National de la Recherche Scientifique" (Belgium). SF and the Si isotope measurements at STAiG were supported by NERC grant NE/R002134/1 to PS; PS would also like to cite the support of a Carnegie Trust Research Incentive Grant, which helped the setup of various isotope techniques in the St Andrews Isotope Geochemistry (STAiG) laboratories. F. von Blanckenburg and the Helmholtz Association are thanked for infrastructure support at GFZ. At the University of Bristol, C. D. Coath is thanked for laboratory support and the European Research Council is acknowledged for funding (ICY-LAB, grant agreement 678371). Hélène Mariot is thanked for her careful maintenance of the CEREGE clean lab. Manuel Henry and Jérôme Chmeleff are thanked for maintaining the GET clean lab and MC-ICP-MS facilities in good working order. J. Schlegel and J. Buhk are acknowledged for laboratory support at GFZ. The authors declare that there is no conflict of interest.

Data availability Statement

The authors confirm that the data supporting the findings of this study are available within the article and its Supplementary Materials. Detailed datasets are available from the corresponding author upon reasonable request.