Interplay between geostrophic vortices and inertial waves in precession-driven turbulence

Rotating turbulence is an ubiquitous phenomenon in a broad context ranging from astrophysical and geophysical flows^1–4 to industrial applications.^5,6 Understanding the impact of rotation on the turbulence dynamics is far from trivial due to the complexity of the nonlinear processes involved. In general, when a fluid is subjected to rotational motion, the nonlinear interactions are affected by the Coriolis force whose strength is quantified by the Rossby number, Ro (the ratio of the advection rate to the rotation rate), and the Reynolds number, Re (the ratio of the advection rate to the viscous dissipation rate). If the Coriolis force is strong enough, the formation of coherent columnar vortices occurs inside the fluid flow. This phenomenon has been observed in experimental campaigns for several systems such as oscillating grids,⁷ decaying turbulence,^8,9 forced turbulence,^10–12 and turbulent convection.^13,14 Also, numerical simulations have been instrumental to analyze such tendency in a myriad of cases,^15–21 making use of large-eddy simulations,^22–24 and even turbulence models.^25,26

The emergence of columnar vortices aligned along the flow rotation axis is accompanied by inertial waves, which are inherent to rotating fluids. Their frequency magnitude ranges between zero and twice the rotation rate Ω of the objects.²⁷ The dynamics and mutual couplings between these two basic types of modes largely depend on Re and Ro. The results of the asymptotic analysis at $Ro≪1$ indicate that three-dimensional (3D) inertial waves and two-dimensional (2D) vortices are essentially decoupled and evolve independently: vortices undergo inverse cascade, while the wave energy cascades forward through resonant wave interactions in the regime of weakly nonlinear inertial wave turbulence.^28–31 However, this picture does not carry over to moderate Rossby numbers $Ro∼0.1$⁠, where the situation is much more complex, since 3D inertial waves (the so-called fast modes) and 2D vortices (also called slow modes) can coexist and be dynamically coupled. In this case, asymptotic analysis cannot be used and more complex mathematical models have been proposed to explain the geostrophic vortices–wave interaction, such as the quartic instability³² or near-resonant instability.³³ However, progress can be made mainly by numerical simulations. Several works are devoted to the study of these two manifolds and their interactions for forced rotating turbulence^29,34–36 and also for convective and rotating turbulence.^37,38

Indeed, other forcing mechanisms have been shown to be characterized by this interplay of vortices and waves, such as elliptical instabilities^39–43 and tidal forcing.^44–48 In this respect, the precession-driven dynamics represents a possible candidate for the development of both 3D waves with embedded 2D vortices,^36,49 but so far these studies did not investigate a wide range of governing parameters. Other works were devoted to the stability analysis of the precession flows.⁵⁰ 

The modified local Cartesian model of a precession-driven flow was proposed by Mason and Kerswell⁵¹ and later used by Barker⁴⁵ to study its nonlinear evolution. In the first paper, rigid and stress-free axial boundaries in the vertical direction were used, while in the second paper an unbounded precessional flow was considered in the planetary context, employing the decomposition of perturbations into shearing waves. In this paper, we follow primarily the approach of Barker,⁴⁵ who analyzed the occurrence of vortices, function of the precession parameter (Poincaré number), including energy spectrum and dissipation properties. The main advantage of the local model is that it allows high-resolution study of linear and nonlinear dynamical processes in precession-driven flows, which is much more challenging in global models. Also, this model allows to focus only on the dynamics of the bulk flow itself avoiding the complications due to boundary layers. This is important for gaining a deeper understanding of perturbation evolution in unbounded precessional flows and then, comparing with the global simulations, for pinning down specific effects caused by boundaries.

In this paper, we continue this path and investigate in detail the underlying dynamical processes in the turbulence of precessional flow in the local model. We decompose perturbations into 2D and 3D manifolds and analyze their dynamics and interplay in Fourier space. Our main goal is to address and clarify several key questions: (i) how the presence and properties of columnar vortices depend on the precession strength and the Reynolds number (here defined as the inverse Ekman number), (ii) what are the mechanisms for the formation of the columnar vortices in precession-driven flows, in particular, how their dynamics are affected by precessional instability of inertial waves, that is, if there are effective nonlinear transfers (NL) (coupling) between vortices and the waves; (iii) what are the dominant nonlinear processes (channels) in this vortex–wave system, that is, the interaction of 2D–3D modes or 2D–2D modes (inverse cascade); (iv) in terms of total shell-average spectral analysis, what type of cascades (inverse, forward) occur and what kind of spectra characterize precessional flows.

The paper is organized as follows: In Sec. II, the local model and governing equations in physical and Fourier space are presented, and numerical methods introduced. Section III presents general evolution of the volume-averaged kinetic energy and dynamical terms as well as flow structure. In Sec. II B, we investigate the nonlinear dynamics of 2D vortices and 3D inertial waves and nonlinear interaction between them in Fourier space. In this section, we also characterize turbulent dissipation as a function of precession parameter Po. Discussions and the future perspectives are presented in Sec. V.

We consider a precessional flow in a local rotating Cartesian coordinate frame (also referred to as the “mantle frame” of a precessing planet) in which the mean total angular velocity of fluid rotation $Ω=Ωez$ is directed along the z axis. In this frame, the equations of motion for an incompressible viscous fluid take the form (see a detailed derivation in Refs. 45 and 51):

