# Mean flow generation via non-resonant interactions in two-dimensional forced stratified turbulence

Antonin Zoppi, Paul Billant, Leo Demaine

► To cite this version:

Antonin Zoppi, Paul Billant, Leo Demaine. Mean flow generation via non-resonant interactions in two-dimensional forced stratified turbulence. Journal of Fluid Mechanics, In press. hal-05361192

HAL Id: hal-05361192

https://hal.science/hal-05361192v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

# Mean flow generation via non-resonant interactions in two-dimensional forced stratified turbulence

Antonin Zoppi¹, Paul Billant ¹†, Leo Demaine¹

¹LadHyX, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

(Received xx; revised xx; accepted xx)

We perform numerical simulations of two-dimensional strongly stratified flows in a square periodic domain (y, z) forced by a steady mode with vorticity of the form $sin(k_{yf}y) sin(k_{zf}z)$, where ($k_{yf}, k_{zf}$) are fixed wavenumbers. It is shown that such deterministic forcing can lead to a transition to turbulence and the emergence of horizontal layers (so-called vertically sheared horizontal flows, VSHF) similarly as for random stochastic forcing. The flow characteristics are studied depending on the Froude and Reynolds numbers. Furthermore, the mechanisms of layers formation are disentangled. Triadic instabilities first lead to the growth of pairs of wavevectors that resonate with each of the four forced wavevectors. Quadratic interactions between these resonant modes and the forcing drive also the growth of several non-resonant modes at the same growth rate. Since the forcing comprises the wavevectors $±(k_{yf}, k_{zf})$ and their mirror symmetric with respect to the horizontal $±(k_{yf}, -k_{zf})$, there exist enslaved/bound modes with the same horizontal wavenumber and different vertical wavenumbers. Hence, the quadratic interactions between the latter modes force a second generation of modes among which some are VSHF. Their growth rate is twice the growth rate of the primary resonant modes. Such mechanism is similar to resonant quartets (Newell 1969; Smith & Waleffe 1999). When the forcing is restricted to only the two wavevectors $±(k_{yf}, k_{zf})$, the second generation of enslaved/bound modes have all a non-zero horizontal wavenumber. However, further quadratic interactions can force VSHF. Thus, horizontal layers also emerge but with a growth rate equal to the number of quadratic interactions times the growth rate of the primary instability.

26

**Key words:**

27

MSC Codes (Optional) Please enter your MSC Codes here

28

## 1. Introduction

29

30

31

Quasi-horizontal layers are ubiquitous in stably stratified flows (Riley & Lelong 2000; Caulfield 2021). In the case of numerical simulations of forced stratified turbulence in periodic domains, the layers take the form of the so-called “shear modes” or “vertically

<footnote>† Email address for correspondence: billant@ladhyx.polytechnique.fr</footnote>
sheared horizontal flows" (VSHF), i.e. horizontal flows purely uniform along the horizontal
32
and varying only along the vertical (Smith 2001; Smith & Waleffe 2002). These shear modes
33
usually grow slowly over a long time scale when using periodic boundary conditions both
34
in three-dimensional (Herring & Métais 1989; Laval et al. 2003; Waite & Bartello 2004;
35
Lindborg 2006; Brethouwer et al. 2007; Augier et al. 2015; Maffioli et al. 2016; Herbert
36
et al. 2016) and two-dimensional numerical simulations of stratified turbulence (Smith
37
2001; Remmel et al. 2014; Kumar et al. 2017; Fitzgerald & Farrell 2018a,b; Linares 2020).
38
Numerical simulations are often interrupted before shear modes become too dominant or, in
39
some cases, they are even artificially suppressed. Several mechanisms have been proposed to
40
explain their origin such as waves interactions (Smith 2001; Smith & Waleffe 2002; Remmel
41
et al. 2014), viscous effects (Lindborg 2006) or statistical instability of homogeneous stratified
42
turbulence (Fitzgerald & Farrell 2018a,b).
43

Similar mean flows varying in only one or two directions are observed in a variety of other
44
contexts such as rotating turbulence (Smith & Waleffe 1999), stratified-rotating turbulence
45
(Smith & Waleffe 2002) and turbulence on the β-plane (Rhines 1975; Smith & Waleffe 1999;
46
Galperin & Read 2019; Lemasquerier et al. 2021). All these fluids can sustain plane waves in
47
a infinite medium and the mean flow corresponds to a zero frequency wave or mode. However,
48
it is well-known that such slow wave or mode cannot be directly generated by resonant triad
49
interactions involving two other waves since the coefficient of interaction vanishes for zero
50
frequency (see Phillips (1968), McEwan et al. (1972) and Lelong & Riley (1991) for the case
51
of stratified fluids). Other plausible mechanisms such as near-resonant triads and resonant
52
quartets of waves have been therefore invoked to explain a transfer into a zero frequency wave
53
(Newell 1969; Smith & Waleffe 1999). In particular, it has been shown that a geostrophic
54
flow in rotating flows can be generated by a near-resonant instability of a single inertial
55
wave (Le Reun et al. 2020) or by a secondary quartetic instability of the triadic instability
56
in turbulence driven by a wave forcing (Brunet et al. 2020). Alternatively, the emergence
57
of a mean flow in geostrophic turbulence or two-dimensional stratified turbulence has been
58
explained in terms of a zonostrophic instability that breaks the statistical homogeneity of
59
turbulence (Srinivasan & Young 2012; Fitzgerald & Farrell 2018a,b; Farrell & Ioannou
60
2019). This framework shows indeed that, statistically, there is a positive feedback between
61
the fluctuating turbulent Reynolds stresses and the mean flow.
62

In contrast to the latter statistical approach which relies on a stochastic forcing of the
63
turbulence, Lemasquerier et al. (2021) have studied experimentally and numerically β-plane
64
flows generated by a steady forcing. Interestingly, they have shown that the resulting flows are
65
also turbulent and exhibit a zonal mean flow. This suggests that the mechanisms responsible
66
for the emergence of a mean flow could be understood in the simple case of a deterministic
67
forcing. With this in mind, we study herein two-dimensional stratified flows generated by a
68
steady forcing with a single wavenumber in each direction. This differs from previous studies
69
of two-dimensional stratified turbulence that were randomly forced within a given range of
70
small or large wavenumbers (Smith 2001; Kumar et al. 2017; Fitzgerald & Farrell 2018a;
71
Linares 2020). Despite this difference, we will see that the present steady forcing leads also
72
to a transition to a turbulent flow that exhibits a power law spectrum and shear modes. The
73
advantage of a deterministic forcing compared to a random forcing is that it will enable
74
us to disentangle the chronology of events leading to the transition to turbulence and the
75
emergence of shear modes.
76

Such evolution towards turbulence is not unexpected since it has been shown that the
77
breaking of a single two-dimensional internal wave can lead to stratified turbulence (Bouruet-
78
Aubertot et al. 1995, 1996; Benielli & Sommeria 1998; Staquet & Sommeria 2002; Koudella
79
& Staquet 2006; Dauxois et al. 2018). Furthermore, while a single monochromatic internal
80
plane progressive wave can not induce a mean flow, an inhomogeneous internal wave or the
81
interaction between several monochromatic waves may generate a mean flow (Sutherland 2010; Semin et al. 2016; Dauxois et al. 2018; Galmiche et al. 2000). The quasi-biennial oscillation of the wind in the stratosphere is one of the most famous example of mean flow induced by waves (Plumb & McEwan 1978; Wedi & Smolarkiewicz 2006; Semin et al. 2018). In the present configuration, we will see that the mean flow is induced by non-resonant interactions between different modes similarly as in Galmiche et al. (2000).

The paper is organized as follows. §2 presents the governing equations, the forcing with two pairs of wavevectors $\pm(k_{yf}, \pm k_{zf})$ and the numerical methods. A typical numerical simulation is first described in §3. The main effects of the control parameters (forcing amplitude, buoyancy frequency and viscosity) are then investigated in §4. In §5, we unravel the mechanisms responsible for the transition to turbulence and the generation of an horizontal mean flow. Numerical simulations for various wavenumbers of the forcing are presented and analysed in §6. Finally, §7 describes a simulation in which the forcing has been truncated to a single pair of wavevectors $\pm(k_{yf}, k_{zf})$. Turbulence and shear modes arise also but through a mechanism slightly different from the one evidenced in §5. The conclusions are drawn in §8.

## 2. Problem formulation and numerical methods

We use the two-dimensional incompressible Navier-Stokes equations under the Boussinesq approximation:

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla \left(\frac{p}{\rho_0}\right) + b\mathbf{e}_z + \nu\Delta\mathbf{u} + \mathbf{F}, \quad (2.1a) $$

$$ \frac{\partial b}{\partial t} + \mathbf{u} \cdot \nabla b = -N^2 w + \kappa \Delta b, \quad (2.1b) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2.1c) $$

where $\mathbf{u} = (v, w)$ is the velocity field in cartesian coordinates $\mathbf{x} = (y, z)$, $p$ the pressure, $b = -g\rho'/\rho_0$ the buoyancy with $g$ the gravity, $\rho'$ the density perturbation and $\rho_0$ the reference density. $\mathbf{e}_z$ is the vertical unit vector pointing upward, $\mathbf{F}$ the forcing, $N = \sqrt{-(g/\rho_0)\partial\bar{\rho}/\partial z}$ the Brunt-Väisälä frequency with $\bar{\rho}(z)$ the mean density profile. $v$ is the viscosity and $\kappa$ the diffusivity of the stratifying agent. The Schmidt number $S_c = v/\kappa$ will be always set to unity in the following.

The forcing is solenoidal with an horizontal vorticity $q = (\nabla \times \mathbf{F}) \cdot \mathbf{e}_x$ taken as a steady mode with a single wavenumber in each direction:

$$ q(y, z) = q_0 \sin(k_{yf} y) \sin(k_{zf} z), \quad (2.2) $$

where $q_0$ is the amplitude and $(k_{yf}, k_{zf})$ the wavenumbers of the forcing, that will be assumed positive without loss of generality. Alternatively, (2.2) can be written as the superposition of four wavevectors: $q = (q_0/4) (\exp(i\mathbf{k}_{f-} \cdot \mathbf{x}) + \exp(-i\mathbf{k}_{f-} \cdot \mathbf{x}) - \exp(i\mathbf{k}_{f+} \cdot \mathbf{x}) - \exp(-i\mathbf{k}_{f+} \cdot \mathbf{x}))$, where

$$ \mathbf{k}_{f+} = (k_{yf}, k_{zf}), \quad \mathbf{k}_{f-} = (k_{yf}, -k_{zf}). \quad (2.3) $$

The subscripts + and - will be convenient in the following to distinguish the different forced wavevectors.

If the velocity $\mathbf{u}$ and buoyancy $b$ are initially null, the response to (2.2) can be determined asymptotically for small amplitude $q_0$ as is done by Thorpe (1968) and Bouruet-Aubertot
et al. (1995) for a standing wave. At leading order, this yields:

$$124 \qquad b = q_0 \left[ \beta + e^{-ct} (a \sin \Omega t - \beta \cos \Omega t) \right] \cos(k_{yf} y) \sin(k_{zf} z) + O(q_0^2), \qquad (2.4a)$$

$$125 \qquad \psi = q_0 \left[ \gamma + e^{-ct} (\mu \sin \Omega t - \gamma \cos \Omega t) \right] \sin(k_{yf} y) \sin(k_{zf} z) + O(q_0^3), \qquad (2.4b)$$

where $\psi$ is a streamfunction such that $v = -\partial\psi/\partial z$, $w = \partial\psi/\partial y$ and

$$\Omega = \sqrt{\frac{N^2 k_{yf}^2}{k_f^2} - \frac{1}{4} \left(1 - \frac{1}{S_c}\right)^2 \nu^2 k_f^4}, \quad c = \frac{1}{2} \left(1 + \frac{1}{S_c}\right) \nu k_f^2, \quad \beta = \left[ k_{yf} \left( 1 + \frac{\nu^2}{N^2 S_c} \frac{k_f^6}{k_{yf}^2} \right) \right]^{-1}, \\ \gamma = -\frac{\nu k_f^2 \beta}{N^2 S_c k_{yf}}, \quad a = -\frac{\beta c}{\Omega}, \quad \mu = -\frac{\beta k_{yf}}{\Omega k_f^2} \left[ 1 - \frac{\nu^2 k_f^4}{N^2 k_{yf}^2 S_c} \left( 1 - \frac{1}{S_c} \right) \right],$$

where $k_f^2 = k_{yf}^2 + k_{zf}^2$. The parameters $\Omega$ and $c$ corresponds to the real and imaginary parts of the frequency of internal waves in presence of dissipative effects. Since $S_c = 1$, $\Omega$ recovers the dispersion relation in the inviscid limit. At the next order, the buoyancy is of the form $q_0^2 f(t) \sin 2k_{zf} z$ whereas $\psi$ is of the form $q_0^3 \sin(k_{yf} y) [g(t) \sin(3k_{zf} z) + h(t) \sin(k_{zf} z)]$, where $f$, $g$ and $h$ are functions of time (not given here). Hence, the horizontally averaged buoyancy is non-zero at order $q_0^2$ with a vertical wavenumber $2k_{zf}$. This comes from the advection of buoyancy in (2.1b). A similar feature occurs for a standing internal wave (Thorpe 1968; Bouruet-Aubertot et al. 1995). However, by noticing that (2.1) is invariant under the change $y \to -y$, $b \to b$, $\psi \to -\psi$, it can be deduced that the streamfunction should be an odd function of $y$ implying that the horizontally averaged streamfunction is zero at any order in $q_0$. In other words, there is no mean horizontal flow directly forced by (2.2).

The solution (2.4) shows that $b$ and $\psi$ oscillate at the frequency $\Omega$ but become steady as $t \to \infty$ since the time-dependent terms decays due to dissipation. The streamfunction is then proportional to $q_0 \gamma \propto q_0 \nu / N^2$ whereas the buoyancy scales as $q_0$. An example of the evolution of $b$ and $\psi$ will be presented at the end of this section. In any case, in the following numerical simulations, we will see that the flow will depart from (2.4) because the initial velocity field will not be strictly null but will consist of a divergence-free random noise of amplitude 0.01.

