2. Simulations of decaying turbulence beneath surface waves

The incompressible, inviscid momentum equation in the presence of a background vorticity $\boldsymbol{\varOmega }$

(2.1) \begin{equation} {\partial }_t {\boldsymbol{u}} + \left ( {\boldsymbol{u}} \,{\boldsymbol{\cdot}}\, {\boldsymbol{\nabla }} \right ) {\boldsymbol{u}} + \boldsymbol{\varOmega } \times {\boldsymbol{u}} + {\boldsymbol{\nabla }} p = 0 , \qquad \text{with} \qquad {\boldsymbol{\nabla }} \,{\boldsymbol{\cdot}}\, {\boldsymbol{u}} = 0 , \end{equation}

describes both rotating flows without surface waves and also wave-averaged flows beneath steady, horizontally uniform surface-wave fields (see for example Craik & Leibovich Reference Craik and Leibovich1976; Holm Reference Holm1996; Suzuki & Fox-Kemper Reference Suzuki and Fox-Kemper2016; Wagner et al. Reference Wagner, Chini, Ramadhan, Gallet and Ferrari2021). In the wave-averaged case, $p$ and $\boldsymbol{u}$ in (2.1) represent the Eulerian-mean kinematic pressure and Lagrangian-mean velocity, respectively. The Lagrangian-mean velocity is the ‘total’ velocity responsible for the advection of mass, momentum, tracers and particles, and may be decomposed into an Eulerian-mean component plus ‘Stokes drift’ correction (van den Bremer & Breivik Reference van den Bremer and Breivik2018; Vanneste & Young Reference Vanneste and Young2022, for example). The wave-averaged momentum equation may be derived by asymptotic expansion leveraging a time-scale separation between rapidly oscillating surface waves and the slower evolution of $\boldsymbol{u}$ (Craik & Leibovich Reference Craik and Leibovich1976; Leibovich Reference Leibovich1980; Vanneste & Young Reference Vanneste and Young2022). The appearance of surface waves in (2.1) as an effective Coriolis force is discussed by Holm (Reference Holm1996).

For rotating flows, the background vorticity $\boldsymbol{\varOmega }$ in (2.1) takes the form

(2.2) \begin{equation} \boldsymbol{\varOmega }_{\textit{rotating}} = f {\boldsymbol{\hat{z}}} , \end{equation}

where $f$ is the Coriolis parameter and ${\boldsymbol{\hat{z}}}$ is the axis of rotation. For wave-averaged flows beneath a uniform surface-wave field propagating in the $x$ direction, the background vorticity $\boldsymbol{\varOmega }$ is alternatively

(2.3) \begin{equation} \boldsymbol{\varOmega }_{w\textit{a}v\textit{es}} = - {\boldsymbol{\nabla }} \times \boldsymbol u^{{S}} = - {\partial }_z {u^{{S}}} {\boldsymbol{\hat{y}}} , \end{equation}

where $\boldsymbol u^{{S}} = {u^{{S}}}(z) {\boldsymbol{\hat{x}}}$ is the surface-wave Stokes drift (Craik & Leibovich Reference Craik and Leibovich1976; Suzuki & Fox-Kemper Reference Suzuki and Fox-Kemper2016; van den Bremer & Breivik Reference van den Bremer and Breivik2018; Vanneste & Young Reference Vanneste and Young2022). We refer to ${\boldsymbol{\nabla }} \times \boldsymbol u^{{S}}$ – the curl of Stokes drift and therefore the curl of pseudomomentum Andrews & McIntyre (Reference Andrews and McIntyre1978) – as the ‘pseudovorticity’. (Stokes drift and pseudomomentum are identical for surface waves (McIntyre Reference McIntyre1981). Conceptually and more generally, ‘Stokes drift’ is the difference between the Eulerian-mean and Lagrangian-mean momentum, while ‘pseudomomentum’ is an apparent, wave-affiliated momentum that appears in the mean momentum budget.) In the context of (2.1), the main difference between rotating turbulence and turbulence beneath steady surface waves is the spatial structure of ${\partial }_z {u^{{S}}}$ : in shallow water ${\partial }_z {u^{{S}}} \sim z$ to leading order, while ${\partial }_z {u^{{S}}}$ decays exponentially for deep water waves.

2.1. Coherent structures in rotating and wave-averaged turbulence

To illustrate the similarity of rotating turbulence and turbulence beneath surface waves we conduct large eddy simulations of (2.1) with implicit grid-scale dissipation in a unit cube with three choices for $\boldsymbol{\varOmega }$

