Near-resonant instability of geostrophic modes: beyond Greenspan's theorem

Abstract

We explore the near-resonant interaction of inertial waves with geostrophic modes in rotating fluids via numerical and theoretical analysis. When a single inertial wave is imposed, we find that some geostrophic modes are unstable above a threshold value of the Rossby number $kRo$ based on the wavenumber and wave amplitude. We show this instability to be caused by triadic interaction involving two inertial waves and a geostrophic mode such that the sum of their eigenfrequencies is non-zero. We derive theoretical scalings for the growth rate of this near-resonant instability. The growth rate scaled by the global rotation rate is proportional to $(kRo)^2$ at low $kRo$ and transitions to a $kRo$ scaling for larger $kRo$. These scalings are in excellent agreement with direct numerical simulations. This instability could explain recent experimental observations of geostrophic instability driven by waves.

1. Introduction

Rotating turbulent flows are ubiquitous in geophysical and astrophysical systems such as stellar interiors, planetary cores, oceans and atmospheres. In a large number of numerical simulations and experiments (see the review by Godeferd & Moisy (2015)), rotating turbulence is observed to develop a strong anisotropy and to spontaneously form vortices that are invariant along the rotation axis. The latter correspond to a first-order balance between the Coriolis force and pressure gradients and are called ‘geostrophic modes’. Yet, the systematic observation of strong geostrophic modes is at odds with various evidence suggesting that rotating turbulence could also be dominated by inertial waves that are sustained by the restoring action of the Coriolis force. Recent numerical (Le Reun etal. 2017) and experimental (Le Reun, Favier & Le Bars 2019; Brunet, Gallet & Cortet 2020) studies have shown that injecting energy in waves solely creates a turbulent state comprising of inertial waves only when the forcing amplitude is sufficiently small, i.e. a discrete version of inertial wave turbulence (Galtier 2003). It is only at larger forcing amplitudes that a secondary instability leads to the classical geostrophic-dominated turbulence. Asymptotic theories describing rotating turbulence in the limit of vanishing forcing amplitude and dissipation also suggest that waves could dominate over geostrophic modes in such a regime (Bellet etal. 2006; Sagaut & Cambon 2018). Hence, although bi-dimensionalisation in the form of geostrophic eddies has been commonly observed, it may not be the only equilibrium state of rotating turbulence, be it at moderate (Yokoyama & Takaoka 2017; Favier, Guervilly & Knobloch 2019) as well as small (van Kan & Alexakis 2019) forcing amplitudes. In addition, the nature of the forcing seems fundamental in determining the equilibrium state of rotating turbulence.

These results altogether call for a better understanding of the fundamental processes by which waves give rise to balanced geostrophic modes. The studies of Le Reun etal. (2019) and Brunet etal. (2020) suggest that such a transfer occurs through an instability. Although wave-to-wave interactions are primarily governed by triadic resonance (Vanneste 2005; Bordes etal. 2012), they cannot account for wave-to-geostrophic transfers (Greenspan 1969), at least in the asymptotic limit of vanishing velocity amplitude and dissipation. Several alternative mechanisms, outside the framework of Greenspan's theorem, have been proposed. Four-mode interactions can transfer energy from waves to geostrophic flows, either directly (Newell 1969; Smith & Waleffe 1999) or through an instability mechanism (Kerswell 1999; Brunet etal. 2020). The growth rate of such an instability scales as $Ro^2$, with $Ro$ the dimensionless wave amplitude or Rossby number. It develops over longer time scales than triad-type interactions between waves. The other inviscid mechanism that has been proposed to account for wave-geostrophic transfer is quasi-resonant triadic interaction (Newell 1969; Smith & Waleffe 1999), that is, a triad between waves whose frequencies do not exactly satisfy the resonance condition (Bretherton 1964). Their presence and their role in the bi-dimensionalisation of rotating turbulence have been assessed by several numerical studies (Smith & Lee 2005; Alexakis 2015; Clark di Leoni & Mininni 2016). While it has been shown that such triads can directly transfer energy from two pre-existing waves to geostrophic modes, we show that this transfer can arise spontaneously through an instability mechanism. More precisely, we show with direct numerical simulations (DNS) and theoretical analysis that there exists a linear mechanism by which a single inertial wave drives exponential growth of geostrophic modes through near-resonant triadic interaction.

2. The stability of a single inertial wave

2.1. Governing equations and numerical methods

Let us consider an incompressible fluid rotating at rate $\Omega \boldsymbol {e}_z$. We investigate the stability of a single plane inertial wave with wavevector $\boldsymbol {k}$ and eigenfrequency $\omega_{k}^{s_k}$. Its amplitude is proportional to $\boldsymbol {U}_w$ with

(2.1)\begin{equation} \boldsymbol{U}_w = \boldsymbol{h}_{\boldsymbol{k}}^{s_k} \exp \textrm{i}(\boldsymbol{k \cdot x}- \omega_{k}^{s_k} t) + \textrm{c.c.} \end{equation}

Here, $\boldsymbol {h}_{\boldsymbol {k}}^{s_k}$ is a helical mode, that is, when $\boldsymbol {k}$ is not parallel to the axis of rotation $\boldsymbol {e}_z$ (Cambon & Jacquin 1989; Waleffe 1992)

(2.2)\begin{equation} \boldsymbol{h}_{\boldsymbol{k}}^{s_k} = \frac{1}{\sqrt{2}} \left(\frac{ \left(\boldsymbol{k} \times \boldsymbol{e}_z \right) \times \boldsymbol{k}}{\vert(\boldsymbol{k} \times \boldsymbol{e}_z) \times \boldsymbol{k} \vert } + \textrm{i}\, s_k \frac{\boldsymbol{k} \times \boldsymbol{e}_z}{\vert \boldsymbol{k} \times \boldsymbol{e}_z \vert} \right), \end{equation}

where $s_k = \pm 1$ is the sign of the helicity of the plane wave. If $\boldsymbol {k}$ is parallel to $\boldsymbol {e}_z$, $\boldsymbol{h}_{\boldsymbol{k}}^{s_k} = \boldsymbol {e}_x + \textrm {i}\, s_k \boldsymbol {e}_y$. Here, $\boldsymbol {U}_w$ automatically satisfies the incompressibility condition since $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {U}_w \propto \boldsymbol {k} \boldsymbol {\cdot } \boldsymbol {h}_{\boldsymbol {k}}^{s_k} = 0$ and $\boldsymbol {U}_w$ satisfies the linearised rotating Euler equation,