The Froude number $F_h$ and Reynolds number $Re$ are defined thanks to the characteristics of the forcing as follows

$$148 \qquad F_h = \frac{\sqrt{q_0}}{N}, \qquad Re = \frac{\sqrt{q_0 l_0^2}}{\nu}, \qquad (2.5)$$

with $l_0 = 2\pi/k_{yf}$ the horizontal wavelength of the forcing.

The equations (2.1a-2.1c) are solved by means of a pseudo spectral method with periodic boundary conditions (Deloncle 2007; Deloncle et al. 2008). Time stepping is performed by a fourth-order Runge-Kutta scheme. Aliasing is removed by using a truncation radius of 2/3.

The numerical domain has been taken square with sizes $L_y = L_z = 40$ with a resolution $N_y = N_z = 512$. In most of the simulations, the wavenumbers of the forcing have been set to $k_{yf} = k_{zf} = 4k_0$ where $k_0 = 2\pi/L_y$ is the lowest non-zero wavenumber. In other words, the forcing has four wavelengths in each direction. However, different forcing wavenumbers will be investigated in §6.

The Brunt-Väisälä frequency, the viscosity and the amplitude of the forcing have been varied in the ranges: $1.5 \le N \le 2.5$, $0.006 \le \nu \le 0.01$ and $0.083 \le q_0 \le 0.91$ (table 1). The corresponding Froude and Reynolds numbers are respectively small and large, in the ranges: $0.16 < F_h < 0.54$, $3580 < Re < 15710$.
<table><caption>Table 1: Summary of the parameters of the main simulations. Re and Fh are the Reynolds and Froude numbers, respectively, as defined in (2.5). v is the viscosity, N the Brunt-Väisälä frequency, q0 the amplitude of the forcing, (k_yf, k_zf) the wavenumbers of the forcing and Tf = 2π/Ω. kmax is the maximum wavenumber and η is the Kolmogorov length scale. Fht, Ret and Reb are, respectively, the turbulent Froude number, Reynolds number and buoyancy Reynolds number defined in (2.6).</caption><thead><tr><th>Re</th><th>F<sub>h</sub></th><th>ν</th><th>N</th><th>q<sub>0</sub></th><th>(k<sub>yf</sub>, k<sub>zf</sub>)</th><th>T<sub>f</sub></th><th>k<sub>max</sub>η</th><th>F<sub>ht</sub></th><th>Re<sub>t</sub></th><th>Re<sub>b</sub></th></tr></thead><tbody><tr><td>3590</td><td>0.16</td><td>0.008</td><td>1.75</td><td>0.083</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>5.2</td><td>0.015</td><td>331</td><td>0.08</td></tr><tr><td>7180</td><td>0.33</td><td>0.008</td><td>1.75</td><td>0.33</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>2.3</td><td>0.044</td><td>965</td><td>1.88</td></tr><tr><td>8840</td><td>0.40</td><td>0.008</td><td>1.75</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>1.8</td><td>0.060</td><td>1308</td><td>4.68</td></tr><tr><td>9520</td><td>0.43</td><td>0.008</td><td>1.75</td><td>0.58</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>1.7</td><td>0.068</td><td>1480</td><td>6.78</td></tr><tr><td>11920</td><td>0.54</td><td>0.008</td><td>1.75</td><td>0.91</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>1.5</td><td>0.076</td><td>1680</td><td>9.87</td></tr><tr><td>11790</td><td>0.40</td><td>0.006</td><td>1.75</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>1.5</td><td>0.058</td><td>1680</td><td>5.58</td></tr><tr><td>7070</td><td>0.40</td><td>0.01</td><td>1.75</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.08</td><td>2.1</td><td>0.061</td><td>1070</td><td>4.02</td></tr><tr><td>8840</td><td>0.47</td><td>0.008</td><td>1.5</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>5.92</td><td>1.9</td><td>0.066</td><td>1240</td><td>5.36</td></tr><tr><td>8840</td><td>0.35</td><td>0.008</td><td>2</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>4.44</td><td>2.0</td><td>0.046</td><td>1140</td><td>2.39</td></tr><tr><td>8840</td><td>0.31</td><td>0.008</td><td>2.25</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>3.95</td><td>1.9</td><td>0.043</td><td>1220</td><td>2.32</td></tr><tr><td>8840</td><td>0.28</td><td>0.008</td><td>2.5</td><td>0.50</td><td>(4k<sub>0</sub>, 4k<sub>0</sub>)</td><td>3.55</td><td>2.5</td><td>0.028</td><td>880</td><td>0.70</td></tr><tr><td>15710</td><td>0.40</td><td>0.008</td><td>1.75</td><td>0.50</td><td>(3k<sub>0</sub>, 3k<sub>0</sub>)</td><td>5.08</td><td>1.8</td><td>0.050</td><td>2000</td><td>4.96</td></tr><tr><td>8840</td><td>0.40</td><td>0.008</td><td>1.75</td><td>0.50</td><td>(4k<sub>0</sub>, 3k<sub>0</sub>)</td><td>4.49</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>8840</td><td>0.40</td><td>0.008</td><td>1.75</td><td>0.50</td><td>(4k<sub>0</sub>, 5k<sub>0</sub>)</td><td>5.75</td><td>2.2</td><td>0.048</td><td>1050</td><td>2.41</td></tr><tr><td>5660</td><td>0.40</td><td>0.008</td><td>1.75</td><td>0.50</td><td>(5k<sub>0</sub>, 4k<sub>0</sub>)</td><td>4.60</td><td>3.9</td><td>0.025</td><td>353</td><td>0.22</td></tr></tbody></table>

For comparison purpose, table 1 also indicates the turbulent Froude, Reynolds and buoyancy Reynolds numbers defined as follows:

$$F_{ht} = \frac{\epsilon_m^{1/3}}{Nl_0^{2/3}}, \quad Re_t = \frac{\epsilon_m^{1/3} l_0^{4/3}}{\nu\epsilon}, \quad Re_b = \frac{\epsilon_m}{\nu N^2}, \qquad (2.6)$$

where $\epsilon_m = \max(\epsilon)$ with $\epsilon$ the dissipation rate of kinetic energy averaged over the domain. These definitions are equivalent to those of Brethouwer et al. (2007) and Maffioli (2019) based on the root-mean square of the horizontal velocity fluctuations U instead of $l_0$ if the Taylor's relation $\epsilon_m \sim U^3/l_0$ holds (Lindborg 2006). As seen in table 1, all the simulations have a Froude number $F_{ht} = O(10^{-2})$ and most of them have a buoyancy Reynolds number larger than unity indicating that they correspond to the strongly stratified turbulent regime: $F_{ht} \ll 1$, $Re_b > 1$ (Brethouwer et al. 2007).

The numerical parameters always satisfy the criterion $k_{max}\eta > 1$ (table 1) required for direct numerical simulations, where $k_{max} = N_y k_0 / 2$ is the maximum wavenumber and $\eta = \nu^{3/4} / \epsilon_m^{1/4}$ the Kolmogorov scale (de Bruyn Kops & Riley 1998). The time step $\Delta t$ has been set to $10^{-3}$ to guarantee the Courant-Friedrichs-Lewy condition. The convergence of the simulations has been tested by running a simulation with the resolution $N_y = N_z = 1024$ for one set of parameters and by checking that it gives similar results as with $N_y = N_z = 512$.

In order to check and illustrate the asymptotic solution (2.4), a control simulation has been first performed with a zero initial velocity and buoyancy fields for $F_h = 0.4$ and $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$. Figure 1 shows the evolution of the amplitudes $\hat{b}$ and $\hat{v}$ (solid lines) of the Fourier modes of the buoyancy and the horizontal velocity component, i.e.

$$b = \sum_k \hat{b}(t, \mathbf{k}) e^{ik_y y + ik_z z}, \quad v = \sum_k \hat{v}(t, \mathbf{k}) e^{ik_y y + ik_z z}, \qquad (2.7)$$

that are non-zero in the wavenumber bands: $0 \le k_y \le 12k_0$, $0 \le k_z \le 12k_0$. For clarity, a
<figure><img src="image_2.png" /><figcaption>Figure 1: Evolution of the amplitudes of the modes in the wavenumber ranges $0 \le k_y \le 12k_0$, $0 \le k_z \le 12k_0$ of the (a) buoyancy $b$ and (b) horizontal velocity $v$ in a simulation without random noise at $t=0$ for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$ (solid lines). The dashed lines represent the asymptotic solution (2.4). The y-axis corresponds to natural logarithm. A moving average over a time interval $\Delta t = 5$ has been applied in order to filter the oscillations at frequency $\Omega$.</figcaption></figure>

moving average over a time interval $\Delta t = 5$ has been applied in order to smooth oscillations with period $T_f = 2\pi/\Omega = 5.08$. Note that the curves do not vanish at $t=0$ because of this moving average. The amplitude of the Fourier mode ($k_y, k_z$) = ($4k_0, 4k_0$) of the leading order asymptotic solution (2.4) is also represented by orange dashed lines and is in excellent agreement with the corresponding Fourier mode in the simulation (blue solid lines). As can be seen in figure 1(a), the dominant Fourier mode of $b$ converges quickly to the steady value $q_0\beta/4$ (see 2.4a). In contrast, according to (2.4b), the dominant mode of $v$ (figure 1b) decreases slowly by dissipation from the initial value $q_0\mu k_{zf}/4$, which is of order unity, towards the steady amplitude $q_0\gamma k_{zf}/4$. The latter value is small since it is proportional to the viscosity $\nu$ and is reached only after $t=800$. The slow modulations visible on this curve and other ones are due to the combined effect of the moving average and the fact that the fields are saved with the sampling period $\Delta t_s = 1$. The other Fourier modes that are generated at the next orders in the expansion with $q_0$ are all much smaller except the mode ($0, 8k_0$) for $b$ (red line in figure 1a) which corresponds to the second order term in (2.4a) ($q_0^2 f(t) \sin(2k_{zf}z)$). Figure 1(b) confirms that the forcing does not generate any mode for $v$ with zero horizontal wavenumber and non-zero vertical wavenumber, i.e. VSHF.

## 3. Description of a typical simulation

### 3.1. Global evolution

We first describe the main features of the simulation for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$, corresponding to the dimensional parameters $q_0 = 0.5$, $\nu = 0.008$, $N = 1.75$ (table 1). We remind that the velocity field is now initialized by a small random noise. Figure 2 displays the horizontal velocity component $v$ at different times. At early stages (figure 2(a)), the flow has exactly the same shape as the forced flow (2.4). At $t = 410$, this regular structure becomes modulated and exhibits inclined bands (figure 2(b)) indicating the emergence of modes different from the forced ones. Later on, at $t = 560$ (figure 2(c)), the velocity field becomes irregular at small scales although some periodicities are still present at large scales. A series of eight horizontal layers with alternate sign then arise (figure 2(d) at $t = 1100$). The vertical wavelength of this layered pattern is the same as the one of the forcing. However, by
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 2: Horizontal velocity field $v$ for $Re = 8840$, $F_h = 0.4$ and $k_{yf} = k_{zf} = 4k_0$ at different times : (a) $t = 1$, (b) $t = 410$, (c) $t = 560$, (d) $t = 1100$, (e) $t = 2100$, (f) $t = 4630$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 3: One-dimensional spectra of kinetic energy $E_k$ (red lines) and potential energy $E_p$ (blue lines) for $Re = 8840$, $F_h = 0.4$ and $k_{yf} = k_{zf} = 4k_0$ averaged over the time interval $2100 \le t \le 2200$. The solid and dashed lines show the horizontal and vertical spectra, respectively. The purple and green lines represent the $k^{-3}$ and $k^{-5/3}$ slopes, respectively. The vertical dashed line indicates the forced wavenumber $k_{yf} = k_{zf}$.</figcaption></figure>

212 $t = 2100$ (figure 2(e)) and $t = 4630$ (figure 2(f)), some layers have disappeared while others have intensified and thickened so that they are no longer regularly alternated.

214 Such flow structure and evolution strongly resemble two-dimensional stratified turbulence generated by random forcing (Smith 2001; Kumar et al. 2017; Fitzgerald & Farrell 2018a; Linares 2020). The turbulent character of the flow is further evidenced in figure 3 which displays one-dimensional horizontal and vertical spectra of kinetic and potential energies around $t = 2100$. As observed in the aforementioned studies of two-dimensional stratified turbulence, the spectra exhibit power-law behaviours with a slope closer to -3 (purple line) than -5/3 (green line). We can notice that the vertical kinetic energy spectrum around the forced wavenumber is larger than the horizontal kinetic energy spectrum reflecting the presence of the horizontal layers. At high wavenumbers, the vertical and horizontal spectra are closer indicating that velocity fluctuations at small scales are more isotropic.

224 In order to analyse the emergence and global evolution of the mean flow, the velocity component $v$ has been decomposed into its mean part $\langle v \rangle_y$ (corresponding to shear modes $\hat{v}(t, k_y = 0, k_z)$ in Fourier space) and its fluctuations $v'$:

$$v = \langle v \rangle_y + v', \quad (3.1)$$

228 where $\langle \cdot \rangle_a$ denotes the average along the direction $a$. Figure 4 first displays a space-time diagram of $\langle v \rangle_y(z, t)$. The layers become clearly visible around $t \approx 600$ and start to be strongly irregular after $t \approx 2000$. Figure 5 shows the root mean square (rms) of each of the three quantities $v$, $\langle v \rangle_y$, $v'$, where the rms of a quantity $X$ is defined as $X_{rms} = \sqrt{\langle X^2 \rangle_{y,z}}$. The rms of the vertical velocity $w$ and of the scaled buoyancy $b/N$ are also plotted with purple and orange lines, respectively. Until $t \approx 400$, $v_{rms}$ is small and follows the solution (2.4). However, an exponential increase of $v_{rms}$ occurs subsequently. During this exponential phase, $v_{rms}$ consists initially only in fluctuations, i.e. the mean flow $\langle v \rangle_{y-rms}$ is zero. Somewhat later, around $t \approx 500$, it starts to grow exponentially also. However, the inset in figure 5 shows that the growth rate of $\langle v \rangle_{y-rms}$ is twice the one of $v'_{rms}$, that has been estimated as
<figure><img src="image_10.png" /><figcaption>Figure 4: Space-time diagram of the horizontal mean flow $\langle v \rangle_y(z, t)$ for $Re = 8840$, $F_h = 0.4$ and $k_{yf} = k_{zf} = 4k_0$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 5: Temporal evolution of the root mean square of the total horizontal velocity $v$, the mean $\langle v \rangle_y$, the fluctuations $v'$, the vertical velocity $w$ and the scaled buoyancy $b/N$ for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$. The inset shows a close-up view of the curves of $\langle v \rangle_y$-rms, $v'_\text{rms}$ and $w_\text{rms}$ plotted in natural logarithm coordinate in the range $200 \le t \le 700$. The black straight lines represent the slopes of pure exponential growth at rates $\sigma$ and $2\sigma$.</figcaption></figure>

