2.1. The wave–vortex decomposition and the potential vorticity field

The solutions of the linear equations obtained by setting $\epsilon =0$ in (2.3) consist of inertia-gravity waves and a divergence-free vortical mode that is in geostrophic balance. The equations governing these two fields, denoted with subscripts ‘W’ and ‘V’ for waves and vortical mode, are

(2.4a)$$\begin{gather} \frac{\partial{ {{ \boldsymbol{v}}_{\scriptscriptstyle{W}}} }}{\partial{t}} + \boldsymbol{f} \times {{ \boldsymbol{v}}_{\scriptscriptstyle{W}}} + \boldsymbol{\nabla} {{h}_{\scriptscriptstyle{W}}} = 0 , \end{gather}$$

(2.4b)$$\begin{gather}\frac{\partial{{{h}_{\scriptscriptstyle{W}}}}}{\partial{t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} {{ \boldsymbol{v}}_{\scriptscriptstyle{W}}}=0, \end{gather}$$

(2.5a)$$\begin{gather} \boldsymbol{f} \times {{ \boldsymbol{v}}_{\scriptscriptstyle{V}}} + \boldsymbol{\nabla} {{h}_{\scriptscriptstyle{V}}} = 0, \end{gather}$$

(2.5b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} {{ \boldsymbol{v}}_{\scriptscriptstyle{V}}} =0. \end{gather}$$

Inertia-gravity waves in (2.4) possesses the dispersion relationship $\omega =\sqrt {f^2 + k^2}$, where $\boldsymbol {k}=(k_x, k_y)$ is the wavenumber vector and $k = \vert \boldsymbol {k} \vert$ is the wavenumber vector magnitude. A special case of the above decomposition to be noted is the non-rotating limit, $f=0$. In this limit, from (2.5a) we see that ${{h}_{\scriptscriptstyle {V}}} = 0$, implying that the vortical flow has no height field fluctuations with respect to the constant mean height.

Subtracting $f$ times (2.4b) from the curl of (2.4a) gives us

(2.6a,b) \begin{equation} \frac{\partial{}}{\partial{t}} \left(\zeta_{\scriptscriptstyle{W}} - fh_{\scriptscriptstyle{W}} \right)= 0, \quad \zeta = \boldsymbol{\hat{z}} \boldsymbol{\cdot} \boldsymbol{\nabla} \times \boldsymbol{v}. \end{equation}

The linear potential vorticity, $\zeta - f h = {{ \zeta }_{\scriptscriptstyle {V}}} - f {{h}_{\scriptscriptstyle {V}}}$, is therefore completely associated with the vortical flow and waves carry no potential vorticity, i.e. ${{ \zeta }_{\scriptscriptstyle {W}}} - f {{h}_{\scriptscriptstyle {W}}} = 0$. Notice that in the non-rotating limit, $f=0$, waves carry no vorticity, $\zeta$, and all the vorticity is associated with the vortical mode.

The linear wave-balance decomposition given in (2.4) and (2.5) has been used in a wide variety of previous works. For example, Warn (Reference Warn1986) used the decomposition to study statistical equilibrium of rotating shallow water, Majda & Embid (Reference Majda and Embid1998) used it to rigorously derive the quasi-geostrophic equation, Remmel & Smith (Reference Remmel and Smith2009) used it to derive intermediate wave-balance models, Ward & Dewar (Reference Ward and Dewar2010) used it to examine scattering of waves by the vortical mode and Thomas & Yamada (Reference Thomas and Yamada2019) used it to examine energy exchanges between waves and balanced flows. We will use the same decomposition to separate inertia-gravity waves from balanced flows in this work.

2.2. Energy expressions

The RSWE model (2.3) conserves energy, with the energy equation being

(2.7)\begin{equation} \frac{{\rm d}{E^{\epsilon}}}{{\rm d}{t}} = \frac{{\rm d}{}}{{\rm d}{t}} \int_{\boldsymbol{x}} \left(\frac{1}{2} {\boldsymbol{v}}^2 + \frac{1}{2} c^2 h^2 + \epsilon \frac{1}{2} h \boldsymbol{v}^2 \right) {\rm d} \boldsymbol{x} =0. \end{equation}

Total energy, $E^{\epsilon }$, is an integral conserved quantity of (2.3) and can be expressed in physical space and spectral space as

(2.8)\begin{align} E^{\epsilon} &= \int_{\boldsymbol{x}} \left \{ \frac{1}{2} (1 + \epsilon h) \left( u^2 + v^2 \right) + \frac{1}{2} h^2 \right \} {\rm d} \boldsymbol{x} \nonumber\\ &=\sum_{\boldsymbol{k}} \left \{ \left \vert \widehat{\left\{ u \sqrt{ 1 + \epsilon h } \right\} } (\boldsymbol{k}, t) \right\vert^2 + \left \vert \widehat{ \left\{ v \sqrt{ 1 + \epsilon h } \right\} } (\boldsymbol{k}, t) \right\vert^2 + \left\vert \hat{h} (\boldsymbol{k}, t) \right\vert^2 \right \} = \sum_{k} \hat{E}^{\epsilon} (k, t), \end{align}

where $\hat {g} (\boldsymbol {k})$ denotes the Fourier transform of $g(\boldsymbol {x})$ and the last summation above is over $k = \vert \boldsymbol {k} \vert$ which is obtained by integrating over angles in wavenumber space in polar coordinates. For $\epsilon \ll 1$, the cubic energy expression (2.8) can be approximated by the following quadratic energy expression:

(2.9)\begin{align} E &= \int_{ \boldsymbol{x}} \left \{ \frac{1}{2} \left( u^2 + v^2 \right) + \frac{1}{2} h^2 \right \} {\rm d} \boldsymbol{x} = \sum_{\boldsymbol{k}} \left \{ \left \vert \hat{u} (\boldsymbol{k}, t) \right \vert^2 + \left\vert \hat{v} (\boldsymbol{k}, t) \right\vert^2 + \left \vert \hat{h} (\boldsymbol{k}, t) \right\vert^2 \right \} \nonumber\\ &= \sum_{k} \hat{E} (k, t). \end{align}

Notice that the quadratic energy (2.9) is the sum of quadratic square terms while the cubic energy (2.8) contains cubic nonlinear terms in addition to quadratic terms. However, the cubic terms in (2.8) are $O(\epsilon )$ and therefore only a small perturbation to the leading-order quadratic terms. In all our numerical integrations we found that the quadratic and cubic energies were close to each other: see tables 1 and 2 for specific energy values for different flows.

Table 1. Parameters for flows where the vortical mode was not initialized.

Table 2. Parameters for flows where the vortical mode was initialized.

The linear wave-balance decomposition based on (2.4) and (2.5) provides us with a vortical field that is in geostrophic balance for $f \neq 0$ while linear inertia-gravity waves form the flow field that is orthogonal to the balanced field such that quadratic energy (2.9) is the sum of these two quantities i.e. $E={{E}_{\scriptscriptstyle {V}}} + {{E}_{\scriptscriptstyle {W}}}$, where $E$, ${{E}_{\scriptscriptstyle {V}}}$ and ${{E}_{\scriptscriptstyle {W}}}$ denote total, balanced and wave energies, respectively. The cubic energy expression $E^{\epsilon }$ on the other hand does not decompose into the exact sum of wave and vortical energy. In this work we will primarily use the linear wave-balance decomposition and the quadratic energy expression to quantify the two different flow fields.

2.3. Numerical integration details

We numerically integrated (2.3) with $\epsilon =0.1$ using a dealiased pseudospectral scheme in the periodic domain $(x, y) \in [0, 2 {\rm \pi}]^2$. The time integrations were performed by using the fourth-order Runge-Kutta method, and numerical convergence in time was ensured by gradually reducing the time step until the results did not depend on the time step. Following this procedure, a time step of $6.25 \times 10^{-4}$ was seen to be required for the time integrations. We used the two-thirds dealiasing and denote the maximum wavenumber obtained after dealiasing by $k_{max}$.

We forced waves at large scales and maintained the energy of the wave field in the wavenumber band $k \leqslant 5$, such that the forced wavenumbers had the same energy at all times. This particular forcing scheme, whose implementation details are explained in Appendix A, has the advantage that it does not force energy into the system at a prespecified rate. If the flow undergoes a forward energy cascade resulting in the transfer of energy to smaller scales, the forcing would proportionally inject energy into the low wavenumber band to maintain the large-scale energy lost to small scales. If there is no transfer of energy to smaller scales, the forcing would shut off and not inject energy into the system. Such a forcing scheme is beneficial since it injects energy into large scales based on the energy transfer rate from large scales to small scales and not based on a preset injection rate. We specifically used this forcing scheme in the present study since we do not a priori know the rate of downscale energy transfer at different rotation rates, and therefore did not want to externally enforce an energy injection rate.

To analyse the turbulent cascade of waves across an inertial range, it would be ideal to have numerical solutions of the inviscid equations. Of course, this is quite impractical since the energy cascading to smaller scales and eventually reaching grid scale has to be removed by some sort of dissipation mechanism. In this regard, a pragmatic solution is to have as much long inviscid inertial range as possible at affordable spatial resolutions, thereby restricting dissipation to close to the grid scale. Adding the usual Laplacian dissipation operator was seen to be ineffective for this purpose, since it considerably limited the inviscid inertial range. We therefore added bi-Laplacian hyperdissipative terms of the form $-\nu \nabla ^4 (u, v, h)$, with hyperviscosity $\nu =5.54 \times 10^{-10}$, to the right-hand sides of (2.3). A comparison between physical and spectral space features of solutions corresponding to these two dissipation operators are given in Appendix B, indicating that the inertial range is severely restricted when using the Laplacian dissipation term. We therefore used the bi-Laplacian hyperdissipative term for the present study.

It is also noteworthy that having dissipation act on the $\boldsymbol {v}$ equation in (2.3a), similar to that implemented by the previous works of Augier et al. (Reference Augier, Mohanan and Lindborg2019), Yuan & Hamilton (Reference Yuan and Hamilton1994), Polvani et al. (Reference Polvani, McWilliams, Spall and Ford1994) and Farge & Sadourny (Reference Farge and Sadourny1989) for example, results in a positive definite expression for the dissipation of the quadratic energy $E$ given in (2.9), while an $O(\epsilon )$ extra term appears in the dissipation equation of the cubic energy expression $E^\epsilon$, given in (2.8). Despite the presence of the $O(\epsilon )$ perturbation term, we found these two different dissipation rates to be quite close in the $\epsilon \ll 1$ regime we are set in: see dissipation values given in tables 1 and 2. The dissipation term expressions and their comparisons are given in Appendix C. Augier et al. suggests an alternate strategy where dissipation is added to the $h \boldsymbol {v}$ equation instead of the $\boldsymbol {v}$ equation, since this would lead to a positive definite cubic energy dissipation expression. However, it is in principle difficult to associate a physical justification for either of these different strategies of adding dissipation to the equations, since hyperdissipation is essentially an artificially introduced mechanism to remove energy reaching grid scale and thereby extending the inviscid inertial range. We therefore persisted with the commonly used strategy of adding dissipation to the $\boldsymbol {v}$ equation and not the $h \boldsymbol {v}$ equation, as in the above-mentioned previous studies.

Based on multiple numerical integrations we confirmed that the resolution $k_{max} = 1024$, which corresponds to 3072 grid points in the $x$- and $y$-directions, generated solutions that were insensitive to further increase in spectral resolution. Specifically, this resolution ensured that the details of the turbulent transfer across the inertial range and the associated statistics described in the following sections were robust and did not change on further increase in resolution. All the results presented in this paper were therefore generated with the spectral resolution $k_{max} = 1024$, which corresponds to $3072^2$ grid points.

In the numerical solutions, inertia-gravity waves steepened and formed strong gradients. We will be using the term ‘shock’ to describe these strong gradients, following the terminology used by the previous studies mentioned above. These shocks have finite thickness due to the presence of viscosity, and are therefore not sharp jump discontinuities that would form the solutions of the completely inviscid equations.

3.1. Dominance of linear dynamics and weakly nonlinear interactions

Since the turbulent dynamics we investigate are set in the weakly nonlinear regime, we expect the flow to be dominated by the linear modes. Here we will compare the weak nonlinearity of the flow, starting with the energy distribution across scales. As discussed earlier with specific comparisons given in tables 1 and 2, although the exact energy expression (2.8) is cubic, in the weakly nonlinear regime with $\epsilon \ll 1$, we found the quadratic and cubic energies to be indistinguishable. A more detailed example comparison across wavenumber space is shown in figure 4(a) that plots the energy spectrum of the flow computed using the quadratic and cubic expressions for the $f=1$ case. Notice how the two spectra overlap, indicating that the quadratic expression approximates the cubic expression quite well.

Figure 4. The figure shows (a) $\hat {E}^{\epsilon } (k)$ and $\hat {E} (k)$ and (b) $k \overline {\hat {R}^\epsilon } (k)$ and $k \bar {\hat {R}} (k)$ for $f=1$.

We will next look at the transfer function distribution across scales. Recall (3.1) that indicates that the energy $\hat {E} (k, t)$ contained at wavenumber $k$ is transferred to other wavenumbers with the transfer being captured by the transfer term $\hat {R} (k, t)$. Along similar lines as those used to arrive at (3.1), we could develop an equation for the time evolution of the cubic energy expression, $\hat {E}^{\epsilon } (k, t)$, and let us identify the corresponding energy transfer term by $\hat {R}^\epsilon (k, t)$. Figure 4(b) shows a comparison between the time-averaged transfer functions $\overline {\hat {R}^\epsilon } (k)$ and $\bar {\hat {R}} (k)$: notice how the two transfer terms overlap in the inertial range. Some difference is seen in the dissipation range, although this difference is exaggerated in the figure since we are comparing $k \bar {\hat {R}} (k)$ and $k \overline {\hat {R}^\epsilon } (k)$ and in the dissipation range as $k$ approaches $10^3$. We therefore conclude that the transfer function based on the quadratic energy expression is a good approximation to the exact energy transfer function based on the cubic energy expression, particularly so in the inertial range. Along similar lines, table 1 lists the domain-integrated quadratic and cubic energy, along with the energy dissipation associated with the quadratic and cubic energy. The comparison reveals that the quadratic quantities closely approximate the cubic variables.

We next examine the linearity of the wave field. Recall that the dispersion relationship of inertia-gravity waves in (2.4) is $\omega =\sqrt {f^2 + k^2}$. For all the numerical integrations we performed, we saved the time series of flow variables in wavenumber space and performed a Fourier transform in time, providing us with the flow variables in the frequency–wavenumber space. Figure 5 shows a representative example illustration: the $x$-velocity component, u, in frequency–wavenumber space for $f=1$ and $f=3$. The white line shows the linear waves’ dispersion relationship. We find that the low wavenumber region shows some spread towards higher frequencies, primarily due to the time-variable external forcing acting at low wavenumbers. Nevertheless, observe that the velocity component distribution aligns well with the dispersion relationship in the inertial range, away from the forcing region. It is also noteworthy that the velocity distribution in figure 5(b) decays much more rapidly than in figure 5(a) with increasing $\omega$ and $k$, concomitant with the results shown in figure 3 that reveal weakened waves’ forward energy flux and reduction in wave energy at small scales with increasing rotation rate.

Figure 5. The figure shows the $\omega$–$k$ space of x-velocity component, $u$, for (a) $f=1$ and (b) $f=3$. The dispersion relationship of the waves, $\omega =\sqrt {f^2 + k^2}$, is marked by the white line.

3.2. Scale-locality of energy transfers

Having confirmed the dominance of linear modes and weakly nonlinear interactions, we return to analysing waves’ energy transfer across scales. Our previous discussion revealed that waves’ energy undergo a forward flux. However, it is not clear if the waves’ forward energy flux is composed of a scale-local cascade, as is the case in 3-D homogenous isotropic turbulence (HIT). To explore the locality of turbulent transfers, we ignore the forcing and dissipative terms and write the energy transfer equation in spectral space as

(3.3)\begin{equation} \frac{\partial{\hat{E} }}{\partial{t}} = \sum_{p} \sum_{q} \hat{S} (k \vert p, q, t) = \sum_{p} \hat{T} (k \vert p, t) = \hat{R} (k, t). \end{equation}

The above equation is similar to (3.1), but with a few intermediate steps introduced. Here $\hat {S} (k \vert p, q, t)$ is the triadic nonlinear term in spectral space responsible for all energy transfer between wavenumbers $\boldsymbol {p}$, $\boldsymbol {q}$ and $\boldsymbol {k}$ subject to the triadic constraint $\boldsymbol {k}=\boldsymbol {p}+ \boldsymbol {q}$. Summing $\hat {S} (k \vert p, q, t)$ over $q$ gives us the transfer function $\hat {T} (k \vert p, t)$, which captures transfers between wavenumbers $k$ and $p$. Finally, summing over $p$ gives us the net energy transfer at a specific wavenumber, $\hat {R}$, which appears in (3.1).

The locality of turbulent transfers can be gauged by examining the behaviour of the $\hat {T}(k \vert p, t)$ term across the inertial range. A transfer is identified as local if most of the transfer term $\hat {T} (k \vert p, t)$ is concentrated in the neighbourhood of wavenumbers $k \sim p$. This means that the transfer at a specific wavenumber $k$ is due to wavenumbers $p$ in the neighbourhood of $k$ such that $k/p \sim 1$. On the other hand, if wavenumbers $p$ significantly away from $k$ such that $k/p \gg 1$ or $k/p \ll 1$ contribute towards the transfer at wavenumber $k$, it is termed a non-local transfer (Brasseur & Wei Reference Brasseur and Wei1994; Lesieur Reference Lesieur2008). We computed $\hat {T}(k \vert p, t)$ across the inertial range and found it to be highly fluctuating in time. We therefore time-averaged $\hat {T}(k \vert p, t)$ over a long enough time window to obtain $\bar {\hat {T}} (k \vert p)$, which has temporal fluctuations removed. A series of examples of $\bar {\hat {T}} (k \vert p)$ is shown in figure 6 for a specific wavenumber in the middle of the inertial range for different rotation rates. We obtained a similar result for other wavenumbers in the inertial range. Observe from the figure that $\bar {\hat {T}}$ is highly localized, with negative values for $k/p < 1$ and positive values for $k/p > 1$. The localized transition of $\bar {\hat {T}}$ from negative to positive values is a reflection of the localized forward energy cascade transferring energy from large scales to small scales. It is interesting that Augier et al. (Reference Augier, Mohanan and Lindborg2019) speculated the possibility of non-local energy transfer from large to small scales by shallow water gravity waves (see discussions right below their figure 3), while our analysis indicates that the downscale energy transfer of waves is scale-local.

Figure 6. The figure shows $\bar {\hat {T}} (k \vert p)$ for four different rotation rates for $k=200$, in the middle of the inertial range. The different panels are (a) $f=0$; (b) $f=1$; (c) $f=2$; (d) $f=3$. Observe that the magnitude of $\bar {\hat {T}} (k \vert p)$ reduces with increasing $f$, as a consequence of weakening of the transfers at higher rotating rates.

Although the above analysis reveals that the turbulent forward cascade of waves in RSWE is scale-local, non-local modes can be involved in the transfer. To see this, recall that in (3.3) $\hat {S}(k \vert p, q, t)$ was summed over all $q$ wavenumbers to obtain $\hat {T} (k \vert p, t)$. Although the lengths of wavenumbers $k$ and $p$ can be comparable, $q$ can be of really small length resulting in $q \ll k \sim p$. Consequently, although the turbulent transfer is scale-local, there is nothing preventing widely separated or non-local wavenumbers affecting the turbulent transfer. To examine the possibility of non-local modes being involved in the transfer, we summed $\hat {S} (k \vert p, q, t)$ over selected wavenumbers $q$ as follows.

For a triad $(k, p, q)$ contributing to the turbulent transfer, we define $\alpha = \text {max}(k, p, q)/\text {min}(k, p, q)$ as the ratio of the maximum to minimum wavenumber. In hydrodynamic turbulence a triad is considered to be local if $\alpha \leqslant 2$ and non-local if $\alpha > 2$ (Lesieur Reference Lesieur2008). We used $\alpha$ to compute the partial transfer function $\hat {T}(k \vert p; \alpha \leqslant m, t )$ as

(3.4)\begin{equation} \hat{T}(k \vert p; \alpha \leqslant m , t) = \sum_{q \ \text{such that } \alpha \leqslant m } \hat{S} (k \vert p, q, t), \end{equation}

where $m$ is a number we set so that the summation on the right-hand side above can be computed subject to the constraint that only triads that fall within a certain value of length ratios are considered. For example, setting $m=2$ would give us the transfer function considering only local triads, $\hat {T}(k \vert p; \alpha \leqslant 2, t)$.

We computed $\bar {\hat {T}} (k \vert p; \alpha \leqslant m)$, the time-averaged $\hat {T}(k \vert p; \alpha \leqslant m, t )$, for different $m$ values across different rotation rates and a plot that summarizes the main results of our analysis is shown in figure 7 for $f=1$. Notice that when we set $\alpha \leqslant 2$ (red curve in figure 7), which constrains the triads to local modes alone, we do not capture any part of the complete transfer function (black curve in figure 7). Setting $\alpha \leqslant 17$, which allows for non-local triads, results in the green curve in figure 7: notice that the green curve captures approximately 50 % of the maxima and minima of the total transfer function shown by the black curve. Finally, further increasing the non-locality of the modes and setting $\alpha \leqslant 117$ leads to the blue curve in figure 7, which captures most of the total transfer function shown by the black curve. This analysis reveals that non-local modes are crucial for obtaining the complete transfer function. Although $\bar {\hat {T}}$ itself shows highly localized behaviour, as seen in figure 6, non-local triads form an important part of $\bar {\hat {T}}$. We therefore conclude that the turbulent downscale transfer of waves in RSWE is a result of the interaction of non-local triads, although the transfer itself is scale-local. The same conclusion, i.e. local transfer due to non-local triads is a well established result in 3-D HIT (Domaradzki & Rogallo Reference Domaradzki and Rogallo1990; Yeung & Brasseur Reference Yeung and Brasseur1991; Ohkitani & Kida Reference Ohkitani and Kida1992; Waleffe Reference Waleffe1992; Zhou Reference Zhou1993; Brasseur & Wei Reference Brasseur and Wei1994; Domaradzki & Carati Reference Domaradzki and Carati2007; Eyink & Aluie Reference Eyink and Aluie2009; Cardesa, Vela-Martin & Jimenez Reference Cardesa, Vela-Martin and Jimenez2017). Our analysis reveal that the same feature is exhibited by the turbulent forward energy cascading inertia-gravity waves in RSWE.

Figure 7. The figure shows $\bar {\hat {T}}(k \vert p)$ and its partial representations for $f = 1$. The black curve is the same curve that appears in figure 6(b), although the $x$-axis range is reduced here compared with figure 6(b) for better visualization of the different curves.

3.3. Finite non-zero time delay in the cascade

As discussed before, the locality results seen in figure 6 required time averaging to remove fast fluctuations. An instantaneous snapshot of $\hat {T}(k \vert p, t)$ did not give clear insights on locality of transfers, although the time-averaged or the slow part of $\hat {T}$ did. This is an indication that the slow energy transfers happen by scale-local transfers, although fast fluctuations associated with energy transfers do not necessarily adhere to locality constraints. Fast fluctuations in energy, transfer functions and fluxes were seen to be a generic feature of the turbulent transfers we examined. An example is shown in figure 8(a) where the time series of transfer function $\hat {R} (k, t)$ for different wavenumbers is plotted in an arbitrary time interval for $f=1$. Notice that the amplitude of oscillations is quite high for low wavenumbers, e.g. $k=20$, while they reduce with increasing wavenumbers and is relatively negligible at high wavenumbers, e.g. $k=200$.

Figure 8. Panel (a) shows the variation of $\hat {R} (k, t)$ with time for different $k$ in an arbitrary time interval. Panels (b,d) show the correlation between $\hat {R} (k_1, t)$ and $\hat {R} (k_2, t + \tau )$, while (c,e) show the correlation between the slow part of the transfer functions, $\hat {R}_{s} (k_1, t)$ and $\hat {R}_{s} (k_2, t + \tau )$, respectively. Panel ( f) shows the variation of the maximum correlation curve, i.e. the white lines in (c,e) above, across the inertial range. Here (a)  $f=1$; (b) $f=1$; (c) $f=1$; (d) $f=3$; (e) $f=3$.

Due to scale-local transfers, energy from large scales will not be able to reach small dissipative scales instantaneously by means of a direct transfer. Instead, a certain finite time is to be expected for energy to get transferred from large to small scales. One way to confirm this non-zero time delay is to compute the time-lag correlation between time series’ of domain-integrated energy and dissipation, as is often done in 3-D HIT investigations (see e.g. Cardesa et al. Reference Cardesa, Vela-Martin, Dong and Jimenez2015). Energy and dissipation may be considered as large-scale and small-scale variables, respectively, which means that their time correlation will inform us about the presence of a delay between the two signals. We computed the time-lag correlation between energy and dissipation for different rotation rates and for all cases found that the maximum correlation between the time series occur with a time lag, hinting at the delay in energy reaching dissipative scales from large domain scales. Here we skip a detailed discussion of these figures, the figures being similar to those seen in 3-D HIT studies, see for example figure 2 in Cardesa et al. (Reference Cardesa, Vela-Martin, Dong and Jimenez2015). For a more illuminating illustration of the time delay in the cascade process, we will now discuss results of an investigation we undertook inspired by the study of Khurshid, Donzis & Sreenivasan (Reference Khurshid, Donzis and Sreenivasan2021) for 3-D HIT, where the time delay between the transfer functions at different wavenumbers is examined.

We computed the time-lag correlation of $\hat {R}$ across different wavenumbers, i.e. the correlation between $\hat {R} (k_1, t)$ and $\hat {R} (k_2, t + \tau )$, where $k_1$ was chosen to be 20, at large scales and a little away from the forced wavenumbers, and we increased $k_2$ starting from $k_1$ and went across the inertial range up to $k_2=600$, where dissipation starts to become non-negligible. Figures 8(b) and 8(d) show the correlation as a function of $k_2$ and the time lag $\tau$ for $f=1$ and $f=3$, respectively. Scanning through these figures reveals that not much correlation is seen across the inertial range. Following this, we performed a running time average on the transfer function $\hat {R}$ at different wavenumbers to remove the high frequency oscillatory part of the transfer function. The averaging window was chosen by trial and error and approximately 10 domain-scale wave periods was seen to be a sufficient time window for removing fast fluctuations and obtain $\hat {R}_{s}$, the slow component of the transfer function. We then computed the correlation between slow part of the transfer functions at different wavenumbers, $\hat {R}_{s} (k_1, t)$ and $\hat {R}_{s} (k_2, t + \tau )$, and the results for $f=1$ and $f=3$ are shown in figures 8(c) and 8(e). Notice that there is a clear correlation between the slow transfer functions across the inertial range, while there was no noticeable correlation between the full transfer functions. In figures 8(c) and 8(e), the maximum correlation curve is marked by a white line. These maximum correlation curves are separately replotted in figure 8( f). Notice that the time lag $\tau$ increases from $k_1=20$ and gradually reaches a saturation value in the inertial range. This behaviour confirms that the slow part of the transfer function is correlated across scales with a time lag, an indication of the time delay associated with the forward energy cascade. It is noteworthy that these results which we obtained for inertia-gravity waves in RSWE are synonymous with the results obtained for 3-D HIT by Khurshid et al. (Reference Khurshid, Donzis and Sreenivasan2021); see for instance their figures 1, 4 and 5. Therefore, just like in 3-D HIT, the turbulent downscale cascade of wave energy takes a certain time to reach dissipation scales from the domain scale, a result that is to be expected in hindsight since the downscale transfer is scale-local.

3.4. Temporal intermittency of energy transfers

As is clear from the discussions in §§ 3.2 and 3.3, the slow downscale energy transfer is scale-local and is associated with a finite time lag. The fast fluctuations on the other hand do not show noticeable correlations across scales and do not respect locality of transfers. Nevertheless, the fast fluctuations are associated with several features closely connected to temporal intermittency of transfers. Specifically, the temporal intermittency in transfers were seen to have close connections to physical structures and their time evolution, which we detail below.

For $f=1$, figure 9(a) shows the time series of spatial maxima of pointwise dissipation, $D^\epsilon$, and pointwise negative divergence, $- \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}$, for an arbitrary time interval. To plot these variables on the same figure we subtracted the mean of the time series from the instantaneous values and then divided the difference by the mean. On scanning figure 9(a), we notice that dissipation and negative divergence are in phase; notice that both of them attain high values at almost the same instances in time. In figure 9(b), we plotted the time series of spatial maxima of dissipation, i.e. the same black curve from figure 9(a), and the time series of energy flux at $k=200$, a wavenumber in the middle of the inertial range. Similar to the other curves in figure 9, we subtracted the mean flux from the instantaneous flux values and normalized by the mean flux to get the blue flux curve plotted in figure 9(b). The three different time series shown in figure 9 indicate that the variables fluctuate intermittently with time, with different local maxima. For convenient viewing we used a grey shading for different regions in figure 9(b) that contains local maxima of the variables. A reader who carefully examines figure 9(b) will notice that the local maxima of flux precedes the local maxima of dissipation; this feature can be seen in each grey shaded region of the figure, starting with the first maxima happening after $t = t_{1} + 2.5$. The flux attaining a maxima before the maxima of dissipation and negative divergence is another explicit illustration of the results discussed in § 3.3; i.e. the finite non-zero time associated with energy cascading across the inertial range to the dissipative scales.

