Matthieu Leclair

Keshav Raja

Chantal Staquet

Nonlinear reflection of a two-dimensional finite-width internal gravity wave onto a slope

 for pure plane waves. The aim of this theory is to determine the conditions under which the incident and the reflected waves form a resonant triad with the second harmonic wave resulting from their interaction. Thorpe's theory leads to an indeterminacy of the second-order amplitude at resonance. In waiving this indeterminacy, we show that the second-order amplitude has a finite behavior at resonance, increasing linearly from the slope. We investigate the influence of background rotation and find similar results with a weaker growth rate. We then adapt the theory to the case of an incident plane wave of finite-width. In this case, nonlinear interactions are confined to the area where the incident and reflected finite-width waves superpose implying that the amplitude of the second-order wave is bounded at resonance. We find good agreement with the results of numerical simulations in a vertical plane as long as the dissipated power of the incident and reflected waves remain smaller than the power transferred to the second harmonic. This is the case for small slope angles. As the slope angle increases, the focusing of the reflected wave enhances viscous effects and dissipation eventually dominates over nonlinear transfer. We finally discuss the relevance of laboratory experiments to assess the validity of these theoretical results.

Introduction

When a plane internal gravity wave impinges on topography, energy transfer to small scales occurs even in a linear regime, as accounted for by a simple geometrical argument proposed by [START_REF] Phillips | Dynamics of the upper ocean[END_REF]. In a stratified rotating fluid, the dispersion relation of internal gravity waves is:

ω 2 = N 2 sin 2 β + f 2 cos 2 β,
(1.1) † Email address for correspondence: chantal.staquet@univ-grenoble-alpes.fr