238 follows:

239

$$ \sigma = \frac{1}{t_2 - t_1} \ln \left( \frac{v'_{\text{rms}}(t_2)}{v'_{\text{rms}}(t_1)} \right), \qquad (3.2) $$

240 where $t_2$ and $t_1$ are two sufficiently distinct times chosen during the regime of pure exponential growth. After $t \ge 800$, $\langle v \rangle_y$-rms becomes larger than the fluctuating part $v'_{\text{rms}}$. The latter saturates after $t \approx 700$ whereas the rms of the mean flow saturates only later, for $t \ge 1200$. The rms of the vertical velocity $w_{\text{rms}}$ is of the same order as $v'_{\text{rms}}$. Hence, the fluctuating turbulent motions seem isotropic despite of the stable stratification. This is probably due to the fact that the forcing is isotropic and is present both in the horizontal and vertical momentum equations (2.1a).

241 242 243 244 245 246 In contrast, the rms of the scaled buoyancy $b_{\text{rms}}/N$ (orange line) exhibits large amplitude oscillations over a very long period around 2000–3000 time units. A first cycle of oscillation
<figure><img src="image_12.png" /><figcaption>Figure 6: Total buoyancy field $b_t = N^2z + b$ for $Re = 8840$, $F_h = 0.4$ and $k_{yf} = k_{zf} = 4k_0$ at different times: (a) $t = 1100$, (b) $t = 2100$.</figcaption></figure>

occurs between $t \approx 1000$ and $t \approx 3800$ while we can see half of the second cycle between $t \approx 3800$ and $t \approx 4700$. Such behaviour can be understood by examining the total buoyancy field $b_t = N^2z + b$ around a minimum and a maximum of $b_{\text{rms}}/N: t = 1100$ (figure 6(a)) and $t = 2100$ (figure 6(b)). The corresponding horizontal velocity field at these two instants were displayed in figures 2(d),2(e).

Near a minimum of $b_{\text{rms}}/N$, the iso-buoyancy lines are approximately equi-spaced along the vertical (figure 6(a)), i.e. there are no large buoyancy fluctuations $b$. In contrast, around a maximum of $b_{\text{rms}}/N$, the mean vertical gradient of $b_t$ is irregular: two regions of low gradient can be seen in the upper part of the domain (figure 6(b)). $b_{\text{rms}}/N$ is then high since the buoyancy $b$ exhibits fluctuations as large as the base buoyancy field $N^2z$. This corresponds to the regions where vortices (looking like a dipole pattern in terms of the velocity component $v$) are visible in the horizontal velocity field instead of horizontal layers (figure 2(e)). Hence, the weakening of the stable stratification in these regions allows for the existence of vortices associated with strong vertical motions. These regions are separated by thin layers where the stratification is intensified and where vortices are absent. Such a structure with weakly stratified regions of strong turbulence alternating with strongly stratified regions with weak turbulence is reminiscent of the mechanism for the emergence of density layers in stratified turbulence proposed by Phillips (1972) and Posmentier (1977).

Later on, $b_{\text{rms}}/N$ is again minimum at $t \approx 4000$ (figure 5). There are then layers almost everywhere (figure 4) although they are irregular. By $t \approx 4700$ (figure 5), a new peak of $b_{\text{rms}}/N$ is observed and again a region where strong vortices are present can be seen around $z = 15-20$ (figure 2(f)). The vertical gradient of total buoyancy $b_t$ is weakened in this region (not shown).

Finally, we stress that if the simulation is repeated for the same parameters but with a different occurrence of random noise, the overall characteristics of the different fields at a given time are identical, but they differ locally. In particular, the location of the layers varies, but their number at a given time is similar.
<figure><img src="image_13.png" /><figcaption>Figure 7: Evolution of the energy dissipations and production as defined in (3.3) for $Re = 8840$, $F_h = 0.4$ and $k_{yf} = k_{zf} = 4k_0$. The different lines correspond to the dissipation of kinetic energy $\epsilon$ (red), the dissipation of potential energy $\chi$ (blue), the energy provided by the forcing $P$ (yellow) and the total dissipation (purple).</figcaption></figure>

## 3.2. Energy analysis

Figure 7 shows the different terms of the total energy equation:

$$ \frac{d(E_k + E_p)}{dt} = -\epsilon - \chi + P, \qquad (3.3) $$

where $E_k = v_{rms}^2/2 + w_{rms}^2/2$, $E_p = b_{rms}^2/(2N^2)$ are the kinetic and potential energies, whereas $\epsilon = \nu \langle \nabla u : \nabla u \rangle_{y,z}$ is the dissipation of kinetic energy averaged over the domain, $\chi = \kappa \langle \nabla b^2 \rangle_{y,z}/N^2$ the dissipation of potential energy and $P = \langle F \cdot u \rangle_{y,z}$ the energy supplied by the forcing. The dissipation of kinetic energy (red line) and potential energy (blue line) are of the same order. The total dissipation $\epsilon + \chi$ (purple line) is almost superposed on the energy production rate $P$ (yellow line). This demonstrates that the flow is close to a statistically stationary state in which all the energy injected is dissipated. The rapid fluctuations of $P$ which are absent in $\epsilon$ and $\chi$ are due to conversion between kinetic and potential energies.

We can also distinguish two regimes in the saturation phase (figure 7). For $800 \le t \le 2500$, the energy production rate $P$ is high while beyond $t \approx 2500$, there is a decrease to a level approximately half the one before. Interestingly, these two time ranges correspond to the regimes of regularly thin alternated layers and of irregular thick layers (figures 2 and 4), respectively. In addition, figure 7 shows that the dissipation rates of potential and kinetic energies are of the same order as expected for strongly stratified turbulence (Maffioli et al. 2016). The mean mixing efficiency is $\Gamma = \chi/\epsilon \approx 0.7$ during the statistically steady regime. Such a large value could come from the fact that the scaled buoyancy $b/N$ is larger than the velocity perturbations (figure 5).

# 4. Parametric study

## 4.1. Effect of the forcing amplitude

The amplitude of the forcing $q_0$ has been varied in the range $0.083 < q_0 < 0.91$ by step of 0.083 or 0.041 for $N = 1.75$, $\nu = 0.008$ and $k_{yf} = k_{zf} = 4k_0$ (table 1). This corresponds to Froude numbers in the range $0.16 \le F_h \le 0.54$. A similar evolution as for the simulation for $q_0 = 0.5$ presented in the previous section, has been observed for each case but with some
<figure><img src="image_14.png" /><figcaption>Figure 8: Growth rates of the rms of the horizontal velocity $v$ (blue circles) and the mean $\langle v \rangle_y$ (green diamonds) as a function of the forcing amplitude $q_0$ for $N = 1.75$, $v = 0.008$ and $k_{yf} = k_{zf} = 4k_0$. The blue and green lines show the growth rate $\sigma$ and $2\sigma$, respectively, computed from (A 9) for each values of $q_0$.</figcaption></figure>

differences. In particular, the transition towards a layered turbulent flow occurs sooner as $q_0$ increases. Indeed, figure 8 shows that the initial growth rate of the rms of $v$ (blue circle) and $\langle v \rangle_y$ (green diamonds) increase linearly with $q_0$. Furthermore, the growth rate of $\langle v \rangle_{y-rms}$ is always twice the one of $v_{rms}$.

In addition, the velocity fluctuations tend to increase with $q_0$. Figure 9(a) shows the root-mean square of the total velocity (black line) averaged over a time interval once the flow is statistically stationary

$$ \bar{u}_{\text{rms}} = \frac{1}{T} \sum_{t=t_0}^{t_0+T} \sqrt{\langle v^2 \rangle_{y,z} + \langle w^2 \rangle_{y,z}}, \quad (4.1) $$

where $T = 1000$ and the starting time $t_0$ is chosen for each $q_0$ such that the rms quantities start to be reasonably statistically stationary after this time. However, it should be stressed that, for $q_0 < 0.4$, the rms can fluctuate largely over a long timescale so that $\bar{u}_{rms}$ depends somewhat on the particular time interval chosen.

We can see that $\bar{u}_{rms}$ first increases continuously and moderately with $q_0$ and, then, increases abruptly between $q_0 = 0.33$ and $q_0 = 0.37$ (figure 9a). For larger values of $q_0$, $\bar{u}_{rms}$ remains approximately constant although a small local decrease can be observed around $q_0 \approx 0.8$. The root-mean squares of each velocity component and of $\langle v \rangle_y$ and $v'$ averaged over the same time range are also plotted in figure 9(a). $\bar{v}_{rms}$ and $\langle \bar{v} \rangle_{y-rms}$ follow approximately the same trend as $\bar{u}_{rms}$ whereas the fluctuations $\bar{v}'_{rms}$ and $\bar{u}_{rms}$ vary a little more smoothly.

The horizontal velocity component $v$ is displayed in figure 10 for various values of $q_0$ at the instant $t = t_0 + T$, i.e. in the statistically steady regime. For low $q_0$ (figure 10(a)), a mean flow component is present but the layers tend to exhibit large and smooth variations that travel along the horizontal direction. Figure 10(a) corresponds to a late time $t = 4043$ and there are only four layers. However, eight layers are present at earlier times as in figure 2. For larger values of the forcing amplitude: $q_0 = 0.33$ (figure 10(b)) and $q_0 = 0.58$ (figure 10(c)), eight regular layers become clearly visible but they tend to be more and more distorted as $q_0$ increases. We can notice that these layers are located at the same vertical coordinates as for $q_0 = 0.5$ but their sign (i.e. direction of $v$) is opposite compared to figure 2(d). Finally, for $q_0 = 0.74$ (figure 10(d)), the layers are not regularly alternated and there exist regions where
<figure><img src="image_15.png" /><figcaption>Figure 9: (a) Averaged root-mean square of the various velocity components as a function of the forcing amplitude $q_0$ for $N = 1.75$ and $v = 0.008$. (b) Averaged root-mean square of the various velocity components divided by $Nl_0$ as a function of the Froude number for $v = 0.008$. The different lines correspond to: $\bar{u}_{rms}$ (black), $\bar{v}_{rms}$ (red), $\bar{w}_{rms}$ (yellow), $\langle v \rangle_{y-rms}$ (blue) and $\bar{v}'_{rms}$ (green). In (b), the solid and dashed lines represent the simulations for variable $q_0$ for $N = 1.75$ and for variable $N$ for $q_0 = 0.5$, respectively. The forcing wavenumbers are $k_{yf} = k_{zf} = 4k_0$.</figcaption></figure>

vortices are present (lower and upper parts) as for $q_0 = 0.5$ at late time (figure 2(f)). Apart from variations in the number and location of layers, which may depend on the particular time $t_0 + T$ chosen, the main observation in figure 10 is that small-scale turbulent fluctuations increase with $q_0$, which is observable at any time in the statistically steady regime, as we have already seen in figure 9a.

## 4.2. Effects of the Brunt-Väisälä frequency and the viscosity

The Brunt-Väisälä frequency has been varied in the range $1.5 \le N \le 2.5$ for $q_0 = 0.5$, $v = 0.008$ and $k_{yf} = k_{zf} = 4k_0$. The corresponding Froude number lie in the range $0.28 \le F_h \le 0.47$. Figure 11(a) shows the root-mean square of the various velocity components averaged over the time interval $[t_0, t_0 + T]$ as in figure 9. The average rms of the mean flow $\langle v \rangle_{y-rms}$ increases abruptly and then decreases as the Froude number $F_h$ increases. In contrast, the average rms of the horizontal velocity fluctuations $\langle v' \rangle_{rms}$ and of the vertical velocity $\langle w' \rangle_{rms}$ increase weakly with $F_h$. Hence, $\bar{u}_{rms}$ and $\bar{v}_{rms}$ evolves with $F_h$ similarly as $\langle v \rangle_{y-rms}$.

The open symbols with dashed lines in figure 9(b) shows the same average rms quantities as in figure 11(a) but now non-dimensionalized by $Nl_0$. We can see that the variations of these non-dimensional rms are then globally consistent with those when $q_0$ is varied (full symbols with solid lines).

Finally, figure 11(b) shows that the various average rms quantities are all almost independent of the Reynolds number when the viscosity is varied in the range $0.006 \le v \le 0.01$ for $q_0 = 0.5$, $N = 1.75$ and $k_{yf} = k_{zf} = 4k_0$. Altogether, the dependencies of $\bar{u}_{rms}$ (and other rms quantities) on the control parameters can be therefore written in the non-dimensional form

$$ \frac{\bar{u}_{rms}}{Nl_0} \approx \mathcal{F}(F_h). \qquad (4.2) $$
<figure><img src="image_16.png" /><figcaption>Figure 10: Snapshots of the horizontal velocity field $v$ in the statistically steady regime for various values of $q_0$: (a) $q_0 = 0.17$, $t = 4043$, (b) $q_0 = 0.33$, $t = 2471$, (c) $q_0 = 0.58$, $t = 2023$, (d) $q_0 = 0.74$, $t = 1968$ for $v = 0.008$, $N = 1.75$ and $k_{yf} = k_{zf} = 4k_0$. In each case, the instant shown corresponds to the time $t_0 + 1000$.</figcaption></figure>

# 5. Mechanisms of transition and of layer formation

## 5.1. Mean flow equation

