2.1 Tidal base flow

We consider a non-rotating planet undergoing tidal deformation due to a moon orbiting at rate $n\boldsymbol{e}_{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. The other limit, which will not 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 $\boldsymbol{U}_{b}^{bulge}$ can be approximated to the following analytical solution (Sridhar & Tremaine Reference Sridhar and Tremaine1992; Kerswell Reference Kerswell2002; Barker Reference Barker2016; Barker, Braviner & Ogilvie Reference Barker, Braviner and Ogilvie2016):

where we have introduced $\unicode[STIX]{x1D6FE}=-n$ , the rotation rate of the fluid in the orbital frame, by analogy with previous studies (Barker & Lithwick Reference Barker and Lithwick2013; Le Reun et al. Reference Le Reun, Favier, Barker and Le Bars2017), and $\unicode[STIX]{x1D6FD}$ the ellipticity of the tidal deformation defined as $\unicode[STIX]{x1D6FD}=(a^{2}-b^{2})/(a^{2}+b^{2})$ (where $a$ and $b$ are the semi-major and semi-minor axes respectively, see figure 1). $\boldsymbol{X}^{b}$ is the position vector whose coordinates are written in the orbital frame rotating with the tidal bulge. The tidal flow (2.1) is an exact steady solution of the incompressible 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 (Sridhar & Tremaine Reference Sridhar and Tremaine1992; Barker Reference Barker2016; Barker et al. Reference Barker, Braviner and Ogilvie2016). The base flow translates to 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 $\unicode[STIX]{x1D6FD}=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}$ (Cébron et al. Reference Cébron, Le Bars, Moutou and Le Gal2012) and is approximately $10^{-7}$ for the Earth’s core and oceans (without taking into account any topographical effect). In the planetary reference frame the base flow $\boldsymbol{U}_{b}$ translates, after coordinate change and velocity transformation, into (Goodman Reference Goodman1993; Barker & Lithwick Reference Barker and Lithwick2013):

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 (2.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 Grannan et al. (Reference Grannan, Favier, Le Bars and Aurnou2017) and Cébron, Maubert & Le Bars (Reference Cébron, Maubert and Le Bars2010)). Conversely, the planetary frame can be more intuitive: in the case of the Earth, 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.

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.

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 $\unicode[STIX]{x1D719}$ , including centrifugal force and tides, the density can be written as a continuous and monotonic function $f$ of $\unicode[STIX]{x1D719}$ so that $\unicode[STIX]{x1D70C}\boldsymbol{g}=-f(\unicode[STIX]{x1D719})\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=-\unicode[STIX]{x1D735}F(\unicode[STIX]{x1D719})$ with $\text{d}F/\text{d}\unicode[STIX]{x1D719}=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 $\unicode[STIX]{x1D6FE}$ (see Ogilvie (Reference Ogilvie2014) 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 $\unicode[STIX]{x1D6FE}$ . 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 Barker & Lithwick (Reference Barker and Lithwick2013). We use the planetary frame of reference, which is considered to be non-rotating. Introducing rotation of this frame at rate $\unicode[STIX]{x1D6FA}\boldsymbol{e}_{\boldsymbol{z}}$ in the model developed hereafter would only require the addition of a Coriolis force in the planetary frame. The base flow would not be modified apart from the fluid rotation rate in the orbital reference frame $\unicode[STIX]{x1D6FE}$ which would then write $\unicode[STIX]{x1D6FE}=\unicode[STIX]{x1D6FA}-n$ (Barker & Lithwick Reference Barker and Lithwick2013). 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.

2.2 The local approach to the dynamics

We aim at studying the incompressible perturbations to the basic flow $\boldsymbol{U}_{b}$ defined in (2.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 $\boldsymbol{X}_{0}(t)$ following the streamlines defined as $\dot{\boldsymbol{X}_{0}}=\boldsymbol{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 Barker & Lithwick (Reference Barker and Lithwick2013) 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 $\boldsymbol{v}^{i}$ the total velocity field in the frame bound to the planet and $\boldsymbol{v}^{c}$ the total velocity field in the frame bound to $\boldsymbol{X}_{0}$ . $\boldsymbol{v}^{i}$ satisfies the following Navier–Stokes equation:

where $\unicode[STIX]{x1D70F}$ stands for time, $\unicode[STIX]{x1D735}_{\boldsymbol{X}}$ for the gradient in the $\boldsymbol{X}=(X,Y,Z)$ coordinates, $P^{i}$ is the pressure and $\unicode[STIX]{x1D70C}$ is the density of the fluid. $\boldsymbol{v}^{i}$ is the total velocity and includes the base flow $\boldsymbol{U}_{b}(\boldsymbol{X})$ and velocity perturbations $\boldsymbol{u}^{i}$ so that $\boldsymbol{v}^{i}=\boldsymbol{U}_{b}(\boldsymbol{X})+\boldsymbol{u}^{i}$ . To transform this equation into the frame in translation bound to $\boldsymbol{X}_{0}$ , we process the following coordinate change:

The corresponding changes in derivatives are $\unicode[STIX]{x1D735}_{\boldsymbol{X}}=\unicode[STIX]{x1D735}_{\boldsymbol{x}}=\unicode[STIX]{x1D735}$ and $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x2202}_{t}-\boldsymbol{U}_{b}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}$ . The velocity measured in the frame bound to $\boldsymbol{X}_{0}$ is $\boldsymbol{v}^{c}=\boldsymbol{v}^{i}-\boldsymbol{U}_{b}(\boldsymbol{X}_{0},t)$ . Transforming the equation (2.3) into this frame yields:

The acceleration term $\unicode[STIX]{x2202}_{t}\boldsymbol{U}_{b}(\boldsymbol{X}_{0})$ is regarded as a volume force. In the frame bound to $\boldsymbol{X}_{0}$ the Navier–Stokes equation reads:

where we have introduced $P^{c}(\boldsymbol{x})=P^{i}(\boldsymbol{X}_{0}+\boldsymbol{x})$ . As $\boldsymbol{v}^{c}=\boldsymbol{v}^{i}-\boldsymbol{U}_{b}(\boldsymbol{X}_{0})$ , it is straightforward that $\boldsymbol{v}^{c}=\unicode[STIX]{x1D63C}\,\boldsymbol{x}+\boldsymbol{u}^{c}$ with $\boldsymbol{u}^{c}=\boldsymbol{u}^{i}=\boldsymbol{u}$ . In the neighbourhood of $\boldsymbol{X}_{0}$ , the perturbed flow satisfies the following equation:

along with the incompressibility condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ .

2.3 Lagrangian effects of the base flow

This paragraph aims at exhibiting the Lagrangian trajectory of the point $M$ at $\boldsymbol{X}_{0}$ to provide a better understanding of the model derived in the preceding paragraphs. The Lagrangian equation $\dot{\boldsymbol{X}_{0}}=\boldsymbol{U}_{b}(\boldsymbol{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:

where $\text{Rot}(-\unicode[STIX]{x1D6FE}t)$ is a rotation matrix of angle $-\unicode[STIX]{x1D6FE}t$ around the $Z$ axis. $\unicode[STIX]{x1D71B}$ and ${\mathcal{E}}$ are defined as follows:

The corresponding trajectories are plotted in figure 2 for $\unicode[STIX]{x1D6FD}=0.5$ , $\unicode[STIX]{x1D6FE}=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 $\unicode[STIX]{x1D6FD}$ limit, it can be shown that, in the planetary frame, the Lagrangian particle rotates around the $Z$ axis at rate $-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FD}^{2}/2$ , the corresponding acceleration being negligible. Around this slow mean rotation corresponding to a Stokes drift with velocity ${\sim}\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FD}^{2}\ell$ with $\ell$ the average distance from the centre of the planetary body, the particle also accomplishes epicycles at a much higher rate $\unicode[STIX]{x1D6FE}/2$ with a displacement of order $\unicode[STIX]{x1D6FD}\ell$ and velocity $\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}\ell$ .

