Thomas Le Reun

Benjamin Favier

Michael Le Bars

Parametric instability and wave turbulence driven by tidal excitation of internal waves

) which has been measured in the oceans. In addition, we find that the mixing efficiency is altered compared to what is computed in the context of DNS of stratified turbulence excited at small Froude and large buoyancy Reynolds numbers and is consistent with a superposition of waves.

Introduction

Tides affect the whole shape of planetary bodies, from their surface oceans to their liquid cores and also their subsurface oceans in the case of icy satellites such as Enceladus ( [START_REF] Thomas | Enceladus's measured physical libration requires a global subsurface ocean[END_REF], Tyler, 2009]). Such a mechanical forcing of the topography of these layers is known for driving a wide variety of bulk flows in planetary interiors (see [Le Bars et al., 2015] for a review). It can directly excite inertial waves ( [Ogilvie andLin, 2004, Goodman andLackner, 2009]), drive zonal flows ( [Tilgner, 2007, Morize et al., 2010, Sauret et al., 2014, Favier et al., 2014]), and generate turbulence from various instabilities. This is the case for instance of the elliptical instability which, from the combination of the elliptical deformation and the planetary rotation, drives the resonance of inertial waves ( [Kerswell, 2002]). Their growth saturates in either sustained turbulence ( [START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF], Le Reun et al., 2017]) or cycles of turbulent bursts in which large scale vortices emerge and inhibit the instability ( [Kerswell, 1999, Schaeffer and Le Dizès, 2010, Barker and Lithwick, 2013, Le Reun et al., 2017]). Tides can therefore inject energy into planetary fluid layers to sustain intense smallscale motions, this energy being then dissipated via viscous friction. This is particularly interesting in the case of planetary cores as it could provide an alternative stirring mechanism to thermal and solutal convection usually invoked to explain dynamo action.

The large scale tidal flow could also couple with buoyancy effects in the case where these planetary interiors are stratified. Many planetary cores are thought to be at least partly stably stratified while maintaining a large scale magnetic field (see [START_REF] Stanley | Effects of an outer thin stably stratified layer on planetary dynamos[END_REF] for a review). For instance, both experimental and numerical determination of liquid iron's thermal conductivity (see [Labrosse, 2015] and references therein) and analysis of the periodicity of the magnetic field fluctuations ( [Buffett, 2014]) point towards the existence of a stably stratifed layer at the top of the Earth's core.

Several routes exist for the excitation of three-dimensionnal turbulent motion within a stably stratified layer by tides or other mechanisms. It is a common issue in physical oceanography where the interaction of the large scale tidal flow with ground topography is known for exciting ([St. Laurent and Garrett, 2002]), focusing ([Maas andLam, 1995, Bajars et al., 2013]) and scattering ( [START_REF] Bühler | Decay of an internal tide due to random topography in the ocean[END_REF]) internal waves which break down into small scale turbulence via triadic resonant interactions [START_REF] Mackinnon | Subtropical catastrophe: Significant loss of low-mode tidal energy at 28.9[END_REF], Bourget et al., 2013, Scolan et al., 2013, Brouzet et al., 2016, Brouzet et al., 2017]). In addition, several studies have strived to examine the resonant excitation of global internal modes by a homogeneous tidal flow without

Moon

Figure 1: Schematic representation of a companion orbiting a planet at rate n. Tidal interactions induce an ellipsoidal deformation of the whole planet, which is supposed uniform: this corresponds to the behaviour of a fully deformable planet.

relying on any small-scale topography. This has been done for instance by [START_REF] Miyazaki | Three?dimensional instability of strained vortices in a stably stratified fluid[END_REF], Kerswell, 1993, Miyazaki, 1993, McWilliams and Yavneh, 1998, Aspden and Vanneste, 2009, Guimbard et al., 2010] with either radial or vertical stratification compared to deformation, but always in the situation where the Coriolis force is of greater or similar influence compared to the buoyancy effects. Whether these resonant instabilities can drive three-dimensionnal turbulence in the bulk of a stratified core or subsurface ocean in a non-rotating case remains to be seen.

In this paper, we derive a local model of fluid planetary interiors, be it a liquid core or a subsurface ocean, which allows to study the idealised limit where stratification completely dominates over rotation, with the stratification axis pointing in any direction relative to that of the tidal deformation. We show in particular that tides excite a parametric subharmonic resonance of internal waves. Moreover, such an idealised local model allows to thoroughly analyse the turbulent saturation of this tidally-driven resonance. It is thus shown hereafter that it drives bulk wave turbulence.

Note that parametric resonances of internal waves driven by large-scale homogeneous forcing have already been investigated in particular experimental setups. [START_REF] Mcewan | Parametric instability of internal gravity waves[END_REF] designed a setup to examine how large-scale internal waves spontaneously generate smaller scale oscillations in a stratified tank designed to mimic the avection by large scales. In close analogy, [Benielli and[START_REF] Benielli | [END_REF]Sommeria, 1998] showed that vertically shaking a stratified fluid leads to a parametric resonance of internal waves (similarly to the classical Faraday instability) whose growth saturates into turbulence. The originality of our study resides in the investigation of the stability of a more realistic homogeneous tidal flow. We also provide a detailed spatio-temporal analysis of the non-linear break-down into turbulence from the primary resonance.

This paper is organised as follows. The first part is devoted to introducing the tidal basic flow and developing a local approach to study its stability following the work of [START_REF] Barker | Non-linear evolution of the tidal elliptical instability in gaseous planets and stars[END_REF] and [Le Reun et al., 2017]. We carefully introduce buoyancy effects under the Boussinesq approximation. With this model, we then investigate the resonance of internal waves via direct numerical simulations and Wentzel-Kramers-Brillouin (WKB) analysis of the local model ( [START_REF] Lifschitz | Local stability conditions in fluid dynamics[END_REF]). We then thoroughly investigate the turbulence resulting from the non-linear saturation of the instability. We show that it is best described as a wave turbulence and study the subsequent mixing and dissipation rates.

2 Local study of the equilibrium base flow

Tidal base flow

We consider a non-rotating planet undergoing tidal deformation due to a moon orbiting at rate ne z as pictured in figure 1. We suppose the tidal deformation to be ellipsoidal and uniform in the whole fluid planetary interior; this assumption corresponds to planetary cores or oceans in between two boundaries whose response to tidal force is the same as that of the fluid's. The other limit, which won't be considered here, is an ocean on top of a non-deformable solid core or inside a non-deformable solid shell for which the ellipticity of the deformation can no longer be considered uniform. Any planetary fluid layer stands between these two limits.

In the reference frame following the tidal deformation, the base flow U bulge b can be approximated to the following analytical solution ( [START_REF] Sridhar | Tidal disruption of viscous bodies[END_REF], Kerswell, 2002, Barker, 2016, Barker et al., 2016]):

U bulge b = γ   0 -1 -β 0 1 -β 0 0 0 0 0     X b Y b Z b   = BX b , (1)
where we have introduced γ = -n, the rotation rate of the fluid in the orbital frame, by analogy with previous studies ([ Barker andLithwick, 2013, Le Reun et al., 2017]), and β the ellipticity of the tidal deformation defined as β = (a 2 -b 2 )/(a 2 + b 2 ) (where a and b are the semi-major and semi-minor axes respectively, see figure 1). X b is the position vector whose coordinates are written in the orbital frame rotating with the tidal bulge.

The tidal flow ( 1) is an exact steady solution of theincompressible Navier-Stokes equation in an ellipsoidal planetary fluid layer with stress-free boundary conditions, under the assumption that the tidal potential is small compared to the body self-gravity potential in its unperturbed spherical shape, and neglecting the bulge self-gravity ( [START_REF] Sridhar | Tidal disruption of viscous bodies[END_REF], Barker, 2016, Barker et al., 2016]). The base flow translates the fact that the streamlines must be lying in the total potential isosurfaces (including gravitational, tidal and rotational contributions) and must match the solid body rotation in the limit β = 0 measured in the orbital frame. They are therefore elliptical and parallel to the (xOy) plane. In planetary fluid layers, the ellipticity of the tidal flow is usually below 10 -3 ( [START_REF] Cébron | Elliptical instability in terrestrial planets and moons[END_REF]) and is about 10 -7 for the Earth's core and oceans (without taking into account any topographical effect). In the planetary reference frame the base flow U b translates, after coordinate change and velocity transformation, into ( [Goodman, 1993, Barker andLithwick, 2013]):

U b = -γβ   sin(2γt) cos(2γt) 0 cos(2γt) -sin(2γt) 0 0 0 0     X Y Z   = A(t)X. (2)
This basic flow is also a solution of the full incompressible Navier-Stokes equations. Both orbital and planetary frames of reference can be used for simulations. The orbital frame, in which the base flow is steady and writes as in equation ( 1), is widely used in global numerical simulations of tidal instabilities for it is numerically easier to keep a constant boundary topography (see for instance [START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF] and [START_REF] Cébron | Tidal instability in a rotating and differentially heated ellipsoidal shell[END_REF]).