Figure 9. Panel (a) shows the time variation of pointwise spatial maximum of dissipation ($D^\epsilon _{max}$) and negative divergence ($- { ( \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v} )}_{min}$). Panel (b) shows the time variation of pointwise spatial maximum of dissipation and energy flux at $k=200$. The inset in (a) shows the zoomed view of the region marked by the dashed blue lines in (a). The black dots in the top-right inset are chosen times at which the physical fields in figure 10 are shown. The first local maximum of the negative divergence happens at $t = t_{1} + 2.81$, this corresponding to the topmost black dot in the inset.

Alongside tracking the temporal fluctuations seen in figure 9, we monitored the details in physical space structures at different times. A key result was the finding that the local maxima of flux, negative divergence and dissipation went hand-in-hand with the intersection of a large number of shocks in physical space. To illustrate this, the inset at the top right-hand side of figure 9(a) shows a small part of the time series that is highlighted by the dashed blue lines in the bottom left of figure 9(a). In the inset we have marked five specific time instances by black dots, two before and two after the local maxima at $t = t_{1} + 2.81$. The physical space structures of the log of the absolute value of divergence and dissipation corresponding to these five time instances are shown in figure 10. A specific region of interest is highlighted by white dashed boxes in figure 10. On carefully examining figure 10 from figure 10(a,b) to figure 10(i,j), especially focusing on the region inside the white dashed boxes, the reader will find that different shocks move across the domain and intersect at specific locations. An incident that stands out is seen in figure 10(e,f), corresponding to $t = t_{1} + 2.81$, where twelve shocks intersect at a location where the maxima of negative divergence and dissipation are attained, these being seen in the inset of figure 9(a). Notice that the moving shocks converge to this point in space before this time instant and diverge from each other after this time instant. We found this feature to be persistent for different local maxima of the curves seen in the grey shaded regions of figure 9. Local maxima of flux, dissipation and negative divergence correspond to a large number of shocks intersecting in physical space, with their intersection being the location of maxima of dissipation and negative divergence.

