2.1. Statistical analysis

Statistics in this paper, as in Hass & Lele (Reference Hass and Lele2026), are based on horizontal and limited time-averaging, denoted by $\langle \boldsymbol{\cdot }\rangle$ . The time-average is restricted over a short enough time where statistics are approximately stationary. This quasistationary approximation is key to the wave–turbulence decomposition introduced below and is verified a posteriori. Based on this average, each dependent variable is decomposed into a mean and fluctuating component. The fluctuations are further divided into a wave and turbulence contribution. The nomenclature for this decomposition is given by

(2.4) \begin{align} h &= \langle h \rangle + h_{\!f}, \;\;\;\; h_{\!f} = \tilde {h} + h' \end{align}

where $h \in \{u_i, p_{nh}, T_{nh}, S_i, S_T, \tau _{ij}, q_{\!j}\}$ , subscript ‘ $f$ ’ denotes a fluctuation about the horizontal and temporal average $\langle h \rangle$ , a tilde $\tilde {h}$ is the wave component and a prime $h'$ is the turbulence. Given the forcing protocol, $\langle S_i \rangle = 0$ and we verified a posteriori there is negligible mean flow in the domain as a result. Because of this, we drop the subscript ‘ $f$ ’ on velocity fluctuations from now on, i.e. $u_i = u_{i,f}$ . The mean temperature is of course non-zero, $\langle T \rangle (z) \ne 0$ . In the following, we refer to ‘waves’ specifically as large-scale sloshing motions, so they do not contribute to fluctuations in $\tau _{ij}$ , $q_{\!j}$ , $S_i$ or $S_T$ , i.e. $\tilde {\tau }_{ij} = \tilde {q}_{\!j} = \tilde {S}_i = \tilde {S}_T = 0$ .

A subscript ‘ $I$ ’ denotes an interfacial quantity which is measured at the location where $\langle T \rangle =0.5$ . The coordinate associated with this location is denoted $z_I$ . Here $\langle T \rangle$ is averaged over the quasistationary period and so $z_I \lt z_i$ where $z_i$ is the initial interface location in (2.3) and reported in table 1. Both $z_i$ and $z_I$ are marked on a number of figures throughout the paper to orient the reader.

Table 1.

Domain size parameters. See figure 1 for definitions. Two mesh resolution levels are considered: medium (M) and fine (F). These represent mesh spacings of $3/128$ and $3/256$ , respectively. Statistics reported in this paper are computed from F-mesh results. However, all F-mesh simulations are initialized from M-mesh runs to efficiently bypass the early transient (propagation of the turbulent front from the forcing region to the density interface). Also, M-mesh simulations are used to obtain long time records for wavenumber–frequency analysis of interfacial waves.

Table 1. Long description

A table comparing domain size parameters for different simulation runs. The table has 3 rows and 11 columns. Column headers are Run I D, L x comma L y, z b, z t, z i, d, L s p, L f, N x comma N y, N z, and Resolution level. Row 1: Run I D, S M 74; L x comma L y, 6; z b, negative 4; z t, 2; z i, blank; d, blank; L s p, 0.6; L f, 1.5; N x comma N y, 512; N z, 512; Resolution level, F. Row 2: Run I D, S M 123 comma S M 124; L x comma L y, 9; z b, negative 5.025; z t, 1.35; z i, negative 1.75; d, 0.125; L s p, 0.6; L f, 1.5; N x comma N y, 384; N z, 272; Resolution level, M. Row 3: Run I D, S M 125 comma S M 126; L x comma L y, blank; z b, blank; z t, blank; z i, blank; d, blank; L s p, blank; L f, blank; N x comma N y, 768; N z, 544; Resolution level, F.

Table 2.

Input and output parameters of the simulations. Here $\textit{Re}$ and $\textit{Fr}$ depend on the reference length and velocity scales, $L^*$ and $U^*$ . More meaningful output parameters are reported as well where $\textit{Re}_t = \textit{Re} k^2/\epsilon$ , $\textit{Fr}_t = \epsilon /Nk$ , $\textit{Re}_b = \textit{Re} \epsilon /N^2$ and $N^2 = \partial _z\langle T \rangle /\textit{Fr}^2$ . Here $z_I$ is defined as the location where $\langle T \rangle = 0.5$ or $z_I=-1.75$ for SM74 (i.e. the initial density interface location). The output parameters are computed from time and planar-averages. The last column documents the amount of time-averaging used to compute statistics in terms of $\tau _{0.7}$ , the eddy turnover time where $\langle T \rangle = 0.7$ (stratified runs) or where $z' = 1.75$ (SM74).

Table 2. Long description

The table has five rows and eight columns. The columns are labeled as follows: Run ID, Re, Fr, Pr, Maximum Re(z), Fr(t), Re(t), Reb(z(t)), and Averaging window t(avg)/t(0.7). The rows present data for different simulation runs: SM74, SM123, SM124, SM125, and SM126. Each row lists specific values for the input and output parameters. For example, Row 1: Run ID, SM74; Re, 100 000; Maximum Re(z), 24 694; Fr(t), 1,659; Averaging window t(avg)/t(0.7), 4.90. Row 2: Run ID, SM123; Re, 100 000; Fr, 14; Pr, 7.2; Maximum Re(z), 25 319; Fr(t), 0.568; Re(t), 425; Averaging window t(avg)/t(0.7), 6.83. Row 3: Run ID, SM124; Re, 100 000; Fr, 7; Pr, 7.2; Maximum Re(z), 25 462; Fr(t), 0.227; Re(t), 99.5; Averaging window t(avg)/t(0.7), 5.80. Row 4: Run ID, SM125; Re, 100 000; Fr, 3.5; Pr, 7.2; Maximum Re(z), 25 359; Fr(t), 0.079; Re(t), 16.3; Averaging window t(avg)/t(0.7), 5.79. Row 5: Run ID, SM126; Re, 100 000; Fr, 1.25; Pr, 7.2; Maximum Re(z), 25 397; Fr(t), 0.009; Re(t), 1.01; Averaging window t(avg)/t(0.7), 2.49.

Plots of statistical quantities are truncated based on sufficient time-averaging. As in Hass & Lele (Reference Hass and Lele2026), at least two eddy turnover times are required to consider a statistical quantity reasonably converged.

2.2. Parameter space

The focus of the present paper is on five simulations: one unstratified and four of varying strengths of stratification. Each stratified simulation is initialized using a hyperbolic tangent temperature profile to approximate two-layer oscillating grid turbulence experiments as shown in figure 1. Also depicted in the schematic are two sponge layers that are used to damp turbulent fluctuations before reaching the domain boundary. (Note, IGWs, a large focus of the present paper, are confined to the interface and so thin sponge layers can be used without concern of reflecting IGWs off the domain boundary.) Dimensions of the computational domain are given in table 1, which refers to the dimensions depicted in figure 1.

Two mesh resolutions are reported in table 1, medium resolution (M) and fine resolution (F). Statistics are computed from F-mesh simulations whereas M-mesh runs are used to obtain long time records to resolve low-frequency interfacial modes. Also, F-mesh simulations are initialized from quasistationary M-mesh simulations to bypass the early transient. See Hass (Reference Hass2025) and Hass & Lele (Reference Hass and Lele2026) for further information including a detailed mesh-convergence study.

The parameter space spanned by the simulations is recorded in table 2. The input parameters result in outputs that give a sense of the overall stability (stratification) and turbulence level in the domain. The turbulent Reynolds number $\textit{Re}_t = \textit{Re} k^2/\epsilon$ is the ratio of inertial-to-viscous forces and characterizes the range of isotropic scales in unstratified turbulence. Here $k\equiv \langle u_i u_i \rangle /2$ is the (non-dimensional) fluctuating kinetic energy and $\epsilon \equiv 2 \langle s_{ij} s_{ij} \rangle /\textit{Re}$ the kinetic energy dissipation rate where $s_{ij} = (\partial _{\!j} u_i + \partial _i u_{\!j})/2$ is the strain rate tensor. The turbulent Froude number $\textit{Fr}_t = \epsilon /Nk$ is the ratio of the buoyancy time scale $\tau _b = 1/N$ to the turbulence time scale $\tau _t = k/\epsilon$ where $N = \sqrt {\partial _z \langle T \rangle } / \textit{Fr}$ is the dimensionless buoyancy frequency. Small $\textit{Fr}_t$ indicates that buoyancy effects are dominant. The Reynolds number $\textit{Re}_t$ is maximum in the forcing region and then relaxes to a nearly constant value for unstratified simulations (Hass & Lele Reference Hass and Lele2026). As defined, $\textit{Re}_t$ is sensitive to interfacial waves and so interfacial turbulence is best characterized by the buoyancy Reynolds number $\textit{Re}_b = \textit{Re} \epsilon /N^2$ which quantifies the time scale separation between buoyancy ( $\tau _b = 1/N$ ) and the viscous time scale $\tau _{\eta } = 1/\sqrt {\textit{Re} \epsilon }$ . (Here $\textit{Re}_b$ is a ratio of a molecular quantity $\epsilon$ and a mean quantity $N^2$ , neither of which are contaminated by large-scale wave motion.) For buoyancy-dominated regions, $\textit{Re}_b$ characterizes the dynamic range of turbulent eddies. For weakly stratified regions, it is not physically meaningful (Mater & Venayagamoorthy Reference Mater and Venayagamoorthy2014) and so we only report $\textit{Re}_b$ at the interface.

3.1. Unforced evolution of the interface

Consider case SM126 ( $\textit{Fr} = 1.25$ ) which is the most strongly stratified simulation (see table 2). We turned the forcing off at a late time when turbulence was fully developed and the system was quasistationary and allowed the system to evolve under its own inertia. Turbulence rapidly decayed, but interfacial modes continued oscillating long after small scale mixing all but shut off. Figure 2(a) shows instantaneous contours of vertical velocity and temperature fluctuations at a time just before turning off the forcing and a much later time. Note the magnitude of the large-scale $w$ oscillations are largely unchanged between the two snapshots whereas small scale turbulence has almost completely dissipated. As noted in § 1, the vertical structure of the temperature modes is much smaller than the velocity modes, which indicates that using a vertical filter to separate waves and turbulence as proposed in Torres et al. (Reference Torres, Klein, D’Asaro, Wang, Thompson, Siegelman, Menemenlis, Rodriguez, Wineteer and Perkovic-Martin2022) is not appropriate for the temperature field.