Conversely, the planetary frame can be more intuitive: in the case of the Earth's, it corresponds to an observer standing on the solid mantle and feeling nearly two tides per day. The base flow in this frame only retains the perturbations induced by tidal deformation. Note that this basic flow is also a solution of the Navier-Stokes equations in the presence of stable stratification in the Boussinesq approximation. Assuming an equilibrium state for which isopycnals are also the surfaces with constant gravitational potential φ, including centrifugal force and tides, the density can be written as a continuous and monotonic function f of φ so that ρg = -f (φ)∇φ = -∇F (φ) with dF/dφ = f . The equilibrium buoyancy term can then be absorbed in the pressure gradient so that this basic situation is purely barotropic. Such a barotropic assumption is valid when the isopycnals can move sufficiently fast to keep track of the orbital motion of the moon and the rotating tidal potential; it is valid for high Brunt-Väisälä frequency compared to the differential rotation frequency γ (see [Ogilvie, 2014] for further discussion). In the opposite regime where the stratification is weaker, the slow motion of the isopycnals should lead to baroclinicity and excite large scale flow consistently trying to restore the alignment between isopycnal and isopotential. This paper is rather focused on small scale instabilities in the regime where the Brunt-Väisälä frequency is larger than γ. We therefore discard any baroclinic situation to keep only the global tidal distortion as a source of instability. Tidally-driven baroclinicity should deserve a study of its own.

In the following, we introduce buoyancy effects under the Boussinesq approximation. We study the tidal instability problem developing a local approach inspired by [START_REF] Barker | Non-linear evolution of the tidal elliptical instability in gaseous planets and stars[END_REF]. We use the planetary frame of reference, which is considered to be non-rotating. Introducing rotation of this frame at rate Ωe z in the model developped hereafter would only require the addition of a Coriolis force in the planetary frame. The base flow would not be modified appart from the fluid rotation rate in the orbital reference frame γ which would then write γ = Ω -n ( [START_REF] Barker | Non-linear evolution of the tidal elliptical instability in gaseous planets and stars[END_REF]). Neglecting planetary rotation is tantamount to assuming the buoyancy effects dominate the dynamics compared to the Coriolis force. Studying the interplay between rotation and stratification would require a separate study.

The local approach to the dynamics

We aim at studying the icompressible perturbations to the basic flow U b defined in (2). Instead of modelling the whole planetary fluid layer, we develop a local model to study this perturbation flow in the neighbourhood of a Lagrangian point M at position X 0 (t) following the elliptical streamlines such that Ẋ0 = U b . This local model will prove particularly convenient to include buoyancy effects as the stratification can be assumed to be locally uniform and linear around the tracked point. It has already been proposed in [START_REF] Barker | Non-linear evolution of the tidal elliptical instability in gaseous planets and stars[END_REF] to investigate the elliptical instability in rotating constant density planetary and stellar interiors. We reproduce here its derivation to carefully introduce buoyancy effects.

Let us call v i the total velocity field in the frame bound to the planet and v c the total velocity field in the frame bound to X 0 . v i satisfies the following Navier-Stokes equation:

∂ τ v i + v i • ∇ X v i = - 1 ρ ∇P i + ν∇ 2 X v i (3)
Figure 2: Trajectory of an initially square patch of fluid for an ellipticity β = 0.5 represented in the orbital frame (left) in which the bulge has a stationnary shape, and in the planetary frame (right). The four configurations correspond to the same instants. The solid black line materialises the trajectory of the center of the patch. This picture highlights the periodic stretching and shearing undergone by the lagrangian parcel.

where τ stands for time, ∇ X for the gradient in the X = (X, Y, Z) coordinates, P i is the pressure and ρ is the density of the fluid. v i is the total velocity and includes the base flow U b (X) and velocity perturbations u i so that v i = U b (X) + u i . To transform this equation into the frame in translation bound to X 0 , we process the following coordinate change:

x = X -X 0 (t) t = τ . ( 4
)
The corresponding change in derivatives is

∇ X = ∇ x = ∇ and ∂ τ = ∂ t -U b • ∇ x . The velocity measured in the frame bound to X 0 is v c = v i -U b (X 0 , t).
Transforming the equation (3) into this frame yields:

∂ t v c + ∂ t U b (X 0 ) + v c • ∇v c = - 1 ρ ∇P i (X 0 + x) + ν∇ 2 v c . (5)
The acceleration term ∂ t U b (X 0 ) is regarded as a volume force. In the frame bound to X 0 the Navier-Stokes equation reads:

∂ t v c + v c • ∇v c = - 1 ρ ∇P c -∂ t U b (X 0 ) + ν∇ 2 v c (6)
where we have introduced

P c (x) = P i (X 0 + x). As v c = v i -U b (X 0 ), it is straightforward that v c = Ax + u c with u c = u i = u.
In the neighbourhood of X 0 , the perturbed flow satisfies the following equation:

∂ t u + A(t)x • ∇u + A(t)u + u • ∇u = - 1 ρ ∇P c -∂ t U b (X 0 ) + ν∇ 2 u (7)
along with the incompressibility condition ∇ • u = 0.

Lagragian effects of the base flow

This paragraph aims at exhibiting the Lagrangian trajectory of the point M at X 0 to provide a better understanding of the model derived in the preceding paragraphs. The Lagrangian equation Ẋ0 = U b (X 0 ) can be solved analytically and the position of M at any time can be related to the initial position (X 0i , Y 0i , Z 0i ) by the following relation:

X 0 = Rot(-γt)    cos(̟t) -E sin(̟t) 0 1 E sin(̟t) cos(̟t) 0 0 0 1      X 0i Y 0i Z 0i   (8)
where Rot(-γt) is a rotation matrix of angle -γt around the Z axis. ̟ and E are defined as follows:

̟ = γ 1 -β 2 and E = 1 + β 1 -β . ( 9
)
The corresponding trajectories are plotted in figure 2 for β = 0.5, γ = 1 and an initial condition (1, 0, 0). The equivalent trajectory in the orbital frame (consistently aligned with the moon) is also indicated for comparison.

In the small β limit, it can be shown that, in the planetary frame, the Lagrangian particle rotates around the Z axis at rate -γβ 2 /2, the corresponding acceleration being negligible. Around this slow mean rotation corresponding to a Stokes drift with velocity ∼ γβ 2 ℓ with ℓ the average distance from the centre of the planetary body, the particle also accomplishes epicycles at a much higher rate γ/2 with a displacement of order βℓ and velocity βγℓ.

To materialise the local effects of the basic flow, we also plot in figure 2 the trajectories of four points forming a square pattern around the tracked point. It can be noticed that this pattern is stretched and sheared during an epicycle and rotates as the particle moves around the Z axis. Note that the slow mean rotation around the Z axis is exaggerated in figure 2 because of the very high ellipticity. This periodic stretching and shearing may drive internal wave parametric instability, in close analogy to the elliptical instability in rotating fluids ( [Kerswell, 2002]).

Introducing buoyancy in the local Navier-Stokes equation

The local model introduced above is particularly useful when introducing buoyancy effects, be it due to temperature or solute concentration. In a sufficiently small patch, the background stratification can be assumed to be uniform, while being mostly radial in global planetary layers. Let us call S(t) the background active scalar gradient such that the total scalar field T can be written as:

T = T 0 + S(t) • x + ϑ ( 10
)
where ϑ is the scalar fluctuation around the background stratification. As it will be pointed out later, it is required to include a temporal dependence in the stratification to account for periodic stretching induced by the background base flow defined in (2) and illustrated in figure 2. The buoyancy effects are first introduced in equation ( 7), in the Boussinesq approximation, via the volume effective gravitational force ρ 0 (g-∂ t U b (X 0 ))(1-α(T -T 0 )) where α is an isobaric thermal expansion coefficient. As we assume the basic state to be barotropic, the term ρ 0 (g -∂ t U b (X 0 ))(1 -αS • x) can be absorbed in a modified pressure Π. Including buoyancy effects, the equations (7) now reads:

∂ t u + A(t)x • ∇u + A(t)u + u • ∇u = -∇Π -α [ g -∂ t U b (X 0 ) ] ϑ + ν∇ 2 u . (11)
Note that the field u also satisfies the incompressibility condition ∇ • u = 0. Two possible instability sources are worth considering. In the RHS of (11), the Lagrangian advection of the studied patch translates into an effective gravity with varying intensity. Such a forcing has already been shown to parametrically excite internal waves with a growth rate proportional to the oscillating acceleration amplitude ([McEwan andRobinson, 1975, Benielli andSommeria, 1998]). The base flow is also coupled in the left hand side to the velocity perturbation; in the analogue context of rotating flows, the induced tidal stretching and shearing is well known for triggering parametric excitation of a pair of inertial waves. This instability mechanism has also been studied in the context of strained vortices with a stratification aligned with the background vorticity ( [START_REF] Miyazaki | Three?dimensional instability of strained vortices in a stably stratified fluid[END_REF], Miyazaki, 1993, Kerswell, 1993, McWilliams and Yavneh, 1998, Aspden and Vanneste, 2009, Guimbard et al., 2010]). Whether a similar stirring mechanism occurs in purely stratified fluids, i.e. with no background vorticity, has never been investigated to the best of our knowledge. In this paper, we propose to drop the gravity-driven parametric excitation to focus on tidal stretching and shearing effects.