where U is the velocity in this frame, ρ is the spatially uniform density, and P is the modified pressure equal to the sum of thermal pressure and the centrifugal potential. The last two terms on the left-hand side in the brackets are the Coriolis and the Poincaré forces, respectively, and $ϵ(t)=Po(cos(Ωt),−sin(Ωt),0)T$ is the precession vector with Po being the Poincaré number characterizing the strength of the precession force. The last term on the right-hand side (RHS) is the second part of the precession force with vertical shear, which is the main cause of hydrodynamic instability in the system, referred to as the precessional instability.⁵² ν is the constant kinematic viscosity.

The basic precessional shear flow in this local frame represents an unbounded horizontal flow with a linear shear along the vertical z axis and oscillating in time t, that is, $Ub=(Ubx(z,t),Uby(z,t),0)$ with the components given by^45,51

where $r=(x,y,z)$ is the local position vector. Our local model deals with perturbations to this basic flow, $u=U−Ub$⁠, for which from Eq. (1) we obtain the governing equation:

where the last three terms on the RHS are related to precession and proportional to Po. The flow field is confined in a cubic box with the same length L in each direction, $Lx=Ly=Lz=L$ (Fig. 1). In other words, both horizontal and vertical aspect ratios of the box are chosen to be equal to one in this paper. Varying these aspect ratios affects the linearly unstable modes that can be excited in the flow and the properties of the vortices.³⁹ 

Below we use the non-dimensionalization of the variables by taking $Ω−1$ as the unit of time, box size L as the unit of length, $ΩL$ as the unit of velocity, $ρL2Ω2$ as the unit of pressure, and perturbation kinetic energy density $E=ρu2/2$⁠. Thus, in our units, $Ω=L=ρ=1$ and every variable is normalized in terms of them. The key parameters governing a precession-driven flow are the Reynolds number (inverse Ekman) defined as

and the Poincaré number Po introduced above. Note that this definition of Re differs from that typically used in forced turbulence theory.^29,31

Our main goal is to perform the spectral analysis of precession-driven turbulence in Fourier (wavenumber k-) space in order to understand dynamical processes (energy injection and nonlinear transfers) underlying its sustenance and evolution. To this end, following Refs. 39 and 45, we decompose the perturbations into spatial Fourier modes (shearing waves) with time-dependent wavevectors $k(t)$⁠,

where $f≡(u,P)$⁠, and the Fourier transforms of these fields are $f¯≡(u¯,P¯)$⁠. Here, the summation is taken over a discrete set of wavenumbers in the cubic box, $ki=niΔki$⁠, where $i=(x,y,z)$ and $Δki=2π/Li=2π/L$ is the grid cell size in Fourier space, which is the same in each direction, with integer $ni=0,±1,±2,…,±Ni/2$⁠. In decomposition (4), the wavevector of modes oscillates in time,

about its constant average value $⟨k(t)⟩=(kx0,ky0,kz0)$ due to the periodic time variation of the basic precessional flow $Ub$⁠. Substituting Eq. (4) into Eq. (3) and taking into account the above non-dimensionalization, we obtain the following equations governing the evolution of velocity amplitude

Note that the wavevector $k(t)$ as given by expression (5) satisfies the ordinary differential equation

and as a result, the last term on the RHS of Eq. (3) related to the basic flow has disappeared when substituting (4) into it. The term $Q(k,t)$ on the RHS of Eq. (6) represents the Fourier transform of the nonlinear advection term $u·∇u=∇·(uu)$ in the original Eq. (3) and is given by the convolution^35,53

where the indices $(m,n)=(x,y,z)$⁠. This term describes the net effect of nonlinear triadic interactions (transfers) among a mode k with two others $k−k′$ and $k′$ and thus plays a key role in turbulence dynamics.

Multiplying both sides of Eq. (6) by the complex conjugate of spectral velocity $u¯∗$⁠, the contribution from Coriolis and part of the Poincaré force in the total kinetic energy of a mode cancel out, since they do not do any work on the flow, $u¯∗·(2ez×u¯−2ε(t)×u¯)=0$⁠, and as a result, we obtain the equation for the (non-dimensional) spectral kinetic energy density $E=|u¯|2/2$ in Fourier space as

The pressure term also cancels out since $u¯∗·k(t)P¯=0$⁠. Thus, the RHS of Eq. (10) contains three main terms:

• Injection

  $A≡−12[u¯∗(Mu¯)+u¯(Mu¯)∗],$

  which is of linear origin, being determined by the matrix M, that is, by the precessing background flow and describes energy exchange between the perturbations and that flow. If A > 0, kinetic energy is injected from the flow into inertial wave modes, and hence, they grow, which is basically due to precessional instability,^45,51,52,54 whereas at A < 0 modes give energy to the flow and decay.

• Nonlinear transfer

  $NL≡12[u¯∗Q+u¯Q∗]$
  describes transfer (cascade) of spectral kinetic energy among modes with different wavenumbers in Fourier space due to nonlinearity. The net effect of this term in the spectral energy budget summed over all wavenumbers is zero, that is,

  $∑kNL(k,t)=0,$

  which follows from vanishing of the nonlinear advection term in the total kinetic energy equation integrated in physical space. Thus, the main effect of the nonlinear term is only to redistribute energy among modes that is injected from the basic flow due to A, while keeping the total spectral kinetic energy summed over all wavenumbers unchanged. Although the nonlinear transfers NL produce no net energy for perturbations, they play a central role in the turbulence dynamics together with the injection term A. The latter is thus the only source of new energy for perturbations drawn from the infinite reservoir of the background precessional flow. Due to this, below we focus on these two main dynamical terms—linear injection and nonlinear transfer functions, compute their spectra, and analyze how they operate in Fourier space in the presence of precession instability using the tools of Refs. 53 and 55.