(2.4) \begin{equation} \boldsymbol{\varOmega }_{\textit{rotating}} = \frac {1}{4} {\boldsymbol{\hat{z}}}, \qquad \boldsymbol{\varOmega }_{w\textit{a}v\textit{es}} = -\frac {1}{2} z {\boldsymbol{\hat{y}}}, \quad \text{and} \quad \boldsymbol{\varOmega }_{\textit{isotropic}} = 0 . \end{equation}

Above, ${z}/{2}$ is the pseudovorticity associated with a shallow water wave in $z \in [0, 1]$ with Stokes drift $\boldsymbol u^{{S}} \approx ((1 + z^2/2) / 2){\boldsymbol{\hat{x}}}$ . Note that turbulence beneath shallow water waves is homeomorphic to $\beta$ -plane turbulence Rhines (Reference Rhines1975) – with ${\boldsymbol{\varOmega }}_{\beta\text{-}plane} = \beta y {\boldsymbol{\hat{z}}}$ – modulo on the $xz$ plane instead of on the $xy$ plane and with the Stokes drift curvature ${\partial }_z^2 {u^{{S}}}$ playing the role of $\beta$ .

The simulation domain is horizontally periodic in $x, y$ , and we use free-slip boundary conditions at the top and bottom boundaries in $z$ . Implicit dissipation of grid-scale kinetic energy is provided by an upwind-biased, nominally ninth-order weighted, essentially non-oscillatory advection scheme (Pressel et al. Reference Pressel, Mishra, Schneider, Kaul and Tan2017; Shu Reference Shu2020; Silvestri et al. Reference Silvestri, Wagner, Campin, Constantinou, Hill, Souza and Ferrari2024). Our implicit dissipation method allows us to simulate decaying turbulence that undergoes a transition from an isotropic to a nearly two-dimensional dynamics.

Following previous isotropic homogeneous turbulence studies (e.g. Orszag & Patterson (Reference Orszag and Patterson1972)), we produce an initial condition for the three simulations by conducting a preliminary simulation initialised with kinetic energy spectrum

(2.5) \begin{equation} \frac {1}{2} | {\boldsymbol{\hat{u}}} |^2 \sim | {\boldsymbol{K}} |^2 \exp \big \{-2 \left ( | {\boldsymbol{K}} | / K_i \right )^2 \big \} , \end{equation}

where ${\boldsymbol{\hat{u}}}$ is the Fourier transform of $\boldsymbol{u}$ , $\boldsymbol{K}$ is the Fourier wavenumber vector and $K_i = 32 \times 2\pi$ is the $32^{\text{}}$ nd wavenumber in the domain. The amplitude of the preliminary initial condition is scaled to produce an initial root-mean-square vorticity

(2.6) \begin{equation} \sqrt {\int | \boldsymbol{\omega } |^2 {\textrm{d}} V} \ {\stackrel {\textit{def}}{=}} \ \omega _{\textit{rms}}(t=-t_0) = 1000 , \quad \text{where} \quad \boldsymbol{\omega }\ {\stackrel {\textit{def}}{=}}\ {\boldsymbol{\nabla }} \times {\boldsymbol{u}} = \xi {\boldsymbol{\hat{x}}} + \eta {\boldsymbol{\hat{y}}} + \zeta {\boldsymbol{\hat{z}}} . \end{equation}

The initial simulation is run for a duration $t_0$ until $t=0$ , defined as the time when the mean-square vorticity has decayed to $\omega _{\textit{rms}}(t=0) = 10$ . The velocity field is then saved to disk to be used as an initial condition in subsequent runs starting from $t=0$ . The simulations are conducted with Oceananigans (Ramadhan et al. Reference Ramadhan, Wagner, Hill, Campin, Churavy, Besard, Souza, Edelman, Ferrari and Marshall2020; Wagner et al. Reference Wagner2025), which discretises (2.1) with a finite volume method. Scripts that reproduce simulations in this paper are stored on GitHub; see the data availability statement.

We focus first on the evolution of the relative vorticity $\boldsymbol{\omega }$ defined in (2.6). Figure 1 shows vorticity components for the three cases after $t=1000$ time units: figure 1(a) shows vertical vorticity $\zeta$ , while figures 1(b) and 1(c) show the horizontal vorticity $\eta$ . Figures 1(a) and 1(b) for rotating and wave-affected turbulence, respectively, both exhibit the formation of coherent structures and relatively greater vorticity levels than the unorganised, small-amplitude isotropic vorticity in figure 1(c). Figure 2 is similar, except that figure 2(a) shows $\zeta$ in rotating turbulence in the $xy$ -plane, while figures 2(b) and 2(c) show $\eta$ in the $xz$ plane for wave-averaged and isotropic turbulence, respectively.