Figure 2.

Trajectory of an initially square patch of fluid for an ellipticity $\unicode[STIX]{x1D6FD}=0.5$ represented in the orbital frame (a) in which the bulge has a stationary shape, and in the planetary frame (b). The four configurations correspond to the same instants. The solid black line materialises the trajectory of the centre of the patch. This picture highlights the periodic stretching and shearing undergone by the Lagrangian parcel.

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 Reference Kerswell2002).

2.4 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 $\boldsymbol{S}(t)$ the background active scalar gradient such that the total scalar field $T$ can be written as:

where $\unicode[STIX]{x1D717}$ 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.2) and illustrated in figure 2.

The buoyancy effects are first introduced in (2.7), in the Boussinesq approximation, via the volume effective gravitational force $\unicode[STIX]{x1D70C}_{0}(\boldsymbol{g}-\unicode[STIX]{x2202}_{t}\boldsymbol{U}_{b}(\boldsymbol{X}_{0}))(1-\unicode[STIX]{x1D6FC}(T-T_{0}))$ where $\unicode[STIX]{x1D6FC}$ is an isobaric thermal expansion coefficient. As we assume the basic state to be barotropic, the term $\unicode[STIX]{x1D70C}_{0}(\boldsymbol{g}-\unicode[STIX]{x2202}_{t}\boldsymbol{U}_{b}(\boldsymbol{X}_{0}))(1-\unicode[STIX]{x1D6FC}\boldsymbol{S}\boldsymbol{\cdot }\boldsymbol{x})$ can be absorbed in a modified pressure $\unicode[STIX]{x1D6F1}$ . Including buoyancy effects, equation (2.7) now reads:

Note that the field $\boldsymbol{u}$ also satisfies the incompressibility condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$ .

Two possible instability sources are worth considering. On the right-hand side of (2.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 & Robinson Reference McEwan and Robinson1975; Benielli & Sommeria Reference Benielli and Sommeria1998). The base flow is also coupled on 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 (Miyazaki & Fukumoto Reference Miyazaki and Fukumoto1992; Kerswell Reference Kerswell1993; Miyazaki Reference Miyazaki1993; McWilliams & Yavneh Reference McWilliams and Yavneh1998; Aspden & Vanneste Reference Aspden and Vanneste2009; Guimbard et al. Reference Guimbard, Le Dizès, Le Bars, Le Gal and Leblanc2010). 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 comparing 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 Lagrangian particle during an epicycle: $\unicode[STIX]{x1D6FD}\ell \unicode[STIX]{x1D6FE}^{2}$ . Following Benielli & Sommeria (Reference Benielli and Sommeria1998), the corresponding growth rate $\unicode[STIX]{x1D70E}_{g}$ should scale like $\unicode[STIX]{x1D70E}_{g}\sim \unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D6FD}\ell \unicode[STIX]{x1D6FE}^{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 $\unicode[STIX]{x1D70E}_{e}$ should then scale like $\unicode[STIX]{x1D70E}_{e}\sim \unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FE}$ (see equation (4) in Barker & Lithwick (Reference Barker and Lithwick2013)). As a consequence, $\unicode[STIX]{x1D70E}_{e}/\unicode[STIX]{x1D70E}_{g}\sim g/(\ell \unicode[STIX]{x1D6FE}^{2})$ . As $\unicode[STIX]{x1D6FE}$ 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, $\unicode[STIX]{x1D6FE}=2\unicode[STIX]{x03C0}/(1\text{day})$ ; at the core–mantle boundary $\ell \sim 3\times 10^{3}$ km and $g\sim 10~\text{m}~\text{s}^{-2}$ . The ratio $\unicode[STIX]{x1D70E}_{e}/\unicode[STIX]{x1D70E}_{g}$ is about $3\times 10^{2}$ at the core–mantle boundary and $1.5\times 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.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.

2.5 The buoyancy equation and time dependence of stratification

The Navier–Stokes equation (2.11) is coupled to the advection–diffusion equation for the scalar $T$ :

where $\unicode[STIX]{x1D705}$ is a diffusivity coefficient assumed to be constant. Plugging the assumption (2.10) in (2.12) leads to the following modified advection–diffusion equation:

Assuming that in the equilibrium state there is no perturbation, i.e. $(\boldsymbol{u},\unicode[STIX]{x1D717})=(\mathbf{0},0)$ , compels the time evolution of the stratification vector $\boldsymbol{S}(t)$ to follow:

where $\unicode[STIX]{x1D63C}^{\text{T}}$ stands for the transpose of $\unicode[STIX]{x1D63C}$ given in (2.2). The periodic stretching and shearing induced by the base flow, as represented in figure 2, impacts the local background density profile. The equation (2.14) can be solved analytically to obtain that an initial stratification $\boldsymbol{S}_{0}=(S_{0x},S_{0y},S_{0z})$ evolves in the following way:

where ${\mathcal{E}}$ , $\text{Rot}(\unicode[STIX]{x1D6FE}t)$ and $\unicode[STIX]{x1D71B}$ are the same as defined in (2.8). A typical time evolution of $\boldsymbol{S}(t)$ is pictured in figure 3. The initial stratification can arbitrarily be set in the $(xMz)$ plane. It is then fully parameterised by the angle $s$ such that $\boldsymbol{S}_{0}=S_{0}(\sin s,0,\cos s)=S_{0}\tilde{\boldsymbol{s}}_{0}$ ; it represents the mean latitude at which the tracked patch is located.

Figure 3.

(a) 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 $(Mxyz)$ are in translation around the $Z$ axis. (b) Time evolution of the $x$ and $y$ components of the stratification. Note that the slow rotation of $\boldsymbol{S}$ is due to the translative motion of the Lagrangian particle around $Z$ .

Note that to avoid the spontaneous appearance of baroclinic instability, gravity has to change its direction to stay aligned with the buoyancy gradient. This is consistent with the fact that the point $\boldsymbol{X}_{0}$ is following an elliptical streamline (in the bulge frame) 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 (2.14):

2.6 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/\unicode[STIX]{x1D6FE}$ and $L\unicode[STIX]{x1D6FE}$ . The initial stratification amplitude $S_{0}$ can be used to build a buoyancy scale $LS_{0}$ . With those definitions, the dimensionless dynamics satisfies:

where we have introduced the dimensionless Brunt–Väisälä frequency $N$ such that $N^{2}=\unicode[STIX]{x1D6FC}gS_{0}/\unicode[STIX]{x1D6FE}$ . $Re$ is the Reynolds number $L^{2}\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D708}$ and $Pr$ is the Prandtl number $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ . $\tilde{\boldsymbol{S}}(t)$ and $\boldsymbol{e}_{s}$ are defined as follows:

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.

2.7 The local model in the orbital frame

The same analysis can be performed in the orbital frame tracking the elliptical deformation rotation at rate $n\boldsymbol{e}_{z}$ . The set of equations obtained is very similar to (2.17) except that $\unicode[STIX]{x1D63C}(t)$ must be replaced by the matrix $\unicode[STIX]{x1D63D}$ defined in (2.1) and a Coriolis acceleration $2n\boldsymbol{e}_{z}\times \boldsymbol{u}$ must be added to the left-hand side of the momentum conservation. In addition, the time evolution of the stratification vector reads $\tilde{\boldsymbol{S}}(t)=\hat{\unicode[STIX]{x1D64D}}(t)\tilde{\boldsymbol{s}}_{0}$ where $\hat{\unicode[STIX]{x1D64D}}$ has been introduced in (2.15). The two frames of reference are equivalent; nevertheless we prefer the planetary frame as it allows us to introduce planetary rotation with the more intuitive addition of a Coriolis acceleration. In the non-stratified case, this frame allows us to clearly identify the inertial waves frequencies (Barker & Lithwick Reference Barker and Lithwick2013; Le Reun et al. Reference Le Reun, Favier, Barker and Le Bars2017). We have therefore considered the model in this frame to facilitate future works concerned with the interplay between rotation and stratification.

2.8 Direct numerical simulations in a shearing box

The dynamics of the perturbations $\boldsymbol{u}$ to the equilibrium state can be simulated via a decomposition of $\left\{\boldsymbol{u},\unicode[STIX]{x1D6F1},\unicode[STIX]{x1D717}\right\}$ into Kelvin modes such that:

In close analogy to the model developed in § 2.5 where stratification is found to be time dependent, evolving the wave vectors $\boldsymbol{k}$ accounts for the periodic shearing induced by the base flow $\boldsymbol{U}_{b}$ . For all $\boldsymbol{k}$ , the equations (2.17) are equivalent to:

A similar derivation was first carried out by Barker & Lithwick (Reference Barker and Lithwick2013) to study the tidally driven elliptical instability in rotating non-stratified fluids. Note that the time evolution of $\boldsymbol{k}(t)$ is the same as that of the stratification (2.14). This shearing box development allows us to simulate the perturbations to the base flow with efficient pseudo-spectral methods, as originally devised by Rogallo (Reference Rogallo1981), for instance implemented in the Snoopy code initiated by Lesur & Longaretti (Reference Lesur and Longaretti2005) that we use here. Apart from the time evolution of the wave vectors which is known analytically, the code solves the equations (2.20) with a fourth-order Runge–Kutta method and applying a $2/3$ -rule for dealiasing the nonlinear terms.

3.1 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 nonlinearities. This is done with the Snoopy code mentioned in § 2.8 which solves the equations (2.20). The Reynolds number $Re$ is usually set between $10^{6}$ and $10^{7}$ while the Prandtl number $Pr$ 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 broadband noise with $k/(2\unicode[STIX]{x03C0})$ 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.

3.2 WKB and Floquet analysis

Along with solving the full problem, we examine the linear inviscid limit of (2.17) via a Wentzel–Kramers–Brillouin analysis (Lifschitz & Hameiri Reference Lifschitz and Hameiri1991). It is easier to perform the stability analysis in the orbital frame, following the rotation of the elliptical deformation, where the base flow matrix $\unicode[STIX]{x1D63D}$ does not depend on time (see § 2.7 for the corresponding change in (2.17)) .

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

where $\unicode[STIX]{x1D702}$ is a small parameter accounting for the quick wave-like spatial variations of the phase $\unicode[STIX]{x1D719}$ compared to the secular evolution of amplitudes $\boldsymbol{a}$ , $p$ and $\unicode[STIX]{x1D6E9}$ .

Plugging the assumption (3.1) into (2.17) in the linear inviscid limit and performing a series expansion in $\unicode[STIX]{x1D702}$ leads to the following set of equations (Lifschitz & Hameiri Reference Lifschitz and Hameiri1991):

The equation on $\boldsymbol{{\mathcal{K}}}$ can be solved analytically: it follows the same time evolution as $\tilde{\boldsymbol{S}}(t)$ , $\boldsymbol{{\mathcal{K}}}=\hat{\unicode[STIX]{x1D64D}}(t)\boldsymbol{{\mathcal{K}}}_{0}$ where $\hat{\unicode[STIX]{x1D64D}}(t)$ has been defined in the time evolution of stratification (2.15) and $\boldsymbol{{\mathcal{K}}}_{0}$ is an initial condition vector.