As a preliminary step to understand how a mean flow can be generated, it is instructive to consider the mean flow equation

$$ \frac{\partial \langle v \rangle_y}{\partial t} = \mathcal{R} + \frac{1}{Re} \frac{\partial^2 \langle v \rangle_y}{\partial z^2}, \qquad (5.1) $$

that can be obtained by applying the mean operator $\langle \cdot \rangle_y$ to the momentum equation (2.1a) in the horizontal direction. The averages of all the other terms of (2.1a) vanish due to the periodic boundary conditions. The source of the mean flow is the Reynolds stress divergence given by

$$ \mathcal{R} = - \frac{\partial \langle v'w \rangle_y}{\partial z}. \qquad (5.2) $$

In order to see which modes can contribute to $\mathcal{R}$, it is interesting to consider a stream-
<figure><img src="image_17.png" /><figcaption>Figure 11: (a) Averaged root-mean square of the various velocity components as a function of the Froude number $F_h$ for variable $N$ for $q_0 = 0.5$ and $v = 0.008$. (b) Averaged root-mean square of the various velocity components as a function of the Reynolds number $Re$ for variable $v$ for $q_0 = 0.5$, $N = 1.75$ and $k_{yf} = k_{zf} = 4k_0$. The different lines correspond to: $\bar{u}_{rms}$ (black), $\bar{v}_{rms}$ (red), $\bar{w}_{rms}$ (yellow), $\langle \bar{v} \rangle_{y-rms}$ (blue) and $\bar{v}'_{rms}$ (green). The forcing wavenumbers are $k_{yf} = k_{zf} = 4k_0$.</figcaption></figure>

function $\psi'$ of the velocity fluctuations ($v'$, $w$) consisting in the sum of a pair of modes with wavevectors $\mathbf{k}_A$ and $\mathbf{k}_B$:

$$ \psi' = \psi_A(t)e^{i(k_y A y + k_z A z)} + \psi_B(t)e^{i(k_y B y + k_z B z)} + \text{c.c.}, \quad (5.3) $$

where $\psi_A$, $\psi_B$ are their amplitudes. This leads to the Reynolds stress divergence:

$$ \begin{aligned} \mathcal{R}(\mathbf{k}_A, \mathbf{k}_B) = & (k_{zA}k_{yB} + k_{zB}k_{yA}) \left[ \delta(k_{yA} - k_{yB})(k_{zA} - k_{zB}) \left( i\psi_A\psi_B^* e^{i(k_{zA}-k_{zB})z} + \text{c.c.} \right) \right. \\ & \left. - \delta(k_{yA} + k_{yB})(k_{zA} + k_{zB}) \left( i\psi_A\psi_B e^{i(k_{zA}+k_{zB})z} + \text{c.c.} \right) \right], \end{aligned} \quad (5.4) $$

where $\delta$ is the Dirac function. Here and hereinafter we use the shorthand notation $\mathcal{R}(\mathbf{k}_A, \mathbf{k}_B) = \mathcal{R}(\mathbf{k}_A, \mathbf{k}_B, \psi_A, \psi_B)$. As shown by Galmiche et al. (2000) for two internal waves, the expression (5.4) shows that the interaction of two modes can contribute to the Reynolds stress divergence, i.e. produce a mean flow, only if they have equal absolute horizontal wavenumbers $|k_{yA}| = |k_{yB}|$ and different absolute vertical wavenumbers $|k_{zA}| \neq |k_{zB}|$. This implies the well-known result that the self-interaction of a mode can not generate a mean flow. In addition, we recover the fact that the modes A, B and the mean flow can not form an unstable resonant triad since the temporal resonance condition $\omega_A = \omega_B$ and $|k_{yA}| = |k_{yB}|$ imply $|k_{zA}| = |k_{zB}|$ so that $\mathcal{R}(\mathbf{k}_A, \mathbf{k}_B) = 0$ (Phillips 1968; McEwan et al. 1972). These simple conditions are valid at any time and for arbitrary amplitude of the modes. They should be kept in mind in the next section where we will determine how a mean flow can emerge.

## 5.2. Instability mechanism

In order to understand the first instability leading to a transition towards turbulence, we now look at the dominant Fourier modes $\hat{v}$ of the horizontal velocity (see (2.7)). Figure 12 first shows the amplitudes of the Fourier modes that are above a threshold in the wavenumbers space $(k_y, k_z)$ at $t = 380$ for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$. The grid represents the discrete wavenumbers resolved in the numerical simulations, i.e. $(nk_0, pk_0)$ where $n$ and
<figure><img src="image_18.png" /><figcaption>Figure 12: Amplitudes of the Fourier modes $\hat{v}$ of the horizontal velocity at $t = 380$ for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$. The colour scale is given in natural logarithm. The grid represents the discrete wavenumbers resolved in the numerical simulations. The wavevectors labelled $k_{1+}$, $-k_{2-}$, $k_{f+}$ are the resonant wavevectors whereas the first non-resonant modes, $k_1^{(i)}$ and $k_2^{(i)}$ with $i = \{1, 2, 3\}$, that are linearly generated by the interactions with the forced wavevectors $k_{f-}$, $-k_{f-}$, $-k_{f+}$, are labelled in red and green, respectively. The wavevectors labelled in black correspond to a sample of modes that are non-linearly generated.</figcaption></figure>

$p$ are integers. Only the quadrant of positive wavenumbers is displayed since the three other quadrants are nearly symmetric and similar. Several modes with an even vertical wavenumber up to $20k_0$ are significant at that time. They are located at $k_z = 4nk_0$ or $k_z = 2 + 4nk_0$, where $n$ is a positive integer, if the horizontal wavenumber $k_y$ is an even or odd number of $k_0$, respectively. The labels of the modes in figure 12 will be explained later.

Figure 13(a) shows the time evolution of all the modes (grey lines) in the ranges: $0 \le k_{yf} \le 12k_0$, $0 \le k_{zf} \le 12k_0$. Like in figure 1, a moving average over a time interval $\Delta t = 5$ has been applied in order to smooth rapid oscillations. The three dominant modes are highlighted in color. The mode $(4k_0, 4k_0)$ (blue line) is the one which is directly forced. It initially follows the same evolution as in figure 1(b). The two other modes have wavenumbers $(k_y, k_z) = (3k_0, 6k_0)$ (green line) and $(k_0, 2k_0)$ (orange line). They are initially very small but they grow exponentially (figure 13(a)) and become dominant for $t \gtrsim 350$. The three other plots of figure 13 are identical except that different subdominant modes are plotted in color. They will be explained later below.

The origin of this instability can be understood from the resonant interaction theory of weak internal waves (McEwan 1971; McEwan et al. 1972; LeBlond & Mysak 1981; Staquet & Sommeria 2002; Dauxois et al. 2018). Following this theory, we consider the wavevector $k_{f+} = (k_{yf}, k_{zf})$ of the solution (2.4) at leading order and two linear internal waves with wavevectors $k_{1+}$ and $k_{2+}$ and frequency $\omega_{1+}$ and $\omega_{2+}$. They can form a resonant triad when
<figure><img src="image_19.png" /><figcaption>Figure 13: Evolution of the modes in the wavenumber ranges [0,12k₀]x[0,12k₀] for Fh = 0.4, Re = 8840 and k_yf = k_zf = 4k₀. The y axis corresponds to natural logarithm. In each of the four plots, all the modes are plotted in grey lines but different modes are highlighted in color. In (a), the three dominant modes k₁+, k₂-, kᵢ⁺ are displayed in color. In (b), the modes k₁⁽ᵢ⁾ and k₂⁽ᵢ⁾ with i = {1, 2, 3}, resulting from the interactions with kᵢ⁻, -kᵢ⁻, -kᵢ⁺ are plotted in colors. In (c), the three shear modes k_vf, k_v1, k_v2 and a sample of other modes, labelled in black in figure 12, that grow at a rate 2σ are plotted in color. In (d), the other shear modes are displayed in color. A moving average over a time interval Δt = 5 has been applied in order to filter the oscillations at frequency Ω.</figcaption></figure>

409 the following conditions are met:

$$410 \qquad k_{1+} = k_{f+} + k_{2+}, \qquad (5.5a)$$

$$4\frac{1}{2} \qquad \omega_{1+} = \omega_{f+} + \omega_{2+}, \qquad (5.5b)$$

413 where ωⱼs for j = {1, 2} and s = {+, −} are given by the dispersion relation:

$$414 \qquad \omega_{js} = \pm N \frac{k_{yjs}}{\|k_{js}\|}. \qquad (5.6)$$

415 The sign convention chosen in (5.5a) is arbitrary since the amplitudes of any pair of wavevectors **k** and −**k** are the complex conjugate of each other so that the fields are real.

416 In contrast with the classical theory of resonant wave interactions, the primary wave will be considered here to be steady, ωᵢ⁺ = 0, since the dominant term of the solution (2.4), i.e. the leading term in q₀ of **b**, becomes quickly steady (see figure 1a).

420 After some manipulations, the conditions of spatial and temporal resonances (5.5) can be
<figure><img src="image_20.png" /><figcaption>Figure 14: Components of the wavevector $\mathbf{k}_{1+}$ (red dashed line) satisfying the resonance conditions (5.5) for $\mathbf{k}_{f+} = (4k_0, 4k_0)$. The red solid line delimits the unstable region predicted by (A 9) for $q_0 = 0.5$, $\nu = 0.008$ and $N = 1.75$. The grid represents the discrete wavenumbers resolved in the numerical simulations. The colored vectors represent the unstable triad ($\mathbf{k}_{f+}, \mathbf{k}_{1+}, \mathbf{k}_{2+}$).</figcaption></figure>

421 reduced to:

$$k_{z1+} = \frac{k_{zf} k_{y1+}}{2k_{y1+} - k_{yf}}. \qquad (5.7)$$

423 Without loss of generality, we can assume $k_{y1+} > k_{yf}/2$ since the complementary interval 424 corresponds simply to the change ($\mathbf{k}_{1+}, \mathbf{k}_{2+}) \to (-\mathbf{k}_{2+}, -\mathbf{k}_{1+})$ in (5.5). The solution (5.7) is 425 shown by a red dashed line in figure 14. In addition, the solid red contour delimits the region 426 where the triad is unstable according to the growth rate (A 9) obtained from a multiscale 427 asymptotic analysis for small Froude number $F_h$ and large Reynolds number $Re$ (appendix 428 A). This analysis does not consider only exact resonance but also quasi-resonance, i.e. when 429 (5.5b) is not satisfied exactly but only with a detuning of order $O(F_h^2)$. This is the reason why 430 the instability band in figure 14 has a finite width around the exact resonance curve. We can 431 see that the sole discrete wavevector resolved in the numerical simulations which is within 432 the unstable region is $\mathbf{k}_{1+} = (3k_0, 6k_0)$. Actually, it coincides with the dashed red curve 433 meaning that the resonance is exact. The complementary wavevector of the resonant triad is 434 then $\mathbf{k}_{2+} = (-k_0, 2k_0)$ as represented in green in figure 14. The wavevector $\mathbf{k}'_{1+} = (6k_0, 3k_0)$ 435 (associated to $\mathbf{k}'_{2+} = (2k_0, -k_0)$) lies also on the curve of exact resonance but it is stable. 436 The stability analysis in appendix A shows indeed that a necessary condition for instability 437 is that $k_{y1+}$ and $k_{y2+}$ should be of opposite sign. Note that the width and upper $k_z$ limit of 438 the unstable domain depend on the parameters ($q_0, \nu, N$) (see A 9).

439 There is therefore a single unstable triad associated to the wavevector $\mathbf{k}_{f+}$:

$$\mathbf{k}_{f+} = (4k_0, 4k_0), \quad \mathbf{k}_{1+} = (3k_0, 6k_0), \quad \mathbf{k}_{2+} = (-k_0, 2k_0). \qquad (5.8a)$$

441 The equivalent triad for the wavevector $\mathbf{k}_{f-}$ (see 2.3) can be simply obtained by the sign 442 change $k_z \to -k_z$:

$$\mathbf{k}_{f-} = (4k_0, -4k_0), \quad \mathbf{k}_{1-} = (3k_0, -6k_0), \quad \mathbf{k}_{2-} = (-k_0, -2k_0). \qquad (5.8b)$$
<figure><img src="image_21.png" /><figcaption>Figure 15: Representation of the wavevectors of the triad with subscript + (solid lines) and of the non-resonant wavevectors (dashed lines) generated by the interaction with the other wavevectors of the base solution ($-\mathbf{k}_{f+}, \mathbf{k}_{f-}, -\mathbf{k}_{f-}$). (a) the non-resonant wavevectors</figcaption></figure>

$\mathbf{k}_{2+}^{(1)}, \mathbf{k}_{2+}^{(2)}, \mathbf{k}_{2+}^{(3)}$ are indicated by green dashed vectors. (b) the non-resonant wavevectors

$\mathbf{k}_{1+}^{(1)}, \mathbf{k}_{1+}^{(2)}, \mathbf{k}_{1+}^{(3)}$ are represented by red dashed vectors.

There exist two other triads with opposite wavevectors but they are linked to their counterparts above due to the relation between the amplitudes of $\mathbf{k}$ and $-\mathbf{k}$ wavevectors. Altogether, these symmetries imply that there is a pair of resonant wavevectors in each quadrant of the wavenumber space ($k_y, k_z$). This is the reason why we can focus our attention to only the upper right quadrant as done in figures 12 and 13

The wavevectors $\mathbf{k}_{1+}$ and $-\mathbf{k}_{2-}$ correspond exactly to the first two Fourier modes that grow exponentially in figure 13(a) (green and orange lines). Furthermore, the growth rate computed from the asymptotic formula (A 9) in appendix A is $\sigma_r = 0.0227$ in very good agreement with the slope of 0.023 of the orange and yellow lines in figure 13(a).

So far, we have considered only the interaction terms between the forced wavevectors and the perturbations that satisfy the resonance conditions (5.5). However, many other non-resonant modes are also generated as illustrated in figure 15 in the case of the triad with subscript +. The green dashed vectors in figure 15(a) represent the non-resonant wavevectors that are generated by the interactions of $\mathbf{k}_{1+}$ with $(-\mathbf{k}_{f+}, -\mathbf{k}_{f-}, \mathbf{k}_{f-})$:

