R D H U I N

Noise generation in the solid Earth, oceans, and atmosphere, from non-linear interacting surface gravity waves in finite depth

Keywords: Hydrodynamic noise, Surface gravity waves, Air/sea interactions 1234546789 28 8A26BC8CD6 EF

Oceanic pressure measurements, even in very deep water, and atmospheric pressure or seismic records, from anywhere on Earth, contain noise with dominant periods between 3 and 10 seconds, that is believed to be excited by ocean surface gravity waves. Most of this noise is explained by a nonlinear wave-wave interaction mechanism, and takes the form of surface gravity waves, acoustic or seismic waves. Previous theoretical works on seismic noise focused on surface (Rayleigh) waves, and did not consider finite depth effects on the generating wave kinematics. These finite depth effects are introduced here, which requires the consideration of the direct wave-induced pressure at the ocean bottom, a contribution previously overlooked in the context of seismic noise. That contribution can lead to a considerable reduction of the seismic noise source, which is particularly relevant for noise periods larger than 10 s. The theory is applied to acoustic waves in the atmosphere, extending previous theories that were limited to vertical propagation only. Finally, the noise generation theory is also extended beyond the domain of Rayleigh waves, giving the first quantitative expression for sources of seismic body waves. In the limit of slow phase speeds in the ocean wave forcing, the known and well-verified gravity wave result is obtained, which was previously derived for an incompressible ocean. The noise source of acoustic, acoustic-gravity and seismic modes are given by a mode-specific amplification of the same wave-induced pressure field near the zero wavenumber.

Introduction

Ocean waves generate noise in a wide range of acoustic frequencies f s . The upper end of the spectrum, f s > 100 Hz, is dominated by wave breaking and associated bubbles [START_REF] Knudsen | Underwater ambient noise[END_REF], whereas the lower frequency part, nominally f s < 2 Hz is mostly expected to be caused by the nonlinearity of the hydrodynamic equations, on which we focus here. The general sound generation by fluid flows was described theoretically by [START_REF] Lighthill | On sound generated aerodynamically[END_REF]. [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF] showed how seismic waves can be generated by the same process, with noise radiating along the Earth's crust in the form of Rayleigh waves. That theory was extended to random waves by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF], and later cast in the more general framework of wave-wave interactions [START_REF] Hasselmann | Feynman diagrams and interaction rules of wave-wave scattering processes[END_REF]. Work on compressible flows has also been extended to the study of tsunamis. In that context, [START_REF] Okal | Seismic parameters controlling far-field tsunami amplitudes: A review[END_REF] discussed the compressibility effect on gravity modes, which we will call here 'acoustic-gravity modes', and the gravity effect in seismic 'pseudo-Rayleigh' waves, that we will refer to as 'Rayleigh' modes. Interest in seismic noise has risen sharply over the last few years with the enforcement of the Comprehensive Nuclear Test Ban Treaty, and motivated by the work of [START_REF] Shapiro | High-resolution surface-wave tomography from ambient seismic noise[END_REF] who demonstrated that seismic noise correlation could provide a unique monitoring method for the properties of the solid Earth.