• Viscous dissipation

  $D≡−2k2ReE$

  is negative definite and describes the dissipation of kinetic energy due to viscosity.

In the present section, we follow a widely used approach in the theory of rotating anisotropic turbulence^{29,31,35–37,39} and decompose the flow field into 2D and 3D modes in Fourier space to better characterize this anisotropy between horizontal and vertical motions. This choice is motivated by the observation of two main types of perturbations: vortices, which are essentially 2D structures, and 3D inertial waves in rotating turbulent flows with external forcing such as libration, elliptical instability,^39,42 precession,^36,45 and other artificial types of forcing concentrated at a particular wavenumber.^34,35,56 The 2D vortical modes, also called slow (geostrophic) modes, have dominant horizontal velocity over the vertical one and are almost uniform, or aligned along the $z$ axis; that is, their wavenumber parallel to this axis is zero k_z = 0. This slow manifold is also referred to as 2D and three-component (2D3C) field in the literature, since it varies only in the horizontal (x, y)-plane perpendicular to the rotation axis, but still involves all three components of velocity with the horizontal one being dominant. On the other hand, 3D inertial wave modes, called fast (with nonzero frequency $ω=±2Ωkz/k$⁠) modes, have comparable horizontal and vertical velocities and vary along z axis; that is, parallel wavenumber is nonzero $kz≠0$⁠.⁵⁷ Specifically, these two mode manifolds are

and the spectral velocities for 2D vortices and 3D inertial waves can be defined as

Note that the definition of 2D manifolds here differs from the Taylor–Proudman problem since it does not necessarily have vanishing vertical flow. Indeed, velocity for both the 3D and 2D modes can be decomposed in turn into respective horizontal $u¯h=(u¯x,u¯y,0)$ and vertical $u¯z$ components.

Using the general Eq. (10), separate equations can be written for 2D and 3D mode spectral energies defined, respectively, as $E2D=|u¯2D|2/2$ and $E3D=|u¯3D|2/2$⁠,^35,39,58

Since injection A and dissipation D terms are of linear origin, they act for 2D and 3D modes separately, that is,

However, the nonlinear transfers can act only among modes, which lie, respectively, within the slow or the fast manifold, that is, nonlinear interactions separately among 2D–2D modes (vortex–vortex),

and among 3D–3D modes (wave–wave),

(long bars denote Fourier transforms) as well as between these two manifolds, that is, nonlinear cross interactions/couplings between 2D and 3D modes (vortex–wave),^35,59

In this case, in a triad, 2D modes can receive/lose energy via nonlinear interaction of two 3D modes with opposite signs of k_z, while 3D modes can receive/lose energy via interaction of only 2D and another 3D mode (interaction between 2D and 2D modes obviously cannot feed 3D modes).

The 2D–2D and 3D–3D nonlinear transfers are each conservative, $∑kNL2D2D=0$ and $∑kNL3D3D=0$⁠, whereas the cross-transfer terms $NL2D3D$ and $NL3D2D$ are not, but their sum is conservative, as the net effect of these terms summed over all wavenumbers, as it should be, are equal in absolute value but have opposite signs, that is,

Below we analyze the action of these injection and transfer terms for different Po and Re (Table I).

In the following, we will mostly use shell averages of these spectra, which are defined in the standard way as³¹ 

for each spectral quantity $f≡(E,A,NL,D)$⁠. For 3D modes $k∈Ψ3D$⁠, the summation is done over spherical shells with radius $k=(kx2+ky2+kz2)1/2$ and width $Δk=2π/L$⁠, while for 2D modes, having $k∈Ψ2D$⁠, the summation is over rings in the $(kx,ky)−$ plane with radius $k=(kx2+ky2)1/2$ and the same width $Δk$⁠. When we plot spectra for 2D and 3D quantities, we implicitly assume each depends on its respective wavenumber magnitude k and ignore any anisotropy in each spectrum.

We solve Eq. (3) using the pseudo-spectral code SNOOPY,^60,61 which is a general-purpose code, solving hydrodynamic (HD) and magnetohydrodynamic (MHD) equations, including shear, rotation, weak compressibility, and several other physical effects. The Fourier transforms are computed using the FFTW3 library. Nonlinear terms are computed using a pseudospectral algorithm with antialiasing 3/2 rule. The original version of the code has been modified⁴⁵ to include precessional forcing, and hence, variables are decomposed in terms of shearing waves with periodically time-varying wavevector (Eq. 5) due to the basic shear flow $Ub$ induced by precession. In this way, the shearing-periodic boundary conditions in the local domain (which are in fact fully periodic in the frame co-moving with the basic flow) are naturally satisfied in the code. We stress that the inclusion of precession-related terms in the code has been already tested through a linear stability analysis.⁴⁵ Indeed, the last three terms on the RHS of Eq. (3) are linear; therefore, the comparison against the growth rate of precession instability predicted by linear stability theory^52,54 (as in fact done in Sec. 4.2 of Barker,⁴⁵ see his Fig. 1) is a sufficient test.