$$ \begin{aligned} \mathbf{k}_{2+}^{(1)} &\equiv \mathbf{k}_{1+} + \mathbf{k}_{f+} = (7k_0, 10k_0), & \mathbf{k}_{2+}^{(2)} &\equiv \mathbf{k}_{1+} + \mathbf{k}_{f-} = (7k_0, 2k_0), \\ \mathbf{k}_{2+}^{(3)} &\equiv \mathbf{k}_{1+} - \mathbf{k}_{f-} = (-k_0, 10k_0), & \end{aligned} \quad (5.9) $$

respectively. Similarly, The red dashed vectors in figure 15(b) show the non-resonant wavevectors produced by the interactions between $\mathbf{k}_{2+}$ and $(-\mathbf{k}_{f+}, -\mathbf{k}_{f-}, \mathbf{k}_{f-})$:

$$ \begin{aligned} \mathbf{k}_{1+}^{(1)} &\equiv \mathbf{k}_{2+} - \mathbf{k}_{f+} = (-5k_0, -2k_0), & \mathbf{k}_{1+}^{(2)} &\equiv \mathbf{k}_{2+} - \mathbf{k}_{f-} = (-5k_0, 6k_0), \\ \mathbf{k}_{1+}^{(3)} &\equiv \mathbf{k}_{2+} + \mathbf{k}_{f-} = (3k_0, -2k_0), & \end{aligned} \quad (5.10) $$

respectively.

The same occurs for the triad with subscript – (5.8b) and those with opposite wavevectors. Hence, like for the resonant wavevectors $\mathbf{k}_{1+}$ and $\mathbf{k}_{2+}$, any of the non-resonant wavevectors (5.10) and (5.9) has symmetric counterparts in the three other quadrants. In the figure 12 showing only the upper right quadrant of the wavenumber space, the corresponding
471 wavevectors are labelled $\mathbf{k}_1^{(1)}$, $\mathbf{k}_1^{(2)}$, $\mathbf{k}_1^{(3)}$, $\mathbf{k}_2^{(1)}$, $\mathbf{k}_2^{(2)}$, $\mathbf{k}_2^{(3)}$ i.e. the sign and subscripts, + or -, are omitted in order to simplify the notation. The evolutions of their amplitudes are highlighted in figure 13(b). Even if their amplitudes are lower than those of the resonant modes, they all grow exponentially at rate $\sigma$ like the resonant modes. Indeed, since these non-resonant modes are driven by the interaction between the resonant modes $\mathbf{k}_{1s}$, $\mathbf{k}_{2s}$ with $s = \{+, -\}$ and the forced wavevectors, they should grow at the same rate as the resonant modes. In addition, their frequency is not the natural one but a forced one: i.e. the same as the wavevectors $\mathbf{k}_{1s}$ and $\mathbf{k}_{2s}$, i.e. $\pm\omega = \omega_{1s} = \omega_{2s}$ with $s = \{+, -\}$.

479 Among these eight modes generated by each triad, it is easy to see that some of them have the same horizontal wavenumber but different vertical wavenumbers. This is the condition for their quadratic interactions to generate shear modes, or equivalently, to contribute to the Reynolds stress divergence (5.4). For example, we see in figure 15(a) that the horizontal components of $\mathbf{k}_{2+} = (-k_0, 2k_0)$ and $\mathbf{k}_{2+}^{(3)} = (-k_0, 10k_0)$ are equal by construction. Hence, the quadratic interactions between $\mathbf{k}_{2+}$ and $\mathbf{k}_{2+}^{(3)}$, or between $\mathbf{k}_{2-}$ and $\mathbf{k}_{2+}^{(3)}$, will generate the shear modes

$$k_{vf} \equiv k_{2+}^{(3)} - k_{2+} = k_{f+} - k_{f-} = (0, 2k_{zf}) = (0, 8k_0), \qquad (5.11)$$

$$k_{v1} \equiv k_{2+}^{(3)} - k_{2-} = k_{1+} - k_{1-} = (0, 2k_{z1+}) = (0, 12k_0). \qquad (5.12)$$

489 The interactions between $\mathbf{k}_{2+}^{(1)}$ and $\mathbf{k}_{2+}^{(2)}$ in figure 15(a), or between $\mathbf{k}_{2+}^{(1)}$ and $\mathbf{k}_{2-}^{(2)}$, will produce the same shear modes. Similar considerations in figure 15(b) show that the interactions of $\mathbf{k}_{1+}$ with $\mathbf{k}_{1+}^{(3)}$, $\mathbf{k}_{1+}^{(1)}$ with $\mathbf{k}_{1+}^{(2)}$, $\mathbf{k}_{1-}$ with $\mathbf{k}_{1+}^{(3)}$ or $\mathbf{k}_{1+}^{(1)}$ with $\mathbf{k}_{1-}^{(2)}$ will produce the shear modes $\mathbf{k}_{vf}$ and

$$k_{v2} \equiv k_{1+}^{(3)} - k_{1-} = k_{2+} - k_{2-} = (0, 2k_{z2+}) = (0, 4k_0). \qquad (5.13)$$

494 Interchanging the subscripts + an - in the above relations leads simply to a sign change. The growth rate of these three shear modes is $2\sigma$ and their frequencies is $[-2\omega, 0, 2\omega]$ since they result from a quadratic interaction between modes growing at rate $\sigma$ with frequency $\pm\omega$. However, the forced response should be highest at zero frequency since this is the natural frequency of shear modes.

499 The generation of the three shear modes: $(\mathbf{k}_{vf}, \mathbf{k}_{v1}, \mathbf{k}_{v2})$, can be also interpreted in the framework of the quartet resonant mechanism (Newell 1969; Smith & Waleffe 1999). For example, the forced wavevector $\mathbf{k}_{f+}$ can be written in terms of its symmetric with respect to the horizontal axis, $\mathbf{k}_{f-}$, and the shear mode $\mathbf{k}_{vf}: \mathbf{k}_{f+} = \mathbf{k}_{f-} + \mathbf{k}_{vf}$. Hence, the resonance conditions (5.5) can be rewritten in terms of a quartet resonance:

$$\mathbf{k}_{1+} = \mathbf{k}_{f-} + \mathbf{k}_{vf} + \mathbf{k}_{2+}, \qquad (5.14a)$$

$$\omega_{1+} = \omega_{f-} + \omega_{vf} + \omega_{2+}, \qquad (5.14b)$$

507 where $\omega_{vf} = 0$. A similar substitution with the symmetric with respect to the horizontal could be applied to the two other wavevectors, i.e. $\mathbf{k}_{1+} = \mathbf{k}_{1-} + \mathbf{k}_{v1}$ or $\mathbf{k}_{2+} = \mathbf{k}_{2-} + \mathbf{k}_{v2}$. Then, (5.5) can be also rewritten in terms of a quartet resonance for these two cases.

508 The evolutions of the shear modes $\mathbf{k}_{v2} = (0, 4k_0)$, $\mathbf{k}_{vf} = (0, 8k_0)$ and $\mathbf{k}_{v1} = (0, 12k_0)$ are highlighted in color in figure 13(c) whereas all the other shear modes $(0, k_z)$ in the band $0 \le k_z \le 12k_0$ are displayed in color in figure 13(d). We can see that the three shear modes $\mathbf{k}_{v2}$, $\mathbf{k}_{vf}$ and $\mathbf{k}_{v1}$ grow exponentially at the rate $2\sigma$ as expected from the above reasoning. 512 513 514 They start to be amplified only after $t \ge 300$ clearly showing that they are a by-product of the triadic instabilities. This explains why the rms of $\langle v \rangle_y$ grows at a rate twice the rms of $v$ (inset of figure 5). In contrast, the other shear modes grow only later, $t \ge 400$, at a different growth
517 rate (figure 13(d)). Figure 13(c) shows also the evolutions of a sample of other significant
518 modes labelled in black in figure 12 that grow at rate $2\sigma$.

519 More generally, all the modes with wavenumbers such that $(k_y, k_z) = ((2n + 1)k_0, (4p + 520 2)k_0)$, where $n$ and $p$ are integers, grow with a growth rate $\sigma$. Indeed, these wavevectors 521 correspond to $q$ successive interactions of $k_{2+}$ with the forced wavevectors, giving $(k_y, k_z) = 522 k_{2+} + 4qk_0(\pm 1, \pm 1)$. In contrast, the modes such that $(k_y = 2nk_0, k_z = 4pk_0)$ grow at rate $2\sigma$. 523 These modes are generated by one interaction between $k_{1+}$ and $k_{2+}$ followed by $q$ successive 524 interactions with the forced wavevectors, yielding $(k_y, k_z) = k_{1+} + k_{2+} + 4qk_0(\pm 1, \pm 1)$. 525

526 The other modes, i.e. those with an odd vertical wavenumber and those with $(k_y, k_z) = 527 ((2n + 1)k_0, 4pk_0)$ or $(k_y, k_z) = (2nk_0, (4p + 2)k_0)$ grow only later at a different growth 528 rate, probably due to a secondary triadic instability of modes amplified by the first triadic 529 instability. The shear modes plotted in figure 13(d) belong to this family.

530 In summary, triadic instabilities with the four forced wavevectors $\pm k_{f+}$, $\pm k_{f-}$ of the 531 base solution (2.4) lead first to the exponential amplification of the internal waves with 532 wave-vector $\pm k_{1+}$, $\pm k_{1-}$ and $\pm k_{2+}$, $\pm k_{2-}$. The quadratic interactions of these resonant waves 533 with the forced wavevectors drive also the growth at the same rate of several non-resonant 534 modes: $\pm k_{j+}^{(i)}$, $\pm k_{j-}^{(i)}$ for $i = \{1, 2, 3\}$ and $j = \{1, 2\}$. These are not internal waves since their 535 frequency is the same as the resonant waves. Then, the quadratic interactions between all 536 these modes further generate many other modes growing at a double growth rate. Among 537 them, we have shown that there are the shear modes: $k_{vf} = (0, 2k_{zf})$, $k_{v1} = (0, 2k_{z1+})$ and 538 $k_{v2} = (0, 2k_{z2+})$. 539

540 We emphasize that this scenario is observed whatever the amplitude of the forcing $q_0$ in 541 the range investigated. In particular, the asymptotic growth rate (A 9) (blue line in figure 8) is 542 in excellent agreement with the growth rate of $v_{rms}$ for all the values of $q_0$ investigated (blue 543 circles). Similarly, the growth rate of $\langle v \rangle_{y-rms}$ (green diamonds in figure 8) corresponds 544 always to twice the asymptotic growth rate (A 9) (green line). 545

543 5.3. Reynolds stress

544 In the previous section, we have shown that the pairs of wavevectors $(k_{1s}, k_{1r}^{(3)})$, $(k_{2s}, k_{2r}^{(3)})$, 545 $(k_{1s}^{(1)}, k_{1r}^{(2)})$ or $(k_{2s}^{(1)}, k_{2r}^{(2)})$, where $s, r$ stand for any pair among the labels $\{+, -\}$, have the 546 same horizontal wavenumber but different vertical wavenumbers. Therefore, their interaction 547 can contribute to the Reynolds stress divergence (5.2) and be responsible for the generation 548 of the mean flow. Among these pairs of modes, the first two are expected to be dominant 549 initially since they each involve one resonant wavevector, $k_{1s}$ or $k_{2s}$, whose amplitudes are 550 larger than those of the non-resonant wavevectors $k_{1s}^{(i)}$, $k_{2s}^{(i)}$, for $s = \{+, -\}$ and $i = \{1, 2, 3\}$. 551 Hence, we have computed the Reynolds stress divergence by taking into account only the 552 former pairs of modes. More precisely, the total Reynolds stress $\mathcal{R}$ is approximated as the 553 sum of their individual contributions obtained with (5.4):

$$
\begin{equation}
\begin{split}
\mathcal{R} ={}& \mathcal{R} \left( k_{1+}, k_{1+}^{(3)} \right) + \mathcal{R} \left( k_{1+}, k_{1-}^{(3)} \right) + \mathcal{R} \left( k_{1-}, k_{1-}^{(3)} \right) + \mathcal{R} \left( k_{1-}, k_{1+}^{(3)} \right) \\
& + \mathcal{R} \left( k_{2+}, k_{2+}^{(3)} \right) + \mathcal{R} \left( k_{2+}, k_{2-}^{(3)} \right) + \mathcal{R} \left( k_{2-}, k_{2-}^{(3)} \right) + \mathcal{R} \left( k_{2-}, k_{2+}^{(3)} \right).
\end{split}
\tag{5.15}
\end{equation}
$$

557 Using the mode amplitudes extracted from the simulations, the value of $\mathcal{R}$ predicted by (5.15) 558 has been then compared to the Reynolds stress divergence computed from the full velocity 559 fields with (5.2).

560 An excellent agreement has been observed in the period of pure exponential growth of the 561 resonant and non-resonant modes, i.e. for $t \le 420$ (see figure 16a for $t = 380$), confirming 562 the reasoning carried out in §5.2. However, later on (figure 16b), the approximation (5.15) no
<figure><img src="image_22.png" /><figcaption>Figure 16: Reynolds stress divergence $\mathcal{R}$ at times (a) $t = 380$ and (b) $t = 460$ for $Re = 8840$, $F_h = 0.4$ and $k_{yf} = k_{zf} = 4k_0$. The blue line represents $\mathcal{R}$ determined from the formula (5.2) while the orange line corresponds to $\mathcal{R}$ computed with the approximation (5.15).</figcaption></figure>

563 longer agrees with the Reynolds stress divergence computed directly from the DNS because
564 many other modes than considered in (5.15) contribute significantly to $\mathcal{R}$.

6. Effect of the wavenumbers of the forcing