where ω is the wave frequency, N is the Brunt-Väisälä frequency of the fluid, assumed to be constant, β is the angle that the group velocity makes with the horizontal and f is the Coriolis parameter (f.i. [START_REF] Gill | Atmosphere-Ocean dynamics[END_REF]. For an incident plane wave propagating against the slope, conservation of frequency upon reflexion implies that the reflected wave is focused whatever the slope angle α of the topography at the reflexion location. Conservation of along-slope wavenumber and energy at reflexion imply that the wavenumber and amplitude of the reflected wave are larger than those of the incident wave, by a factor sin(β + α)/ sin |β -α|. If β > α, the reflected wave propagates (away from and) against the slope and reflection is said sub-critical (see figure 1); for β < α, the reflected wave propagates along the slope and the reflection is super-critical. For β close to α, referred to as critical incidence, focusing leads to strongly nonlinear processes close to the boundary (f.i. [START_REF] Mcphee-Shaw | Boundary layer intrusions from a sloping bottom: a mechanism for generating intermediate nepheloid layers[END_REF][START_REF] Chalamalla | Turbulence during the reflection of internal gravity waves at critical and near-critical slopes[END_REF]) and organized vortex structures sometimes called bores are then observed (f.i. [START_REF] Hosegood | Near-bed solibores over the continental slope in the faroe-shetland channel[END_REF].

Away from critical incidence, the interaction between the incident and the reflected waves leads to the generation of higher harmonic motions of frequencies nω such that 2 n < ω/N and amplitude A n , where A is the incident wave amplitude, with a much weaker impact on the fluid. However, as [START_REF] Thorpe | On the reflection of a train of finite-amplitude internal waves from a uniform slope[END_REF] (hereafter referred to as TH87) showed it theoretically, the incident and the reflected plane waves can form a resonant triad with a second harmonic wave resulting from their interaction, opening the possibility of breaking and mixing away from the slope. This result is valid only when the slope is inclined (α = 0) as no harmonic motion is produced when the incident plane wave reflects on a flat surface [START_REF] Thorpe | On the shape of progressive internal waves[END_REF]). In [START_REF] Thorpe | On the interactions of internal waves reflecting from slopes[END_REF], the computation of the resonance condition is extended to a rotating fluid.

However, in natural flows, internal gravity waves do not propagate as pure plane waves. This is the case for the internal tide, which propagates as beams from a submarine mountain or a continental shelf. (The internal tide is the internal wave field resulting from the interaction of the barotropic tide with submarine topography.) This is also the case, to some extent, for lee waves, which propagate in the lee of a mountain and most often, above it. In laboratory experiments as well, the generation of internal gravity waves always leads to wave packets of finite widths, for practical reasons. Even the wave generator originally designed by [START_REF] Gostiaux | A novel internal waves generator[END_REF] to model pure plane waves produces a wave packet containing a finite number of wavelengths in the direction normal to the propagation of the packet. For a large number of wavelengths, say four and above, the incident wave structure may be assumed to be close to that of a pure plane wave, at least away from the boundaries of the packet; when the packet contains one wavelength, the wave structure is closer to that of a beam. The wave packet produced by this wave generator will be referred to as a finite-width plane wave hereafter.

When a beam incident on a surface interacts with the reflected beam, harmonic motions are produced, whether the surface is flat or not [START_REF] Tabaei | Nonlinear effects in reflecting and colliding internal wave beams[END_REF]. A single beam is indeed a solution of the inviscid Boussinesq equations [START_REF] Tabaei | Nonlinear internal gravity wave beams[END_REF] while the superposition of two beams is not. This accounts for harmonics to be generated close to the boundary, where the beams interact, and this result holds for finite-width plane waves. The generation of higher harmonics by beams (or by finite-width plane waves) reflecting either on a flat surface or on an inclined boundary has been addressed in several numerical and experimental works performed in a two-dimensional (or quasitwo-dimensional) vertical geometry. [START_REF] Gostiaux | Quantitative laboratory observations of internal wave reflection on ascending slopes[END_REF] studied from laboratory experiments the structure of the higher harmonics resulting from the reflexion of a finitewidth plane wave on a flat surface. In [START_REF] Pairaud | Generation of harmonics and sub-harmonics from an internal tide in a uniformly stratified fluid: numerical and laboratory experiments[END_REF], the structure and evolution of an internal tide beam emitted from a two-dimensional topography is analysed from laboratory experiments and numerical simulations. It is shown that the reflection of the wave beam on the flat boundary at the foot of the slope leads to the generation of harmonic beams, consisting of free and trapped waves. [START_REF] Rodenborn | Harmonic generation by reflecting internal waves[END_REF] considered a beam reflecting on a simple slope and investigated empirically the conditions that lead to the largest amplitude of the second harmonic wave.

The purpose of the present paper is to examine the validity of the theory by TH87 when the conditions of a resonant triad involving an incident and a reflected wave of finite width and their second-order harmonic wave are satisfied. In a first part, we revisit TH87 theory for the case of a pure plane wave reflecting on a uniform slope. It is indeed unclear whether this theory predicts a finite amplitude of the second harmonic wave at resonance because the expression of this amplitude involves an indeterminacy (the numerator and denominator both vanish), which leads to the common inference in the literature that this amplitude is unbounded at resonance. In the second part, predictions of the theory of TH87 are determined for the case of an incident plane wave of finite width. The theoretical predictions are compared to the results of numerical simulations in a vertical plane, consistent with the resonant triad geometry, for parameters of laboratory experiments.

The plan of the paper is as follows. In section 2, the inviscid theory of TH87 for a pure plane wave incident onto a uniform slope is presented and extended; we show that, at resonance, the amplitude of the second harmonic plane wave actually grows linearly with distance from the slope. The predictions of the theory when the incident plane wave is of finite width are presented in section 3. We performed numerical simulations to estimate the validity of these predictions and the results are discussed in section 4. In section 5, the theoretical model is used again, in which viscous and diffusive effects are introduced in order to compare the transfer of energy from the incident and reflected waves to the second harmonic wave with the energy dissipated by molecular effects. A summary and conclusions are given in section 6.

Nonlinear reflection of a pure plane wave

Governing equations

The governing equations are the Navier-Stokes equations in the Boussinesq approximation for a rotating fluid with uniform stratification:

∂ ∂t u + (u • ∇)u = - 1 ρ 0 ∇p -f e z × u + b e z , (2.1a) ∂ ∂t b + u • ∇b + N 2 w = 0 , (2.1b) ∇ • u = 0 , (2.1c) u • n| slope = 0 .
(2.1d ) u = (u, v, w) is the velocity vector, w being the vertical velocity component, p is the non-hydrostatic pressure and b is the buoyancy; ρ 0 is a reference density and e z and n are unit vectors along the vertical axis pointing upward and normal to the slope, respectively. The parameters f and N have been defined in the Introduction. The last equation expresses the impermeability condition at the boundary.

We shall assume that wave propagation occurs in the vertical (x, z) plane and that the problem is invariant in the y-direction. Therefore all three velocity components depend on x and z (and on time) only and the incompressibility condition (2.1c) reduces to ∂u/∂x + ∂w/∂z = 0. As usual, this condition is automatically satisfied if a stream function ψ is introduced. Expressing the velocity vector as u = (∂ψ/∂z, v, -∂ψ/∂x) and combining the equations for u and w into a single equation for the vorticity ∂w/∂x -∂u/∂z = -∇ 2 ψ, equations (2.1) become

∂ ∂t ∇ 2 ψ + ∂ ∂x b -f ∂ ∂z v = J(ψ, ∇ 2 ψ) , (2.2a) ∂ ∂t v + f ∂ ∂z ψ = J(ψ, v) , (2.2b) ∂ ∂t b -N 2 ∂ ∂x ψ = J(ψ, b) , (2.2c) u • n| slope = 0 , (2.2d )
where J designates the Jacobian operator. Equation (2.2a) is rewritten by taking the time derivative and by replacing the fields ∂v/∂t and ∂b/∂t by their expression from equations (2.2b) and (2.2c), respectively. This yields

∂ 2 ∂t 2 ∇ 2 ψ + N 2 ∂ 2 ∂x 2 ψ + f 2 ∂ 2 ∂z 2 ψ = ∂ ∂t J(ψ, ∇ 2 ψ) - ∂ ∂x J(ψ, b) + f ∂ ∂z J(ψ, v) . (2.
3)

The set of equations (2.2b), (2.2c), (2.2d) and (2.3) will be referred to as equations (E) hereafter.

Let us rewrite this set of equations in a rotated coordinate system (x , z ) with x and z being the along-and normal-slope coordinates respectively. Following TH87, we use the convenient abbreviations c α , s α , c β and s β , standing for cos α, sin α, cos β and sin β, respectively. Thus (x , z ) = (xc α + zs α , zc α -xs α ), α being the slope angle. In this coordinate system, for a pure plane wave with frequency ω and along-and normal-slope wavenumbers k and nk, respectively, the dispersion relation (1.1) becomes

ω 2 = N 2 s 2 β + f 2 c 2 β = N 2 (ns α -c α ) 2 1 + n 2 + f 2 (s α + nc α ) 2 1 + n 2 . (2.4)
Solving this equation for n yields two roots

n I = s β c β -s α c α s 2 β -s 2 α , and n R = - s β c β + s α c α s 2 β -s 2 α , (2.5)
which are associated with the incident and reflected wave, respectively. Note that these roots are no longer defined as in TH87 as we have chosen to scale the normal-slope wave number by the along-slope wave number k. This turns out being handier for some further calculations as k is conserved through the reflection process.

Equations (E) are made nondimensional using the amplitude of the incident wave velocity, denoted U , as a velocity scale, its wavelength λ as a length scale and 1/N as a time scale. These scales come into play in the nondimensional equations of motions through a Froude number Fr = U/λN and through the parameter f /N . This Froude number is related to the wave steepness s (defined such that s = 1 when overturning occurs) by s = 2πFr / tan β.

(2.6)

Both F r and s can be interpreted as a non dimensional amplitude of the incident wave field.

When expressed in the rotated coordinate system and made nondimensional, equations (E) become, keeping the same notations for the nondimensional variables as well as for the nondimensional parameters N and f , to keep track of the Brunt-Väisälä frequency (otherwise equal to 1) in the equations:

∂ 2 ∂t 2 ∇ 2 ψ + N 2 c α ∂ ∂x -s α ∂ ∂z 2 ψ + f 2 s α ∂ ∂x + c α ∂ ∂z 2 ψ = Fr ∂ ∂t J(ψ, ∇ 2 ψ) -c α ∂ ∂x -s α ∂ ∂z J(ψ, b) + f s α ∂ ∂x + c α ∂ ∂z J(ψ, v) , (2.7a) ∂ ∂t v + f s α ∂ ∂x + c α ∂ ∂z ψ = Fr J(ψ, v) , (2.7b) ∂ ∂t b -N 2 c α ∂ ∂x -s α ∂ ∂z ψ = Fr J(ψ, b) , (2.7c) ∂ ψ ∂x z =0 = 0 .
(2.7d )

In the following, Fr is considered as a small parameter (i.e. F r 1), namely the wave dynamics are weakly nonlinear. The ψ, b and v fields can thus be decomposed into a multiple scale expansion as

ψ = ψ 1 + Fr ψ 2 + o(Fr 2 ), b = b 1 + Fr b 2 + o(Fr 2 ) and v = v 1 + Fr v 2 + o(Fr 2
). Introducing these expansions in the equations of motion (2.7), the resulting equations can be solved order by order, by matching terms multiplied by the same power of F r. At each order the solution satisfies the same linear operator with a right hand side determined by lower order nonlinear terms. Identifying terms multiplied by Fr 0 , the first-order fields ψ 1 , v 1 and b 1 are thus found to satisfy the homogeneous linear wave equations. Matching terms multiplied by Fr 1 shows that the second-order fields ψ 2 , v 2 and b 2 are solutions of the same linear operator forced by nonlinear terms solely involving the first-order fields. The detailed calculations associated with this general principle are now presented.

First order solution

The linear equation satisfied by ψ 1 is given by

∂ 2 ∂t 2 ∇ 2 ψ 1 + N 2 c α ∂ ∂x -s α ∂ ∂z 2 ψ 1 + f 2 s α ∂ ∂x + c α ∂ ∂z 2 ψ 1 = 0 (2.8a)
with the boundary condition

∂ ψ 1 ∂x z =0 = 0. (2.8b)
Looking for a pure plane wave solution with wavenumbers (k, nk) and frequency ω, the wave parameters are found to satisfy the dispersion relation (2.4), implying that the solution ψ 1 is a linear superposition of waves with wavenumbers (k, n R k) and (k, n I k),

with n I and n R given by (2.5). Imposing the boundary condition (2.8b) yields [START_REF] Phillips | Dynamics of the upper ocean[END_REF]:

ψ 1 = a 1 [sin (kx + n I kz -ωt) -sin (kx + n R kz -ωt)] .
(2.9)

The first-order stream function is therefore the sum of a wave incident onto the slope, which we denote for simplicity as (k, n I k, ω), and of its reflected counterpart (k, n R k, ω).

The wave amplitude a 1 is arbitrary. The stream function ψ 1 is also defined up to a constant additional phase, chosen to be 0 in the present case, which does not imply any loss of generality.

Once ψ 1 is known, v 1 and b 1 are inferred from the equations

∂ ∂t v 1 + f s α ∂ ∂x + c α ∂ ∂z ψ 1 = 0 (2.10) and ∂ ∂t b 1 -N 2 c α ∂ ∂x -s α ∂ ∂z ψ 1 = 0 (2.11) which yields v 1 = a 1 f k ω [(s α + c α n I ) sin (kx + n I kz -ωt) -(s α + c α n R ) sin (kx + n R kz -ωt)]
(2.12) and

b 1 = a 1 N 2 k ω [(c α -s α n R ) sin (kx + n R kz -ωt) -(c α -s α n I ) sin (kx + n I kz -ωt)] .
(2.13)

Second order solution

Matching terms proportional to Fr in (2.7) gives the equations for the second-order component. The equation for ψ 2 is

∂ 2 ∂t 2 ∇ 2 ψ 2 + N 2 c α ∂ ∂x -s α ∂ ∂z 2 ψ 2 + f 2 s α ∂ ∂x + c α ∂ ∂z 2 ψ 2 = ∂ ∂t J(ψ 1 , ∇ 2 ψ 1 ) -c α ∂ ∂x -s α ∂ ∂z J(ψ 1 , b 1 ) + f s α ∂ ∂x + c α ∂ ∂z J(ψ 1 , v 1 )
(2.14a) with the boundary condition

∂ ψ 2 ∂x z =0 = 0. (2.14b)
The solution of equations (2.14), which generalizes TH87 solution to a rotating fluid, is:

ψ 2 =3a 2 1 k 2 ω s 2 β c 2 β s α c α (s 2 β -s 2 α )D sin 2kx + m 2 kz -2ωt -sin 2kx + (n I + n R )kz -2ωt - a 2 1 k 2 ω s β c β s 2 β -s 2 α sin (n I -n R )kz (2.15) with D = N 2 s 2 β 4s 4 β -7s 2 β s 2 α + 4s 2 α -s 2 β + f 2 c 2 β 4s 4 β -7s 2 β s 2 α + 3s 2 α (2.16)
and

m 2 = -2s α c α (N 2 -f 2 ) + 2 4s 2 β N 4 (1 -4s 2 β ) + 4c 2 β f 4 (1 -4c 2 β ) + N 2 f 2 (3 -32s 2 β c 2 β ) N 2 (4s 2 β -s 2 α ) + f 2 (4c 2 β -c 2 α )
.

(2.17) Details of this calculation can be found in appendix A.1 where the full second order solution is derived.

The second-order stream function (2.15) is the sum of an oscillatory component (for α = 0) and a steady component. The oscillatory component, referred to as the second harmonic wave in the introduction, is made of two terms, associated with the waves (2k, m 2 k, 2ω) and (2k, (n I + n R )k, 2ω). The former term is the solution of the homogeneous equation associated with equation (2.14a) so that the wave (2k, m 2 k, 2ω) satisfies the dispersion relation and radiates energy away from the boundary. This term is called a "free wave" by TH87. As detailed in appendix A.1, (2.15) is obtained by adding to this homogeneous solution a particular solution of the complete equation, referred to as a "forced" wave by TH87, and by prescribing the homogeneous solution parameters through the boundary condition (2.14b). The forced wave becomes a free wave when (2k, (n I + n R )k, 2ω) satisfies the dispersion relation, namely when n I + n R = m 2 ; in this case, the incident, reflected and second harmonic waves form a resonant triad. A sketch of the interaction between the incident and reflected waves for conditions away from resonance is displayed in figure 1.

The steady component of (2.15) is an Eulerian mean current opposing the Stokes

x[m] x' z[m] z' Cg i k i Cg r k r Cg 2 k 2 0 0.5 1 1.5 2 2.5 3 -0.8 -0.6 -0.4 -0.2 0 Figure 1.
Sketch of the interaction between the incident and reflected waves in the (x, z) plane.

The incident wave propagates in a vertical plane and has a finite width in the direction normal to propagation. This finite width contains four wavelengths in the sketch. The wave vector and group velocity of the incident wave are denoted k i = (k, nI k) and Cg i . Those for the reflected wave are kr = (k, nRk) and Cgr. The blue and green solid lines represent phase lines of the incident and the reflected waves, respectively, and the interaction area between these waves is highlighted by a yellow triangle. The red solid lines in the interaction area represent phase lines of the second harmonic "forced wave" (2k, (nI + nR)k, 2ω) and red dashed lines correspond to the second harmonic "free wave" (2k, m2 k, 2ω) able to propagate outside the interaction area (at resonance the red solid and dashed lines coincide, namely the forced wave becomes a free wave). The slope angle is α and β is the angle of incidence. The coordinates (x , z ) refer to the rotated frame of reference such that x lies along the slope. The domain dimensions, expressed in meters, are those used in the numerical simulations.

drift associated with the first-order solution ((2.9),(2.12),(2.13)), as stated by TH87 and [START_REF] Thorpe | On the interactions of internal waves reflecting from slopes[END_REF]. Hence the total Lagrangian mean flow, which is the sum of the Eulerian mean flow and of the Stokes drift (Longuet-Higgins 1969), is zero (see appendix A.3). Indeed, because of the sloping boundary, there cannot be any horizontal mass transport in the present two-dimensional configuration, implying that the Lagrangian mean flow must vanish.

We now consider the oscillatory component, denoted ψ h 2 . The denominator of its amplitude goes to 0 for values of (α, β) satisfying D = 0. For f = 0, these values span the intervals [0, 8.21 • ] for α and [0, 30 • ] for β, these ranges shrinking to empty intervals as f /N increases and reaches 0.5. As stated in TH87 and detailed in appendix A.2, the condition D = 0 also corresponds to the case of a resonant triad between the primary incident and reflected waves and the second harmonic free wave. Since m 2 = n I + n R in this case, as just discussed, the oscillatory component in (2.15) is an indeterminate form 0/0 as D → 0. One aim of the present paper is to remove this indeterminacy to clarify the behavior of the oscillatory component at resonance. For this purpose, we rewrite the oscillatory component of ψ 2 as

ψ h 2 = 6a 2 1 k 2 ω s 2 β c 2 β s α c α s 2 β -s 2 α sin(δkz ) D cos 2kx + (m 2 -δ)kz -2ωt (2.18) where δ = 1 2 (m 2 -(n I + n R )). For α = 0, one can show that δ ∼ D→0 D 6s α c α s 2 β -s 2 α ω 2 , (2.19) so that sin(δkz ) D → D→0 kz 6s α c α s 2 β -s 2 α ω 2
(2.20) (see appendix A.2 for details). Hence, as D → 0, the second harmonic wave becomes

ψ h 2 = a 2 1 k 2 ω s 2 β c 2 β (s 2 β -s 2 α ) 2 kz cos 2kx + (n I + n R )kz -2ωt .
(2.21)

The amplitude of the second harmonic wave is thus a periodic function of the distance from the slope z whose period 2π/δk goes to infinity as (α, β) approaches the resonance condition D = 0. In this case, the amplitude grows linearly from the slope and is therefore unbounded. This linear growth is consistent with the behavior of a resonantly forced harmonic oscillator, whose amplitude grows linearly with time. This unbounded amplitude is not physical and is due to the weakly nonlinear character of the theory: the first-order solution is supposed to be not modified by the transfer of energy to the second-order solution. Hence, at resonance, the second harmonic wave is continuously fed during its propagation by the unaltered primary wave.

Nonlinear reflection of a plane wave of finite width

The results presented in section 2.3 have been obtained for a pure plane wave. As discussed in the Introduction, this type of wave is commonly modelled in laboratory experiments by a wave generator, which creates in practice a wave packet containing a finite number of wavelengths in the direction normal to the propagation of the packet. The interaction area between the incident and the reflected waves is now spatially bounded and is represented by the yellow triangle in figure 1. In the present section, we compute the amplitude of the second-order wave amplitude when resonant conditions are met for this generation method. The incident wave field is modelled theoretically as a plane wave with a finite number of wavelengths in the direction normal to wave propagation, with no envelop. In section 4, we compare the theoretical predictions with numerical simulations for an incident wave with various widths.

We recall that, to be valid, the theory requires the amplitude of F r ψ 2 to stay small with respect to the amplitude of ψ 1 . Since we focus on the oscillatory component of ψ 2 , namely ψ h 2 , we scale the amplitude of F r ψ h 2 with that of ψ 1 . More precisely, we scale Fr ψ h 2 by ψ 1 where norm . is the infinity norm in space and time. This norm is defined by ψ = max

(x,t)∈Ω×R +
|ψ|, where Ω is the interaction area between the incident and reflected waves. We thus compute the ratio

R th = Fr ψ h 2 ψ 1 , (3.1)
the th superscript standing for theoretical. At resonance, the amplitude of ψ h 2 should reach a maximum value at the largest distance from the slope in the interaction area, which coincides with the top of the yellow triangle. Our purpose here is (i) to compute R th as a function of the parameters of the incident wave and (ii) to investigate the conditions under which the theory is strictly valid (namely R th 1).

Let us compute ψ 1 and ψ h 2 . If one rewrites ψ 1 in the same way as ψ h 2 , i.e.

ψ 1 = 2a 1 sin n I -n R 2 kz cos kx + n I + n R 2 kz -ωt , (3.2)
the expressions of ψ 1 and ψ h 2 are given by

ψ 1 = 2a 1 max 0 kz kh sin n I -n R 2 kz (3.3a) and ψ h 2 = 6a 2 1 k 2 ω s 2 β c 2 β s α c α s 2 β -s 2 α max 0 kz kh sin(δkz ) D , (3.3b)
where h is the height of the interaction triangle normal to the slope. Introducing n λ , the number of wavelengths contained in the primary wave, and thanks to simple trigonometric calculations, one can show that

kh = 2π sin(β -α) sin(β + α) sin(2β) n λ = 2π (s 2 β -s 2 α ) sin(2β) n λ = 2π n I -n R n λ . (3.4)
This result implies that 0.5 (n I -n R ) kz πn λ for 0 kz kh. Since n λ 1, it ensures that max 0 kz kh {sin (0.5 (n I -n R ) kz )} = 1 and leads to ψ 1 = 2a 1 . Note that the first equality of (3.4) is in agreement with [START_REF] Thorpe | On the reflection of internal wave groups from sloping topography[END_REF].

The expression of ψ h

2 cannot be simplified in the same way, especially when considering the situation of (α, β) close to the resonant triad case where δ → 0. For simplicity, we introduce

M = max 0 kz kh sin(δkz ) D =        kh/(6s α c α (s 2 β -s 2 α ) ω 2 ) if D = 0 1/D if D = 0 and δkh > π/2
sin(δkh)/D if D = 0 and δkh π/2, (3.5) assuming α = 0. A measure of the ratio between the second harmonic wave and the first-order solution is thus given by

R th = 3Fr a 1 k 2 ω s 2 β c 2 β s α c α s 2 β -s 2 α M. (3.6)
For the sake of completeness, we provide the expressions of the nondimensional quantities a 1 , k and ω with the scaling defined in section 2.1:

a 1 = 1/2π, k = 2π sin(α + β), ω = s 2 β + (f /N ) 2 c 2 β , (3.7)
while N and f become equal to 1 and f /N , respectively. In the following, the parameter f /N is denoted as γ. With this scaling, the expression of R th becomes:

R th = 6πFr s 2 β + γ 2 c 2 β sin (β + α) s 2 β c 2 β s α c α sin (β -α) M. (3.8)
R th is proportional to Fr and, at resonance (D = 0), to the number of wavelengths n λ as well (through M and equation (3.4)). This implies that, as n λ → ∞, this ratio may well become much larger than 1 since F r and n λ are independent parameters. Hence the internal wave reflection problem is highly nonlinear at resonance, which figures 2 and 3 below attest.

Figure 2 displays R th in a (α, β) diagram for γ = 0 (no rotation), F r = 0.005 and for four different values of n λ equal to 1, 2, 4 and 8. These values, except for n λ = 8, are those of the numerical simulations reported in section 4. The counterpart diagram for γ = 0.2 is displayed in figure 3. Both figures show that, even if the Froude number is quite low and n λ is at most equal to 8, R th reaches in all cases values close to unity at resonance. Figure 2 also shows that, for small values of n λ , the location of resonant triads (marked with a dashed curve) differs from the location where R th reaches a maximum value, both locations superimposing as n λ grows. In other words, we recover TH87 prediction in the limit of pure plane waves. The latter statement is confirmed in figure 8, further discussed in section 4, where R th (blue curve) is plotted versus α for n λ = 1, 2 and 4 and for β = 18.9 • : for n λ 2, R th displays a maximum which is all the more pronounced n λ is larger; this maximum is reached for a value of α that approaches the theoretical value, equal to 7.848 • , at which resonance occurs when n λ increases. Thus, the second harmonic resonance in the sense of a pronounced local maximum for R th in the (α, β) plane only arises for high values of n λ which are not likely to occur in reality.

Figure 3 shows that rotation does not qualitatively modify the resonance process, only making this maximum lower than in the non-rotating case. This behavior can be explained with equation (2.21): the influence of rotation effects in the second-order stream function amplitude at resonance only comes into play through the frequency ω in the denominator of the expression of ψ h 2 . All parameters being kept the same, adding rotation increases ω and therefore lowers the second-order amplitude.

Comparison with numerical simulations

In order to estimate the validity of the theoretical predictions exposed in the previous section, two-dimensional numerical simulations have been performed. As just discussed, including rotation does not qualitatively change the results so that the non-rotating case is considered in the present section. All quantities below are dimensional but we keep the same notation as before for the variables and parameters, for simplicity. In the following, the term first harmonic wave refers to the superposition of the incident and reflected waves.

Figure 2. Second-order to first-order stream function norms R th , defined by (3.8), as a function of the slope angle α and incidence angle β, for increasing width of the incident wave defined by the number of wavelengths n λ it contains. For all panels, the Froude number Fr = 0.005 and f = 0 (no rotation). The dashed curve corresponds to the resonance condition D = 0 for a pure plane wave (n λ → +∞), with D defined by equation (2.16). The crosses refer to numerical simulations reported in section 4. The domain is limited to β 30 • because the second harmonic free wave becomes evanescent above this range. Values of R th are also not plotted when approaching the critical case α = β which is beyond the scope of this article.

Numerical configuration

The simulations were performed with the numerical model NHOES (Non Hydrostatic Ocean model for Earth Simulator) which solves the free-surface non-hydrostatic Boussinesq equations in a Cartesian coordinate system [START_REF] Aiki | A numerical study on the successive formation of meddy like lenses[END_REF]. For the present study, equations (2.1) were solved with a viscous term ν∇ 2 u added to the right-hand-side of equation (2.1a) and a diffusive term κ∇ 2 b to that of equation (2.1b); ν and κ are the kinematic viscosity and the diffusivity of buoyancy, respectively. The equation of state is linear and only depends on salinity. The choice of salinity rather than temperature is dictated by laboratory experiments conducted in parallel to the present study on the Coriolis platform at Grenoble in which density variations are created by a vertical profile of salt concentration. The viscosity is thus equal to ν = 10 -6 m 2 s -1 , namely that of water, while the diffusivity is set to 1.49 10 -9 m 2 s -1 since the Prandtl number of salt dissolved in water is about 700. No sub-grid scale parametrisation is used. The background stratification is linear, with a value of the Brunt-Väisälä frequency equal to 0.46 s -1 and, as said above, the Coriolis frequency f is set to 0. These values and all values of the physical and geometrical parameters below are those of the laboratory experiments.

The numerical setup is sketched in figure 1. The domain dimensions are 3 m in the horizontal direction and 0.8 m in the vertical direction. The bottom boundary condition is of the free-slip type and an implicit free-surface boundary condition is imposed at the top of the domain. The origin of the vertical axis, oriented positively upwards, is taken at the free surface.

An internal gravity wave of finite vertical width is continuously forced at the left boundary of the numerical domain (x = 0) by adding a forcing function to the righthand-side of the momentum equation for the u-component (see figure 4):

F x (0, z, t) = U E(z) (1 -e -5t T ) ω sin(k z z -ωt) + 5t T e -5t T cos(k z z -ωt) (4.1a) with E(z) = 1 2 1 + tanh 3 2 tan - π 2 + 1 2 min(2π, max(0, π(n λ + 1) -k z |z -z 0 |))
(4.1b) where U is the velocity amplitude of the forcing, ω the wave frequency (and T the wave period) and k z is the vertical wave number. This forcing generates the u-component of the incident plane wave. The smooth envelope in the z-direction E(z) is centered at z 0 = -0.42 m. The value of the frequency is equal to 0.149 rad s -1 implying that the propagation angle of the incident wave β = 18.9 • . The value of the vertical wavelength λ z is 0.125 m implying that the wavelength λ is equal to 0.118 m; four wavelengths are therefore permitted at most in the incident wave. The slope angle α varies from 0 • to 12 • . The value of α at which resonance theoretically occurs in this case, referred to as α res , is equal to 7.848 • . The horizontal velocity amplitude is set to U = 3 10 -4 m s -1 . With these values, the Froude number of the incident wave Fr = 5 10 -3 1 which satisfies the condition for resonant interaction to occur. The corresponding value of the incident wave steepness is s = 0.09, using equation (2.6).

A sponge layer of width 0.5 m is added at the right boundary in order to prevent wave reflection. The horizontal and vertical resolutions are equal to 1.25 10 -3 m and 0.625 10 -3 m, respectively. The latter value has been chosen so that the height h of the interaction triangle between the incident and reflected waves, which decreases as α increases, contains about 20 grid points for the largest value of α we consider (equal to 12 • ). This vertical resolution implies that the vertical wavelength λ z is described by 200 grid points.

Results of the numerical simulations

As indicated by the crosses in figure 2, numerical simulations have been performed for three different widths of the incident wave (n λ = 1, 2 and 4) and 7 slope angles, equal to 0 • , 2 • , 4 • , 6 • , α res = 7.848 • , 10 • and 12 • . Each simulation has been carried out over 40 wave periods.

The off-slope velocity (normal to the slope) w = -∂ψ/∂x , filtered at either the incident wave frequency, denoted w 1 , or twice this frequency, denoted w 2 , is displayed for n λ = 2 and n λ = 4 in figures 5 and 6, respectively. The harmonic filtering has been performed over the last 8 periods of the simulations. The choice of the w variable is dictated by the forthcoming comparison with the theoretical predictions presented in the previous section.

First harmonic wave

We first consider the incident wave component w 1 displayed in figures 5a and 6a. Note that w 1 is proportional to k by definition, namely to | k| sin(α + β). Before entering the interaction area, the amplitude of w 1 grows with α, since | k| and β are kept fixed in all simulations. Figures 5b and6b show that very weak amplitude second-order harmonic waves are radiated from the forcing boundary of the incident wave due to the modulation of this wave by the envelope function (4.1b).

The incident and reflected waves superpose in the interaction area and the resulting stream function is given by equation (3.2). This equation implies that the ω-filtered interaction pattern does not vary along x and has a sinusoidal dependency in z , as confirmed by figures 5a and 6a, with wavenumber k(n I -n R )/2. Figures 5a and 6a however display a major feature, which is not accounted for in the theory. The theory assumes indeed that the off-slope velocity amplitude of the reflected wave is uniform and equal to that of the incident wave to satisfy the impermeability condition at the sloping boundary. While this is the case in the interaction area, at least close to the sloping boundary, the figures actually show that this amplitude is much weaker than that of the incident wave outside the interaction area, all the more so the slope angle is larger. This behavior has two origins. The main one arises from the focusing of the reflected wave and is due to molecular effects. The wavelength of the reflected wave, equal to λ sin(β -α)/ sin(β + α), is indeed smaller than the incident wavelength λ so that the viscous (or diffusive) time scale is smaller for the reflected wave as well. When scaled by the incident wave period, the viscous time scale decays from 16.8 to 0.9 as α increases from 0 • to 12 • . Dissipation thus becomes an important effect in the dynamics of the reflected wave as α increases. The second effect is due to energy transfer to higher harmonic waves. Figures 5b and6b show that, as α increases, the amplitude of the second harmonic wave becomes no longer small compared to the amplitude of the incident wave. There is therefore a significant energy flux from the first harmonic wave to the second harmonic one which also reduces the amplitude of the reflected wave. Theoretical estimates of the power of the first-harmonic wave lost by dissipation and by energy transfer to the second-order harmonic are provided in section 5.

Second harmonic wave

As discussed in the Introduction, the interaction of a finite-width wave incident of a 6b show indeed that a (weak amplitude) second harmonic wave propagates from the interaction area for α = 0. In the present case, only harmonics 2 and 3 can be generated (with 3 ω = 0.447 very close to N ) so that energy is dominantly transferred to the second harmonic wave.

405

When α is non zero, the second-order wave field ψ h 2 is the sum of a forced and a free wave, as discussed in section 2.3. The forced wave exists only in the interaction area between the incident and reflected waves. While being also generated inside the interaction area, the free wave can radiate away from this area. At resonance, the forced and free waves coincide (namely m 2 = n I + n R , see equation (2.15)) and the amplitude of ψ h 2 grows linearly with the normal to the slope z , the maximum value being reached at the top of the interaction area, for z = h.

Figures 5b and 6b

show that, for α = 0, the second-order free wave is clearly visible outside the interaction area, with amplitude equal to that at the boundary of this area. This amplitude is largest for α = α res , as expected. It is damped by dissipative effects along its travel to the upper boundary where it reflects. Focusing now on the interaction area, figures 5b and 6b show that the amplitude of w 2 does not vary along the slope, as predicted by equation (2.18). This is visible for α = α res and α = 12 • , and not so clearly for α = 4 • , consistent with figure 2 showing that the second-order amplitude is much weaker for α = 4 • than for the former angles. The behavior of w 2 normal to the slope at resonance is analysed in the next section.

Comparison with the theoretical predictions

According to equation (2.21), in the absence of fluid viscosity, the amplitude of ψ h 2 evolves linearly with z at resonance. To assess the validity of this theoretical prediction, this linear law is compared to the behavior normal to the slope of the amplitude of w 2 inside the interaction area (0 z h) for the simulation with α = α res and n λ = 4 (see figure 7a). A poor agreement is obtained. Since the height h is fairly well resolved in the numerical simulation (by more than 100 grid points), this discrepancy may be due to the inviscid character of the theory. The numerical simulation was therefore rerun with a ten times lower viscosity, namely ν = 10 -7 m 2 s -1 . The comparison is displayed in figure 7b: the agreement is much better, which confirms the prediction that the second-order harmonic wave at resonance evolves linearly with z inside the interaction area in the inviscid limit.

The validity of the theory can be further assessed by comparing the ratio R th defined by (3.1) with its numerical counterpart, denoted R num . Rewriting R th as

R th = Fr 2 2k ψ h 2 k ψ 1 = Fr 2 ∂ x ψ h 2 ∂ x ψ 1 , (4.2)
R num can be defined as, in dimensional form

R num = 1 2 w 2 w 1 . (4.3)
The norm of the first and second harmonic off-slope velocities w 1 and w 2 are defined by the maximum of their respective amplitude over the interaction area. R th is displayed in figure 8 for n λ = 1, 2 and 4. As discussed in section 3, R th displays a local maximum for n λ 2, which grows with n λ (while being upper-bounded) and is reached for a value of α approaching α res as n λ increases.

R num is compared to R th in figure 8 (star symbols). For α = 0, as expected, their values do not match since no harmonic wave is produced by the interaction of the plane incident and reflected waves (so that R th = 0) while a harmonic wave is generated when these waves are of finite thickess (implying that R num = 0). When α has a non-zero value, figure 8 shows that the simulations are in good agreement with the theory for the lowest angles only, up to α = α res for n λ = 1 and up to α = 4 • for n λ = 2 and n λ = 4. This can be explained as follows. The largest theoretical amplitude of the first harmonic wave is always reached in the simulations, but in the bottom part of the interaction area where the reflected wave has not been damped yet. The agreement with the theory (namely with R th ) therefore requires the theoretical maximum for the second harmonic wave to be reached in the simulation. This theoretical maximum occurs at the top of the interaction area. In the simulation, w 2 will also be reached at that location if the amplitude of the reflected wave has hardly decayed in the interaction area, namely if the focusing effect on that wave is moderate enough. This occurs for small slope angles.

For larger slope angles, the amplitude of the reflected wave has already weakened when reaching the top of the interaction area, leading to a weaker amplitude of the second harmonic wave at that location than predicted.

The major role of viscosity in damping the reflected wave inside the interaction area is attested by comparing the results of the theoretical predictions with those of simulations with ν = 10 -7 m 2 s -1 , for n λ = 4 and α 6 • . Results are plotted again in figure 8 (triangle symbols). The agreement with the theory strongly improves, confirming that the dissipative damping of the focused reflected wave is responsible for the discrepancy with the theory.

Theoretical estimate of the energy budget

In order to estimate the relative importance of the damping processes of the first harmonic wave, we determine in this section the power lost in the interaction area Ω by that wave due either to energy transfer to the second harmonic wave or to dissipation and compare both powers to the incident energy flux. All results are obtained from the weakly nonlinear theory presented in section 2 in which dissipative effects are introduced, with primary and second-order wave solutions derived in section 3.

We start from the general governing equations (2.1a) and (2.1b), which we write in The Froude number of the incident wave is Fr = 5 10 -3 , the angle of incidence is β = 18.9 • and an incident wave with three different widths is considered: n λ = 1, 2 and 4. The solid blue line is the theoretical ratio R th defined by equation (3.1) (expression (3.8) is plotted here). Its numerical counterpart R num defined by equation ( 4.3) is represented with symbols: stars for simulations with ν = 10 -6 m 2 s -1 and triangles for ν = 10 -7 m 2 s -1 . The dashed line indicates the value of α for which resonance is predicted by the theory of TH87.

nondimensional form:

∂ ∂t u + ∇p + f e z × u -b e z = -Fr (u • ∇)u , (5.1a) ∂ ∂t b + N 2 w = -Fr u • ∇b . (5.1b)
At first order, equations (5.1) become

∂ ∂t u 1 + ∇p 1 + f e z × u 1 -b 1 e z = 0 , (5.2a)
∂ ∂t b 1 + N 2 w 1 = 0 .
(5.2b)

These equations are analogous to equations (2.8a), (2.10) and (2.11) of section 2.2, now written in terms of the full velocity and buoyancy fields.

The usual operation u 1 • (5.2a) + (b 1 /N 2 )(5.2b) yields the conservation equation for the total energy of the first harmonic wave:

∂ ∂t E 1 + ∇ • (p 1 u 1 ) = 0 , (5.3)
with

E 1 = 1 2 u 2 1 + b 2 1 /N 2 .
Following the same steps at second order yields the energy equation

∂ ∂t E 2 + ∇ • (p 2 u 2 ) = -Fr 2 ((u 1 • ∇)u 1 ) • u 2 + 1 N 2 (u 1 • ∇b 1 )b 2 .
(5.4)

The right-hand-side of the equation is the transfer term from the first harmonic wave to the second order solution. This second order solution contains the second harmonic wave plus an Eulerian mean flow. This mean flow does not radiate any wave so that, when the transfer term from the first harmonic to the second harmonic wave power density is seeked for, it does not come into play in this term. The latter term, denoted tp 12 , can therefore be written as

tp 12 = -Fr 2 ((u 1 • ∇)u 1 ) h • u h 2 + 1 N 2 (u 1 • ∇b 1 ) h b h 2 , (5.5)
where the superscript h stands for the harmonic part. Note that for consistency, the nonlinear term involving the first harmonic wave is also written with the superscript h (despite it contains only wave components). The transfer term of power involved in the energy budget of the second harmonic wave over the interaction area Ω is

T P 12 = Ω tp 12 dx dz , (5.6)
where . designates the average operator over one wave period. We refer the reader to appendix B for the expression of this term and detailed calculations.

This transfer term has been determined in the weakly non-linear framework, hence with the assumption that the first harmonic wave is negligibly influenced by the second harmonic generation. In the same way, one can determine an estimate of the power of the first harmonic wave lost by viscous dissipation (diffusive effects being very small) by adding a viscous term to the RHS of equation (5.2a) and assuming that the inviscid first-order solution still holds. In this framework, the first harmonic power density lost by dissipation is given by

dp 1 = F r Re ∇u 1 : ∇u 1 = F r Re |∇u 1 2 | + |∇v 1 2 | + |∇w 1 2 | , (5.7)
where Re = U λ/ν is the Reynolds number. As above, when the energy budget of the primary wave over Ω is considered, an average operator over a primary wave period and over Ω should be applied to dp 1 leading to

DP 1 = Ω dp 1 dx dz (5.8)
(see appendix B for the expression of this term and detailed calculations).

We finally need to compare the power of the primary wave lost into nonlinear transfer T P 12 and dissipation DP 1 to the incident energy flux

P I = ∂Ω p I u I • dn = ∂Ω E I c gI • dn (5.9)
where the I subscript indicates the incident part of the first harmonic wave, c g is the group velocity, ∂Ω designates the boundary of the interaction area and dn the outgoing unit normal vector on ∂Ω. The second equality in equation (5.9) holds only for a plane monochromatic wave. The reader is again referred to the appendix for its expression.

The scaled powers T P 12 /P I , DP 1 /P I and their sum are shown in figure 9 for two different values of fluid viscosity, ν = 10 -6 m 2 s -1 and ν = 10 -7 m 2 s -1 . It should first be noted that both T P 12 and DP 1 are overestimated as their expressions have been derived under the assumption that the first harmonic wave is not influenced by second harmonic transfer and dissipation (in particular, the lower amplitude of the reflected wave is not taken into account). The first harmonic wave therefore behaves as an infinite energy source for the transfer and dissipation terms so that the scaled power ratios can exceed 1. Figure 9 shows that, whatever the value of ν, the dissipated power of the primary wave increases with α because of the focusing of the reflected wave. For ν = 10 -6 m 2 s -1 (first row), the power lost by the first harmonic wave is dominated by dissipation whatever n λ , even when resonance occurs (at which the transfer term displays a maximum). By contrast, when the value of fluid viscosity is lowered to ν = 10 -7 m 2 s -1 (second row of figure 9), transfers to the second harmonic wave are much higher than power lost by dissipation, up to α 12 • or so. For a larger angle, the nonlinear transfer term decreases while the impact of focusing on the reflected primary wave increases, promoting dissipation.

These power estimates confirm the previous statement that dissipation of the first harmonic wave is the main process accounting for the discrepancy between R th and R num for ν = 10 -6 m 2 s -1 , due to the focusing of the reflected wave.

These results also imply that two-dimensional laboratory experiments (with water) cannot evaluate the validity of TH87 theory, for a finite-width incident wave for the value of the primary wave Froude number we consider (F r = 0.005). The reason is that the Reynolds number of the primary wave Re = U λ/ν, equal to 35 in the simulations, is too low. Decreasing ν by a factor 10 in the Reynolds number, as done in the simulations, is equivalent to increasing U by a factor 10, suggesting that favorable conditions to test TH87 theory could be achieved in the laboratory experiments. However, with F r increasing by a factor 10 as well, a strongly nonlinear regime is obtained since the incident wave steepness becomes equal to 0.92. The conditions to test Thorpe's theory for a finitewidth wave are therefore not met either. (Increasing the wavelength by a factor 10 would ensure both a higher Reynolds number and a lower Froude number but is not feasible in practice.)

Summary and conclusion

The purpose of this paper is to address the nonlinear reflection of a finite-width internal gravity wave incident onto a uniform slope, away from critical incidence, in a twodimensional vertical plane. The incident wave propagates against the slope, leading to focusing of the reflected wave. The paper has three objectives: (i) to revisit the inviscid and weakly nonlinear theory of TH87 when resonant conditions between an incident pure plane wave, the reflected wave and the second-order harmonic waves are met; at resonance indeed, the amplitude of the second-order stream function is expressed as an indeterminate form leading to the common inference in the literature that this amplitude diverges; (ii) to apply this theory to an incident wave of finite width in the direction normal to wave propagation, as produced by the wave generator device commonly used in laboratory experiments to model plane wave dynamics; (iii) to compare the latter theoretical predictions with results of two-dimensional numerical simulations of a wave of finite width.

We show that the indeterminacy at resonance can be waived and that the amplitude of the second-order stream function is a linear function of the distance from the slope. This amplitude is therefore unbounded for a pure plane wave in the weakly nonlinear limit.

For an incident wave of finite width, which we model theoretically as a pure plane wave . Energy budget for the first harmonic in the interaction area Ω (yellow triangle in figure 1) as obtained from the weakly nonlinear theory in which dissipative effects have been added. The molecular viscosity is equal to ν = 10 -6 m 2 s -1 in the first row and to ν = 10 -7 m 2 s -1 in the second row. The green and red curves represent the power averaged over one period transferred to the second harmonic and dissipated by viscous effects, respectively, as a function of the slope angle α; each power is scaled by the incident energy flux in Ω. The blue curve is the sum of the scaled powers. The fixed parameters are the incidence angle β = 18.9 • , the Froude number Fr = 5 10 -3 and the Reynolds number Re = 35 used in the numerical simulations.

with a finite number of wavelengths, the generation of second-order harmonic waves is limited to the area where the incident and reflected waves superpose. The amplitude of the second-order stream function is therefore now bounded. At resonance, its maximum value scaled by that of the first-order wave increases with the number of wavelengths n λ and becomes close to 1 for n λ 4 even for a Froude number of the incident wave as small as 0.005. This implies that the reflection problem at resonance is highly nonlinear. We also showed that rotation does not qualitatively change the results, the second-order amplitude being weaker when rotation is imposed. Numerical simulations of an incident wave of finite width were thus performed in the non-rotating case to estimate the validity of the theoretical predictions, for parameters of laboratory experiments performed in parallel to the present study. Various slope angles, still well below critical incidence, and incident wave widths were considered in the numerical simulations.

The numerical simulations show that the amplitude of the reflected wave is damped for a sloping boundary, all the more as the slope angle is larger. For the highest angle we consider (α = 12 • with β = 18.9 • ) almost no reflected wave leaves the interaction region. The reason lies in the focusing of that wave, which enhances viscous dissipation. As a consequence, the numerical simulations show a good agreement with the theoretical results for low slope angles only, well below critical incidence. This comparison is quantified by introducing dissipative effects in the weakly non linear theory. We show that, as the slope angle increases, the dissipated power of the primary wave becomes larger than the power transferred to the second harmonic wave, accounting for the discrepancy with the theory. The agreement strongly improves when the viscosity is divided by a factor 10, consistent with this discrepancy being due to the dissipative damping of the focused reflected wave. These results imply that the validity of the theoretical predictions cannot be assessed by (quasi-) two-dimensional laboratory experiments in a vertical plane, because of the too low value of the incident wave Reynolds number (equal to 35 in the numerical simulations). Decreasing the value of the viscosity by a factor 10 in the simulations is equivalent to increase the wave amplitude by the same factor in the laboratory experiments to keep the same Reynolds number. However the now ten times larger Froude number implies that the incident wave would no longer be of weak amplitude in the experiments (its steepness being equal to 0.95), so that resonant conditions cannot be met.

The three-dimensional laboratory experiments conducted in parallel to the present two-dimensional numerical simulations are not appropriate either to test the theoretical predictions. Indeed a Lagrangian mean flow is induced in the interaction area as a result of nonlinear and dissipative effects, whose amplitude can be as large as that of the incident primary wave due to dissipative cumulative effect [START_REF] Grisouard | Large scale energy transfer from an internal gravity wave reflecting on a simple slope[END_REF]. As a result, the incident wave frequency is shifted through Doppler effect by the mean flow and refracted (with the phase lines bending toward the horizontal), which deeply modify the incident primary wave geometry. This mean flow does not occur in a two-dimensional vertical plane because of the topography, which breaks the horizontal homogeneity required for the mean flow to develop. (This mean flow results from the same generation process as that induced by a single wave beam, as observed by [START_REF] Bordes | Experimental observation of a strong mean flow induced by internal gravity waves[END_REF] in their laboratory experiments, and theoretically modelled by [START_REF] Kataoka | On three-dimensional internal gravity wave beams and induced large-scale mean flows[END_REF].) Threedimensional numerical simulations should therefore be performed for any comparison to be performed with the laboratory experiments. This joint approach has been carried out [START_REF] Raja | Forcing of a mean flow during the reflection of a three dimensional internal wave beam[END_REF]) and the analysis of the resulting nonlinear wave-topography interaction will be reported in a subsequent paper.