The resolution of the various simulations is listed in Table I. As it is seen, for the largest Re and Po, we accordingly increase the number of grid points. This type of resolutions is consistent with the ones typically used in forced rotating turbulence studies.^36,45,62 For further details regarding the numerical accuracy, see  Appendix A.

The simulations performed in this paper for different pairs of (Po, Re) are listed in Table I. All the runs start with small solenoidal random noise perturbations of the velocity on top of the basic flow composed of wavenumbers $k/2π∈[1,21]$ with random phases and the same amplitude for each k, so that in physical space this yields the root mean square (RMS) of the perturbed velocity $⟨u2⟩1/2=1.12×10−4$⁠. We have checked that the results do not depend on different realizations of the random noise, because shortly after the evolution starts the most unstable mode emerges and determines subsequent time development.

In Fig. 2, we plot the time evolution of the volume-averaged kinetic energy, which is equal to the sum of energies of 2D and 3D modes over all wavenumbers, $⟨E⟩=∑kE=∑k(E2D+E3D)$⁠, for several precession parameters Po and at an intermediate Reynolds number $Re=104.5$ (see  Appendix B for the effect of varying Re). In the early linear regime, the energy grows exponentially corresponding to dominant 3D inertial waves being excited first by the precessional instability^45,51,52,54 (see inset in Fig. 2). In the given range of Po, the growth rate of the precession instability increases with Po. After about several hundreds of orbital times, the exponential growth attains sufficient amplitude for nonlinearity to come into play and cause the instability to saturate with higher amplitudes and shorter saturation times for larger Po. Depending on Po, the saturated states are qualitatively different, exhibiting statistically steady turbulence at higher $Po≳0.1$ with only mild fluctuations, whereas strong quasi-periodic oscillations (bursts) are observed at small $Po≲0.1$ with more than an order of magnitude variations in the kinetic energy. This temporal behavior of the volume-averaged kinetic energy in the nonlinear state of the precessional instability with Po is consistent with previous related local studies.^36,45 Below we interpret this temporal evolution of the total kinetic energy in terms of the dynamics of 2D vortices and 3D waves and their interplay.

A first comparison between the 2D and 3D mode dynamics is shown in Fig. 3 where the evolution of the total kinetic energies for 2D modes, $⟨E2D⟩$⁠, and 3D modes, $⟨E3D⟩$⁠, is plotted for three precession parameters. For larger $Po≳0.1$ the energy of 2D modes is more than one order of magnitude larger than that of 3D modes. However, the saturated value of $⟨E3D⟩$ tends to increase more than that of $⟨E2D⟩$ with increasing Po, implying that the waves, as they should be, are more affected and intensified by precession rather than the vortices. The 2D vortices are linearly stable against precession instability and hence cannot grow due to the latter.⁵² They are driven and energetically supplied by waves via nonlinear transfers,⁴⁵ which will be examined in detail below using the spectral analysis. On the other hand, for the lowest precession parameter Po = 0.075, corresponding to the bursty regime, the energy of 3D waves periodically dominates over the 2D vortical mode energy during the growth (burst) phase (referred to as State 1). In this burst phase, waves excited by the precessional instability lose their energy to 2D vortices due to nonlinearity. As a consequence, the energy of the waves drops by about an order of magnitude (affected additionally by viscous dissipation) relative to the 2D mode energy (State 2). After that it starts to increase again due to precessional instability, closing the cycle. Although 2D mode energy also decreases at this stage, it does so much slower, on a viscous time.⁴⁵ This cyclic behavior of both components is remarkable, indicating the quasi-periodic nature of evolution due to weak precessional forcing (⁠$Po≲0.1$⁠), which is relevant to astrophysical and geophysical regimes.^45,63,64 We explore this behavior in more detail in the spectral analysis section below.

Figure 4 shows the structure of the vertical component of vorticity, $ωz=(∇×u)z$⁠, in physical space well after the saturation for the above regimes of weak, moderate, and strong precessions. The top row shows the two different states in the Po = 0.075 case with bursts. In State 1, the energy of 3D waves dominates over that of 2D vortical modes (Fig. 3), and consequently, small-scale structures, varying along the $z$ axis, are more pronounced in the vorticity field. By contrast, in State 2, where the wave energy quickly decays afterward and 2D modes dominate, only vertically uniform columnar vortical structures aligned with the rotation z axis are present. At larger Po = 0.2 and Po = 0.3 shown, respectively, in bottom left and right panels of Fig. 4, the nonlinear states consist of vortices embedded in 3D waves, coexisting at all times. At medium Po = 0.2, two columnar vortices aligned with the rotation axis with opposite vorticity (cyclonic/anticyclonic) are clearly seen in the small-scale waves, whose strength with respect to vortices has increased compared to that in the above bursty regime. At even higher Po = 0.3, the contribution of 3D wave energy is somewhat larger (Fig. 3), and therefore, small-scale turbulent wave structures are more pronounced with respect to a single 2D vortex.

