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

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 provided that the wave amplitude, or Rossby number $Ro$, is sufficiently large. We show this instability to be caused by triadic interaction involving two inertial waves and a geostrophic mode such that the sum of their eigen frequencies is non-zero. We derive theoretical scalings for the growth rate of this near-resonant instability which is proportional to $Ro^2$ at low $Ro$ and transitions to a $Ro$ law at moderate $Ro$. These scalings are in excellent agreement with direct numerical simulations. This instability could explain recent experimental observations of geostrophic instability driven by waves.


Introduction
Rotating turbulent flows are ubiquitous in geo-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 geostrophic modes, that is, columnar eddies parallel to the axis of rotation and corresponding to a first-order balance between the Coriolis force and the pressure gradient. Inertial waves, sustained by the restoring effect of the Coriolis force, have also been observed (Yarom & Sharon 2014;Clark di Leoni et al. 2014;Campagne et al. 2015). Recent numerical (Le Reun et al. 2017) and experimental (Le Reun et al. 2019;Brunet et al. 2019) 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). At larger forcing amplitudes, a secondary instability leads to the classical geostrophic-dominated turbulence. Hence, although bidimensionalisation in the form of geostrophic eddies has been commonly observed, it may not be the only equilibrium state of rotating turbulence, as pointed out for instance by Yokoyama & Takaoka (2017), Favier et al. (2019) and van Kan & Alexakis (2019).
The emergence of balanced geostrophic flows out of waves remains a challenging problem. Although the wave-to-wave interactions occur primarily through triadic resonance (Bordes et al. 2012), the latter 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-modes interactions can transfer energy from waves to geostrophic flows, either directly (Newell 1969;Smith & Waleffe 1999) or through an instability mechanism (Kerswell 1999;Brunet et al. 2019). The growth rate of such an instability scales like Ro 2 , with Ro the dimensionless wave amplitude or Rossby number. It develops over longer timescales that 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 has 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 transfer directly 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 that a single inertial wave may be unstable and drive exponential growth of geostrophic modes using direct numerical simulations (DNS) and theoretical analysis.

Governing equations and numerical methods
Let us consider an incompressible fluid rotating at rate Ωe z . We investigate the stability of a single plane inertial wave with wave vector k and eigen frequency ω k . Its amplitude is proportional to U w with h s k k is a helical mode, that is, when k is not parallel to the axis of rotation e z (Smith & Waleffe 1999;Bellet et al. 2006) where s k = ±1 is the sign of the helicity of the plane wave. If k is parallel to e z , h k = e x + is k e y . U w automatically satisfies the incompressibility condition since ∇ · U w ∝ k · h s k k = 0. U w satisfies the linearised rotating Euler equation, provided that ω s k k and k are related by the dispersion relation of inertial waves θ being the angle between the wavevector k and the rotation axis e z (ranging from 0 to π), and k = |k|. We solve for the time evolution of perturbations u to the wave U w maintained at a constant amplitude via the following set of equations where time is scaled by Ω −1 and length by the domain size L. We have introduced the Ekman number E = ν/(L 2 Ω), ν 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 {u, π} is determined with pseudo-spectral methods, that is, {u, π} is decomposed into a truncated sum of Fourier modes {û q ,π q } e iq·x . A wave-vector q writes 2π/L (n x , n y , n z ) where n x,y,z are integers varying from 0 to N . In the following, N = 96; higher resolutions have been tested and yield the exact same results. The temporal dynamics of the modesû 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 kRo rather than Ro alone.