Figure 1.

Vorticity in simulations decaying turbulence with $512^3$ finite volume cells in a unit cube after $1000$ time units. (a) Vertical relative vorticity $\zeta = {\boldsymbol{\hat{z}}} \,{\boldsymbol{\cdot}}\, ( {\boldsymbol{\nabla }} \times {\boldsymbol{u}} )$ in rotating turbulence with Coriolis parameter $f = 1/4$ , (b) horizontal relative vorticity $\eta = {\boldsymbol{\hat{y}}} \,{\boldsymbol{\cdot}}\, ( {\boldsymbol{\nabla }} \times {\boldsymbol{u}} )$ in turbulence beneath surface waves with Stokes shear ${\partial }_z {u^{{S}}} = - z / 2$ and (c) $\eta$ in isotropic turbulence.

Figure 2.

As figure 1 but showing (a) $\zeta$ for rotating turbulence in the $xy$ -plane; (b) $\eta$ for turbulence beneath surface waves in the $xz$ -plane and (c) $\eta$ for isotropic turbulence in the $xz$ -plane.

2.2. Zonation and the analogy of wave-averaged turbulence with beta-plane turbulence

One of the most striking aspects of the similarity between $\beta$ -plane turbulence and turbulence beneath surface waves is the development of ‘zonal jets’ – coherent, alternating jets in the direction of surface-wave propagation (and perpendicular to the pseudovorticity direction). This similarity allows us, therefore, to borrow intuition on structure formation in $\beta$ -plane turbulence (for example, Huang & Robinson Reference Huang and Robinson1998; Farrell & Ioannou Reference Farrell and Ioannou2007; Srinivasan & Young Reference Srinivasan and Young2012, Constantinou, Farrell & Ioannou Reference Constantinou, Farrell and Ioannou2014; for an overview see Constantinou Reference Constantinou2015; Farrell & Ioannou Reference Farrell, Ioannou, Galperin and Read2019; Marston & Tobias Reference Marston and Tobias2023).

To illustrate this, we consider two additional cases with background vorticity

(2.7) \begin{equation} \boldsymbol{\varOmega }_{\textit{deep}} = - \frac {1}{4} \textrm{e}^{8 (z - 1)} {\boldsymbol{\hat{y}}}. \qquad \text{and} \qquad \boldsymbol{\varOmega }_{w\textit{eak}} = - \frac {1}{8} z {\boldsymbol{\hat{y}}} . \end{equation}

These and subsequent simulations use $384^3$ finite volume cells to reduce their computational expense for the purpose of performing longer simulations to $t = 10^4$ .

Figure 3 shows time series of the $y$ -momentum $v$ , which plays the role that vertical velocity plays in rotating turbulence. The right side of figure 3 shows vertical profiles of the horizontally averaged $x$ -momentum

(2.8) \begin{equation} U(z, t)\ {\stackrel {\textit{def}}{=}} \int u \, \textrm{d} x \, \textrm{d} y . \end{equation}

These $U$ -profiles exhibit the development of depth-alternating jets, which are most prominently exhibited by the zig-zagging structure in the medium waves case on the right side of figure 3. To our knowledge, depth-alternating jets have not been observed in shear-free wave-modified turbulence. Comparing the slices of $v$ for medium and weak waves reveals how the medium-strong waves induce a strong inverse cascade, more coherent vortices and fewer small-scale motions than the weak waves.

Figure 3.

The evolution of cross-wave momentum $v$ in the $xz$ -plane at $t=20$ for (a) deep waves, (b) medium waves, and (c) weak waves, and at $t=400$ (e–g). The horizontally-averaged momentum $u$ is shown for all three cases in panel (d) at $t=20$ and panel (h) at $t=400$ . Dashed lines in the top right plot show the Stokes drift profile (normalised) for each case. The ‘medium’ waves case uses the background vorticity $\boldsymbol{\varOmega }_{w\textit{a}v\textit{es}}$ in (2.4). Both the $u$ profiles and the light grey ‘zero lines’ are spaced apart by the scale $\delta u = 10^{-2}$ . The $u$ -profiles at right show the development of depth-alternating jets, including a counter-wave surface jet for all cases. The $v$ -slices are similar between the three cases at early times, but at later times exhibit strong, localised wave impacts in their respective regions of significant Stokes shear. Note that $v$ plays the role that vertical velocity plays in rotating turbulence. Note that the Rhines scale may be estimated as $\ell _R = 2 \pi \sqrt {U / | {\partial }_z^2 {u^{{S}}} |}$ . With $ |{\partial }_z^2 {u^{{S}}} | = {1}/{2}$ and $U \approx 10^{-2}$ for the medium waves case, for example, we estimate $\ell _R \approx 0.8$ , consistent with the jet spacing apparent on the right panels.