(2.3)\begin{equation} \partial_t \boldsymbol{U}_w + 2 \boldsymbol{\Omega} \times \boldsymbol{U}_w ={-} \boldsymbol{\nabla} {\rm \pi}_w, \end{equation}

provided that $\omega_{k}^{s_k}$ and $\boldsymbol {k}$ are related by the dispersion relation of inertial waves

(2.4)\begin{equation} \omega_{k}^{s_k} = 2 s_k \Omega \frac{k_z}{k} = 2 s_k \Omega \cos \theta, \end{equation}

$\theta$ being the angle between the wavevector $\boldsymbol {k}$ and the rotation axis $\boldsymbol {e}_z$ (ranging from $0$ to ${\rm \pi}$), and $k = \vert \boldsymbol {k} \vert$. We solve for the time evolution of perturbations $\boldsymbol {u}$ to the wave $\boldsymbol {U}_w$ maintained at a constant amplitude via the following set of equations:

(2.5a)\begin{gather} \partial_t \boldsymbol{u} + Ro \left(\boldsymbol{U}_w \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} +\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U}_w\right) + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + 2 \boldsymbol{e}_z \times \boldsymbol{u} ={-} \boldsymbol{\nabla} {\rm \pi} + E \nabla^2 \boldsymbol{u}, \end{gather}

(2.5b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}

where time is scaled by $\Omega ^{-1}$ and length by the domain size $L$. We have introduced the Ekman number $E = \nu /(L^2 \Omega )$, $\nu$ being the kinematic viscosity, and an input Rossby number $Ro$ quantifying the dimensionless amplitude of the plane wave.

Equations (2.5) are solved numerically in a triply periodic cubic box using the code ‘Snoopy’ (Lesur & Longaretti 2005). The dynamics of the perturbation flow $\{ \boldsymbol {u}, {\rm \pi} \}$ is determined with pseudospectral methods, that is, $\{ \boldsymbol {u}, {\rm \pi} \}$ is decomposed into a truncated sum of Fourier modes $\{ \hat {\boldsymbol {u}}_{\boldsymbol {q}},\hat {{\rm \pi} }_{\boldsymbol {q}} \} \textrm {e}^{ \textrm {i}\, \boldsymbol {q} \boldsymbol {\cdot } \boldsymbol {x}}$. A wavevector $\boldsymbol {q}$ writes $2 {\rm \pi} /L ( n_x, n_y, n_z )$ where $n_{x,y}$ are integers varying from $-N$ to $N$ and $n_z$ is an integer ranging from $0$ to $N$ because of Hermitian symmetry. In the following, $N = 96$; higher resolutions have been tested and yield the exact same results. The temporal dynamics of the modes $\hat {\boldsymbol {u}}_{\boldsymbol {q}}$ is solved using a third-order Runge–Kutta method. Note that the size of the box $L$ is artificial and we thus expect our results to depend on the intrinsic Rossby number based on the imposed wavelength $k Ro$ rather than $Ro$ alone.

2.2. Numerical results

Keeping the Ekman number to $E = 10^{-6}$, two simulations of the stability of the inertial wave $\boldsymbol {k} = 2 {\rm \pi} [ 4,0,8 ]$ (with $s_k = 1$) are performed at low ($Ro = 2.83 \times 10^{-3}$) and moderate ($Ro = 2.83 \times 10^{-2}$) wave amplitudes. They are both initiated with a random noise comprising wavenumbers ranging between $0$ and $15 {\rm \pi}$, with very small initial amplitude. The use of spectral methods allows the separation of kinetic energy of the perturbation flow $\boldsymbol {u}$ into a two-dimensional component $\mathcal{E}_{\text{2-D}}$, accounting for all modes $\boldsymbol {q}$ with $q_z=0$, and a complementary three-dimensional component $\mathcal {E}_{\text{3-D}}$. In addition, performing a Fourier transform in space and time allows the projection of the kinetic energy in the sub-space of the dispersion relation of inertial waves to draw the spatiotemporal spectrum $\mathcal {E}(\theta , \omega )$ (Yarom & Sharon 2014; Le Reun etal. 2017). Note that, in these maps, $\theta$ is restricted to $[0, {\rm \pi} /2]$ since the flow is real and only the wavevectors $\boldsymbol {q}$ with $q_z \geq 0$ are simulated due to the Hermitian symmetry. Moreover, the spectra are symmetrised with respect to $\omega =0$ and the maps are shown as a function of $\vert \omega \vert$ to be more compact.

The kinetic energy time series and maps are shown for both simulations in figure 1. At low imposed wave amplitude (figure 1a,b), three-dimensional perturbations dominate the growth of the instability. The energy map displays two spots aligned on the dispersion relation with negative frequencies $\omega _{1,2}$ such that $\omega_{k}^{s_k} + \omega _1 +\omega _2 = 0$ which is indicative of several waves undergoing triadic resonance with the imposed mode. The growth of two-dimensional modes is delayed and their growth rate is approximately two times larger than the rate of three-dimensional modes. Removing the nonlinear term $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}$ (see (2.5)) suppresses the growth of two-dimensional modes which are thus not unstable themselves, at least on the time scale of the growth of waves. In fact, we find that $\mathcal{E}_{\text{2-D}} \simeq 3 \times 10^{-2} (k \mathcal {E}_{\text{3-D}})^2$ (see figure 1a) in the growth phase, which suggests that two-dimensional modes’ growth is due to direct forcing by nonlinear interaction of two growing waves involved in triadic resonances, with close frequencies and opposed vertical wavenumbers. This mechanism corresponds to the direct excitation of geostrophic modes by two waves identified by Newell (1969) and Smith & Waleffe (1999).