The regime diagram in Fig. 5 summarizes the properties of all the runs for different pairs (Po, Re). The colored dots represent the ratio of the time- and volume-averaged energy of 2D vortices, $⟨E2D⟩$⁠, to the total energy of all the modes, $⟨E⟩$⁠, in the statistically steady turbulent state, as shown in Figs. 2 and 3. The empty points represent the cases where the energy drops to negligible value meaning that the local flow in the box is stable against precessional instability. The colors show that at given Re, the fraction of 2D mode energy vs total energy initially increases with Po when $Po≲0.1$⁠, then reaches a maximum at medium precessions $Po∼0.1$ and decreases at larger $Po≳0.1$⁠. The maximum shifts toward smaller Po with increasing Re. We will carry out the analysis distinguishing the latter three groups in Sec. IV.

Having analyzed the time-development of the mode energies, next in Fig. 6 we plot the evolution of the volume-averaged (or equivalently integrated in Fourier space) dynamical terms in Eqs. (13) and (14), that is, the energy injection, $⟨A⟩=∑kA$⁠, for 3D waves (top) and 2D vortical (bottom) modes together with nonlinear transfer terms between them, $⟨NL⟩=∑kNL$ (middle). To obtain a better visualization and a clear trend of temporal evolution, we have smoothed these volume-averaged terms additionally over short time intervals around each time moment, thereby removing fast temporal oscillations and getting meaningful time averages. For 3D waves, energy injection occurs due to the precession instability and hence increases with increasing precession strength Po, while the 2D modes are excited by the nonlinear interaction among these waves via $NL3D2D$ term, which also increases with Po. This term is overall positive in time, implying transfer of energy from 3D waves to 2D vortices (see also Ref. 39). The excited 2D modes in turn also extract energy from the basic flow via positive $⟨A2D⟩>0$ term. Note that $⟨A2D⟩$ is not necessarily zero, unlike for the 2D linear modes in the background precessional flow. As a result, the evolution of $⟨A2D⟩$ is determined by the nonlinear transfer term $⟨NL2D3D⟩$ and hence follows the latter, as it is seen in Fig. 6 where the peaks of both these functions nearly coincide. For Po = 0.2, $⟨A2D⟩≈⟨A3D⟩$ (red), while for Po = 0.3, $A3D>A2D$ (orange curves), indicating that the precession instability feeds the waves, while the waves in turn feed vortices via nonlinear cascade. Below we will see how this process occurs scale by scale in Fourier space.

In Fig. 4, we have shown the total vorticity field including both 2D vortices and 3D waves. To better visualize these fields, we computed the inverse Fourier transforms from $u¯2D$ and $u¯3D$ and show respective energy densities in physical space in the saturated regime in Fig. 7. The left panels of this figure show energy of 2D modes, where now we can clearly distinguish two vortices for Po = 0.2 (top row) and a single vortex for Po = 0.3 (bottom row). The middle panels show the small-scale 3D mode energies in the $(x,y)−$ plane at the central height (z = 0) of the box. There is a noticeable difference between the Po = 0.2 and Po = 0.3 cases: for Po = 0.2, we observe larger-scale wave structures, whereas for Po = 0.3, the wave field is more fluctuating and rich in smaller scales, implying that increasing precession parameter intensifies first of all 3D waves and indirectly vortices due to their nonlinear coupling with the former. Note also that the 3D wave structures are concentrated around the vortices—a feature observed experimentally in a precessing spherical container.⁶⁵ The right panels show the vertical structure of 3D mode energy in the $(y,z)−$ plane at the center of vortices (located at $x=−0.33$ and $x=−0.2$⁠, respectively, for Po = 0.2 and 0.3, which are marked with a green dotted line in the middle row). Again, the Po = 0.3 case shows more fluctuating behavior with fine scales surrounding the column. So, the main dynamical picture consists of the coexisting columnar (geostrophic) vortices and waves whose magnitudes and length scales depend on the precession strength.

So far the study has been mainly conducted in the physical (coordinate) space. However, a deeper insight into the precession-driven turbulence dynamics can be gained by investigating the dynamical processes—energy injection, nonlinear transfers, and viscous dissipation—in Fourier space, where a much richer dynamical picture unfolds and becomes more accessible to analysis. Following the approach of Refs. 35, 53, and 55, we individually compute and visualize the injection A, viscous dissipation D, and various nonlinear transfer NL terms entering spectral energy Eqs. (13) and (14) in Fourier space using the simulation data, and analyze their interplay in different regimes with respect to the precession parameter identified above.

As we have seen above, the case with weak precession forcing is characterized by cyclic bursts, where the system alternates between State 1 and State 2 and the purpose of this section is to understand the underlying mechanisms of this behavior. With this goal, we analyze and compare the dynamics in two different intervals shown in Fig. 3: $Δt1$ corresponding to State 1, when the energy of 3D wave modes initially increases due to precessional instability, while the energy of 2D modes is still at its minimum, and $Δt2$ corresponding to State 2, when both 3D and 2D mode energies drop.