To support this approximation, we suggest to compare the order of magnitude of the expected growth rate of both instabilities. In the case of excitation by gravity variations, the amplitude of the forcing is proportional to the acceleration of a Lagragian particle during an epicycle: βℓγ 2 . Following [START_REF] Benielli | Excitation and breaking of internal gravity waves by parametric instability[END_REF]], the corresponding growth rate σ g should scale like σ g ∼ γ(βℓγ 2 )/g where g is the mean intensity of gravity. On the other hand, if the coupling between the base flow and the velocity perturbation acts as in the elliptical instability in rotating fluids, the growth rate σ e should then scale like σ e ∼ βγ (see eq. ( 4) in [START_REF] Barker | Non-linear evolution of the tidal elliptical instability in gaseous planets and stars[END_REF]). As a consequence, σ e /σ g ∼ g/(ℓγ 2 ). As γ is at most comparable to the spin rate of the considered planetary body, this ratio equivalently compares the self-gravity of the body to the centrifugal acceleration. It should then always be large to ensure self-cohesion. For instance, for the Moon-Earth system, γ = 2π/(1day); at the Core-Mantle boundary ℓ ∼ 3 × 10 3 km and g ∼ 10 m.s -2 . The ratio σ e /σ g is about 3 × 10 2 at the Core-Mantle boundary and 1.5 × 10 2 at the surface of the Earth, which justifies dropping the varying gravity forcing term. Note however that in the case of a confined layer above a non-deformable core, the elliptical base flow, which is then no longer described by (2), would create large-amplitude lateral flows whose contribution to the effective gravity would not necessarily be negligible. Although of interest for instance for the Earth's oceans, we do not consider the latter case here to focus on a fully deformable object.

The buoyancy equation and time dependence of stratification

The Navier-Stokes equation ( 11) is coupled to the advection-diffusion equation for the scalar T :

∂ t T + (A(t)x + u) • ∇T = κ∇ 2 T ( 12
)
where κ is a diffusivity coefficient assumed to be constant. Plugging the assumption (10) in ( 12) leads to the following modified advection-diffusion equation:

dS dt • x + ∂ϑ ∂t + A(t)x • S + A(t)x • ∇ϑ + S • u + u • ∇ϑ = κ∇ 2 ϑ . ( 13
)
Assuming that in the equilibrium state there is no perturbation, i.e. (u, ϑ) = (0, 0), compels the time evolution of the stratification vector S(t) to follow:

dS dt = -A T (t)S (14)
where A T stands for the transpose of A given in equation ( 2). The periodic stretching and shearing induced by the base flow, as represented in figure 2, impacts the local background density profile. The equation ( 14) can be solved analytically to obtain that an initial stratification S 0 = (S 0x , S 0y , S 0z ) evolves in the following way:

S(t) = Rot(-γt)    cos(̟t) - 1 E sin(̟t) 0 E sin(̟t) cos(̟t) 0 0 0 1      S x0 S y0 S z0   = Rot(-γt) R(t)S 0 (15)
where E, Rot(γt) and ̟ are the same as defined in (8). A typical time evolution of S(t) is pictured in figure 3.

The initial stratification can arbitrarily be set in the (xM z) plane. It is then fully parameterised by the angle s such that S 0 = S 0 (sin s, 0, cos s) = S 0 s0 ; it represents the mean latitude at which the tracked patch is located. Note that to avoid the spontaneous appearance of baroclinic instability, the gravity has to change its direction to stay aligned with the buoyancy gradient. This is consistent with the fact that the point X 0 is following an elliptical streamline included in an equipotential surface of the total gravitational field including tidal force and self-gravitation of the unperturbed state; the gravitational field must remain perpendicular to streamlines. We arbitrarily choose to keep the gravitational field amplitude constant throughout time. Calling g the gravity intensity, the final set of local equations is, in addition to ( 14):

     ∂u ∂t + A(t)x • ∇u + A(t)u + u • ∇u = -∇Π + αg S S ϑ + ν∇ 2 u ∂ϑ ∂t + A(t)x • ∇ϑ + u • ∇ϑ = -S(t) • u + κ∇ 2 ϑ . (16)

Conclusion on the equations and the local model

We will consider hereafter a patch of typical size L. Typical time and velocity scales are then given by 1/γ and Lγ. The initial stratification amplitude S 0 can be used to build a buoyancy scale LS 0 . With those definitions, the dimensionless velocity-temperature dynamics satisfies:

         ∂u ∂t + A(t)x • ∇u + A(t)u + u • ∇u = -∇Π + N 2 ϑe s + 1 Re ∇ 2 u ∂ϑ ∂t + A(t)x • ∇ϑ + u • ∇ϑ = -S(t) • u + 1 P r Re ∇ 2 ϑ ∇ • u = 0 (17)
where we have introduced the dimensionless Brunt-Väisälä frequency N such that N 2 = αgS 0 /γ. Re is the Reynolds number L 2 γ/ν and P r is the Prandtl number ν/κ. S(t) and e s are defined as follows:

S(t) = Rot(-γt) R(t)S 0 /S 0 = Rot(-γt) R(t)s 0 and e s = S(t) / S(t) (18)
This set of equations is particularly convenient as it reduces a global problem with a non-trivial ellipsoidal geometry to a local one in Cartesian coordinates with uniform stratification. It retains all the key ingredients to understand the homogeneous dynamics of tidally-forced flows in stratified layers while avoiding to account for boundary layers.

The local model in the orbital frame

The same analysis can be performed in the orbital frame tracking the elliptical deformation rotation at rate ne z . The set of equations obtained is very similar to (17) except that A(t) must be replaced by the matrix B defined in (1) and a Coriolis acceleration 2ne z × u must be added to the left hand side of the momentum conservation. In addition, the time evolution of the stratification vector reads S(t) = R(t)s 0 where R has been introduced in (15). The two frames of reference are equivalent; nevertheless we prefer the planetary frame as it allows to introduce planetary rotation with the mere intuitive addition of a Coriolis acceleration. In the non stratified case, this frame allows to clearly identify the inertial waves frequencies ( [Barker andLithwick, 2013, Le Reun et al., 2017]). We have therefore considered the model in this frame to facilitate future works concerned with the interplay between rotation and stratification.

-1

0 1 S x (t) -1 0 1 S y (t)

Direct numerical simulations in a shearing box

The dynamics of the perturbations u to the equilibrium state can be simulated via a decomposition of {u, Π, ϑ} into Kelvin modes such that:

{u, Π, ϑ} = k ûk (t), Πk (t), θk (t) e ik(t).x . ( 19
)
In close analogy to the model developed in section 2.5 where stratification is found to be time dependent, evolving the wave vectors k accounts for the periodic shearing induced by the base flow U b . For all k, the equations ( 17) are equivalent to:

               dk dt = -A T (t) k = -A(t) k dû k dt = -A(t)û k -ik Πk + N 2 θk e s - k 2 Re ûk -(u • ∇u) k d θk dt = -S(t) • ûk - k 2 Re P r ûk -(u • ∇ϑ) k . ( 20
)
A similar derivation was first carried out by [START_REF] Barker | Non-linear evolution of the tidal elliptical instability in gaseous planets and stars[END_REF] to study the tidally-driven elliptical instability in rotating non-stratified fluids. Note that the time evolution of k(t) is the same as that of the stratification ( 14). This shearing box development allows to simulate the perturbations to the base flow with efficient pseudo-spectral methods, as originaly devised by [Rogallo, 1981], for instance implemented in the Snoopy code initiated by [START_REF] Lesur | On the relevance of subcritical hydrodynamic turbulence to accretion disk transport[END_REF]] that we use here. Apart from the time evolution of the wave vectors which is known analytically, the code solves the ODEs (20) with a fourth order Runge-Kutta method and applying a 2/3-rule for dealiasing the non-linear terms.

3 Stability analysis

Direct Numerical Simulations

We first investigate the stability of the elliptical base flow by performing direct numerical simulations of the full problem, including viscosity and non-linearities. This is done with the Snoopy code mentioned in section 2.8 which solves the equations (20). The Reynolds number Re is usually set between 10 6 and 10 7 while the Prandtl number P r is kept constant at 1. The resolution used is up to 96 modes in each direction in a square box of size L which is used as a length scale. The simulations are initiated by a broad-band noise with k/(2π) ranging from 4 to 20. The time evolution of the volume-averaged kinetic energy of the fluctuations is tracked until an exponential phase is reached from which a growth rate is derived.