Figure 1. Kinetic energy time series (a and c) and heat map of $\log (\mathcal {E}(\theta ,\omega ))$ (b and d) resulting from two numerical simulations of the stability of the inertial wave $\boldsymbol {k} = 2 {\rm \pi} [4,0,8 ]$ at $Ro = 2.83 \times 10^{-3}$ and $Ro= 2.83 \times 10^{-2}$. In panels (a) and (c), the labels indicate the slope of the best fit for the exponential growth. In panel (b) and (d), the plain line materialises the dispersion relation of inertial waves and the horizontal dashed line the frequency of the imposed wave ($\omega_{k}^{s_k} \simeq 1.78$). For the spectral energy maps, the temporal Fourier transforms have been performed until $t= 800$ for panel $(b)$ and $t= 100$ for panel (d). In panel (b), we have indicated the extremal frequencies of the two energy locations.

At larger wave amplitude (figure 1c,d), this picture changes: two- and three-dimensional modes grow at the same rate from the start of the simulation. The vertical vorticity field of the growing perturbation and its vertical average are displayed in figure 2(a,b). The frequency of two-dimensional modes ($\theta ={\rm \pi} /2$) is close to $0$ according to figure 1(d), they are thus geostrophic, that is, slow as well as invariant along the $z$ axis. As it may be seen in figure2(b), the geostrophic unstable flow contains the wavevector $\boldsymbol {p} = 2 {\rm \pi} [0,5,0 ]$ (and $-\boldsymbol {p}$ for hermiticity), which grows along with other three-dimensional structures whose energy location in the $(\theta , \omega )$ space is reminiscent of triadic resonant instability, but with significantly more spreading. Note that the properties of the transition and the instability we report are robust to changes in the aspect ratio of the box, which discards any spurious effect of the discretisation (Smith & Lee 2005).

Figure 2. (a) Three-dimensional vertical vorticity field of the growing perturbation at $Ro = 2.83 \times 10^{-2}$ and (b) its vertical average. The fields are normalised by their maximum values. (c) Growth rate of the geostrophic modes $\boldsymbol{p}= 2 {\rm \pi} [0, p_{y}, 0 ]$ with $p_{y} \in \{ 1,3,5 \}$ as a function of the imposed wave amplitude $Ro$. The vertical line materialises $k Ro = 1$. The lines joining the markers are used to facilitate the identification of each curve. (d) Heat map of $\log (\mathcal {E}(\theta ,\omega ))$ in the case $p_{y} = 3$ and $Ro = 1.3 \times 10^{-2}$. The same lines as in figure 1(b,d) are reported and a vertical line materialises the angle of the modes closing the triad, that is $-\boldsymbol {k} \pm \boldsymbol {p}$.

To further characterise the growth of geostrophic modes, we carry out simulations with the same imposed wave, but we use as initial condition only one wavevector, $\boldsymbol {p} = 2 {\rm \pi} [0,p_{y},0 ]$ with $p_{y} \in \{ 1,3,5 \}$ and $s_{p} = \pm 1$, instead of random noise. The Ekman number is set to $10^{-8}$ to discard effects of viscosity in the mode growth. The growth rate of geostrophic modes is reported in figure 2(c) where $Ro$ is systematically varied. At large wave amplitude, the growth rate increases linearly, but it goes to zero at a finite value of $Ro \sim 10^{-2}$. This value is very large compared with the viscous damping $k^2 E \sim 10^{-5}$ which plays no role here, hence suggesting an inviscid mechanism.

On the kinetic energy map $\mathcal {E}(\theta ,\omega )$ (figure 2d) two particular spots appear, one at the location of geostrophic modes and the other at the angle of the vectors closing the triads between $\boldsymbol {k}$ and $\pm \boldsymbol {p}$, which is indicative of triad-type interactions. This may seem a priori in contradiction with the result of Greenspan (1969). Yet, our results are outside Greenspan's theorem framework in two aspects. First, our study is necessarily limited to finite Rossby number. Second, the frequency of the closing mode $-\boldsymbol {k} \pm \boldsymbol {p}$ (figure 2d) is outside its eigenfrequency. In fact, the discrepancy between the observed frequencies and eigenfrequencies of the closing mode corresponds to the sum of the eigenfrequencies of the three modes, ${\rm \Delta} \omega _{kpq}$. Moreover, as $Ro$ is decreased, the fastest growing mode $\boldsymbol {p}_g$ in figure 2(c) goes from $p_y = 5$ with ${\rm \Delta} \omega _{kpq}=0.22$ to $p_y = 1$ with ${\rm \Delta} \omega _{kpq}=0.01$: the triad is nearly resonant and draws closer to exact resonance. Two mechanisms may be considered to explain the growth of geostrophic modes. One hypothesis is that the growth of the geostrophic modes is a by-product of a classical resonance between, say, $\boldsymbol {k}$ and $-\boldsymbol {q}$ that forces $\boldsymbol {q}$ (by Hermitian symmetry) and hence $\boldsymbol {p}$. However, in this case, the mode closing the wave triad would have an eigenfrequency $-\omega_{k}^{s_k} + \omega_{q}^{s_q} = {\rm \Delta} \omega_{kqp} -2 \omega_{k}^{s_k}$ which is outside the range of inertial waves for small ${\rm \Delta} \omega _{kpq}$ since $\omega_{k}^{s_k} \simeq 1.78$. Instead, our results point towards a near-resonant triadic instability that transfers energy from an inertial wave to a $z$-invariant geostrophic mode.

3. Near resonance of geostrophic modes: theoretical approach

3.1. The low Rossby number limit

To provide theoretical insight into the inertial-wave destabilisation observed in the numerical simulations, we turn to linear stability analysis. The flow is $\boldsymbol {U} = \boldsymbol {u} + Ro\, \boldsymbol {U}_w$ where $Ro \,\boldsymbol {U}_w$ is the base flow ($Ro$ being the finite Rossby number) and $\boldsymbol {u} \ll Ro\, \boldsymbol {U}_w$ is now an infinitesimal perturbation. We consider a Cartesian domain with periodic boundary conditions and we decompose $\boldsymbol {U}$ into a superposition of plane waves $\boldsymbol {h}_{p}^{s_p} \exp \textrm {i} (\,\boldsymbol{p} \boldsymbol {\cdot } \boldsymbol {x}-\omega _{p}^{s_p} t)$ with time-dependent amplitudes $b_{p}^{s_p} (t)$. The Euler equation governing $\boldsymbol {U}$ is then equivalent to a set of ordinary differential equations governing the amplitudes $b_{p}^{s_p}$ (Smith & Waleffe 1999)