The shell-averaged spectra of the linear injection A and dissipation D terms for 2D and 3D modes as well as nonlinear transfer terms for 2D–2D, $NL2D2D$⁠, for 2D–3D, $NL2D3D$ and $NL3D2D$⁠, and for 3D–3D $NL3D3D$ mode interactions in these two states (also averaged over $Δt1$ and $Δt2$ time intervals, respectively) are shown in Fig. 8. The basic dynamical picture in this regime is the following. In State 1 (top row), the most important contribution is due to $A3D$ (blue), which injects energy into waves from the basic flow due to the initial development of precessional instability. This reaches a large peak at $k=8π$ whose value is positive and larger than the comparable effects of viscous dissipation $D3D$ (black) and transfer $NL2D3D$ (red), which are both negative reaching a minimum at the same wavenumber. The effect of nonlinear transfers among waves, $NL3D3D$ (green), is relatively small at this time. This also implies that the viscosity is already important at the injection scale; that is, there is not a good scale separation (inertial range) between the injection and viscous scales. Nevertheless, $A3D>0$ is sufficiently large to overcome both these negative (sink) terms and give rise to wave growth in State 1. Since $NL2D3D<0$ at the injection wavenumbers, its counterpart for 2D modes $NL3D2D>0$⁠, indicating that the waves nonlinearly transfer their energy to and amplify 2D vortices, but at lower wavenumbers near the peak $k=2π$ of this term. These vortices additionally receive some energy from the basic flow due to the positive $A2D$ (blue) term. However, the dissipation $D2D<0$ (black curve in top left panel) for vortices is quite high, prevailing over the positive $NL3D2D$ (red) and $A2D$⁠, and as a result, vortices do not yet grow at these times.

The nonlinear transfers between waves and vortices, $NL3D2D$ and $NL2D3D$⁠, increase by absolute value (but retain their signs) with time and already in State 2 the mostly negative $NL2D3D$⁠, together with dissipation $D3D<0$⁠, dominates the positive injection $A3D$ (bottom right plot in 8). As a result, wave energy quickly drops in State 2 (see also Fig. 3). On the other hand, the 2D vortices, which now receive much more energy from waves via the term $NL3D2D>0$⁠, also develop an inverse cascade themselves described by $NL2D2D$ (bottom left panel). This is manifested in the emergence of large-scale vortices in physical space in State 2 (top right plot of Fig. 4). The injection $A2D$ is relatively small/insignificant at these times. However, dissipation $D2D$ is still larger than the nonlinear replenishment by $NL3D2D$⁠, and consequently, the energy of vortices slowly decreases too (which is consistent with Fig. 3). Once vortices have become weak enough, the waves can grow again due to the precessional instability and close the cycle loop. Thus, we can conclude that the bursts are caused by a quasi-periodic behavior of 3D dynamical terms which in State 1 with $∑k(A3D+NL3D+D3D)>0$ leading to wave energy amplification, whereas in State 2 with $∑k(A3D+NL3D+D3D)<0$ leading to energy decay.

In this subsection, we present a similar analysis of the dynamical processes in Fourier space, focusing on the regime of large Po = 0.2 and 0.3 where the saturated state is characterized by a quasi-steady turbulence, where both 3D and 2D mode energies evolve in time with only mild fluctuations in contrast to the small-Po regime displaying quasi-periodic bursts (Fig. 3). In this section, we keep the Reynolds number fixed at $Re=104.5$ to focus on the impact of increasing precession on the spectral properties and dynamical balances of the turbulence.

Figure 9 shows the shell-averaged kinetic energy spectra of 2D vortices and 3D waves divided further into horizontal, $Eh=(|u¯x|2+|u¯y|2)/2$⁠, and vertical, $Ez=|u¯z|2/2$⁠, components at Po = 0.2 (left panel) and 0.3 (right panel). The time average has been done over $Δt≈1000$ in the saturated state.⁶⁶ In this figure, the gray vertical line marks the location of the Zeman wavenumber $kΩ$ defined in the presence of the energy injection due to the precessional instability as $kΩ=(Ω3/⟨A2D+A3D⟩)1/2$ (equal to $1/⟨A2D+A3D⟩1/2$ in non-dimensional units), where $⟨A2D+A3D⟩$ is the volume-averaged total injection term introduced. This definition of $kΩ$ differs from the usual one used in the rotating turbulence theory^31,34,67 in that the energy injection rate, ε, due to an external forcing is replaced here by the injection due to the instability, but these definitions are consistent in steady state.

The most remarkable aspect is the different shape and scaling of energy spectra for the 2D and 3D modes. The 2D mode energy dominates over the 3D one at small wavenumbers $k≲10$ where it increases with decreasing k, reaching a maximum at the largest box scale, with its horizontal component being about an order of magnitude larger than the vertical one. This corresponds to large-scale horizontal vortical motions in physical space, as is seen in Figs. 4 and 7. At higher $10≲k≲kΩ$⁠, the horizontal and vertical components are comparable in the $E2D$ spectrum and its slope is close to $k−3$⁠, exhibiting the same power-law dependence of rotating geostrophic 2D turbulence,^{20,29,35,36,68,69} which does not appear to change with precession parameter Po.

The energy spectrum of 3D waves, $E3D$⁠, has a peak at larger $k≈15$ than that of $E2D$ (which approximately coincides with the peak of injection $A3D$ in Fig. 12). $E3D$ decreases then at lower wavenumbers, while at higher wavenumbers, $15≲k≲kΩ$ follows a scaling $∼k−2±0.5$⁠, which has been typically observed in forced rotating turbulence of inertial waves.^{30,31,34,36,59} However, in contrast to these papers using a forcing in a very narrow wavenumber band, we do not prescribe the forced wavenumber a priori; instead, the injection wavenumbers are determined by the background flow itself through the precessional instability and extend over a broad range (see below). Precession influences the scaling exponent of the $E3D$ spectrum: its slope seems to become shallower with increasing Po, as is seen in Figs. 9 (compare left and right panels) and 14 below. Like for 2D mode energy, also for 3D mode energy, horizontal and vertical components are comparable at higher wavenumbers, but the horizontal one dominates at lower wavenumbers.