Recent numerical models based on the Longuet-Higgins-Hasselmann (LH-H) theory for Rayleigh wave generation have shown good agreement of modeled seismic noise spectra with observations [START_REF] Kedar | The origin of deep ocean microseisms in the north Atlantic ocean[END_REF][START_REF] Ardhuin | Ocean wave sources of seismic noise[END_REF]. It is still unclear whether most of the uncertainties on the modeled noise level can be attributed to errors in the seismic sources, associated with a poorly constrained directional distribution of the ocean surface wave spectrum, or to errors in the seismic propagation. Seismic observations have also revealed body waves (e.g. [START_REF] Koper | On the composition of Earth's short-period seismic noise field[END_REF][START_REF] Landès | Origin of deep ocean microseisms by using teleseismic body waves[END_REF][START_REF] Hillers | Global oceanic microseism sources as seen by seismic arrays and predicted by wave action models[END_REF], for which no complete theory has been proposed to date. [START_REF] Vinnik | Sources of microseismic P waves[END_REF] did propose a theory for compressional (P ) waves, but he did not consider the important effect of the water layer.

Further, the generation of Love waves, which are surface shear waves polarized in the horizontal direction, is not well understood. These Love waves are particularly important for frequencies below 0.02 Hz or above 1 Hz (e.g. [START_REF] Bonnefoy-Claudet | The nature of noise wavefield and its applications for site effects studies a literature review[END_REF][START_REF] Nishida | Background Love and Rayleigh waves simultaneously generated at the Pacific ocean floors[END_REF][START_REF] Kurrle | The horizontal hum of the Earth: A global background of spheroidal and toroidal modes[END_REF]. The low frequency Love waves may be excited by the direct action of long surface gravity waves (known as infragravity waves), on a sloping bottom [START_REF] Fukao | Seafloor topography, ocean infragravity waves, and background Love and Rayleigh waves[END_REF], in a way similar to the generation of primary microseisms described by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF]. Here we narrow the scope of our investigation and only consider non-linear wave-wave interactions, which leads to noise with frequencies double that of the surface gravity waves.

Finally, the level of acoustic noise has also been explored as a potential source of information on the poorly known directional spectrum of short gravity waves [START_REF] Tyler | Wave directional spectra from synthetic aperture observations of radio scatter[END_REF]. [START_REF] Farrell | What do deep sea pressure fluctuations tell about short surface waves?[END_REF], 2010) and [START_REF] Duennebier | Wind, waves, and acoustic background levels at station aloha[END_REF] showed a large variability of the spectral level in the frequency range 0.1-50 Hz that is clearly related to the sea state. Their interpretation of the data, following previous studies of underwater noise (e.g. [START_REF] Hughes | Estimates of underwater sound (and infrasound) produced by nonlinearly interacting ocean waves[END_REF][START_REF] Lloyd | Underwater sound from surface waves according to the Lighthill-Ribner theory[END_REF][START_REF] Kibblewhite | Wave-wave interactions, microseisms, and infrasonic ambient noise in the ocean[END_REF], is based on sound generated by waves in an unbounded ocean. Although this simplified approach is reasonable for high frequency noise, the neglected reflection from the seafloor and subsequent reverberations in the bounded ocean may strongly amplify the lower frequency resonant modes.

Given the renewal of interest in seismic and acoustic noise, we found it appropriate to revisit Hasselmann's theory. We thus illustrate, correct, and add a few missing aspects. These corrections include important terms for intermediate and shallow water that have not been considered before for the compressible conditions, although they were verified in the incompressible limit [START_REF] Herbers | Nonlinear wave interactions and high-frequency seafloor pressure[END_REF]. The compressible equations of motion are used to derive a consistent solution in terms of gravity, acoustic and seismic modes, including both surface Rayleigh waves and compressional (P) or shear (S) body waves. This solution for body waves has not been given before.

Building on the seismo-acoustic paradigm proposed by [START_REF] Arrowsmith | The seismoacoustic wavefield: a new paradigm in studying geophysical phenomena[END_REF], the basic idea of the present paper is that all modes of motion can be excited by ocean surface gravity waves (OSGW) of any frequency. For a given pair of interacting frequencies f and f ′ , the frequency of the generated noise is f +f ′ , and the different types of waves are only distinguished by their phase velocity, or equivalently by their horizontal wavenumber K which is the norm of the vector sum of the wave numbers, K = k + k ′ of the interacting OSWG, as shown in figure 1. This type of wave-wave interaction is one of the lowest order interactions [START_REF] Hasselmann | Feynman diagrams and interaction rules of wave-wave scattering processes[END_REF].

In the physical space, one such interaction excites waves with horizontal wavelength L h = 2π/K, and different vertical patterns in the atmosphere, ocean and crust, due to the very different speeds for compression waves in these three media. Figure 2 shows how the same forcing can give almost vertical-propagating waves in the atmosphere, waves that propagate almost horizontally in the ocean, and evanescent waves in the crust. The analysis of a pair of interacting wave-trains is key to our interpretation of the different noise modes. The broad-banded wave spectrum of ocean waves results in the superposition of all possible pairs of OSGW waves trains, and thus all possible noise waves propagating or evanescent, that radiate in all directions. This broad spectrum allowed [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF] to consider all interactions at once, replacing the wave-induced forcing by an equivalent point source exerted on the sea surface. However, that latter approach is only valid for noise wavelengths much larger than those of the interacting OSWG. Table 1 points directly to the main theoretical results of the present paper, and all symbols and notations are listed in tables 2-3.

Water waves theory and noise sources

Here we give only a brief derivation of the solution, which is a straightforward compressible extension to the solution given by [START_REF] Hasselmann | On the non-linear energy transfer in a gravity wave spectrum, part 1: general theory[END_REF]. Water column motions are expanded in powers of the sea surface slope with a linear motion for which compressible effects may be neglected [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF], and a second order motion with pressure p 2 and velocity potential φ 2 . It is for this second order motion that we are Figure 2. Left, schematic of second-order pressure field, in colours, associated with double-frequency noise, forced by the interaction of a single pair of directionally opposing monochromatic wave trains on the sea surface that form periodic groups. The individual waves in the group advance very slowly, at the OSWG phase speed. However, in this particular case, the group as a whole travels slightly faster than the sound speed in water, and slower than the crust shear and compression velocity. As a result, the crust elastic waves are evanescent while ocean and atmospheric waves propagate also in the vertical. In the ocean, the superposition of upward and downward waves (dashed lines) gives the vertical mode structure. The number of waves in the group was reduced for visibility, and for the same reason the amplitudes of the sea surface and bottom elevations have been exaggerated as well as the pressure fluctuation in the atmosphere relative to those in the ocean. Right, schematic of ocean waves with a relatively broad spectrum, giving rise to the interaction of all possible pairs of wave trains and noise radiation in all directions. extending previous results. Following again Longuet-Higgins (1950) we will see that φ 2 approximately obeys the same linear acoustic wave equation as φ 1 , so that the forcing of the double-frequency noise is only due to boundary conditions. In contrast to [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF], we consider here the general finite depth expression for the surface pressure forcing p 2,surf , and in the boundary condition at the ocean bottom the additional forcing p 2,bot that accounts for the Bernoulli effect of the near-bed orbital wave motion. The theory presented here extends the second order finite depth theory of [START_REF] Hasselmann | On the non-linear energy transfer in a gravity wave spectrum, part 1: general theory[END_REF] to a compressible ocean and elastic seafloor.

Equations of motion

We decompose the water density into a mean value ρ w and a perturbation ρ ≪ ρ w . Neglecting stratification effects due to temperature and salinity, the fluctuations in pressure p and water density ρ are related by an equation of state, which involves the speed of sound in water α w (Lighthill 1978, eq. 32).

dp dt = α 2 w dρ dt , (2.1)
We assume that the motion is irrotational so that the velocity field is given by the gradient of the velocity potential φ. This assumption, which implies a frictionless interior ocean, is well supported by the observed weak attenuation of swells propagating across ocean basins [START_REF] Ardhuin | Observation of swell dissipation across oceans[END_REF]) and local comparisons of the observed wave orbital motion with second order wave theory (e.g. [START_REF] Herbers | Field observations of orbital velocities and pressure in weakly nonlinear surface gravity waves[END_REF]. Vorticity effects can be important in the vicinity of the surface and bottom boundary layers (e.g. [START_REF] Longuet-Higgins | Mass transport in the boundary layer at a free oscillating surface[END_REF], or in the presence of sheared currents (e.g. [START_REF] Peregrine | Interaction of water waves and currents[END_REF]). The former effect has little influence on the pressure field, which is our primary interest, and we will not consider here the effects of currents. Because of the importance of the apparent gravity acceleration g, which defines the vertical axis, we separate horizontal and vertical components using vectors and gradient operators in the horizontal plane, e.g. u = ∇φ = (∂φ/∂x, ∂φ/∂y) and w = ∂φ/∂z. The conservation of mass of sea water is

dρ dt = -ρ w ∇ 2 φ -ρ w ∂ 2 φ ∂z 2 . (2.2) Equations (2.2)-(2.1) can be combined to eliminate ρ, dp dt = -ρ w α 2 w ∇ 2 + ∂ 2 ∂z 2 φ (2.3)
From eq. (2.1), the water density is only a function of the pressure. The two unknowns p and φ are also related by the momentum conservation equation, with can be cast in the form of Bernoulli's equation (see e.g. [START_REF] Lamb | Hydrodynamics[END_REF]

ζ = ζ 1 + ζ 2 , φ = φ 1 + φ 2 .
(2.8a,b)

Linear solution

We consider the case of a constant depth h. Compressibility effects in the linear solution are negligible for our purpose (see Longuet-Higgins 1950, eq. 123), so that we may use

ζ 1 = k,s Z s 1,k e i(k•x-sσt) , φ 1 ≃ k,s -i sg σ cosh(kz + kh) cosh(kh) Z s 1,k e i(k•x-sσt) (2.9a,b)
where k is the norm of the horizontal wavenumber vector k, s is a sign index equal to -1 or 1, so that s = 1 corresponds to waves propagating in the direction of the vector k, and s = -1 corresponds to the opposite direction. The radian frequency σ is given by the dispersion relation for linear waves (de Laplace 1776) , σ = gk tanh(kh), (2.10)

giving the group speed

C g = ∂σ ∂k = σ k 1 2 + kh sinh(2kh)
.

(2.11)

The negligible near-surface forcing

Eliminating p between (2.3) and (2.4), we obtain the acoustic wave equation

α 2 w ∇ 2 + ∂ 2 ∂z 2 φ = d dt ∂φ ∂t + 1 2 |∇φ| 2 + ∂φ ∂z 2 , (2.12)
where the gravity term has been removed by our approximation of a constant mean density. That term is usually neglected eventually in the solution [START_REF] Stoneley | The effect of the ocean on Rayleigh waves[END_REF]. The expression for the C(t) term in eq. (2.4), is given by eq. ( 28) in [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF] and it is actually zero provided that there are no standing waves, i.e. we do not consider the case of waves trains of exactly equal frequency and opposite direction. For broad wave spectra we may indeed neglect these contributions since the measure of such pairs of wave components is zero while the measure of nearly opposite waves is finite. Our eq. (2.12) corresponds to eq. ( 130) in [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF], who showed that the non-linear terms yield a contribution (his F terms) that is of the order of g/(2α 2 w k) times the other terms. Even for very long waves with a wavelength of 50 km, this factor is 0.02, so that we shall neglect this source of wave forcing. This leads to a linear wave equation for the second order velocity potential

∂ 2 φ 2 ∂ 2 t ≃ α 2 w ∇ 2 + ∂ 2 ∂z 2 φ 2 , for -h z 0 (2.13)

Surface forcing

Because the compressibility only affects the mass conservation equation, it does not modify the kinematic and dynamic boundary conditions. These are given by [START_REF] Hasselmann | On the non-linear energy transfer in a gravity wave spectrum, part 1: general theory[END_REF]. The unknown ζ is eliminated from the linear terms by adding ∂(2.4)/∂t evaluated with (2.5), and g× (2.7). This combination of kinetic and dynamic boundary condition gives an equation for the velocity potential to second order in the wave slope, valid at z = 0. Keeping the lowest order non-linear terms we have,

∂ 2 ∂t 2 + g ∂ ∂z φ 2 = -g ∇φ 1 • ∇ζ 1 + ζ 1 ∂ 2 φ 1 ∂ 2 z - 1 2 ∂ ∂t |∇φ 1 | 2 + ∂φ 1 ∂z 2 + 2ζ 1 ∂ 2 φ 1 ∂t∂z .
(2.14)

We now give the explicit form of the right hand side using eqs. (2.9). As noted by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF], we may rewrite (2.14) as

∂ 2 ∂t 2 + g ∂ ∂z φ 2 = - 1 ρ w ∂ p 2,surf ∂t .
(2.15)

Namely, our problem is equivalent to the effect of a pressure field p 2,surf applied at z = 0. Using the linear solution (2.9),

p 2,surf = ρ w k,s,k ′ ,s ′ D z (k, s, k ′ , s ′ ) Z s 1,k Z s ′ 1,k ′ e iΘ(k,k ′ ,s,s ′ )
(2.16)

with the phase function of interacting wave trains defined by

Θ(k, k ′ , s, s ′ ) = [(k + k ′ ) • x -(sσ + s ′ σ ′ ) t]
(2.17)

and coupling coefficient D z (k, s, k ′ , s ′ ) given by eq. (A 1). The shape of this pressure pattern is well understood when the full sum in eq. (2.16) is reduced to only two interacting deep water wave trains with amplitudes a and a ′ and slightly different frequencies σ and σ ′ , traveling in opposite directions, as shown in figure 3. For kh ≫ 1 we may use D z (k, 1, -k, 1) = -2σ 2 and D z (k, 1, k, 1) = 0. Defining K = k -k ′ , the equivalent surface pressure is given by the sum interaction,

p 2,surf = -2ρ w σ 2 aa ′ {cos [Kx + 2σt]} .
(2.18)

The corresponding difference interaction yields a short wave, with wavenumber 2k, that will not be considered here (see e.g. [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF]). We note that the wave forcing p 2,surf is out of phase with -u 2 /2, the part of the Bernoulli head which comes from the horizontal velocities. Indeed, as the wave trains propagate in opposite directions their velocity partially cancel where the surface elevations add up, thus causing a pattern of higher pressure under the groups of high waves and lower pressure under the lulls, that is opposite to the familiar Bernoulli effect of set-down under groups of uni-directional waves [START_REF] Longuet-Higgins | Radiation stresses and mass transport in surface gravity waves with application to 'surf beats[END_REF]. As a result the double-frequency perturbations of -u 2 /2 have a sign opposite to p 2,surf . That effect will be very important for waves in shallow water. In deep water, the contribution of -w 2 /2 must also be considered. Using σ ≃ σ ′ , we have

u = σa cos(kx -σt) -σ ′ a ′ cos(kx + σt) (2.19) w = σa sin(kx -σt) -σ ′ a ′ sin(kx + σt) (2.20) -(u 2 + w 2 )/2 = aa ′ σ 2 cos [Kx + 2σt] -σ 2 a 2 + a ′2 /2.
(2.21)

Eq. (2.18) generalizes the result given by Longuet-Higgins (1950) for equal wave periods. The standing wave studied by [START_REF] Longuet-Higgins | A theory of the origin of microseisms[END_REF] is thus obtained, somehow paradoxically, as the limit of wave groups that travel at an infinite speed, but that are infinitely long.

For wave directions nearly opposite, instead of exactly opposite, this first term can propagate in any horizontal direction, given by the direction of K = k + k ′ .

Noise spectrum and finite depth effects

Both LH50 and H63 used the deep-water approximation

p 2,surf ≃ ρ w |∇φ| 2 + ∂φ ∂z 2 valid for kh ≫ 1, (2.22) 12222222222222223222222222222222422222222222222252222222222222226222222222222227122222222222273 123456 -2 -2 897AB2C -2 0 2 891 D292732C DE29277A662C D29252CF2229231112C 3 3 0 2 89A62C Figure 3
. Generation of supersonic wave groups (in blue) by the superposition of two opposing deep water monochromatic wave trains of nearly equal periods T and T ′ . The curves show the surface elevation of the individual wave trains (black and red) or their combination (blue). The group with length 12 km propagates in the same direction as the wave train with the shortest wavelength. If the wave trains are not exactly directionally opposing, the group propagates in the direction of the vector K = k + k ′ . The red and black dots are attached to the wave crest of each train, and move 100 times slower than the blue group.

instead of the more complex but more general form given by eq. (2.16). We will illustrate the important differences between these two expressions for kh < 1, by considering the forcing of very long components. Indeed, sound waves in the ocean have velocities in excess of 1.4 km s -1 . As a result, regardless of the vertical wavenumber l, the horizontal wavenumber vector is relatively small K = (ω 2 s /α 2 w -l 2 ) < ω s /α w . Thus the acoustically noisy wave interactions verify K ≪ k, which gives k ′ ≃ -k, and f ≃ f ′ with f s ≃ 2f . We may thus focus on the estimation of the spectrum of p 2,surf at K ≃ 0.

For this purpose we introduce the spectral density of this surface pressure in the three spectral dimensions (K x , K y , f s ). Using the Fourier amplitude p 2,surf (K, f s ) of the forcing pressure p 2,surf , with wavenumber vector K and frequency f s ,

F p2,surf (K, f s ) = 2 lim |dK|→0,dfs→0 | p 2,surf (K, f s )| 2 dK x dK y df s , (2.23)
The factor 2 in the expression makes this a single-sided spectrum, with non-zero values only for f s 0. This spectral density is in a three-dimensional spectral space, with S.I. units of Pa 2 m 2 s; it is denoted F p3D in [START_REF] Ardhuin | Ocean wave sources of seismic noise[END_REF]. Using (2.16) the surface pressure spectrum can be expressed in terms of quadratic products of the (linear) sea surface elevation spectrum

E(k x , k y ) = 2 lim |dk|→0 Z + 1,k 2 dk x dk y (2.24)
with a coupling coefficient from (A 1) that simplifies for K ≃ 0 to

D z (k, 1, -k, 1) = -2σ 2 1 + 1 4 sinh 2 (kh)
.

(2.25)

To transform the spectra to frequency-direction space we use the Jacobian transfor-mation

E(f, θ) = 2πk C g E(k x , k y ).
(2.26)

We now introduce the directional distribution M such that E(f, θ) = E(f )M (f, θ), and we define the directional integral

I(f ) = π 0 M (f, θ)M (f, θ + π)dθ. (2.27)
With these notations we finally obtain

F p2,surf (K ≃ 0, f s ) ≃ ρ 2 w D 2 z E(f, θ)E(f, θ + π) C 2 g dk x dk y k 2 4π 2 df s ≃ ρ 2 w g 2 f s E 2 (f )I(f ) tanh 2 (kh) 1 + 2kh sinh(2kh) 1 + 1 4 sinh 2 (kh) 2 (2.28)
In deep water (kh ≫ 1), the equivalent surface pressure is p 2,surf ≃ ρ w (u 2 1 + w 2 1 ), and its spectrum [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF][START_REF] Ardhuin | Ocean wave sources of seismic noise[END_REF])