(3.1)\begin{equation} \frac{\mathrm{d} b_{p}^{s_p}}{\mathrm{d} t}= \sum_{\boldsymbol{k} + \boldsymbol{p} + \boldsymbol{q} = \boldsymbol{0}} \sum_{s_k, s_q ={\pm} 1 } C_{pkq}^{s_p s_k s_q}\, b_{k}^{s_k*}\, b_{q}^{s_q*} \exp ( \textrm{i}\, {\rm \Delta} \omega_{k p q} t ), \end{equation}

with

(3.2a,b)\begin{align} C_{pkq}^{s_p s_k s_q} \equiv \frac{1}{2} (s_{q} q - s_k k) \, \boldsymbol{h}_{p}^{s_p*} \boldsymbol{\cdot} \left(\boldsymbol{h}_{k}^{s_k *} \times \boldsymbol{h}_{q}^{s_q *} \right) \quad \mbox{and} \quad {\rm \Delta} \omega_{k p q} \equiv \omega_{k}^{s_k} + \omega_q^{s_q} + \omega_p^{s_p}, \end{align}

and where $k,p = \vert \boldsymbol {k} \vert , \vert \boldsymbol {p} \vert$. Here, ${\rm \Delta} \omega _{k p q}$ is the sum of the eigenfrequencies of the three modes involved in the triad. In general, maximum energy transfer between the three modes is ensured when the oscillations due to the detuning in the right-hand side of (3.1) are cancelled, that is, when ${\rm \Delta} \omega _{k p q} \rightarrow 0$. If all three modes $\boldsymbol {k}$, $\boldsymbol {p}$ and $\boldsymbol {q}$ are inertial waves, this leads to the well known mechanism of triadic resonance (Bretherton 1964; Vanneste 2005; Bordes etal. 2012). When the mode $\boldsymbol {k}$ is imposed with an amplitude $Ro$ and helicity $s_k$, $\boldsymbol {p}$ and $\boldsymbol {q}$ grow exponentially with a rate proportional to $k Ro$. This picture is changed when one of the modes, say $\boldsymbol {p}$, is geostrophic (i.e. $\omega _p^{s_p} = 0$ and $p_z = 0$, regardless of $s_p$). The spatial interaction condition $\boldsymbol {k} + \boldsymbol {p} + \boldsymbol {q} = \boldsymbol {0}$ forces $k_z = - q_z$ and the resonance condition becomes

(3.3)\begin{equation} {\rm \Delta} \omega_{k p q} = k_z \left( \frac{s_k}{k} - \frac{s_q}{q} \right) = \omega_{k}^{s_k} \frac{1}{q} \left( q - \frac{s_q}{s_k} k \right) \rightarrow 0. \end{equation}

It may be fulfilled only when $s_q = s_k$ and $q - k \rightarrow 0$. At exact resonance, these conditions impose the location of $\boldsymbol {p}$ on a circle centred on $-\boldsymbol {k}_{\perp } = -[k_x,k_y ]$ with radius $\vert \boldsymbol {k}_{\perp }\vert$. Moreover, in the governing equation for $\dot {b}_p^{s_p}$, the slow oscillation terms involve a coupling coefficient $C_{p k q}^{s_p s_k s_k} \propto k - q \rightarrow 0$. At exact resonance, the coupling coefficient vanishes: there is no energy transfer from waves to geostrophic modes, as proved by Greenspan (1969).

Nevertheless, wave-to-geostrophic transfer is still possible when the detuning ${\rm \Delta} \omega _{kpq}$ is small but non-zero (Newell 1969; Smith & Waleffe 1999; Alexakis 2015), and we investigate instabilities of geostrophic modes driven by this mechanism. Let us assume that the wave $\boldsymbol {k}$ with helicity $s_k$ is imposed with a small constant amplitude $Ro$, and that $\boldsymbol {p}$ is geostrophic. To infer from (3.1) the time evolution of the geostrophic mode amplitude, we proceed to an asymptotic expansion using a two-time method involving a fast time $\tau = t$ and a slow time $T$. The hierarchy between them must be a power of $kRo$, the intrinsic Rossby number based on the imposed wavelength. Because the wave-to-geostrophic transfer coefficient $C_{ p k q }^{s_p s_k s_k}$ vanishes as ${\rm \Delta} \omega _{kpq} \rightarrow 0$, we find via a heuristic analysis, detailed in appendix A, that $T = (kRo)^2 t$, instead of $(k Ro) \,t$ for classical wave triads. In addition, it imposes the condition that the amplitude of $\boldsymbol {p}$ be smaller by a factor $kRo$ compared with the mode closing the triad $\boldsymbol {q} = -\boldsymbol {k} -\boldsymbol {p}$. The amplitudes of the modes interacting with the imposed wave ($\,\boldsymbol {p}$ and $\boldsymbol {q}$ with both helicity signs $s_{p,q}$) are thus expanded as

(3.4)\begin{equation} \left.\begin{array}{c@{}} b_q^{s_q} = (kRo) B_{q 1}^{ s_q} (T,\tau) + (kRo)^3 B_{q 2}^{s_q} (T,\tau),\\ b_p^{s_p} = (kRo)^2 B_{p 1}^{s_p} (T,\tau) + (kRo)^4 B_{p 2}^{s_p} (T,\tau), \end{array}\right\} \end{equation}

where the $B_i^j$ are all $O(1)$. The hierarchy between orders is imposed by the need to match the slow time derivative of the leading order with the next order in the multiple scale expansion. Since the slow derivation introduces a factor $(kRo)^2$, there must be a $(kRo)^2$ hierarchy between the first and second orders.