Figure 10. The figure shows the $\log$ of the absolute value of the divergence field, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}$, (a,c,e,g,i) and the dissipation field, $D^\epsilon$, (b,d,f,h,j) at different times, corresponding to the times marked by black dots in the inset of figure 9(a). The dashed white boxes in the panels above highlight a specific region where a large number of shocks meet at $t = t_{1} + 2.81$, this being the time at which pointwise negative divergence and pointwise dissipation attain the maximum values, once again as seen from the top-most black dot at $t = t_{1} + 2.81$ in the inset of figure 9(a). Shocks move towards this intersection point at earlier times and shocks move away from this intersection point at later times, as can be seen by examining the dashed white box at different times from (a,b) to (i,j). Here (a,b) $t = t_{1} + 2.45$; (c,d) $t = t_{1} + 2.62$; (e,f) $t = t_{1} + 2.81$; (g,h) $t = t_{1} + 2.95$; (i,j) $t = t_{1} + 3.46$.

The above-described features were also seen for the other rotation rates and a key finding was that the level of intermittency increased with increasing rotation rate. To illustrate this, figure 11 shows the spectral flux for the $f=1$ case and $f=3$ case, the former shown in figure 11(a) and is the same blue curve that appears in figure 9(b). Notice that while we see multiple flux bursts from the mean for the $f=1$ case, only one distinguished burst in flux is seen for the $f=3$ case in figure 11(b) in the same time interval. Furthermore, the single burst seen in $f=3$ case is of much higher amplitude than that of the $f=1$ case. We noticed this to be a distinguishing feature of the turbulent transfers with increasing rotation rates: higher rotation rates were associated less frequent but larger amplitude bursts in flux, this being followed by multiple shocks intersecting in the domain and high dissipation and negative divergence values at the intersection points.