F p2,surf (K ≃ 0, f s ) = ρ 2 w g 2 f s E 2 (f )I(f ). (2.29)
As shown in figure 4, this spectral density of surface pressure is, in finite depth, up to four times smaller than the deep water approximation obtained from eq. ( 2.22) and used by [START_REF] Webb | The earth's 'hum' is driven by ocean waves over the continental shelves[END_REF] and [START_REF] Tanimoto | Excitation of normal modes by non-linear interaction of ocean waves[END_REF]. This is because the surface pressure actually combines two terms. One is ρ w (u 2 1 + w 2 1 )/2, from the momentum conservation equation -which here takes the form of the Bernoulli equation (2.4) -and the other, equal in magnitude in deep water, comes from the non-linearity of the surface boundary condition (2.7). As kh goes to zero, both the latter term and the vertical velocity contribution to the Bernoulli pressure are small compared with ρ w u 2 1 /2, thus reducing the surface pressure forcing by a factor four.

This effect is not important for the dominant microseismic peak, generated by waves of period about 10 s, but it may be important for the much longer waves, known as hum, driven by long surface gravity waves. However, in that case, one should also consider the direct action of the wave-induced pressure on the bottom.

Additional bottom forcing

To evaluate noise in the water column, both the second order pressure forcing applied at the surface and bottom have to be taken into account. This includes the well known Bernoulli effect of a pressure drop in response to an increase in velocity. At the bottom, this wave-induced pressure is -ρ w u 2 1 /2 which is exactly out of phase with p 2,surf . In the limit kh → 0, it cancels the source of noise that would have resulted from the surface forcing alone. In physical terms, for kh = 0, we have the same momentum balance at the sea surface and bottom, and thus the same pressure perturbations, which are zero because the pressure just below the surface has to match the atmospheric pressure.