566 So far, the wavenumbers of the forcing were fixed to $(k_{yf}, k_{zf}) = (4k_0, 4k_0)$. Here,
567 we now present four simulations with the forcing (2.2) with different forcing wavenum-
568 bers: $(k_{yf}, k_{zf}) = (3k_0, 3k_0)$, $(k_{yf}, k_{zf}) = (4k_0, 3k_0)$, $(k_{yf}, k_{zf}) = (5k_0, 4k_0)$ and
569 $(k_{yf}, k_{zf}) = (4k_0, 5k_0)$ for $q_0 = 0.5$, $N = 1.75$ and $\nu = 0.008$. These four examples
570 are typical of the different behaviours that have been observed in a series of simulations
571 performed systematically in the wavenumber ranges: $2k_0 \le k_{yf} \le 5k_0$, $2k_0 \le k_{zf} \le 5k_0$.

572 Figure 17 first displays the wavevector $\mathbf{k}_{1+}$ (red dashed line) satisfying the resonance 573 condition (5.5) for each forcing wavevector $\mathbf{k}_{f+}$. As in figure 14, the red solid contour 574 delimits the unstable region according to (A 9). The unstable discrete wavenumbers are 575 highlighted by filled circles whose color indicates the growth rate value (colorbar). For 576 $\mathbf{k}_{f+} = (3k_0, 3k_0)$ (figure 17a), the wavevector $\mathbf{k}_{1+} = (2k_0, 6k_0)$ (red symbol) satisfies 577 exactly the resonance condition and is unstable. The wavevector $\mathbf{k}_{1+} = (2k_0, 7k_0)$ (green 578 symbol) is also weakly unstable but it is only quasi-resonant, i.e. it is inside the unstable 579 band without lying on the curve of exact resonance. For $\mathbf{k}_{f+} = (4k_0, 3k_0)$ (figure 17b), none 580 of the discrete wavenumbers are unstable.

581 For $\mathbf{k}_{f+} = (4k_0, 5k_0)$ (figure 17c), none of the discrete wavenumbers are exactly resonant 582 but there exists two wave-vectors that are quasi-resonant: $\mathbf{k}_{1+} = (3k_0, 7k_0)$ (dark blue symbol) 583 and $\mathbf{k}_{1+} = (3k_0, 8k_0)$ (blue symbol) although the former is only weakly unstable.

584 Finally, for $\mathbf{k}_{f+} = (5k_0, 4k_0)$ (figure 17d), there is one unstable exact resonance for 585 $\mathbf{k}_{1+} = (3k_0, 12k_0)$ and four unstable quasi-resonance for wavevectors with the horizontal 586 wavenumber $k_{y1+} = 3k_0$ and various vertical wavenumbers in the range $10k_0 \le k_{z1+} \le 14k_0$. 587 Among all these unstable modes, $\mathbf{k}_{1+} = (3k_0, 11k_0)$ is dominant whereas $\mathbf{k}_{1+} = (3k_0, 12k_0)$ 588 is slightly subdominant despite it is an exact resonance. The three other wavevectors are 589 somewhat less unstable.

590 We emphasize that the results in figure 17 refer to the amplitude of the forcing $q_0 = 0.5$.
591 If $q_0$ were varied, quasi-resonant wave-vectors could appear or disappear, because the width
<figure><img src="image_23.png" /><figcaption>Figure 17: Wavevector $k_{1+}$ satisfying the resonance conditions (5.5) (red dashed line) for (a) $k_{f+} = (3k_0, 3k_0)$, (b) $k_{f+} = (4k_0, 3k_0)$, (c) $k_{f+} = (4k_0, 5k_0)$ and (d) $k_{f+} = (5k_0, 4k_0)$ for $q_0 = 0.5$, $v = 0.008$, $N = 1.75$. The red solid line delimits the unstable region according to (A 9). The blue vector represents the forced wavevector $k_{f+}$. The grid represents the discrete wavenumbers resolved in the numerical simulations. Those that are unstable are highlighted by a circle whose color indicates the growth rate level (colorbar).</figcaption></figure>

of the unstable region around the exact resonance curve increases with the amplitude of the forced solution (2.4) (see (A 9)). Thus, quasi-resonant wave-vectors should also appear in figure 17b if $q_0$ were increased beyond a critical value. Conversely, the quasi-resonant wave-vectors in figure 17c should disappear if $q_0$ were sufficiently decreased. In contrast, the exact resonances are expected to be less sensitive to variations in $q_0$.

Figure 18 displays space-time diagram of the mean horizontal flow $\langle v \rangle_y(z, t)$ for the three unstable cases. When $(k_{yf}, k_{zf}) = (4k_0, 3k_0)$ (not shown), the flow simply evolves according to the solution (2.4) and there are no transition to turbulence and emergence of mean horizontal flows in agreement with the absence of triadic instability (figure 17b).

We see that for $(k_{yf}, k_{zf}) = (3k_0, 3k_0)$ (figure 18(a)) and $(k_{yf}, k_{zf}) = (4k_0, 5k_0)$ (figure 18(b)), i.e. when a single pair of resonant modes is dominant, the magnitude of the mean flow and the layers dynamics resemble the one for $(k_{yf}, k_{zf}) = (4k_0, 4k_0)$ (figure 4). In contrast, for $(k_{yf}, k_{zf}) = (5k_0, 4k_0)$, the mean flow is much less prominent and the layers tend to fade away after their onset. For the latter case, figure 19 further shows that the root-mean square of the mean horizontal flow is smaller than the one of the fluctuation $v'_{rms}$ in contrast to figure 5.

The time evolution of all the modes have been also investigated for each case. For $(k_{yf}, k_{zf}) = (3k_0, 3k_0)$ and $(k_{yf}, k_{zf}) = (4k_0, 5k_0)$, i.e. when there is a single pair of dominant resonant wavevectors, the dynamics is very similar to the case $(k_{yf}, k_{zf}) = (4k_0, 4k_0)$ and the emerging shear modes have been observed to be $k_{vf} = (0, 2k_{zf})$, $k_{v1} = (0, 2k_{z1+})$ and $k_{v2} = (0, 2k_{z2+})$ for both cases in agreement with §5.2. These two cases are therefore not shown. In contrast, for $(k_{yf}, k_{zf}) = (5k_0, 4k_0)$, the dynamics differ because there exist two pairs of dominant resonant modes: $k_{1+} = (3k_0, 11k_0)$, $k_{2+} = (-2k_0, 7k_0)$ and $k_{1+} = (3k_0, 12k_0)$, $k_{2+} = (-2k_0, 8k_0)$. As seen in figure 20(a), these four modes grow approximately at the same rate in agreement with the asymptotic stability analysis (figure 17d). According to (5.10) and (5.9), the six non-resonant modes generated by the first triad
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>Figure 18: Space-time diagrams of the horizontal mean flow $(v)_y(z, t)$ for (a) $(k_{yf}, k_{zf}) = (3k_0, 3k_0)$, (b) $(k_{yf}, k_{zf}) = (4k_0, 5k_0)$ and (c) $(k_{yf}, k_{zf}) = (5k_0, 4k_0)$ for $q_0 = 0.5$, $v = 0.008$, $N = 1.75$.</figcaption></figure>

618 are

619

$$ \mathbf{k}_{1+}^{(1)} = (-7k_0, 3k_0), \quad \mathbf{k}_{1+}^{(2)} = (-7k_0, 11k_0), \quad \mathbf{k}_{1+}^{(3)} = (3k_0, 3k_0), $$

620

$$ \mathbf{k}_{2+}^{(1)} = (8k_0, 15k_0), \quad \mathbf{k}_{2+}^{(2)} = (8k_0, 7k_0), \quad \mathbf{k}_{2+}^{(3)} = (-2k_0, 15k_0), \quad (6.1) $$

622

while for the second one, they are

623

$$ \mathbf{k}_{1+}^{(1)} = (-7k_0, 4k_0), \quad \mathbf{k}_{1+}^{(2)} = (-7k_0, 12k_0), \quad \mathbf{k}_{1+}^{(3)} = (3k_0, 4k_0), $$

624

$$ \mathbf{k}_{2+}^{(1)} = (8k_0, 16k_0), \quad \mathbf{k}_{2+}^{(2)} = (8k_0, 8k_0), \quad \mathbf{k}_{2+}^{(3)} = (-2k_0, 16k_0). \quad (6.2) $$

626

The shear modes that are generated independently by each of the two triads are: $\mathbf{k}_{vf} = (0, 8k_0)$, $\mathbf{k}_{v1} = (0, 22k_0)$, $\mathbf{k}_{v2} = (0, 14k_0)$ and $\mathbf{k}_{vf} = (0, 8k_0)$, $\mathbf{k}_{v1} = (0, 16k_0)$, $\mathbf{k}_{v2} = (0, 24k_0)$, respectively. They are highlighted in color in figure 20(a).

627

628

629

However, the mutual interactions between the resonant and non-resonant modes of the two triads offer the possibility to generate several additional shear modes since they have the same

630
<figure><img src="image_27.png" /><figcaption>Figure 21: Same as figure 4 but with the forcing (7.1).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 22: Same as figure 5 but with the forcing (7.1).</figcaption></figure>

steady regime turns out to be lower for $(k_{yf}, k_{zf}) = (5k_0, 4k_0)$ than when there is a single dominant triadic instability. The understanding of this paradoxical effect is left for future studies.

## 7. Forcing with a single pair of wavenumber

In the previous sections, we have seen that a vertically sheared horizontal mean flow can be generated by the quadratic interactions between two modes that are driven by two different triads that are symmetric with respect to the horizontal. This is possible since the forcing (2.2) contains symmetric wavevectors: $\pm \mathbf{k}_{f+}$ and $\pm \mathbf{k}_{f-}$.

The next question is: what happens if the forcing does not contain wavevectors symmetric with respect to the horizontal? The previously observed mechanism is a priori no longer possible and the question arises as to whether a mean flow can still occur. In order to answer this question, we have truncated the forcing to

$$q(y, z) = -\frac{q_0}{2} \cos(k_{yf} y + k_{zf} z), \qquad (7.1)$$

so that the forced wavevectors are now only $\pm \mathbf{k}_{f+}$.

Figure 21 shows a space-time diagram of a simulation with such forcing for the same
657 parameters as the simulation described in §3: $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$.
658 We can see that layers do arise around the same time as before but their dynamics is more
659 erratic and their distribution on the vertical is less uniform than in figure 4. Furthermore,
660 they wander along the vertical and can disappear more abruptly.

661 As seen in figure 22, the rms of the horizontal mean flow $\langle v \rangle_y$ is again largely dominant
662 relative to the rms of the fluctuations $v'$. However, it exhibits larger fluctuations than in figure
663 5. We can also remark that the average levels of the rms of $v'$ and $w$ are approximately half
664 the ones in figure 5. Interestingly, the inset in figure 22 shows that $\langle v \rangle_{y-\text{rms}}$ grows initially
665 exponentially with a growth rate that is fourth times the one of $v'_{\text{rms}}$, instead of twice in
666 figure 5.

667 In order to understand these results, figure 23 shows the amplitudes of the Fourier modes
668 of $v$ at $t = 550$. We can see that only the discrete wavenumbers aligned along oblique lines
669 $k_z/k_0 = 12n-2k_y/k_0$, where $n$ is an integer, are non-zero at that time. The first line (i.e. $n = 1$)
670 correspond to wavevectors $k = k_{f+} + pk_{2+}$, where $p$ is an integer, generated by $p-1$ quadratic
671 interactions of $k_{2+} = (-k_0, 2k_0)$ with itself and one interaction with the forced wavevector
672 $k_{f+} = (4k_0, 4k_0)$. Hence, the first interaction $p = 1$ leads to the resonant wavevector
673 $k_{1+} = (3k_0, 6k_0)$ whereas subsequent non-resonant interactions, $p = 2, 3, 4, \dots$, force non-
674 resonant wavevectors with a decreasing horizontal wavenumber: $(2k_0, 8k_0)$, $(k_0, 10k_0)$,
675 (0, 12k_0) and so on.

676 As seen in figure 24(a), the growth rate of these successive forced modes is $p\sigma$, where
677 $\sigma$ is the growth rate of the resonant modes. Hence, four successive nonlinear interactions
678 can generate a non-resonant mode with a zero horizontal wavenumber that grows at a rate
679 $4\sigma$. This explains why $\langle v \rangle_{y-\text{rms}}$ in figure 22 grows exponentially with a growth rate fourth
680 times the one of $v'_{\text{rms}}$. Similarly, the modes along the next oblique lines ($n > 1$) can be
681 produced from those of the first line by additional $n-1$ non-linear interactions with $k_{1+}$:
682 $k = k_{f+} + pk_{2+} + (n-1)k_{1+}$. A given line $n$ will then produce the shear mode: (0, 12nk$_0$).
683 Figure 24(b) shows that all the other shear modes, i.e. those different from (0, 12nk$_0$), grow
684 only later and at a different growth rate, probably as a consequence of a secondary triadic
685 instability of the modes produced by the primary instability.

686 By writing the shear mode $k_v = (0, 12k_0)$ as follows $k_v = k_{1+} + 3k_{2+}$, the resonance
687 condition (5.5) can be rewritten in terms of a sextet resonance:

$$k_v = k_{f+} + k_{2+} + k_{2+} + k_{2+} + k_{2+}, \qquad (7.2a)$$

$$\omega_v = \omega_{f+} + \omega_{2+} \pm \omega_{2+} \pm \omega_{2+} \pm \omega_{2+}, \qquad (7.2b)$$

691 where $\omega_v = 0$ and the sign of the frequencies $\pm\omega_{2+}$ in (7.2b) can be chosen so that the
692 temporal resonance condition is satisfied.

693 In summary, shear modes are generated by quadratic interactions between enslaved non-
694 resonant modes Newell et al. (1993); Manneville (2014) as before even if the precise
695 mechanism is different from the one described in the previous sections. For arbitrary forcing
696 wavenumber ($k_{yf}$, $k_{zf}$) and pair of resonant wavevectors $k_1 = (n_1, m_1)k_0$, $k_2 = (n_2, m_2)k_0$,
697 this mechanism should generate a shear mode after $p = -n_2/M$ interactions of the wavevector
698 $k_1$ and $q = n_1/M$ interactions of the wavevector $k_2$, where $M$ is the greatest common divisor
699 of $(n_1, n_2)$, since

$$pk_{1y} + qk_{2y} = 0. \qquad (7.3)$$