Figure 11. Figure shows the time variation of the spectral flux at $k=200$ for (a) $f = 1$ and (b) $f = 3$. Panel (a) shows the same blue curve in figure 9(b), replotted above for a convenient comparison.

We conclude this section by noting that the turbulent wave transfers are in general intermittent in time, with the temporal intermittent behaviour being correlated with distinguished features in the physical domain and with the temporal intermittency of the signals increasing with increase in the rotation rate.

4.1. Wave-balance energy exchanges

Applying the wave-balance decomposition to RSWE (2.3) with the dissipation term on the right-hand side, we get the evolution equation for the vortical flow variables. Taking the Fourier transform of those equations we can form the energy of the vortical flow as

(4.5)\begin{equation} \frac{\partial{ {{\hat{E} }_{\scriptscriptstyle{V}}} (\boldsymbol{k}, t) }}{\partial{t}} = {{ \hat{R}}_{\scriptscriptstyle{V}}} (\boldsymbol{k}, t) - {{ \hat{D}}_{\scriptscriptstyle{V}}} ( \boldsymbol{k}, t), \end{equation}

where the quadratic vortical energy is given by ${{ \hat {E} }_{\scriptscriptstyle {V}}} (\boldsymbol {k}, t) = \vert {{ \hat {u}}_{\scriptscriptstyle {V}}} (\boldsymbol {k}, t) \vert ^2 + \vert {{ \hat {v}}_{\scriptscriptstyle {V}}} ( \boldsymbol {k}, t) \vert ^2 + \vert {{ \hat {h}}_{\scriptscriptstyle {V}}} ( \boldsymbol {k}, t) \vert ^2$ and the corresponding dissipation is given by ${{\hat D}_{\scriptscriptstyle {V}}} (\boldsymbol {k}, t) = \nu \boldsymbol {k}^4 {{ \hat {E} }_{\scriptscriptstyle {V}}} (\boldsymbol {k}, t)$. Additionally, ${{ \hat {R}}_{\scriptscriptstyle {V}}}$ is the nonlinear transfer term that is responsible for energy exchanges. This term is similar to the ${{ \hat {R}}_{\scriptscriptstyle {W}}}$ term in (4.1), although unlike in (4.1) here we do not further decompose the term into different triadic contributions. Summing the above equation over all wavenumbers gives us the domain-integrated energy equation for the vortical flow as