Numerical results
Keeping the Ekman number to E = 10 −6 , two simulations of the stability of the inertial wave k = 2π [4, 0, 8] (with s k = 1) are performed at low (Ro = 2.83×10 −3 ) and moderate (Ro = 2.83 × 10 −2 ) wave amplitudes. They are both initiated with a random noise comprising wavevectors ranging between 0 and 15π, with very small initial amplitude. The use of spectral methods allows separating the kinetic energy of the perturbation flow u into a two-dimensional component E g , accounting for all modes q with q z = 0, and a complementary three-dimensional component E 3d . In addition, performing Fourier transform in space and time allows projecting the kinetic energy of u in the sub-space of the dispersion relation of inertial waves to draw spatio-temporal spectrum E(θ, ω) (Yarom & Sharon 2014;Le Reun et al. 2017). Note that, in these maps, θ is restricted to [0, π/2] since the flow is real and only the wavevectors q with q z 0 are simulated due to the Hermitian symmetry.
The kinetic energy time series and maps are shown for both simulations in figure 1. At low imposed wave amplitude (figure 1 a-b), three-dimensional perturbations dominate the growth of the instability. The energy map displays two spots aligned on the dispersion relation with frequencies ω 1,2 such that ω 1 + ω 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 twice larger than the rate of three-dimensional modes. The geostrophic modes are not unstable themselves, at least on the timescale of the growth of waves. Their growth is due to direct forcing by non-linear interaction of two growing waves involved in triadic resonances, with close frequencies and and opposed vertical wavenumbers, hence the twice larger growth rate. This mechanism corresponds to the direct excitation of geostrophic modes by two waves identified by Newell (1969) and Smith & Waleffe (1999).
At larger wave amplitude (figure 1c-d), this picture changes: two-and threedimensional 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 2a-b. The geostrophic unstable flow contains the wavevector p g = 2π [0, 5, 0] (and −p g for hermiticity), which grows along with other three-dimensional structures whose energy location in the (θ, ω) space is reminiscent of triadic resonant instability, but with significantly more spreading.
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, p = 2π [0, p y , 0] with p y ∈ {1, 3, 5} and s p = ±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, too large to be attributed to viscosity, hence suggesting an inviscid mechanism.
On the kinetic energy map E(θ, ω) (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 k and ±p g , 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 −k ± p (figure 2d) is outside its eigen frequency, indicating that the triad is not exactly resonant. Our results therefore point towards a near-resonant triadic instability involving geostrophic modes.
3. Near-resonance of geostrophic modes: theoretical approach 3.1. The low Rossby number limit In Cartesian domains with periodic boundary conditions, the total flow U = u+Ro U w is decomposed into a superposition of plane waves h sp The Euler equation governing U is then equivalent to a set of ordinary differential equations governing the amplitudes b sp p (Smith & Waleffe 1999): and ∆ω kpq ≡ ω s k k +ω sq q +ω sp p , (3.2) and where k, p = |k|, |p|. ∆ω kpq is the sum of the eigen frequencies 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 ∆ω kpq → 0. If all three modes k, p and q are inertial waves, this leads to the well known mechanism of triadic resonance (Bretherton 1964;Vanneste 2005;Bordes et al. 2012). When the mode k is imposed with an amplitude Ro and helicity s k , p and q grow exponentially with a rate proportional to Ro. This picture is changed when one of the modes, say p, is geostrophic (i.e. ω sp p = 0 and p z = 0, regardless of s p ). The spatial interaction condition k + p + q = 0 forces k z = −q z and the resonance condition becomes It may be fulfilled only when s q = s k and q − k → 0. In the governing equation forḃ sp p , the slow oscillation terms involve a coupling coefficient C sps k s k pkq ∝ k − q → 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 ∆ω 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 k with helicity s k is imposed with a small constant amplitude Ro, and that 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 τ = t and a slow time T . Because the wave-to-geostrophic transfer coefficient C sps k s k pkq vanishes as ∆ω kpq → 0, we find via a heuristic analysis that T = Ro 2 t, instead of Ro t for classical wave triads. In addition, it imposes the amplitude of p to be smaller by a factor Ro compared to the mode closing the triad q = −k − p. The amplitudes of the modes interacting with the imposed waves (p and q with both helicity signs s p,q ) are thus expanded as where the B j i are all O(1). Injecting the ansatz (3.4) into (3.1) and cancelling the fast secular growth of higher order terms gives the governing equations for the leading order coefficients. We find that triads may be unstable only when ∆ω kpq , and thus the wave-togeostrophic transfer coefficient C sps k s k pkq , are O(Ro 2 ). This condition may be fulfilled only for the modes p with s p = ±1 and q with s q = s k . Their amplitudes are then governed by the following equations, Ro 2 T , (3.5) the rescaled quantities C sps k s k pkq /Ro 2 and ∆ω kpq /Ro 2 being O(1). These equations have exponentially growing solutions and the complex growth rate of the instability is Note that the product of coupling coefficients C sps k s k pkq C s k s k sp * qkp is real. The expression of the real part of growth rate, σ k (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 ∆ω kpq and the coupling coefficient C sps k s k pkq are both O(Ro 2 ), the maximum geostrophic growth rate remains O(Ro 2 ) at most.
The Ro 2 scaling governing the maximum growth rate of geostrophic near-resonance is found quantitatively by expanding the frequency detuning, the transfer coefficients and then σ k (p; Ro) in the neighbourhood of exact resonance. As illustrated in figure 3a, the exactly resonant modes p 0 are located on a circle of centre −k ⊥ and radius k ⊥ , with p ⊥ = [k x , k y ]. We thus introduce p = p 0 + Ro 2 kδp where δp is a O(1) dimensionless vector in the geostrophic plane. Consider the closing modes q 0 = −(k+p 0 ) and q = −(k+p), then q = q 0 −Ro 2 kδp. The imposed mode k is left unperturbed. Since ∆ω kpq = C sps k s k p0kq0 = 0 at exact resonance, the leading order of the frequency detuning and the transfer coefficient are found from (3.2) and (3.3) to be O(Ro 2 ). The perturbation of the wavevectors is thus consistent with the asymptotic expansion. At leading order, where we have used that q 0 = k and s q = s k . In the growth rate σ k (p; Ro), the product of the coupling coefficient is Therefore, at leading order in powers of Ro, the growth rate is where X ≡ (q 0 · δp)/k and C k (p 0 ) is the sum in the right hand side of (3.8). When C k (p 0 ) > 0, the growth rate reaches an optimum at X = C k k 2 /(2ω 2 k ) with value k 2 C k (p 0 )/|4ω k |, which remains to be maximised over all exactly resonant wavevectors p 0 . The coefficient C k (p 0 ) is shown in figure 3b for several wavevectors k with different frequencies ω k . When plotted against p 0 /k ⊥ , all the curves C k collapse on a master curve that reaches a maximum value of 1 at p 0 → 0, regardless of the helicity sign s k of the imposed wave. In the low Rossby number limit, the maximal growth rate is then (3.10) We retrieve that kRo, 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 (kRo) 2 upper envelope. More details are given below in section 3.3 where we compare the law (3.10) to exact computation of the growth rate curves from (3.6).

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 a Ro law. This is proven directly by multiplying (2.5) by u and integrating over the fluid domain V , thus giving: where · n denotes the L n -norm. In virtue of Hölder's inequality (Gallet 2015), This upper bound applies in particular to the low Ro scaling (3.10), which thus holds up to kRo ∼ 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.

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 k = 2π [4, 0, 8], as functions of the Rossby number. We notice that each growth rate curve follows a Ro scaling at sufficiently large Ro, which is consistent with the asymptotic expansion carried out in section 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 a Ro 2 law that perfectly matches the theoretical prediction (3.10). This instability is thus fundamentally different from fourmodes interactions for which each growth rate curve follows a Ro 2 law (Kerswell 1999;Brunet et al. 2019). Besides, in agreement with the upper bound (3.12), we observe the Ro 2 maximum growth rate law to break down above Ro ∼ 1. Beyond that, the nearresonance growth rate derived from (3.6) remains below the upper bound and follows a Ro law. However, additional instabilities (shear,centrifugal, etc.) may also drive the growth of geostrophic modes in this regime.

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 ∆ω kpq as Ro → 0. Hence, discretisation implies the existence of a finite value of Ro below which the near-resonant instability vanishes. Regarding viscous effects, at finite but low Ekman number, i.e. when (kRo) 2 k 2 E, the low Rossby number scaling is unaltered since the near-resonant modes p and q have at most similar wavenumbers to k. The O(Ro) upper bound on the growth rate is also unaltered by the inclusion of viscosity.

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 figure  2a and (3.6) predicts a frequency ∆ω kpq /2 = 0 for the geostrophic modes, while it is 0 according to 2b. We must include in our model that not only k, but also −k (with the same helicity sign s) is imposed due to Hermitian symmetry. Consider two triads 2)) for wavevectors p restricted to |px| < 2π and |py| < 12π. 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 k + p + q = 0 and −k + p + q = 0 (see figure 4) with helicity signs s q = s q = s and s p = ±1. As shown schematically in figure 4a, the exact resonance circles associated to ±k coincide at p = 0. In a neighbourhood of this point (the dark grey intersection in figure 4a), a geostrophic mode q is in near-resonance with both ±k with both detunings ∆ω kpq and ∆ω −kpq small. This is all the more important that p = 0 corresponds to the largest near-resonant growth rate as Ro → 0.
We thus consider two triads k + p + q = 0 and −k + p + q = 0 (see figure 4) with helicity signs s q = s q = s k and s p = ±1. In a very similar way to the two-time asymptotic expansion of section 3.1, we can retrieve the amplitude equations governing the slowly varying envelopes B sp p1 and B s Q1 , with K = ±k, Q = −K − p, and s p = ±1. Similarly to (3.5), the four amplitude equations are Ro 2 T . (4.1) The characteristic polynomial P(σ) of this set of four linear differential equations in terms of the Rossby number Ro and the detunings ∆ω K is: with S K and P 0 two real coefficients defined as (4.3) In general, the growth rate σ k (p; Ro) is found by numerical computation of the roots of P. Nevertheless, an estimate of the maximum growth rate can be obtained in the limit |p x | |p y |, which is relevant to our DNS. In this case, symmetries impose ∆ω kpq = −∆ω −kpq , S k = S −k and P 0 = S 2 k . The polynomial P then has a purely real root, σ k (p; Ro) = 2Ro 2 S k (p) − ∆ω 2 kpq (4.4) and the frequency of the growing geostrophic mode is then 0. We show in figure 4b the excellent agreement between this theoretical law and the DNS data of figure 2a. An expansion similar to section 3.1 reveals that the maximum growth rate follows the exact same 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 figure 4c. We compare our theoretical predictions with the geostrophic growth rate extracted from DNS initiated with a large-scale noise, as in section 2, down to Ro = 5 × 10 −3 As shown in figure 4c, the numerical growth rate coincides with the maximum nearresonant growth rate, both in the small and moderate Rossby numbers regimes. For the two lowest 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 wavedriven 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). Ro = 2.83 × 10 −2 (kRo 1.6) corresponds to the transition zone where the geostrophic growth rate is O(Ro), as the wave-only triadic resonances (see figure 4c). Below, the latter grow much faster than the geostrophic near-resonant instability and appear to dominate.