It is seen in Fig. 9 that the observed power-law scalings of both 2D and 3D mode energy spectra occur at $k<kΩ$ and therefore are strongly influenced by rotation and precession, deviating from the classical Kolmogorov $k−5/3$ spectrum. However, as it is seen in this figure, with increasing precession strength Po, the Zeman wavenumber $kΩ$ (gray vertical lines) decreases; that is, the effect of rotation becomes increasingly weaker for smaller and smaller k. We will see below that in this case the energy spectrum at $k>kΩ$ indeed approaches the Kolmogorov spectrum.

To see the structure of spectra of energy and dynamical terms in the quasi-steady precessional turbulence, in Figs. 10 and 11 we show the time-averaged spectra of the kinetic energy E, energy injection A, and the total nonlinear transfer NL in two different horizontal $(kx,ky)−$planes: at k_z = 0 for the 2D modes and at $kz=2π/Lz$ for the 3D wave modes. We have chosen $kz=2π/Lz$ because it corresponds to the maximum injection along k_z axis, for which therefore the precession instability reaches the largest growth rate in the box.⁵⁰ The most striking observation is the anisotropic nature of the 2D manifold in Fourier space (top row), for $A2D$ and $NL2D$⁠, and hence, for the kinetic energy spectrum, $E2D$ determined by the joint action of these terms, which are all localized at smaller wavenumbers, with a clear inclination toward the k_x axes. The injection term $A2D$ is always positive, implying some energy injection from the basic flow into vortices; however, $NL2D$ is negative at the same wavenumbers for Po = 0.2, but changes sign at Po = 0.3. As a result, the dynamical balances for 2D modes are different for these two values of Po, which will be discussed below. In contrast, the 3D manifold exhibits a quasi-isotropic distribution (bottom row, similarly at larger $kz>2π/Lz$ not shown here) for both Po = 0.2 and 0.3, whose range extends over larger wavenumbers than that of 2D quantities. Notice the impact of Po on the 3D manifold: for Po = 0.2 a weak preferential direction (anisotropy) is present in the center of the wave plane, which tends to become more isotropic with increasing Po. Comparing the $A3D$ and $NL3D$⁠, the first injection term, which is due to the precessional instability, is always positive and appreciable at $5<|kx|,|ky|<25$ (yellow/red area), while the second nonlinear term is negative (blue) and also appreciable at these wavenumbers. The similar shape of these two functions in Fourier space and comparable absolute values imply that these two processes are in balance: 3D modes receive energy from the precessional background flow predominantly in the range $5<|kx|,|ky|<25$⁠, while nonlinearity, counteracting injection at these wavenumbers, transfers this energy to other 3D and 2D modes with different wavenumbers.

From Figs. 10 and 11 showing the distribution of total nonlinear terms $NL2D$ and $NL3D$ in Fourier space, one cannot establish specifically what kind of transfer mechanisms operate, that is, whether the cascades inside a given manifold are direct or inverse or if there are transfers of energy between these two manifolds, since these terms encapsulate nonlinear interactions among all kinds of modes. To get insight into the details of linear (energy injection) and nonlinear cascade processes in the precessional turbulence, in Fig. 12, we show the shell- and time-averaged spectra of all the dynamical terms—injection A, nonlinear NL, and dissipation D terms—entering Eqs. (13) and (14) again for Po = 0.2 and Po = 0.3, as we did for the bursty case Po = 0.075 in the above section. The precession instability injects energy into 3D waves, which is described by positive $A3D$ (blue, right panels). Unlike the case of a forcing localized about certain wavenumber,³⁵ the injection due to the instability extends over a range of wavenumbers, reaching a peak at $k=6π$⁠, and increases with increasing Po. Since it is a quasi-steady state, the energy injection is balanced by nonlinear transfers, $A3D+NL3D3D+NL2D3D≈0$⁠, at those dynamically active wavenumbers where $A3D>0$ is appreciable (the role of viscous dissipation $D3D$ (black, right panels) is not important at these wavenumbers). Therefore, being negative at those injection wavenumbers, $NL3D3D<0$ (green, left panels) and $NL2D3D<0$ (red, left panels) drain energy from the active 3D modes there and transfer it, respectively, to smaller-scale 3D waves due to positive $NL3D3D>0$ at k > 40 (forward/direct cascade) and to 2D vortical modes. The latter process is mediated by positive $NL3D2D>0$ (red, left panels) at k > 18, peaking at the same $k=6π$⁠. At these wavenumbers, 2D–2D transfer term is negative $NL2D2D<0$ (green, left panels), with a minimum also at $k=6π$⁠, and causes inverse cascade of 2D mode energy to even smaller wavenumbers, where it is positive $NL2D2D>0$ and reaches a maximum at $k=2π$ that corresponds to the largest box scale. These vortices draw some energy from the basic flow as well due to $A2D$ term, which has a peak at the same $k=2π$ as $NL2D2D>0$⁠. It is seen that $A2D$ decreases, whereas $NL2D2D$ increases with Po. Notice that in the basic shear flow $Ub$⁠, idealized steady 2D vortices would be linearly stable against precessional instability (i.e., $A2D=0$⁠),⁵² so the positive $A2D>0$ in this case can be attributed to the fact that these vortices are being nonlinearly generated by the waves. The 2D–3D nonlinear interaction term $NL2D3D$ is positive at large wavenumbers k > 40, redistributing part of energy of 2D large-scale vortices back to smaller-scale 3D waves (forward cascade). Thus, as it is shown from Fig. 12 (left panels), similar to that for 3D modes, also for 2D modes, there is a balance among production of these modes by 2D–3D transfers, energy extraction, 2D–2D transfers, and viscous dissipation, $A2D+NL2D2D+NL3D2D+D2D≈0$⁠. Note that viscous dissipation for 2D and 3D modes has completely different behaviors (compare black curves in left and right panels). It is stronger and more significant for the 2D vortices: $D2D$ has a clear minimum at low wavenumber $k=2π$⁠, which coincides with the peak of $A2D$ and $NL3D2D$⁠, and counteracts these terms, indicating a dissipative nature of the vortices. On the other hand, for 3D modes viscosity is important only at higher k > 40; that is, small scales are dissipative, in contrast to that in the bursty regime, where viscous and injection scales coincide (right panels in Fig. 8).

