Lapo Boschi

Fabrizio Magrini

Fabio Cammarano

Mark Van Der Meijde

On seismic ambient noise cross-correlation and surface-wave attenuation

We derive a theoretical relationship between the cross correlation of ambient Rayleigh waves (seismic ambient noise) and the attenuation parameter α associated with Rayleigh-wave propagation. In particular, we derive a mathematical expression for the multiplicative factor relating normalized cross correlation to the Rayleigh-wave Green's function. Based on this expression, we formulate an inverse problem to determine α from cross correlations of recorded ambient signal. We conduct a preliminary application of our algorithm to a relatively small instrument array, conveniently deployed on an island. In our setup, the mentioned multiplicative factor has values of about 2.5 to 3, which, if neglected, could result in a significant underestimate of α. We find that our inferred values of α are reasonable, in comparison with independently obtained estimates found in the literature. Allowing α to vary with respect to frequency results in a reduction of misfit between observed and predicted cross correlations.

Introduction

A number of theoretical and experimental studies have proved that the cross correlation between two recordings of a diffuse surface-wave field approximately coincides with the surfacewave Green's function associated with the two points of observation. This is relevant to the field of seismology, since recorded seismic ambient signal has been found to mostly consist of seismic surface waves, and empirical Green's functions are now routinely retrieved from the cross correlation of seismic ambient noise [e.g. Boschi and Weemstra, 2015, and references therein]. Most authors have been able to derive the medium's velocity from the phase of the reconstructed Green's functions; this resulted in successful applications of ambient-noise theory to imaging and monitoring [see the reviews by, e.g., [START_REF] Campillo | Seismic imaging and monitoring with ambient noise correlations[END_REF][START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]. The amplitude of the Green's function in principle provides complementary information on the medium's anelastic properties; but it tends to be less accurately reconstructed by cross correlation. Initial attempts to constrain surface-wave attenuation from ambient noise [e.g. [START_REF] Prieto | Anelastic earth structure from the coherency of the ambient seismic field[END_REF][START_REF] Harmon | Distribution of noise sources for seismic interferometry[END_REF][START_REF] Lawrence | Attenuation tomography of the western United States from ambient seismic noise[END_REF][START_REF] Weemstra | Seismic attenuation from recordings of ambient noise[END_REF] were based on the assumption (questioned by [START_REF] Weaver | On the amplitudes of correlations and the inference of attenuations, specific intensities and site factors from ambient noise[END_REF]) that the "lossy" Green's function be simply the product of the elastic Green's function times an exponential damping term. [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF] validated mathematically this assumption, but both he and Harmon et al.

[2010] emphasized the difficulty of constraining earth's attenuation whenever the ambient field is not perfectly diffuse. The study of [START_REF] Weemstra | On estimating attenuation from the amplitude of the spectrally whitened ambient seismic field[END_REF] additionally showed that data-processing techniques typically used in ambient-noise literature, such as whitening or time-domain normalization [e.g. [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF], could also affect attenuation estimates.

Whitening and/or normalization, however, are necessary to avoid that localized events of relatively large amplitude, like earthquakes, obscure random noise and thus bias cross correlations.

It is the purpose of this study to introduce a new normalization criterion, which is in practice similar to whitening, but is derived directly from the reciprocity theorem as stated, e.g., by [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]; i.e., it does not follow from data-processing considerations, but from the physics of ambient-noise cross correlation. The relationship we obtain (eq. (30) in sec. 2.3.2), and on whose basis we formulate an inverse problem to determine attenuation, can be summarized simplistically by stating that cross correlations are normalized against the power spectral density of emitted noise. An approximate equation expressing source power spectral density in terms of recorded data is also found (eqs. ( 27) and (29), sec. 2.3.2). To achieve all this, we assume the same theoretical framework of, e.g., [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], describing Love-and Rayleigh-wave as combinations of two-dimensional membrane waves and "radial eigenfunctions." Our treatment (sec. 2) involves an independent derivation of the surface-wave Green's function in a lossy medium (sec. 2.1 and app. B).

Importantly, the fact that the medium is lossy requires that noise sources be uniformly distributed over the surface of the earth (and not just over all azimuths) for the Green's function to be reconstructed by noise cross-correlation [START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF][START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF]. In addition, our formulation is strictly valid only if the spectrum of signal emitted by ambient-noise sources is laterally homogeneous (sec. 2.3). Because seismic ambient noise is mostly originated by coupling between the solid earth and the oceans, it is reasonable to expect that both requirements will be best approximated by deploying instruments on an island. We accordingly test our algorithm on two years of continuous data recorded by an array of fourteen broadband stations in Sardinia. We explore two different parameterizations of attenuation in the area of interest, i.e. constant (sec. 3.2.1) vs. frequency-dependent (sec. 3.2.2) attenuation parameter. In both cases, we identify an attenuation model that minimizes data misfit. We discuss the resulting models in light of independent studies of Rayleigh-wave attenuation.

Theory 2.1 Green's problem for a lossy membrane

Let us first summarize the treatment of surface-wave theory given by [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], based on earlier studies by [START_REF] Tanimoto | Modelling curved surface wave paths: membrane surface wave synthetics[END_REF] and [START_REF] Tromp | Variational principles for surface wave propagation on a laterally heterogeneous Earth-III. Potential representation[END_REF]. We shall work in the frequency (ω) domain and employ the Fourier-transform convention of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]. It is convenient to first introduce the Rayleigh-(u R ) and Love-wave

(u L ) displacement Ansätze u R (x 1 , x 2 , x 3 , ω) = U (x 3 , ω)x 3 φ R (x 1 , x 2 , ω) + V (x 3 , ω)∇ 1 φ R (x 1 , x 2 , ω),
(1)

u L (x 1 , x 2 , x 3 , ω) = W (x 3 , ω)(-x 3 × ∇ 1 )φ L (x 1 , x 2 , ω), (2)
respectively, where x 1 , x 2 , x 3 are Cartesian coordinates, with the x 3 -axis perpendicular to Earth's surface (which we assume to be flat) and oriented downward; the unit-vectors x 1 , x 2 ,

x 3 are parallel to the Cartesian axes; ∇ 1 denotes the surface gradient

x 1 ∂ ∂x 1 + x 2 ∂ ∂x 2 .
The functions U (x 3 , ω), V (x 3 , ω) and W (x 3 , ω) control the dependence of surface-wave amplitude on depth, and the functions φ R (x, ω) and φ L (x, ω) are respectively dubbed Rayleigh-and Love-wave "potentials". Upon substituting expressions (1) and ( 2) into the frequency-domain displacement equation, it is found that φ R and φ L can be determined through the Helmholtz' equations

∇ 2 1 φ L,R (x 1 , x 2 , ω) + ω 2 c 2 L,R φ L,R (x 1 , x 2 , ω) = f (x 1 , x 2 , ω), (3)
where c L,R (ω) denote the value of Rayleigh-or Love-wave phase velocity at frequency ω and f is a generic forcing term.

This study is limited to recordings of seismic ambient noise. Because seismic noise has been shown to essentially amount to surface waves, it is safe (provided that signals related to large or nearby earthquakes are excluded) to assume that eqs. ( 1) and (2) correctly describe the corresponding ground displacement. Furthermore, for the sake of simplicity we only consider vertical-component recordings, i.e., the

x 3 -component of u R , or u R,3 = U (0, ω)φ R (x 1 , x 2 , ω), (4)
with x 3 =0 as long as recordings are made at Earth's surface. It is inferred upon multiplying eq. ( 3) by U (0, ω) that the following vertical-displacement equation holds,

∇ 2 1 u R,3 (x 1 , x 2 , ω) + ω 2 c 2 R u R,3 (x 1 , x 2 , ω) = U (0, ω)f (x 1 , x 2 , ω). (5)
Eqs.

(3) and ( 5) coincide with the equation governing the displacement of a lossless, stretched membrane [e.g. [START_REF] Kinsler | Fundamentals of Acoustics[END_REF]. Following [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], we define the Green's function G 2D as the membrane response to impulsive (Dirac δ(x 1 )δ(x 2 )) initial velocity at the reference-frame origin, with zero initial displacement and zero forcing term f ; it follows (app. A.1) that, for Rayleigh waves,

G 2D (x 1 , x 2 , ω) = - iP 4 √ 2πc 2 R H (2) 0 ωx c R , ( 6
)
where i denotes the imaginary unit, H

(2) 0 the zeroth-order Hankel function of the second kind [e.g. Abramowitz and Stegun, 1964, eq. (9.1.4)], P accounts for the physical dimensions of G 2D as discussed in app. A.1, and x = x 2 1 + x 2 2 is the distance between (x 1 , x 2 ) and the "source."

Following [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF], we now make the assumption that surface-wave attenuation can be accounted for by replacing eq. ( 5) with a damped membrane equation, i.e. introducing an additional forcing term, or "loss term proportional to, and oppositely directed from, the velocity of the vibrating element" [Kinsler et al., 1999, sec. 4.6]. It is convenient to denote 2α c R the proportionality factor between force and velocity, with the "attenuation coefficient"

α coinciding with that of [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF]. In the frequency domain, the resulting displacement equation reads

∇ 2 1 u(x 1 , x 2 , ω) + ω 2 c 2 -i 2αω c u(x 1 , x 2 , ω) = U (0, ω)f (x 1 , x 2 , ω) (7)
[e.g. [START_REF] Kinsler | Fundamentals of Acoustics[END_REF][START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF], where all unnecessary subscripts have been dropped.

It is inferred by comparing eqs. ( 5) and ( 7) that the expression ( 6) is still a solution of ( 7), if the real ratio ω/c in its argument is replaced by the complex number ω 2 c 2 -2iαω c [START_REF] Kinsler | Fundamentals of Acoustics[END_REF][START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF][START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF][START_REF] Weemstra | On the attenuation of the ambient seismic field: Inferences from distributions of isotropic point scatterers[END_REF]; the two-dimensional, damped

Green's function therefore reads

G d 2D (x 1 , x 2 , ω) = - iP 4 √ 2πc 2 H (2) 0 x ω 2 c 2 - 2iαω c . (8)
We show in app. B that as long as attenuation is relatively weak, i.e. α ω/c [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF], and provided that frequency is high and/or the effects of near-field sources are negligible, expression (8) can be reduced to the more convenient, approximate form

G d 2D (x 1 , x 2 , ω) ≈ - iP 4 √ 2πc 2 H (2) 0 ωx c e -αx ≈ G 2D (x 1 , x 2 , ω)e -αx ; (9)
in other words, if ωx c 1 and α ω/c, the lossy-membrane Green's function G d 2D can be roughly approximated by the product of the lossless two-dimensional Green's function G 2D times a damping term e -αx ; we verified this result via numerical tests (Fig. 1).

Reciprocity theorem for a lossy membrane

We extend the reciprocity-theorem derivation of [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], who only considered lossless media, to the case of a lossy membrane. The procedure that follows is similar to that of [START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF], which, however, is limited to lossy three-dimensional media. In analogy with [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], let us introduce a vector v = -1 iω ∇ 1 u, such that

∇ 1 u + iωv = 0. ( 10
)
Figure 1: Comparison between the exact (red) and approximate (blue) formulae (eqs. ( 8) and ( 9), respectively) for the lossy-membrane Green's function. The real parts of the Green's function are shown on the left (panels a, c, e, g), the imaginary parts on the right (b, d, f and h). Green's functions corresponding to interstation distances of 50 km (a through d) and 200 km (e through h), and α=2 × 10 -5 m -1 (a, b, e and f) and α=5 × 10 -5 m -1 (c, d, g and h) are shown. The approximation is good when α is small and ω is high, and decays with growing α and decreasing ω.

Upon substituting (10) into the damped-membrane displacement equation ( 7),

∇ 1 • (-iωv) + ω 2 c 2 -i 2αω c u = U (0, ω)f. (11)
After some algebra,

∇ 1 • v + i ωκ c 2 u -q = 0, (12)
where κ = 1i 2αc ω for brevity, and q= i ω U (0, ω)f to be consistent with [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF].

In the frequency domain, q has units of squared time over distance. [START_REF] Kinsler | Fundamentals of Acoustics[END_REF] show that, in the stretched-membrane model, this forcing can be thought of as proportional to "pressure divided by surface density".)

Eqs. ( 10) and ( 12) are similar to eqs. ( 24) and ( 25) of [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], except that the real number ω c 2 in eq. ( 25) of [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF] is replaced here by the complex ωκ c 2 .

Following [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], consider an area S on the membrane, bounded by the closed curve ∂S; let q A (x 1 , x 2 , ω), u A (x 1 , x 2 , ω) and v A (x 1 , x 2 , ω) denote a possible combination of the fields q, p and v co-existing at (x 1 , x 2 ) in S and ∂S. A different forcing q B would give rise, through eqs. ( 10) and ( 12), to a different "state" B, defined by u B (x 1 , x 2 , ω) and

v B (x 1 , x 2 , ω).
The reciprocity theorem is obtained by combining the left-hand sides of eqs.

(10) and ( 12) as follows,

S d 2 x [(10) A • v * B + (10) * B • v A + (12) A u * B + (12) * B u A ] = 0, (13)
where x=(x 1 , x 2 ), d 2 x = dx 1 dx 2 , and * denotes complex conjugation. (10) A is short for the expression one obtains after substituting u = u A (x, ω) and v = v A (x, ω) into the left-hand side of eq. ( 10), etc. Following the same procedure as in sec. 3 of [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], we find

S d 2 x [∇ 1 (u A • v * B ) + ∇ 1 (u * B • v A )] + iω c 2 (κ -κ * ) S d 2 x u A u * B = S d 2 x (q A u * B + q * B u A ) .
(14) Notice that the second term at the left-hand side of eq. ( 14) would be 0 if the membrane were lossless (α=0, and therefore κ=κ * ); it is easy to see that in that case ( 14) is equivalent eq. ( 31) of [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF].

The divergence theorem allows to reduce the first surface integral at the left-hand side of (14) to a line integral along ∂S. Following [START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF], we consider the particular case where surface integration is over the entire two-dimensional space R 2 , i.e. the area S is infinite. This is relevant to seismic ambient-noise applications, where receiver arrays are typically deployed within a relatively small area, receiving signal from "noise" sources that are distributed with (approximately) equal probability over the entire surface of the globe.

Then, for attenuating media the wave field vanishes exponentially at infinity, and the integral along ∂S accordingly vanishes. We are left with

4α c R 2 d 2 x u A u * B = R 2 d 2 x (q A u * B + q * B u A ) . ( 15
)
Consider now the states A and B associated with impulsive forcing terms q A =F δ(x -x A )

and q B =F δ(x -x B ), respectively. x A and x B are two arbitrary point-source locations, and the factor F accounts for the physical dimensions of q (recall that δ(x) has dimensions of one over squared distance). It follows from app. A, eq. (A.22), that the corresponding membrane

displacements are u A = iω F c 2 P G d 2D (x, x A , ω) and u B = iω F c 2 P G d 2D (x, x B , ω), respectively.
Substituting into eq. ( 15), F simplifies out and

2αωc P R 2 d 2 x G d 2D (x, x A , ω)G d * 2D (x, x B , ω) = -[G d 2D (x A , x B , ω)], (16)
which is the sought reciprocity theorem. Eq. ( 16) is similar, e.g., to eq. ( 39) of [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF], with one fundamental difference: the integral in ( 16) is not over a closed curve containing the receiver pair (as in [START_REF] Boschi | A simple method for earthquake location by surfacewave time-reversal[END_REF]), but over the entire real plane. As shown in the following, this implies that, for the lossy Green's function G d 2D to be accurately reconstructed by seismic interferometry, sources should be uniformly distributed over space, rather than azimuth, and both in the near and far field of the receivers [START_REF] Snieder | Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves[END_REF][START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF][START_REF] Weemstra | On the attenuation of the ambient seismic field: Inferences from distributions of isotropic point scatterers[END_REF].

Cross-correlation amplitude as a constraint for attenuation

We shall next (sec. 2.3.1) use eq. ( 16) to establish a relationship between the cross correlation of ambient surface-wave signal ("noise") recorded at two locations x A , x B , and the imaginary part of the Green's function

G d 2D (x A , x B , ω).
Based on this relationship, in sec. 2.3.2 an inverse problem will be formulated, having attenuation α as unknown parameter, and the cross correlation of recorded noise as data. In this endeavour, it is assumed that ambient noise can be represented by a distribution of point sources of constant spatial density, emitting at random times (i.e., random phase in the frequency domain) but, in analogy with the numerical study of [START_REF] Cupillard | On the amplitude of surface waves obtained by noise correlation and the capability to recover the attenuation: a numerical approach[END_REF], all sharing the same spectral amplitude.

The assumption that the spectral amplitude of noise sources be constant across the globe is based on the idea that Rayleigh-wave noise on earth is generated by the coupling, at the ocean bottom, between oceans and the solid earth. It has been shown that, while local effects play a role, the resulting spectrum has maxima determined by the main frequencies of ocean waves (i.e., primary and secondary microseisms at 0.05-0.12 and 0.1-0.25 Hz, respectively), independent of location [e.g. [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF][START_REF] Ardhuin | Modelling secondary microseismic noise by normal mode summation[END_REF][START_REF] Hillers | Global oceanic microseism sources as seen by seismic arrays and predicted by wave action models[END_REF].

We accordingly consider our assumption to be valid at least as a rough approximation of the real world.

Noise cross-correlation and the Green's function

The vertical-component, Rayleigh-wave displacement associated with a noise "event" can be thought of as the time-domain convolution of G d 2D and a source time function. In the frequency domain, convolution is replaced by product, and the signal emitted at a point x and recorded at, say, receiver x A reads h(ω)G d 2D (x A , x, ω)e iωφ , with h(ω) and φ denoting the amplitude and phase of the emitted signal. It follows that the ambient noise recorded at x A can be written

s(x A , ω) = h(ω) N S j=1 G d 2D (x A , x j , ω)e iωφ j , (17)
where N S denotes the total number of sources, and the index j identifies the source; as anticipated, it is assumed that the unitless, frequency-domain amplitude h(ω) is approximately the same for all noise sources. Based on eq. ( 17), the cross correlation of noise recorded at two receivers x A , x B can be written

s(x A , ω)s * (x B , ω) =|h(ω)| 2   N S j=1 G d 2D (x A , x j , ω)e iωφ j   N S k=1 G d * 2D (x B , x k , ω)e -iωφ k =|h(ω)| 2 N S j=1 G d 2D (x A , x j , ω)G d * 2D (x B , x j , ω) + N S j=1 N S k=1,k =j G d 2D (x A , x j , ω)G d * 2D (x B , x k , ω)e iω(φ j -φ k ) . (18)
The phases φ 1 , φ 2 , φ 3 , . . . are assumed to be random (uniformly distributed between 0 and 2π); it follows that the cross correlations of signals emitted by different sources, i.e. the second term at the right-hand side of (18) ("cross terms"), can be neglected if noise is recorded over a sufficiently long time, or if a sufficiently large amount of uniformly distributed sources are present [e.g. [START_REF] Weemstra | On estimating attenuation from the amplitude of the spectrally whitened ambient seismic field[END_REF]Boschi and Weemstra, 2015, App. D]; then

s(x A , ω)s * (x B , ω) ≈ |h(ω)| 2 N S j=1 G d 2D (x A , x j , ω)G d * 2D (x B , x j , ω). (19)
It is convenient to transform the sum at the right-hand side of eq. ( 19) into an integral; said ρ the surface density of noise sources, which we assume to be constant,

s(x A , ω)s * (x B , ω) ≈ ρ |h(ω)| 2 R 2 d 2 x G d 2D (x A , x, ω)G d * 2D (x B , x, ω). ( 20
)
It will be noticed that the integral in ( 20) is over the entire real plane: this follows from the assumption, made in sec. 2.1, that sources be uniformly distributed over R 2 . Dividing both sides by ρ|h(ω)| 2 , we find from eq. ( 20) that

R 2 d 2 x G d 2D (x A , x, ω)G d * 2D (x B , x, ω) ≈ 1 ρ|h(ω)| 2 s(x A , ω)s * (x B , ω), (21)
and finally, substituting eq. ( 21) into (16),

s(x A , ω)s * (x B , ω) ≈ - |h(ω)| 2 P ρ 2αωc [G d 2D (x A , x B , ω)]. (22)
Eq. ( 22) is the "lossy" counterpart of, e.g., eq. ( 65) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF].

Let us emphasize, again, that eq. ( 22) was obtained under the assumption that sources are uniformly distributed over the entire real plane R 2 . It follows from eq. ( 22) that, in the absence of attenuation, i.e. α=0, the cross correlation of ambient signal at its left-hand side is divergent. This is why in ambient-noise literature, whenever attenuation is neglected, the assumption is made that sources are uniformly distributed with respect to azimuth, rather than in space: for instance, the mentioned eq. ( 65) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] results from a uniform distribution of sources along a circle that surrounds the receivers.

Like eq. ( 65) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], eq. ( 22) can be used as the basis of inverse problems: the seismic observations at its left-hand side are "inverted" to constrain unknown parameters contained in the theoretical formula at its right-hand side. Importantly, surfacewave phase velocity c(ω) can be determined from eq. ( 22) without knowledge of the factor |h(ω)| 2 P ρ (that is to say, of the power spectral density, surface density and intensity of the noise sources); Ekström et al. [2009] show that this amounts to identifying the values of c and ω for which the left-hand side of eq. ( 22) is zero (i.e., the "zero crossings" of the reconstructed Green's function). The factor |h(ω)| 2 P ρ becomes relevant if the Green's function's amplitude is to be accurately reconstructed, which is necessary if one wants to determine attenuation.

From noise cross-correlation to attenuation: the inverse problem

We show in the following how eq. ( 22) can be manipulated to formulate an inverse problem with α as unknown parameter, without neglecting |h(ω)| 2 P ρ. Let us start by expressing |h(ω)| 2 as a function of noise data.

It follows from eq. ( 17) that the power spectral density of the ambient signal recorded at a location x can be written

|s(x, ω)| 2 =|h(ω)| 2   N S j=1 G d 2D (x, x j , ω)e iωφ j   N S k=1 G d * 2D (x, x k , ω)e -iωφ k =|h(ω)| 2 N S j=1 |G d 2D (x, x j , ω)| 2 + N S j=1 N S k=1,k =j G d 2D (x, x j , ω)G d * 2D (x, x k , ω)e iω(φ j -φ k ) ; (23)
then, if one neglects cross terms based on the same arguments as above,

|s(x, ω)| 2 ≈|h(ω)| 2 N S j=1 |G d 2D (x, x j , ω)| 2 ≈|h(ω)| 2 N S j=1 |G d 2D (r j , ω)| 2 , (24)
where we have introduced the source-receiver distance r j =|x -x j |, to emphasize the fact that the value of G d 2D at a given point depends on its distance from the source, but not on the absolute locations of source and receiver. We next transform the sum at the right-hand side of eq. ( 24) into an integral over source-receiver distance; let us first replace the summation over sources with a summation over N D distance bins, i.e.,

|s(x, ω)| 2 ≈ |h(ω)| 2 N D k=1 N k |G d 2D (r k , ω)| 2 , ( 25
)
where N k denotes the number of sources at distances between r k and r k+1 from the receiver, and it is assumed that

G d 2D (r k , ω) ≈ G d 2D (r, ω)
as long as r k ≤ r ≤ r k+1 (which will be the case as long as the increment δr=r k+1 -r k is small). The area of the annulus centered at the receiver and bounded by the circles of radii r k and r k+1 is approximately 2πr k δr. It follows that N k ≈ 2πr k ρδr, and

|s(x, ω)| 2 ≈ 2πρ|h(ω)| 2 N D k=1 δr r k |G d 2D (r k , ω)| 2 ≈ 2πρ|h(ω)| 2 ∞ 0 dr r |G d 2D (r, ω)| 2 ≈ ρP 2 |h(ω)| 2 16c 4 ∞ 0 dr r H (2) 0 ωr c 2 e -2αr , (26)
where we have replaced G d 2D with its leading term, according to eq. ( 9). We have not been able to find a closed-form solution for the integral at the right-hand side of eq. ( 26); let us denote

I(α, ω, c) = ∞ 0 dr r H (2) 0 ωr c 2 e -2αr . (27)
Then, solving eq. ( 26)

for |h(ω)| 2 , |h(ω)| 2 ≈ 16c 4 ρP 2 I(α, ω, c) |s(x, ω)| 2 . ( 28
)
Eq. ( 28) stipulates that the power-spectral density of emitted ambient noise can be obtained from the power-spectral density of signal recorded at any receiver x, by application of a simple filter (provided that the surface density ρ and "intensity" P of noise sources are known). Since it was assumed that the function h is the same for all source-receiver vectors

x, the right-hand side of ( 28) can be replaced by an average over all available receivers, which we denote < |s(x, ω)| 2 > x :

|h(ω)| 2 ≈ 16c 4 ρP 2 I(α, ω, c) < |s(x, ω)| 2 > x ; (29)
this is irrelevant from a purely theoretical perspective, but useful when processing real data,

as averaging over all receivers will reduce effects that are not accounted for in our theoretical formulation, i.e. structural heterogeneities, dependence of the source time function on source location, nonuniformities in noise source distribution, etc.

Substituting (29) into ( 22), we find after some algebra that

s(x A , ω)s * (x B , ω) < |s(x, ω)| 2 > x ≈ c ω I(α, ω, c) 2 π J 0 ω|x A -x B | c e -α|x A -x B | α , (30)
where J 0 denotes the zeroth-order Bessel function of the first kind [e.g. [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]; importantly, the intensity P and surface density ρ of noise sources have canceled out, and only α and c remain to be determined. A nonlinear inverse problem can then be formulated, with α as its only unknown parameter. In practice, the dispersion curve c(ω) is first inverted for, ideally through a robust method that bypasses amplitude information and only accounts for phase [e.g. [START_REF] Cupillard | On the amplitude of surface waves obtained by noise correlation and the capability to recover the attenuation: a numerical approach[END_REF][START_REF] Kästle | Two-receiver measurements of phase velocity: cross-validation of ambient-noise and earthquake-based observations[END_REF]. We discuss in sec. 3.2 how a cost function to be minimized can then be introduced, based on eq. ( 30).

It might be noticed that a closed-form expression for α can be derived from the above treatment. Let us rewrite eq. ( 30) after replacing x A and x B with the locations of two other stations in our array, denoted x C and x D . We next divide eq. ( 30) by the equation so obtained, and find

s(x A , ω)s * (x B , ω) s(x C , ω)s * (x D , ω) ≈ J 0 (ω|x A -x B |/c) J 0 (ω|x C -x D |/c) e α(|x C -x D |-|x A -x B |) , (31)
which can be solved for α to obtain the sought formula

α(ω) = log [s(x A , ω)s * (x B , ω)] [J 0 (ω|x C -x D |/c)] [s(x C , ω)s * (x D , ω)] [J 0 (ω|x A -x B |/c)] 1 |x C -x D | -|x A -x B | , ( 32
)
where log denotes the natural logarithm. We have found that application of eq. ( 32) to our database does not lead to stable results, and, at the present stage, have not pursued this approach further. Eq. ( 32) might be of interest in the presence of a more diffuse ambient field, or a higher number of receivers allowing, e.g., for averaging over different azimuths as in [START_REF] Prieto | Anelastic earth structure from the coherency of the ambient seismic field[END_REF].

3 Application to Sardinian data set

Data cross correlation

Recordings of earthquakes are characterized by amplitudes much larger than those of truly diffuse, "ambient" signal, and can accordingly bias cross correlations [e.g., [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements[END_REF]. After subdividing seismic recordings into relatively short time intervals, some authors minimize this bias by identifying intervals where anomalously large displacements are recorded, to then exclude them from cross correlation. An alternative solution consists of cross-correlating separately the segments of seismic recording associated with each time interval; then, "partial" cross correlations so obtained can be normalized independently, usually by whitening, before being summed.

We follow here the latter approach, but, instead of whitening, normalize by the power spectral density < |s(x, ω)| 2 > x , as in the left-hand side of eq. ( 30). As explained in sec. 2.3.2, < |s(x, ω)| 2 > x is averaged over all stations x, and is proportional, through eq. ( 29), to the To verify that possible biases introduced by anomalous events are indeed suppressed by PSD-normalization, we also attempted to remove them directly from the data. We define as "event" the Fourier transform of one day of recording, and flag it as anomalous (an "outlier")

Figure 4: Power spectral density < |s(x, ω)| 2 > x (red line) used to normalize the observed cross-spectra as in eq. ( 30). The power spectral density obtained after removing from the data all 24 hour-long time intervals where anomalous events (sec. 3.1) were recorded (blue line) is also shown.

if the maximum amplitude exceeds a certain value. After testing various criteria to identify outliers, we decided to follow the Interquartile Range rule (IQR) [START_REF] Tukey | Exploratory data analysis[END_REF], with outlier constant set to 5, which we applied separately on 2-month-long subsets of the entire data set.

The power spectral density < |s(x, ω)| 2 > x obtained after removing the so defined outliers is also shown in Fig. 4 for the sake of comparison. We show in Fig. 5 the cross correlation of signals recorded over more than one year at several stations, before and after removing outliers, as described. After conducting the same test on several other station pairs, we conclude that the two approaches result in practically coincident results. We prefer PSDnormalization as it involves no arbitrary choices, e.g. in the definition of outlier. In Fig. 6 the results of PSD-normalizing and whitening the same cross correlations are compared.

Discrepancies are, again, very small, which confirms the stability of our results.

Figs. 5 and 6 also show that the imaginary parts of our observed cross correlations can be relatively large. This is in contrast with eq. ( 30), stipulating that the imaginary part of the frequency-domain cross correlation should be zero (or, equivalently, the causal and anticausal parts of the time-domain cross correlation should coincide [e.g. App. B of Boschi and Weemstra, 2015]). The same limitation affects many other seismic ambient-noise studies, where the conditions that the field be diffuse and that sources be uniformly distributed over space are not exactly met. To quantify the azimuthal bias of recorded ambient signal, we narrow-band-pass filter and inverse-Fourier-transform the frequency-domain cross-correlation associated with each station pair; we then compute, in the time domain, the signal-to-noise ratio (SNR) of both causal and anticausal parts of each cross correlation, defined as the ratio of the maximum signal amplitude to the maximum of the trailing noise [e.g. [START_REF] Yang | Characteristics of ambient seismic noise as a source for surface wave tomography[END_REF][START_REF] Kästle | Two-receiver measurements of phase velocity: cross-validation of ambient-noise and earthquake-based observations[END_REF]. The results are shown in Fig. 7, where it is apparent that, because our array is small, its azimuthal sampling is limited (sampling is particularly poor around the East-West direction). The sampled azimuths appear however to be characterized, on average, by similar values of SNR, at least at periods ≥5s, indicating a relatively isotropic ambient field: compare, e.g., with Fig. 12 of [START_REF] Kästle | Two-receiver measurements of phase velocity: cross-validation of ambient-noise and earthquake-based observations[END_REF], which was obtained by applying the exact same procedure to a larger and denser array.

In the following, we shall further reduce the unwanted effects of azimuthal bias, by implicitly averaging over all station pairs (and therefore all available azimuths) as only one model of α is sought that fits cross correlation data for all station pairs. Authors that process data from larger/denser arrays often also group station pairs in distance bins, and for each distance bin take an average over all azimuths [e.g. [START_REF] Prieto | Anelastic earth structure from the coherency of the ambient seismic field[END_REF][START_REF] Weemstra | Seismic attenuation from recordings of ambient noise[END_REF], but this is not feasible in our case owing to the limited size of our array.

Dispersion and attenuation parameters

For each station pair i, j, a dispersion curve c ij =c ij (ω) is derived via the frequency-domain method of Ekström et al. normalization. Importantly, in both cases the frequency range over which the dispersion curve is defined changes depending on the station pair; as a general rule, it is hard to constrain its low-frequency end if stations are relatively close to one another.

Attenuation parameter as a scalar constant

Taking the squared modulus of the difference of the left-and right-hand sides of eq. ( 30), and summing over all frequency samples ω k and station pairs i, j, the cost function

i,j k s(x i , ω k )s * (x j , ω k ) < |s(x, ω k )| 2 > x - 2c ij (ω k ) ω k √ 2π 1 I[α, ω k , c ij (ω k )] J 0 ω k |x i -x j | c ij (ω k ) e -α|x i -x j | α 2 (33)
is obtained.

The right-hand side of eq. ( 30) is, through the Bessel function J 0 , an oscillatory function of ω. The value of the attenuation parameter α, however, only affects its envelope, and not its oscillations, with respect to ω. Following other authors who estimated attenuation on the basis of ambient-noise cross correlation [e.g., [START_REF] Prieto | Anelastic earth structure from the coherency of the ambient seismic field[END_REF], we accordingly define the Figure 9: Cost C 1 as defined by eq. ( 34), as a function of the scalar attenuation parameter α.

cost function

C 1 (α) = i,j k env s(x i , ω k )s * (x j , ω k ) < |s(x, ω k )| 2 > x -env 2c ij (ω k ) ω k √ 2π 1 I[α, ω k , c ij (ω k )] J 0 ω k |x i -x j | c ij (ω k ) e -α|x i -x j | α 2 , ( 34
)
where env denotes the envelope function, which we implement by fitting a linear combination of splines [e.g. [START_REF] Press | Numerical Recipes in Fortran[END_REF] to the maxima of the absolute value of its argument.

The attenuation parameter α in this case is a scalar value, independent of both frequency and location. We show in Fig. 9 how C 1 varies as a function of α; a clear minimum is identified at α=3.03 × 10 -5 m -1 , which is our preferred attenuation model in this scenario. We next use this value for α to evaluate numerically the right-hand side of eq. ( 30), and compare the results with the normalized cross-correlation at the left-hand side; this is illustrated in Fig. 10 for four different station pairs. The modeled Green's functions have their zeroes at approximately the same frequencies as the measured ones, indicating that dispersion curves derived as described above are reliable. At short interstation distances, the found value of α also results in a relatively good fit of observed amplitude at most frequencies. The fit deteriorates with increasing distance, indicating that the assumption that α be constant does not honour the actual complexity of the medium.

As an additional test, we found the cost function defined by expression (33) (i.e., no envelopes are taken) to be similar to C 1 , with a less prominent but well defined minimum at α=2.75 × 10 -5 m -1 . The similarity of this estimate of α with the one based on function C 1

C 2 (α,ω)-min[C 2 (α,ω)] max[C 2 (α,ω)]-min[C 2 (α,ω)]
, where min[C 2 ] and max[C 2 ] denote the minimum and maximum values of C 2 for all sampled values of α and ω. The stepwise trend of the minima of C 2 is correlated with the stepwise growth (also as a function of ω) of the number of station pairs for which cross-correlation data are available.

suggests that this result is robust.

Frequency-dependent attenuation parameter

We next allow the attenuation parameter α to change as a function of ω; in practice, we evaluate the cost function

C 2 (α, ω) = i,j env s(x i , ω)s * (x j , ω) < |s(x, ω)| 2 > x -env 2c ij (ω) ω √ 2π 1 I[α(ω), ω, c ij (ω)] J 0 ω|x i -x j | c ij (ω) e -α|x i -x j | α 2 , (35)
shown in Fig. 11, after normalization, as a function of both α and ω. Since both terms within the square brackets in eq. ( 35) are close to a Bessel function of ω, it is not surprising that their difference has an oscillatory behaviour with respect to ω; because only a discrete and limited set of interstation distances are available from our data set, this effect is not canceled by summation, as is apparent from Fig. 11. Fig. 11 also shows that C 2 (α, ω) has a single, well-defined minimum at all frequencies, resulting in the α(ω) curve of Fig. 12.

The actual values of the minima of C 2 (α, ω), without normalization, are shown in Fig. 12.

C 2 decreases with growing ω, meaning that relatively high frequencies are better fit than low frequencies. Fig. 13 also shows that the amplitude fit between observed cross correlations and modeled Green's functions is worse for large interstation distances.

In analogy with sec. 3.2.1, we also evaluated an alternative cost function, where the difference of observed and theoretical, normalized cross correlation is computed without extracting The corresponding values of α are similar to those obtained based on C 2 ; we do not show or discuss them here in the interest of brevity.

Cost function as a weighted sum

Attenuation models of secs. 3.2.1 and 3.2.2 achieve a systematically worse fit of cross correlations between faraway as opposed to nearby stations (see examples in Figs. 10 and13).

To some (minor) extent, this bias might stem from the error involved in the far-field/highfrequency approximation discussed in sec. 2.1, causing a fictitious loss of amplitude of the theoretical cross correlation at large interstation distances (Fig. 1). More importantly, it might result from the fact that, by geometrical spreading, cross-correlation amplitude decreases with growing interstation distance; assuming the relative misfit on cross-correlation amplitude to be independent of distance, pairs of faraway stations are then systematically associated with smaller absolute errors and thus contribute less to the cost functions C 1 and C 2 . We attempt to reduce this effect by replacing the sum in C 2 with a weighted sum,

C 3 (α, ω) = i,j w(|x i -x j |) env s(x i , ω)s * (x j , ω) < |s(x, ω)| 2 > x -env 2c ij (ω) ω √ 2π 1 I[α(ω), ω, c ij (ω)] J 0 ω|x i -x j | c ij (ω) e -α|x i -x j | α 2 , (36)
where w(|x i -x j |) = |x i -x j | e and e is the Euler number. We selected this weighting scheme after a suite of preliminary tests, where the weight w was chosen to coincide in turn with different powers (from square root to fourth power) of interstation distance. We show C 3 (α, ω) in Fig. 14, and the corresponding best-fitting values of α in Fig. 15. Interestingly, it is apparent from Fig. 15 that the so obtained function α(ω) spans a smaller range of values than its counterpart discussed in sec. 3.2.2; α is also generally smaller, resulting in larger amplitude of modeled cross correlations (Fig. 16) at all interstation distances.

We are unable to determine a unique function α(ω) that results in a comparably good fit of cross-correlation amplitude for all station pairs: observed amplitudes tend to be overestimated by our "model" at short interstation distances, and underestimated at large interstation distances. This effect might point to possible lateral heterogeneities of α that our data set is too limited to constrain; it could also be associated with the error inherent in ambientnoise-based reconstruction of the Green's function, when the seismic ambient field (as in most practical applications) is not truly diffuse.

Discussion and conclusions

The main purpose of this study was to clarify some aspects of the relationship between the cross correlation of seismic ambient noise and surface-wave attenuation (attenuation parameter α or quality factor Q). It is known that this relationship is complicated by the need to process ambient-noise cross correlation data so as to reduce as much as possible the bias introduced by anomalous high-amplitude events (earthquakes). This is often achieved by subdividing continuous traces into shorter time intervals, which are then whitened and crosscorrelated separately before being "stacked" together; but whitening has a complicated effect on the mentioned noise-attenuation relationship [START_REF] Weemstra | On estimating attenuation from the amplitude of the spectrally whitened ambient seismic field[END_REF]. We develop here a different normalization criterion, with practical effects similar to whitening, but obtained by simply manipulating the reciprocity theorem without any additional data processing. This results in the relationship (30), which is strictly valid provided that sources of seismic ambient noise be uniformly distributed over R 2 , that their phases be random and uncorrelated, and that the spectrum of noise sources be spatially uniform (sec. 2.3.2). Our experimental setup, consisting of a small array deployed on an island, is chosen to approximately satisfy these requirements.

Eq. ( 30) involves a proportionality factor, relating normalized cross-correlation and the product J 0 (ω∆/c)e -α∆ (with ∆ denoting interstation distance), that had been neglected in previous studies. Compare, e.g., with eq. ( 7) of [START_REF] Prieto | Anelastic earth structure from the coherency of the ambient seismic field[END_REF] or eq. ( 1) of Lawrence and

Prieto [2011]. The need to account for such a factor was pointed out theoretically by Tsai

[2011], while [START_REF] Harmon | Distribution of noise sources for seismic interferometry[END_REF] and [START_REF] Weemstra | Seismic attenuation from recordings of ambient noise[END_REF] introduced it as a free parameter of their inversions. Similar to [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF], we have derived an analytical expression for it, and evaluate it numerically based on our estimates for α and phase velocity c(ω). Fig. 17 shows that the factor in question is of the order of unity, which would explain the success of [START_REF] Prieto | Anelastic earth structure from the coherency of the ambient seismic field[END_REF], [START_REF] Lawrence | Attenuation tomography of the western United States from ambient seismic noise[END_REF] and others in inferring reasonable values for α.

On the basis of eq. (30), we explored several possible definitions of cost function (secs.

3.2.1 through 3.2.3), quantifying the misfit between observed and modeled cross-correlation amplitudes. We first assumed the attenuation parameter α to be constant, independent of frequency and position. We next identified a frequency-dependent α = α(ω) model that minimizes the misfit for all receiver-receiver pairs at the same time. Since the cost function involves a sum over station pairs, we finally introduced a cost function where the contribution of each station pair was weighted differently depending on interstation distance.

To compare quantitatively the misfit achieved by different models, let us introduce the [START_REF] Tukey | Exploratory data analysis[END_REF] of the distributions at (a), (b) and (c). misfit function

m ij = k s(x i , ω k )s * (x j , ω k ) < |s(x, ω k )| 2 > x - 2c ij (ω k ) ω k √ 2π 1 I[α(ω k ), ω k , c ij (ω k )] J 0 ω k |x i -x j | c ij (ω k ) e -α(ω k )|x i -x j | α(ω k ) 2 , (37)
associated to the station pair ij. We implement eq. ( 37) for each model (corresponding to the cost functions C 1 through C 3 ) and for each station pair ij, and visualize the resulting values of m ij in Fig. 18, in the form of one histogram per model.

It is apparent from Fig. 18, and could be anticipated from a visual comparison of Fig. 10 with Fig. 13, that allowing α to vary with respect to ω results in an overall improvement of fit with respect to constant-α models. Minimizing the cost function C 3 , on the other hand, results in an increase in the misfit m ij with respect to C 2 , as nearby station pairs tend to contribute to m ij more than faraway ones (see discussion of the cost function C 2 in sec.

3.2.3). We found by a Kolmogorov-Smirnoff test [e.g. [START_REF] Press | Numerical Recipes in Fortran[END_REF] that the probability that the three histograms in Fig. 18 correspond to the same statistical distribution is always 3%, and always 1% for the histogram in Fig. 18b (C 2 in Fig. 18d). This suggests that the improvement in fit achieved by the cost function C 2 is significant. Only a limited number of samples (station pairs) is available, however, and a more reliable statistical analysis should be conducted in the future on a larger database. In addition, the level of complexity (number for example

∇ 2 1 f - 1 c 2 ∂ 2 f ∂t 2 = 0 (A.
(x, t) = 1 D t 0 dζ f (x, t -ζ; ζ) = c 2 P t 0 dζ R 2 d 2 x G 2D (x, x , t -ζ)η(x , ζ) (A.19)
solves the general inhomogeneous problem (A.7)-(A.9) for any forcing term η(x, t).

A.4 Green's problem as inhomogeneous equation

We next consider the membrane-displacement eq. ( 5); it is convenient to choose the forcing term U (0, ω)f (x 1 , x 2 , ω) = -iωF δ(x -x S ) so that q = F δ(x -x S ) in sec. 2.2. In the time domain the resulting equation reads In the literature, problem (A.7)-(A.9) with an impulsive forcing term is often presented as the Green's problem, and its solution (A.21), (A.22) as the Green's function. It is important to realize that this is not mathematically the same as our Green's function G 2D . Other definitions of Green's function likewise exist. In practice, the phrase "Green's function" always refers to the response of a medium to an impulsive excitation; but this definition is ambiguous, and, for any given medium, the "impulsive response" can be defined mathematically in a virtually unlimited variety of ways.

∇ 2 1 u(x, t) - 1 c 2 ∂ 2 ∂t 2 u(x, t) = -F δ(x -x S )δ (t), (A.
Appendix B Approximate expression for the damped-membrane

Green's function

We show in the following that, as long as attenuation is weak, the lossy-membrane Green's function

G d 2D (x 1 , x 2 , ω) = - iP 4 √ 2πc 2 H
(2) 0

x ω 2 c 2 -2iαω c (copy of eq. 8)

can be approximated by the product of the lossless Green's function G 2D and a term that decays exponentially with source-receiver distance. Two different arguments are provided in the next two sections.

B.1 Stationary-phase approach

If z is small, √ 1 + iz ≈ 1 + iz/2 is accurate to first order in z. It is then reasonable to write the argument of H

(2) 0 in eq. ( 8)

x ω 2 c 2 - 2iαω c = ωx c 1 - 2iαc ω = ωx c (1 -εi), (B.1)
where ε ≈ cα/ω 1.

Substituting expression (B.1) into the integral representation of H

(2) 0 , e.g. eq. ( 10.9.11) of DLMF,

H (2) 0 x ω 2 c 2 - 2iαω c = i π +∞ -∞ dt e -i cosh(t)[ ωx c (1-iε)] = i π +∞ -∞
dt e -i ωx c cosh(t) e -ε ωx c cosh(t) .

(B.2)

V. Tsai (personal communication, 2013) first derived eq. (B.2), and solved the integral on its right-hand side via the stationary-phase approximation. The right-hand side of (B.2) has indeed the same form as the one-dimensional stationary-phase integral, e.g. eq. (A1) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], except for the fact that the non-oscillatory term e -ε ωx c cosh(t)

depends not only on the integration variable t, but also on x. Because this term is small in the vicinity of the (only) stationary point t=0, we assume that the stationary-phase formula still applies; considering that cosh coincides with its second derivative, and that cosh(0) = 1, it follows from eq. (A2) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF] that +∞ -∞ dt e -i ωx c cosh(t) e -ε ωx c cosh(t) ≈ 2 e -ε ωx c e -i( ωx c -π 4 ) πc 2ωx , (B.3) with the factor 2 at the right-hand side to account for the fact that the stationary point t=0

is not one of the integration limits. Consequently,

H

(2) 0

x ω Let us emphasize, again, that this simplification relies not only on the weak-scattering approximation, but also on the far-field approximation for H

(2) 0 .

Figure 2 :

 2 Figure 2: (a) Geographical locations of receivers (black triangles with station names) over the island of Sardinia. (b) Distribution of interstation distances for all station pairs in our deployment; the mean and median of the distribution are 114.45 km, 104.95 km, respectively. Acronyms starting with the letters UT identify stations deployed by our team.

  At the end of June, 2016, our team has deployed an array of broadband seismic stations (Trillium Nanometrics 120s posthole broadband stations) around Sardinia, as shown in Fig.2a. This temporary deployment was complemented by three permanent stations belonging to the Italian MN and IV networks. Except for UT001 and UT011, stations recorded continuously for 24 months. Station UT001 recorded from June 2016 until November 2017; it was then removed and redeployed at location UT011, where it recorded November 2017 to September 2018. We next explain how ambient recordings of displacement(vertical component only) were cross correlated to one another, to determine first a set of Rayleigh-wave dispersion curves, and then, according to sec. 2.3.2, the attenuation parameter α.

Figure 3 :

 3 Figure 3: Five examples of power spectral density < |s(x, ω)| 2 > x , computed for five consecutive six-hour long intervals. The computation is repeated for each season, with panels a, b, c and d showing the results obtained in summer, autumn, winter and spring, respectively. Recordings used for this examples start on (a) July 8, 2017 at 8:49AM; (b) October 5, 2017, at 2:49AM; (c) January 20, 2018, at 2:49AM; April 14, 2018, at 8:49PM.

  actual power spectral density |h(ω)| 2 of ambient noise. 6-hour long non-overlapping time windows are normalized independently before being summed. We shall refer to this procedure as PSD-(power spectral density) normalization. Examples of the factor < |s(x, ω)| 2 > x are shown in Fig. 3 for twenty different time windows, sampling all four seasons. The cumulative power spectral density < |s(x, ω)| 2 > x is shown in Fig. 4.

Figure 5 :

 5 Figure 5: Real (a, c, e) and imaginary (b, d, f) parts of the cross correlations of the entire available recordings at stations UT.006 and UT.009 (a, b), IV.AGLI and IV.DGI (c, d), UT.002 and UT.003 (e, f). The associated interstation distances are 231, 97 and 57 km, respectively. Cross correlations of the entire recorded traces are PSD-normalized (red lines); alternatively (blue), outliers are removed from the traces prior to cross-correlating, as discussed in sec. 3.1. PSD-normalization is preferred throughout the rest of this study.

Figure 6 :

 6 Figure6: Same as Fig.5, but PSD-normalized (red lines) and whitened cross correlations (blue) are now compared. In both cases, the entire available seismic records are cross correlated, without attempting to identify and remove outliers. For each interstation distance ∆, the Bessel function J 0 (ω∆/c) (yellow), to which the real parts of cross correlations should be proportional according to (30), is also shown.

Figure 7 :

 7 Figure7: Estimates of the SNR of our cross correlations, as a function of station-pair azimuth, at periods of 3 to 15s, as specified. The value of SNR determines the length of the red segments, while their orientation coincides with the station-pair azimuth, with 0 • corresponding to the North, 90 • to the East, etc. The signal-to-noise ratios associated with both the causal and anticausal parts of the cross correlations are computed and plotted; for each station pair, the causal and anticausal segments point in opposite directions. In practice, longer segments should point to the direction where most seismic ambient signal comes from.

Figure 8 :

 8 Figure 8: Phase velocity curves retrieved from the whitened (red lines) and PSD-normalized (blue) cross-spectra, for station pairs (a) UT.006 -UT.009 (interstation distance ∼ 231 km), (b) UT.002 -UT.004 (∼ 152 km), (c) IV.AGLI -IV.DGI (∼ 97 km), and (d) UT.002 -UT.003 (∼ 57 km).

Figure 10 :

 10 Figure 10: Comparison of normalized data (red lines) and model (blue), i.e. left-and righthand side of eq. (30), after substituting α=3.03 × 10 -5 m -1 , as explained in sec. 3.2.1 (i.e., inversion via the cost function C 1 ). As in Fig. 8, panels a, b, c and d correspond to station pairs UT.006-UT.009, UT.002 -UT.004, IV.AGLI-IV.DGI and UT.002-UT.003, respectively, with interstation distances decreasing from 231 to 57 km.

Figure 11 :

 11 Figure 11: Cost function C 2 (α, ω) (sec. 3.2.2) shown (after normalization) as a function of the attenuation parameter α and frequency ω. We normalize C 2 according to the formula

Figure 12 :

 12 Figure 12: Attenuation parameter α (red dots, scale on the left) and corresponding values of the cost function C 2 (α, ω) (blue dots, scale on the right), both plotted as functions of frequency ω.

Figure 13 :

 13 Figure 13: Same as Fig. 10, but the blue curves are obtained by substituting into eq. (30) the values of α(ω) obtained by minimizing the cost function C 2 of sec. 3.2.2.

Figure 14 :

 14 Figure 14: Same as Fig. 11, but the cost function C 3 is evaluated, where contributions of different station pairs are weighted differently (sec. 3.2.3) according to interstation distance.

Figure 15 :

 15 Figure 15: Same as Fig. 12, but the values α(ω) (red dots) that minimize at each ω the cost function C 3 (sec. 3.2.3), and the corresponding values of C 3 (blue dots) are shown.

Figure 16 :

 16 Figure 16: Same as Figs. 10 and 13, but the blue curves are obtained by substituting into eq. (30) the values of α(ω) obtained by minimizing the cost function C 3 (sec. 3.2.3).

Figure 17 :

 17 Figure 17: Proportionality factor

Figure 18 :

 18 Figure 18: Number (vertical axis) of station pairs for which the misfit m ij falls within a given interval (horizontal axis), for models of α resulting from the minimization of cost functions (a) C 1 , (b) C 2 , and (c) C 3 . (d) Box plots[START_REF] Tukey | Exploratory data analysis[END_REF] of the distributions at (a), (b) and (c).

  20)where δ (t) denotes the derivative of the delta function. Assuming zero inital displacement and velocity, eq. (A.20) is a particular case of problem (A.7)-(A.9). A solution is found bysubstituting η(x, t) = -F δ(x -x S )δ (t) into eq. (A.19), G 2D (x, x , t -ζ)δ(x -x S )δ ζ) G 2D (x, x S , t -ζ) = F c 2 P ∂ ∂t G 2D (x, x S , t).

  the frequency domain, this maps tou(x, ω) = iω F c 2 P G 2D (x, x S , ω). (A.22)

  the leading term at the right-hand side coincides with the right-hand side of eq. (B.5).

  Importantly, this approximation is only valid for large x, i.e. the presence of sources in the near field of the receivers, which is necessary to reconstruct the Green's function according to sec. 2.2, introduces a possibly significant error.Making use, again, of the Taylor expansion√ 1 + iz ≈ 1 + iz/2 in the argument of G d 2D ,If one substitutes into (B.6) the far-field and/or high-frequency approximation for H Expression (B.7) can be simplified if one considers that, for small z, (1iz) -1 2 ≈ 1 + iz/2is accurate to first order in z; after so expanding the square root,

	2 c 2 -	2iαω c		≈	2c πωx	c -π 4 ) e -αx e -i( ωx	(B.4)
							≈ H 0 (2)	ωx c	e -αx ,
	where we have replaced ε with cα/ω. Substituting eq. (B.4) into (8),
	G d 2D (x 1 , x 2 , ω) ≈ -	iP √ 4 2πc 2	H 0 (2)	ωx c	e -αx	(B.5)
	G d 2D (x 1 , x 2 , ω) = -	iP √ 4 2πc 2	H 0 (2)	ωx c	1 -	2iαc ω	(B.6)
		≈ -	iP √ 4 2πc 2	H 0 (2)	ωx c	-iαx .
								(2) 0 [e.g.
	Boschi and Weemstra, 2015, eq. (C5)],
	G d 2D (x 1 , x 2 , ω) ≈ -	4	iP √ 2πc 2	2 c -iα x π ω	e -i( ωx c -π 4 ) e -αx	(B.7)
	G d 2D (x 1 , x 2 , ω) ≈ -	iP √ 4 2πc 2		1 + i	αc 2ω	2 π ωx

B.2 Taylor-series approach

Acknowledgments

Our many exchanges with Julien de Rosny, Emanuel Kaestle, Piero Poli, Victor Tsai and Kees Weemstra were very beneficial to this study. We are also grateful to Nicholas Harmon, one anonymous reviewer and the editor Ana Ferreira for carefully reading and commenting our manuscript. LB is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement 641943 (WAVES network)

remains to be determined; it is beyond the scope of our current contribution.

It would be useful to compare our observations with independent estimates of α in the region of interest, but, to our knowledge, no studies of surface-wave attenuation in Sardinia have been published so far. Comparison with global surface-wave attenuation literature suggests that our estimates of α (on the order of 10 -5 m -1 in the period range 2-10s according to Figs. 12 and 15) are relatively large. Studies based on earthquake data suggest values of α on the order of 10 -6 m -1 , but quickly growing, with decreasing period, within the period range of interest to us [START_REF] Mitchell | Anelastic structure and evolution of the continental crust and upper mantle from seismic surface wave attenuation[END_REF][START_REF] Romanowicz | Deep earth structure: Q of the earth from crust to core[END_REF]. Surface waves at those periods are perhaps best sampled by ambient-noise cross correlations; Prieto et al.

[2009] find α≈6.4 × 10 -6 m -1 from seismic ambient noise at 5s period, consistent with the laterally-varying values obtained by [START_REF] Lawrence | Attenuation tomography of the western United States from ambient seismic noise[END_REF], and not far from the values proposed here; [START_REF] Lin | On the reliability of attenuation measurements from ambient noise cross-correlations[END_REF] estimate α≈1 × 10 -6 m -1 or lower. Those studies neglect the proportionality factor shown here in Fig. 17, which might partially explain the discrepancy.

Both [START_REF] Harmon | Distribution of noise sources for seismic interferometry[END_REF] and [START_REF] Weemstra | Seismic attenuation from recordings of ambient noise[END_REF] account for this effect, although in a different way than done here; [START_REF] Harmon | Distribution of noise sources for seismic interferometry[END_REF] find α≈5 × 10 -7 m -1 at 7.5s period, while [START_REF] Weemstra | Seismic attenuation from recordings of ambient noise[END_REF] obtain estimates of α actually larger than ours. [START_REF] Viens | Retrieving impulse response function amplitudes from the ambient seismic field[END_REF] likewise fit ambient-noise surface-wave data with α=1.4 × 10 -5 m -1 in the period range 3-10s, consistent with our estimates.

As values of α obtained from different methods are compared, one should keep in mind the significant uncertainties associated with the many practical issues quantified, e.g., by [START_REF] Menon | On the apparent attenuation in the spatial coherence estimated from seismic arrays[END_REF]. Estimates of surface-wave attenuation might be affected by the presence, in seismic ambient noise, of body-wave signal not accounted for by the theory [e.g. [START_REF] Gerstoft | PKP wave microseisms observed from distant storms[END_REF]. Differences in the terrains where the data were collected also play a role. This preliminary application, limited to a small data set, demonstrates that our new algorithm leads to reasonable estimates of α. Future applications to denser instrument arrays, with a more thorough account of heterogeneity in source distribution, are likely to benefit more significantly from the theoretical improvements that we have introduced.

Appendix A Green's functions of the scalar wave equation (homogeneous lossless media)

We describe in the following two definitions of the Green's function that are commonly found in the literature. In one case (sec. A.1), the Green's function is obtained by prescribing nonzero initial velocity at the source; initial displacement is zero and no other forcing is present. In another case (sec. A.4), both initial displacement and velocity are everywhere zero, but an impulsive force is applied at the source. Following [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF],

we adhere throughout this study to the former definition, but in sec. 2.2 implicitly make use of the latter. Through the mathematical tools provided in secs. A.2 and A.3, a relationship between the two definitions is obtained; this relationship is employed in sec. 2.2.

A.1 Green's problem as homogeneous equation

In analogy with [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], we call Green's function G 2D = G 2D (x, x S , t) (with t denoting time) the solution of

with initial conditions

i.e. an impulsive source at x S . Only "causal" solutions, vanishing when t < 0, are relevant.

The scalar constant P serves to remind us of the physical dimensions of G 2D ; for instance, one can think of (A.1) as the displacement equation for a membrane, in which case ∂G 2D ∂t (x, x S , 0)

is the initial velocity, and P has dimensions of cubed distance over time (recall that, in two dimensions, δ has dimensions of one over squared distance).

Boschi and Weemstra [2015] show that, in the time domain, the "Green's problem" (A.1)-

where H denotes the Heaviside function. This corresponds to

in the frequency domain, with H

(2) 0 denoting the zeroth-order Hankel function of the second kind. When ωx/c 1, the expression (A.5) can be replaced by the far-field/high-frequency approximation

It might be noticed upon comparing our expression (A.5) for G 2D (x, ω) with that of, e.g., [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF], that the membrane-wave Green's function given in that study is proportional to the zeroth-order Hankel function of the first kind: this apparent discrepancy is explained by the fact that the Fourier-transform convention assumed by [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF] differs from ours (compare eq. ( 4) of [START_REF] Tsai | Understanding the amplitudes of noise correlation measurements[END_REF] and eq. (B2) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF], and consider eq. ( E16) of [START_REF] Boschi | Stationary-phase integrals in the cross-correlation of ambient noise[END_REF]).

A.2 Duhamel's principle

Consider the initial-value problem

with η an arbitrary forcing term. Suppose that a solution v(x, t) to the following, similar homogeneous problem can be found: This result is known as Duhamel's principle [e.g. [START_REF] Hildebrand | Advanced Calculus for Applications[END_REF][START_REF] Strauss | Partial Differential Equations[END_REF].

A.3 General initial condition

Once the Green's function associated with a given differential equation is found, it can be used to solve rapidly more general initial-value problems based on the same equation. Consider


Emmanuel Franck

Laurent Navoret

Semi-implicit two-speed Well-Balanced relaxation scheme for Ripa model

Keywords: relaxation method, shallow-water model, semi-implicit scheme

In this paper, we propose a semi-implicit well-balanced scheme for the Ripa model based on a two-speed relaxation. The method both preserves equilibria and has an implicit step that reduces to the inversion of a constant Laplacian. Numerical simulations show that the scheme well capture low-Froude flows.

Introduction

To discretize two-scale hyperbolic problems with a good accuracy, several methods have been proposed and are generally based on semi-implicit schemes to prevent us from the fast scale stringent stability condition. In [START_REF] Coulette | Implicit time schemes for compressible fluid models based on relaxation methods[END_REF]5], such a scheme have been proposed and it is based on very recent relaxation method and a dynamical splitting. This method allows to adapt the time discretization to the regime and the implicit step just reduces to the inversion of a discrete Laplacian with constant coefficient. In this paper, we adapt this numerical scheme to the Ripa model with topography, which describes shallow-water flows with horizontal temperature gradients. In the Ripa model, the fast scale dynamics are perturbative gravity waves around an equilibrium and the slow scale dynamics is the convection. Here equilibria are balance between the pressure gradient and the topography source term. In addition to the implicit treatment of the perturbative waves, the scheme has to preserve the equilbria so to prevent the generation of spurious waves around them. This is the so-called wellbalanced property.

We consider the one-dimensional Ripa model [START_REF] Berthon | Well-balanced scheme to capture non explicit steady-states: Ripa model Mathematics of Computation[END_REF]:

Emmanuel Franck, Laurent Navoret INRIA Grand-Est and IRMA Strasbourg, 7 Rue René Descartes, 67000, Strasbourg, France, e-mail: emmanuel.franck@inria.fr, laurent.navoret@math.unistra.fr

1        ∂ t h + ∂ x (hu) = 0, ∂ t (hu) + ∂ x (hu 2 + p(h, Θ)) = -ghΘ∂ x z, ∂ t (hΘ) + ∂ x (hΘu) = 0, (1)
where h(x, t) is the water height, u(x, t) the velocity, Θ(x, t) the temperature and z(x) the topography and the pressure law is given by: p(h, Θ) = gΘ 1 2 h 2 , with g the gravity constant. This system is hyperbolic and has three characteristic speeds: Σ = {uc, u, u + c}, with c = ghΘ. We introduce the Froude number Fr = u/c, which describes the ratio between the two velocity scales. We are interested into perturbations of stationary solutions. These solutions are balance between pressure force and source term: ∂ x p = -ghΘ∂ x z. Hereafter we consider the following three families of equilibria, with zero velocity:

       u = 0, Θ = cst, h + z = cst,          u = 0, z = cst, Θ h 2 2 = cst,        u = 0, h = cst, z + h 2 ln(Θ) = cst.
(

) 2
The first aim is to write a scheme that preserves these steady states. Indeed, if the scheme does not preserve the steady states, spurious velocity and pressure modes would appear and would destroy the accuracy of the scheme for small velocity (low Froude regime). The second aim is to capture the dynamics near an equilibrium with an acceptable cost. For example, we consider the following perturbation O(Fr) of the first steady equilibria with Fr 1. In that case, the perturbation has a small amplitude but moves with a large propagation speed of order O(1/Fr). Therefore, implicit schemes are usually required to filter these small fast waves.

Two-speed relaxation system

To simplify the implicit treatment of the dynamics, we propose a relaxation model that linearizes the fast scale associated to the gravity waves. We introduce two additional unknows Π(x, t) and v(x, t) and consider the following extended hyperbolic system

                   ∂ t h + ∂ x (hv) = 0, ∂ t (hu) + ∂ x (huv + Π) = -gh∂ x z, ∂ t (hΘ) + ∂ x (hΘv) = 0, ∂ t Π + v ∂ x Π + h m λ 2 ∂ x v = 1 ε p(h, Θ) -Π ∂ t v + v ∂ x v + 1 h m ∂ x Π = -h h m gΘ∂ x z + 1 ε v -u (3)
with h m > 0 and λ > 0 are constant relaxation parameters and where ε > 0 is the relaxation parameter. This system is an approximation of (1) in the limit ε tends to zero. This result can be shown formally.

Proposition 1 As ε → 0, the relaxation system ( 3) is consistent at first order in ε with

       ∂ t h + ∂ x (hu) = ε ∂ x ( β(∂ x p + gh∂ x z)), ∂ t (hu) + ∂ x (hu 2 + Θ g 2 h 2 ) = -hgΘ∂ x z + ε ∂ x (u β(∂ x p + gh∂ x z)) + ε∂ x (γ∂ x u) ∂ t (hΘ) + ∂ x (hΘu) = ε∂ x (Θ β(∂ x p + gh∂ x z)) , ( 4
) with β = h h m -1 , γ = h m λ 2 -hc 2 .
The proof is based on a classical Chapman-Enskog expansion detailed for the Euler system case in [5]. Note that an equilibrium given by ∂ x p = -ghΘ∂ x z and u = 0 is still a steady state of the first order approximation (4) of the relaxation system. It is not the case for all the relaxation models like the Jin-Xin relaxation used in [START_REF] Coulette | Implicit time schemes for compressible fluid models based on relaxation methods[END_REF]. This property is necessary, but not sufficient, to obtain a well-balanced scheme since we will discretize the relaxation system (3) and not the original one. In [5], the entropy stability of ( 4) is analysed. This computation adapted to our model gives the following stability conditions: β ≥ 0, γ ≥ 0.

Semi-implicit scheme

The structure of the relaxation system (3) enables us to devise a semi-implicit scheme with a simple implicit part. This scheme is based on a splitting method between the different time-scale dynamics contained in (3). We split system (3) into the convection part (C), the gravity waves part (W) and the relaxation part (R):

(C)                    ∂ t h + ∂ x (hv) = 0, ∂ t (hu) + ∂ x (huv + F 2 Π) = -F 2 ghΘ∂ x z, ∂ t (hΘ) + ∂ x (hΘv) = 0, ∂ t Π + v∂ x Π + h m λ 2 ∂ x v = 0 ∂ t v + v∂ x v + F 2 h m ∂ x Π = -F 2 h h m gΘ∂ x z (W )                  ∂ t h = 0, ∂ t (hu) + (1 -F 2 ) (∂ x Π + hg∂ x z) = 0, ∂ t hΘ = 0 ∂ t Π + (1 -F 2 )h m λ 2 ∂ x v = 0 ∂ t v + (1 -F 2 ) 1 h m ∂ x Π + h h m g∂ x z = 0 (R) ∂ t Π = 1 ε (p(h, Θ) -Π) , ∂ t v = 1 ε (u -v) ,
where

F = max F min , min u √ hΘg
, 1 is an estimation of the global Froude number.

Introduced in [START_REF] Iampietro | A Mach-Sensitive Implicit-Explicit Scheme Adapted to Compressible Multi-scale Flows[END_REF], this dynamic splitting allows to adapt itself to the dynamics. In practice, the relaxation step (R) is treated as a projection: Π = p(h, Θ) and v = u. The key point is the discretization of the waves part (W) with an implicit solver. First, we note that ∂ t h = 0 so that the topography source term can be treated explicitly. The construction of the implicit scheme is based on the following remark: the equations on v and Π form a linear independent system. We can thus discretize all the equations of (W) with an implicit Euler scheme and then we get an implicit elliptic problem on Π n+1 by plugging the expression of v n+1 in the Π equation. This elliptic problem is linear and has constant coefficient. After solving this problem, we obtain Π n+1 and then v n+1 and (hu) n+1 can be computed with an explicit cost. For the spatial discretization, we use a classical finite volume scheme for the elliptic part and the central flux for the first derivative. The final algorithm writes:

• Step 1: solve Π n+1 j -(1 -F 2 ) 2 ∆t 2 λ 2 Π n+1 j+1 -2Π n+1 j + Π n+1 j-1 ∆x 2 = Π n j -∆t(1 -F 2 )λ 2 v n j+1 -v n j-1 2∆x (1 -F 2 ) 2 ∆t 2 λ 2 1 ∆x S n j+ 1 2 -S n j-1 2 , with S n j+ 1 2 = h n j+ 1 2 Θ n j+ 1 2 z j+1 -z j ∆x , (5)
where the quantities h j+ 1 2 and Θ j+ 1 2 will be define below.

• Step 2: compute v n+1 j = v n j -(1 -F 2 ) ∆t h m Π n+1 j+1 -Π n+1 j-1 2∆x -(1 -F 2 ) ∆t h m g 2 S n j+ 1 2 -S n j-1 2 , (hu) n+1 j
= (hu) n j -∆t(1 -F 2 ) Π n+1 j+1 -Π n+1 j-1 2∆x - g∆t 2 (1 -F 2 ) S n j+ 1 2 -S n j-1 2
.

Well-balanced fluxes

The convective part (C) is discretized with a first order explicit finite volume scheme. The numerical flux for the transport terms is constructed so that it preserves the steady states and this will provide the value of h j+ 1 2 and Θ j+ 1 2

. Following [START_REF] Tiam Kapen | A New Flux Splitting Scheme Based on Toro-Vazquez and HLL Schemes for the Euler Equations[END_REF][START_REF] Toro | Flux splitting schemes for the Euler equations[END_REF], the idea consists in splitting the flux and then using a specific numerical fluxes for each part. Here, we propose to decompose the flux term in three parts:

∂ x (hv), ∂ x (huv + F 2 Π) ∂ x (hΘv) v∂ x Π + h m λ 2 ∂ x v v∂ x v + F 2 h m ∂ x Π = ∂ x hv huv + F 2 Π hΘv 0 0 + 0 0 0 v∂ x Π v∂ x v + ∂ x 0 0 0 h m λ c v 1 h m F 2 Π = ∂ x F c + F * nc + ∂ x F a
For the non-conservative part F * nc , we use a non-conservative upwind scheme. For the linear acoustic part F a , we propose a modified acoustic flux that takes into account spatial variation due to the source term using the Jin-Levermore method [START_REF] Jin | Numerical Schemes for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF].

The linear acoustic part writes:

∂ t Π + h m λ 2 ∂ x v = 0, ∂ t v + F 2 1 h m ∂ x Π = 0,
and can be diagonalized as follows:

     ∂ t h m λ F v -Π -F λ∂ x h m λ F v -Π = 0, ∂ t h m λ F v + Π + F λ∂ x h m λ F v + Π = 0.
To define the intermediate value of Π and v in the fluxes at node x j+ 1 2 , we thus consider the upwinded quantities:

       h m λv F -Π j+ 1 2 = h m λ F v(x + j+ 1 2 ) -Π(x + j+ 1 2
),

h m λv F + Π j+ 1 2 = h m λ F v(x - j+ 1 2 ) + Π(x - j+ 1 2 ). (6)
Following [START_REF] Jin | Numerical Schemes for Hyperbolic Conservation Laws with Stiff Relaxation Terms[END_REF] and as already done for the Euler gravity system in [START_REF] Franck | Finite Volume Scheme with Local High Order Discretization of the Hydrostatic Equilibrium for the Euler Equations with External Forces[END_REF], we precise these formula by considering the possible spatial variation of Π at equilibria:

Π(x j ) ≈ Π(x - j+ 1 2 ) - ∆x 2 ∂ x Π(x j+ 1 2 ) ≈ Π(x - j+ 1 2 ) + ∆x 2 h(x j+ 1 2 )Θ(x j+ 1 2 )g∂ x z(x j+ 1 2
).

We define similarly Π(

). Then plugging these values in ( 6) and considering

v(x + j+ 1 2 ) = v(x j+1 ) and v(x - j+ 1 2 ) = v(x j )
, we obtain at the discrete level:

       h m λv F -Π j+ 1 2 = h m λ F v j+1 -Π j+1 + ∆x 2 h j+ 1 2 Θ j+ 1 2 g z j+1 -z j ∆x , h m λv F + Π j+ 1 2 = h m λ F v j + Π j -∆x 2 h j+ 1 2 Θ j+ 1 2 g z j+1 -z j ∆x . (7)
Then we obtain the following intermediate values for the modified acoustic fluxes:

       v * j+ 1 2 = 1 2 v j+1 + v j -F 2h m λ c Π j+1 -Π j + gh j+ 1 2 Θ j+ 1 2 (z j+1 -z j ) , Π * j+ 1 2 = 1 2 Π j+1 + Π j -h m λ c 2F v j+1 -v j .
Finally, for the convection part, we use an upwind scheme at the velocity v * for the pressure term. To ensure the well-balanced property, we consider the following discretization for the source terms [START_REF] Franck | Finite Volume Scheme with Local High Order Discretization of the Hydrostatic Equilibrium for the Euler Equations with External Forces[END_REF]:

S j = -g 1 2 S j+ 1 2 + S j-1 2 . ( 8
)
where S j+ 1 2 are defined in (5).

Proposition 2 Considering that u j = 0, v j = 0 and Π j = 1 2 gΘ j h 2 j , the schemes for the convective part (C) and for the wave part (W) are well-balanced for the three type of steady states:

Θ j = cst, h j + z j = cst, z j = cst, Θ j (h j ) 2 /2 = cst, h j = cst, z j + h j ln(Θ j ) = cst, ( 9
)
if we choose in the scheme

h j+ 1 2 = 1 2 (h j + h j+1 ), Θ j+ 1 2 =      Θ j+1 -Θ j ln(Θ j+1 )-ln(Θ j ) , if Θ j+1 Θ j , Θ j , if Θ j+1 = Θ j .
Indeed, with this choice, then v * )/∆x = (Π n j+1 -Π n j-1 )/(2∆x) and the soure term. This follows from the choice [START_REF] Berthon | Well-balanced scheme to capture non explicit steady-states: Ripa model Mathematics of Computation[END_REF]. We do not detail the computation here. For the implicit scheme for (W ), the preservation of these steady states results from this balance between discret gradient pressure and the source term and also from the centered discretization.

Numerical results

First, we consider a classical test case for well-balanced schemes. The initial data is taken as a steady state and we compare the classical Rusonov scheme with the semi-implicit two-speed well-balanced scheme (SI two-speed WB) on a large time interval [0, T f ]. The CFL condition of our scheme is given by ∆t ≤ ∆x max x |u(t, x) + F (t) h(t, x)Θ(t, x)g| .

The parameter F (t) depends on an important parameter F min . If u = 0, there is still a small part of the gravity wave dynamics which is explicit and generates a CFL condition: ∆t ≤ ∆x/(F min hΘg). Hence, decreasing this parameter F min increases the time step. We remark also that when F min = 1, the scheme is explicit. The initial data are given by u 0 (x) = v 0 (x) = 0, Π 0 (x) = 1 2 Θ 0 (x)h 0 (x) 2 with the following topography z(x), initial water height h 0 (x) and temperature Θ 0 (x):

(ST1) z(x) = 0.1 + G x 0 ,σ (x), h 0 (x) = 8.0 -z(x), Θ 0 (x) = 1, (ST2) z(x) = 1, h 0 (x) = 1.0 + 0.2G x 0 ,σ (x), Θ 0 (x) = 1 gh 0 (x) 2 , (ST3) z(x) = x(1 -x), h 0 (x) = 1, Θ 0 (x) = 2e -x(1-x) .

Conclusion

In this paper, we propose a semi-implicit relaxation scheme for the Ripa model, that is well-adapted to treat the low-Froude regime. Indeed, we are able to take very large time step compared to the gravity waves time scale and the accuracy of the scheme is independent of the Froude number. The two-speed relaxation allows to have an implicit step with a constant linear Laplacian. Additionally the scheme is able to preserve non-trivial steady states with a very good accuracy, with small errors for very large time steps. However the scheme is not able to treat wet/dry transitions. The 2D extension have been performed for the Euler system in [5]. In the future, we propose to extend the method to flows around equilibria MHD system.

j+ 1 2 for

 2 the quantities h, hu, hΘ and we use Π * j+ 1 2

j+ 1 2 = 2 in

 22 0 for all the steady-states. Consequently, all the transport terms at velocity v * j+ 1 (C) vanish. It remains to prove that the balance between the pressure term (Π *

Fig. 1

 1 Fig.1Left: explicit Rusanov scheme; In green the initial data. In red the solution on a semi-coarse grid (1200 cells), in blue the solution on a fine grid (12000 cells). Right: SI two-speed WB; in green the initial data. In red the solution on a coarse grid (600 cells), in blue the solution on a semi-coarse grid (4800 cells).

) and σ = 0.06. We consider g = 1, N c = 200 cells and T f = 20. The results are given in the table below.

The results show that the SI two-speed WB is exactly well-balanced in the explicit version (F min = 1). In the semi-implicit case, it is more complicated. At the theoretical level the scheme is well-balanced, but the implicit part generates very small errors close to machine precision and these errors are slowly propagated. The scheme does not preserve the steady states exactly. However, we remark that with F min = 0.005, the errors are very small (between 1.0E -11 and 1.0E -13 ) with a time step 200 larger than the one of the explicit schemes. In the following, we will show that the scheme well capture the flow around steady states and that these pertubative errors does not deteriorate the results contrary to classical non well-balanced schemes. Then we introduce a perturbation in the (ST3) test-case, h 0 (x) = 1+0.001G 0.8,σ (x), on the domain is [0, 3]. In Figure [START_REF] Tiam Kapen | A New Flux Splitting Scheme Based on Toro-Vazquez and HLL Schemes for the Euler Equations[END_REF], we compare the explicit Rusanov scheme with the SI two-speed WB scheme. We observe that the non-WB explicit Rusanov scheme is not able to capture the physical perturbation: the scheme creates a numerical perturbation larger that the physical one (in red), for 1200 cells, or with the same size (in blue) with 12000 cells. The SI two-speed WB scheme does not create numerical perturbations and capture correctly the physical perturbation with a coarser grid. The results are better with 600 cells and the SI two-speed WB scheme that with 12000 cells with the non-WB Rusanov scheme. Additionally, the results obtained with the SI two-speed WB scheme are given with a 10 times larger time step than the explicit one. We cannot increase too much the time step here since the implicit part of the schemes would create numerical diffusion on the gravity wave part.