To find the equations governing the leading-order coefficients, we follow the method of Bretherton (1964) and inject the ansatz (3.4) into (3.1). First, the leading-order coefficients are found to be independent of $\tau$ and their long-time evolution with $T$ is determined at next order. As noted by Bretherton (1964), for fast triads with detuning larger than $O((kRo)^2)$, the imposed wave only drives fast and bounded oscillations of the second-order terms, and the leading-order terms must be zero to avoid secular growth. It is only for slow triads, i.e. when ${\rm \Delta} \omega _{kpq} = O((kRo)^2)$, that the exponential term in (3.1) drives slow oscillations that contribute to secular growth. Such a condition on the detuning is consistent with the numerics: the fastest growing modes at $Ro = 1 \times 10^{-2}$ ($kRo \sim 0.6$) in figure 2(c) is $\boldsymbol {p} = 2{\rm \pi} [0,3,0]$ for which ${\rm \Delta} \omega _{kpq} \simeq 0.1 \sim 0.3 (k Ro)^2$. The condition ${\rm \Delta} \omega _{kpq} = O((kRo)^2)$ may be fulfilled only for the modes $\boldsymbol {p}$ with $s_p = \pm 1$ and $\boldsymbol {q}$ with $s_q = s_k$ and, as noted earlier, the wave-to-geostrophic transfer coefficient $C_{ p k q }^{s_p s_k s_k}$ is also $O((kRo)^2)$. Cancelling the secular growth terms at second order gives the following amplitude equations :

(3.5a,b)\begin{equation} \partial_T B_{q 1}^{s_q} = \sum_{s_p} C_{qkp}^{s_k s_k s_p}\, B_{p1}^{s_p \,*}\, \textrm{e}^{\textrm{i}\, \frac{{\rm \Delta} \omega_{kpq}}{(kRo)^2} T }\quad \mbox{and}\quad \partial_T B_{p 1}^{s_p} = \frac{C_{pkq}^{s_p s_k s_k}}{(kRo)^2} \, B_{q1}^{s_q\,*}\, \textrm{e}^{\textrm{i}\, \frac{{\rm \Delta} \omega_{kpq}}{(kRo)^2} T }, \end{equation}

the rescaled quantities $C_{ p k q }^{s_p s_k s_k}/(kRo)^2$ and ${\rm \Delta} \omega _{kpq}/(kRo)^2$ being $O(1)$. These equations have exponentially growing solutions and the complex growth rate of the instability normalised by the rotation rate is

(3.6)\begin{equation} \sigma = \textrm{i}\, \frac{{\rm \Delta} \omega_{kpq}}{2} +\frac{1}{2} \sqrt{ 4 \textstyle\sum_{s_p} C_{pkq}^{s_p s_k s_k} C_{qkp}^{s_k s_k s_p\, *} Ro^2 - {\rm \Delta} \omega_{kpq}^2 } \equiv \textrm{i}\,\frac{{\rm \Delta} \omega_{kpq}}{2} + \sigma_k(\,\boldsymbol{p};Ro). \end{equation}

Note that the product of coupling coefficients $C_{pkq}^{s_p s_k s_k} C_{qkp}^{s_k s_k s_p\, *}$ is real. The expression of the real part of growth rate, $\sigma _k(\,\boldsymbol{p};Ro)$, is consistent with the numerical findings: for a given near-resonant triad, it drops to zero at a finite value of $Ro$ and it is proportional to $Ro$ at larger $Ro$. However, in the small Rossby number limit, because the detuning ${\rm \Delta} \omega _{kpq}$ and the coupling coefficient $C_{ p k q }^{s_p s_k s_k}$ are both $O((kRo)^2)$, the maximum geostrophic growth rate remains $O((kRo)^2)$ at most.

The $(kRo)^2$ scaling governing the maximum growth rate of geostrophic near resonance is found quantitatively by expanding the frequency detuning, the transfer coefficients and then $\sigma _k(\,\boldsymbol{p};Ro)$ in the neighbourhood of exact resonance. The fact that the growth rate is non-zero only close to the exact resonance is illustrated qualitatively in figure 3(a). We thus introduce $\boldsymbol {p} = \boldsymbol {p}_0 + (kRo)^2 \boldsymbol {\delta } \boldsymbol {p}$ where $\boldsymbol {\delta }\boldsymbol {p}$ is an $O(1)$ vector in the geostrophic plane. Consider the closing modes $\boldsymbol {q}_0 = -(\boldsymbol {k} +\boldsymbol {p}_0)$ and $\boldsymbol {q} = -(\boldsymbol {k} +\boldsymbol {p})$, then $\boldsymbol {q} = \boldsymbol {q}_0 - (kRo)^2 \boldsymbol {\delta } \boldsymbol {p}$. The imposed mode $\boldsymbol {k}$ is left unperturbed. Since ${\rm \Delta} \omega _{kpq} =C_{p_0kq_0}^{s_p s_k s_k} = 0$ at exact resonance, the leading orders of the frequency detuning and the transfer coefficients are found from (3.2a,b) and (3.3) to be $O((kRo)^2)$. The perturbation of the wavevectors is thus consistent with the asymptotic expansion. At leading order,

(3.7a,b)\begin{equation} \frac{{\rm \Delta} \omega_{kpq} }{(kRo)^2} \simeq s_k \omega_{k}^{s_k} \frac{\boldsymbol{q}_0 \boldsymbol{\cdot}\boldsymbol{\delta} \boldsymbol{p} }{k^2} \quad\mbox{and} \quad \frac{C_{pkq}^{s_p s_k s_k}}{(kRo)^2} \simeq \frac{1}{2} s_k \left(\boldsymbol{h}_{p_0}^{s_p\,*} \boldsymbol{\cdot} \left( \boldsymbol{h}_k^{s_k\,*} \times \boldsymbol{h}_{q_0}^{s_k\,*} \right)\right) \frac{\boldsymbol{q}_0 \boldsymbol{\cdot}\boldsymbol{\delta} \boldsymbol{p}}{k}, \end{equation}

where we have used that $q_0 = k$ and $s_q = s_k$. In the growth rate $\sigma _k(\,\boldsymbol{p};Ro)$, the product of the coupling coefficients is

(3.8)\begin{equation} \sum_{s_p} C_{pkq}^{s_p s_k s_k} C_{qkp}^{s_k s_k s_p\, *} \simeq \frac{1}{4} (k \,Ro)^2 k^2 \frac{\boldsymbol{q}_0 \boldsymbol{\cdot} \boldsymbol{\delta} \hat{ \boldsymbol{p}}}{k^2} \sum_{s_p} \left(1 - s_p \frac{p_0}{k} \right) \left| \boldsymbol{h}_{p_0}^{s_p}\boldsymbol{\cdot} \left( \boldsymbol{h}_k^{s_k} \times \boldsymbol{h}_{q_0}^{s_k} \right) \right|^2. \end{equation}

Therefore, at leading order in powers of $Ro$, the growth rate is