Appendix A. Weakly nonlinear second-order calculations

A.1. Computing the second-order solution

In this section we solve the linear partial differential equation (2.14). Using (2.9), (2.12) and (2.13), the right hand side of (2.14a) becomes

-12a 2 1 k 4 ω s 2 β c 2 β s α c α s 2 β -s 2 α 3 sin(ϕ I +ϕ R )+4a 2 1 k 4 N 2 s 2 α + f 2 c 2 α ω s 3 β c 3 β s 2 β -s 2 α 3 sin(ϕ I -ϕ R ) . (A 1)
where the simplified notations ϕ I = (kx + n I kz -ωt) and ϕ R = (kx + n R kz -ωt)

have been introduced. As expected, nonlinear terms computed for the first-order solution result in an oscillatory term of frequency 2ω and a steady term. A particular solution of (2.14a) is therefore sought in the form

ψ part 2 = a h 2 sin(ϕ I + ϕ R ) + a s 2 sin(ϕ I -ϕ R ) , ( A 2)
the h and s superscripts standing for harmonic and steady, respectively. The first term in (A 2) corresponds to the forced wave in TH87. Substituting expression (A 2) in (2.14a) yields

a h 2 = 3a 2 1 k 2 ω D s 2 β c 2 β s α c α s 2 β -s 2 α and a s 2 = - a 2 1 k 2 ω s β c β s 2 β -s 2 α , ( A 3)
where