Conclusion
We have described a geostrophic instability excited by a single inertial wave, by the means of numerical simulations and theoretical analysis. We have proved that this is instability is driven by near-resonant triadic interaction and derived its theoretical growth rate. It follows a (kRo) 2 and a kRo laws at respectively small and moderate wave amplitude Ro, k being the wavenumber. The near-resonant instability completes the picture proposed in Brunet et al. (2019) where another inviscid geostrophic instability based on two imposed modes and four-mode interaction is detailed. Although of different nature, both instabilities have a Ro 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 Ro times smaller than the waves in the present mechanism. On the other hand, the near-resonant instability achieves a growth rate of order Ro for moderately low Ro. Our analysis may thus explain the experimental results of Le Reun et al. (2019): they found at moderate Rossby number a geostrophic instability with a Ro growth rate, which thus matches the moderate Ro law derived here. Our work also paves the way for new studies in rotating turbulence, in particular regarding the conditions for the bi-dimensionalisation of these flows. Since wave-wave triadic interactions grow at a rate O(Ro) while wave-geostrophic interaction have a O(Ro 2 ) growth rate, we may speculate that, at sufficiently low forcing amplitude, 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 et al. (2017Reun et al. ( , 2019 and Brunet et al. (2019), but a systematic investigation of their existence in various realisations of rotating turbulence remains to be carried out. These studies will hopefully be facilitated by our quantification of wave-driven geostrophic instabilities.