For any value of kh, the coupling coefficient given by eq. (2.25) differs from the full second order coefficient for the bottom pressure (e.g. eq. 4 in [START_REF] Herbers | Wind-wave nonlinearity observed at the sea floor. part I: forced-wave energy[END_REF], which also involves the Bernoulli head (the bracket in eq. 2.4). However, that extra term is also relevant to the generation of seismic noise due to the bottom boundary condition that couples the solid crust to the water column. Indeed, the second-order pressure perturbation at the bottom writes,

p 2 (-h) = -ρ w ∂φ 2 ∂t + p 2,bot , (2.30)
where the Bernoulli head contribution to the pressure can be expressed from the first order wave amplitudes,

p 2,bot = ρ w k,s,k ′ ,s ′ D pb (k, s, k ′ , s ′ , z = -h) Z s 1,k Z s ′ 1,k ′ e iΘ(k,k ′ ,s,s ′ ) , (2.31)
with a coupling coefficient D pb given by eq. (A 2). We may interpret the bottom pressure (2.30) as the sum of the surface forcing p 2,surf transmitted to the bottom by φ 2 , and a direct effect of the Bernoulli head at the bottom which is an additional forcing p 2,bot that partly cancels p 2,surf .

We shall see in the next section that the forcing term for seismic noise is p 2,surf + cos(lh) p 2,bot , with l K ≪ k the vertical wavenumber in the water. For shallow water gravity waves, kh ≪ 1 and thus cos(lh) ≃ 1 so that the effective forcing term becomes p 2,surf + p 2,bot , which equals the bottom pressure in the incompressible limit. The shallow water asymptote of the spectrum of this total forcing term is very different from the surface pressure only. Compared to eq. (2.28), the 1 + 0.25/ sinh 2 (kh)

2 factor is now replaced by 1. For kh ≪ 1, this is a factor (kh) 4 /16 smaller, as shown in figure 4. This asymptote is relevant for the hum, the noise with periods larger than 30 s, which is believed to be driven by long (infragravity) surface gravity waves [START_REF] Webb | The earth's 'hum' is driven by ocean waves over the continental shelves[END_REF]. The source of this hum is attenuated by several orders of magnitude on the continental shelves and not amplified according to the deep water approximation given by eq. (2.22).

3. From surface pressure spectrum to noise spectra

Dispersion relation and modes

The problem of noise generation by waves has been reduced to that of noise generation by an equivalent surface pressure field p 2,surf (x, y, t), with the possible addition of a bottom pressure field p 2,bot (x, y, t) for finite depths. From a statistical point of view, this equivalent surface pressure is fully represented by its spectrum, F p2,surf (K, f s ) where K = k+k ′ is the sum of the interacting wavenumbers and f s = f ±f ′ is the sum or difference of the interacting frequencies. Subjected to this surface forcing, our linearized wave equation (2.13) will have linear solutions. In particular, any propagating or evanescent solution will take the following form

φ 2 ∝ exp[i(K x x + K y y + lz -ω s t)],
(3.1) which, substituted in the wave equation, gives the dispersion relation

-ω 2 s + α 2 w K 2 + l 2 = 0, (3.2)
where ω s = 2πf s . Both ω s and K = |K| are imposed by the forcing, so that the magnitude of the complex vertical wavenumber l is given by . Compared to deep water, the source of seismic noise power in finite water depth is amplified by a factor tanh 2 (kh) [1 + 2kh/sinh(2kh)], when accounting for both the surface and bottom forcing. This is very different from the approximation that considers the surface forcing only (leading to eq. 2.28), even more so when using its deep water approximation (eq. 2.22) which gives an additional factor 4 difference for kh ≪ 1.

l = K ω 2 s K 2 α 2 w -1, (3.3) which yields φ 2 = Ce il(z+h) + De -il(z+h) e iΘ , Θ = (K x x + K y y -ω s t) (3.
where C and D are the bottom amplitudes of the upward and downward propagating waves, determined by the surface and bottom boundary conditions. There are two classes of solutions. Those for which ω s /K α w and thus l is imaginary: these are 'acoustic-gravity ' modes with an amplitude that decays exponentially from the surface. For the shortest components, we have ω s /K ≪ α w and thus l = ik, corresponding to the incompressible limit in which only gravity is important. The other class of solutions, for which ω s /K α w and l is real, are the acoustic modes that propagate along both vertical and horizontal dimensions, as illustrated in figure 5.

The equations of elasticity in the crust yield the same wave equation, one for compressional and another for shear motions, with α w replaced by the crusts compressional wave speed α c and shear wave speed β s , respectively (e.g. [START_REF] Aki | Quantitative seismology[END_REF]. Since we have α w < β s < α c , we can distinguish four different regimes depending on the value of the horizontal phase speed ω s /K, relative to these three velocities (fig. 1).

In the family of gravity-acoustic modes, there is no resonance in the forcing, namely there are no free waves with K = k + k ′ and f s = f ± f ′ and thus the motion is locked with the forcing wave groups, as verified by [START_REF] Herbers | Nonlinear wave interactions and high-frequency seafloor pressure[END_REF]. The acoustic modes propagate obliquely to the vertical away from the surface, and interact with the ocean bottom (figure 5). The seismic Rayleigh waves modes propagate horizontally, combining an acoustic-like motion in the water and evanescent elastic waves in the crust, illustrated . Schematic of wave groups and forced acoustic and seismic wave motion. For a given wave group period T , the horizontal wavelength L = 2π/K can be larger than the acoustic wavelength in the water due to the oblique incidence of the sound waves. The superposition of two obliquely propagating sound waves (arrows) forms a mode pattern that propagates horizontally at a supersonic speed. The acoustic wavelength is La = 2π/ka = 2π/ √ K 2 + l 2 = αw/T . Both vertical and horizontal wavelengths are larger. For readability, the wave and bottom amplitude is not to scale, and we have reduced the number of waves in the group from 104 to 10. Other than that, the angles are preserved. The configuration shown here corresponds to the conditions for maximum amplification of mode 1 (see below), with a vertical wavelength to water depth ratio of 0.75. For a water depth h = 4400 m depth, this corresponds to L = 7.7 km and a seismic frequency fs = 0.29 Hz.

by figure 2. The coupling with the bottom motion selects a few resonant modes that dominate the solution, with an energy that grows linearly with the propagation distance [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF]. For wavenumbers that allow a vertical propagation in the crust, we obtain body waves. These waves can be compression or shear (P or S) waves, and these two types occur in overlapping ranges of K.

Acoustic noise in an unbounded ocean

In order to illustrate the different types of solutions, it is interesting to evaluate the solution for an unbounded ocean, in which sound waves are radiated from the surface only. The velocity field and associated pressure fluctuations are

φ 2 = 1 ρ w iω s p 2 (K, f s ) -ω 2 s + igl e i[-lz+Θ(k,k ′ ,s,s ′ )] dKdf s (3.5) p 2 = p 2 (K, f s ) 1 -igl/ω 2 s e i[-lz+Θ(k,k ′ ,s,s ′ )] dKdf s (3.6)
where p 2 has been obtained using the linearized version of eq. (2.4). The measured pressure signal is the sum of the linear pressure p 1 , the second-order wave pressure p 2 given by eq. (3.6), and the Bernoulli correction p 2,B given by

p 2,B (z) = ρ w k,s,k ′ ,s ′ D pb (k, s, k ′ , s ′ , z) Z s 1,k Z s ′ 1,k ′ e iΘ(k,k ′ ,s,s ′ ) . (3.7)
We note that p 2,bot defined in eq. (2.31) is equal to p 2,B (z = -h).

We shall neglect g|l|/ω 2 s , which is bounded by the ratio between the deep water gravity and sound speeds, which is less than 0.1 for wave periods less than 180 s. We express the velocity potential as a sum of propagating (acoustic, l real) and evanescent (acousticgravity, l imaginary) modes,

φ 2 = φ 2,p + φ 2,e .
(3.8)

We get the frequency spectrum of the propagating modes by integrating over the inner regions of the wavenumber space (labelled P+S, S and R in figure 1),

F p2,p (f s ) = K<ωs/αw F p2,surf (K, f s )dK.
(3.9)

For this range of wavenumbers |k -k ′ | < K < ω s /α w , and using the relations ω s ≃ 4πf and, (for small |f -

f ′ |), |k-k ′ | ≃ 2π|f -f ′ |/C g ≃ 8π 2 f |f -f ′ |,
we obtain an upper bound for the frequency difference |f -f ′ | < g/(2πα w ) which is close to= 0.001 Hz. Typical ocean wave spectra have a relative frequency half-width σ f /f that is between 0.03 for swells and and 0.07 for wind-seas [START_REF] Hasselmann | Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project[END_REF], so that E(f ) ≃ E(f ′ ) is a good approximation for the interactions that drive long wavelength pressure fluctuations. The wave spectrum is thus broad enough for us to evaluate F p2,surf at K = 0 using eq. (2.29), and take it out of the integral in eq. (3.9). The acoustic spectrum simplifies to

F p2,p (f s ) = πω 2 s α 2 w ρ 2 w g 2 f s E 2 (f )I(f ). (3.10)
This is identical to the expression given by [START_REF] Lloyd | Underwater sound from surface waves according to the Lighthill-Ribner theory[END_REF].

Gravity noise in an unbounded ocean

The pressure associated with acoustic-gravity modes is the other part of the integral in (3.9), for K > ω s /α w . The imaginary wave number l gives a vertical attenuation of the power spectrum by a factor e -2|l|z . With that attenuation we may, for large enough depths, assume that only modes with K ≪ k contribute to the result, so that we may take F p2,surf (K, f s ) ≃ F p2,surf (K = 0, f s ), and take it out of the integrand. This approximation is valid only up to a maximum wave number K max that is a small fraction of k, K max = ǫk. For numerical applications we used ǫ = 0.2.

With this approximation we have,

F p2,e (f s , z) = F p2,surf (K = 0, f s )2π Kmax ωs/αw Ke 2|l|z dK = F p2,surf (K = 0, f s )2π Kmax 0 |l|e 2|l|z d|l| = π 2z 2 ρ 2 w g 2 f s 1 -e 2zKmax E 2 (f )I(f ) (3.11)
A previous investigation by [START_REF] Cox | Cartesian diver observations of double frequency pressure fluctuations in the upper levels of the ocean[END_REF] included an extra factor (1 + zK max ) in front of the exponential term e 2zKmax , because they neglected compressibility effects. That term , however, is negligible in the upper part of the water column, and their observations collected within 100 to 290 m of the surface in 4000 m depth, are thus not affected by this small compressibility correction.

As shown in figure 6, the oceanic pressure signals can be dominated by linear gravity waves down to depths of a few hundred meters. When looking at the double frequency band, linear waves may only dominate in the top 100 m. At these frequencies, the acoustic-gravity modes have the most important contribution between about 100 to 500 m, provided that E 2 (f )I(f ) is large enough. Propagating modes should dominate only beyond about 1000 m in the case of an unbounded ocean, or only 300 m, when accounting for the reverberation in a finite depth ocean, assuming a typical tenfold amplification for a sea floor with realistic elastic properties †. These depths will be reduced in the case of surface gravity waves with periods shorter than the 15 s swells example shown in figure 6.

Atmospheric noise source: microbaroms

The source of noise in the atmosphere can also be derived with the same formalism, as an alternative to the Green functions used by [START_REF] Waxler | The radiation of atmospheric microbaroms by ocean waves[END_REF]. Indeed, we may consider the atmospheric motion to be irrotational, so that the equations of motion are identical in the atmosphere and in an unbounded ocean, with the only difference that the atmospheric density is ρ a and the atmospheric sound speed is α a . The second-order velocity potential takes the form,

φ 2,a ∝ exp[i(K x x + K y y + l a z -ω s t)] for z > 0, (3.12) with l a = ω 2 s α 2 a -K 2 . (3.13)
Because ρ w /ρ a ≃ 1000, the air motion has only a small O(ρ w /ρ a ) local influence on the water motion, so that the solutions derived earlier for the water motion remain valid in the presence of air. The air motion, with a velocity potential φ a also obeying eq. (2.13) is fully determined from the water motion via the kinematic boundary conditions on the air and water-sides of the interface (2.7),

∂φ a ∂z - ∂φ ∂z ≃ ∇ (φ a -φ) • ∇ζ -ζ ∂ 2 (φ a -φ) ∂ 2 z at z = 0. (3.14)
From the first order potential in the air (e.g. [START_REF] Waxler | The radiation of atmospheric microbaroms by ocean waves[END_REF])

φ 1 = k,s is g σ Z s 1,k e -kz e i(k•x-sσ) (3.15)
we obtain the second order potential,

∂φ 2,a ∂z = ∂φ 2 ∂z + k,s,k ′ ,s ′ D za (k, s, k ′ , s ′ , z) Z s 1,k Z s ′ 1,k ′ e iΘ(k,k ′ ,s,s ′ ) (3.16)
and a new coupling coefficient

D za (k, s, k ′ , s ′ , z) = - 2isg σ (kk ′ + k • k ′ ) .
(3.17)

We note that for k ′ = -k, D za = 0, so that the long-wavelength motion with K ≪ k simplifies to ∂φ 2a ∂z ≃ ∂φ 2 ∂z at z = 0.

(3.18) consistent with the result given by [START_REF] Posmentier | A theory of microbaroms[END_REF] for the interaction of monochromatic wave trains, and in disagreement with a factor 8 correction proposed by [START_REF] Arendt | Acoustic radiation by ocean surface waves[END_REF]. This gives a pressure spectrum for the propagating atmospheric waves,

F p2,ap (f s ) = K<Kmax ρ 2 a l 2 ρ 2 w l 2 a F p2,surf (K, f s )dK = R(K max ) πω 2 s α 2 w ρ 2 a g 2 f s E 2 (f )I(f ). (3.19) with the non-dimensional factor R(K max ) = 2 α 2 w ω 2 s Kmax 0 l 2 l 2 a KdK == 2 α 2 w ω 2 s Kmax 0 l 2 l 2 a KdK. (3.20)
In order to avoid the singularity for l a = 0, and atmospheric ducting effects not represented here, we take K max = ω s /(2α a ) which restricts the acoustic propagation directions to be within 30 degrees from the vertical. In that case we have R(K max ) ≃ 0.54 instead of R(K max ) = 0.25 with the vertical propagation approximation of [START_REF] Waxler | The radiation of atmospheric microbaroms by ocean waves[END_REF], which replaces the l 2 /l 2 a factor in the integral by its value α 2 a /α 2 w for K = 0. Other than that, our expression is consistent with their low Mach number asymptote, i.e. σ/k ≪ α a (Waxler & Gilbert 2006, eq. 61). The present theory also allows the estimation of the evanescent wave components given by wavenumbers K > ω s /α a .

Noise in a finite depth ocean

For large depths compared to the OSGW wavelength, kh ≫ 1, the finite depth has little effect on the evanescent modes except for a doubling of the motion amplitude near the bottom, as the vertical profiles of the form exp(Kz) are replaced by cosh(Kz + Kh)/ cosh(Kh). This is similar to the finite depth effect on linear wave motions. However, the propagating modes radiated by the surface will now undergo multiple reflections at the bottom and sea surface, as shown in figure 5. The oceanic acoustic field is tightly coupled to elastic waves in the crust through these reflections.

One of the greatest complications induced by the presence of a bottom is the heterogeneity of the sediment and rock layers below the water column. The natural layering of the crust has a strong influence on the sound reflection and the nature of the seismic modes (e.g. [START_REF] Latham | Seismic measurements on the ocean floor[END_REF][START_REF] Abramovici | Diagnostic diagrams and transfer functions for oceanic wave-guides[END_REF]). These effects will not be considered here, and we follow exactly the theoretical setting of [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF].

Elastic wave theory

For simplicity, we assume here that the ocean bottom is a uniform and semi-infinite solid, with constant density ρ s and compression and shear wave velocities α c and β given from the Lame coefficients λ and µ,

α 2 c = λ + 2µ ρ s , β 2 = µ ρ s . (4.1a,b)
Assuming that the crust is a Poisson solid, µ = λ and α c = √ 3β. The free wave problem in these conditions was solved by [START_REF] Stoneley | The effect of the ocean on Rayleigh waves[END_REF]. Here we consider the forced problem treated by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF] with a forcing by a pressure field p 2,surf at the sea surface, but now generalized to an additional wave-induced bottom pressure p 2,bot .

The equations of motion in the water column are unchanged from the previous section, but they are now coupled to the elastic motions of the crust. Crustal motions can be separated in an irrotational part with a velocity potential φ c and a rotational part with a stream function ψ, both solutions to Laplace's equation. With a wave source at the surface, a horizontal propagating wave of phase Θ = Kx -ω s t, implies that φ and ψ are either decaying or propagating downwards. They must therefore take the following form, φ c = Ae m(z+h) e iΘ , ψ = Be n(z+h) e iΘ .

(4.2a,b)

Both vertical wavenumbers m and n are given by the Fourier transform of eq. (2.13),

m = K 2 - ω 2 s α 2 c , n = K 2 - ω 2 s β 2 (4.3a,b)
where the sound speed in water has been replaced by the compression and shear velocities.

For K > ω s /β s , m and n are real and both compression and shear waves are evanescent.

For ω s /α c < K < β s the compression wave is evanescent but there is a shear (S) wave that propagates through the crust.

The constants A and B have dimensions of m 2 /s and are determined by the boundary conditions at the ocean bottom.

Horizontal and vertical ground displacements are given by the real parts of

ξ x = KAe m(z+h) + inBe n(z+h) e iΘ /ω s (4.4) ξ z = -imAe m(z+h) + KBe n(z+h) e iΘ /ω s . (4.5)
Hooke's law of elasticity gives

τ zz = λ ∂ξ x ∂x + ∂ξ z ∂z + 2µ ∂ξ z ∂z , τ xz = µ ∂ξ x ∂z + ∂ξ z ∂x (4.6a,b)
The zero tangential stress on the ocean bottom τ xz (z = -h) = 0 yields the following relationship between A and B, which is typical of Rayleigh waves,

B = 2iKm n 2 + K 2 A. (4.7)
Thus, in addition to the unknown water-side amplitudes C and D of the velocity potential at the bottom, we have one more unknown, the compression wave amplitude A on the solid side.

The three equations that relate C, D and A are: the combined kinematic and dynamic boundary condition (2.15), and the bottom continuity of normal velocity

∂φ 2 ∂z = ∂φ c ∂z + ∂ψ ∂x at z = -h (4.8)
and normal stress,

-τ zz (-h) = p(-h) = ρ w ∂φ 2 (-h) ∂t + p 2,bot . (4.9)
For waves in intermediate or shallow water, i.e. kh < π, the Bernoulli term p 2,bot that was not considered by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF] should be included. We thus obtain the linear system of equations

-ω 2 s ρ w + gil e ilh C + -ω 2 s ρ w -gil e -ilh D = iω s ρ w p 2,surf (K, f s ) (4.10) qA -ilC +ilD = 0 (4.11) rA -iω s ρ w C -iω s ρ w D = -p 2,bot (K, f s ) (4.12) with q = mω 2 s ω 2 s -2K 2 β 2 (4.13) r = i ω s ρ s - 4β 4 K 2 mn ω 2 s -2K 2 β 2 + ω 2 s -2K 2 β 2 = i ω s -ρ s m 2 α 2 c + λk 2 + 4µ K 2 mn n 2 + K 2 (4.14)
Since we are in the range where l < ω s /α w we may neglect again g|l|/ω 2 s , which is less than 0.1 for OSGW periods less than 180 s. We rewrite these equations in matrix form,

M [A, C, D] T = [-p 2,surf (K, f s ), 0, -p 2,bot (K, f s )] T (4.15) with M =   0 iρ w ω s e ilh iρ w ω s e -ilh q -il il r -iω s ρ w -iω s ρ w   (4.16)
The general solution of eq. ( 4.15) is the sum of one particular forced solution and the general solution of the homogeneous system, without the right hand side forcing, i.e. the free waves. This combination of free and forced waves is fully determined by the initial conditions.

The forced solution is readily expressed using the determinant of the system det(M ) = -ρ w ω s le ilh r + ω s ρ w e -ilh q -le -ilh r -ω s ρ w e ilh q (4.17 (4.18) in the form

) = 2ρ w (ω 2 s -2K 2 β 2 ) lρ s cos(lh) 4β 4 K 2 mn -ω 2 s -2K 2 β 2 2 -ρ w m sin(lh)ω 4 s ,
A = 2lρ w ω s p 2,surf (K, f s ) + cos(lh) p 2,surf (K, f s ) det(M ) (4.19a) C = -i (lr + qω s ρ w ) p 2,surf (K, f s ) + qe -ilh p 2,bot (K, f s ) det(M ) (4.19b) D = i (qω s ρ w -lr) (qω s ρ w -lr) p 2,surf (K, f s ) -qe -ilh p 2,bot (K, f s ) det(M ) . (4.19c)
As detailed below, this determinant vanishes for the pairs (ω s , K) that fall on the dispersion relation of Rayleigh waves. But the solution for random waves can always be obtained by integrating across this singularity, following [START_REF] Hasselmann | On the non-linear energy transfer in a gravity wave spectrum, part 1: general theory[END_REF].

For kh ≫ 1 we may neglect p 2,bot and,

C = p 2,surf (K, f s ) ρ w ω s ir ′ -q ′ 2r ′ cos(lh) -2q ′ sin(lh) (4.20a) D = p 2,surf (K, f s ) ρ w ω s ir ′ + q ′ 2r ′ cos(lh) -2q ′ sin(lh) (4.20b) with q ′ = ρ w ω s q, r ′ = ilr. (4.21a,b)

Acoustic noise in a bounded ocean

Taking typical values of the water and crust density and sound speeds gives r ′ /q ′ > 70 for the free modes that are significantly generated by the waves (i.e. c j > 0.1, as defined below by eq. 4.34). We can thus consider that q ′ /r ′ ≪ 1, which gives C ≃ D ≃ p 2,surf (K, f s )/[2ρ w ω s (cos(lh) + q ′ /r ′ )]. The velocity potential and pressure in the water are given by,

φ 2 ≃ i p 2,surf (K, f s ) cos[l(z + h) -q ′ /r ′ ] ρ w ω s (cos(lh) -q ′ /r ′ ) e iΘ , p 2 ≃ i p 2,surf (K, f s ) cos[l(z + h)] -q ′ /r ′ (cos(lh) -q ′ /r ′ ) e iΘ . (4.22a,b)
The small but finite factor q ′ /r ′ ensures that the solution remains finite as a small fraction of the acoustic energy is radiated into the crust, otherwise the acoustic energy would accumulate in the water column. The pressure oscillations are thus maximum for resonant frequencies such that the ratio of the water depth and vertical wavelength lh/(2π) is 1/4, 3/4, 5/4 ... We note that the vertical wavelength 2π/l is always greater than the acoustic wavelength 2π/ √ K 2 + l 2 = f s α. As a result, the resonant frequencies are shifted to higher values compared to the vertical resonant condition that is given by f s hα w = 1/4 ...

We may now integrate the pressure spectrum for all acoustic wavenumbers to find, again, the frequency spectrum,

F p2,ap (f s ) = 2πρ 2 w g 2 f s E 2 (f )I(f ) ωs/αw 0 cos(lz + lh) -q ′ /r ′ cos(lh) -q ′ /r ′ 2 KdK. (4.23)
We now have a depth dependence of the sound spectrum. At the surface it is equal to the widely used unbounded ocean value [START_REF] Lloyd | Underwater sound from surface waves according to the Lighthill-Ribner theory[END_REF], but in the water column it can be strongly amplified at depths where cos(lh) approaches 0, which includes the near-bottom region where all the resonant modes have an anti-node.

K (rad/m) K a K 1 K R Re(tan(lh)) Im(tan(lh)) RHS 123456748 Figure 7
. Wavenumbers of seismic and acoustic modes for a fixed frequency, which are the solutions of eq. (4.24). The right hand side of that equation is abbreviated as RHS. Here we used h = 5000 m, αw = 1500 m/s, β = 3200m/s, ρw = 1000 kg/m 3 and ρs = 2500 kg/m 3 . This graph corresponds to a seismic frequency fs = 0.263 Hz, for which fsh/αw = 0.88. Here the determinant has three real roots, one acoustic mode Ka for which l = 0 (horizontal propagation) and two seismic modes K0 and K1. The number of these roots increases with the frequency. One new solution appears every time 2πfs becomes larger than ω 0 s,j defined by eq. (4.29). For reference we also give the wavenumber for which the right hand side of eq. (4.24) is zero, which is the Rayleigh wavenumber KR in the absence of the ocean layer (h = 0).

Rayleigh waves

In order to simplify the algebra, we consider in this section waves in deep water and neglect bottom forcing. As a result, this section brings no new results compared to [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF], but the properties are discussed in more detail and more explicit expressions are presented that will be used later in the derivation of new solutions for other types of seismic waves.

For a fixed frequency ω s , there is at least one wavenumber K for which det(M ) = 0. This condition defines the dispersion relation of the Rayleigh modes [START_REF] Stoneley | The effect of the ocean on Rayleigh waves[END_REF],

tan (lh) = lρ s mρ w × 4β 4 K 2 mn -ω 2 s -2K 2 β 2 2 ω 4 s (4.24)
with the fundamental mode corresponding to the largest K value. Figure 7 illustrates the family of modes for a give frequency.

The phase speeds of the Rayleigh modes vary continuously from the shear wave velocity β, in the limit n = 0 where the shear waves transition from propagating to evanescent, to the sound speed in water α w , in the limit l = 0 where the acoustic modes become evanescent in the water (figure 8.a). This variation of the phase speed has an inflexion point close to the phase speed of Rayleigh waves without the water layer, corresponding to a maximum in group speed (figure 8.b).

We may now use eq. (4.19) to obtain the ground displacement amplitude as a function of the amplitude of the sea surface pressure, with the transfer function

δ ≡ ξ z (-h) = G(K x , K y , ω s ) p 2,surf (K x , K y , ω s ) (4.25) f s h / 1 w 1234564755869AB4C DEF764755869AB4C 1000
G = 2iρ w lm(K 2 -n 2 ) (n 2 + K 2 ) det(M ) . (4.26)
An example of this transfer function is shown in figure 9, in which the dispersion relation of the Rayleigh modes appear as the narrow red bands around the singularities of G.

To evaluate the complete family of Rayleigh wave solutions with dispersion relations ω s,j (K), where j is the mode number, we need to examine the nature of these singularities. We re-write the determinant of the system as,

det(M ) = 2ilρ s ρ w cos(lh) (ω 2 s -2K 2 β 2 ) 4β 2 k 2 mn + ω 2 s -β 2 K 2 -ω 4 s (1 + tan(lh)ρ w m/ρ s l) .
(4.27) For most frequencies, the singularities are simple, allowing a Taylor expansion of G in this form

G(K, ω s ) ≃ G ′ (K j (ω s )) ω 2 s -ω 2 s,j (K)
.

(4.28)

However, for each mode j, there is a critical frequency

ω 0 s,j = 1 h arctan(- l m ρ s ρ w ) + jπ (4.29)
for which both tan(lh)ρ w m/ρ s l = 1, and ω 2 s -K 2 β 2 = 0. These critical frequencies ω 0 s,j exist for modes j > 0. These frequencies are those for which new Rayleigh modes appear, in a way similar to the Love waves discussed by Aki & Richards (2002, figure 7.3). The singularities at (K, ω s ) = (ω 0 s,j /β, ω 0 s,j ) are not simple and for these we have G(ω s , k) ≃ G ′ (K j (ω s ))/ K 2 -K 2 j (ω s ). Solutions for the vertical displacement at the top of the crust have a spectral amplitude which is linearly related to the equivalent surface pressure amplitude,

δ(K x , K y , ω s ) = G(K x , K y , ω s ) p 2,surf (K x , K y , ω s ).
(4.30)

For the simple singularities of G, we may write G(K, ω s ) = G ′ (K, ω s,j )/(ω 2 s -ω 2 s,j ) + O(ω s -ω s,j ). Taking initial conditions δ(t = 0) = 0 and ∂δ/∂t = 0 gives the full solution (Hasselmann 1962, eq. 3.2) of the j th Rayleigh mode response

δ(K x , K y , ω s ) = G ′ (K, ω s,j )e iKx 1 ω 2 s,j -ω 2 s e -iωst - 1 2ω s,j
e iωs,jt ω s,j -ω s + e iωs,j t ω s,j + ω s , (4.31) as a function of K, where ω s,j (K) is frequency of the j th Rayleigh mode. For a forcing that varies slowly on the scale of the seismic period 2π/ω s , and provided that the forcing spectrum is continuous in spectral space, this gives the rate of change of the ground displacement spectrum given by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF],