WKB and Floquet analysis

Along with solving the full problem, we examine the linear inviscid limit of equations ( 17) via a Wentzel-Kramers-Brillouin analysis ( [START_REF] Lifschitz | Local stability conditions in fluid dynamics[END_REF]). It is easier to perform the stability analysis in the orbital frame, following the rotation of the elliptical deformation, where the base flow matrix B does not depend on time (see paragraph 2.7 for the corresponding change in equations ( 17)) .

Figure 4: Schematic diagram defining the relevant angles to describe the relative position between initial stratification s0 and wave-vector K 0 . The angle between s0 and the axis z is s defined in figure 3.

We then assume that the velocity, pressure and buoyancy fluctuations around the equilibrium state may be written as follows:

{u, Π, ϑ} = {a, p, Θ} e i φ(x, t) η ( 21
)
where η is a small parameter accounting for the quick wave-like spatial variations of the phase φ compared to the secular evolution of amplitudes a, p and Θ.

Plugging the assumption ( 21) into (17) in the linear inviscid limit and performing a series expansion in η lead to the following set of equations ( [START_REF] Lifschitz | Local stability conditions in fluid dynamics[END_REF]):

                   K = ∇φ dK dt = -B T K da dt = 2 KK T K 2 -I Ba -2 KK T K 2 -I (e z × a) -N 2 KK T K 2 -I Θe s dΘ dt = -S(t) • a . ( 22
)
The equation on K can be solved analytically: it follows the same time evolution as S(t), K = R(t)K 0 where R(t) has been defined in the time evolution of stratification (15) and K 0 is an initial condition vector. At the lowest order (β = 0), the linear operators B and e z × • are equal. Since the shearing and stretching effects are entirely due to the ellipticity of streamlines, the stratification S and wave-vector K have a purely rotating motion at rate γ. Taking the time derivative of the last equation in ( 22) gives a second order differential equation:

d 2 Θ dt 2 -N 2 (s 0 • K 0 ) 2 K 2 0 -1 Θ = 0 ( 23
)
where s0 is the unit vector defining the initial stratification direction, as defined in (18). It is then convenient to introduce ξ, the angle between the initial stratification and wave-vector as represented in figure 4. Θ, and consequently the velocity component sensitive to the stratification a • S, oscillate at a frequency λ = ±N sin ξ.

At the lowest order, the internal wave dispersion relation is retrieved.

To further analyse the linear growth of the instability and to obtain quantitative growth rates, we perform the Floquet analysis of the system (22) ( [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF]). This can be done since the vectors K, S and e s all oscillate at the same frequency ̟, defined in (9). The linear operator RHS of ( 22) therefore oscillates with a period 2π/̟. Knowing the time evolution of K, we solve the differential equations over the vector (a, θ) from t = 0 to t = T = 2π/̟. The initial condition is the identity matrix I. The final value (a(T ), θ(T )) for each initial condition corresponds to a monodromy matrix Φ and the growth rate σ of the instability is then related to its eigenvalues µ i such that ( [Bender andOrszag, 1978, Cébron et al., 2014]):

σ = ̟ 2π max {ln µ i } or σ = ̟ max {ln µ i } (24)
where • refers to a growth computed per tidal cycles. As the resonant wave-vector is a priori unknown, the Floquet analysis is performed for different K 0 . This initial wave vector is parameterised by its norm K 0 (which As it will be shown hereafter, the growth rate is linearly growing with the ellipticity. We therefore normalise the growth rate by β which allows to identify the limits of the Floquet resonance tongue. The latter converges towards N sin ξ = γ = 1 for β → 0 showing that the resonance is subharmonic. Note that the higher growth rate area is always above the dashed line N sin ξ = 1.

does not play any role in the inviscid limit), the angles s, ξ and α represented in figure 4, such that:

K 0 = K 0   sin ξ cos α cos s + sin s cos ξ sin ξ sin α -sin ξ cos α sin s + cos s cos ξ   . (25)
Such a parameterisation merely comes from the expression of K 0 in the spherical coordinates (α, ξ) with a polar axis (M, s0 ) (see figure 4). With s used as a control parameter, resonance is found exploring the values of σ in the (ξ, α) space.

The range of angles ξ to explore is non-trivial. We show in figure 5 the maximum (respectively to α) growth rate computed with the Floquet theory as a function of ξ and β. The non-zero growth rate area delimits Floquet tongues which stretch towards N sin ξ = γ = 1 as β goes to 0. The resonant waves are therefore parametrically excited close to half the frequency of the forcing flow. Note however that the Floquet tongues are not symmetric around N sin ξ = 1: the areas with maximum growth rates are always slightly above this line. As a consequence, to compute the theoretical maximum growth rate, we explore a range of angle ξ around arcsin(1/N ) with a tolerance of order β.

Note that the theoretical growth rate could also have been analytically computed via a multi-scale analysis where βt is the slow time. This would have given the asymptotic resonant values of α for β → 0. However, the complexity of the first order operator respective to β is such that the problem might be intractable. This might be due to the low degree of symmetry as the angle between the orbital plane and the stratification axis is arbitrary.

Comparison between DNS and linear stability analysis

We perform direct numerical simulations setting β = 0.05, N = 1.5, s = 90 • with a Reynolds number Re = 10 6.5 . A first try was initiated from a broadband white noise with k/(2π) ∈ [4,20]. The kinetic energy displays an exponential growth but snapshots reveal several entangled growing modes. To better quantify the growth rate and modes selection, we then restrict the broad-band noise to three intervals 1:[5, 10], 2:[10, 14] and 3:[14, 20] which allows to isolate growing modes with approximately the same wavelength. The kinetic energy corresponding to these three initial conditions is shown in figure 6 with snapshots of the buoyancy perturbation field ϑ during the exponential growth phase. In each case, the kinetic energy is exponentially growing and the buoyancy field bears a plane wave structure, confirming that the instability mechanism is based on wave resonance. These DNS results allow to calculate the viscous growth rate σv in tidal units. The inviscid growth rate σth (expressed in tidal cycles) is then obtained by substracting the viscous damping of the growing mode, i.e. σth = σv + 2πk 2 Re -1 with k the wave number of the mode (as σth is in tidal units, a 2π factor must be added to the damping rate).

These results are then compared to the theoretical inviscid growth rate σf given by the Floquet analysis of equations ( 22). The map of σ f (α, ξ) is displayed in figure 7 where the location of the growing mode for each DNS is highlighted by a black dot associated to the corresponding σf . We first note that, both in DNS and theory, the angle ξ satisfies the condition N sin ξ = 1 with a tolerance smaller than β, as it was expected from qualitative arguments developed in the preceding paragraph. In addition, the theoretical growth rate σf is close to the growth rate σth measured in DNS with a relative difference less than 2% . With this very good agreement between DNS and the linear theory, we can now analyse the dependence of the growth rate on the control parameters using rapid linear theory only.

Linearity with the ellipticity β

The amplitude of the periodic stretching and shearing responsible for the parametric excitation of internal waves is proportional to the ellipticity β. Another way to validate the DNS and the linear stability analysis is to examine the consequent expected linearity in β of the growth rate as in the case of tidally-driven elliptical instability ([Le Dizès, 2000, Kerswell, 2002, Grannan et al., 2017]). As shown in figure 7 (right panel), the theoretical maximal growth rate σf inferred from Floquet theory is very well described by a linear function over several orders of magnitude, in the case with N = 1.5 and s = 90 • . The growth rate computed from DNS is in addition very close to this theoretical line.

Dependence on the stratification angle s

To illustrate the dependence on the latitude or equivalently the stratification angle s (see figures 3 and 4), we show in figure 8 the maximum theoretical growth rate of the instability for fixed N = 1.5. The main conclusion is that this instability can be triggered at any latitude in a planetary fluid layer. For this Brunt-Väisälä frequency, the growth is optimal between roughly 50 • and 60 • . The mode selection of the parametric instability changes with latitude. Unlike the mode selected in the case N = 1.5 and s = 90 • where α ≃ π/2 (see figure 7), the mode selected at s = 60 • lies in the plane (xM z), i.e. α ≃ π (see figure 4).

Dependence on the Brunt-Väisälä frequency N

The growth rate of the instability is also a function of the Brunt-Väisälä frequency. As shown in figure 9, it tends to a limit value when N ≫ 1. The consequence is that the instability can be triggered at any Brunt-Väisälä frequency provided that it is larger than 1 in tidal units. Note that for large N , the selected modes' wave vectors draw closer to the stratification direction as sin ξ ∼ 1/N .

Conclusion