3. Phenomenological model for the evolution of kinetic energy

We turn to the evolution of the domain-averaged kinetic energy

(3.1) \begin{equation} k(t)\ {\stackrel {\textit{def}}{=}} \int \frac {1}{2} \big ( u^2 + v^2 + w^2 \big ) {\textrm{d}} V . \end{equation}

Figure 4.

The decay of kinetic energy $k(t)$ in isotropic, rotating and surface-wave-modulated turbulence. Solid lines show kinetic energy normalised by it’s initial value, $k / k_0$ , computed from large eddy simulations. Dashed lines show $k / k_0$ given by (3.7), which solves the phenomenological two-equation system in (3.2)–(3.3). The initial dissipation rate $\epsilon _0$ in (3.7) is estimated from the numerical solution; specifically, we evaluate (3.8) at $\Delta t$ , substitute $k(t=\Delta t)$ from the numerical solution and solve for $\epsilon _0$ . We use $a = 11/6$ and estimate $b$ from (3.10). For the rotating case, $\varOmega = 1/4$ yields $b = 0.033$ . For the four surface-wave-modulated cases we use $b = 0.036$ , where $\varOmega \mathop {=}\limits ^{\textit{def}} \int {\partial }_z {u^{{S}}} {\textrm{d}} z = (1/4, 1/8, 1/16, 1/32 )$ . Note that using $a = 1.75$ along with commensurate adjustments to $b$ matches the simulation data even more closely.

Figure 4 plots time-series of the normalised kinetic energy $k(t) / k(t=0)$ for the three cases presented in figure 1, as well as three additional cases that use $\boldsymbol{\varOmega } = S z {\boldsymbol{\hat{y}}}$ with $S = (1/2, 1/8, 1/16)$ . Figure 4 illustrates another common feature to rotating turbulence and turbulence beneath surface waves: a suppression of the kinetic energy dissipation rate, such that at long times the kinetic energy levels off. We theorise that the suppression of dissipation is linked to the formation of coherent structures: rather than cascading to the grid scale, where dissipation can occur, kinetic energy accumulates in much larger-scale coherent structures. Similar to the dynamics of two-dimensional turbulence McWilliams (Reference McWilliams1984), the concentration of kinetic energy in coherent structures suppresses the forward cascade of kinetic energy.

Inspired by Bardina et al. (Reference Bardina, Ferziger and Rogallo1985), we model the kinetic energy $k(t)$ with the phenomenological two-equation system

(3.2) \begin{align} \frac {\textrm{d}}{\textrm{d} t} k &= -\epsilon , \end{align}

(3.3) \begin{align} \frac {\textrm{d}}{\textrm{d} t} \epsilon &= - a \frac {\epsilon ^2}{k} - b \varOmega \epsilon , \\[8pt] \nonumber \end{align}

where $\epsilon (t)$ describes the dissipation rate of $k$ , $\varOmega$ is a characteristic scale for $\boldsymbol{\varOmega }$ and $a$ and $b$ are $O(1)$ free parameters. Equation (3.2) follows directly from $\int {\boldsymbol{u}} \,{\boldsymbol{\cdot}}\, (2.1){\textrm{d}} V$ after accounting for the effect of the implicit numerical dissipation in our simulations. (Alternative models that include surface-wave-associated source terms for turbulent kinetic energy, such as the one proposed by Axell (Reference Axell2002), are inconsistent with (2.1)). The catalytic nature of rotation or surface waves may also be deduced from the fact that $\boldsymbol{\varOmega }$ does not appear in (3.2). The first term in (3.3) models the destruction of $\epsilon$ on the turbulent time scale $\tau = k / \epsilon$ , since $\textrm{d} \epsilon / \textrm{d} t \sim - \epsilon / \tau$ when $\varOmega \epsilon \ll \epsilon / \tau$ .