701 Then, the shear mode $(0, pk_{1z} + qk_{2z})$ will grow at rate $(p+q)\sigma$. An additional condition to
702 have a steady and large forced response is that $(p+q)$ should be even. This ensures that the
703 forcing of shear modes can contain one contribution with zero frequency, i.e. corresponding
704 to their natural frequency.

705 Finally, it is worth to point out that even if the shear modes can grow twice faster by
<figure><img src="image_29.png" /><figcaption>Figure 23: Amplitudes of the Fourier modes $\hat{v}$ of the horizontal velocity with the forcing (7.1) at $t = 550$ for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$. The colour scale is given in natural logarithm. The wavevectors labelled $k_{1+}, \hat{k}_{f+}$ are the resonant wavevectors. They drive the other wavevectors that are all non-resonant. The evolutions of the labelled wavevectors are displayed in figure 24.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 24: Evolution of the amplitudes of the modes (grey lines) in the wavenumber ranges $[0,24k_0]$x$[0,24k_0]$ for $F_h = 0.4$, $Re = 8840$ and $k_{yf} = k_{zf} = 4k_0$ with the forcing (7.1). In (a), the modes labelled in figure 23 are plotted in color. In (b), all the shear modes $(0, k_z)$ with $0 < k_z \le 11k_0$ are highlighted in color. The y axis corresponds to natural logarithm. A moving average over a time interval $\Delta t = 5$ has been applied in order to filter rapid oscillations.</figcaption></figure>

the sextet mechanism compared to the quartet mechanism, the former mechanism was not visible in the simulations with two pairs of forced wave-vectors ($\pm k_{f+}, \pm k_{f-}$). Comparing figure 24a and figure 13c, we can see that the initial amplitude of the common shear mode $(0, 12k_0)$ which is amplified by the two mechanisms, is much smaller when the forcing consists of a single pair of forced wave-vectors: $\exp(-22) \approx 3.10^{-10}$ than with two pairs:
711 $\exp(-17) \approx 4.10^{-8}$ when it starts to grow exponentially. The exponential growth begins also
712 later in the first case: $t_i \approx 300$ (figure 24a) than in the second case: $t_i \approx 200$ (figure 13c).

713 Hence, in figure 13c, the potential growth at rate $4\sigma$ is masked by the growth at rate $2\sigma$ since
714 the latter begins before and from a larger amplitude.

715 **8. Conclusion**

716 We have performed direct numerical simulations of the dynamics of a two-dimensional
717 strongly stratified flow in a square periodic domain forced in the momentum equation by a
718 steady mode with vorticity of the form $\sin(k_{yf}y)\sin(k_{zf}z)$. Various values of the forcing
719 wavenumbers have been investigated in the ranges: $2k_0 \le k_{yf} \le 5k_0$, $2k_0 \le k_{zf} \le 5k_0$
720 but we have first investigated in details the values $k_{yf} = k_{zf} = 4k_0$. For this case, starting
721 from small random noise velocity perturbations, the flow evolves spontaneously towards
722 turbulence and exhibits a horizontal mean flow varying along the vertical. The initial growth
723 of the horizontal mean flow is exponential with a growth rate twice the one of the velocity
724 fluctuations. Initially, the horizontal layers are rather regularly distributed with a thickness
725 corresponding typically to half a wavelength of the forcing. However, later on, their size
726 tends to increase and their distribution become irregular. Such evolution has been found
727 to occur for a wide range of the forcing intensity and for different values of the buoyancy
728 frequency and viscosity. The root-mean squares of the horizontal mean flow and of the
729 velocity fluctuations non-dimensionalized by $Nl_0$ depend mostly on the Froude number $F_h$
730 and not on the Reynolds number $Re$. The flows have similar kinetic and potential energy
731 spectra as two-dimensional strongly stratified turbulence generated by random stochastic
732 forcing (Smith 2001; Kumar et al. 2017; Fitzgerald & Farrell 2018a; Linares 2020).

733 Thanks to the use of a deterministic forcing, we have been able to provide clear evidences
734 for the mechanisms leading to the emergence of layers and the transition to turbulence. The
735 detailed analysis of the evolutions of the different Fourier modes of the horizontal velocity
736 shows that a triadic instability involving the forced wavevectors ($\pm k_{yf}$, $\pm k_{zf}$) first occurs.
737 It is slightly different from the classical triadic instability of weak internal waves because
738 the forced mode is not oscillating at its natural frequency but is steady. As is well-known,
739 a mean flow can not be generated directly by such triadic instability. However, the growth
740 of the resonant internal waves is accompanied by the amplification of six non-resonant
741 enslaved modes at the same growth rate for each triad. Then, the quadratic interactions
742 between the latter modes from two different triads involving forced wavevectors that are
743 symmetric with respect to the horizontal can generate a mean horizontal flow varying along
744 the vertical. Accordingly, these shear modes grow exponentially at a rate twice those of the
745 primary resonant internal waves as observed in the numerical simulations. Such mechanism
746 is equivalent to the resonant quartet mechanism proposed by Newell (1969) and Smith
747 & Waleffe (1999) for inertial waves in rotating turbulence and Rossby waves in $\beta$-plane
748 turbulence. However, resonant quartets were only suspected to be one possible mechanism
749 among others, such as near-resonances, to explain the emergence of a mean flow. Here,
750 the use of deterministic forcing has allowed us to provide clear evidences for its role. It
751 should be also noted that resonant quartets differ significantly from the quartet instability
752 evidenced by Brunet *et al.* (2020) in experiments on wave-driven rotating turbulence. The
753 latter instability is a secondary instability of the triadic instability when some wave-vectors
754 parallel to the primary resonant wave-vectors and to the forced wave-vectors have a sufficient
755 amplitude. They can emerge in these experiments because it is the frequency of the forcing
756 that is imposed, i.e. the direction of the forced wave-vector is fixed, but not its norm. In
757 contrast, resonant quartets are not a secondary instability of the triadic instabilities but are a
758 non-linear mechanism associated to the triadic instabilities leading to the growth of enslaved
759 non-resonant modes. We expect that the same phenomenon could occur if an internal wave
760 were forced instead of a steady mode.

761 The same scenario has been observed for most other values of the forcing wavenumbers
762 $(k_{yf}, k_{zf})$. However, there exist values for which no triadic instability occurs. Accordingly,
763 the flow then simply evolves towards the forced solution (2.4). One case has been also
764 observed where two distinct pairs of internal waves are unstable with close growth rates. In
765 addition to the previous mechanism, some shear modes can be then generated via the mutual
766 interactions between resonant and non-resonant modes of the two triads. Paradoxically, the
767 mean flow has been observed to be weaker than the fluctuations in the statistically steady
768 regime for this case.

769 In a second step, we have restricted the forcing to only the form $\cos(k_{yf}y + k_{zf}z)$ such
770 that there are no forced wavevectors that are symmetric with respect to the horizontal. In
771 this case, the previous mechanism is no longer possible. However, the numerical simulations
772 for $k_{yf} = k_{zf} = 4k_0$ still exhibit the emergence of a mean flow but growing exponentially
773 at a growth rate fourth times the one of the resonant internal waves. The analysis of the
774 evolutions of the different Fourier modes show that four successive quadratic interactions
775 between non-resonant and resonant modes lead to the generation of a mean flow at a growth
776 rate four times the growth rate of the primary triadic instability. This can be interpreted in
777 terms of a sextet resonance. Thus, resonant quartet is only a special case of a more general
778 mechanism by which shear modes are generated via forced non-resonant modes.

779 For arbitrary values of $(k_{yf}, k_{zf})$, a similar scenario should be always possible provided
780 that a triadic instability occurs. However, the number $n$ of successive quadratic interactions
781 needed for the generation of shear modes will depend on each case. Accordingly, the growth
782 rate of the shear modes will be $n\sigma$, where $\sigma$ is the growth rate of the primary triadic
783 instability.

784 Therefore, for the two types of forcing studied herein, shear modes, i.e. VSHF, are
785 spontaneously generated by the interactions between non-resonant modes driven by triadic
786 instabilities. Such uprising of enslaved modes is able to efficiently channel energy from
787 the resonant modes towards the mean flow. As pointed out by Remmel et al. (2014),
788 this mechanism is absent in weak turbulence theories since only resonant interactions are
789 taken into account. Nevertheless, the link between the present simulations with a steady
790 monochromatic forcing and those using random forcing remains to be made, as shear modes
791 do not appear to grow exponentially in the latter simulations. This could be due to the fact
792 that all wave vector directions in a given wave number band are forced randomly, so that
793 a given wave vector is only forced intermittently. However, further investigations would be
794 needed to prove this.

795 Furthermore, the interactions between non-resonant modes explain the initial emergence of
796 shear modes but other mechanisms may be responsible for their maintenance and subsequent
797 evolution such as the statistical instability of homogeneous stratified turbulence (Fitzgerald
798 & Farrell 2018a,b). Hence, we believe that these different mechanisms are complementary.

799 In the future, it would be interesting to study along the same lines three-dimensional
800 stratified turbulence and different types of turbulence such as two-dimensional $\beta$-plane
801 turbulence or rotating turbulence. Three-dimensional stratified turbulence is expected to be
802 more complicated than the present configuration since shear modes are not constrained to a
803 single horizontal direction and there exist potential vorticity modes.

804 **Acknowledgements.** We thank M. Le Bars, B. Favier and B. Gallet for insightful discussions, V. Toai for technical assistance and the anonymous referees for their comments that have improved the manuscript.

806 **Funding.** A part of this work was performed using HPC resources from GENCI-IDRIS (Grant 2022-AD012A07419R1).
809 Appendix A. Stability analysis

810 To study the stability of the base solution (2.4), it is first convenient to non-dimensionalize
811 (2.1) by taking 1/N and l₀ as time and length units. Accordingly, the streamfunction and
812 buoyancy are non-dimensionalized by Nl₀² and N²l₀, respectively. Then, (2.1) can be rewritten
813 in the form

$$
\frac{\partial \tilde{\Delta} \tilde{\psi}}{\partial \tilde{t}} + J(\tilde{\psi}, \tilde{\Delta} \tilde{\psi}) = \frac{\partial \tilde{b}}{\partial \tilde{y}} + \frac{F_h}{Re} \tilde{\Delta}^2 \tilde{\psi} + \frac{q}{N^2}, \quad (A 1a)
$$

$$
\frac{\partial \tilde{b}}{\partial \tilde{t}} + J(\tilde{\psi}, \tilde{b}) = -\frac{\partial \tilde{\psi}}{\partial \tilde{y}} + \frac{F_h}{ReSc} \tilde{\Delta} \tilde{b}, \quad (A 1b)
$$

817 where $J(.,,.) = \partial./\partial\tilde{y}\partial./\partial\tilde{z} - \partial./\partial\tilde{z}\partial./\partial\tilde{y}$ is the Jacobian and tilded quantities are non-dimensional.
818

819 Then, the base solution (2.4), taken at leading order in $q_0$ for $t \to \infty$ and non-
820 dimensionalized as above, can be decomposed as the sum of four modes:

$$
\begin{equation}
\begin{pmatrix} \tilde{\psi} \\ \tilde{b} \end{pmatrix} = \begin{pmatrix} \tilde{\psi}_+ \\ \tilde{b}_+ \end{pmatrix} e^{(i\tilde{k}_{yf} \tilde{y} + i\tilde{k}_{zf} \tilde{z})} + \text{c.c.} + \begin{pmatrix} \tilde{\psi}_- \\ \tilde{b}_- \end{pmatrix} e^{(i\tilde{k}_{yf} \tilde{y} - i\tilde{k}_{zf} \tilde{z})} + \text{c.c.}, \tag{A 2}
\end{equation}
$$

822 where $\tilde{\psi}_+ = -\tilde{\psi}_- = -q_0\gamma/(4Nl_0^2)$, $\tilde{b}_+ = -\tilde{b}_- = -iq_0\beta/(4N^2l_0)$. This shows that the non-
823 dimensional base buoyancy and streamfunction scale as $\tilde{b}_s = O(F_h^2)$ and $\tilde{\psi}_s = O(F_h^3/Re)$
824 with $s = \{+, -\}$. By assuming the distinguished scaling $Re = \tilde{R}e/F_h$, we can write $\tilde{b}_s =$
825 $F_h^2 B_s$ and $\tilde{\psi}_s = F_h^4 \Psi_s$, where capital quantities are assumed to be of order unity. The total
826 perturbations are scaled similarly: $\tilde{b} = F_h^2 B$ and $\tilde{\psi} = F_h^4 \Psi$.
827

827 Therefore, the stability problem can be solved asymptotically with $F_h^2$ as a small parameter
828 following weakly interacting waves theories (LeBlond & Mysak 1981; Bouruet-Aubertot
829 et al. 1995; Staquet & Sommeria 2002; Dauxois et al. 2018). The time and the infinitesimal
830 perturbations are expanded with the Froude number as follows: $\partial/\partial\tilde{t} = \partial/\partial t_0 + F_h^2\partial/\partial t_1 + \dots$
831 and

$$
\begin{equation}
\begin{pmatrix} \psi' \\ b' \end{pmatrix} = \begin{pmatrix} \psi'_0 \\ b'_0 \end{pmatrix} + F_h^2 \begin{pmatrix} \psi'_1 \\ b'_1 \end{pmatrix} + \dots \tag{A 3}
\end{equation}
$$

833 Introducing the infinitesimal perturbations (A 3) into (A 1) leads at leading order to

$$
\frac{\partial \tilde{\Delta} \psi_0'}{\partial t_0} = \frac{\partial b_0'}{\partial \tilde{y}}, \qquad (A 4a)
$$

$$
\frac{\partial b'_0}{\partial t_0} = - \frac{\partial \psi'_0}{\partial \tilde{y}}. \qquad (A 4b)
$$

837 The solution is taken in the form