Note that the strength of all the dynamical processes depicted in Fig. 12 increases with increasing Po. In all cases, the peaks of energy injections into 3D wave modes due to the precessional instability are concentrated at wavenumbers smaller than the corresponding Zeman wavenumber, $kin,3D=6π<kΩ$ (gray vertical lines in Fig. 9). As we have seen in the left panels of Fig. 12, this injection affects the 2D–3D nonlinear transfer $NL3D2D$⁠, which describes driving of 2D vortices by 3D wave modes, and also has a maximum at the same k_in. This is in agreement with the general condition for the upscale/inverse energy cascade of 2D vortices toward wavenumbers smaller than the injection one, that is, $k<kin$⁠, in rotating turbulence.³⁵ However, in the present case of precessional driving, in the inverse cascade regime, the energy spectrum of 2D modes is slightly steeper than the $k−3$ slope (Fig. 9), which is usually observed in the same regime in a purely rotating case.

Overall, the above-described processes of nonlinear transfers inside the 3D manifold, inside the 2D manifold, and the coupling between these two manifolds are consistent with previous spectral analysis of turbulence dynamics under rotation and an imposed external forcing.^31,35,59 In particular, in Fig. 12 we observe the split (simultaneous inverse and forward) cascade of energy typical of rotating turbulence as demonstrated in those papers, that is, forward cascade of 3D wave mode energy to high wavenumbers (small scales) due to $NL3D3D$ and $NL2D3D$⁠, and inverse cascade of 2D modes to small wavenumbers (large scales) due to $NL2D2D$⁠.

To confirm the overall type (direction) of the nonlinear cascades inferred above on the basis of the nonlinear transfers NL as a function of k, we also analyze shell-to-shell flux of the total energy defined as³¹ 

Figure 13 shows the resulting flux function F(k) and, for reference, the wavenumbers $kin,2D=2π$ and $kin,3D=6π$ at which the injection terms $A2D$ and $A3D$⁠, respectively, reach their maximum (Fig. 12). The gray line in this figure represents the Zeman wavenumber $kΩ$⁠. The shape of the total fluxes is qualitatively similar for Po = 0.2 and 0.3 and indeed display split/dual cascade: they are positive, F > 0, at $k>kin,2D$ with a maximum value around $k≈50$⁠, indicating a forward cascade of energy, and negative, F < 0, at small wavenumbers $k<kin,2D$⁠, indicating inverse cascade. These forward/inverse cascade regimes deduced from the behavior of the energy flux function F(k) in fact confirm those found above based on the behavior of the transfer functions in Fig. 12. Specifically, the forward cascade at $k>kin,2D$ is related to the transfer of 3D wave mode energy to higher k, while the inverse cascade at $k<kin,2D$ is related to the transfer of 2D vortical mode energy to smaller k.

In the following, we draw conclusions on the properties of 3D wave modes, which are the ones directly influenced and driven by the precession instability. We have already seen clear indications that these modes exhibit characteristics of isotropicity, direct cascade, and decreasing the wavenumber range where the rotation is substantially dominant, that is, decreasing the Zeman wavenumber $kΩ$⁠, with increasing precession intensity. In order to confirm and generalize these concepts, we run another simulation for quite a high precession parameter Po = 0.5 and with the same $Re=104.5$ to check this trend.

Figure 14 shows the spectra of 3D mode energy for the three precession parameters. This time we compensated $E3D$ spectra with the Kolmogorov spectrum $k−5/3$ to better see whether the energy spectrum approaches the Kolmogorov one. Indeed, it is seen in this figure that increasing Po, the Zeman wavenumber decreases and the compensated spectrum at $k>kΩ$ becomes gradually flatter, indicating approach to the Kolmogorov one $k−5/3$ already at Po = 0.5, that is, the regime of isotropic homogeneous turbulence. Thus, we showed that the rotating-dominated range of wavenumbers $k<kΩ$ is narrowed as Po increases because of dramatic decrease in Zeman wavenumbers (e.g., $kΩ=103$ for Po = 0.2 reducing to $kΩ=38$ for Po = 0.5). This is also reflected in the increase in the Rossby number at the injection scales $Ro=(A2Dkin,2D2+A3Dkin,3D2)1/3/Ω$⁠. Specifically, we have $Ro≈0.176$ for Po = 0.2, $Ro≈0.24$ for Po = 0.3 and Ro = 0.32 for Po = 0.5.