D = N 2 s 2 β 4s 4 β -7s 2 β s 2 α + 4s 2 α -s 2 β + f 2 c 2 β 4s 4 β -7s 2 β s 2 α + 3s 2 α . (A 4)
This expression is in agreement with Thorpe (1997) (equation ( 7)).

Since the frequency and along-slope wavenumber of the second-order harmonic motions are 2ω and 2k, respectively, the solution of the homogeneous equation must be of the form

ψ 0 2 = a 0 2 sin 2kx + m 2 kz -2ωt + ϕ 0 2 , ( A 5)
called the free wave in TH87. To ensure that this wave satisfies the dispersion relation and radiates energy away from the boundary, m 2 must be defined by

m 2 = -2s α c α (N 2 -f 2 ) + 2s β 4s 4 β N 4 (1 -4s 2 β ) + 4c 4 β N 4 (1 -4c 2 β ) + N 2 f 2 (3 -32s 2 β c 2 β ) N 2 (4s 2 β -s 2 α ) + f 2 (4c 2 β -c 2 α )
.

(A 6) Combining the particular and homogeneous solutions (A 2) and (A 5), the final solution of (2.14) is given by

ψ 2 =3a 2 1 k 2 ω s 2 β c 2 β s α c α (s 2 β -s 2 α )D sin 2kx + m 2 kz -2ωt -sin 2kx + (n I + n R )kz -2ωt - a 2 1 k 2 ω s β c β s 2 β -s 2 α sin (n I -n R )kz .
(A 7) where the constant parameters a 0 2 = -a h 2 and ϕ 0 2 = 0 have been determined through the boundary condition (2.14b).