At the lowest order ( $\unicode[STIX]{x1D6FD}=0$ ), the linear operators $\unicode[STIX]{x1D63D}$ and $\boldsymbol{e}_{z}\times \cdot$ are equal. Since the shearing and stretching effects are entirely due to the ellipticity of streamlines, the stratification $\tilde{\boldsymbol{S}}$ and wave vector $\boldsymbol{{\mathcal{K}}}$ have a purely rotating motion at rate $\unicode[STIX]{x1D6FE}$ . Taking the time derivative of the last equation in (3.2) gives a second-order differential equation:

where $\tilde{\boldsymbol{s}}_{0}$ is the unit vector defining the initial stratification direction, as defined in (2.18). It is then convenient to introduce $\unicode[STIX]{x1D709}$ , the angle between the initial stratification and wave vector as represented in figure 4. $\unicode[STIX]{x1D6E9}$ , and consequently the velocity component sensitive to the stratification $\boldsymbol{a}\boldsymbol{\cdot }\tilde{\boldsymbol{S}}$ , oscillate at a frequency $\unicode[STIX]{x1D706}=\pm N\sin \unicode[STIX]{x1D709}$ . At the lowest order, the internal wave dispersion relation is retrieved.

Figure 4.

Schematic diagram defining the relevant angles to describe the relative position between initial stratification $\tilde{\boldsymbol{s}}_{0}$ and wave vector $\boldsymbol{{\mathcal{K}}}_{0}$ . The angle between $\tilde{\boldsymbol{s}}_{0}$ and the axis $z$ is $s$ defined in figure 3.

To further analyse the linear growth of the instability and to obtain quantitative growth rates, we perform the Floquet analysis of the system (3.2) (Bender & Orszag Reference Bender and Orszag1978). This can be done since the vectors $\boldsymbol{{\mathcal{K}}}$ , $\tilde{\boldsymbol{S}}$ and $\boldsymbol{e}_{s}$ all oscillate at the same frequency $\unicode[STIX]{x1D71B}$ , defined in (2.9a,b ). The linear operator on the right-hand side of (3.2) therefore oscillates with a period $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D71B}$ . Knowing the time evolution of $\boldsymbol{{\mathcal{K}}}$ , we solve the differential equations over the vector $(\boldsymbol{a},\unicode[STIX]{x1D703})$ from $t=0$ to $t=T=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D71B}$ . The initial condition is the identity matrix $\unicode[STIX]{x1D644}$ . The final value $(\boldsymbol{a}(T),\unicode[STIX]{x1D703}(T))$ for each initial condition corresponds to a monodromy matrix $\unicode[STIX]{x1D731}$ and the growth rate $\unicode[STIX]{x1D70E}$ of the instability is then related to its eigenvalues $\unicode[STIX]{x1D707}_{i}$ such that (Bender & Orszag Reference Bender and Orszag1978; Cébron, Vantieghem & Herreman Reference Cébron, Vantieghem and Herreman2014):

where $\tilde{\cdot }$ refers to a growth computed per tidal cycles. As the resonant wave vector is a priori unknown, the Floquet analysis is performed for different $\boldsymbol{{\mathcal{K}}}_{0}$ . This initial wave vector is parameterised by its norm ${\mathcal{K}}_{0}$ (which does not play any role in the inviscid limit), the angles $s$ , $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D6FC}$ represented in figure 4, such that:

Such a parameterisation merely comes from the expression of $\boldsymbol{{\mathcal{K}}}_{0}$ in the spherical coordinates $(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D709})$ with a polar axis $(M,\tilde{\boldsymbol{s}}_{0})$ (see figure 4). With $s$ used as a control parameter, resonance is found exploring the values of $\tilde{\unicode[STIX]{x1D70E}}$ in the $(\unicode[STIX]{x1D709},\unicode[STIX]{x1D6FC})$ space.

Figure 5.

Map of the growth rate computed via the Floquet theory $\tilde{\unicode[STIX]{x1D70E}}_{f}$ as a function of the ellipticity $\unicode[STIX]{x1D6FD}$ and the angle $\unicode[STIX]{x1D709}$ defined in figure 4 and (3.5) for (a)  $N=1.5$ and $s=90^{\circ }$ and (b) $N=4$ and $s=90^{\circ }$ . As will be shown hereafter, the growth rate is linearly growing with the ellipticity. We therefore normalise the growth rate by $\unicode[STIX]{x1D6FD}$ which allows us to identify the limits of the Floquet resonance tongue. The latter converges towards $N\sin \unicode[STIX]{x1D709}=\unicode[STIX]{x1D6FE}=1$ for $\unicode[STIX]{x1D6FD}\rightarrow 0$ showing that the resonance is subharmonic. Note that the higher growth rate area is always above the dashed line $N\sin \unicode[STIX]{x1D709}=1$ .

The range of angles $\unicode[STIX]{x1D709}$ to explore is non-trivial. We show in figure 5 the maximum (respectively to $\unicode[STIX]{x1D6FC}$ ) growth rate computed with the Floquet theory as a function of $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D6FD}$ . The non-zero growth rate area delimits Floquet tongues which stretch towards $N\sin \unicode[STIX]{x1D709}=\unicode[STIX]{x1D6FE}=1$ as $\unicode[STIX]{x1D6FD}$ 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 \unicode[STIX]{x1D709}=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 $\unicode[STIX]{x1D709}$ around $\arcsin (1/N)$ with a tolerance of order $\unicode[STIX]{x1D6FD}$ .

Note that the theoretical growth rate could also have been analytically computed via a multi-scale analysis where $\unicode[STIX]{x1D6FD}t$ is the slow time. This would have given the asymptotic resonant values of $\unicode[STIX]{x1D6FC}$ for $\unicode[STIX]{x1D6FD}\rightarrow 0$ . However, the complexity of the first-order operator respective to $\unicode[STIX]{x1D6FD}$ 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.

3.3 Comparison between direct numerical simulations and linear stability analysis