The linear stability analysis examined theoretically via WKB analysis and Floquet theory is quantitatively consistent with the results of direct numerical simulations. This first study has two consequences. It allows to assert that provided the dissipation is low enough, a parametric excitation of internal waves can be excited in a planetary fluid layer undergoing homogeneous tidal deformations. This instability can be triggered at any latitude. However, the mode selection seems difficult to predict as it depends on latitude. At least it is confirmed that the growing waves are selected because their frequencies are close to half the forcing frequency 2γ. In addition, this first linear study validates the use of local direct numerical simulations under a shearing box approximation with time-varying wave-vectors, as this method is in excellent agreement with the linear WKB-Floquet theory. 1: Input parameters and measured statistical properties of the flow for each simulation. k res is the principal wave number of the resonant modes emerging during the growth phase. u rms is the rms velocity computed from the mean of the kinetic energy. The Froude number F r is computed as u rms /(N λ res ). ε k is the saturation dissipation rate defined as -Re -1 ∂ i u j ∂ i u j summed over the whole box. The output Reynolds Re o number and the buoyancy Reynolds number R are respectively defined in equations ( 28) and ( 29). When errorbars are given, they correspond to the variance of the quantity over the total duration of the saturation phase. Note that N = 4 and Re = 10 6.5 is intermittently turbulent and the Reynolds number had to be pushed up to 10 7 to observe a sustained turbulence.

Non-linear saturation of the instability

To further analyse this tidally-driven instability of internal waves, we now focus on its non-linear saturation. This regime is the most relevant to the understanding of natural systems, in particular to comprehend the dissipation rate of the input tidal energy and the turbulent mixing in the oceans, or to study the possible existence of dynamo action in stably stratified planetary cores. Although we cannot address here all those issues, we strive to exhibit the key features of the non-linear saturation of this tidally-driven instability as a basis for future works. We performed many simulations, all with an ellipticity β = 5 × 10 -2 , to explore the influence of the Reynolds number Re, the Brunt-Väisälä frequency N , and the latitude or stratification angle s in the low forcing intensity and low dissipation regime relevant to geophysics. The Prandtl (or Schmidt) number P r is kept constant at P r = 1. These simulations are all summed up in table 1 where the input parameters are referenced along with the main output statistical quantities.

Sustained instability and turbulence

As in the linear stability analysis, the non-linear fate of the instability is mainly tracked via the time evolution of the total kinetic energy in the Lagrangian box. Figure 10 shows that once the instability has reached saturation, the kinetic energy is maintained throughout time for the considered parameters. As indicated by figure 10, changing the Reynolds number from 10 6 to 10 7 does not influence the time-averaged value of the [150,400]. The black solid line materialises a k -3 power law. The dashed-dotted spectra is computed during the growth phase to show that the energy is primarily injected in a narrow band of wave numbers through the instability mechanism.

10 0 10 1 10 2 k/(2π) 10 -8 10 -6 E(k) k -3 Angle s 0 • 45 • 90 • 10 0 10 1 10 2 k/(2π) 10 -8 10 -6 E θ (k) k -3
Re = 10 6.5 N = 1.5 Re = 10 6.5 N = 1.5 10 0 10 1 10 2 k/(2π) , 45 • , 90 • } at Re = 10 6.5 and N = 1.5. Bottom: the same quantities for s = 90 • with N ∈ {1.5, 2, 4}. The Reynolds number is 10 6.5 except for N = 4 where it had to be increased to 10 7 to observe sustained turbulence. kinetic energy provided that N and s are kept constant. We note though that the variations around the mean energy level are larger for the lowest values of the Reynolds number. This is reminiscent of cyclic resonance and turbulence breakdown often occurring in systems close to the instability threshold. It has been observed for instance for the elliptical instability in rotating fluids ( [START_REF] Grannan | Experimental study of global-scale turbulence in a librating ellipsoid[END_REF], Favier et al., 2015]) or in the case of parametrically excited internal waves in a Faraday instability setup ( [START_REF] Benielli | Excitation and breaking of internal gravity waves by parametric instability[END_REF]).

10 -9 10 -7 10 -5 E(k) k -3 N 1.5 2.0 4.0 10 0 10 1 10 2 k/(2π) 10 -9 10 -7 10 -5 E θ (k) k -3 s = 90 • s = 90 •
To examine whether the saturation flow is influenced by stratification, we compute the Froude number based on the resonant wavelength and the rms velocity. The resonant wavelength λ res = 2π/k res is an output parameter resulting from the mode selection during the growth phase (see table 1), and the rms value of the velocity in the saturated phase is computed from the saturation time-average of the kinetic energy E k such that u rms = 2E k . The Froude number F r is therefore defined as

F r = u rms λ res N . ( 26
)
The values of u rms and F r are all referenced in table 1. As F r = O(10 -2 ) for all the simulated configurations, we conclude that the background stable stratification strongly affects the saturated flow.

In addition, the saturation of the instability is associated with the spontaneous creation of small scales. The isotropic power spectrum of both velocity (E(k)) and buoyancy (E θ (k)) are shown in figures 10 and 11. They are computed as:

E(k) = k≤|k|<k+1 |û k | 2 and E θ (k) = k≤|k|<k+1 | θk | 2 . ( 27
)
In the high Reynolds number limit, they converge towards a k -3 power spectrum, independently of N and s. Note that a similar velocity power spectrum has been observed by [Brouzet, 2016] in the close context of turbulence driven by a forced internal wave attractor. It has also been measured at low Froude number by [START_REF] Rorai | Stably stratified turbulence in the presence of large-scale forcing[END_REF] in a stratified turbulence randomly forced at large scale. From the excitation of a few unstable internal waves, this instability mechanism manages to create sustained stratified turbulence and smaller scales. In addition, the apparent equipartition between velocity and buoyancy power spectra points towards a dynamics dominated by internal waves.

To better characterise the turbulent flow resulting from the saturation of the instability, we introduce two dimensionless parameters. We compute an output Reynolds number Re o based on the rms velocity and the resonant wavelength such that:

Re o = Re u rms λ res .

(

) 28
With this output Reynolds number, we can also compute the buoyancy Reynolds number R defined as ( [Brethouwer et al.,200

R = Re o F r 2 . ( 29
)
It compares a scale ℓ b beyond which buoyancy effects are negligible and a scale ℓ v beyond which viscosity affects the flow [START_REF] Godoy-Diana | Vertical length scale selection for pancake vortices in strongly stratified viscous fluids[END_REF], Brethouwer et al., 2007]). For instance, in the context of classical stratified turbulence, it compares the so-called Ozmidov and the Kolmogorov length scales ( [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF]).

The output Reynolds number is a O(100-1000) but because the Froude number is small, the buoyancy Reynolds number is at most O(1) (see table 1). This means that although the flow is turbulent, there is no significant range of scales where inertia dominates over buoyancy: all the non-dissipative scales are affected by the background stratification. This is drastically different from recent studies on forced stratified turbulence, which are mostly focused on the R ≫ 1 regime (see for instance [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF], Bartello and Tobias, 2013, Maffioli et al., 2016, Maffioli and Davidson, 2016]). Lastly, the instantaneous dissipation rate associated with this type of turbulence is shown in figure 10 and its mean saturation value is given in table 1. They are computed in our dimensionless units as ε k = -Re -1 ∂ i u j ∂ i u j > 0 where • is a volume averaging operator. Figure 12 sums up the evolution of this dissipation with the Reynolds number at constant N and s. The dissipation rate ε k is a decreasing function of the Reynolds number in the considered range of parameters and this decay is shallower than a Re -1 decrease. This is an additional signature of the development of turbulence as it indicates that the velocity gradients become steeper as the input Reynolds number is decreased. However, in the present range of parameters accessible with reasonable computing time, no saturation of ε k at high Re is reached. As a conclusion, at large Reynolds number, the flow resulting from the saturation of this tidally-driven instability is developing over a wide range of spatial scales from an initial resonance dominated by a most unstable wavelength. At a given Brunt-Väisälä frequency, this turbulence develops at any latitude. Typical snapshots of this turbulent state can be found in figure 13.

Internal wave turbulence

In this paragraph, we aim at thoroughly characterising the structures generated by the non-linear saturation of the initially unstable waves. In simulations and experiments of stratified turbulence, the emergence of The insert shows the same amplitude compensated by ω 2 with ω the frequency to highlight a ω -2 power law consistent with the high-frequency branch of the Garrett and Munk spectrum ( [START_REF] Garrett | Internal Waves in the Ocean[END_REF]).

layerwise structures in which the flow is quasi-two-dimensional is frequently observed. These so-called "pancakes modes" correspond to the quasi-static limit of the internal waves dynamics (i.e. ξ → 0 and ω → 0); threedimensional motion comes through shear instability between those layers (see e.g. [Billant andChomaz, 2001, Brethouwer et al., 2007]). Conversely, the turbulence excited by internal wave attractors leads to a different situation where the turbulence is a cascade of triadic resonances between the excited waves and a swarm of daughter waves. It results in an internal wave turbulence ( [START_REF] Brouzet | Energy cascade in internal-wave attractors[END_REF]).