∂F δ (K x , K y ) ∂t = S DF (K x , K y ) = π |G ′ | 2 2ω 2 s,j F p2,surf (K x , K y , ω s ). (4.32)
The other discrete singularities, at (K, ω s ) = (ω 0 s,j /β, ω 0 s,j ), are associated with 'conical' or 'head' waves [START_REF] Aki | Quantitative seismology[END_REF], for which the vertical wavenumber n = 0, and that propagate along the ocean-crust interface. These horizontally propagating shear waves are generated with an evanescent compression wave. The singularity is integrable over the two spectral dimensions ω s and K.

Although it looks like only the resonant forcing contributes to the solution, it is in fact the near-resonant forcing (ω s ≃ ω s,j ) that builds up the seismic noise, because the exact resonant terms have a zero measure in spectral space. Indeed, a purely resonant forcing would give an amplitude that increases linearly with time, and an energy that increases like t 2 . The linear growth of energy in time can be interpreted as an effect of the narrowing with time of the frequency bandwidth in which the interaction is significant. This is a general property of wave-wave interactions (see also [START_REF] Hasselmann | Feynman diagrams and interaction rules of wave-wave scattering processes[END_REF].

We express the source of seismic noise with the rate of increase of the variance of δ per unit of propagation distance. It is

S DF (ω s ) = K(s)S DF (K x , K y ) U 2 = 4π 2 f s c 2 j β 5 ρ 2 s F p2,surf (K x , K y , ω s ) (4.33)
where U is the group speed of the seismic waves, and c j is a dimensionless coefficient that depends on ω s h/α w and the seismic mode index j, shown in figure 10,

c 2 j = β 5 ρ 2 s K j U 2 j 2πω s π G ′ j 2 2ω 2 s . (4.34)
We note that a missing 2π in eq. ( 5) of [START_REF] Ardhuin | Ocean wave sources of seismic noise[END_REF] has been corrected here.

A very rough simplification can be obtained by taking U and K independent of j and h. Then S DF (f s ) can be taken as a sum of all Rayleigh modes in the form

S DF (f s ) = 4π 2 f s β 5 ρ 2 s ∞ i=0 c 2 j F p (K ≃ 0, f s ). (4.35)
Values of c j are obtained from eq. (4.34), and shown in figure 10. We may propagate these sources of seismic waves in a vertically symmetric earth model, neglecting all three-dimensional propagation effects, and parameterizing seismic wave scattering and dissipation with a uniform quality factor Q. Under these assumptions, the spectral density of the vertical ground displacement at z = -h and at the longitude λ O . Dimensionless coefficients cj that amplify the wave-induced pressure into ground displacement. The maxima of cj correspond to quarter-wavelength resonance typical of organ pipes, except that here the sound waves propagate obliquely in the water column, which is why the maxima are at values of fsh/αw which are not exactly at 1/4, 3/4..., but shifted to higher frequencies by a factor l/ (K 2 +l 2 ). The amplitudes of the peaks depend on the impedance ratio of the sea water and crust. Hence the peak amplitude increases with ρsβ/(ρwαw). For example, β = 2800 m s -1 , gives a maximum of 0.88 for c0 instead of 1.03 here with β = 3000 m s -1 .

and latitude φ

O F δ (λ O , φ O , f s ) = π/2 -π/2 2π 0 S DF (f s ) R E sin ∆ e -2πfs∆RE/(UQ) (R 2 E sin φ ′ O dλ ′ O dφ ′ O ) (4.36)
with R E the Earth radius, and U the seismic group velocity. The term (R 2

E sin φ ′ O dλ ′ O dφ ′ O )
is the Earth surface area element. The denominator (R E sin ∆) is the geometrical spreading factor for wave energy that follows geodesics on the sphere (e.g. [START_REF] Kanamori | Use of long-period surface waves for rapid determination of earthquake-source parameters[END_REF], replacing the distance (R E ∆) used in flat Earth models (e.g. [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF].

The Rayleigh waves thus generated propagate away like free modes, those that exist without the local forcing. For these free modes with a monochromatic ground displacement δ(x, y, t) the surface pressure is constant so that C = -D and the velocity potential and pressure in the water take the form,

φ 2 (z) = ω s l cos(lh) sin(lz)δ(x, y, t) (4.37a,b) p 2 (z) = iρ w ω 2 s l cos(lh) sin(lz)δ(x, y, t). (4.37a,b)
For varying water depths, and on land, one may assume that the seismic energy is propagated along a ray and apply the refraction coefficient given by [START_REF] Hasselmann | A statistical analysis of the generation of microseisms[END_REF].

Seismic P and S waves

Unlike the Rayleigh waves that grow resonantly in time owing to the trapping in the ocean/crust waveguide, the P and S waves, for which the vertical wavenumbers m or n are complex, radiate into the earth's interior, and their energy level is given directly by its value at the source. In particular, for K < ω s /α c we have propagating P waves with a velocity potential amplitude A given by eq. (4.19). In the particular case K = 0, which correspond to standing waves, we have only have P waves, no S waves or Rayleigh waves, and these propagate exactly along the vertical axis. The frequency spectrum of vertical ground displacements at (z = -h) can be evaluated directly with (4.30) because G has no singularity in this range of wavenumbers,

F δ,P (f s ) = f s E 2 (f )I(f ) ρ 2 w g 2 ρ 2 s β 4 s c 2 P (4.38)
with a non-dimensional coefficient c P ,

c 2 P = 2π ωs/αc 0 4l 2 m 2 ρ 2 s β 4 s ω 2 s det 2 (M )
KdK.

(4.39)

A similar expression can be written for S waves, and both are illustrated in figure 11,

F δ,S (f s ) = f s E 2 (f )I(f ) ρ 2 w g 2 ρ 2 s β 4 s c 2 S . (4.40)
However, in the range of wavenumbers where S waves exist, k < ω s /β, there can also be evanescent P waves, and the system can approach the singularity for ω s = ω s,j and k = ω s /β. We evaluated numerically the coefficient

c 2 S = 2π ωs/β 0 4l 2 m 2 k 2 ρ 2 s β 4 s ω 2 s (n 2 + k 2 ) det 2 (M )
KdK.

(4.41)

Due to the typically three times stronger attenuation of S waves compared to P waves in the Earth mantle (e.g. [START_REF] Anderson | Q of the Earth[END_REF][START_REF] Pasyanos | A simultaneous multiphase approach to determine p-wave and s-wave attenuation of the crust and upper mantle[END_REF], we will now focus on P waves only, which should dominate in the far field of the noise source.

For the estimation of the spectrum recorded outside of a source area, it is more convenient to express the local seismic source as a function of the horizontal propagation angle θ, and the vertical take-off angle ϕ. For P waves, this gives,

F δ,P (f s , θ, ϕ) = f s E 2 (f )I(f ) 1 ρ 2 s β 4 s c 2 P,ϕ sin ϕ (4.42)
with the non-dimensional coefficient c P,ϕ defined by

c 2 P,ϕ = 4l 2 m 2 ρ 2 s β 4 s ω s α c det 2 (M ) ∂K ∂ϕ , (4.43)
which is the normalized source per unit solid angle Ω, so that the average over the half space of downward directions Ω -is

c 2 P = 2π 0 π/2 0 c 2 P,ϕ sin ϕdϕdθ = Ω - c 2 P,ϕ dΩ (4.44)
as defined by eq. (4.39). It is noteworthy that the distribution of the P -wave energy with the take-off angle is very close to the one given by a small disk pushing at the top of a uniform half space, as given by [START_REF] Miller | On the partition of energy between elastic waves in a semi-infinite solid[END_REF] and used by [START_REF] Vinnik | Sources of microseismic P waves[END_REF], although it also varies with the non-dimensional water depth f s h/α w . The only missing item in the work by [START_REF] Vinnik | Sources of microseismic P waves[END_REF] is the very strong amplification of the motion for resonant frequencies associated with the water layer. Due to the large impedance contrast at the water-crust interface the relative amplification of P waves is one order of magnitude stronger than for Rayleigh waves. We thus expect a much tighter correspondence of the strong seismic noise sources with the water depths that correspond to a maximum amplification.

Observable P wave spectra

We will now finish our analysis of these body waves by estimating the ground motion due to P and Rayleigh waves as a function of the distance from the source, an application of practical interest. For a seismic station or hydrophone in the ocean, the incoming energy at the receiver will amplified by multiple bottom and surface reflections in the water column. However, this P -wave signal from remote sources is likely to be dwarfed by locally generated noise.

For a land-based station, we may assume that P waves arrive directly from the source area. There are many other seismic wave phases that have undergone multiple reflections at the surface, these are called P P , P P P ..., or at any inner boundary of the Earth, like P kP phases that have reflected off the mantle-core interface [START_REF] Aki | Quantitative seismology[END_REF]. These can be treated exactly like the direct P phase.

On arrival at the receiver, these waves are totally reflected at the Earth surface at z = 0, which doubles the ground motion, so that the ground displacement is given by the integral of four times the incoming spectral densities over the directions θ and ϕ, which can be replaced by an integral over the source positions (λ S , φ S ). The transformation from the ray parameters (θ, ϕ) to the geographical coordinate can be obtained approximately for any type of seismic wave using travel time tables (e.g. [START_REF] Snoke | Traveltime tables for iasp91 and ak135[END_REF]), which also provide the travel time τ .

We may now express the ground displacement due to P waves at the observing station of coordinates (λ O , φ O ), as a function of the same quantity at the location of sources, as given by eq. (4.42)

F δ,P (λ O , φ O , f s ) = Ω + 4F δ,P (λ, φ, f s , θ, ϕ)e -2πfsτ (λS ,φS)/Q dΩ, (4.45)
with dΩ an element of solid angle that corresponds to the ensemble of rays arriving from an Earth surface element around the sources located at (λ S , φ S ). This sum may also be transformed as an integral over the ocean surface by properly mapping Ω to (λ S , φ S ). The elementary solid angle dΩ is zero for the so-called shadow zones, the regions for which there is no P -wave ray that connects to the observing station. For a single phase of seismic waves, this ensemble of rays has a half-banana shape, as illustrated on figure 13.

From our calculations, we expect that P waves will dominate the signal at large distances from the source. The exact location where P -wave levels overtakes Rayleigh-wave levels depends on the relative attenuation of the two types of waves. With a realistic Q = 2000 for the P waves, and Q = 400 for the Rayleigh waves, figure 14 shows that it occurs at an epicentric angle of 40 • , which is a distance of 4400 km, consistent with the observations reported by [START_REF] Vinnik | Sources of microseismic P waves[END_REF] using Kazakhstan array data.

Conclusions

We have shown how the same physical process, the interaction of ocean surface gravity wave (OSGW) trains, can produce a wide variety of noises, in the atmosphere, ocean and Earth's crust, that can be classified according to their horizontal phase speed. The slowest noises in the ocean are acoustic-gravity waves that dominate pressure records at depths less than about one tenth of the acoustic wavelength. These acoustic-gravity waves cannot exist in the absence of OSGWs and are thus confined to the region of active wave forcing. Intermediate phase speeds correspond to Rayleigh waves that contain most of the energy of the seismic modes for distances less than about 4000 km from the source.

We corrected previously published asymptotic behaviour for very long period noise (T > 30 s). In particular we find that the sources of this long period noise are attenuated on the continental shelves, consistent with previous studies of forced gravity wave motion. This finding supports a spatial distribution of these sources outside of the continental

Figure 1 .

 1 Figure1. (a) The vertical evanescent or propagating nature of the noise field in the solid and liquid layers is defined by the horizontal phase speeds relative to the distinct values of the sound speed in the ocean (αw), and the shear (β) and compression (αc) speeds in the crust. From slow to fast, there are the acoustic-gravity (AG) domain, the Rayleigh (R) wave domain, and two body wave domains (S only, and P and S together) (b) For any fixed frequency, the four domains correspond to 4 concentric regions in the wavenumber plane. For three selected noise frequencies generated by OSGW in the infragravity, dominant and high frequency ranges of the forcing wave field, the limiting wavelengths between the four domains are indicated, using αw = 1.5 km/s, β = 3.2 km/s, αc = 5.54 km/s. One example of interaction (black vectors) is shown with two gravity wave modes that interact to generate a Rayleigh wave.

  Figure4. Compared to deep water, the source of seismic noise power in finite water depth is amplified by a factor tanh 2 (kh) [1 + 2kh/sinh(2kh)], when accounting for both the surface and bottom forcing. This is very different from the approximation that considers the surface forcing only (leading to eq. 2.28), even more so when using its deep water approximation (eq. 2.22) which gives an additional factor 4 difference for kh ≪ 1.

  Figure5. Schematic of wave groups and forced acoustic and seismic wave motion. For a given wave group period T , the horizontal wavelength L = 2π/K can be larger than the acoustic wavelength in the water due to the oblique incidence of the sound waves. The superposition of two obliquely propagating sound waves (arrows) forms a mode pattern that propagates horizontally at a supersonic speed. The acoustic wavelength is La = 2π/ka = 2π/ √ K 2 + l 2 = αw/T . Both vertical and horizontal wavelengths are larger. For readability, the wave and bottom amplitude is not to scale, and we have reduced the number of waves in the group from 104 to 10. Other than that, the angles are preserved. The configuration shown here corresponds to the conditions for maximum amplification of mode 1 (see below), with a vertical wavelength to water depth ratio of 0.75. For a water depth h = 4400 m depth, this corresponds to L = 7.7 km and a seismic frequency fs = 0.29 Hz.

Figure 6 .

 6 Figure6. Example of (a) directional wave spectrum and (b) resulting profiles of the different contributions to the pressure fluctuations in the ocean, assuming infinite water depth. The ratio of double frequency to linear wave contributions depends on the amplitude of the waves and on the directional spectral shape, because all double frequency contributions are proportional to E 2 (f )I(f ). This directional spectrum was estimated with a numerical wave model, and corresponds to the loudest noise event recorded at the ocean bottom seismometer H2O, on 31 May 2002 at 25 N, 136 W. This unusual spectrum has large wave energies in opposite directions, radiated from a North Pacific storm and Hurricane Alma (This event is analyzed in detail by[START_REF] Obrebski | How moderate sea states can generate loud seismic noise in the deep ocean[END_REF].

Figure 8 .

 8 Figure 8. Phase (top) and group (bottom) speeds of Rayleigh waves as a function of the dimensionless water depth. Calculations used αw = 1470 m s -1 , β = 3000m s -1 , ρw = 1000 and ρs = 2500 kg m -3 . Values without the water layer are also indicated with the blue line. In that case Rayleigh waves are not dispersive.

Figure 9 .

 9 Figure 9. Magnitude of the transfer function G(K, fs) illustrating the singularities along the dispersion relations of the free Rayleigh modes with wavenumbers K = Kn(fs), that correspond to det(M ) = 0. The oblique dashed lines correspond to phase speeds equal to αc, β, and αw, and separate the four domains body waves (P+S), mixed body and evanescent waves (S), Rayleigh waves (R) and acoustic-gravity modes (A-G).

  Figure10. Dimensionless coefficients cj that amplify the wave-induced pressure into ground displacement. The maxima of cj correspond to quarter-wavelength resonance typical of organ pipes, except that here the sound waves propagate obliquely in the water column, which is why the maxima are at values of fsh/αw which are not exactly at 1/4, 3/4..., but shifted to higher frequencies by a factor l/ (K 2 +l 2 ). The amplitudes of the peaks depend on the impedance ratio of the sea water and crust. Hence the peak amplitude increases with ρsβ/(ρwαw). For example, β = 2800 m s -1 , gives a maximum of 0.88 for c0 instead of 1.03 here with β = 3000 m s -1 .

Figure 11 .

 11 Figure 11. Non-dimensionless coefficients cP and cS that amplify the wave-induced pressure into ground displacement associated with P and S waves.

Figure 12 .

 12 Figure12. Dimensionless coefficients cP,ϕ that amplify the wave-induced pressure into ground displacement. The maximum for a zero take-off angle corresponds to vertically propagating compression waves, and the compression waves that propagate along the crust have a vanishing amplitude.

Figure 13 .

 13 Figure13. Schematic of rays for seismic P waves radiated from a point source (red dashed lines) with a directional distribution cP,ϕ (solid red line), and ensemble of rays received at a given station from an extended noise source (yellow half-banana).

Figure 14 .

 14 Figure14. Estimates of the rms vertical ground displacement associated with Rayleigh or P waves, as a function of the epicentric angle ∆, for a source of intensity F p2,surf (K = 0, f )df = 4.2 × 10 4 hPa 2 m 2 over a 330 by 330 km square, assuming an attenuation factor Q = 2000 for the P waves[START_REF] Pasyanos | A simultaneous multiphase approach to determine p-wave and s-wave attenuation of the crust and upper mantle[END_REF], with travel times given by the ak135 reference Earth model[START_REF] Snoke | Traveltime tables for iasp91 and ak135[END_REF].

  This gives, for z = -h,C = -D z (sσ + s ′ σ ′ ) 2 g [gK tanh(Kh) -(sσ + s ′ σ ′ ) 2 ] + D pb (z = -h)

Table 1 .

 1 Summary of theoretical expressions for the gravity, acoustic and seismic noise fields.

		infinite water depth finite water depth h
	Surface forcing field	eq. (2.29)	eq. (2.28)
	Bottom forcing field	n.a.	eq. (2.31)
	gravity-like modes	eq. (3.11)	not given
	acoustic spectrum	eq. (3.10)	eq. (4.23)
	microbaroms	eq. (3.19)	eq. (3.19)
	seismic source (R)	n.a.	eq. (4.32)
	seismic source (P)	n.a.	eq. (4.38)
	seismic source (S)	n.a.	eq. (4.40)
	seismic spectrum (R)	n.a.	eq. (4.36)
	seismic spectrum (P)	n.a.	eq. (4.45)

Table 2 .

 2 critical values of ωs for which new modes appear when ωs increases Table of arab and greek symbols Appendix B. Definition of symbols Rayleigh wave mode number, counted from 0 k and k ′ vector wavenumbers of surface gravity waves K horizontal vector wavenumbers of acoustic or seismic waves l, la, m and n vertical wavenumbers for φ2,φ2,a , φc, and ψ p pressure p 2,surf wave-induced forcing at the sea surface eq. (2.16) p2,B Bernoulli head pressure eq. (3.7) p 2,bot wave-induced forcing at the bottom, p 2,bot = p2,B(z = -h) eq. (2.31) M matrix of the linear system of equations Q seismic quality factor (i.e. damping coefficient)

	Symbol	meaning	where
	1 and 2	indices denoting first and second order motions
	αw	sound speed in water
	αc	compression wave speed in crust
	αa	sound speed in atmosphere
	β	shear wave speed in crust
	∆	angular distance
	δ ζ	vertical ground displacement at the top of the crust, δ = ξz(z = -h) elevation of the sea surface
	λ and µ	Lame elastic coefficients for the crust
	λO and φO	longitude and latitude of the observation location
	λS and φS	longitude and latitude of a noise source
	ρ	perturbation of density
	ρw	mean water density
	σ	radian frequency of surface gravity waves σ = 2πf
	τ	stress tensor
	Θ	phase function of the seismic or acoustic waves
	φ, φa, φc	velocity potentials in the water, atmosphere, crust
	ψ	stream function in the crust
	ϕ	take-off angle for seismic body waves
	ξ	displacement of particles
	ωs	radian frequency of noise ωs = 2πfs
	ω 0 s,j

Table 3 .

 3 Table of roman notations

† This amplification depends not only on the impedance ratio of the water and crust, which defines the amplification coefficients cj derived below, but also on the seismic attenuation coefficient Q, which is discussed in section 4. Realistic calculations following[START_REF] Ardhuin | Ocean wave sources of seismic noise[END_REF] typically give a factor 10 to 20 amplification of the sound in the water column due to the bottom elasticity.

Appendix A. Coupling coefficients

Using the coupling coefficient D given by Hasselmann (1962, eq. 4.3) for the velocity potentials, our coupling coefficient for the elevation amplitudes is

In the bottom pressure, the additional term arising from the orbital velocity has a coupling coefficient

The relationship with the coupling coefficient C given by Herbers & Guza (1991, their eq. 4) for the bottom pressure, expressed in meters of water, is given by solving eq. (2.15) for φ 2 , and then rewriting Bernoulli's equation (2.4), as


M Naaim

Y D Rand

DENSE AVALANCHE FRICTION COEFFICIENTS INFLUENCE OF NIVOLOGICAL PARAMETERS

High altitudes, heavy snow, and sharp mountainous relief of Mont Blanc massif make the upper part of the Arve valley (Chamonix, France) one of the most threatened areas by snow avalanches. 115 paths from this area were included in the avalanches national observatory, 26 of them have experienced between 20 and 100 events over a period of almost 100 years, among these avalanches 729 are posterior to 1958. Recently, Meteo France reconstructed, for all the French massifs, the daily evolution of the snowpack properties since 1958 producing a consistent data set with available measurements, advanced meteorological models and snow physics. These two data sources were combined to calibrate an avalanche dynamics model including sub models for friction, erosion, entrainment and deposition. First topographic profiles of the 26 well documented paths were built using an accurate laser scan model. Then, the snow conditions were determined from the reconstructed snow packs for each of the 729 avalanches. Finally the full range of friction coefficients was scanned and a numerical simulation was performed for each pair of friction parameters and thus the corresponding run-out altitude determined. Only the pairs of parameters, for which the run out altitude is found close enough to the observed one (±5m), was retained. Statistical methods such as CPA and stepwise were used to investigate correlations between the obtained friction coefficients and the snow properties. Concerning the static friction coefficient, an increasing tendency with the temperature and the density was evidenced as well as a decreasing tendency with the liquid water content and the snow depth. 1

INTRODUCTION

The propagation of snow avalanches displays highly complex features. Longitudinally, the avalanche is made of a front followed by a core and a tail. The flow may present a strong vertical stratification of velocity and density. In addition, the nivological occurrence conditions and the morphology of avalanche paths vary widely making the avalanche dynamics, and consequently its modeling, a very difficult matter.

AVALANCHE MODELING

Despite this complexity several simple models were proposed during the last century to capture the main features of snow avalanches (see [START_REF] Harbitz | A survey of computational models for snow avalanche motion[END_REF]. The majority has focused on the dense flowing part.

* Corresponding author:

Mohamed Naaim, ETNA-IRSTEA, 2 Rue de la papeterie, BP 76, 38402 Saint Martin d'Hères, France.

: 33 4 76762722. email : mohamed.naaim@irstea.fr One of the first models for dense snow avalanche has been proposed by Voellmy in 1955. This onedimensional model has assimilated the avalanche to a sliding snow block, submitted to the sum of a Coulomb friction and a dynamic drag. This last is proportional to the square of the velocity. Under quasi similar assumptions, several other authors have developed various models that predict velocities along the path and run-outs of the center of mass of the avalanche. To account for the flow depth variations with time and space during an avalanche event, the hydraulic framework has been adopted and a second generation of models, based on depth-averaged hydraulic framework, has emerged. It incorporates the two-parameter friction model of Voellmy and several other terms to deal with erosion and deposition processes and with vertical (active / passive) pressure distribution.

In addition to the velocities and run-outs, this second class of models allows predicting the flow depth along the path as well as the depth of the deposit. These models have been developed in one and two dimensional by numerous authors.

The model used in the following concerns the dense avalanche flows and belongs to the hydraulic models class. Because the ratio of the height h of an avalanche to its length L, is small enough ( ≪ 1), the shallow water formalism has been extended to dense avalanches. The shear, (and consequently the resulting dissipation), is mainly located at the interface between the moving snow, and the snow or the soil at rest at the base of the flow, in accordance with the measured velocity profiles, published in the literature [START_REF] Dent | Density, velocity and friction measurements in dry snow avalanche[END_REF] ).

In this model the avalanche path is represented by its longitudinal profile z(x) where z is the altitude and x the longitudinal horizontal coordinate. The "wetted" cross-section S(x) is given as a function of the flowing snow depth h ( ( ) = ( )ℎ ),

where k(x) and α are parameters representing the local traverse morphology of the path, parameters to be calibrated from the local crosswise section.

Model equations

The equations of the model are the conservation of mass and the conservation of momentum. If θ is the local slope angle and t the time, governing equations are given by the following partial differential equations:

+ = (1)
+ + ! ℎ" = cos & (tan & -μ) (2)
According to the formalism developed in [START_REF] Naaim | Dense snow avalanche modeling: flow, erosion, deposition and obstacle effects[END_REF], a specific distributed erosion and deposition model is used to determine ϕ.

All of the existing models require the specification of values for the friction parameters. Because the physics of flowing snow is still poorly understood, this cannot be achieved with reference to rheological measurements. An empirical calibration is required. This paper summarizes the methodology and the main results of a simple calibration of a hydraulic type model over 729 historical events and addresses the correlation between the calibrated friction coefficients, and the nivological parameters of the involved snow.

Detailed developments and detailed results could be found in the version submitted to CRST.

DATA AND PROCEDURE OF CALIBRATION

Avalanche historical release and run-out altitudes

In the national French avalanche observatory, 115 avalanches paths of Chamonix valley are monitored. The number of events per path is ranged from 5 to 100 and the observation period is ranged from 50 to 100 years. Among many various data, the date, the release altitude and the run-out altitude of each event are recorded. The altitudes of the trigger and the run-out are determined from in-situ observation. The snow avalanches are recorded only if it reaches a predefined threshold.

The magnitude of the recorded avalanches can then be considered moderate to large. In this list, the pair of numbers between brackets attached to the name of each path, correponds to its number of event posterior to 1958 and its total events number.

The longitudinal and traverse profiles of each avalanche path have been determined from 1 m digital terrain model obtained using a laser scan technology.

Avalanche historical nivological data

The avalanche dynamics model used in this study requires accurate knowledge, at the date of the event, of the distributed snow mantle and its vertical structuring. Whether it's remote (satellite, radar) or classic (polls/surveys), currently no operational network allows the knowledge or the reconstruction of the snow cover and its structure at a spatial scale of an avalanche path and at an hourly time scale, especially in times of imminent release of avalanches, where the snow conditions changes significantly. To overcome this difficulty, the snow center of Meteo-France, developed an integrated tool combining assimilation of meteorology data and snow pack models. The spatial scale of this chain of models is the mountain massif (approximately 500 km 2

) and the time scale is few hours. In this chain each massif is represented by several altitude slices of 300m, in six exposures (North, East, South East, South, Southwest and West) and the avalanche diagnostic is made for two slopes 20° and 40°. This modeling allows taking into account the significant effects of altitude, exposure and sloping which are the effects that contribute significantly in the structural differences of the snowpack. This chain (SCM) made of Safran, Crocus and Mepra, was submitted to several validations through comparisons between the snow profiles measured to the simulated ones. Its results in terms of avalanche activity were successfully compared to the observed natural avalanche activity. Recently, runs since 1958 until today were made to build two data sets daily on all the Alps and Pyrenees [START_REF] Durand | Reanalysis of 44 years of climate in the French Alps (1958-2002): methodology, model validation, climatology, and trends for air temperature and precipitation[END_REF]). This simulation campaign produced a consistent snowpack data set with available measurements, advanced meteorological models and snow physics. An important aspect to be highlighted here to support the use of these data in our study (scale of path and valley) is that almost the entire surface of the Mont Blanc massif is occupied by the Chamonix Valley and the magnitude of the events analyzed here are medium to large for which the weather can be considered much more homogenous.

Among the SCM chain data, one has the total snow depth (h t ), the snow density (ρ), the average temperature (T), the average liquid water content (Lwc), and the burial depth of the weakest layer (h), the cohesion (τ c ). These quantities have been generated for each avalanche date, for the main orientation of each path and for all the three hundred meters from 900 m to 3100 m altitude. The initial conditions let say release altitude and initial snow pack (the snow depth above the weakest layer, the total snow depth and the cohesion τ c (x)) was defined. The run out altitude, used as calibration criterion, is known. So, one has all the ingredients to perform a calibration as long as one considers constant the friction coefficients for each event.

The calibration procedure has been as simple and free as possible. One was scanned the range of the static friction parameter from 0.1 to 0.7 in steps of 0.01 and the range of the inertial friction coefficient from 500 to 1500 in steps of 50 m.s -2 . Therefore for each pair (µ 0 , ξ), a numerical simulation was performed and the corresponding run-out altitude was determined and the gap to the actually observed altitude was calculated. The pairs of parameters that minimize this gap (and for which the gap is lower than 10 m) was retained. 874800 simulations have been performed. The simulation campaign has taken 2890 hours on a twelve processors parallel computer.

Inertial friction coefficient

The set of inertial friction coefficients ξ obtained has been analyzed. No obvious correlation was detected neither with temperature nor with the density. The only outcome that could be drawn up from this simulation campaign was the statistical distribution of ξ. As shown in the Figure 3, the statistical distribution of ξ can be approximated by a Gaussian distribution (mean 940 m.s -2 and standard deviation 220 m.s -2

).

Proceedings, 2012 International Snow Science Workshop, Anchorage, Alaska

Static friction coefficients

For each path, one studied the eventual correlations between the obtained static friction coefficients and the corresponding nivological quantities. The stepwise method and the principal component analysis were used to this end. Several trends were observed. For example, the temperature and the total liquid water content have a significant influence on the static friction coefficients for five paths and the initial snow depth have a significant influence for 13 paths. However, no systematic trend was observed due to the number of dependant nivological quantities and the limited number of data per path. One then grouped all the data in a single file and one analyzed, one by one, the influence of the depth, of the density, of the total liquid water content and of the temperature. For each one, the dataset was sampled into 6 to 14 classes of equal amplitude and for each class, the mean and the standard deviation ware calculated for the static friction coefficient and for the nivological parameter studied. The results are presented and commented her after. Before drawing any conclusions, it is essential to keep in mind the level of uncertainty when manipulating field observations. Another important matter to recall is related to the various approximations and assumptions we have adopted for modeling both the dynamics of dense avalanches and the reconstruction of the snowpack. An additional important aspect to recall is the statistical character of data obtained. The coefficients are averaged over significant number of data. The conclusions concern average trends.

This study revealed or confirmed several existing knowledge concerning the avalanches friction parameters dependency to nivological quantities. The main result is the linear increase of the static friction coefficient with the temperature and with the density, even if for this last the trend is observed only for density below 200 (kg.m -3

) (see figures 4) .

Another important result concerns the case of wet snow avalanches (figure 5). The results showed that the liquid water content has no significant effect except for high values, for which the static friction coefficient decreases significantly and may reach the same low value obtained for very low temperature. The last conclusion is related to the snow depth. The results obtained in this work showed a decrease of the static friction coefficient with the snow depth (figure 5). Finally no evident trend or correlation between the inertial friction coefficient and the nivological quantities could be found.

  and many others). In the model used here the friction at the base is represented by the two parameters standard model of Voellmy. The total friction coefficient µ (μ( ) = μ + ) is a sum of a Coulomb friction μ and a dynamic drag coefficient proportional to the square of the Froude number, ( = / ℎ) ratio of the mean velocity U to the gravitational wave velocity ℎ ), and inversely proportional to the inertial friction coefficient ξ (m.s -2

Figure 1 :

 1 Figure 1: Example of path topographic profile and historical run-out altitudes (Taconnaz path) We have focused on the 26 paths well documented. This means paths that have (1) a continuous observation, (2) a good data quality (few obviously wrong data), and (3) the number of events exceeding 20. These paths are listed in the following : Le petit Cervin (24/27), La Fillaz (4/23), Les Vorgeats (28/61), Orthaz (51/101), Les Plants (44/75), Brevent (33/63), Pecheux (40/64), Favre (39/65), Entremene (37/68), Rocher (13/25), Vouillours (38/63), Etrangleur (23/51), Lays (23/44), Egolerons (32/50), Epinettes (20/29), Pichet (21/31), Greppon (21/34), Chapeau (28/43), Grand-Chantet (36/63), Nant-Montets (22/34), Lyapet (33/56), Griaz (30/79), Bourgeat (30/82), Arandellys (7/37), Taconnaz (40/88), Nant-Freynet (12/27).In this list, the pair of numbers between brackets attached to the name of each path, correponds to its number of event posterior to 1958 and its total events number.

Figure 2 :

 2 Figure 2 : Example of nivological data profile 4 BACK ANALYSIS OF HISTORICAL DATA

Figure 3

 3 Figure 3 : ξ ξ ξ ξ statistical distribution

Figure 4 :

 4 Figure 4: mean static friction coefficient µ according to the mean initial snow temperature and mean density

Figure 5 :

 5 Figure 5 : mean static friction coefficient according to the mean initial snow depth and liquid water content 5 CONCLUSIONS The calibration of a two parameters hydraulic avalanche dynamics model has been performed on 729 historical avalanches occurred at Chamonix valley during the last 54 years. All the studied avalanches have reached a predefined alert threshold, criteria to be recorded in the national data base. The magnitude of these avalanches can then be considered as moderate to large. The weather in this valley, especially during avalanche times, has been considered homogeneous and the model initial conditions have been extracted from the reconstructed snow mantle by the nivo meteorological chain SCM. The friction parameters ranges were scanned and for each pair (µ 0 ,ξ) the run-out altitude was determined. The run out altitude served as criteria to calibrate the friction parameters. The correlations of parameters that reproduce historical avalanches with the other available quantities let say snow density, snow depth, snow temperature and snow liquid water content were searched for. The range of each nivological quantity has been split into classes of equal amplitudes. For each class the mean and the standard deviation of the static and inertial friction coefficient have been determined.

Proceedings, 2012 International Snow Science Workshop, Anchorage, Alaska