(3.9)\begin{equation} 4 \sigma_k(\,\boldsymbol{p};Ro)^2 = ( \mathcal{C}_k(\,\boldsymbol{p}_0) X - {\omega_{k}^{s_k}}^2 X^2) (kRo)^4, \end{equation}

where $X \equiv (\boldsymbol {q}_0 \boldsymbol {\cdot }\boldsymbol {\delta } \boldsymbol {p})/k^2$ and $\mathcal {C}_k(\,\boldsymbol{p}_0)$ is the sum in the right-hand side of (3.8). When $\mathcal {C}_k(\,\boldsymbol{p}_0) >0$, the growth rate reaches an optimum at $X = \mathcal {C}_k /(2 {\omega_{k}^{s_k}}^2)$ with value $(kRo)^2 \mathcal {C}_k (\,\boldsymbol{p}_0) /\vert 4 \omega_{k}^{s_k}\vert$, which remains to be maximised over all exactly resonant wavevectors $\boldsymbol {p}_0$. The coefficient $\mathcal {C}_k(\,\boldsymbol{p}_0)$ is shown in figure 3(b) for several wavevectors $\boldsymbol {k}$ with different frequencies $\omega_{k}^{s_k}$. When plotted against $p_0/k_{\perp }$, all the curves $\mathcal {C}_k$ collapse on a master curve that reaches a maximum value of $1$ at $p_0 \rightarrow 0$, regardless of the helicity sign $s_k$ of the imposed wave. As $(kRo) \rightarrow 0$, the unstable geostrophic mode becomes large-scale and stems from the interaction between $\boldsymbol {k}$ and $\boldsymbol {q} \simeq -\boldsymbol {k}$. In the low Rossby number limit, the maximal growth rate is then

(3.10)\begin{equation} \sigma_{k}^{max} (Ro) = \frac{1}{4} \frac{(k Ro)^2}{\vert \omega_{k}^{s_k} \vert}. \end{equation}

We confirm that $k Ro$, the Rossby number based on the wavelength, is the relevant parameter to describe the growth rate of the instability. While each geostrophic mode follows the law (3.6), all the growth rate curves lie below a $(k Ro)^2$ upper envelope. More details are given below in § 3.3 where we compare the law (3.10) to the exact computation of the growth rate curves from (3.6).

Figure 3. (a) Map of the growth rate of geostrophic modes at $Ro=7.5\times 10^{-3}$ computed from (3.6). The imposed wave is $2 {\rm \pi} [ 4,0,8]$. The white dashed circle locates the exactly resonant geostrophic modes. The colour scale gives the amplitude of the growth rate. Where it is maximum $(p_x=0 , p_y \simeq \pm 2)$, the detuning is approximately $0.04 \sim 0.24 (kRo)^2$. (b) Plot of $\mathcal {C}_k (\,\boldsymbol{p}_0)$ on the exact resonant circle against $p_0$ normalised by the horizontal wavenumber $k_{\perp }$ for several wavevectors $\boldsymbol {k}$ with different frequencies $\omega_{k}^{s_k}$. The curve is the same regardless of the imposed helicity sign $s$. (c) Growth rate curves of the geostrophic modes as a function of the Rossby number. The geostrophic modes are sampled over 15 circles whose centres are the same as the exact resonance circle with five points on each circle. The line colour codes the frequency detuning $\vert {\rm \Delta} \omega_{kpq} \vert$. The upper envelope is compared with the law (3.10) and the upper bound (3.12).

3.2. The moderate to large Rossby number regime

The $Ro^2$ law governing the growth rate in the small Rossby number limit cannot hold as the wave amplitude is increased since the growth rate has an upper bound following an $Ro$ law. This is proved directly by multiplying (2.5) by $\boldsymbol {u}$ and integrating over the fluid domain $V$, thus giving

(3.11)\begin{equation} \sigma =\frac{1}{2} \frac{\mathrm{d} \ln \left\lVert\boldsymbol{u}\right\rVert_2^2 }{\mathrm{d} t} ={-} \frac{Ro}{\left\lVert u \right\rVert^2_2} \int_V \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U}_w \boldsymbol{\cdot} \boldsymbol{u} \end{equation}

where $\left \lVert \boldsymbol {\cdot }\right \rVert _n$ denotes the $L_n$-norm. In virtue of Hölder's inequality (Gallet 2015),

(3.12)\begin{equation} \vert \sigma \vert \leq Ro \left\lVert\boldsymbol{\nabla} \boldsymbol{U}_w\right\rVert_{\infty} \leq k Ro. \end{equation}

This upper bound applies in particular to the low $Ro$ scaling (3.10), which thus holds up to $k Ro \sim 1$ at most. Note that, beyond that point, non-triad type instabilities (shear, centrifugal) may add to near resonance in driving the dynamics of the flow excited by the maintained wave.

3.3. Comparison with exact computation

In figure 3, we sample the growth rate curves given by (3.6) of many geostrophic modes in near-resonant interaction with $\boldsymbol {k} = 2 {\rm \pi} [4,0,8]$, as functions of the Rossby number. We notice that each growth rate curve follows an $Ro$ scaling at sufficiently large $Ro$, which is consistent with the asymptotic expansion carried out in § 3.1. However, because for all modes the growth rate vanishes at a finite value of $Ro$ (which decreases to 0 close to exact resonance), the upper envelope is an $Ro^2$ law that perfectly matches the theoretical prediction (3.10). This instability is thus fundamentally different from four-modes interactions for which each growth rate curve follows a $Ro^2$ law (Kerswell 1999; Brunet etal. 2020). Besides, in agreement with the upper bound (3.12), we observe the $Ro^2$ maximum growth rate law to break down above $k Ro \sim 1$. Beyond that, the near-resonant growth rate derived from (3.6) remains below the upper bound and follows an $Ro$ law. However, additional instabilities (shear, centrifugal, etc.) may also drive the growth of geostrophic modes in this regime.

3.4. Finite size and viscous effects