To determine which scenario is relevant here, we propose to map the energy in the same representation as [START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF], [START_REF] Brouzet | Energy cascade in internal-wave attractors[END_REF] and [Le Reun et al., 2017], i.e. to project the energy in the spectral space along the temporal frequency ω and the angle ξ between the stratification direction S and the wave vector of a mode. This allows to determine where the energy is located around the dispersion relation of internal waves. This energy map is in the present case straightforward to draw thanks to the spectral nature of our simulations. The flow is indeed known through the velocity in spectral space ûk (t). They can be sorted by their angles ξ to obtain the quantity û(ξ, t) defined as:

û(ξ 0 , t) = k,α, ξ∈[ξ0,ξ0+∆ξ] ûk (t) with k = k   sin ξ cos α cos s + sin s cos ξ sin ξ sin α -sin ξ cos α sin s + cos s cos ξ   , ( 30
)
where α is an azimutal angle, ξ is a polar angle respective to the stratification axis as defined in figure 4 and k is the norm of the wave vector k. ∆ξ is a given tolerance to assume the angle ξ of a mode is equal to ξ 0 . A time Fourier transform is applied to û(ξ, t) to finally get û(ξ, ω). With these definitions, we remind at this point that the dispersion relation of internal waves at a given k(k, ξ, α) is ω 2 = N 2 sin 2 ξ.

The result of such a process is shown in figures 14 and 15. The most striking feature is the coincidence between the main energy locations and the dispersion relation of internal waves, similarly to [START_REF] Yarom | Experimental observation of steady inertial wave turbulence in deep rotating flows[END_REF] and [Le Reun et al., 2017] in rotating fluids, to [START_REF] Brouzet | Energy cascade in internal-wave attractors[END_REF] in stratified fluids, and to [START_REF] Aubourg | Nonlocal Resonances in Weak Turbulence of Gravity-Capillary Waves[END_REF] in the case of gravity-capillary waves. It confirms that, in the saturation phase, the non-linear interactions between the growing modes give rise to a cascade of daughter internal waves. In figure 14, it can be noticed at low Reynolds number that only a few modes emerge in the non-linear saturation. Increasing the Reynolds number leads to filling continuously the dispersion relation. Note that as energy is injected into the resonant modes only and as the Froude number is always small, the only way to create new waves is via a cascade of triadic resonances.

Secondary locations of the energy mirroring the dispersion relation of internal waves can be noticed in figures 14 and 15. Their frequencies match the relation ω = 2 -N sin ξ and are therefore associated to non-linear and non-resonant interactions between the waves of frequency ±N sin ξ and the base flow of frequency ±2.

The filling of the dispersion relation depends though on the Brunt-Väisälä frequency (figure 15). When N is increased, modes with frequency around or below γ seem to be more excited via nonlinear interactions than modes with frequency between γ and N , at least for the Reynolds numbers considered here (see N = 4 in figure 15).

In order to quantify more precisely how frequencies are excited, we propose to focus on temporal spectra obtained via the local acquisition of the three components of the velocity at several points. As shown in

E(k) k -3 N = 1.5 Re = 10 6.5 s = 45 • k ⊥ k 10 0 10 1 10 2 k/(2π) 10 -11 10 -9 10 -7 10 -5 E(k) k -3 N = 4 Re = 10 7 s = 90 • k ⊥ k
Figure 18: Anisotropic velocity power spectra as functions of k ⊥ and k for left: N = 1.5, Re = 10 6.5 and s = 45 • and right: N = 4, Re = 10 7 and s = 90 • . On each plot, the black solid line materialises a k -3 power law for comparison. figure 16, and as expected theoretically for an internal wave turbulence, there are no significant fluctuations beyond the Brunt-Väisälä frequency N . Below this frequency the excited modes are homogeneously distributed down to frequencies which are an order of magnitude smaller than both N and γ (see figure 16 left). These frequencies correspond to the lower branch (i.e. small ξ ) of the dispersion relation observed in figures 14 and 15. When N is increased, i.e. when a scale separation appears between the forcing frequency γ and N , the energy contained in the higher frequencies (N sin ξ > 1) follows a ω -2 power law. Such a trend is reminiscent of oceanographic measurements of the velocity which display a similar ω -2 power law in the range of frequencies above the tidal forcing and which is interpreted as a signature of internal wave turbulence ( [START_REF] Garrett | Space-Time scales of internal waves[END_REF], Garrett and Munk, 1975, Garrett and Munk, 1979, Levine, 2002]). To provide a definitive comparison, it would be necessary to increase N while keeping a turbulent saturation, which requires large computational time as the Reynolds number must also be increased.

As a conclusion, the tidally-driven parametric instability of internal waves saturates in a state reminiscent of "internal wave turbulence". The sustained, broadband frequency and small-scale saturation flow is composed of non-linearly interacting internal waves, although the non-linearity is weak compared to the effects of the background stratification.

Anisotropy and decoupling

Stratified turbulence has often been studied in the perspective of the emergence of layerwise, low-frequency structures, leading to a strong anisotropy and decoupling between horizontal and vertical variations. This paragraph aims at comparing the internal wave turbulence previously identified to the classical theories of stratified turbulence in the high buoyancy Reynolds number regime developed in particular by [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF], [Lindborg, 2006] and [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF].

First, figures 14 and 15 do not indicate any energy accumulation in the layerwise structures around ξ = 0 and ω = 0. To support this assertion, we show in figure 17 the kinetic energy spectrum for the layerwise modes at ξ = 0 and for the rest of the flow. As in the two cases the stratification is along the x axis, the layerwise modes are easily identified in the spectral space as their wave number k is such that k z , k y = 0. As it can be noticed, at all scales, the slow modes are subdominant. This could be confusing as layers perpendicular to the stratification can be noticed in figure 13, they however do not correspond to slow modes since they are not exactly invariant along the axes perpendicular to the stratification.

Although layerwise structures are ubiquitous in stratified turbulence excited by a random forcing or by large scale vortices, forcing waves leads to a completely different state with low energy transfers towards those particular modes. The present picture is reversed in the close context of rotating turbulence excited by inertial waves: geostrophic vortices happen to grow up to taking over the whole dynamics in absence of specific dissipative process and to strongly back-react on wave propagation ( [Barker andLithwick, 2013, Le Reun et al., 2017]). This could be linked to the fundamental mathematical difference between layerwise modes in stratified turbulence and geostrophic vortices in rotating turbulence ( [Cambon, 2001]), which leads in particular to the absence of inverse cascade in the purely stratified case ( [START_REF] Marino | Inverse cascades in rotating stratified turbulence: Fast growth of large scales[END_REF], Herbert et al., 2016]).

This result regarding slow modes has strong implications in the anisotropy of the turbulent flow. In the classical theory of stratified turbulence, the velocity power spectrum is an anisotropic function of k ⊥ = k sin ξ and k = k cos ξ. The velocity power spectrum integrated over k , i.e. along the stratification axis, E(k ⊥ ), follows a Kolmogorov-like power law k Lindborg, 2006, Brethouwer et al., 2007]). However, in the case of the internal wave turbulence presented above, both spectra E(k ⊥ ) and E(k ) follow the same power law close to k -3