Interestingly, as noted by [START_REF] Wunsch | Note on some reynolds stress effects of internal waves on slopes[END_REF] in a similar weakly nonlinear study of internal waves encountering a shoaling region, the system becomes degenerate when considering the steady problem in the presence of background rotation. The fields v s 2 and b s 2 are indeed linked through the equation

∂ ∂x b s 2 -f ∂ ∂z v s 2 = J(ψ 1 , ∇ 2 ψ 1 ) s (A 10)
where the subscript s on the right-hand-side refers to the steady component, but remain undetermined.

A.2. Determining the behavior at resonance

In this section, we analyse the behaviour of ψ 2 as D → 0. The point is that D appears in the denominator of (2.15) and the numerator also vanishes when D → 0. One can indeed show that

m 2 = n I + n R ⇔ N 2 (4s 2 β -s 2 α ) + f 2 (4c 2 β -c 2 α ) D = 0 , ( A 11)
so that D = 0 implies that the incident (k, n I k, ω), reflected (k, n R k, ω) waves and the second harmonic free wave (2k, m 2 k, 2ω) form a resonant triad.

A way to remove this indeterminacy is to reformulate the expression of ψ 2 as (2.18). The following expression of δ = 1 2 (m 2 -(n I + n R )) is then obtained (from where (A 11) where the h superscript stands for the harmonic component of the corresponding term as opposed to its constant one. This can be rewritten as

pt 12 = Fr 2 J(ψ 1 , ∂ z ψ 1 ) h ∂ z ψ h 2 + J(ψ 1 , v 1 ) h v h 2 + J(ψ 1 , ∂ x ψ 1 ) h ∂ x ψ h 2 + 1 N 2 J(ψ 1 , b 1 ) h b h 2 .
(B 2) Using the expressions of the first and second harmonic solutions determined in appendix A.1, one can show that the time averaged transfer power density is given by

pt 12 = 3Fr 2 a 4 1 k 6 2ω s 4 β c 4 β s α c α (s 2 β -s 2 α )D 2s α c α (N 2 -f 2 ) -(N 2 (2s 2 β + s 2 α ) + f 2 (2c 2 β + c 2 α ))m 2 sin(2δkz ) . (B 3)
The transfer power density has been determined in the weakly non-linear framework, hence with the implied hypothesis that the first harmonic wave is negligibly influenced by the second harmonic generation. In the same way, one can determine an estimate of the power density by viscous effects (diffusive ones being very small) by adding a viscous term to the RHS of equation (5.2a) but assuming that the inviscid first harmonic solution (2.9), (2.12) and (2.13) still holds. Equation (5.2a) thus takes the form

∂ ∂t u 1 = RHS + St -1 ∇ 2 u 1 , ( B 4)
where St = Re/F r is the Stokes number. Taking the dot product of eq. (B 4) by u 1 the viscous term can be rewritten as the sum of an energy diffusion term and an energy dissipation one using the identity

∇ 2 u 1 • u 1 = 1 2 ∇ 2 u 1 -∇u 1 : ∇u 1 , ( B 5)
so that the dissipated power density is

pd 1 = St -1 ∇u 1 : ∇u 1 = St -1 ∇u 1 2 + ∇v 1 2 + ∇w 1 2 . (B 6)
Finally, using the original first harmonic solution yields the time averaged dissipated power density

pd 1 = St -1 a 2 1 k 4 2 1 + f 2 ω 2 c 2 β sin 2 (β + α) sin 2 (β -α) + sin 2 (β -α) sin 2 (β + α) -2(s 2 β -c 2 β ) s 2 β -1 + f 2 ω 2 c 2 β cos ((n I -n R )kz ) . (B 7)
Having determined an expression of these power densities and in order to build the energy budget, we now need to integrate them over the interaction area. Thanks to simple geometric considerations in the interaction triangle, one can calculate the integrals, The transfer and dissipated powers within the interaction area need finally to be compared to the incident energy flux. A well-known result of the linear internal wave theory is that the energy flux can be expressed as P u = Ec g . E and c g are respectively the mechanical energy density and group velocity whose expressions for the incident wave

Figure 3 .

 3 Figure 3. Same as figure 2, for a ratio of the Coriolis frequency to the Brunt-Väisälä frequency equal to 0.2.

Figure 4 .

 4 Figure 4. Envelop function E(z) given by (4.1b). The vertical width of the incident wave is equal to n λ λz where λz is the vertical wavelength. The width of the wave is defined by the number of wavelengths n λ contained in the envelope.

Figure 5 .

 5 Figure5. Amplitude of the off-slope velocity component in mm s -1 filtered at the first harmonic frequency ω (right column) and at the second harmonic frequency 2 ω (left column) for four different slope angles α (the angle of incidence β is constant and equal to 18.9 • ). The incident and reflected waves are delineated with dashed lines. The unit of the horizontal and vertical axes is in m. The envelop of the incident wave contains two wavelengths (n λ = 2).

Figure 6 .

 6 Figure 6. Same as figure 5 for an incident wave envelop containing 4 wavelengths (n λ = 4).

Figure 7 .

 7 Figure 7. Amplitude of the second harmonic offslope velocity component at resonance as a function of the coordinate normal to the slope (z ) in the interaction area for n λ = 4, α = αres, for two different values of fluid viscosity (a) ν = 10 -6 m 2 s -1 , and (b) ν = 10 -7 m 2 s -1 . The dashed line is the theoretical prediction.

Figure 8 .

 8 Figure8. Comparison of the weakly nonlinear theory presented in section 3 with numerical simulations. The Froude number of the incident wave is Fr = 5 10 -3 , the angle of incidence is β = 18.9 • and an incident wave with three different widths is considered: n λ = 1, 2 and 4. The solid blue line is the theoretical ratio R th defined by equation (3.1) (expression (3.8) is plotted here). Its numerical counterpart R num defined by equation (4.3) is represented with symbols: stars for simulations with ν = 10 -6 m 2 s -1 and triangles for ν = 10 -7 m 2 s -1 . The dashed line indicates the value of α for which resonance is predicted by the theory of TH87.

Figure 9

 9 Figure9. Energy budget for the first harmonic in the interaction area Ω (yellow triangle in figure1) as obtained from the weakly nonlinear theory in which dissipative effects have been added. The molecular viscosity is equal to ν = 10 -6 m 2 s -1 in the first row and to ν = 10 -7 m 2 s -1 in the second row. The green and red curves represent the power averaged over one period transferred to the second harmonic and dissipated by viscous effects, respectively, as a function of the slope angle α; each power is scaled by the incident energy flux in Ω. The blue curve is the sum of the scaled powers. The fixed parameters are the incidence angle β = 18.9 • , the Froude number Fr = 5 10 -3 and the Reynolds number Re = 35 used in the numerical simulations.

  -n R )k 2 , Ω sin(2δkz ) dx dz = n λ π δk 2 1 -sinc 4n λ π δ n I -n R and Ω cos((n I -n R )kz ) dx dz = 1 -cos(n λ π) (n I -n R )k 2 . (B 8)Using these expressions and the non-dimensional parameters given by eq. (3.7), one can express the transfer and dissipated powers in Ω as a function of the angles α and β and the non-dimensional parameters Fr , St, γ and n λ :β s α c α sin(β + α) s 2 β + γ 2 c 2 β sin 3 (β -α) 2s α c α (1 -γ 2 ) -(2s 2 β + s 2 α + γ 2 (2c 2 β + c 2 α )β + α) sin 2 (β -α) + sin 2 (β -α) sin 2 (β + α) -2(1 -cos(n λ π))(s 2 β -c 2 β ) s 2 β -

  -dimensional parameters given by eq. (3.7), the incident energy flux in Ω is thusP I = ∂Ω E I c gI • dn = E I c gI n λ λ = (1 -γ 2 )s β c β 4π s 2 β + γ 2 c 2 β n λ . (B12)