$$
\begin{pmatrix} \psi'_0 \\ b'_0 \end{pmatrix} = \begin{pmatrix} \psi'_{01}(t_1) \\ b'_{01}(t_1) \end{pmatrix} e^{(i\tilde{k}_{y1}\tilde{y} + i\tilde{k}_{z1}\tilde{z} - i\tilde{\omega}_{1t_0})} + \begin{pmatrix} \psi'_{02}(t_1) \\ b'_{02}(t_1) \end{pmatrix} e^{(i\tilde{k}_{y2}\tilde{y} + i\tilde{k}_{z2}\tilde{z} - i\tilde{\omega}_{2t_0})}, \quad (A 5)
$$

839 where $\tilde{\omega}_1$ and $\tilde{\omega}_2$ satisfy the dispersion relation (5.6) in non-dimensional form.
840 At order $F_h^2$, we get

$$ \frac{\partial \tilde{\Delta} \psi'_1}{\partial t_0} - \frac{\partial b'_1}{\partial \tilde{y}} = - \frac{\partial \tilde{\Delta} \psi'_0}{\partial t_1} + \frac{1}{\tilde{R}e} \tilde{\Delta}^2 \psi'_0, \qquad (A 6a) $$

$$ \frac{\partial b'_1}{\partial t_0} + \frac{\partial \psi'_1}{\partial \tilde{y}} = - \frac{\partial b'_0}{\partial t_1} - J(\psi'_0, B) + \frac{1}{\tilde{R}e S_c} \tilde{\Delta} b'_0. \qquad (A 6b) $$

844 If the resonance conditions (5.5) are satisfied for the Fourier components associated to $B_+$, the compatibility conditions to be able to solve these equations are

$$ \frac{\partial \psi'_{01}}{\partial t_1} = - \frac{1}{2\tilde{R}e} \tilde{k}_1^2 \left(1 + \frac{1}{S_c}\right) \psi'_{01} - \frac{\tilde{\omega}_1 \tilde{\alpha}}{2\tilde{k}_{y1}} \psi'_{02} B_+ e^{2i\delta_{\omega}t_1}, \qquad (A 7) $$

$$ \frac{\partial \psi'_{02}}{\partial t_1} = - \frac{1}{2\tilde{R}e} \tilde{k}_2^2 \left(1 + \frac{1}{S_c}\right) \psi'_{02} + \frac{\tilde{\omega}_2 \tilde{\alpha}}{2\tilde{k}_{y2}} \psi'_{01} B_+^* e^{-2i\delta_{\omega}t_1}, \qquad (A 8) $$

849 where $\tilde{k}_j = \|\tilde{k}_j\|$ and $\tilde{\alpha} = \tilde{k}_{yf} \tilde{k}_{z1} - \tilde{k}_{zf} \tilde{k}_{y1} = \tilde{k}_{yf} \tilde{k}_{z2} - \tilde{k}_{zf} \tilde{k}_{y2}$ because of (5.5a). In (A 7-A 8), it has been assumed that the resonance condition (5.5b) could be not satisfied 850 exactly by introducing the detuning parameter $\tilde{\delta}_{\omega} = (\tilde{\omega}_1 - \tilde{\omega}_2)/2$. By assuming that this 851 parameter is small such that $\tilde{\delta}_{\omega} = \delta_{\omega} F_h^2$ with $\delta_{\omega} = O(1)$, we have rewritten the phase shift 852 in terms of the slow time $t_1$: $\tilde{\delta}_{\omega}t_0 = \delta_{\omega}t_1$. 853

854 By writing $\psi'_{01} = \psi_{01} \exp(\tilde{\sigma}t_1 + i\delta_{\omega}t_1)$ and $\psi'_{02} = \psi_{02} \exp(\tilde{\sigma}t_1 - i\delta_{\omega}t_1)$ and by rescaling 855 the growth rate $\sigma = NF_h^2\tilde{\sigma}$ and all the other variables, we obtain in dimensional form, i.e. in 856 terms of untilded quantities

$$ \sigma = -\frac{\nu}{4} \left(1 + \frac{1}{S_c}\right) \left(k_1^2 + k_2^2\right) \pm \frac{1}{4} \sqrt{\left(\nu \left(1 + \frac{1}{S_c}\right) (k_1^2 - k_2^2) + 2i(\omega_1 - \omega_2)\right)^2 - \frac{4\alpha^2 |b_+|^2 \omega_1 \omega_2}{N^4 k_{y1} k_{y2}}}. \qquad (A 9) $$

858 The destabilization being due to the last term of (A 9), we see that a necessary condition for 859 instability when the resonance is exact, $\omega_1 = \omega_2$, is that $k_{y1}$ and $k_{y2}$ be of opposite sign so 860 that this term is positive.

861 The same expression is obtained when considering $(b_-, \psi_-)$ except the changes $k_{zf} \rightarrow -k_{zf}$ and $b_+ \rightarrow b_-$. 862

REFERENCES

863 AUGIER, P., BILLANT, P. & CHOMAZ, J.-M. 2015 Stratified turbulence forced with columnar dipoles: numerical study. *J. Fluid Mech.* **769**, 403–443.

864 BENIELLI, D. & SOMMERIA, J. 1998 Excitation and breaking of internal gravity waves by parametric instability. *J. Fluid Mech.* **374**, 117–144.

865 BOURUET-AUBERTOT, P., SOMMERIA, J. & STAQUET, C. 1995 Breaking of standing internal gravity waves through two-dimensional instabilities. *J. Fluid Mech.* **285**, 265–301.

866 BOURUET-AUBERTOT, P., SOMMERIA, J. & STAQUET, C. 1996 Stratified turbulence produced by internal wave breaking: two-dimensional numerical experiments. *Dyn. Atmos. Oceans* **23** (1-4), 357–369.

871 BRETHOUWER, G., BILLANT, P., LINDBORG, E. & CHOMAZ, J.-M. 2007 Scaling analysis and simulation of strongly stratified turbulent flows. *J. Fluid Mech.* **585**, 343–368.

872 BRUNET, M., GALLET, B. & CORTET, P.-P. 2020 Shortcut to geostrophy in wave-driven rotating turbulence: the quartetic instability. *Phys. Rev. Lett.* **124** (12), 124501.

873 DE BRUYN KOPS, S. M. & RILEY, J. J. 1998 Direct numerical simulation of laboratory experiments in isotropic turbulence. *Phys. Fluids* **10** (9), 2125–2127.

874 CAULFIELD, C. P. 2021 Layering, instabilities, and mixing in turbulent stratified flows. *Annu. Rev. Fluid Mech.* **53**, 113–145.
879 DAUXOIS, T., JOUBAUD, S., ODIER, P. & VENAILLE, A. 2018 Instabilities of internal gravity wave beams. *Annu. Rev. Fluid Mech.* **50**, 131–156.

880 DELONCLE, A. 2007 Three dimensional instabilities in stratified fluids. PhD thesis, Ecole Polytechnique, Palaiseau, France.

881 DELONCLE, A., BILLANT, P. & CHOMAZ, J.-M. 2008 Nonlinear evolution of the zigzag instability in stratified fluids: a shortcut on the route to dissipation. *J. Fluid Mech.* **599**, 229–239.

885 FARRELL, B. F. & IOANNOU, P. J. 2019 Statistical state dynamics: A new perspective on turbulence in shear flow. *Zonal Jets Phenomenology, Genesis, and Physics. Ed. Boris Galpirin and Peter L. Read.* Cambridge University Press 2019, pp. 380–400.

888 FITZGERALD, J. G. & FARRELL, B. F. 2018a Statistical state dynamics of vertically sheared horizontal flows in two-dimensional stratified turbulence. *J. Fluid Mech.* **854**, 544–590.

890 FITZGERALD, J. G. & FARRELL, B. F. 2018b Vertically sheared horizontal flow-forming instability in stratified turbulence: Analytical linear stability analysis of statistical state dynamics equilibria. *J. Atmos. Sci.* **75** (12), 4201–4227.

893 GALMICHE, M., THUAL, O. & BONNETON, P. 2000 Wave/wave interaction producing horizontal mean flows in stably stratified fluids. *Dyn. Atmos. Oceans* **31** (1-4), 193–207.

895 GALPERIN, B. & READ, P. L. 2019 *Zonal jets: Phenomenology, genesis, and physics*. Cambridge University Press.

897 HERBERT, C., MARINO, R., ROSENBERG, D. & POUQUET, A. 2016 Waves and vortices in the inverse cascade regime of stratified turbulence with or without rotation. *J. Fluid Mech.* **806**, 165–204.

898 HERRING, J. R. & MÉTAIS, O. 1989 Numerical experiments in forced stably stratified turbulence. *J. Fluid Mech.* **202**, 97–115.

901 KOUDELLA, C.R. & STAQUET, C. 2006 Instability mechanisms of a two-dimensional progressive internal gravity wave. *J. Fluid Mech.* **548**, 165–196.

903 KUMAR, A., VERMA, K. M. & SUKHATME, J. 2017 Phenomenology of two-dimensional stably stratified turbulence under large-scale forcing. *J. Turbul.* **18** (3), 219–239.

905 LAVAL, J.-P., MCWILLIAMS, J. C. & DUBRULLE, B. 2003 Forced stratified turbulence: successive transitions with reynolds number. *Phys. Rev. E* **68** (3), 036308.

907 LE REUN, T., GALLET, B., FAVIER, B. & LE BARS, M. 2020 Near-resonant instability of geostrophic modes: beyond greenspan's theorem. *J. Fluid Mech.* **900**, R2.

909 LEBLOND, P. H. & MYSAK, L. A. 1981 *Waves in the Ocean*. Elsevier.

910 LELONG, M.-P. & RILEY, J. J. 1991 Internal wave—vortical mode interactions in strongly stratified flows. *J. Fluid Mech.* **232**, 1–19.

912 LEMASQUERIER, D., FAVIER, B. & LE BARS, M. 2021 Zonal jets at the laboratory scale: hysteresis and rossby waves resonance. *J. Fluid Mech.* **910**, A18.

914 LINARES, M. C. 2020 Numerical study of 2d stratified turbulence forced by internal gravity waves. PhD thesis, Université Grenoble Alpes.

916 LINDBORG, E. 2006 The energy cascade in a strongly stratified fluid. *J. Fluid Mech.* **550**, 207–242.

917 MAFFIOLI, A. 2019 Asymmetry of vertical buoyancy gradient in stratified turbulence. *J. Fluid Mech.* **870**, 266–289.

919 MAFFIOLI, A., BRETHOUWER, G. & LINDBORG, E. 2016 Mixing efficiency in stratified turbulence. *J. Fluid Mech.* **794**, R3.

921 MANNEVILLE, P. 2014 *Dissipative Structures and Weak Turbulence*. Academic Press.

922 McEWAN, A. D. 1971 Degeneration of resonantly-excited standing internal gravity waves. *J. Fluid Mech.* **50** (3), 431–448.

924 McEWAN, A. D., MANDER, D. W. & SMITH, R. K. 1972 Forced resonant second-order interaction between damped internal waves. *J. Fluid Mech.* **55** (4), 589–608.

926 NEWELL, A. C. 1969 Rossby wave packet interactions. *J. Fluid Mech.* **35** (2), 255–271.

927 NEWELL, A. C., PASSOT, T. & LEGA, J. 1993 Order parameter equations for patterns. *Annu. Rev. Fluid Mech.* **25**, 399–453.

929 PHILLIPS, O. M. 1968 The interaction trapping of internal gravity waves. *J. Fluid Mech.* **34** (2), 407–416.

930 PHILLIPS, O. M. 1972 Turbulence in a strongly stratified fluid—is it unstable? *Deep Sea Res.* **19** (1), 79–81.

931 PLUMB, R. A. & McEWAN, A. D. 1978 The instability of a forced standing wave in a viscous stratified fluid: A laboratory analogue of the quasi-biennial oscillation. *J. Atmos. Sci.* **35** (10), 1827–1839.

933 POSMENTIER, E. S. 1977 The generation of salinity finestructure by vertical diffusion. *J. Phys. Oceanogr.* **7** (2), 298–300.
935 REMMEL, M., SUKHATME, J. & SMITH, L. M. 2014 Nonlinear gravity-wave interactions in stratified turbulence. *Theor. Comput. Fluid Dyn.* **28**, 131–145.

936 RHINES, P. B. 1975 Waves and turbulence on a beta-plane. *J. Fluid Mech.* **69** (3), 417–443.

938 RILEY, J. J. & LELONG, M.-P. 2000 Fluid motions in the presence of strong stable stratification. *Annu. Rev. Fluid Mech.* **32** (1), 613–657.

940 SEMIN, B., FACCHINI, G., PÉTRÉLIS, F. & FAUVE, S. 2016 Generation of a mean flow by an internal wave. *Phys. Fluids* **28** (9).

942 SEMIN, B., GARROUM, N., PÉTRÉLIS, F. & FAUVE, S 2018 Nonlinear saturation of the large scale flow in a laboratory model of the quasibiennial oscillation. *Phys. rev. lett.* **121** (13), 134502.

944 SMITH, L. M. 2001 Numerical study of two-dimensional stratified turbulence. *Contemporary Mathematics* **283**, 91–106.

946 SMITH, L. M. & WALEFFE, F. 1999 Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence. *Phys. Fluids* **11** (6), 1608–1622.

948 SMITH, L. M. & WALEFFE, F. 2002 Generation of slow large scales in forced rotating stratified turbulence. *J. Fluid Mech.* **451**, 145–168.

950 SRINIVASAN, K. & YOUNG, W. R. 2012 Zonostrophic instability. *J. Atmos. Sci.* **69** (5), 1633–1656.

951 STAQUET, C. & SOMMERIA, J. 2002 Internal gravity waves: from instabilities to turbulence. *Annu. Rev. Fluid Mech.* **34** (1), 559–593.

953 SUTHERLAND, B. R. 2010 *Internal gravity waves*. Cambridge university press.

954 THORPE, S. A. 1968 On standing internal gravity waves of finite amplitude. *J. Fluid Mech.* **32** (3), 489–528.

955 WAITE, M. L. & BARTELLO, P. 2004 Stratified turbulence dominated by vortical motion. *J. Fluid Mech.* **517**, 281–308.

957 WEDI, N. P. & SMOLARKIEWICZ, P. K. 2006 Direct numerical simulation of the plumb–mcewan laboratory analog of the qbo. *J. Atmos. Sci.* **63** (12), 3226–3252.