-5/3 ⊥ , while conversely E(k ) ∝ k -3 ([
,⊥ , as it can be seen in figure 18. The situation is even reversed as E(k ⊥ ) is slightly below E(k ) at large k while it is expected to be dominant in classical stratified turbulence. This result shows that there is no decoupling between the horizontal and vertical variations, as observed in the frequently studied high buoyancy Reynolds number regime, which is coherent with the fact that the turbulent state considered here is a superposition of many internal waves propagating in multiple directions in a quasi-isotropic manner.

The spectra displayed in figure 18 suggest that the gradients in the direction perpendicular to the stratification are less steep than they should be if the turbulence was due to shear instability between layerwise modes. To investigate whether shear instabilities are possible in the saturated flow, we compute the local Richardson number defined as:

Ri(x, t) = N 2 1 + dϑ dz s (x) 2 dv ⊥ dz s (x) 2 . ( 31
)
where z s is a linear coordinate along the stratification axis and v ⊥ is the velocity component perpendicular to the stratification direction. It compares the local Brunt-Väisälä frequency, including buoyancy fluctuations, with the shearing rate along the stratification direction. Linear stability analysis indicates that a sheared stratified flow is unstable when Ri < Ri c = 1/4. Following [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF], we compute for several Reynolds numbers, at N = 1.5 and s = 45 • , the PDFs of the local Richardson number. As shown in figure 19, for the two lower input Reynolds numbers, there is no event likely to create shear instabilities. The buoyancy Reynolds number R being smaller than one, this is coherent with the remark of [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF] that in the low buoyancy Reynolds number regime, there should be no disturbances of Kelvin-Helmholtz type in the flow. The picture seems to change at the highest Reynolds number (10 7 ), which corresponds to our most extreme simulation where rare events with Ri < 1/4 are observed. The appearance of rare unstable events could be reminiscent of a transition towards a high buoyancy Reynolds number regime, which is further discussed in the concluding section of the present paper. Still we conclude that in the regime we explore in the present paper, the internal wave turbulence is mostly stable to shear instabilities and is unable to drive strong overturning events.

Mixing

At this stage, we know that tides are able to amplify buoyancy perturbations over a background stratification, and that this amplification saturates into an internal wave turbulence. We would like to quantify then how this turbulent state mixes the buoyancy perturbations, i.e. how it irreversibly converts potential energy perturbations into background potential energy ( [START_REF] Peltier | Mixing efficiency in stratified shear flows[END_REF]). Following the work of [START_REF] Lindborg | Vertical dispersion by stratified turbulence[END_REF], [START_REF] Salehipour | Diapycnal diffusivity, turbulent Prandtl number and mixing efficiency in Boussinesq stratified turbulence[END_REF] and [Maffioli et al., 2016], we propose to quantify the mixing via a coefficient Γ defined as:

Γ = ε p ε k (32)
where ε p = -N 2 (Re P r) -1 (∇ϑ) 2 quantifies the diffusion of the buoyancy perturbations and ε k is the kinetic energy dissipation defined earlier. This coefficient Γ was originally introduced to evaluate how turbulence induces an effective diapycnal diffusivity K = Γε k /N 2 ( [Osborn, 1980, Salehipour et al., 2016]). The expression given for ε p can be retrieved by considering that it is a potential energy dissipation ( [START_REF] Lindborg | Vertical dispersion by stratified turbulence[END_REF]).

In the derivation of the energy equations from (17), the buoyancy equation must be multiplied by N 2 ϑ to obtain the same energy transfer from velocity to buoyancy, which finally yields to the definition given earlier to ε p . 10 6 , 10 6.5 , 10 7 . PDFs are computed from snapshots of the buoyancy and velocity fields, the result presented here is the ensemble average of all the PDFs computed in the saturation phase. The number of samples is usually between 10 and 20. Maffioli et al., 2016]) as a function of the input Reynold number. The limit value in the small Froude regime obtained by [Maffioli et al., 2016] is given as a reference. Right: evolution of the mixing coefficient with N for s = 90 • and Re = 10 6.5 except at N = 4 where Re = 10 7 .

In their study, [Maffioli et al., 2016] found that, forcing a turbulence in a stratified fluid with vortices aligned with stratification, at low Froude and high buoyancy Reynolds numbers the mixing coefficient Γ converges towards 0.33. In figure 20 (left), we display the mixing coefficient Γ at N = 1.5 as a function of the input Reynolds number Re with the limit value found by [Maffioli et al., 2016] given as a reference. Despite Re is increased from 10 6 up to 10 7 we do not observe any variation of the mixing coefficient. Instead, Γ remains constant around 1, well above the limit reference value. The evolution of the mixing coefficient Γ with the Brunt-Väisälä frequency N is also computed for s = 90 • . As it can be noticed in figure 20 (right), again, Γ remains constant and around 1.

This result can be inferred from a very simple model assuming the flow is only a superposition of low amplitude internal waves with weak non-linear interactions. A single wave of frequency ω and wave vector k, {u, ϑ, Π} = {u k0 , ϑ k0 , Π k0 } e i(k•x-ωt) with ω 2 = N 2 sin ξ, must obey the following linear inviscid set of equations:

∂ t u = -∇Π + N 2 ϑe s ∂ t ϑ = -e s • u ⇒ -iω u k0 = -kΠ k0 + N 2 ϑ k0 e s -iω ϑ k0 = -e s • u k0 (33)
where e s is the stratification direction unit vector. We wish then to compute the volume averaged dissipation associated with the wave field, which is merely the sum of each single wave dissipation. For one wave only,

u * • ∇ 2 u = -k 2 |u k0 | 2 and ϑ * ∇ 2 ϑ = -k 2 |ϑ k0 | 2 .
To compute Γ, we need to relate |u k0 | 2 to |ϑ k0 | 2 which can be done for instance applying k × (k × •) to the velocity equation in (33). We then obtain the exact balance |u k0 | 2 = N 2 |ϑ k0 | 2 (which does not apply at ω = 0). Thus we find with the following simple scaling for the mixing coefficient:

Γ = 1 P r (34)
where P r is the Prandtl or Schmidt number. We retrieve for our simulations at P r = 1 that Γ = 1. The numerical result Γ = 1 should therefore be regarded as an additional signature of internal wave turbulence.

To conclude, internal wave turbulence offers a picture completely different from the classical stratified turbulence at high buoyancy Reynolds number. The flow being a superposition of low to moderate amplitude waves, the mixing coefficient is different compared to a situation where the most energetic structures are the non-propagative layerwise modes. Note that although the mixing coefficient is increased, the consequent turbulent diapycnal diffusivity should still be lower than what is measured in high buoyancy Reynolds number regime, essentially because the forcing introduced here and the associated dissipation rates are small.

Conclusion

Throughout this paper, we have shown with an idealised local Lagrangian model that tidal flows are able to drive bulk turbulence in stratified planetary fluid layers. This turbulence is driven by the parametric subharmonic resonance of unstable internal waves. The latter continuously feeds a cascade of daughter waves to create a flow which bears signatures of internal wave turbulence, in particular the focalisation of the energy along the dispersion relation of internal waves. Such a turbulent flow has already been characterised in an experimental setup designed to mimic the effects of tides on a particular topography ( [START_REF] Brouzet | Energy cascade in internal-wave attractors[END_REF]). We claim from our results that internal wave turbulence can take place homogeneously in a whole fluid layer undergoing tidal distortion, provided that the latter's amplitude is large enough to overcome dissipation and that the Brunt-Väisälä frequency N is larger than the tidal frequency γ.

In addition, our local approach provides an efficient way to numerically investigate the detailed and possibly universal properties of weakly forced internal wave turbulence in low dissipation regimes. Despite recent experimental ( [START_REF] Brouzet | Energy cascade in internal-wave attractors[END_REF]) and theoretical [START_REF] Gamba | On the wave turbulence theory for stratified flows in the ocean[END_REF]) advances, this particular state of stratified turbulence remains challenging and difficult to be compared to the classical theory of wave turbulence ( [START_REF] Zakharov | Kolmogorov Spectra of Turbulence I: Wave Turbulence[END_REF], Nazarenko, 2011]). As in rotating turbulence, this is essentially due to the anisotropy of the dispersion relation (see for instance the discussion in [Brouzet, 2016]), the role played by nearresonant interactions and the non-linear interaction with non-propagative modes ( [Cambon, 2001, Galtier, 2003, Smith and Lee, 2005, Bellet et al., 2006, Scott, 2014, Gelash et al., 2017, Gamba et al., 2017]). Although our model is introduced in a targeted geophysical context, it could be used to test universal internal wave turbulence models or closure.

Future work will strive to introduce rotation, which is also a key ingredient to planetary fluid dynamics. In particular, in the limit where buoyancy effects still dominate over the Coriolis force, it should be interesting to study the consequence of its introduction on the filling of the dispersion relation and the subsequent lowfrequency cut-off. As it can be noticed in figure 16, for N = 4, there is an energy accumulation at frequencies ω ∈ [0.1, 1]. If the rotation rate were to be in this range, would the energy accumulate in the lowest frequency modes, i.e. the layerwise structures which were never observed to develop in our simulations? If so, tidal flows would convey energy into modes which could then undergo shear instabilities and therefore drive more intense Figure 21: Schematic relative distribution of the low and high buoyancy Reynolds number R regimes as a function of the input Reynolds number Re and the ellipticity of the deformation β. The area covered by geophysical regimes and our simulations is indicated in particular to highlight the fact that the low R regime is also relevant to planetary layers in the case of bulk wave turbulence excited by tides.

turbulence with enhanced dissipation rate and mixing. Whether this turbulence is sustained or lead to the temporary inhibition of the instability which feeds it remains to be seen.

Moreover, it would be interesting to investigate the persistence of the results found here in the regime of high or low Prandtl (or Schmidt) number, which are both relevant to geophysical fluid dynamics. Note though that it is already known from the experiments of [START_REF] Brouzet | Energy cascade in internal-wave attractors[END_REF] that internal wave turbulence can be excited in salted water, i.e. at high Schmidt number. As linear internal waves are characterised by energy equipartition (see paragraph 4.4), we should expect that the resonant energy transfer towards small scale is inhibited as soon as either viscosity or diffusion balances non-linear advective transfer. What happens to the non-dissipated quantity and how it interacts with the larger scale waves beyond this cut-off remains an open question. The Prandtl (or Schmidt) number should not play any significant role in the large scale behaviour of the flow.

Lastly, the type of turbulence resulting from the saturation of the tidally-driven instability occurs in a regime of high Reynolds, low Froude and low buoyancy Reynolds (R) numbers, which leads to a completely different picture compared to the high R regime frequently studied and branded as the regime relevant to geophysical flows. In particular, the mixing coefficient is increased in the regime we describe, and is coherent with the scaling Γ = 1/P r that we have derived theoretically assuming the flow is a superposition of linear internal waves only. This result is an additional signature of wave turbulence. However, this enhanced mixing coefficient may not result in an increase in the turbulent diapycnal diffusivity as the forcing and the dissipation are small.