We perform direct numerical simulations (DNSs) setting $\unicode[STIX]{x1D6FD}=0.05$ , $N=1.5$ , $s=90^{\circ }$ with a Reynolds number $Re=10^{6.5}$ . A first try was initiated from a broadband white noise with $k/(2\unicode[STIX]{x03C0})\in [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 broadband noise to three intervals 1  $[5,10]$ , 2  $[10,14]$ and 3  $[14,20]$ which allows us 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 $\unicode[STIX]{x1D717}$ 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 us to calculate the viscous growth rate $\tilde{\unicode[STIX]{x1D70E}}_{v}$ in tidal units. The inviscid growth rate $\tilde{\unicode[STIX]{x1D70E}}_{th}$ (expressed in tidal cycles) is then obtained by subtracting the viscous damping of the growing mode, i.e. $\tilde{\unicode[STIX]{x1D70E}}_{th}=\tilde{\unicode[STIX]{x1D70E}}_{v}+2\unicode[STIX]{x03C0}k^{2}Re^{-1}$ with $k$ the wavenumber of the mode (as $\tilde{\unicode[STIX]{x1D70E}}_{th}$ is in tidal units, a $2\unicode[STIX]{x03C0}$ factor must be added to the damping rate).

Figure 6.

Growth of the instability for $\unicode[STIX]{x1D6FD}=5\times 10^{-2}$ , $N=1.5$ and $s=90^{\circ }$ initiated from broadband white noise with $k/(2\unicode[STIX]{x03C0})$ 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–linear coordinates (shown as dash-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.

(a) Map of the growth rate $\tilde{\unicode[STIX]{x1D70E}}_{f}$ computed with the Floquet analysis as a function of the angles $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D6FC}$ . The black dots correspond to the location of the growing modes observed in figure 6 for which the theoretical $\tilde{\unicode[STIX]{x1D70E}}_{f}$ growth rate is given. The dashed black lines highlight the marginal stability. (b) Comparison between systematic calculation of the maximum growth rate from the Floquet analysis and DNS results. The growth rates $\tilde{\unicode[STIX]{x1D70E}}_{f}$ computed with Floquet theory are all aligned on $\tilde{\unicode[STIX]{x1D70E}}_{f}=1.65\unicode[STIX]{x1D6FD}$ .

These results are then compared to the theoretical inviscid growth rate $\tilde{\unicode[STIX]{x1D70E}}_{f}$ given by the Floquet analysis of equations (3.2). The map of $\unicode[STIX]{x1D70E}_{f}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D709})$ is displayed in figure 7 where the location of the growing mode for each DNS is highlighted by a black dot associated with the corresponding $\tilde{\unicode[STIX]{x1D70E}}_{f}$ . We first note that, both in DNS and theory, the angle $\unicode[STIX]{x1D709}$ satisfies the condition $N\sin \unicode[STIX]{x1D709}=1$ with a tolerance smaller than $\unicode[STIX]{x1D6FD}$ , as it was expected from qualitative arguments developed in the preceding paragraph. In addition, the theoretical growth rate $\tilde{\unicode[STIX]{x1D70E}}_{f}$ is close to the growth rate $\tilde{\unicode[STIX]{x1D70E}}_{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.

3.4 Linearity with the ellipticity $\unicode[STIX]{x1D6FD}$

The amplitude of the periodic stretching and shearing responsible for the parametric excitation of internal waves is proportional to the ellipticity $\unicode[STIX]{x1D6FD}$ . Another way to validate the DNS and the linear stability analysis is to examine the consequent expected linearity in $\unicode[STIX]{x1D6FD}$ of the growth rate as in the case of tidally driven elliptical instability (Le Dizès Reference Le Dizès2000; Kerswell Reference Kerswell2002; Grannan et al. Reference Grannan, Favier, Le Bars and Aurnou2017). As shown in figure 7(b), the theoretical maximal growth rate $\tilde{\unicode[STIX]{x1D70E}}_{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^{\circ }$ . The growth rate computed from DNS is in addition very close to this theoretical line.

3.5 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^{\circ }$ and $60^{\circ }$ . The mode selection of the parametric instability changes with latitude. Unlike the mode selected in the case $N=1.5$ and $s=90^{\circ }$ where $\unicode[STIX]{x1D6FC}\simeq \unicode[STIX]{x03C0}/2$ (see figure 7), the mode selected at $s=60^{\circ }$ lies in the plane $(xMz)$ , i.e. $\unicode[STIX]{x1D6FC}\simeq \unicode[STIX]{x03C0}$ (see figure 4).

Figure 8.

(a) Theoretical and numerical growth rates with $\unicode[STIX]{x1D6FD}=5\times 10^{-2}$ and $N=1.5$ varying the stratification angle $s$ . (b) Corresponding map of the growth rate $\tilde{\unicode[STIX]{x1D70E}}_{f}$ computed with Floquet theory for the case $s=60^{\circ }$ . The black dashed lines highlight the marginal stability.

3.6 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\gg 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 \unicode[STIX]{x1D709}\sim 1/N$ .

Figure 9.

(a) Evolution of the maximum growth rate for $\unicode[STIX]{x1D6FD}=5\times 10^{-2}$ at three different stratification angles; a few DNSs have been carried out to confirm the agreement between theory and numerics. (b) Map of the growth rate computed from the Floquet analysis in the case $N=4$ with $s=90^{\circ }$ . The black dashed lines highlight the marginal stability.

4.1 Sustained instability and turbulence

Table 1.

Input parameters and measured statistical properties of the flow for each simulation. $k_{res}$ is the principal wavenumber of the resonant modes emerging during the growth phase. $u_{rms}$ is the root-mean-square (r.m.s.) velocity computed from the mean of the kinetic energy. The Froude number $Fr$ is computed as $u_{rms}/(N\unicode[STIX]{x1D706}_{res})$ . $\unicode[STIX]{x1D700}_{k}$ is the saturation dissipation rate defined as $-Re^{-1}\left\langle \unicode[STIX]{x2202}_{i}u_{j}\unicode[STIX]{x2202}_{i}u_{j}\right\rangle$ summed over the whole box. The output Reynolds $Re_{o}$ number and the buoyancy Reynolds number ${\mathcal{R}}$ are respectively defined in (4.3) and (4.4). When error bars 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.

Figure 10.

(a,b) Time evolution of the kinetic energy and the dissipation $\unicode[STIX]{x1D700}_{k}$ when $Re$ is varied while keeping constant $s=45^{\circ }$ and $N=1.5$ . The time is normalised by the theoretical growth rate $\unicode[STIX]{x1D70E}=\tilde{\unicode[STIX]{x1D70E}}_{f}$ . (c,d) Corresponding power spectra of the velocity $E(k)$ and of the buoyancy $E_{\unicode[STIX]{x1D703}}(k)$ ; they are averaged for $\tilde{\unicode[STIX]{x1D70E}}t\in [150,400]$ . The black solid line materialises a $k^{-3}$ power law. The dash-dotted spectra are computed during the growth phase to show that the energy is primarily injected in a narrow band of wavenumbers through the instability mechanism.

Figure 11.

(a,b) Power spectra of the velocity $E(k)$ and of the buoyancy $E_{\unicode[STIX]{x1D703}}(k)$ for $s\in \left\{0^{\circ },45^{\circ },90^{\circ }\right\}$ at $Re=10^{6.5}$ and $N=1.5$ . (c,d) The same quantities for $s=90^{\circ }$ with $N\in \left\{1.5,2,4\right\}$ . The Reynolds number is $10^{6.5}$ except for $N=4$ where it had to be increased to $10^{7}$ to observe sustained turbulence.

As in the linear stability analysis, the nonlinear 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 kinetic energy provided that $N$ and $s$ are kept constant. We note however 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 (Grannan et al. Reference Grannan, Le Bars, Cébron and Aurnou2014; Favier et al. Reference Favier, Grannan, Le Bars and Aurnou2015) or in the case of parametrically excited internal waves in a Faraday instability set-up (Benielli & Sommeria Reference Benielli and Sommeria1998).

To examine whether the saturation flow is influenced by stratification, we compute the Froude number based on the resonant wavelength and the r.m.s. velocity. The resonant wavelength $\unicode[STIX]{x1D706}_{res}=2\unicode[STIX]{x03C0}/k_{res}$ is an output parameter resulting from the mode selection during the growth phase (see table 1), and the r.m.s. value of the velocity in the saturated phase is computed from the saturation time average of the kinetic energy $\overline{E_{k}}$ such that $u_{rms}=\sqrt{2\overline{E_{k}}}$ . The Froude number $Fr$ is therefore defined as

The values of $u_{rms}$ and $Fr$ are all referenced in table 1. As $Fr=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_{\unicode[STIX]{x1D703}}(k)$ ) are shown in figures 10 and 11. They are computed as:

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 (Reference Brouzet2016) in the close context of turbulence driven by a forced internal wave attractor. It has also been measured at low Froude number by Rorai, Mininni & Pouquet (Reference Rorai, Mininni and Pouquet2015) 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 r.m.s. velocity and the resonant wavelength such that:

With this output Reynolds number, we can also compute the buoyancy Reynolds number ${\mathcal{R}}$ defined as (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007):

It compares a scale $\ell _{b}$ beyond which buoyancy effects are negligible and a scale $\ell _{v}$ beyond which viscosity affects the flow (Godoy-Diana, Chomaz & Billant Reference Godoy-Diana, Chomaz and Billant2004; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). For instance, in the context of classical stratified turbulence, it compares the so-called Ozmidov and the Kolmogorov length scales (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). The output Reynolds number is $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 ${\mathcal{R}}\gg 1$ regime (see for instance Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), Bartello & Tobias (Reference Bartello and Tobias2013), Maffioli, Brethouwer & Lindborg (Reference Maffioli, Brethouwer and Lindborg2016), Maffioli & Davidson (Reference Maffioli and Davidson2016)). 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 $\unicode[STIX]{x1D700}_{k}=-Re^{-1}\left\langle \unicode[STIX]{x2202}_{i}u_{j}\unicode[STIX]{x2202}_{i}u_{j}\right\rangle >0$ where $\left\langle \cdot \right\rangle$ 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 $\unicode[STIX]{x1D700}_{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 $\unicode[STIX]{x1D700}_{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.

Figure 12.

Mean value of the dissipation rate $\unicode[STIX]{x1D700}_{k}$ as a function of the Reynolds number for $N=1.5$ and $s=45^{\circ }$ . The error bars account for the standard deviation of $\unicode[STIX]{x1D700}_{k}$ during the saturation phase.

Figure 13.

Typical snapshots of the buoyancy field (a,c) and the $y$ component of the vorticity (b,d) in the saturated phase for a,b: $N=1.5$ , $s=45^{\circ }$ and $Re=10^{7}$ ; c,d: $N=4$ , $s=90^{\circ }$ and $Re=10^{7}$ .

4.2 Internal wave turbulence

In this paragraph, we aim at thoroughly characterising the structures generated by the nonlinear saturation of the initially unstable waves. In simulations and experiments of stratified turbulence, the emergence of 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. $\unicode[STIX]{x1D709}\rightarrow 0$ and $\unicode[STIX]{x1D714}\rightarrow 0$ ); three-dimensional motion comes through shear instability between those layers (see e.g. Billant & Chomaz (Reference Billant and Chomaz2001), Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007)). 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 (Brouzet et al. Reference Brouzet, Ermanyuk, Joubaud, Sibgatullin and Dauxois2016).