Finite size domain translates into discretisation of the modes. Since the exactly resonant geostrophic modes lie on a finite radius circle, they cannot be approached with arbitrarily low detuning ${\rm \Delta} \omega _{kpq}$ as $Ro \rightarrow 0$. Hence, discretisation implies the existence of a finite value of $Ro$ below which the near-resonant instability vanishes (Bourouiba 2008). Regarding viscous effects, at finite but low Ekman number, i.e. when $(k Ro)^2 \gg k^2 E$, the low Rossby number scaling is unaltered since the near-resonant modes $\boldsymbol {p}$ and $\boldsymbol {q}$ have at most similar wavenumbers to $\boldsymbol {k}$. The $O(Ro)$ upper bound on the growth rate is also unaltered by the inclusion of viscosity.

4. A refined model: the double near-resonant triad

Although promising, the model of the previous section needs to be refined to fully account for numerical results. It misses by a factor two the growth rate curves of figure2(c) and (3.6) predicts a frequency ${\rm \Delta} \omega _{kpq} /2 \neq 0$ for the geostrophic modes, while it is $0$ according to figure 2(d). We must include in our model that not only $\boldsymbol {k}$, but also $- \boldsymbol {k}$ (with the same helicity sign $s$) is imposed due to Hermitian symmetry. Consider two triads $\boldsymbol {k} + \boldsymbol {p} + \boldsymbol {q} = \boldsymbol {0}$ and $-\boldsymbol {k} + \boldsymbol {p} + \boldsymbol {q}' = \boldsymbol {0}$ (see figure 4) with helicity signs $s_q = s_{q'} = s_k$ and $s_p = \pm 1$. As shown schematically in figure 4(a), the exact resonance circles associated with $\pm \boldsymbol {k}$ coincide at $\boldsymbol {p} = 0$. In a neighbourhood of this point (the dark grey intersection in figure 4a), a geostrophic mode $\boldsymbol {p}$ is in near resonance with both $\pm \boldsymbol {k}$ with both detunings ${\rm \Delta} \omega _{kpq}$ and ${\rm \Delta} \omega _{-kpq'}$ small. It is all the more important that $\boldsymbol {p} = \boldsymbol {0}$ corresponds to the largest near-resonant growth rate as $Ro \rightarrow 0$.

Figure 4. (a) Schematic cartoon of a geostrophic mode $\boldsymbol {p}$ in near resonance with both imposed modes $\pm \boldsymbol {k}$ at the same time, based on the map of figure 3(a). The $\perp$ indices denote the horizontal component of wavevectors. (b) Comparison between theory (4.4) and the direct numerical simulations of figure 2(c). The imposed wave is $\boldsymbol {k} = 2 {\rm \pi} [4,0,8]$ and the geostrophic modes are $\boldsymbol {p} = 2 {\rm \pi} [0, p_{y}, 0]$, $p_{y}$ being given in the legend. (c) Samples of the growth rate curves $\sigma _k(\,\boldsymbol{p};Ro)$ of modes $\boldsymbol {p}$ interacting with $\pm \boldsymbol {k} = 2 {\rm \pi} [4,0,8 ]$ including simple (dark grey) and double (light grey) triad mechanisms. The double triad growth rate curves are determined by finding the roots of $\mathcal {P}$ (see (4.2)) for wavevectors $\boldsymbol {p}$ restricted to $\vert p_x\vert < {\rm \pi}$ and $\vert p_y \vert < 12{\rm \pi}$. We recall the $Ro^2$ law (3.10) and the growth rate upper bound (3.12). The dots represent the geostrophic growth rate found in the DNS with imposed wavevector $\boldsymbol {k} = 2 {\rm \pi} [4,0,8 ]$ (DNS 1), $\boldsymbol {k} = 1.5\times 2 {\rm \pi} [4,0,8 ]$ (DNS 1.5) and $\boldsymbol {k} = 2\times 2 {\rm \pi} [4,0,8 ]$ (DNS 2).

In a very similar way to the two-time asymptotic expansion of § 3.1, we can retrieve the amplitude equations governing the slowly varying envelopes $B_{p1}^{s_p}$ and $B_{Q1}^s$, with $\boldsymbol {K}= \pm \boldsymbol {k}$, $\boldsymbol {Q} = -\boldsymbol {K} - \boldsymbol {p}$ (i.e. $\boldsymbol{Q} = \boldsymbol{q}$ or $\boldsymbol{q}'$) and $s_p = \pm 1$. Similarly to (3.5a,b), the four amplitude equations are

(4.1)\begin{equation} \left. \begin{array}{c} \partial_T B_{Q1}^{s_k} = \displaystyle \sum_{s_p ={\pm} 1} C_{Q K p}^{ s_k s_k s_p} \, B_{p1}^{s_p\,*} \, \textrm{e}^{\textrm{i}\, \frac{{\rm \Delta} \omega_{KpQ}}{(kRo)^2} T} , \\ \partial_T B_{p1}^{s_p} = \displaystyle \sum_{\boldsymbol{K} ={\pm} \boldsymbol{k}} \dfrac{C_{p K Q}^{ s_p s_k s_k}}{(kRo)^2} \, B_{Q1}^{s_k\,*} \, \textrm{e}^{\textrm{i}\, \frac{{\rm \Delta} \omega_{KpQ}}{(kRo)^2} T} . \end{array} \right\} \end{equation}

The characteristic polynomial $\mathcal {P}(\sigma )$ of this set of four linear differential equations in terms of the Rossby number $Ro$ and the detunings is

(4.2)\begin{equation} \mathcal{P}(\sigma) = \left( \sigma^2 - \textrm{i} \,{\rm \Delta} \omega_{kpq} \sigma - Ro^2 S_k \right) \left( \sigma^2 - \textrm{i} \,{\rm \Delta} \omega_{{-}kpq'} \sigma - Ro^2 S_{{-}k} \right) - Ro^4 P_0, \end{equation}

with $S_{K}$ and $P_0$ two real coefficients defined as

(4.3a,b)\begin{align} S_K = \sum_{s_p} C_{pKQ}^{s_p s_k s_k} C_{QKp}^{s_k s_k s_p\, *} \quad\mbox{and}\quad P_0 = \left(\sum_{s_p} C_{q k p}^{ s_k s_k s_p} C_{p -k q'}^{ s_p s_k s_k} \right) \left(\sum_{s_p} C_{q' -k p}^{ s_k s_k s_p} C_{p k q}^{ s_p s_{{+}k} s_k}\right). \end{align}

In general, the growth rate $\sigma _k(\,\boldsymbol{p};Ro)$ is found by numerical computation of the roots of $\mathcal {P}$. Nevertheless, an estimate of the maximum growth rate can be obtained in the limit $\vert p_x \vert \ll \vert p_y \vert$, which is relevant to our DNS. In this case, symmetries impose ${\rm \Delta} \omega _{kpq} = -{\rm \Delta} \omega _{-kpq'}$, $S_k = S_{-k}$ and $P_0 = S_{k}^2$. The polynomial $\mathcal {P}$ then has a purely real root,

(4.4)\begin{equation} \sigma_k (\,\boldsymbol{p}; Ro) = \sqrt{2 Ro^2 S_k (\,\boldsymbol{p}) - {\rm \Delta} \omega_{kpq}^2} \end{equation}

and the frequency of the growing geostrophic mode is $0$. We show in figure 4(b) the excellent agreement between this theoretical law and the DNS data of figure 2(a). An expansion similar to § 3.1 reveals that the maximum growth rate follows the exact same $(kRo)^2$ law as in the case of the single triad (3.10). This is confirmed by the systematic computation of the growth rate in the double triad case, as shown in figure4(c). Note that although this instability is driven by two imposed modes, it remains different from the four-mode interaction mechanism detailed in Brunet etal. (2020). The latter features intermediate, non-resonant modes that are absent in the double triad mechanism. Moreover, the near-resonant growth rate (4.4) is proportional to $Ro$ when $Ro$ is sufficiently large, whereas it always follows an $Ro^2$ law in the case of four-mode interaction.

We compare in figure 4(c) our theoretical predictions with the geostrophic growth rate extracted from DNS initiated with a large-scale noise, as in § 2, down to $Ro = 5 \times 10^{-3}$. We use the same imposed wavevector as previously ($\boldsymbol {k} = 2 {\rm \pi} [4,0,8 ]$) but also 1.5 and 2 times longer wavevectors to confirm that the growth rate is a function of the intrinsic Rossby number $k Ro$. Such a dependence is expected since in the infinite-domain limit, $L \rightarrow \infty$, $L$ becomes irrelevant and $1/k$ is the only remaining length scale, and the associated Rossby number is $k Ro$. As shown in figure 4(c), the numerical growth rate coincides with the maximum near-resonant growth rate in the low $Ro$ regime but also in the moderate $Ro$ regime where it is proportional to $kRo$. Note that for the three lowest $k Ro$ points, the noise has been implemented on two-dimensional modes ($\,p_z = 0$) to facilitate the isolation of the geostrophic instability. It delays the growth of two-dimensional modes with non-zero frequency by direct forcing at the lowest values of $Ro$. Despite their rapid growth, the latter are subdominant in the saturation of wave-driven flows, and do not prevent the long-term growth of unstable geostrophic modes under the mechanism examined here.

Lastly, our analysis allows us to understand the transition in the stability of the geostrophic flow observed in the numerical study (see figure 1). Here, $Ro = 2.83 \times 10^{-2}$ ($k Ro \simeq 1.6$) corresponds to the transition zone where the geostrophic growth rate is $O(kRo)$, as the wave-only triadic resonances (see figure 4c). When the Rossby number is decreased by one order of magnitude, as shown in figure 4(c), the growth rate of geostrophic modes scales as $(kRo)^2$, i.e. a factor $k Ro$ smaller than the growth rate of wave-only triadic resonances. This is why the geostrophic instability is not observed in the simulation at $Ro = 2.83 \times 10^{-3}$.

5. Conclusion

By means of numerical simulations and theoretical analysis, we have described a new instability mechanism by which inertial waves excite $z$-invariant geostrophic modes. We have proved that this instability is driven by near-resonant triadic interaction and derived its theoretical growth rate. When normalised by the global rotation rate $\Omega$, the growth rate follows a $(kRo)^2$ law and a $kRo$ law at small and moderate wave amplitude, respectively, $k$ being the imposed wavenumber. It translates into $\Omega ^{-1}$ and $\Omega ^0$ laws, respectively, for the dimensional growth rate. The near-resonant instability completes the picture proposed by Brunet etal. (2020) where another inviscid geostrophic instability based on two imposed modes and four-mode interaction is detailed. Although of different nature, both instabilities have a $(kRo)^2$ growth rate in the limit of small Rossby numbers which makes them possibly difficult to distinguish. On the one hand, in the linear growth phase, the eigenmode of the four-mode instability consists of waves and geostrophic flow of comparable amplitude, whereas the geostrophic flow is $kRo$ times smaller than the waves in the present mechanism. On the other hand, the near-resonant instability achieves a growth rate of order $kRo$ for $kRo \gtrsim 1$. Our analysis may thus explain the experimental results of Le Reun etal. (2019): they found at moderate Rossby number a geostrophic instability with a growth rate proportional to $Ro$, which thus matches the moderate Rossby number law derived here.

Our work paves the way for new studies dealing with the energy transfers from waves to geostrophic modes in rotating fluids, in particular in rotating turbulence. Firstly, it remains to be seen under which conditions the near-resonant instability may excite slow nearly geostrophic modes with small but non-zero frequencies. These modes have been proved to be important for the development of anisotropy in rotating turbulence (Smith & Lee 2005) even in the asymptotic limit of small Rossby numbers (van Kan & Alexakis 2019). With a heuristic analysis similar to the one developed in § 3.1 and appendix A, we can predict again a $(kRo)^2$ growth rate, but a detailed analysis is difficult and remains to be done. Lastly, the instability we describe in the present article could play an important role in the rotating turbulence dynamics that remains to be fully deciphered. Since wave–wave triadic interactions grow at rate $O(kRo)$ while wave–geostrophic interactions have a $O((kRo)^2)$ growth rate, we may speculate that, at sufficiently low forcing amplitudes, rotating turbulence may remain purely three-dimensional as the wave amplitudes would remain below the threshold of geostrophic instabilities. Such states have been observed by Le Reun etal. (2017), Le Reun etal. (2019) and Brunet etal. (2020), but a systematic investigation of their existence in various realisations of rotating turbulence remains to be carried out. In the future, our study could help bridge the gap between finite Rossby number experiments and simulations of turbulence (Godeferd & Moisy 2015) and asymptotic models of rotating turbulence (Galtier 2003; Bellet etal. 2006; van Kan & Alexakis 2019).