(4.6)\begin{equation} \frac{{\rm d}{ {{E}_{\scriptscriptstyle{V}}} (t) }}{{\rm d}{t}} = {{R}_{\scriptscriptstyle{V}}} ( t) - {{D}_{\scriptscriptstyle{V}}} ( t). \end{equation}

As explained earlier, we used the equilibrated fields from $t=90\unicode{x2013}170$ to compute all the transfers and statistics described so far. During this same interval, we time-integrated the right-hand side terms in (4.6) to get the change in the vortical energy during the time interval $t=90\unicode{x2013}170$. For example, we found this energy change per unit area for the $f=1$ case to be $\delta {{E}_{\scriptscriptstyle {V}}}/A = 0.0117$, almost two orders lower than the vortical energy per unit area, which is close to 1. This drop in vortical energy is reported in the last column of table 2 for different rotation rates. From the table we see that the energy exchange with the wave field and vortical dissipation is negligible, leading to vortical energy changing very little during our numerical integrations. This is why we did not have to force the vortical field in our numerical integrations, although we forced the wave field.

The above discussion points out how weak energy exchange is between waves and vortical modes in the regime considered here for RSWE, where wave and vortical modes are comparable in strength. In the past, a broad set of studies have examined wave-balance energy exchanges in different contexts using asymptotic models, two-dimensional and 3-D models (Gertz & Straub Reference Gertz and Straub2009; Xie & Vanneste Reference Xie and Vanneste2015; Taylor & Straub Reference Taylor and Straub2016; Wagner & Young Reference Wagner and Young2016; Rocha, Wagner & Young Reference Rocha, Wagner and Young2018; Thomas & Yamada Reference Thomas and Yamada2019; Thomas & Arun Reference Thomas and Arun2020; Thomas & Daniel Reference Thomas and Daniel2020; Taylor & Straub Reference Taylor and Straub2020; Xie Reference Xie2020; Thomas & Daniel Reference Thomas and Daniel2021; Zhang & Xie Reference Zhang and Xie2023). These studies were primarily inspired by the need to identify whether different wave fields such as internal tides, near-inertial waves and internal wave continuum could act as an energy sink for oceanic balanced flow and the results and their consequences are reviewed in detail in Thomas (Reference Thomas2023). Based on the results from these studies, noticeable energetic interactions between waves and balanced flow of comparable strength is seen only in the case of small Burger number near-inertial waves. Other wave fields, like internal tides and internal wave continuum are seen to have negligible energy exchange with the balanced flow in regimes when the two fields have comparable strengths.

Based on the scaling we used in (2.2a–e), the Burger number, $Bu$, for the flow regimes is $Bu =(c_0/fL)^2 \sim O(1)$, this being the case for low mode internal tides. As a result, amongst the broad set of studies that have examined wave-balance energy exchanges mentioned above, energy exchanges examined by Thomas & Yamada (Reference Thomas and Yamada2019) in the $Bu \sim O(1)$ regime using a two-vertical-mode model that is similar to RSWE, where waves did not form shocks, is the closest to the energy exchange discussions in our present work. In that study, it was seen that in parameter regimes where waves and balanced flow have comparable strengths, the energy exchange between the two fields was negligible. In the present study, we explored interactions between waves and the vortical field in the regime where the two fields had comparable energy levels, since this is the regime where vortical scattering of waves is strongest and our goal was to examine how the vortical flow affects the waves’ downscale cascade. If we made the vortical energy weaker than waves, there would be less scattering of waves by the vortical flow and less modification of the waves’ cascade by the vortical flow. Consequently, we conclude by noting that the primary role of the vortical mode in the regime where wave and balanced energy levels are comparable is to scatter waves to smaller scales, with negligible energy exchange with the waves.