To determine which scenario is relevant here, we propose to map the energy in the same representation as Yarom & Sharon (Reference Yarom and Sharon2014), Brouzet et al. (Reference Brouzet, Ermanyuk, Joubaud, Sibgatullin and Dauxois2016) and Le Reun et al. (Reference Le Reun, Favier, Barker and Le Bars2017), i.e. to project the energy in the spectral space along the temporal frequency $\unicode[STIX]{x1D714}$ and the angle $\unicode[STIX]{x1D709}$ between the stratification direction $\tilde{\boldsymbol{S}}$ and the wave vector of a mode. This allows us 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 $\hat{\boldsymbol{u}}_{\boldsymbol{k}}(t)$ . They can be sorted by their angles $\unicode[STIX]{x1D709}$ to obtain the quantity $\hat{\boldsymbol{u}}(\unicode[STIX]{x1D709},t)$ defined as:

where $\unicode[STIX]{x1D6FC}$ is an azimuthal angle, $\unicode[STIX]{x1D709}$ is a polar angle respective to the stratification axis as defined in figure 4 and $k$ is the norm of the wave vector $\boldsymbol{k}$ . $\unicode[STIX]{x0394}\unicode[STIX]{x1D709}$ is a given tolerance to assume the angle $\unicode[STIX]{x1D709}$ of a mode is equal to $\unicode[STIX]{x1D709}_{0}$ . A time Fourier transform is applied to $\hat{\boldsymbol{u}}(\unicode[STIX]{x1D709},t)$ to finally get $\hat{\boldsymbol{u}}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D714})$ . With these definitions, we remind at this point that the dispersion relation of internal waves at a given $\boldsymbol{k}(k,\unicode[STIX]{x1D709},\unicode[STIX]{x1D6FC})$ is $\unicode[STIX]{x1D714}^{2}=N^{2}\sin ^{2}\unicode[STIX]{x1D709}$ .