Acknowledments

We thank Nicolas Grisouard, Louis Gostiaux and Joël Sommeria for useful discussions and Hidenori Aiki for his help in using the NHOES code. ML was supported by the ANR program and computations were performed at the CINES supercomputer center.

flat surface with the reflected wave generates harmonic waves n ω < N . Figures 5b and is inferred):

This expression enables to find the behaviour of δ as D → 0

and the expression of the second harmonic wave at resonance

A.

Lagrangian and Eulerian mean flows

In this section, we quickly demonstrate a statement of TH87, namely that there is no mass transport associated with the along slope steady current

appearing in (2.15). This Eulerian current is indeed found to be compensated by the Stokes drift associated with the first harmonic wave (where the incident and reflected waves superpose).

Let u 1 = (u 1 , w 1 ) = (∂ψ 1 /∂z , -∂ψ 1 /∂x ) be the velocity field of this first harmonic wave. Assuming that fluid parcel displacements are small compared with the length scale over which u 1 varies, the Stokes drift is defined by (Longuet-Higgins 1969)

and

where the overbar designates the time average over one wave period. Using the expression of ψ 1 given by (2.9) yields u sd 2 = -u s 2 and w sd 2 = 0 . (A 17)

Appendix B. Energy budget in the interaction area

In this appendix we establish the expressions of the different quantities used in the energy budget in the interaction area presented in section 5.

As explained in that section, the power density transferred from the first harmonic to the second one is given by