A regime of high buoyancy Reynolds number turbulence excited by the parametric instability introduced here is possible in the very high Reynolds number limit, but could not be investigated because it is highly demanding in resolution and computational time, or it requires an increase of the ellipticity β to unrealistic values. As seen earlier, events with Richardson number under 1/4 are measured in the most extreme simulation where the dissipation is so small that the buoyancy Reynolds number reaches ∼ 5. It would be interesting to see whether at high R wave turbulence can drive strong over-turning events or not, and how it would impact the mixing coefficient and the turbulent diapycnal viscosity. Still we claim that both regimes should be considered as relevant to geophysical applications due to the specificity of our forcing mechanism favouring weak wave interactions. Indeed, the buoyancy Reynolds number R can be expanded as

R = u 3 rms Re λ res N 2 . ( 35
)
Assuming that the saturation results from the balance between the forcing term A(t)u ∼ βu rms and the nonlinear term u • ∇u ∼ u 2 rms /λ res leads to u rms ∼ βλ res . As a result, the buoyancy Reynolds number goes like:

R ∼ β 3 Re λ res N 2 . ( 36
)
The area with high R lies above a line β ∝ Re -1/3 (N/λ res ) 2/3 . In addition, the instability grows when the forcing overcomes the volume viscous dissipation, i.e. for β (λ 2 res Re) -1 . As a result, in the (β, Re) plane, both regimes are worth considering in the geophysical limit where usually β is smaller than 10 -3 and Re is large, as indicated in figure 21. Note that this discussion is unchanged if we consider the dissipation to be due to solid wall friction, for which the unstable zone lies above the line β ∝ Re -1/2 . Moreover, as indicated in figure 21, the area of small R is extended as N is increased. In future work, it would be interesting to delimit more precisely those two regimes. Note that a possible transition could be approached in our most extreme simulation for which R ∼ 5. Exploring the internal wave turbulence driven at high buoyancy Reynolds number would therefore require increasing the ellipticity and thus the forcing intensity. This, we believe, deserves a study of its own. Lastly, we believe the results presented here should not change as the ellipticity is lowered provided that R < 1 and the flow is unstable. In addition, when three waves of frequencies (ω 1 , ω 2 , ω 3 ) exchange energy via triadic resonance, the resonance condition on frequency must be satisfied with a tolerance O(F r) i.e ω 1 ± ω 2 ± ω 3 = O(F r) due to detuning by larger scales advection (see [START_REF] Smith | On near resonances and symmetry breaking in forced rotating flows at moderate Rossby number[END_REF] for a discussion in the analog context of inertial waves in rotating flows). As u rms scales like βλ res , decreasing β corresponds to decreasing F r and therefore to more exact resonances. The only significant change, we believe, is a thiner focalisation of the energy along the dispersion relation of internal waves.

Figure 3 :

 3 Figure 3: Left: schematic picture linking stratification seen from a global and a local view. The angle s between stratification and (OZ) is then a proxy for the latitude at which the patch is located. Note that the local axes (M xyz) are in translation around the Z axis. Right: time evolution of the x and y components of the stratification. Note that the slow rotation of S is due to the translative motion of the Lagrangian particle around Z.

NNFigure 5 :

 5 Figure 5: Map of the growth rate computed via the Floquet theory σf as a function of the ellipticity β and the angle ξ defined in figure 4 and equation (25) for left: N = 1.5 and s = 90 • and right: N = 4 and s = 90 • .As it will be shown hereafter, the growth rate is linearly growing with the ellipticity. We therefore normalise the growth rate by β which allows to identify the limits of the Floquet resonance tongue. The latter converges towards N sin ξ = γ = 1 for β → 0 showing that the resonance is subharmonic. Note that the higher growth rate area is always above the dashed line N sin ξ = 1.

Figure 6 :

 6 Figure 6: Growth of the instability for β = 5 × 10 -2 , N = 1.5 and s = 90 • initiated from broad-band white noise with k/(2π) in three different intervals 1:[5, 10], 2:[10, 14] and 3:[14, 20]. The viscous growth rate is obtained by a linear fitting in log-lin coordinates (shown as dashed-dotted lines). The corresponding snapshot of the buoyancy perturbation field is given for each experiment with an indication of the stratification direction. It is used to determine the growing mode.

Figure 7 :Figure 8 :

 78 Figure7: Left: map of the growth rate σf computed with the Floquet analysis as a function of the angles ξ and α. The black dots correspond to the location of the growing modes observed in figure6for which the theoretical σf growth rate is given. The dashed black lines highlight the marginal stability. Right: comparison between systematic calculation of the maximum growth rate from the Floquet analysis and DNS results. The growth rates σf computed with Floquet theory are all aligned on σf = 1.65β.

Figure 9 :

 9 Figure 9: Left: evolution of the maximum growth rate for β = 5 × 10 -2 at three different stratification angles; a few simulations have been carried out to confirm the agreement between theory and numerics. Right: map of the growth rate computed from the Floquet analysis in the case N = 4 with s = 90 • . The black dashed lines highlight the marginal stability

Figure 10 :

 10 Figure 10: Top: Time evolution of the kinetic energy and the dissipation ε k when Re is varied while keeping constant s = 45 • and N = 1.5. The time is normalised by the theoretical growth rate σ = σf . Bottom: corresponding power spectra of the velocity E(k) and of the buoyancy E θ (k); they are averaged for σt ∈[150, 400]. The black solid line materialises a k -3 power law. The dashed-dotted spectra is computed during the growth phase to show that the energy is primarily injected in a narrow band of wave numbers through the instability mechanism.

Figure 11 :

 11 Figure11: Top: power spectra of the velocity E(k) and of the buoyancy E θ (k) for s ∈ {0 • , 45 • , 90 • } at Re = 10 6.5 and N = 1.5. Bottom: the same quantities for s = 90 • with N ∈ {1.5, 2, 4}. The Reynolds number is 10 6.5 except for N = 4 where it had to be increased to 10 7 to observe sustained turbulence.

Figure 12 :

 12 Figure 12: Mean value of the dissipation rate ε k as a function of the Reynolds number for N = 1.5 and s = 45 • . The errorbars account for the standard deviation of ε k during the saturation phase.

Figure 13 :

 13 Figure 13: Typical snapshots of the buoyancy field (a. and c.) and the y component of the vorticity (b. and d.) in the saturated phase for a., b.: N = 1.5, s = 45 • and Re = 10 7 ; c., d.: N = 4, s = 90 • and Re = 10 7 .

Figure 14 :Figure 15 :Figure 16 :

 141516 Figure14: Map of the kinetic energy as a function of the frequency of the modes and the angle ξ between the stratification and the wave vector of a mode for Re = 10 6 , 10 6.5 and 10 7 at N = 1.5 and s = 45 • kept constant. The Fourier transforms are performed for σt ∈ [200, 400] and the energy is normalised by its maximum value. The range of wave number k over which the transform is processed goes from k = 3 to k = 100; changing these boundaries does not affect the map provided the most energetic scales (k 20) are included. The horizontal line represents the frequency of the first excited modes, the plain line gives the dispersion relation of the internal waves ω = N sin ξ and the dashed line locates the modes due to non-resonant non-linear interactions between the tidal basic flow and the internal waves resulting in frequencies 2 -N sin ξ. The Snoopy code computes the time evolution of half the spectral space as the fields are real, the angle ξ is therefore between 0 and π -π/4 = 3π/4.

  Figure 17: Kinetic energy spectrum of the layerwise slow modes defined by ξ = 0 and the rest of the flow for left: N = 1.5, Re = 10 6.5 and right: N = 4, Re = 10 7 at s = 90 • .

Figure 19 :

 19 Figure 19: Probability density function of the Richardson number for s = 45 • and N = 1.5 with Re ∈10 6 , 10 6.5 , 10 7 . PDFs are computed from snapshots of the buoyancy and velocity fields, the result presented here is the ensemble average of all the PDFs computed in the saturation phase. The number of samples is usually between 10 and 20.

Figure 20 :

 20 Figure20: Left: mixing coefficient Γ as defined by([Maffioli et al., 2016]) as a function of the input Reynold number. The limit value in the small Froude regime obtained by[Maffioli et al., 2016] is given as a reference. Right: evolution of the mixing coefficient with N for s = 90 • and Re = 10 6.5 except at N = 4 where Re = 10 7 .

Acknowledgement: We acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 681835-FLUDYCO-ERC-2015-CoG). We also acknowledge support from IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) for computational time on Turing (Projects No. 100508 and 100614) and from the HPC resources of Aix-Marseille Université (Projects No.15b011 and 16b020) financed by the project Equip@Meso (No. ANR-10-EQPX-29-01) of the program Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.