Figure 14.

Map of the kinetic energy as a function of the frequency of the modes and the angle $\unicode[STIX]{x1D709}$ 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^{\circ }$ kept constant. The Fourier transforms are performed for $\tilde{\unicode[STIX]{x1D70E}}t\in [200,400]$ and the energy is normalised by its maximum value. The range of wavenumber $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\lesssim 20$ ) are included. The dotted horizontal line represents the frequency of the first excited modes, the plain line gives the dispersion relation of the internal waves $\unicode[STIX]{x1D714}=N\sin \unicode[STIX]{x1D709}$ and the dashed line locates the modes due to non-resonant nonlinear interactions between the tidal basic flow and the internal waves resulting in frequencies $2-N\sin \unicode[STIX]{x1D709}$ . The Snoopy code computes the time evolution of half the spectral space as the fields are real, the angle $\unicode[STIX]{x1D709}$ is therefore between $0$ and $\unicode[STIX]{x03C0}-\unicode[STIX]{x03C0}/4=3\unicode[STIX]{x03C0}/4$ .

Figure 15.

Map of the kinetic energy as a function of the frequency of the modes and the angle $\unicode[STIX]{x1D709}$ for $N=1.5$ , $2$ and $4$ with $s=90^{\circ }$ . The Reynolds number is $10^{6}$ except for $N=4$ for which it was increased to $Re=10^{7}$ to observe sustained turbulence. The Fourier transforms are performed for $\tilde{\unicode[STIX]{x1D70E}}t\in [150,400]$ and the energy is normalised by the maximum value. Again, secondary and mirroring locations of the energy corresponding to non-resonant and nonlinear interaction of the waves and the base flow can be noticed.

Figure 16.

Temporal spectrum resulting from the acquisition of 1024 local velocity signals of all three components. The plain vertical line marks the upper boundary of internal waves and the dashed vertical line highlights the excitation frequency. (a) Temporal spectrum for $N=1.5$ , $s=45^{\circ }$ and $Re=10^{7}$ . (b) Temporal spectrum for $N=4$ , $s=90^{\circ }$ and $Re=10^{7}$ . The inset shows the same amplitude compensated by $\unicode[STIX]{x1D714}^{2}$ with $\unicode[STIX]{x1D714}$ the frequency to highlight a $\unicode[STIX]{x1D714}^{-2}$ power law consistent with the high-frequency branch of the Garrett and Munk spectrum (Garrett & Munk Reference Garrett and Munk1979).

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 Yarom & Sharon (Reference Yarom and Sharon2014) and Le Reun et al. (Reference Le Reun, Favier, Barker and Le Bars2017) in rotating fluids, also Brouzet et al. (Reference Brouzet, Ermanyuk, Joubaud, Sibgatullin and Dauxois2016) in stratified fluids and to Aubourg & Mordant (Reference Aubourg and Mordant2015) in the case of gravity–capillary waves. It confirms that, in the saturation phase, the nonlinear 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 nonlinear 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 $\unicode[STIX]{x1D714}=2-N\sin \unicode[STIX]{x1D709}$ and are therefore associated with nonlinear and non-resonant interactions between the waves of frequency $\pm N\sin \unicode[STIX]{x1D709}$ and the base flow of frequency $\pm 2$ .

The filling of the dispersion relation depends however on the Brunt–Väisälä frequency (figure 15). When $N$ is increased, modes with frequency around or below $\unicode[STIX]{x1D6FE}$ seem to be more excited via nonlinear interactions than modes with frequency between $\unicode[STIX]{x1D6FE}$ 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 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 $\unicode[STIX]{x1D6FE}$ (see figure 16 a). These frequencies correspond to the lower branch (i.e. small $\unicode[STIX]{x1D709}$ ) 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 $\unicode[STIX]{x1D6FE}$ and $N$ , the energy contained in the higher frequencies ( $N\sin \unicode[STIX]{x1D709}>1$ ) follows a $\unicode[STIX]{x1D714}^{-2}$ power law. Such a trend is reminiscent of oceanographic measurements of the velocity which display a similar $\unicode[STIX]{x1D714}^{-2}$ power law in the range of frequencies above the tidal forcing and which is interpreted as a signature of internal wave turbulence (Garrett & Munk Reference Garrett and Munk1972, Reference Garrett and Munk1975, Reference Garrett and Munk1979; Levine Reference Levine2002). 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 nonlinearly interacting internal waves, although the nonlinearity is weak compared to the effects of the background stratification.

4.3 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 Billant & Chomaz (Reference Billant and Chomaz2001), Lindborg (Reference Lindborg2006) and Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007).

Figure 17.

Kinetic energy spectrum of the layerwise slow modes defined by $\unicode[STIX]{x1D709}=0$ and the rest of the flow for (a) $N=1.5$ , $Re=10^{6.5}$ and (b) $N=4$ , $Re=10^{7}$ at $s=90^{\circ }$ .

Figure 18.

Anisotropic velocity power spectra as functions of $k_{\bot }$ and $k_{\Vert }$ for (a) $N=1.5$ , $Re=10^{6.5}$ and $s=45^{\circ }$ and (b) $N=4$ , $Re=10^{7}$ and $s=90^{\circ }$ . On each plot, the black solid line shows a $k^{-3}$ power law for comparison.