The second term in (3.3) models the suppression of kinetic energy dissipation – or alternatively, the growth of the correlation length due to the coalescence of coherent structures – by the background vorticity $\varOmega$ . The relative importance of ‘intrinsic’ destruction of $\epsilon$ and background-vorticity-induced destruction of $\epsilon$ is measured by the non-dimensional number $\varOmega \tau = \varOmega k / \epsilon$ , whose significance is revisited in § 4. We note that (3.2) is exact – and unchanged whether or not the system is rotating or modulated by surface waves. The modulation of turbulence by surface waves is therefore fundamentally ‘catalytic’, and can only affect $k$ indirectly by changing (3.3).

The free parameter $a$ may be constrained by considering isotropic turbulence with $\varOmega =0$ . In this case we expect the turbulent kinetic energy to decay according to $k \sim t^{-6/5}$ (Saffman Reference Saffman1967), which implies $ ({\textrm{d}}/{\textrm{d} t}) k = - ({6}/{5}) ({k}/{t})$ , and thus in turn, via (3.2), leads to

(3.4) \begin{equation} \epsilon = \frac {6}{5} \frac {k}{t} , \quad \text{so that} \quad \frac {\textrm{d}}{\textrm{d} t} \epsilon = -\frac {66}{25} \frac {k}{t^2} . \end{equation}

Inserting (3.4) into (3.3) yields $a = 11/6$ .

From (3.2)–(3.3) we deduce that the rate of change of the combination $\log \epsilon - a \log k$ is proportional to $b$ , that is

(3.5) \begin{gather} \frac {\textrm{d}}{\textrm{d} t} \left ( \log \epsilon - a \log k \right ) = - b \varOmega . \end{gather}

Integrating (3.5) produces

(3.6) \begin{equation} \frac {\epsilon }{\epsilon _0} =\left ( \frac {k}{k_0} \right )^a \textrm{e}^{- b \varOmega t } , \end{equation}

where $\epsilon _0$ and $k_0$ are the dissipation and kinetic energy at $t = 0$ , respectively. If we insert the expression for the dissipation $\epsilon (t)$ from (3.6) into (3.2) and integrate in time, we obtain

(3.7) \begin{equation} \frac {k}{k_0} = \left [ 1 + \frac {1}{n b \varOmega } \frac {\epsilon _0}{k_0} \big ( 1 - \textrm{e}^{-b \varOmega t} \big ) \right ]^{-n} , \end{equation}

where $n = 1 / (a - 1) = 6/5$ . (See also equation (10) in Bardina et al. (Reference Bardina, Ferziger and Rogallo1985).) Taking the limit $\varOmega \to 0$ (or $t \to 0$ with finite $\varOmega$ ), we obtain the corresponding solution for isotropic turbulence

(3.8) \begin{equation} k_{\textit{isotropic}}(t) = k_0 \left ( 1 + \frac {1}{n} \frac {\epsilon _0}{k_0} t \right )^{-n} , \end{equation}

which yields the expected power law $k \sim t^{-n}$ when $\epsilon _0 t/ n k_0 \gg 1$ .

At long times the isotropic and vortical solutions diverge: $k_{\textit{isotropic}} \to 0$ as $t \to \infty$ , while in the vortical case $k/k_\infty$ limits to the constant

(3.9) \begin{equation} \frac {k_\infty }{k_0} = \left ( 1 + \frac {\epsilon _0}{n b \varOmega k_0} \right )^{-n} , \end{equation}

where $k_\infty = \lim _{t \to \infty } k(t)$ . Equation (3.9) yields a formula for $b$ in terms of $k_\infty$

(3.10) \begin{equation} b \varOmega = \frac {\epsilon _0}{n k_0} \frac {1}{\left ( \frac {k_\infty }{k_0} \right )^{1/n} - 1} , \end{equation}

which we use to diagnose $b$ from our numerical simulations.

The dashed curves in figure 4 show solutions to the two-equation system in (3.2)–(3.3) that corresponds to the solid-line simulated results. Figure 4 shows that a single value of $b = 0.036$ qualitatively describes the evolution of kinetic energy beneath surface waves for a wide range of background vorticity magnitudes $\varOmega$ , where we estimate $\varOmega$ via the $z$ -average of ${\partial }_z {u^{{S}}}$ , that is $\varOmega \approx \int {\partial }_z {u^{{S}}} {\textrm{d}} z$ . For the rotating case, we use $b = 0.033$ and $\varOmega = 1/4$ . While qualitatively excellent considering that $b \approx 0.036$ describes a wide range of conditions, we also find that the phenomenological model overestimates the dissipation rate – and therefore underestimates the kinetic energy $k$ – during the transition between isotropic and background-vorticity-dominated regimes.