To quantify the energy of interfacial waves we integrated the fluctuating kinetic energy below the forcing layer edge and plotted its time evolution in figure 2(b). Overlayed is the integrated dissipation rate (dashed black line), where each quantity has been normalized by its initial value (i.e. the value before turning the forcing off). These domain-integrated quantities are defined as follows:

(3.1) \begin{align} E_k &= \int _{z_b}^{z_{\!f}} \int _0^{L_y} \int _0^{L_x} \frac {u_i u_i}{2} {\rm d}x {\rm d}y {\rm d}z, \\[-12pt]\nonumber \end{align}

(3.2) \begin{align} \epsilon _V &= \int _{z_b}^{z_{\!f}} \int _0^{L_y} \int _0^{L_x} \frac {2 s_{ij} s_{ij}}{\textit{Re}} {\rm d}x {\rm d}y {\rm d}z, \end{align}

where $z_{\!f}$ corresponds to the bottom edge of the forcing region ( $z_{\!f} = -0.75$ in all simulations reported here) and $z_b$ is the bottom boundary of the domain.

Time in figure 2(b) is normalized by $\tau _I$ , the initial eddy turnover time measured at the interface. (This was calculated based on the kinetic energy and dissipation rate measured where $\langle T \rangle =0.5$ before the forcing was turned off. This gives an overestimate of the eddy turnover time since the sloshing modes contain significant kinetic energy. We therefore took 40 % of $k/\epsilon$ to estimate $\tau _I$ in these plots which is the value of the ratio $k'/k$ at $z_I$ as demonstrated in § 5.2.) At five eddy turnover times, the total kinetic energy in the domain has dropped to 7.3 % of its initial value whereas the dissipation rate decreased to 0.056 % demonstrating that the sloshing motions persist long after the majority of the turbulence has been dissipated. If we assume that viscous dissipation is negligible for the sloshing modes then this suggests that the waves contribute approximately 7 % to the integrated kinetic energy. The local contribution near the interface is even greater (approximately 60 % for this case – see § 5.2) highlighting the difficulty in interpreting statistics near the interface.

Figure 2.

System evolution without external forcing. (a) Instantaneous contours of (i,iii) vertical velocity and (ii,iv) temperature fluctuations. Subpanels (i) and (ii) are at the initial condition and (iii) and (iv) are at a later time when almost all of the turbulence has decayed providing an approximate representation of the interfacial sloshing motions. (b) Time history of integrated kinetic energy and the kinetic energy dissipation rate. Each quantity is integrated over the subdomain below the bottom edge of the forcing region and normalized by its value at $t=0$ (subscript `0’).

Further evidence of the sloshing modes comes courtesy of the two-point correlation of quasistationary data (forcing turned ‘on’). We can define an autocorrelation coefficient as

(3.3) \begin{align} c_{ij,k}(r,z) &= \frac {\langle u_i(\boldsymbol{x},t)u_{\!j}(\boldsymbol{x} + r\boldsymbol{e}_k,t) \rangle }{\sqrt {\langle u_i^2(\boldsymbol{x},{t}) \rangle \langle u_{\!j}^2(\boldsymbol{x} + r\boldsymbol{e}_{k},t) \rangle }} \end{align}

where the bracket $\langle \boldsymbol{\cdot }\rangle$ denotes an average in $x$ , $y$ and time (see § 2.1). The correlation of vertical velocity with vertical separation $c_{33,3}$ is particularly sensitive to sloshing motions and an ideal diagnostic for their presence. Figure 3 shows this quantity for the SM126 case at various $z$ -locations. Due to statistical inhomogeneity in $z$ , the correlation is not symmetric and negative $r$ values show the correlation with points below the observation point. Each of the observation points plotted in the figure are above the density interface (see figure 3 b) and so negative $r$ values show the correlation with points in the direction of the interface. The effect of interfacial sloshing is manifest as large plateaus in the negative $r$ direction. As the observation location gets closer to the interface, the magnitude of these plateaus increases, consistent with the picture of large-scale interfacial motions.

Figure 3.

(a) Two-point correlation $c_{33,3}(r,z)$ for SM126 ( $\textit{Fr} = 1.25$ ). Three observer locations above the density interface are shown. Circles mark the location where $\langle T \rangle =0.5$ . (b) Observer locations visualized relative to the mean temperature profile (black curve). Line colours correspond to the legend in (a). For reference, the initial interface location is at $z = -1.75$ .

3.2. Linearized equations as a model for interfacial waves

The essentially inviscid nature of the interface sloshing suggests that oscillating modes may be well-approximated by the linearized equations governing small-amplitude, inviscid perturbations. This approximation becomes exact as $\textit{Fr} \to 0$ . Neglecting nonlinear and viscous/diffusive terms in (2.1), the linear model is given by

(3.4) \begin{align} \frac {\partial u_{\!j}}{\partial x_{\!j}} &= 0, \end{align}

(3.5) \begin{align}\frac {\partial u_i}{\partial t} + \frac {\partial p_{nh}}{\partial x_i} - \frac {T_{nh}}{\textit{Fr}^2} \delta _{i3} &= 0, \end{align}

(3.6) \begin{align} \frac {\partial T}{\partial t} + w\textit{Fr}^2 N^2 &= 0. \end{align}

Where the dependent variables are understood to be small-amplitude perturbations. Here $\textit{Fr}$ and $N^2$ are defined in § 2. The equations can be reduced to a fourth-order partial differential equation for the vertical velocity perturbation $w$ (Gill Reference Gill1982),

(3.7) \begin{align} \frac {\partial ^2 }{\partial t^2}{\nabla} ^2 w + N^2 {\nabla} _h^2 w &= 0 \end{align}

where ${\nabla} _h^2 = {\partial ^2 }/{\partial x^2} + {\partial ^2 }/{\partial y^2}$ is the horizontal Laplacian operator. Assuming normal mode solutions

(3.8) \begin{align} w \sim \hat {w}(z) e^{i(\alpha x + \beta y - \omega t)} + \text{c.c.} \end{align}

results in a linear homogeneous ordinary differential equation (ODE) for $\hat {w}(z)$ ,

(3.9) \begin{align} \frac {d^2 \hat {w}}{{d}z^2} - \kappa _h^2\hat {w} + \frac {1}{\omega ^2}N^2(z)\kappa _h^2 \hat {w} &= 0 \end{align}

Here $\alpha$ and $\beta$ are the wave-vector components in the $x$ and $y$ directions respectively, $\omega$ is the angular frequency, and $\kappa _h = \sqrt {\alpha ^2 + \beta ^2}$ is the horizontal wavenumber. ‘c.c.’ denotes the complex conjugate in (3.8) which was omitted in (3.9) for simplicity and will be dropped from now on.

The depth dependence of $N^2$ is made explicit in equation (3.9) demonstrating that it is valid for arbitrary background density profiles. When the boundary conditions are homogeneous Dirichlet conditions (i.e. $\hat {w} = 0$ at the boundaries, the so-called ‘rigid-lid approximation’ (Gill Reference Gill1982)), this becomes a regular Sturm–Liouville eigenvalue problem, solutions of which yield an infinite set of eigenfunctions $\hat {w}_n$ and eigenfrequencies $\omega _n$ for each horizontal wavenumber $\kappa _h$ , i.e. $\omega _n = \omega _n(\kappa _h, N^2(z))$ and $\hat {w}_n = \hat {w}_n(z; \kappa _h, N^2(z))$ $\forall n \in \{1, 2, \ldots \}$ . Sturm–Liouville theory guarantees that the eigenfrequencies are unique and the eigenfunctions are orthogonal under the weighting function $N^2(z)$ , i.e.

(3.10) \begin{equation} [\hat {w}_n, \hat {w}_m]_{N^2} \equiv \int _z N^2(z) \hat {w}_n(z) \hat {w}_m^*(z) {\rm d}z = \mathcal{W}_n \delta _{nm} \end{equation}

where $\mathcal{W}_n = [\hat {w}_n, \hat {w}_n]_{N^2}$ , $\delta _{nm}$ is the Kronecker delta and, abusing notation, $*$ henceforth denotes the complex conjugate (dimensional quantities will not be used from now on). (Solutions to (3.9) are real-valued, but inner products of complex-valued functions are used below and so the general definition is presented here.) We have ordered the set of eigenfrequencies and eigenfunctions such that $\omega _1 \gt \omega _2 \gt \ldots$ . The eigenfunction $\hat {w}_n$ represents the vertical structure of a wave associated with frequency $\omega _n$ . Note, the no-penetration boundary condition on $w$ , while an approximation for the ocean free surface, is the logical choice for our numerical simulations.

When $N^2$ is an arbitrary function of $z$ , the equations must be solved numerically. This was done using a pseudospectral method: $\hat {w}$ is expanded as a sine series to enforce the no-penetration boundary conditions, i.e.

(3.11) \begin{align} \hat {w}_n(z) &= \sum \limits _{j=1}^{N_z} W_{\!j} \sin {\bigg ( j \pi \frac {z-z_b}{z_t - z_b} \bigg )}. \end{align}

Substitution of this ansatz into (3.9), for each discrete $z_k$ ( $k \in \{1, \ldots , N_z\}$ ) results in a linear system of equations for the eigenfrequencies $\omega _n$ and eigenvectors $\boldsymbol{W}$ , from which the function $\hat {w}_n$ can be reconstructed according to (3.11). The method is pseudo-spectral because $N^2$ was not represented in a sine basis, rather its pointwise values at each $z_k$ were used. This results in a full coefficient matrix which is easily inverted on a personal laptop given the size of our problem.

The mode shapes for the first two modes of $\kappa _h = \sqrt {2} (2\pi /L_x)$ are shown in figure 4 where the eigenmodes $\hat {u}$ and $\hat {T}$ are related to $\hat {w}$ through the linearized equations. Substituting normal mode representations into (3.4)–(3.6), the following equations, known as polarization relations, are easily derived:

(3.12) \begin{align} \hat {T}_n(z) &= \frac {-i \textit{Fr}^2 N^2(z)}{\omega _n} \hat {w}_n(z), \end{align}

(3.13) \begin{align} \hat {p}_n(z) &= \hat {p}(z_b) + \int _{z_b}^z i\omega _n \bigg (1 - \frac {N^2(\sigma )}{\omega _n^2} \bigg ) \hat {w}_n(\sigma ) {\rm d} \sigma, \end{align}

(3.14) \begin{align} \hat {u}_{n}(z) &= \frac {\alpha }{\omega _n}\hat {p}_n(z), \end{align}

(3.15) \begin{align} \hat {v}_{n}(z) &= \frac {\beta }{\omega _n}\hat {p}_n(z), \end{align}

where $\hat {p}_n(z_b)$ is imaginary and determined by solving the Poisson equation for pressure. Because $\hat {w}$ is real-valued, $\hat {p}$ , $\hat {T}$ , $\hat {u}$ and $\hat {v}$ are imaginary, hence the imaginary components of $\hat {T}$ and $\hat {u}$ are plotted in figure 4. The normalization condition employed in this paper is $\max \limits _z \{|\hat {w}_n|\} = 1$ .

Figure 4.

Mode shapes for $\kappa _h L_x/2\pi = \sqrt {2}$ and $n=1$ (solid) and $n=2$ (dashed) for SM126. For reference, the initial interface location is at $z = -1.75$ .

Figure 4. Long description

Panel A: A line graph shows velocity data with the z-axis ranging from −5 to 1. The x-axis represents velocity values from -1 to 2. The graph includes two lines: a solid blue line labeled w-hat and a dashed orange line labeled Imaginary-u-hat. The solid blue line shows a peak around 0 on the x-axis and then decreases sharply. The dashed orange line also peaks around 0 but with a different pattern, showing a more pronounced dip before stabilizing. Panel B: Another line graph displays temperature data with the z-axis ranging from −5 to 1. The x-axis represents temperature values from −10 to 10. This graph features a solid blue line labeled Imaginary-T-hat. The line shows a significant dip around 0 on the x-axis and then stabilizes, with a dashed blue line indicating another component of the data.

3.2.1. Linear model comparison with nonlinear dynamics

We now ask the question: How well do solutions to (3.9) represent interfacial motions seen in the simulations? To help answer this, we compute $E_{ww}(\kappa _h,z,\omega )$ , the vertical kinetic energy wavenumber–frequency spectra for the two most strongly stratified cases SM125 and SM126 ( $\textit{Fr} = 3.5$ and 1.25). Here $E_{ww}$ is defined as

(3.16) \begin{align} E_{ww}(\kappa _h,z,\omega ) = \langle \check {w}(\alpha ,\beta ,z,\omega ) \check {w}^*(\alpha ,\beta ,z,\omega ) \rangle _{\boldsymbol{\kappa }_h} \end{align}

where the averaging is taken over annuli in the $\alpha$ – $\beta$ plane, and $\check {w}$ is the horizontal and time Fourier transform of vertical velocity. The result of this calculation is shown in figure 5 for a $z$ -location below the interface midline where turbulence has been significantly attenuated due to buoyancy. Overlayed in white is the dispersion relation from solving (3.9). Apparently, the interface sloshing is well-approximated by the $n=1$ and $n \in \{1,2\}$ modes for SM125 and SM126, respectively.

Figure 5.

Dispersion relation comparison from the LES data and solutions to (3.9): (a) SM125 ( $\textit{Fr} = 3.5$ ); (b) SM126 ( $\textit{Fr} = 1.25$ ). The coloured contours are the $\kappa _h$ – $\omega$ spectra computed from the LES data and the white lines are the dispersion relations for the first mode and the first two modes (a and b, respectively) from solving (3.9). Spectra are visualized at $z = -1.9664$ which is below the $\langle T \rangle =0.5$ location in both cases: $z_I = -1.8351$ for SM125 and $-1.7648$ for SM126 where $z_I$ is defined as $\langle T \rangle (z_I)=0.5$ .

Figure 6.

Instantaneous eigenfrequencies. Frequencies are plotted for the first four horizontal wavenumbers and the $n=1$ mode over the quasistationary time record. Time is normalized by the eigenperiod.

Figure 6. Long description

The image contains four line graphs, each representing the relationship between normalized frequencies and time for different horizontal wavenumbers and Froude numbers. Panel A: The line graph shows the normalized frequency (ω1/Nmax) over time (t〈ω〉/2π) for κhLx/2π = 1. The x-axis represents time normalized by the eigenperiod, and the y-axis represents the normalized frequency. Different lines represent different Froude numbers (Fr = 14, 7, 3.5, 1.25). Panel B: The line graph shows the normalized frequency (ω2/Nmax) over time for κhLx/2π = √2. The axes and lines are labeled similarly to Panel A. Panel C: The line graph shows the normalized frequency (ω3/Nmax) over time for κhLx/2π = 2. The axes and lines are labeled similarly to Panel A. Panel D: The line graph shows the normalized frequency (ω4/Nmax) over time for κhLx/2π = √5. The axes and lines are labeled similarly to Panel A. Each graph shows how the normalized frequency varies over time for different Froude numbers and horizontal wavenumbers.

3.3. Projection procedure using orthogonal eigen-functions

Solutions to (3.9) are orthogonal, as already discussed. Using (3.12) it is easy to show that $\hat {T}$ is also orthogonal under the weight $1/N^2$ , i.e. $[{\hat {T}}_n, {\hat {T}}_m]_{{1}/{N^2}} = [{\hat {T}}_n, {\hat {T}}_n]_{ {1}/{N^2}} \delta _{nm}$ where the temperature inner product is defined as

(3.17) \begin{align} [{\hat {T}}_n, {\hat {T}}_m]_{\frac {1}{N^2}} &= \int _z \frac {1}{N^2} {\hat {T}}_n {\hat {T}^*}_m {\rm d}z. \end{align}

Using $\hat {w}$ and $\hat {T}$ as a basis, the wave component of velocity and temperature can be represented as eigenfunction expansions,

(3.18a) \begin{align} \boldsymbol{\check {\tilde {u}}}(\alpha , \beta , z, t) &= \sum \limits _n a_n^{(w)}(t; \alpha , \beta ) \boldsymbol{{\hat {u}}}_n(z;\alpha , \beta ), \\[-12pt]\nonumber \end{align}

(3.18b) \begin{align} \check {\tilde {T}}(\alpha , \beta , z, t) &= \sum \limits _n a_n^{(T)}(t; \alpha , \beta ) {\hat {T}}_n(z; \kappa _h), \end{align}

where $\boldsymbol{\check {\tilde {u}}}$ and $\check {\tilde {T}}$ are the horizontal Fourier transforms of $\boldsymbol{\tilde {u}}\equiv (\tilde {u}, \tilde {v}, \tilde {w})$ and $\tilde {T}$ . We explicitly show the parametric dependence of $a_n^{(w)}$ , $a_n^{(T)}$ , $\hat {\boldsymbol{u}}_n$ and $\hat {T}_n$ on the horizontal wavevector components (and/or magnitude) as a reminder that a projection must be done for each Fourier mode resulting in different coefficients for each.

The expansion coefficients $a_n^{(w)}$ and $a_n^{(T)}$ are computed by projecting the LES data onto the eigenbasis:

(3.19) \begin{align} a_n^{(w)}(t;\alpha ,\beta ) &= \frac {[\check {w}, \hat {w}_n]_{N^2}}{[\hat {w}_n, \hat {w}_n]_{N^2}}, \\[-12pt]\nonumber \end{align}

(3.20) \begin{align} a_n^{(T)}(t;\alpha ,\beta ) &= \frac {[\check {T}, \hat {T}_n]_{\frac {1}{N^2}}}{[\hat {T}_n, \hat {T}_n]_{\frac {1}{N^2}}}. \end{align}

Note, $a_n^{(w)}$ is vertical velocity projected onto $\hat {w}_n$ which is then used to expand all velocity components in (3.18a ), thus $a_n^{(w)}$ is a scalar coefficient in the vector expansion. To reiterate, $\check {w}$ and $\check {T}$ in the inner products in (3.19) and (3.20) are the full nonlinear LES data.

The existence of two sets of expansion coefficients to represent the wave fields is natural. We noted previously the vertical structure of the $\hat {w}$ modes, which is set by the mixed layer depth, is much larger than the temperature modes which are proportional to the interface thickness (see figures 2 a and 4) and so both modes are required to fully describe wave dynamics. Furthermore, we show in Appendix B that the $a_n^{(w)}$ coefficients are the kinetic energy mode amplitudes while the $a_n^{(T)}$ are proportional to the potential energy mode amplitudes, so both are required for a total accounting of the energy.

Equations (3.18)–(3.20) are well-defined and consistent mathematical operations. The question is whether the result has anything to do with the linearized equations and, more importantly, with the interfacial waves we posit are modelled by such equations.

4.1. Relating the projection coefficients $a_n^{(w)}$ and $a_n^{(T)}$

We begin by taking the time derivative of $a_n^{(T)}$ (3.20),

(4.1) \begin{align} \dfrac {{\rm d} a_n^{(T)}}{{\rm d}t} = \dfrac {\int _z \dfrac {1}{N^2} \dfrac {\partial T}{\partial t} \hat {T}^*_n {\rm d}z}{\int _z \dfrac {1}{N^2}\hat {T}_n \hat {T}^*_n {\rm d}z} = \dfrac {\int _z \dfrac {1}{N^2} \bigg ( -w \dfrac {\partial \langle T \rangle }{\partial z} \bigg ) \bigg ( \dfrac {i \textit{Fr}^2 N^2}{\omega _n} \hat {w}_n \bigg ) {\rm d}z}{\int _z \dfrac {1}{N^2} \bigg ( \dfrac {-i\textit{Fr}^2 N^2}{\omega _n} \hat {w}_n\bigg ) \bigg ( \dfrac {i\textit{Fr}^2 N^2}{\omega _n} \hat {w}_n\bigg ) {\rm d}z} = -i\omega _n a_n^{(w)}, \end{align}

where the second equality is the result of using (3.6) to replace $\partial _t T$ (assuming linear dynamics) and (3.12) to eliminate $\hat {T}_n$ . Similarly, we show in Appendix B that

(4.2) \begin{align} \frac {{\rm d} a_n^{(w)}}{{\rm d}t} = -i \omega _n a_n^{(T)}. \end{align}

This gives a set of simple first-order ODEs relating the two projection coefficients:

(4.3) \begin{equation} \frac {{\rm d} a_n^{(w)}}{{\rm d}t} = -i \omega _n a_n^{(T)}; \;\;\;\; \frac {{\rm d} a_n^{(T)}}{{\rm d}t} = -i \omega _n a_n^{(w)}. \end{equation}

The system (4.3) can be decoupled by taking another time-derivative of each ODE and using (4.3) to eliminate first-order derivatives:

(4.4) \begin{equation} \frac {d^2 a_n^{(w)}}{{\rm d}t^2} + \omega _n^2 a_n^{(w)} = 0; \;\;\;\; \frac {d^2 a_n^{(T)}}{{\rm d}t^2} + \omega _n^2 a_n^{(T)} = 0. \end{equation}

In words, both sets of coefficients are harmonic oscillators with frequency $\omega _n$ . Therefore, the projection procedure does in fact sift out the eigenmodes of the linearized equations: they have the right spatial structure by using the eigenfunction basis and they have the right temporal dynamics. This assertion is true to the extent that the projected data satisfy the linearized equations which was assumed in the derivation of (4.4). In other words, we can test the validity of this result by computing the left- and right-hand sides of (4.3) for our data and seeing whether the equality holds. Figure 7 shows this result for the first mode of the first horizontal wavenumber and for the two low- $\textit{Fr}$ cases. The coefficients are complex-valued and so only the real component is shown, but similar curves are obtained for the imaginary component.

Table 3.

Error metrics for the curves shown in figure 7.

Table 3. Long description

A table with two rows and three columns comparing error metrics for two cases. The columns are labeled Case, E1W (2π/Lx), and E1T (2π/Lx). The first row represents Case SM125 with values 0.4138 and 0.0381. The second row represents Case SM126 with values 0.1891 and 0.0168.

Figure 7.

Panels (i) and (iii) compare the left- and right-hand sides of $\dot {a}_n^{(w)}=-i\omega _n a_n^{(T)}$ ; panels (ii) and (iv) compare the left- and right-hand sides of $\dot {a}_n^{(T)} = -i\omega _n a_n^{(w)}$ . The time derivative terms (left-hand sides) are plotted as solid blue lines and the dashed red lines correspond to the right-hand sides. Only the real component computed using the $\kappa _hL_x/2\pi = 1$ and $n=1$ mode is shown: (a) SM125 ( $\textit{Fr} = 3.5$ ); (b) SM126 ( $\textit{Fr} = 1.25$ ).

Clearly, the linearized model represents our LES data very well. Interestingly, the equation for $a_1^{(T)}$ has better agreement with the data than that for $a_1^{(w)}$ . This suggests that the projected data satisfies the linearized heat equation better than that for vertical momentum. However, drawing general conclusions from a single curve can be misleading. Defining an error metric enables quantification of the disagreement between all curves. First, we define the complex residuals,

(4.5) \begin{align} r_n^w &\equiv \dot{a}_n^{(w)} + i\omega_n a_n^{(T)} \nonumber\\r_n^T &\equiv \dot{a}_n^{(T)} + i\omega_n a_n^{(w)}. \end{align}

We then define a real-valued error metric that quantifies the disagreement separately for the real and imaginary parts of the curves depicted in figure 7:

(4.6) \begin{align} E_n^w(\kappa_h) &= \bigg \langle{\frac{\sqrt{\langle{[\mathcal{P}(r_n^w)]^2}\rangle_t}}{\langle{| \mathcal{P}(\dot{a}^{(w)}_n)|} \rangle_t}}\bigg \rangle_{\boldsymbol{\kappa}_h, \mathcal{P}} \nonumber\\ E_n^T(\kappa_h) &= \bigg \langle{\frac{\sqrt{\langle{[\mathcal{P}(r_n^T)]^2}\rangle_t}}{\langle{| \mathcal{P}(\dot{a}^{(T)}_n)|}\rangle_t}}\bigg \rangle_{\boldsymbol{\kappa}_h, \mathcal{P}}. \end{align}

Dots denote time differentiation, and the particular average is indicated by a subscript. A subscript ‘ $t$ ’ denotes a time average and ‘ $\boldsymbol{\kappa }_h$ ’ an azimuthal average in the $\alpha$ – $\beta$ plane (one error metric for each horizontal wavenumber). The operator $\mathcal{P} \in \{\text{Re}, \text{Im}\}$ extracts the real or imaginary part, while a subscript $\mathcal{P}$ denotes an arithmetic average over the errors obtained from the two parts, yielding a real-valued error metric.

For the present data, the numerical value is non-negative and approximately bounded above by unity. Table 3 shows the value for each curve in figure 7 to assist in interpreting the numerical value. To get a broader perspective, the metrics are computed for the first 60 wavenumbers and four modes and visualized in figure 8. The metrics corresponding to the curves in figure 7 can be located at the bottom-left-hand corner of each figure. Consistent with our expectation, the error metric is smallest across the board for SM126 compared with SM125. In fact, only SM126 has reasonably low error for the $n=2$ mode consistent with the existence of a second band of spectral peaks seen in figure 5(b). It is clear that the higher-order modes are not relevant to the LES data (high error for larger $n$ ). It is interesting that for low wavenumbers and $n=1$ , $E^T_n \lt E_n^w$ and is quite low for both levels of stratification as seen in the individual curves of figure 7, again, suggesting that the linearized temperature equation is a better model for our projected data than the linearized vertical momentum equation.

Figure 8.

Error metrics quantifying the disagreement with the linearized equations for the first four modes and 60 horizontal wavenumbers: (a) SM125 ( $\textit{Fr} = 3.5$ ); (b) SM126 ( $\textit{Fr} = 1.25$ ).

By only considering the error, we get the impression that the accuracy (and thereby, the utility) of the projection procedure is sensitive to which modes are included in the basis. However, if a particular mode results in large error, but the projection captures negligible energy, then the turbulence has not been affected. This suggests looking at the ratio of the variance of the $n$ th projected component to the total variance of the corresponding field contained in the Fourier modes at each $\kappa _h$ :

(4.7) \begin{align} R_n^f(\kappa _h) &= \bigg \langle \frac {\int _z \langle |\tilde {f}_{n}(\alpha ,\beta ,z,t)|^2 \rangle _{\phi _k} {\rm d}z}{\int _z \langle |\check {f}(\alpha ,\beta ,z,t)|^2 \rangle _{\phi _k} {\rm d}z} \bigg \rangle _t. \\[-12pt]\nonumber \end{align}

Where $f$ generically represents any of the five field variables and $\tilde {f}_{n}$ is the $n$ th term in the eigenfunction expansion (3.18). This quantity is visualized for $w$ and $T$ for each mode–wavenumber pair in figure 9. (The ratios for $u$ , $v$ and $p$ are similar to, or smaller than those for $w$ .) The high-energy modes are correlated with low errors suggesting that the projection procedure is robust to the particular subset of basis functions included. Also striking about figure 9 is that the ratio of projected temperature variance to the total variance is typically higher than the vertical kinetic energy ratio, indicating that the sloshing motions contain a larger portion of the potential energy than the kinetic.

4.2. Determining the set of basis functions

Based on the spectral peaks seen in the wavenumber–frequency spectra in figure 5 it is clear that only the first two modes should be included for SM126 and the first mode for SM125. A clear cutoff in terms of the wavenumber is less obvious but based on the error metrics and energy ratios reported in the previous section, a reasonable choice is the first 20 wavenumbers. This choice is somewhat arbitrary, but the results are not sensitive beyond approximately 10 wavenumbers as shown in Appendix E of Hass (Reference Hass2025).

In light of these observations, we propose the values in table 4 as a reasonable choice for the low-Fr cases.

Table 4.

Projection parameter selection for low-Fr cases.

Table 4. Long description

A table with two rows and three columns comparing projection parameter selections for low-Fr cases. The columns are labeled Case, n max, and κ h,maxL x/2π. The first row under Case is labeled SM125 with values 1 and √40 under n max and κ h,maxL x/2π respectively. The second row under Case is labeled SM126 with values 2 and √40 under n max and κ h,maxL x/2π respectively.

Figure 9.

Variance ratios $R_n^w$ and $R_n^T$ defined by (4.7): (a) SM125 (Fr = 3.5); (b) SM126 (Fr = 1.25).

Figure 9. Long description

Panel A: Two heat maps showing energy ratios for different modes and wave numbers. The x-axis represents (κhLx/2π)² ranging from 1 to 146, and the y-axis represents modes ranging from 1 to 4. The color scale on the right indicates values from 0 to 1, with blue representing lower values and red representing higher values. The left heat map shows Rwn(kh) with higher values concentrated in the lower left corner, indicating a higher energy ratio for lower modes and wave numbers. The right heat map shows Rtn(kh) with a similar pattern but with slightly different distributions. Panel B: Two heat maps showing energy ratios for different modes and wave numbers. The x-axis represents (κhLx/2π)² ranging from 1 to 146, and the y-axis represents modes ranging from 1 to 4. The color scale on the right indicates values from 0 to 1, with blue representing lower values and red representing higher values. The left heat map shows Rwn(kh) with higher values concentrated in the lower left corner, indicating a higher energy ratio for lower modes and wave numbers. The right heat map shows Rtn(kh) with a similar pattern but with slightly different distributions.

4.3. Moderate-Fr regime

Because SM123 and SM124 do not have adequate time records to resolve the low- $\kappa _h$ eigenmodes, we cannot analyse their dynamics with the same rigour as the low-Fr cases. Furthermore, these cases do not satisfy the $\textit{Fr} \ll 1$ limit implied by the linearized equations, and so the linear model is unlikely capable of representing the interfacial dynamics. We confirmed this by computing the error metrics for each case and saw that they are substantial for all modes (see Appendix E in Hass Reference Hass2025).

In spite of the poor detailed correspondence with the linear model, we maintain that a wave–turbulence decomposition is still desirable for two reasons. One, the inhomogeneous longitudinal correlation $c_{33,3}$ still shows a prominent plateau due to organized interfacial sloshing motions as seen in figure 10. Two, when we turn the turbulence forcing off, the interface remains energetically active for a long time despite negligible mixing. In this scenario, however, the nonlinear nature of the waves is visually evident in movies of temperature contours and we observe wave breaking and localized shear instabilities leading to small scale turbulence production (not shown here).

Because of these observations, we used the first nine wavenumbers of the first baroclinic mode as our projection basis to filter out wave-like motions. This choice is conservative, but somewhat arbitrary. See § 4.4.3 in Hass (Reference Hass2025) for more details regarding this selection.

5.1. Qualitative demonstration

Before analysing the statistics, we begin with a qualitative demonstration by looking at instantaneous contours $w$ , $p_{\!f}$ and $T_{\!f}$ . These are shown in figure 11 for the $\textit{Fr} = 1.25$ case. It is visually evident that the protocol is able to extract organized sloshing motions.

Figure 11.

Wave–turbulence decomposition for the $\textit{Fr}=1.25$ case. Columns (i) – (iii) show instantaneous contours of the $w$ , $p$ and $T$ fluctuations, respectively. Panel (a) is the raw LES data. Panel (b) shows the result of the eigenprojection, and (c) is the difference of the two, i.e. turbulence.

5.2. Fluctuating kinetic energy

The total kinetic energy can be written in terms of a wave contribution, turbulence and a cross-correlation: $k = \tilde {k} + k' + 2 \tilde {k}'$ where $\tilde {k} \equiv \langle \tilde {u}_i\tilde {u}_i \rangle /2$ , $k' \equiv \langle u_i'u_i' \rangle /2$ and $\tilde {k}' \equiv \langle u_i' \tilde {u}_i \rangle /2$ . The ratios $k'/k$ and $\tilde {k}/k$ are shown in figures 12(a) and 12(b). The kinetic energy ratios in figure 12(a) are plotted against depth whereas in figure 12(b) depth is scaled by $h_1$ resulting in a significant collapse of the profiles, demonstrating that buoyancy effects become significant for $\textit{Fr}_t \gtrsim 1$ . The collapsed profiles indicate that we can parameterize $k'$ in terms of $k$ since $h_1$ can reliably be determined by using the unfiltered fluctuations in the definition of $\textit{Fr}_t$ (not shown here, but demonstrated in Hass (Reference Hass2025)).

In the lowest- $\textit{Fr}$ case, SM126, nearly 60 % of the kinetic energy just above the interface is due to interfacial waves emphasizing the necessity of handling these in a systematic way. While the ratios of the other cases may seem unexpectedly high, the high ratios occur below the interface midline (where $\langle T \rangle = 0.5$ , marked by circles in figure 12) where we expect the wave contribution to dominate.

Figure 12.

Ratios of turbulent kinetic energy (TKE) $k'$ and wave kinetic energy (WKE) $\tilde{k}$ to the fluctuating kinetic energy $k$ . Open circles mark where $\langle T \rangle =0.5$ and closed circles (panel a only) mark where $\textit{Fr}_t = 1$ . (a) Kinetic energy ratios versus depth. (b) Kinetic energy ratios versus scaled depth.

Figure 12. Long description

The image contains four line graphs. Panel A shows kinetic energy ratios versus depth, while Panel B shows kinetic energy ratios versus scaled depth. Each graph has two y-axes: the left y-axis represents the ratio of TKE to fluctuating kinetic energy, and the right y-axis represents the ratio of wave kinetic energy to fluctuating kinetic energy. The x-axis in Panel A represents depth, and the x-axis in Panel B represents scaled depth. The graphs include data series for different Froude numbers (Fr = 14, Fr = 7, Fr = 3.5, Fr = 1.25), each represented by different colored lines. Open circles mark where the ratio of TKE to fluctuating kinetic energy is 0.5, and closed circles (in Panel A only) mark where the ratio of wave kinetic energy to fluctuating kinetic energy is 0.5. The graphs show how these ratios vary with depth and scaled depth for different Froude numbers.

5.3. Large-scale anisotropy and velocity variance

The ratio of the vertical to horizontal velocity root-mean-square (RMS), known as the isotropy factor $I$ , is a measure of the large-scale anisotropy of a flow where $I=1$ indicates isotropy. This quantity is plotted in figure 13. Dashed lines in figure 13(b) are computed using the full fluctuating fields and solid lines are computed after projecting out the eigenmodes. Overlayed on the plot of $I$ is the unstratified data which sets the baseline for our simulations which is $1.125 \lesssim I \lesssim 1.2$ . Clearly, the effect of the density interface is to suppress vertical velocity relative to horizontal as evidenced by the reduction in $I$ with increasing stratification. However, this effect is significantly more pronounced with the projected data since the sloshing interface acts to increase $w$ relative to its turbulent value (see figure 4.18 in Hass (Reference Hass2025)).

Figure 13.

Isotropy factor, $I = w_{rms}/u_{h,rms}$ . (a) Unfiltered fluctuations. (b) Solid lines are after sloshing modes have been removed and dashed lines are from the unfiltered data (same as solid lines in a). Open circles mark where $\langle T \rangle =0.5$ . Closed circles show the location where $\textit{Fr}_t = 1$ .

5.4. Two-point correlations

Previously, we have seen that the inhomogeneous longitudinal correlation $c_{33,3}$ is sensitive to interfacial sloshing. In this section, we re-evaluate this correlation and one other correlation relevant for length scale estimation namely the longitudinal autocorrelation coefficients defined as

(5.1) \begin{equation} f(r,z) = \frac {1}{2}\big ( c_{11,1}(r,z) + c_{22,2}(r,z) \big ). \end{equation}

Figure 14 shows $c_{33,3}$ at a number of observer locations before (dashed lines) and after (solid lines) the sloshing modes are removed. It is interesting that for nearly all the cases, the large plateau was replaced by a persistent, albeit smaller, negative correlation. The existence of a non-zero plateau suggests that the projection procedure did not completely separate waves and turbulence which is not surprising considering a perfect decomposition is likely unattainable. However, the much improved decorrelation (over the full fluctuating field) makes the prospect of computing a vertical length scale a viable option. Prior to the eigenprojection, we could not estimate a meaningful vertical length scale due to the large plateau.

Figure 14.

Inhomogeneous longitudinal autocorrelation coefficient $c_{33,3}$ at various $z$ -locations. Line colours correspond to the legend in figure 15. Dashed lines are computed from the full fluctuating variables and solid lines are after the eigenmodes have been filtered out.

Figure 14. Long description

The image contains four line graphs side by side, each depicting the inhomogeneous longitudinal autocorrelation coefficient at various locations. The x-axis is labeled ’r’ and the y-axis is labeled ’c333 (r, z)’. Each graph has a range from −2 to 2 on the x-axis and 0 to 1 on the y-axis. The graphs are titled with different Froude numbers: Fr = 14, Fr = 7, Fr = 3.5, and Fr = 1.25. Line colors correspond to a legend, with dashed lines representing full fluctuating variables and solid lines representing data after eigenmodes have been filtered out. Each graph shows a peak at r = 0, with varying patterns and amplitudes for different Froude numbers. The trends indicate how the autocorrelation coefficient changes with different Froude numbers and filtering methods.

Figure 15.

Longitudinal autocorrelation coefficient $f$ . Dashed lines are computed from the full fluctuating variables and solid lines are after the eigenmodes have been filtered out.

A horizontal length scale can be computed from the longitudinal autocorrelation coefficient defined in (5.1). This is plotted in figure 15 for each level of stratification and multiple $z$ -locations. Fortunately, from a turbulence modelling and diagnostic perspective, the horizontal correlation is robust in the presence of interfacial sloshing as evidenced by the close correspondence between solid and dashed curves. The one exception is very close to the density interface in the strongest stratified case (SM126, $\textit{Fr} = 1.25$ ). However, at this location the turbulent Froude number and buoyancy Reynolds numbers are 0.01 and 1.74, respectively, indicating that the flow is likely laminar here and so attempting to define a ‘turbulence’ length scale is erroneous. We therefore conclude that sloshing motions have negligible impact on horizontal two-point correlations, and thereby the horizontal integral length scales.

7.1. Budget equations

Contracting the equation for $u_i$ (2.2b ) with $u_i'$ and averaging results in an equation for the TKE $k'$ . Similarly, contracting (2.2b ) with $\tilde {u}_i$ and averaging yields an equation for WKE $\tilde {k}$ . These are given symbolically as

(7.1)

(7.2)

The fluxes are defined as

(7.3) \begin{align} F'_w &= -\langle wu_i'u_i' \rangle /2, & F'_{\!p} &= -\langle p'w' \rangle , & F'_{\nu } &= \frac {2}{\textit{Re}}\langle u_i' s_{i3}' \rangle - \langle u_i' \tau _{i3} \rangle, \nonumber \\[3pt] \tilde {F}'_{w1} &= -\langle u_i' \tilde {u}_i \tilde {w} \rangle , & \tilde {F}'_{\!p1} &= -\langle \tilde {p}w' \rangle , & \tilde {F}'_{\nu 1} &= \frac {2}{\textit{Re}}\langle u_i' \tilde {s}_{i3} \rangle, \nonumber \\[3pt] \tilde {F}_w &= -\langle w\tilde {u}_i \tilde {u}_i \rangle /2, & \tilde {F}_{\!p} &= -\langle \tilde {p} \tilde {w} \rangle , & \tilde {F}_{\nu } &= \frac {2}{\textit{Re}} \langle {\tilde {u}}_i \tilde {s}_{i3} \rangle, \nonumber \\[3pt] \tilde {F}'_{w2} &= -\langle \tilde {u}_i u_i' w' \rangle , & \tilde {F}'_{\!p2} &= -\langle p' \tilde {w} \rangle , & \tilde {F}'_{\nu 2} &= \frac {2}{\textit{Re}}\langle {\tilde {u}}_i s_{i3}' \rangle - \langle {\tilde {u}}_i \tau _{i3} \rangle. \end{align}

Below, we also employ the symbol $\mathcal{T}$ to denote the $z$ derivative of each flux, i.e. $\mathcal{T}_w' = \partial _z F'_w$ is the turbulent transport of TKE. The buoyancy fluxes and dissipation rates are given by

(7.4) \begin{align} \mathcal{B}' &= -\langle w'b' \rangle , & \tilde {\mathcal{B}}'_1 &= -\langle w'\tilde {b} \rangle, \nonumber \\[3pt] \tilde {\mathcal{B}} &= -\langle {\tilde {w}} \tilde {b} \rangle , & \tilde {\mathcal{B}}'_2 &= -\langle {\tilde {w}} b' \rangle ,\nonumber \\[3pt] \epsilon ' &= \frac {2}{\textit{Re}}\langle s_{ij}'s_{ij}' \rangle - \langle \tau _{ij}s_{ij}' \rangle , & \tilde {\epsilon }'_1 &= \frac {2}{\textit{Re}}\langle s'_{ij} \tilde {s}_{ij} \rangle, \nonumber \\[3pt] \tilde {\epsilon } &= \frac {2}{\textit{Re}} \langle \tilde {s}_{ij} \tilde {s}_{ij} \rangle , & \tilde {\epsilon }'_2 &= \frac {2}{\textit{Re}} \langle s_{ij}'\tilde {s}_{ij} \rangle - \langle \tau _{ij} \tilde {s}_{ij} \rangle .\end{align}

Finally, the unsteady cross-correlations are

(7.5) \begin{align} \tilde {\mathcal{U}}'_{k1} &= \bigg \langle u_i' \frac {\partial \tilde {u}_i}{\partial t} \bigg \rangle , \;\;\; \tilde {\mathcal{U}}'_{k2} = \bigg \langle {\tilde {u}}_i \frac {\partial u_i'}{\partial t} \bigg \rangle. \end{align}

In (7.1) and (7.2), cross-correlations between turbulence and wave quantities are shown in blue to highlight the difference between these and the traditional TKE equation. In addition to cross-correlation terms, there is a term $\mathcal{I}$ which shows up in both equations but with opposite sign, and therefore is responsible for the transfer of kinetic energy between waves and turbulence:

(7.6) \begin{align} \mathcal{I} &= \langle {\tilde {u}}_i {\tilde {u}}_{\!j} s_{ij}' \rangle -\langle u_i'u_{\!j}'\tilde {s}_{ij} \rangle. \end{align}

The form of the term is directly analogous to the interscale transfer term in Wyngaard’s thought experiment during his introduction of LES (Wyngaard Reference Wyngaard2010). It is the difference between the turbulent stress working against wave strain and the wave stress working against turbulent strain. Only when these two processes are unequal can there be a transfer of energy from one component to the other.

For maximum transparency, we also decompose the kinetic energies $k'$ and $\tilde {k}$ into horizontal and vertical components: $k' = 2k_h' + k_w'$ and $\tilde {k} = 2\tilde {k}_h + \tilde {k}_w$ . These are defined as

(7.7) \begin{align} k_h' &= \langle u^{\prime 2} + v^{\prime 2} \rangle /4, & k_w' &= \langle w^{\prime 2} \rangle /2 ,\nonumber \\[3pt] \tilde {k}_h &= \langle {\tilde {u}}^2 + {\tilde {v}}^2 \rangle /4, & \tilde {k}_w &= \langle {\tilde {w}}^{\prime 2} \rangle /2.\end{align}

The budget equations for horizontal kinetic energy are

(7.8)

(7.9)

and vertical kinetic energy

(7.10)

(7.11)

The $\phi$ -terms are pressure-strain correlations responsible for kinetic energy redistribution between vertical and horizontal motions. The remaining terms are componentwise analogues to the TKE budget terms such as dissipation ( $\epsilon _w$ and $\epsilon _h$ ) and fluxes ( $F$ -terms). In the interest of space, all terms are defined in Appendix A. Due to their importance in the following analysis, we define the pressure-strain correlations and wave–turbulence transfer terms here:

(7.12) \begin{align} \phi '_h &= \langle p'( s_{11}' + s_{22}') \rangle /2 ,\nonumber \\ \tilde {\phi }_h &= \langle {\tilde {p}} (\tilde {s}_{11} + \tilde {s}_{22}) \rangle /2, \nonumber \\ \phi _w' &= \langle p' s_{33}' \rangle, \nonumber \\ \tilde {\phi }_w &= \langle {\tilde {p}} \tilde {s}_{33} \rangle ,\nonumber \\ \mathcal{I}_h &= \langle {\tilde {u}} {\tilde {u}}_{\!j} \partial _{\!j} u' + {\tilde {v}} {\tilde {u}}_{\!j} \partial _{\!j} v' \rangle /2 - \langle u' u_{\!j}' \partial _{\!j} {\tilde {u}} + v' u_{\!j}' \partial _{\!j} {\tilde {v}} \rangle /2, & & \nonumber \\ \mathcal{I}_w &= \langle {\tilde {w}} {\tilde {u}}_{\!j} \partial _{\!j} w' \rangle - \langle w' u_{\!j}' \partial _{\!j} {\tilde {w}} \rangle. & & \end{align}

Due to incompressibility $\phi ' \equiv \phi _h' = -\phi'_{w}/2$ and $\tilde {\phi } \equiv \tilde {\phi }_h = -\tilde {\phi }_{w}/2$ .

The fluctuating potential energy is defined as

(7.13) \begin{align} k_{\!p} = \frac {\langle T_{\!f} T_{\!f} \rangle }{2\textit{Fr}^4 N^2} \end{align}

which can further be decomposed via the wave–turbulence decomposition: $k_{\!p} = k_{\!p}' + \tilde {k}_{\!p} + 2\tilde {k}_{\!p}'$ where $k_{\!p}' = \langle T'T' \rangle / 2\textit{Fr}^4 N^2$ is the turbulent potential energy (TPE), $\tilde {k}_{\!p} = \langle \tilde {T} \tilde {T} \rangle /2\textit{Fr}^4N^2$ is the wave potential energy (WPE) and $\tilde {k}'_{\!p} = \langle \tilde {T} T' \rangle /2\textit{Fr}^4N^2$ is the cross-correlation. The potential energy budget equations are derived by multiplying (2.2c ) by $T'/\textit{Fr}^4$ (for the $k_{\!p}'$ equation) and $\tilde {T}/\textit{Fr}^4$ (for the $\tilde {k}_{\!p}$ equation), averaging, and dividing by $N^2$ (consistent with the quasistationary approximation used in computing averages, $N^2$ can be brought into the time-derivative on the left-hand side):

(7.14)

(7.15)

Where the buoyancy variance flux $F_b$ and potential energy destruction rates $\chi$ are defined as

(7.16) \begin{align} F_b' &= -\langle b' b' w \rangle /2 + \frac {1}{\textit{Re} {\textit{Pr}}} \bigg \langle b' \frac {\partial b'}{\partial z} \bigg \rangle -\frac{1}{\textit{Fr}^2} \langle b' q_3' \rangle & \tilde {F}_{b2}' &= - \langle b' {\tilde {b}} {\tilde {w}} \rangle + \frac {1}{\textit{Re} {\textit{Pr}}}\bigg \langle b' \frac {\partial {\tilde {b}}}{\partial z} \bigg \rangle, \nonumber \\[5pt] \tilde {F}_{b} &= - \langle {\tilde {b}} {\tilde {b}} w \rangle /2 + \frac {1}{\textit{Re} {\textit{Pr}}} \bigg \langle {\tilde {b}} \frac {\partial {\tilde {b}}}{\partial z} \bigg \rangle & \tilde {F}'_{b1} &= - \langle {\tilde {b}} b' w' \rangle + \frac {1}{\textit{Re} {\textit{Pr}}}\bigg \langle {\tilde {b}} \frac {\partial b'}{\partial z} \bigg \rangle -\frac{1}{\textit{Fr}^2} \langle {\tilde {b}} q_3' \rangle, \nonumber \\[5pt] \chi ' &= \frac {1}{N^2} \bigg [\frac {1}{\textit{Re} {\textit{Pr}}}\bigg \langle \frac {\partial b'}{\partial x_{\!j}}\frac {\partial b'}{\partial x_{\!j}} \bigg \rangle - \frac {1}{\textit{Fr}^2}\bigg \langle q_{\!j}' \frac {\partial b'}{\partial x_{\!j}} \bigg \rangle \bigg ] & \tilde {\chi }'_2 &= \frac {1}{N^2}\frac {1}{\textit{Re} {\textit{Pr}}}\bigg \langle \frac {\partial b'}{\partial x_{\!j}}\frac {\partial {\tilde {b}}}{\partial x_{\!j}} \bigg \rangle, \nonumber \\[5pt] \tilde {\chi } &= \frac {1}{N^2}\frac {1}{\textit{Re} {\textit{Pr}}}\bigg \langle \frac {\partial {\tilde {b}}}{\partial x_{\!j}}\frac {\partial {\tilde {b}}}{\partial x_{\!j}} \bigg \rangle & \tilde {\chi }'_1 &= \frac {1}{N^2} \bigg [\frac {1}{\textit{Re} {\textit{Pr}}}\bigg \langle \frac {\partial {\tilde {b}}}{\partial x_{\!j}}\frac {\partial b'}{\partial x_{\!j}} \bigg \rangle - \frac {1}{\textit{Fr}_2}\bigg \langle q_{\!j}' \frac {\partial {\tilde {b}}}{\partial x_{\!j}} \bigg \rangle \bigg ] ,\end{align}

where $b\equiv T/\textit{Fr}^2$ is the non-dimensional buoyancy. The wave–turbulence transfer of potential energy is quantified by $\mathcal{I}_{k_{\!p}}$ :

(7.17) \begin{align} \mathcal{I}_{k_{\!p}} &= \frac {1}{N^2}\big [\langle {\tilde {b}} {\tilde {u}}_{\!j}\partial _{\!j} b' \rangle - \langle b' u_{\!j}' \partial _{\!j} {\tilde {b}} \rangle \big ] .\end{align}

Figure 21.

Kinetic energy budgets of the turbulence $k'$ and wave components $\tilde {k}$ : (a) SM123 ( $\textit{Fr} = 14$ ); (b) SM126 ( $\textit{Fr} = 1.25$ ). Dashed horizontal grey lines depict the $\langle T \rangle =0.5$ location and the dash–dot line is where $\textit{Fr}_t = 1$ . All terms are normalized by the total transport of $k'$ . All cross-correlations are shown in dashed lines.

7.2. Kinetic energy

The non-negligible terms in (7.1) are plotted in figure 21 where each has been normalized by the total TKE transport ( $\mathcal{T}' = \mathcal{T}'_w + \mathcal{T}_{\!p}' + \mathcal{T}_{\nu }'$ ). (In every case and for each budget viscous transport is negligible.) Also overlayed in each figure is a grey dashed line corresponding to the $\langle T \rangle =0.5$ location and a grey dash–dot line marking where $\textit{Fr}_t = 1$ , approximately where buoyancy effects first begin to become important. Only budgets for the limiting cases SM123 and SM126 ( $\textit{Fr} = 14$ and 1.25) are shown in this section for brevity, but budgets for the other simulations are included in Appendix A.

Conveniently, it turns out that the cross-correlation terms (blue terms) in (7.1) are nearly in balance with each other and so we can simply analyse the ‘pure’ turbulence terms (this was verified, but not shown here). We do include the cross-correlations as thin dashed lines to demonstrate that, in general, wave–turbulence correlations are non-zero.

There is a balance between total transport and dissipation above the $\textit{Fr}_t=1$ line for each case where $\mathcal{T}_{\!p}' \simeq -\mathcal{T}_w'/3$ , see figure 21. Just below the $\textit{Fr}_t=1$ line, the buoyancy flux begins to be a non-negligible contribution to the overall budget. As we move farther down, the relative contribution of turbulent and pressure transport begins to change noticeably. At this location the covariance terms $\tilde {\mathcal{U}}'$ and $\tilde {\mathcal{B}}'$ also become significant, but again, are in balance with each other.

The WKE budgets have negligibly small terms above the $\textit{Fr}_t=1$ line. Once the budget terms become active, the two wave-pure terms that are substantial are the buoyancy flux $-\langle {\tilde {w}} {\tilde {b}} \rangle$ and pressure transport $\partial _z\langle {\tilde {w}} {\tilde {p}} \rangle$ . For waves governed by linear dynamics $\langle {\tilde {w}} \tilde {b} \rangle$ is identically zero since $\tilde {w}$ and $\tilde {b}$ are $90^{\circ }$ out of phase. However, for the wave-breaking mechanism to be active, nonlinear interactions must be important and so a non-zero buoyancy flux is expected.

It is interesting to note that the unsteady terms, $\partial _t k'$ and $\partial _t \tilde {k}$ , are required to balance the budget for SM123 in the interfacial region (figure 21 a) indicating a breakdown in the quasistationary assumption at this weakest level of stratification.

7.2.1. Horizontal and vertical kinetic energy

The horizontal and vertical kinetic energy budgets for SM123 and SM126 ( $\textit{Fr} = 14$ and 1.25) are shown in figures 22 and 23 with budgets for the remaining cases in Appendix A. Dissipation and transport are in balance above the $\textit{Fr}_t=1$ line for both vertical and horizontal TKE. Below this line, buoyancy begins to tax $k_w'$ at the expense of spatial transport. For high stratification strengths, pressure transport changes sign near the interface becoming a source of TKE similar to the near wall behaviour of wall-bounded shear flows (Pope Reference Pope2001). The $k_h'$ reservoir also taxes $k_w'$ in this region via the pressure-strain correlation which is most evident in figure 23.

Looking at the wave budgets $\tilde {k}_h$ and $\tilde {k}_w$ below the $\textit{Fr}_t=1$ location, energy flows from horizontal motions to vertical via the pressure-strain correlation $\tilde {\phi }$ . Furthermore, the wave–turbulence transfer term is a significant source of energy for horizontal WKE in the interfacial region. $\tilde {k}_w$ is significantly taxed by the wave buoyancy flux, feeding the WPE reservoir as seen in figure 24 in the next section.

Figure 22.

Kinetic energy budgets for SM123 ( $\textit{Fr} = 14$ ). (a) Horizontal kinetic energy. (b) Vertical kinetic energy. Horizontal dash–dot and dashed lines correspond to the location where $\textit{Fr}_t = 1$ and $\langle T \rangle =0.5$ , respectively.

Figure 22. Long description

The image contains four line graphs depicting kinetic energy budgets. Panel A shows the horizontal kinetic energy budget on the left and the modified horizontal kinetic energy budget on the right. Panel B shows the vertical kinetic energy budget on the left and the modified vertical kinetic energy budget on the right. Each graph has the vertical axis labeled ’z’ with units ranging from −2.5 to −1.0. The horizontal axis for the left graphs in both panels ranges from −1 to 1, while the horizontal axis for the right graphs ranges from −0.2 to 0.1 in Panel A and from −0.5 to 0.5 in Panel B. The graphs include multiple lines representing different components of the kinetic energy budget, as indicated by the legend on the right. The legend includes labels such as ∂t k’, ε’, T’w, T’p, −B’Ŭ, I, Ŭk1, Ŭk2, T’w1, T’w2, T’p1, T’p2, −B1, −B2, and Residual, each represented by different colors and line styles. The graphs show various trends and interactions between these components across different depths (z values).

Figure 23.

Kinetic energy budgets for SM126 ( $\textit{Fr} = 1.25$ ). (a) Horizontal kinetic energy. (b) Vertical kinetic energy. Horizontal dash–dot and dashed lines correspond to the location where $\textit{Fr}_t = 1$ and $\langle T \rangle =0.5$ , respectively. Line colours and styles correspond to the legends in figure 22.

7.3. Potential energy and the mixing coefficient $\varGamma$

The $k_{\!p}'$ budgets show a dominant balance between buoyancy flux and destruction (i.e. $\mathcal{B}\simeq \chi$ ) above the $\langle T \rangle =0.5$ line. This simple expression underlies the arguments of Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) and is exact in triply periodic statistically homogeneous and stationary simulations of stratified turbulence. The fact that we have an inherently inhomogeneous flow and yet this balance holds suggests that the scaling relationships from homogeneous stratified shear-free DNS may carry over to our problem, e.g. the $\varGamma \sim \textit{Fr}_t^{-n}$ relationships reported in Maffioli et al. (Reference Maffioli, Brethouwer and Lindborg2016), Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) and Yi & Koseff (Reference Yi and Koseff2022). To test this, we plot a locally defined mixing coefficient $\varGamma (z)=\chi (z)/\epsilon (z)$ versus the turbulent Froude number $\textit{Fr}_t(z) = \epsilon (z)/N(z)k'(z)$ in figure 25. (Note that $\chi \simeq \chi '$ and $\epsilon \simeq \epsilon '$ and so the total destruction rates can be considered.) The functional dependence of each variable on depth $z$ is explicitly written to contrast this coefficient from the volume-averaged $\varGamma$ typically used to parameterize oceanographic measurements.

Figure 24.

Potential energy budgets of the turbulence $k_{\!p}'$ and wave components $\tilde {k}_{\!p}$ : (a) SM123 ( $\textit{Fr} = 14$ ); (b) SM126 ( $\textit{Fr} = 1.25$ ). Dashed horizontal grey lines depict the $\langle T \rangle =0.5$ location and the dash–dot line is where $\textit{Fr}_t = 1$ . All cross-correlations are shown in dashed lines.

Figure 24. Long description

The image contains four line graphs depicting potential energy budgets of turbulence and wave components. The graphs are divided into two rows, each with two panels. Panel A and Panel B in each row show different budgets. Panel A in the top row shows the k’p budget and kp budget, while Panel B in the bottom row shows the k’p budget and kp budget. Each graph has a vertical axis labeled ’z’ and a horizontal axis with different units. The graphs include multiple lines representing different components of the budgets, such as −χ’, χ’, Tw,b/N^2, and B. The lines are color-coded and labeled accordingly. The graphs show various trends and patterns, with lines intersecting and diverging at different points. The dashed horizontal grey lines depict specific locations, and the dashdot line indicates where certain conditions are met. All cross-correlations are shown in dashed lines.

Figure 25.

Mixing coefficient $\varGamma = \chi /\epsilon$ . Symbol colours denote $\textit{Fr} = 14$ , 7, 3.5, and 1.25 in red, yellow, purple, and green, respectively.

It is striking that this inhomogeneous flow possesses two regimes $\textit{Fr}^{-2}$ and $\textit{Fr}_t^0$ for $\textit{Fr}_t \gg 1$ and $\textit{Fr}_t \ll 1$ consistent with the scaling arguments and DNS results of Maffioli & Davidson (Reference Maffioli and Davidson2016), Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019) and Yi & Koseff (Reference Yi and Koseff2022) for homogeneous stratified turbulence. Issaev et al. (Reference Issaev, Williamson and Armfield2023) also computed a local mixing coefficient in an inhomogeneous flow (radiatively heated half-channel) and found similar power laws relating $\varGamma$ and $\textit{Fr}_t$ when intermittency was not accounted for. When averaging was restricted to turbulent-only regions the $\textit{Fr}_t^{-1}$ scaling persisted to the smallest values of $\textit{Fr}_t$ suggesting that the low- $\textit{Fr}_t$ plateau is a signature of intermittency. Whether such a conclusion is supported in our simulations would be interesting to investigate in the future.

Figure 26.

Energy conversion pathways. Arrow colours correspond to the line colours in the budget figures.

In the present case, each point in figure 25 can be interpreted as representing a small horizontal control volume where the influx of $k_{\!p}'$ is roughly balanced by its efflux resulting in no net transport.

The data plotted in figure 25 is restricted to locations above the interface midline (i.e. where $\langle T \rangle =0.5$ ). Below this location, there are rare events where an energetic turbulent eddy penetrates the interface carrying warm fluid into the cold quiescent layer below. In this scenario the buoyancy gradients become larger than the momentum gradients and $\chi \gt \epsilon$ within the eddy. When averaged over the plane, these mixing events result in large values of $\varGamma$ ( $1 \leqslant \varGamma \lt 2$ ) and result in a departure from the scaling seen in figure 25. One such turbulent eddy is visualized in figure 32 for case SM123 ( $\textit{Fr}=14$ ).

The key to connecting our data to homogeneous turbulence simulations is the approximate balance $\mathcal{B}\simeq \chi$ . Other inhomogeneous flows may not possess such a balance and in those cases it seems unlikely that local mixing coefficients would follow the scaling reported here and in previous studies.

7.4. Wave–turbulence coupling and the flow of energy

Based on the observations noted in §§ 7.2 and 7.3, and a close inspection of the six budget equations for the remaining cases in Appendix A, we propose the energy flow diagram in figure 26 as a general description of the interfacial region (i.e. in the neighbourhood of $\langle T \rangle = 0.5$ ). The diagram does not exactly describe any one simulation, but the general path is a reliable approximation for the observed budgets in this region. In particular, cross-correlations do not internally balance for the wave budgets and can be a source or sink. For example, $\tilde {\phi }'_{h2} = \langle p'( \tilde {s}_{11} + \tilde {s}_{22}) \rangle /2$ is a significant source for $\tilde {k}_h$ in the $\textit{Fr} = 1.25$ case (SM126) near the density interface, but not so for the other cases. Clearly, the detailed picture is more complicated than the diagram in figure 26, yet the figure is a general description of the processes involved.

Kinetic energy is injected into the domain via the forcing term. Beneath the forcing layer, this corresponds to feeding the $k_h'$ and $k_w'$ reservoirs through turbulent transport $\mathcal{T}_{w,h}'$ and $\mathcal{T}_{w,w}'$ . Vertical kinetic energy is converted to horizontal and potential energy (via $\phi _w'$ and $\mathcal{B}'$ , respectively) near the density interface as well as dissipated throughout the water column via $\epsilon _w'$ . Horizontal kinetic energy that is not dissipated drives wave motions through $\mathcal{I}_h$ which in turn produces vertical motions by the pressure-strain correlation $\tilde {\phi }$ . Here $\tilde {k}_w$ is fed via the pressure field by the pressure-strain redistribution from $\tilde {k}_h$ as well as the divergence of the wave pressure flux $\tilde {T}_{\!p,w} = -\partial _z\langle {\tilde {w}} {\tilde {p}} \rangle$ . The WPE is primarily driven by the wave buoyancy flux $\tilde {\mathcal{B}}$ which is then transferred to TPE via the transfer term $\mathcal{I}_{k_{\!p}}$ and ultimately destroyed by molecular diffusion $\chi '$ .

To validate the energy flow diagram in figure 26 we plot the integral of each term (arrow in the diagram) normalized by the integral of all sources for each reservoir (e.g. $k_h'$ , $\tilde {k}_h$ , etc.). Specifically, we define the operator $\mathcal{R}[\boldsymbol{\cdot }]$ acting on the budget terms as

(7.18) \begin{align} \mathcal{R} \big[T_i^{({E_{\textit{bdg}}})}\big] &\equiv \dfrac {\int _{z_{co}}^{z_1} T_i^{({E_{\textit{bdg}}})} {\rm d}z}{\int _{z_{co}}^{z_1} \dfrac {1}{2} \sum \limits _{i} |T_i^{({E_{\textit{bdg}}})}| {\rm d}z} \end{align}

where ${E_{\textit{bdg}}} \in \{k_h', k_w', k_{\!p}', \tilde {k}_h, \tilde {k}_w, \tilde {k}_{\!p}\}$ labels the energy budget and the $T_i$ are individual terms in the respective equation (i.e. (7.8)–(7.11) and (7.14)–(7.15)). Here $z_{co}$ is the cutoff coordinate based on sufficient time averaging (see § 2.1 and Hass & Lele (Reference Hass and Lele2026)) and $z_1$ is defined as $\textit{Fr}_t(z_1) \equiv 1$ which can be visualized as a dash–dot line in the budget figures above (and in Appendix A). This is where buoyancy effects begin to influence the large scales. One half the sum of the magnitudes of $T_i$ is a proxy for the total of all sources of energy since the residuals (sum of sources and sinks) are nearly zero (black lines in budget figures). This normalization allows for an objective calculation of the amount of ‘source’ even if individual terms change sign.

These ratios are plotted in figure 27. The takeaway is that the qualitative flow described in figure 26 is supported by the data with one notable exception occurring for simulation SM126 ( $\textit{Fr}=1.25$ ) for the $\tilde {k}_h$ budget (lower left-hand corner of figure 27) where the wave–turbulence transfer term appears to be a sink for horizontal WKE. Close inspection of the horizontal WKE budget in figure 23(a) reveals that near the interface midline $\mathcal{I}_{h}$ changes sign and $\tilde {\phi }_{h2}' \equiv \langle p'(\tilde {s}_{11} + \tilde {s}_{22}) \rangle /2$ becomes the predominant source suggesting that there is still a transfer of energy from horizontal TKE to WKE. This term is difficult to interpret as it does not show up in the $k_h'$ budget and so is not a true ‘transfer’ term. Regardless, the other arrows in figure 26 are supported by the rest of the data in figure 27.

Figure 27.

Integrals of individual budget terms normalized by the total source-term contribution for each energy reservoir (7.18). Symbol colours correspond to arrow colours in figure 26. The horizontal axis is numbered by the simulation Froude number – see table 2.

In general, the flow of energy between turbulence and waves follows $k_w' \to k_h' \to \tilde {k}_h \to \tilde {k}_w \to \tilde {k}_{\!p} \to k_{\!p}'$ . The transfer from waves to turbulence by $\mathcal{I}_{k_{\!p}}$ generally happens below the $\langle T \rangle =0.5$ location. This is clearly seen by plotting the ratio $\mathcal{I}_{k_{\!p}}/\mathcal{B}'$ as a function of depth in figure 28 which quantifies the relative contribution of waves to the TPE reservoir compared with the TKE contribution. The ratio is significant in the neighbourhood $(z-z_{\!f})/H = 1$ and increases with increasing stratification. This is consistent with the experimental observations of Hannoun & List (Reference Hannoun and List1988) and McGrath et al. (Reference Mcgrath, Fernando and Hunt1997) where in the limit of strong stratification, entrainment is due almost entirely to internal wave breaking (i.e. a transfer from wave motions to turbulence). We provide further evidence of the wave breaking mechanism in the following section.

Figure 28.

Ratio of wave-to-turbulence transfer of potential energy to the buoyancy flux $\mathcal{I}_{k_{\!p}}/\mathcal{B}'$ . Depth is normalized by the mixed-layer depth $H$ , defined as the distance from the forcing layer edge $z_{\!f}$ and the location where $\langle T \rangle =0.5$ .

7.5. Wave breaking and entrainment

It has been experimentally demonstrated that the primary entrainment mechanism at high stratification is interfacial waves breaking. The statistical analysis of the previous section reinforces this picture as the ratio $\mathcal{I}_{k_{\!p}}/\mathcal{B}'$ is ${O}(10)$ for the highest level of stratification. The term $\mathcal{I}_{k_{\!p}}$ is a difference of correlations (7.17) and does not prove by itself that waves are breaking and so we present further evidence in this section. Figure 29 shows a sequence in time of temperature contours in the $x$ – $z$ plane to visualize a wave breaking event. Two waves are moving towards each other which is visualized by tracking two interfacial disturbances. As they collide, the wave steepens and appears to break generating small scale turbulence and localized mixing.

Figure 29.

Visualizing a wave breaking event on the interface. Instantaneous contours of temperature are shown in the $x$ – $z$ plane where blue-to-red corresponds to increasing temperature (bounded by 0 and 1). Snapshots proceed from top to bottom; the white boxes track two interfacial disturbances, and the arrows indicate their directions of travel. The time spanned by the snapshots is 0.16 eddy turnover times measured at the interface.

We expect the location of this mixing event to have significant magnitudes of either $\tilde {b} \tilde {u}_{\!j} \partial _{\!j} b'$ , $-b' u_{\!j}' \partial _{\!j} \tilde {b}$ or both. Figures 30 and 31 demonstrate this to be the case where the first figure visualizes pointwise values in the $x$ – $z$ plane and the latter in the $x$ – $y$ plane. Figure 31 demonstrates that at this time instant there are multiple such breaking events.

Figure 30.

Pointwise values of wave–turbulence terms contributing to $\mathcal{I}_{k_{\!p}}$ in the $x$ – $z$ plane for simulation SM126 ( $\textit{Fr}=1.25$ ). A vertical white or black line marks the location of the breaking event and a horizontal line marks the $z$ -plane visualized in figure 31.

Figure 30. Long description

Panel A: A heatmap shows temperature distribution in a fluid. The horizontal axis is labeled ’x’ with a range from 0 to 9. The vertical axis is labeled ’z’ with a range from −2.5 to −1.0. The color scale on the right ranges from 0 to 1, with blue indicating lower temperatures and red indicating higher temperatures. A vertical white line marks the location of a breaking event, and a horizontal line marks the plane visualized in another figure. Panel B: A heatmap displays the product of velocity and buoyancy fluctuation gradients. The horizontal axis is labeled ’x’ with a range from 0 to 9. The vertical axis is labeled ’z’ with a range from −2.5 to −1.0. The color scale on the right ranges from −1 to 1 times 10^−6, with blue indicating negative values and red indicating positive values. A vertical black line marks the location of the breaking event, and a horizontal line marks the plane visualized in another figure. Panel C: A heatmap shows the product of buoyancy and velocity fluctuation gradients. The horizontal axis is labeled ’x’ with a range from 0 to 9. The vertical axis is labeled ’z’ with a range from −2.5 to −1.0. The color scale on the right ranges from −1 to 1 times 10^−6, with blue indicating negative values and red indicating positive values. A vertical black line marks the location of the breaking event, and a horizontal line marks the plane visualized in another figure.

Figure 31.

Pointwise values of wave–turbulence terms contributing to $\mathcal{I}_{k_{\!p}}$ in the $x$ – $y$ plane for simulation SM126 ( $\textit{Fr}=1.25$ ). A vertical line marks the $x$ -location of the breaking event in figure 30 and the horizontal line marks the $y$ -location of the $x$ - $z$ plane visualized in figure 30.

The peak value of $\mathcal{I}_{k_{\!p}}/\mathcal{B}'$ for the weakest stratified case, SM123 is 0.97 (red curve in figure 28) suggesting that wave breaking may be in play for this case as well. Instantaneous $x$ – $z$ contours of $T$ , $\tilde {b}\tilde {u}_{\!j} \partial _{\!j} b'$ and $-b' u_{\!j}' \partial _{\!j} \tilde {b}$ are visualized in figure 32 for this simulation. Here, there is a large energetic turbulent eddy penetrating the interface engulfing quiescent cold fluid as it rebounds due to the stabilizing buoyancy force. This eddy registers large magnitudes of $\tilde {b} \tilde {u}_{\!j} \partial _{\!j} b'$ and yet is clearly not a wave-breaking event. In this case, the turbulent eddy has steep buoyancy gradients due to the relatively cold surrounding fluid. The large value of $|\boldsymbol{\nabla }b'|$ results in moderate values of $ \langle \tilde {b} \tilde {u}_{\!j} \partial _{\!j} b' \rangle$ so long as $\tilde {b} \tilde {u}_{\!j}$ is non-zero.

Figure 32.

Pointwise values of wave–turbulence terms contributing to $\mathcal{I}_{k_{\!p}}$ in the $x$ – $z$ plane for simulation SM123 ( $\textit{Fr}=14$ ). Vertical and horizontal white/black lines mark the location of a turbulent eddy penetrating the interface.

It is clear that while the decomposition proposed in this work provides a meaningful separation between wave and turbulent fluctuations, caution must be taken when interpreting statistics for cases that depart too far from the linear model used to derive the decomposition. It is fascinating that while the physical mechanism of entrainment changes from the $\textit{Fr}=14$ simulation to $\textit{Fr}=1.25$ , the measured entrainment rate continues to follow a simple power law with the Richardson number (Hass & Lele Reference Hass and Lele2026).