First, figures 14 and 15 do not indicate any energy accumulation in the layerwise structures around $\unicode[STIX]{x1D709}=0$ and $\unicode[STIX]{x1D714}=0$ . To support this assertion, we show in figure 17 the kinetic energy spectrum for the layerwise modes at $\unicode[STIX]{x1D709}=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 wavenumber $\boldsymbol{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 to take over the whole dynamics in the absence of specific dissipative process and to strongly back-react on wave propagation (Barker & Lithwick Reference Barker and Lithwick2013; Le Reun et al. Reference Le Reun, Favier, Barker and Le Bars2017). This could be linked to the fundamental mathematical difference between layerwise modes in stratified turbulence and geostrophic vortices in rotating turbulence (Cambon Reference Cambon2001), which leads in particular to the absence of inverse cascade in the purely stratified case (Marino et al. Reference Marino, Mininni, Rosenberg and Pouquet2013; Herbert et al. Reference Herbert, Marino, Rosenberg and Pouquet2016).

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_{\bot }=k\sin \unicode[STIX]{x1D709}$ and $k_{\Vert }=k\cos \unicode[STIX]{x1D709}$ . The velocity power spectrum integrated over $k_{\Vert }$ , i.e. along the stratification axis, $E(k_{\bot })$ , follows a Kolmogorov-like power law $k_{\bot }^{-5/3}$ , while conversely $E(k_{\Vert })\propto k_{\Vert }^{-3}$ (Lindborg Reference Lindborg2006; Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). However, in the case of the internal wave turbulence presented above, both spectra $E(k_{\bot })$ and $E(k_{\Vert })$ follow the same power law close to $k_{\Vert ,\bot }^{-3}$ , as it can be seen in figure 18. The situation is even reversed as $E(k_{\bot })$ is slightly below $E(k_{\Vert })$ 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:

where $z_{s}$ is a linear coordinate along the stratification axis and $\boldsymbol{v}_{\bot }$ 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 Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), we compute for several Reynolds numbers, at $N=1.5$ and $s=45^{\circ }$ , 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 ${\mathcal{R}}$ being smaller than one, this is coherent with the remark of Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) 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.

Figure 19.

Probability density function of the Richardson number for $s=45^{\circ }$ and $N=1.5$ with $Re\in \left\{10^{6},10^{6.5},10^{7}\right\}$ . Probability density functions (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.

4.4 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 (Peltier & Caulfield Reference Peltier and Caulfield2003). Following the work of Lindborg & Brethouwer (Reference Lindborg and Brethouwer2008), Salehipour & Peltier (Reference Salehipour and Peltier2015) and Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016), we propose to quantify the mixing via a coefficient $\unicode[STIX]{x1D6E4}$ defined as:

where $\unicode[STIX]{x1D700}_{p}=-N^{2}(Re~Pr)^{-1}\left\langle (\unicode[STIX]{x1D735}\unicode[STIX]{x1D717})^{2}\right\rangle$ quantifies the diffusion of the buoyancy perturbations and $\unicode[STIX]{x1D700}_{k}$ is the kinetic energy dissipation defined in table 1. This coefficient $\unicode[STIX]{x1D6E4}$ was originally introduced to evaluate how turbulence induces an effective diapycnal diffusivity $K=\unicode[STIX]{x1D6E4}\unicode[STIX]{x1D700}_{k}/N^{2}$ (Osborn Reference Osborn1980; Salehipour et al. Reference Salehipour, Peltier, Whalen and MacKinnon2016). The expression given for $\unicode[STIX]{x1D700}_{p}$ can be retrieved by considering that it is a potential energy dissipation (Lindborg & Brethouwer Reference Lindborg and Brethouwer2008). In the derivation of the energy equations from (2.17), the buoyancy equation must be multiplied by $N^{2}\unicode[STIX]{x1D717}$ to obtain the same energy transfer from velocity to buoyancy, which finally yields the definition given earlier to $\unicode[STIX]{x1D700}_{p}$ .

Figure 20.

(a) Mixing coefficient $\unicode[STIX]{x1D6E4}$ as defined by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) as a function of the input Reynolds number. The limit value in the small Froude regime obtained by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) is given as a reference. (b) Evolution of the mixing coefficient with $N$ for $s=90^{\circ }$ and $Re=10^{6.5}$ except at $N=4$ where $Re=10^{7}$ .

In their study, Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) 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 $\unicode[STIX]{x1D6E4}$ converges towards $0.33$ . In figure 20(a) we display the mixing coefficient $\unicode[STIX]{x1D6E4}$ at $N=1.5$ as a function of the input Reynolds number $Re$ with the limit value found by Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016) given as a reference. Despite the fact that $Re$ is increased from $10^{6}$ up to $10^{7}$ we do not observe any variation of the mixing coefficient. Instead, $\unicode[STIX]{x1D6E4}$ remains constant around $1$ , well above the limit reference value. The evolution of the mixing coefficient $\unicode[STIX]{x1D6E4}$ with the Brunt–Väisälä frequency $N$ is also computed for $s=90^{\circ }$ . As it can be noticed in figure 20(b), again, $\unicode[STIX]{x1D6E4}$ 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 nonlinear interactions. A single wave of frequency $\unicode[STIX]{x1D714}$ and wave vector $\boldsymbol{k}$ , $\left\{\boldsymbol{u},\unicode[STIX]{x1D717},\unicode[STIX]{x1D6F1}\right\}=\left\{\boldsymbol{u}_{\boldsymbol{k}0},\unicode[STIX]{x1D717}_{\boldsymbol{k}0},\unicode[STIX]{x1D6F1}_{\boldsymbol{k}0}\right\}\text{e}^{\text{i}(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}-\unicode[STIX]{x1D714}t)}$ with $\unicode[STIX]{x1D714}^{2}=N^{2}\sin \unicode[STIX]{x1D709}$ , must obey the following linear inviscid set of equations:

where $\boldsymbol{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, $\left\langle \boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}\right\rangle =-k^{2}|\boldsymbol{u}_{\boldsymbol{k}0}|^{2}$ and $\left\langle \unicode[STIX]{x1D717}^{\ast }\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D717}\right\rangle =-k^{2}|\unicode[STIX]{x1D717}_{\boldsymbol{k}0}|^{2}$ . To compute $\unicode[STIX]{x1D6E4}$ , we need to relate $|\boldsymbol{u}_{\boldsymbol{k}0}|^{2}$ to $|\unicode[STIX]{x1D717}_{\boldsymbol{k}0}|^{2}$ which can be done for instance applying $\boldsymbol{k}\times (\boldsymbol{k}\times \boldsymbol{\cdot })$ to the velocity equation in (4.8). We then obtain the exact balance $|\boldsymbol{u}_{\boldsymbol{k}0}|^{2}=N^{2}|\unicode[STIX]{x1D717}_{\boldsymbol{k}0}|^{2}$ (which does not apply at $\unicode[STIX]{x1D714}=0$ ). Thus we find with the following simple scaling for the mixing coefficient:

where $Pr$ is the Prandtl or Schmidt number. We retrieve for our simulations at $Pr=1$ that $\unicode[STIX]{x1D6E4}=1$ . The numerical result $\unicode[STIX]{x1D6E4}=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.