3.1. Time-averaged statistics

Turbulent statistics are presented in figure 2. In this paper the overbar represents the combination of planar ($x$ and $z$) and temporal averaging over time $T$ (table 1), and the primes represent the remaining fluctuating component of a variable: $\rho = \bar {\rho } + \rho '$. The same averaging process is used to calculate root-mean-square and flux variables, for example, $u_{rms}^{2} = \overline {u' u'}$. Plane Poiseuille flow is characterised by steep velocity gradients and high turbulent kinetic energy (TKE) near the wall, and a maximum velocity ${U}_{max}$ at the channel centre. For clarity, data presented in figure 2 are for $y=0\to 1$, given that all flow statistics are symmetric about the centreline apart from $\bar {\rho }$ and $\overline {u'v'}^{+}$, which are antisymmetric (i.e. symmetric with a change of sign). For the unstratified case (U), where density is transported as a passive scalar, the advective scalar (‘buoyancy’) flux ($\overline {\rho 'v'}$) is constant for $y \gtrsim 0.2$ and the momentum flux ($\overline {u'v'}$) and density field vary linearly with the vertical coordinate. Deviations between unstratified (U) and stratified (S) averaged statistics are most clear far from the walls ($y \gtrsim 0.2$).

Figure 2. Profiles of flow statistics, averaged in $x$, $z$ and $t$. Superscript $+$ represents wall-unit normalisation by $u_\tau$, $\delta _\nu = \nu /u_\tau$, and $\rho _\tau = -\Delta \rho \textit {Re}_\tau ^{-1} \textit {Pr}^{-1} \partial _y \bar {\rho }$. Solid lines are Case U data, dashed lines are Case S data. Areas shaded light and dark grey respectively highlight the inner ($y \leqslant 0.2$) and channel core ($y > 0.8$) regions, while the outer region ($0.2 < y \leqslant 0.8$) is unshaded. Areas shaded green represent regions where the phase speed of ‘linear’ internal waves is equal to the background planar and temporally averaged flow velocity: $c = \overline{U}_{mean} - N_{mean} / k_x = \overline{U}$. Here, $N$ is the buoyancy frequency, defined in (3.1a–d), and the subscript ‘mean’ denotes averaging over the channel core ($y > 0.8$); $c$ is calculated for a streamwise wavenumber ($k_x$) range of 2.5 to 5.0. These $k_x$ bounds correspond to 50 % of the maximum $v'$ spectral energy at $y=1$ (as shown in figure 8).

At high levels of stratification three distinct regions in the flow are observed, the inner, outer and core regions. The inner region, $y \leqslant 0.2$ and $y \geqslant 1.8$, comprises the viscosity-affected boundary layer, while the outer region, $0.2 < y \leqslant 0.8$ and $1.2 \leqslant y < 1.8$ comprises the turbulent region where the direct effects of viscosity on $\overline{U}$ are negligible. Finally, the channel core region, $0.8 < y < 1.2$, is where buoyancy dominates the dynamics. This classification is based upon the standard definitions of the inner and outer regions of an (unstratified) turbulent boundary layer (see e.g. Pope Reference Pope2000), with the channel core defining the region of high density gradient. While these labels are used throughout the manuscript it is important to note their qualitative nature; the quantitative bounds between these different regions are (naturally) less well defined in such a transient flow.

The inner region is largely unchanged by stratification when $\textit {Ri}_\tau$ is increased, owing to its high shear and production of turbulence (figure 2). The collapse of near-wall U and S data is expected due to the imposed constant pressure gradients and therefore friction (quantified through $\textit {Re}_\tau$). However, it should be noted that the bulk Reynolds number $\textit {Re}_b$ is significantly increased for Case S (see table 1). The (relatively) high turbulence production leads to a well-mixed region of essentially constant density extending through the outer region of the flow. Conversely, the channel core is characterised by a strong density gradient and low levels of turbulence, evidenced through profiles of $u^{+}_{rms}$ and $w^{+}_{rms}$. Note that $v^{+}_{rms}$ actually increases in this region due to the presence of (internal gravity) waves, as will be shown in spectra and DMD data. This is further illustrated by profiles of the buoyancy flux ($\overline {\rho ' v'}^{+}$) which tends to zero at the channel centre for Case S. This behaviour is explained by the (essentially perfect) decorrelation of $\rho '$ and $v'$ which are out of phase by ${\rm \pi} /2$. Such a phase difference is characteristic of the polarisation relations of linear internal waves, as noted by Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011).

Apart from the Prandtl number $\textit {Pr}$ (determined by molecular properties of the fluid), three other natural non-dimensional quantities for turbulent stratified shear flows are the gradient Richardson number, $\textit {Ri}_g$, the horizontal Froude number, $\textit {Fr}_h$, and the buoyancy Reynolds number, $\textit {Re}_B$:

(3.1a–d)\begin{equation} \textit{Ri}_g = \frac{N^{2}}{\left( \partial \overline{U}/\partial y\right)^{2}}, \quad \textit{Fr}_h = \frac{\bar{\varepsilon}}{N u_{{rms}}^{2}}, \quad \textit{Re}_B = \textit{Re}_\tau \frac{\bar{\varepsilon}}{N^{2}}, \quad N^{2} ={-}\textit{Ri}_\tau \frac{\textrm{d} \bar{\rho}}{\textrm{d} y}, \end{equation}

where $\bar {\varepsilon } = \textit {Re}_\tau ^{-1} \overline {\partial _i u_j'\partial _i u_j'}$ is the time- and planar-averaged dissipation rate of TKE and $N$ is the buoyancy frequency. These quantities capture the relative significance of the time scales of buoyancy, shear, turbulence and viscous dissipation, and are in general functions of $y$ (when constructed from planar and temporal averaging). Their profiles are presented in figure 3, although it should be noted that quantities dependent on the dissipation rate of TKE ($\textit {Fr}_h$ and $\textit {Re}_B$) are only approximate due to the filter-based artificial dissipation. Here, the defining boundaries of the channel core ($0.8 < y < 1.2$) are more evident, as they can be identified to coincide with the sharp rise in $\textit {Ri}_g$. Note in particular that $\textit {Ri}_g > 0.2$ for $0.75 \gtrsim y \gtrsim 1.25$ due to the steep density gradient, although it should also be appreciated that $\textit {Ri}_g$ is singular at $y=1$ due to the zero-valued velocity gradient. Similar trends are observed for both $\textit {Fr}_h$ and $\textit {Re}_B$ which peak in the inner region at $y \approx 0.1$ with $\textit {Fr}_h \approx 0.7$ and $\textit {Re}_B \approx 200$ (note that $\textit {Fr}_h$ is singular at $y=0$). Both quantities decrease monotonically in the outer region until reaching their minimum values in the channel core. A local maximum is present for $\textit {Fr}_h$ at $y=1$ due to the decrease in $u_{rms}$ (figure 2). The ability of a stratified flow to sustain active turbulence is often quantified by a threshold on the buoyancy Reynolds number: $\textit {Re}_B \gtrsim 20$ (Smyth & Moum Reference Smyth and Moum2000). Case S has $\textit {Re}_B < 20$ for $y \gtrsim 0.6$, and $\textit {Re}_B = 0.14$ at the channel centreline. The low $\textit {Re}_B$ and $\textit {Fr}_h$ in the channel core indicate a viscosity-affected (or indeed viscosity-dominated) stratified flow regime, as noted by Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011). Mean profiles indicate (unsurprisingly) that relatively very strong stratification suppresses turbulence in the channel core leading to stable, and close to laminar flow.

Figure 3. Profiles of key dimensionless quantities for Case S, defined in (3.1a–d). Shaded regions are as per figure 2.

Table 1 also reports the Obukhov scale normalised by the channel half-height: $\varLambda = 2\textit {Re}_\tau \textit {Pr}/ \varkappa \textit {Ri}_\tau Nu$, where $\varkappa =0.41$ is the Kármán constant and $Nu = -2\partial _y \bar {\rho }|_w$ is the Nusselt number. The Obukhov scale quantifies the length at which turbulent production due to buoyancy is equal to that due to shear. Case S has $\varLambda > 1$, indicating that shear dominates over buoyancy production in the present simulation.

3.2. Spectra

The dependency of energy spectra on the three-dimensional wavenumber $k_i = (k_x, k_y, k_z)'$ with respective wavelengths $\lambda _i = 2{\rm \pi} /k_i$, and temporal frequency $\omega$, is now investigated. Two-dimensional spatial energy spectra computed on three $y$-normal planes are presented in figures 4 to 6. Energy spectra are reported for fluctuations of the vertical velocity $v'$, density $\rho '$, the momentum flux $u'v'$ and the buoyancy flux $\rho 'v'$, with the latter two statistics omitted from spectra collected at $y=1$ (figure 6), where they are uncorrelated (see figure 2). In each figure, spectra from S data are plotted with shaded contours, while spectra from U data are plotted with solid contour lines.

Figure 4. Two-dimensional premultiplied spatial energy spectra for the (negative) momentum flux (a), vertical velocity (b), density (c) and buoyancy flux (d), at $y^{+} = 100$. Solid contour lines represent U data, shaded contours represent S data. Contour lines are 20 %, 40 %, 60 % and 80 % of respective maximum values.

All spatial spectra of figures 4 to 6 are consistent with findings of Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011). Near the wall, energy spectra are very similar for the stratified and unstratified cases, aside from a reduction of energy associated with large-scale wavelengths in the S data, relative to respective maxima. These suppressed long and wide structures are termed global modes (Del Alamo & Jiménez Reference Del Alamo and Jiménez2003). Conversely, at the channel centre (as shown in figure 6) U and S spectra are markedly different, where Case S has a peak at $(\lambda _x,\lambda _z) \approx (1.5,6)$, significantly larger and more anisotropic than Case U, which peaks at $(\lambda _x,\lambda _z) \approx (1,1)$. The contour lines are also narrower in $\lambda _x$ for Case S. In addition, $v'$ and $\rho '$ energy spectra appear significantly more correlated for Case S, than for the Case U. Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011) argued that such correlation is further evidence of internal waves.

At the intermediate $y=0.75$ (shown in figure 5), two peaks are present for Case S, one of which is due to turbulence, $(\lambda _x,\lambda _z) \approx (0.9,0.5)$, and the other is due to the presence of internal waves, $(\lambda _x,\lambda _z) \approx (1.3,4)$. Interestingly, the peak corresponding to the internal waves has different height dependence for the $v'$ and $\rho '$ spectra. Specifically, the peak in $\rho '$ energy spectra is consistent between $y=1$ and $y=0.75$, while the peak in $v'$ spectra is shifted slightly towards larger streamwise wavelengths at $y=0.75$. Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011) noted that stratification suppresses spanwise and streamwise wavelengths, pushing the peak in $E^{2D}(\lambda _x, \lambda _z)$ associated with turbulence to smaller wavelengths than the unstratified case. Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011) also found that stable stratification damps the tall turbulent structures (global modes), such that dynamics in the outer region are governed by a local scaling based on the Obukhov length $\varLambda (y)$, defined as (Nieuwstadt Reference Nieuwstadt1984)

(3.2)\begin{equation} \varLambda(y) ={-} \frac{1}{\varkappa \textit{Ri}_\tau} \frac{\tau(y)^{3/2}}{q(y)}, \end{equation}

where $\tau (y)$ and $q(y)$ represent the total (dimensionless) momentum and buoyancy fluxes and $\varkappa =0.41$ is the Kármán constant. Here, $\varLambda (y)$ is made dimensionless by the channel half-height $\delta$, and $\varLambda (0)$ corresponds to the values of $\varLambda$ reported in table 1. Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011) found that energy spectra as a function of spanwise and streamwise wavelengths normalised by $\varLambda (y)$ collapsed in the outer region of the flow for their simulations with $\textit {Ri}_\tau = 240 - 960$ if the ratio $y/\varLambda (y)$ was kept constant. This suggests that turbulent processes in the outer layer are governed by local processes, but it is not as yet clear over what range of $y$ this scaling holds, due to the relatively low $\textit {Re}_\tau$ and $\textit {Ri}_\tau$ considered by Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011).

Figure 5. Two-dimensional premultiplied spatial energy spectra for the (negative) momentum flux (a), vertical velocity (b), density (c) and buoyancy flux (d), at $y = 0.75$. Solid contour lines represent U data, shaded contours represent S data. Contour lines are 20 %, 40 %, 60 % and 80 % of respective maximum values.

Figure 6. Two-dimensional premultiplied spatial energy spectra for the vertical velocity (a) and density (b), at $y = 1$. Solid contour lines represent U data, shaded contours represent S data. Contour lines are 20 %, 40 %, 60 % and 80 % of respective maximum values.

Energy spectra as a function of the vertical and streamwise wavenumbers are presented in figure 7. Spectra are calculated using a Hamming window over the interval $y\in [0.2,0.8]$, capturing the turbulent motion in the outer region of the flow. The $x$- and $y$-axes of figure 7 have been scaled to a small range focused on the highest energy wavenumbers, and both halves of the vertical wavenumber spectra are shown here due to asymmetry about $k_y = 0$. Given the energy spectra are premultiplied with streamwise and vertical wavenumbers, energies are expected to be negative for $k_y < 0$ and positive for $k_y > 0$ for $E^{2D}_{vv}k_x k_y$, $E^{2D}_{\rho \rho }k_x k_y$ and $E^{2D}_{\rho v}k_x k_y$ and the opposite is expected for $E^{2D}_{u v}k_x k_y$. Figure 7 indicates that U and S spectra are similar in the region $y\in [0.2,0.8]$, where large differences only arise for the momentum and buoyancy fluxes. Most energy is contained in the lowest wavenumber mode, corresponding to structures with a wavelength equal to the full spatial window, $y\in [0.2,0.8]$. The spectral energy of figure 7 should therefore be interpreted with care, since spectra are clearly affected by the window size and the non-periodic nature of the vertical direction $y$. However, it is argued that this is a physical constraint owing to the channel domain height (or Reynolds number), and these results demonstrate that there is insufficient space for vertical wavenumbers to dominate. Peak energy occurs between a wave vector angle range of $45 - 65^{\circ }$ for all spectra, similar to those observed in the classic experiments of Dohan & Sutherland (Reference Dohan and Sutherland2003), where waves at such an angle were observed propagating away from a turbulent mixed region forced by an oscillating grid. Dohan & Sutherland (Reference Dohan and Sutherland2003) speculated that this preferred propagation of angle was associated with maximisation of the vertical transport of horizontal momentum. Although the angle of wave propagation in this flow simulation is consistent with that observed in their experiments, any argument of direct connection between the key dynamics of the two flows should perhaps be treated with caution due to the absence of large-scale shear in the experiments and the occurrence of the peak at the lowest wavenumber for the simulation data.

Figure 7. Two-dimensional premultiplied energy spectra as a function of vertical and streamwise wavenumbers. Spectra are calculated in the window $y \in [0.2,0.8]$. Contour lines are 20 %, 40 %, 60 % and 80 % of respective minimum (cyan) and maximum (magenta) values. Dashed lines represent equal magnitudes of vertical and streamwise wavenumbers, $k_x = \pm k_y$ ($45^{\circ }$ wave vector angle). Note that spectra are truncated at $k_y=\pm 10$.

A further interesting feature of the vertical-streamwise wavenumber spectra is evident in the spectra for the momentum and buoyancy fluxes, which not only show a new peak in spectra for Case S compared with Case U, but also show it is of a different sign to the other peaks. This new peak occurs at $(k_x,k_y) \approx (20,-40)$, or a wave vector angle of $63^{\circ }$ to the horizontal (figure 7). However, the high streamwise wavenumber of this peak suggests it is not related (at least directly) to the internal waves in the channel core, which have a maximum spectral energy at $\lambda _x \approx 1.5$ or $k_x \approx 4.2$. It is speculated that these high wavenumber downward (away from the core and towards the walls) travelling structures arise from wave breaking at the core edge.

The two-dimensional (2-D) energy spectra are shown as a function of streamwise wavenumber and temporal frequency in figure 8. Two-dimensional spectra are calculated using the method of Welch, with a 50 % overlap Hamming window in the temporal dimension. The window has a length of 4096 timesteps with the full time series containing at least 20 000 snapshots, deemed appropriate through a sensitivity study (not shown). Spectra are subsequently averaged in the spanwise ($z$) direction, and can be directly used to calculate the dispersion relation, $\omega (k_x)$. ${\rm U}$ spectra (figure 8) indicate that all spectral energy closely follows the relation $\omega = \overline{U}_{max} k_x$, where $\overline{U}_{max} k_x$ represents the Doppler shift due to the background flow at the midpoint of the channel $y=1$, a (perhaps) unsurprising result since that is precisely where the spectral data are reported.

Figure 8. Two-dimensional energy spectra as a function of streamwise wavenumber $k_x$ and frequency $\omega$ at $y=1$ for Case U (a,b) and Case S, (c,d). Contours are log-scaled relative energy spectra magnitude. The lines represent different dispersion relations. Subscript max denotes maximum values, and subscript mean denotes the average value in the channel core region ($0.8 < y < 1.2$).

However, the dispersion relation for Case S deviates significantly from $\omega = \overline{U}_{max} k_x$. Indeed, two peaks are observed in the spectral energy on either side of $\omega = \overline{U}_{max} k_x$ for a given $k_x$, with the lower-frequency peak having a spectral energy several orders of magnitude larger than the high-frequency peak. Two additional (model) dispersion relations are shown in figure 8, i.e. $\omega = \overline{U}_{max} k_x \pm N_{max}$ and $\omega = \overline{U}_{mean} k_x \pm N_{mean}$, where the subscript ‘max’ represents the maximum value of a variable (equal to the value at $y=1$) and subscript ‘mean’ represents the spatial average value of a variable only in the channel core region, $0.8 < y < 1.2$. Of course, these are the dispersion relations for linear internal waves with zero vertical wavenumber $k_y$, Doppler shifted by a notional uniform flow equal to either $\overline{U}_{max}$ or $\overline{U}_{mean}$, and maximal magnitude intrinsic frequency (i.e. the local frequency in a fluid at rest, see for example the pedagogical discussion in Bühler Reference Bühler2014) with constant buoyancy frequency $N_{max}$ or $N_{mean}$. It is important to remember that this relationship is highly idealised, and at best is expected to apply only approximately in the real turbulent flow, as both $U$ and $N$ vary significantly even within the central core region. Nevertheless, figure 8 shows that these dispersion relations appear to be relevant, as both ‘forward’ (i.e. with positive intrinsic frequency, and hence positive phase speed relative to the streamwise flow) and ‘backward’ (i.e. with negative intrinsic frequency) propagating internal waves are present in the channel core. As there is (dominant) streamwise velocity throughout the central core region, both types of waves are still convected downstream by the mean flow.

Interestingly, the observed peaks of figure 8 agree with $\omega = \overline{U}_{mean} k_x \pm N_{mean}$, suggesting that internal waves in some sense extend across the entire core-region of the channel and are therefore largely governed by the spatially and temporally averaged buoyancy frequency and velocity in that region ($0.8 < y < 1.2$). It must be remembered that spectral data of figure 8 are extracted only from the very midpoint of the channel at $y=1$.

The highest magnitude waves are contained in streamwise wavenumbers consistent with the 2-D spatial spectra of figure 6, with temporal frequencies of $50 \lesssim \omega \lesssim 150$. Somewhat curiously, the higher extrinsic frequency (i.e. the frequency observed by a stationary observer, using again the nomenclature of Bühler Reference Bühler2014) ‘forward’ propagating mode has a spectral energy several orders of magnitude smaller than the lower extrinsic frequency ‘backward’ propagating mode. This implies that the (streamwise) phase speed $\omega /k_x < \overline{U}_{mean}$ for the dominant low extrinsic frequency mode, which is indeed observed.

The highest energy ($E^{2D}/E^{2D}_{max} \geqslant 0.5$) ‘backward’ propagating internal waves occur between a wavenumber range of $2.5 \lesssim k_x \lesssim 5$ for Case S. The critical levels for these (highest energy) internal waves, where $c = \omega /k_x = \overline{U}_{mean} - N_{mean}/k_x = \overline{U}$, lie in the shaded regions of figures 2 and 3, in the range $0.55 \lesssim y \lesssim 0.75$. The importance of these critical levels will be discussed in § 3.5.

For brevity only 2-D spectra at $y=1$ is reported. Further from the channel centreline spectra for Case S collapse onto the dispersion relation $\omega = U_c k_x$ with $U_c$ representing the planar and temporally averaged velocity $\overline{U}$ at a corresponding $y$ value, suggesting a negligible effect of the stratification. Nevertheless, the data clearly show, surprisingly, that at least the observed periodic structures in the central core region are still well described by linear internal wave theory based around constant streamwise velocity and buoyancy frequency, even though the wave-like structures are undoubtedly finite amplitude, and both the streamwise velocity and buoyancy frequency vary significantly in the vertical. The apparently robust coherence of these waves will be analysed in the following section, using DMD.

3.3. Dynamic mode decomposition

While spectra clearly demonstrate that time-periodic internal-wave-like flow is present in the core region, the 2-D spatial structure of these ‘waves’ is not yet clear. DMD (Schmid Reference Schmid2010) provides a technique to obtain dynamically relevant flow structure from time-resolved data. DMD approximates the flow by

(3.3)\begin{equation} \boldsymbol{q}_k = \sum^{N-1}_{j=0} \boldsymbol{\varPhi}_j \lambda_j^{k-1}, \quad k = 1, \ldots, N-1, \end{equation}

where $\boldsymbol {q}_k$ represents a snapshot of the data at timestep $k$ (including all variables at all points in space, with planar and temporal means subtracted), $\boldsymbol {\varPhi }_j$ represents DMD modes and $\lambda _j$ represents DMD eigenvalues. Alternatively, this can be expressed in terms of the modal growth rate $\sigma _j$ and frequency $\omega _j$:

(3.4)\begin{equation} \boldsymbol{q}_k = \sum^{N-1}_{j=0} \boldsymbol{\varPhi}_j \exp({(\sigma_j - \textrm{i} \omega_j)t_k}), \quad k = 1, \ldots, N-1. \end{equation}

This decomposition is applied to a $z$-normal slice at $z=L_z/2$ and a $y$-normal slice at $y=1$ (detailed in § 2). The dimensionless time period between slice samples is $2\times 10^{-3}$ (20 simulation timesteps), and the spatial resolution of the slices is coarsened by omitting every second point in both directions to reduce memory requirements. For each slice the DMD modes are computed using the method of Schmid (Reference Schmid2010): first, the matrix $\boldsymbol {q}_k$ is constructed from $U$, $V$, $W$ and $\rho$ slice data, with planar and temporal means subtracted. DMD modes are subsequently calculated by first applying a singular value decomposition to $\boldsymbol {q}_k$ before calculation of the DMD modes $\boldsymbol {\varPhi }_j$ and eigenvalues $\lambda _j^{k-1}$. Note that decompositions are calculated without smoothing of the data or truncation of the singular values. Reconstruction of $\boldsymbol {q}_k$ from the DMD results in an $l_2$ norm error of order $1\times 10^{-8}$. Importance of certain modes is assessed by their modal energy, $E_{\varPhi,j} = ||\boldsymbol {\varPhi }_j||^{2}$. As simulations are of fully converged stationary flow all decompositions presented here have growth rates $\sigma _j \approx 0$ which will subsequently not be reported.

Figure 9 shows DMD for Case S ($\textit {Ri} = 480$) on a $z$-normal slice. DMD modal energies have been plotted with 1-D temporal energy spectra, $E^{1D}_{vv}$, obtained at $y=1$. Generally $E_\varPhi$ decreases with increasing $\omega$, seemingly in disagreement with temporal energy spectra (figure 9). This is explained by noting that DMD modal energies are calculated over the full channel height, while spectral energies are only calculated at $y=1$. When narrowing the $\omega$ axis range, figure 9 shows that modal energies, $E_{\varPhi }$, lead to the same local peaks in frequency as the temporal spectra $E_{vv}^{1D}$, corresponding to spatial forcing of the domain length, and the internal waves. The spatially forced local peaks arise due to the finite periodic domain size, where streamwise wavenumbers ($k_x$) in multiples of $2{\rm \pi} /L_x$ are ‘favoured’ over others. The peaks in temporal frequency correspond to these structures passing through the domain at some convective velocity $U_c$. Note that these local peaks arising from spatial forcing are also present in the 2-D spectra of figure 8 but are partially hidden by the logarithmic scaling. Examples of $U$, $V$ and $\rho$ components of $\varPhi$ are also reported in figure 9 and demonstrate that the internal waves are highly coherent, spanning the full height of the core region. The vertical extent of waves suggests a plausible reason for the surprising agreement between the (spatially localised) spectral data from $y=1$ and the linear dispersion relation for internal waves based upon the buoyancy frequency and velocity averaged over the channel core region (figure 8). Further, it is noted that the density and vertical velocity modes are out of phase by ${\rm \pi} /2$, consistent with the polarisation relations for monochromatic internal gravity waves, and indeed the instantaneous visualisations of Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011).

Figure 9. DMD modes for Case S on a $z$-normal slice at $z=L_z/2$. DMD modal amplitudes $E_{\varPhi,j} = ||\boldsymbol {\varPhi }_j||^{2}$ as a function of temporal frequency $\omega _j$ are plotted in (a) with black dots. Panel (b) presents the same data as (a) with axis limits focused on the internal waves. One-dimensional temporal vertical velocity energy spectra, $E_{vv}^{1D}$, calculated at $y=1$, are also plotted with a red line. The spatial structure of the highlighted mode in (a,b) ($\bullet$, blue) is shown in (c–f), where $\varPhi _U$ is the streamwise velocity component, $\varPhi _V$ is the vertical velocity component and $\varPhi _\rho$ is the density component. Streamlines calculated from $\varPhi _U$ and $\varPhi _V$ are reported in (e). The observed structure of this example mode ($\bullet$, blue) is typical of all modes corresponding to local peaks in $E_{vv}^{1D}$, deviating only by the streamwise wavenumber.

The dispersion relation can also be obtained empirically from DMD using $y$-normal slices at $y=1$, an example of which is presented in figure 10. Modal amplitudes are in excellent agreement with temporal energy spectra. Further, all modes are characterised by structures that span the full width of the domain, and have a dominant underlying wavenumber $k_x$. This is demonstrated in figure 10 by presenting spatial spectra in the streamwise direction of two example modes, where the dominant peaks are at wavenumbers of $k_x \approx 2.5$ and $k_x \approx 3.5$ respectively. Note also that there is an additional much smaller peak for both signals at a slightly lower $k_x$ value. These second peaks are approximately three orders of magnitude smaller in spectral energy, but are a consistent feature of modes with $\omega \lesssim 200$. These secondary small-amplitude peaks correspond to the higher-frequency ‘forward’ internal wave-like structures observed in the spatio-temporal spectra of figure 8, while the high-amplitude peaks correspond to the lower-frequency ‘backward’ propagating internal wave-like structures.

Figure 10. DMD modes for Case S on a $y$-normal slice at $y=1$. DMD modal amplitudes $E_{\varPhi,j} = ||\boldsymbol {\varPhi }_j||^{2}$ as a function of temporal frequency $\omega _j$ are plotted in (a) with black dots. One-dimensional temporal vertical velocity energy spectra, $E_{vv}^{1D}$, calculated at $y=1$, are also plotted with a red solid line. The spatial structure of the two highlighted modes in (a) ($\blacksquare$, $\blacktriangleright$) are shown in (b,c), where $\varPhi_V$ is the density component. The observed structure of these two example modes is typical of all modes with local peaks in $E_\varPhi$. Energy spectra of $\varPhi_V$ as a function of streamwise wavenumber $k_x$ are presented in (d) for both Mode 1 ($\blacksquare$) and Mode 2 ($\blacktriangleright$). Note that both signals are bimodal, with a dominant high $k_x$ peak and a weaker low $k_x$ peak.

The spatial spectra of DMD modes at $y=1$ are subsequently used to determine the dispersion relation empirically, by obtaining the $k_x$ values of all notable peaks in the spectral energy (achieved using a peak-finding algorithm). The dispersion relation for Case S is plotted in figure 11. There is clear agreement between the dispersion relation obtained by DMD and 2-D spectra (figure 8). However, the DMD dispersion relation shows the two separate peaks, on either side of $\omega = \overline{U}_{max} k_x$ more clearly. The flow is comprised of a series of internal waves, with peak amplitudes occurring over a temporal frequency range of $50 \lesssim \omega \lesssim 200$. The lower-frequency dominant waves closely follow the dispersion relation $\omega = \overline{U}_{mean} k_x - N_{mean}$, while the higher-frequency waves have a significantly less well-defined dispersion relation that lies between $\omega = \overline{U}_{mean} k_x + N_{mean}$ and $\omega = \overline{U}_{max} k_x + N_{max}$.

Figure 11. Dispersion relation for Case S, obtained from DMD modes at $y=1$. Data are $\varPhi _v$ modes (a) and $\varPhi _\rho$ modes (b). Figures are coloured by spectral energy, normalised by the maximum spectral energy.

Therefore, simulations have so far shown that stratified channel flow leads to three well-defined regions in the flow, i.e. the inner region, the outer region and the channel core region, characterised by a strong density gradient and a relatively low turbulence intensity. Spectra and DMD reveal two sets of structures which are sustained in the channel core, and can be interpreted, at least loosely, as linear internal waves that are highly coherent and span the full core-region height. The dominant mode occurs at a lower frequency, (and hence lower phase speed for given streamwise wavenumber) and has an energy several orders of magnitude higher than the (also present) higher-frequency mode. However, the mechanisms driving and sustaining these periodic wave-like structures are as yet unclear. In particular, spectra in the inner and outer regions show that vertical wavenumbers are negligible in this flow, suggesting that the generation mechanisms for these waves are indeed different from those in the non-sheared grid-turbulence experiments reported by Dohan & Sutherland (Reference Dohan and Sutherland2003), as discussed above. The spectral energy dependence on vertical wavenumbers is found to be dominated by the lowest wavenumber modes (figure 7) suggesting that modes are physically restricted by the channel height (or Reynolds number). It is also important to remember the strong (wall-normal) spatial variation in streamwise velocity and buoyancy frequency in the central core region. This variation is crucial, as in the following section the waves are demonstrated to arise from the sensitivity of the mean profiles to perturbation growth, using two different, although related and complementary, analysis techniques.

3.4. Stochastic forcing of the linearised equations

The first hypothesis that is tested is the proposition that the internal ‘waves’ are generated and sustained by an inhomogeneous incoherent turbulent forcing. This process may be thought of as either a filtering of random forcing from the turbulence in the inner and outer regions, or equivalently a sensitive response of particular structures in the central core region to such forcing. With either viewpoint, the outcome is proposed to be that essentially random turbulence continually forces small perturbations from mean flow profiles over a range of wavenumbers and frequencies, and what emerges are the most favoured or sensitive perturbation structures, which exhibit the strongest ‘response’ (in the sense that the perturbation energy grows). Specifically, it is proposed that the ‘waves’ are generated by the amplification of certain frequencies and the fastest growing modes dominate consistently with the heuristic linear internal wave dispersion relation using the mean streamwise velocity and mean buoyancy frequency within the central core region. This hypothesis is tested by investigating the response of the linearised momentum, buoyancy and continuity equations to white-noise forcing (uncorrelated in $t$ and $x$, but correlated in $y$), which is a highly idealised yet still instructive model for the (assumed) incoherent forcing due to turbulence generated principally in the inner and outer regions of the flow. The stochastically forced linearised Navier–Stokes equations have been successfully adopted to analyse the stability of wall-bounded shear flows (see e.g. Farrell & Ioannou Reference Farrell and Ioannou1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; Hwang & Cossu Reference Hwang and Cossu2010; Jovanović Reference Jovanović2021), with previous work often adopting an approach based on linear dynamical systems analysis (see Jovanović (Reference Jovanović2021) for a recent review of this methodology). This work evolves the stochastically forced linearised Navier–Stokes equations in time directly by introducing a streamfunction formulation. In this way, the governing equations are a direct precursor to the viscous Taylor–Goldstein analysis of § 3.5.

Following Liu, Thorpe & Smyth (Reference Liu, Thorpe and Smyth2012) the dimensionless linearised momentum, buoyancy and continuity equations for this system are

(3.5)\begin{gather} \frac{\partial {u}}{\partial {t}} + \overline{U} \frac{\partial {u}}{\partial {x}} + v \frac{\textrm{d} {\overline{U}}}{\textrm{d} {y}} ={-} \frac{\partial {p}}{\partial {x}} + A_h \frac{\partial^{2} u}{\partial x^{2}} + A_v \frac{\partial^{2} u}{\partial y^{2}} + \frac{\textrm{d} A_v}{\textrm{d} y} \frac{\partial {u}}{\partial {y}}, \end{gather}

(3.6)\begin{gather}\frac{\partial {v}}{\partial {t}} + \overline{U} \frac{\partial {v}}{\partial {x}} ={-} \frac{\partial {p}}{\partial {y}} + b + A_h \frac{\partial^{2} v}{\partial x^{2}} + A_v \frac{\partial^{2} v}{\partial y^{2}} + \frac{\textrm{d} A_v}{\textrm{d} y} \frac{\partial {v}}{\partial {y}}, \end{gather}

(3.7)\begin{gather}\frac{\partial {u}}{\partial {x}} + \frac{\partial {v}}{\partial {y}} = 0, \end{gather}

and

(3.8)\begin{equation} \frac{\partial {b}}{\partial {t}} + \overline{U} \frac{\partial {b}}{\partial {x}} + v \frac{\textrm{d} {\overline{B}}}{\textrm{d} {y}} = K_h \frac{\partial^{2} b}{\partial x^{2}} + K_v \frac{\partial^{2} b}{\partial y^{2}} + \frac{\textrm{d} K_v}{\textrm{d} y} \frac{\partial {b}}{\partial {y}}, \end{equation}

where $u$, $v$, $p$ and $b$ denote streamwise velocity, vertical velocity, pressure and buoyancy perturbations and $\overline{U} = \overline{U}(y)$, and $\overline{B} = \overline{B}(y) = -\textit {Ri}_\tau \bar {\rho }(y)$ represent planar and temporally averaged streamwise velocity and buoyancy, obtained from the full 3-D simulations. Note that this formulation allows for wall-normal variation in both the streamwise velocity and the density (or equivalently buoyancy frequency). The buoyancy perturbation relates to the density perturbation by $b = - \textit {Ri}_\tau \rho$. $A_h$, $A_v$, $K_h$ and $K_v$ represent horizontal and vertical diffusivities of momentum ($A$) and buoyancy ($K$), equal to the sum of both viscous and turbulent diffusion.

The streamfunction $\psi$ is introduced, where $u = \partial \psi / \partial y$ and $v = - \partial \psi / \partial x$, transforming the equations to

(3.9)\begin{equation} \frac{\partial {}}{\partial {t}}\nabla^{2} \psi + \overline{U} \frac{\partial {}}{\partial {x}}\nabla^{2} \psi - \frac{\textrm{d}^{2} {\overline{U}}}{\textrm{d} {y}^{2}}\frac{\partial {\psi}}{\partial {x}} ={-} \frac{\partial {b}}{\partial {x}} + \frac{1}{\textit{Re}_\tau} \nabla^{4} \psi + \mathcal{W}(y), \end{equation}

and

(3.10)\begin{equation} \frac{\partial {b}}{\partial {t}} + \overline{U} \frac{\partial {b}}{\partial {x}} - \frac{\textrm{d} {\overline{B}}}{\textrm{d} {y}}\frac{\partial {\psi}}{\partial {x}} = \frac{1}{\textit{Pr} \textit{Re}_\tau} \nabla^{2} b. \end{equation}

Here, the diffusivities of momentum have been separated into viscous ($1/\textit {Re}_\tau$) and turbulent components, the turbulent part of which is modelled by $\mathcal {W}(y)$, a white-noise function correlated in $y$. The streamfunction formulation allows the white-noise forcing to be imposed while maintaining mass continuity. It is assumed that the buoyancy fluxes have a minor role in the inner region of the boundary layer, such that the turbulent components of $K_h$ and $K_v$ can be neglected while molecular diffusion ($1/\textit {Pr} \textit {Re}_\tau$) is retained. The dimensionless parameters are the same as for the Case S simulation, and so $\textit {Re}_\tau = 550$, $\textit {Ri}_\tau = 480$ and $\textit {Pr} = 1$.

The system of equations is discretised and solved using Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). Equations are solved on a 2-D domain of size $L_x \times L_y = 8 {\rm \pi}\times 2$, discretised using a Fourier basis in the periodic ($x$) direction and a Chebyshev basis in the vertical ($y$) direction, dealiased using the $3/2$ rule. The grid is chosen to match the LES resolution (table 1). At each grid point $\mathcal {W}(y)$ is calculated by a normal distribution around a mean value of $0$ and a standard deviation of $1/\sqrt {\Delta t}$, where $\Delta t$ is the simulation timestep. To better approximate turbulence in the inner region of the flow the normal distribution is multiplied by $\nu _t/\nu _{t,{max}}$, where $\nu _t = -\overline {u'v'}/\partial _y \overline{U}$ quantifies the eddy viscosity of the 3-D nonlinear simulations, such that $\mathcal {W}(y)$ is strongest in the inner region of the flow and zero at the walls and channel centreline (visualised in figure 12). The initial value problem, initially at rest, is integrated in time using a first-order Runge–Kutta scheme for a total simulation time of $T=80$, where the flow reaches a pseudo-steady state at $T \approx 10$.

Figure 12. Snapshots of solutions to the white-noise-forced linearised equations. Variables are the streamwise velocity (a), vertical velocity (b) and buoyancy (c), perturbations, and the white-noise function (d).

A snapshot of the simulation is presented in figure 12. The wave-like structure is clearly seen to develop due to the stochastic forcing away from the channel core. The dispersion relation at the channel centreline is plotted in figure 13. The dominant mode arising in the (real) nonlinear system is clearly well captured by this stochastically forced linearised system (compare with figure 8), indicating that it is insensitive to the particular forcing in the inner/outer regions of the flow. Note that the solutions to the stochastically forced linearised system are particularly insensitive to the form of the forcing term $\mathcal {W}(y)$. Alongside white-noise forcing at the grid scale, further noise functions were tested, including ‘red noise’ with time correlation, uniform distributions without correlation with $\nu _t$ and spatially correlated noise calculated by up-sampling noise generated on a coarser grid. None of these functions made an appreciable difference to the dispersion relation, and are therefore not reported.

Figure 13. Two-dimensional energy spectra as a function of streamwise wavenumber $k_x$ and frequency $\omega$ at $y=1$, for the white-noise-forced system. Contours are log-scaled relative energy spectrum magnitude. The lines represents different dispersion relations. Subscript max denotes maximum values, and subscript mean denotes the average value in the channel core region ($0.8 < y < 1.2$).

These results suggest that wavelike structures arise as a sensitive response of the mean flow profiles to incoherent forces. Inspection of the gradient Richardson number $\textit {Ri}_g$ (figure 3) shows that $\textit {Ri}_g \approx 0.2$, and so it is perhaps unsurprising that perturbations can grow significantly. The mean flow profiles develop as a function of both Reynolds number and Richardson number. The particular conditions that allow waves and turbulence to coexist in the channel arise from the combination of the imposed density gradient and the near-wall turbulence. In the weakly stratified regime the near-wall dynamics are governed by the turbulence generated by near-wall shear, which appears largely unaffected by stratification, due to the imposed friction Reynolds number. This leads to the well-mixed inner and outer regions of the flow. However, the imposed density gradient at the boundaries must be satisfied somewhere in the flow, which subsequently and inevitably leads to the sharp density gradient at the channel core, where the flow is furthest from the near-wall shear. This leads to the well-defined inner, outer and channel core regions of figure 2. The system acts as a spectral filter such that incoherent perturbations nearer the walls cause ringing and lead to coherent periodic wavelike structures in the core.

However, the high-frequency mode obtained through the white-noise-forced system more closely matches the dispersion relation $\omega = \overline{U}_{max} k_x$ rather than a dispersion relation for linear internal waves. This discrepancy will be discussed in the following section, where rather than investigating the linear response to incoherent forcing we consider the linear (normal mode) instability of the mean profiles.

3.5. Solutions to the viscous Taylor–Goldstein equations

Solutions are sought to the viscous Taylor–Goldstein equations (Liu et al. Reference Liu, Thorpe and Smyth2012; Lian, Smyth & Liu Reference Lian, Smyth and Liu2020) which are derived from the linearised momentum, continuity and buoyancy transport equations of § 3.4 by assuming perturbations take the form

(3.11a,b) \begin{equation} v = \mathcal{R}[\hat{v}(y) \exp({\lambda t + \textrm{i} k_x x})], \quad b = \mathcal{R}[\hat{b}(y) \exp({\lambda t + \textrm{i} k_x x})], \end{equation}

where $\lambda = \sigma - \textrm {i} \omega$ is the complex growth rate, $\hat {v}(y)$ and $\hat {b}(y)$ represent complex vertical structure functions and it is assumed the wave vector is aligned with the $x$ direction (${k_y = 0}$). Following Lian et al. (Reference Lian, Smyth and Liu2020), and ensuring dimensional consistency with the full 3-D simulations, the dimensionless viscous Taylor–Goldstein equations are

(3.12)$$\begin{gather} (\lambda + \textrm{i}k_x\overline{U}) \nabla^{2} \hat{v} - \textrm{i}k_x \frac{\textrm{d}^{2} \overline{U}}{{\textrm{d} y}^{2}}\hat{v} = T_w \hat{v} - k_x^{2} \hat{b}, \end{gather}$$

(3.13)$$\begin{gather}(\lambda + \textrm{i}k_x\overline{U}) \hat{b} + \frac{\textrm{d} \overline{B}}{\textrm{d} y}\hat{v} = T_b \hat{b}, \end{gather}$$

with viscous and diffusive operators

(3.14)$$\begin{gather} T_w = \frac{d^{2}}{{\textrm{d} y}^{2}}\left(A_v \frac{\textrm{d}^{2}}{{\textrm{d} y}^{2}} \right) - k_x^{2} \frac{\textrm{d}}{\textrm{d} y} \left[ (A_h + A_v) \frac{\textrm{d}}{\textrm{d} y} \right] + k_x^{4} A_h, \end{gather}$$

(3.15)$$\begin{gather}T_b = \frac{\textrm{d}}{\textrm{d} y} \left(K_v \frac{\textrm{d}}{\textrm{d} y} \right) + k_x^{2} K_h, \end{gather}$$

where $A_h$ and $A_v$ are horizontal and vertical viscosities, and $K_h$ and $K_v$ are horizontal and vertical diffusivities. Although in principle this form of the equations allows for the modelling of turbulence through a turbulent diffusion closure, here, only molecular diffusion is allowed, such that $A_h = A_v = 1/\textit {Re}_\tau$ and $K_h = K_v = 1/\textit {Pr} \textit {Re}_\tau$. Inclusion of eddy viscosity (not shown) actually does not affect the obtained dispersion relation and only has a minor smoothing effect on the structure functions at the critical levels, where $\omega = \overline{U} k_x$. Equations (3.12) and (3.13) can be reconstructed as a generalised differential eigenvalue problem, and solved numerically using finite differencing for given profiles of $\overline{U}(y)$ and $\overline{B}(y)$, and wavenumber $k_x$ (see Lian et al. (Reference Lian, Smyth and Liu2020) for further details). The eigenvalue problem is solved on a uniform grid consisting of $401$ points over a wavenumber range of $k_x = 0.5 - 80$. For each $k_x$ $801$ modes are found with corresponding complex structure functions $\hat {v}$ and $\hat {b}$, and the eigenvalues $\lambda = \sigma - \textrm {i}\omega$.

The dispersion relation $\omega (k_x)$ obtained from the eigenvalue problem is plotted in figure 14, as well as the eigenstructures of the two modes at $k_x = 4$ that have the largest growth rates. The dispersion relation is coloured by $\sigma$, and is negative for all modes such that the system is stable to small perturbations, with maximum growth rates close to but less than $\sigma = 0$. There is good agreement between the (real) nonlinear simulations of figure 8, the stochastically forced linear solutions of figure 13 and the viscous Taylor–Goldstein solutions of figure 14, all of which obtain the dispersion relation $\omega = \overline{U}_{mean} k_x - N_{mean}$ for the dominant ‘backward’ travelling mode. The spatial structure of the two highest growth rate modes for $k_x = 4$, labelled Mode 1 and Mode 2, are also plotted in figure 14. Mode 1 lies on the linear dispersion relation for ‘backward’ travelling internal waves, and shows that there is a critical level at $y \approx 0.7$, where $\omega = \overline{U} k_x$. This approximately corresponds to the point at which $\textit {Ri}_g = 0.2$ in figure 3, before rapidly increasing in the channel core. Further, the spatial structure is qualitatively similar when comparing $v(x,y)$ and $b(x,y)$ for Mode 1 in figure 14 to the stochastically forced instantaneous snapshot of figure 12, and the DMD modes of figure 9. There is therefore a plethora of evidence to suggest that turbulent perturbations ring the system at the edge of the channel core and generate coherent internal waves (figures 2, 3, 8, 11, 13 and 14); the critical levels for the most energetic waves occur in a turbulent region of the flow at $0.55 \lesssim y \lesssim 0.75$ (figures 2 and 5). Any wave breaking in the flow is therefore likely to be induced by turbulent events rather than a saturation of energy at a critical level.

Figure 14. Solutions to the viscous Taylor–Goldstein equations. The dispersion relation $\omega (k_x)$ is plotted in (a), with modes coloured by their growth rate. The lines represent different (linear) dispersion relations, where subscript max denotes maximum values, and subscript mean denotes the average value in the channel core region ($0.8 < y < 1.2$). The two highlighted modes correspond to those with the highest growth rates for $k_x = 4$. Their complex structure eigenfunctions $\hat {v}(y)$ and $\hat {b}(y)$ are reported in (b–e), with real components marked in black, complex components marked in grey, and the origin marked by the dotted black lines. Critical lines are marked in red, where the (streamwise) phase speed $\omega /k_x = \overline{U}$. In addition, the spatial structures $v(x,y)$ and $b(x,y)$, obtained from (3.11a,b) are plotted for one wavelength, in (b–e).

However, neither the stochastically forced system nor the viscous Taylor–Goldstein solutions reproduce the high-frequency modes, $\omega = \overline{U}_{mean} k_x + N_{mean}$ as per the nonlinear, 3-D simulations. Mode 2 of figure 14 actually lies close to the dispersion relation $\omega = \overline{U}_{max} k_x$. The spatial structure, particularly the buoyancy, is much narrower than for Mode 1, thus implying stronger gradients in the structure of this mode, and also a more specific sensitivity of the mode to the particular structure of the (averaged) background flow in the vicinity of the peak in the amplitude of the perturbation. The stronger values of the gradient suggest that Mode 2 may be more sensitive to nonlinearities, while the fact that the (actual) background flow varies in time also suggest that the particular structure of the background flow may not be in the required form conducive to growth of Mode 2. Although speculative, it therefore seems plausible that the actual flow conditions may suppress Mode 2 while allowing Mode 1 to persist. The disagreement between the high-frequency mode obtained from the linearised systems and the nonlinear simulations suggests that the observed ‘forward’ propagating internal waves arise either through some nonlinearity in the system, or possibly a dependence on instantaneous deviations away from the mean flow profiles. Possible reasons for the disagreement will be further discussed in § 4.

3.6. Analysis of instantaneous data

The (apparently essentially linear) mechanisms generating and sustaining the low-frequency mode can now be summarised. The background mean flow profiles, arising from the well-mixed outer region coupled with the steep density gradient in the channel core, are linearly unstable. Wave-like structures also arise as a sensitive filtered response to incoherent turbulence perturbations at the core edge, over a range of wavenumbers and frequencies, generating internal waves which approximately satisfy the linear dispersion relation $\omega = \overline{U}_{mean} k_x - N_{mean}$. The observed appearance of the high-frequency low energy mode in spectra and DMD is less well explained, given that the two sets of linearised equations considered above do not appear to reproduce these waves. In this section, the flow is considered from a different viewpoint, with an investigation of the nonlinear transient flow dynamics, with a particular focus on the interactions between the turbulent outer region and the buoyancy dominated channel core, associated with specific coherent, yet inherently spatio-temporally intermittent, hairpin vortical structures.

Streamwise-wall-normal planar snapshots of the instantaneous local TKE dissipation rate, $\varepsilon = \textit {Re}_\tau ^{-1}\partial _i u_j'\partial _i u_j'$, are plotted in figures 15 and 16 for Case S. The two sets of snapshots track turbulent events at the edge of the channel core in a Lagrangian reference frame, with a horizontal coordinate given by $x' = x - U_c t$ where $U_c$ is the spatially and temporally averaged velocity in the interval $0.7< y<1.3$ ($x' \in [0,L_x]$ due to periodicity). Grey regions indicate an unstable density field, where dense fluid is above light fluid, indicative of wave-breaking events.

Figure 15. Time series of $z$-slice data focused on a turbulent ejection event in a Lagrangian frame of reference, where $t=0$ represents the start of the event. The streamwise coordinate is given by $x' = x - U_c t$ where $U_c$ is the temporally and spatially averaged velocity in the given window. Filled contours are TKE dissipation rate, overlaid by iso-contours of constant density with a difference between them of $0.1$. The thick solid contour represents $\rho = 0$. Negative contours are dashed lines, positive contours are solid lines. The grey regions indicate an unstable density field. The dashed box focuses on a single ejection event.

Figure 16. Time series of $z$-slice data focused on a turbulent ejection event in a Lagrangian frame of reference, where $t=0$ represents the start of the event. The streamwise coordinate is given by $x' = x - U_c t$ where $U_c$ is the temporally and spatially averaged velocity in the given window. Filled contours are TKE dissipation rate, overlaid by iso-contours of constant density with a difference between them of $0.1$. The thick solid contour represents $\rho = 0$. Negative contours are dashed lines, positive contours are solid lines. The grey regions indicate an unstable density field. The dashed box focuses on a single ejection event.

Figure 15 shows a strong ejection of turbulent fluid from the outer region into the channel core, nearly penetrating downwards to the centreline of the flow. The event is highly dissipative, causing local ‘wave breaking’ as it passes through the core edge. All observed wave-breaking events appear to arise due to the same mechanism. Therefore, in this paper, ‘wave breaking’ refers to the over-turning of density contours due to a saturation in energy associated with highly dissipative ejection events passing from the outer region and into the channel core. The bulk of the flow in the channel core has $\varepsilon \approx 0$; all dissipation in the channel core is due to turbulent activity penetrating the core edge. The ejection also perturbs waves as it passes through the core edge.

Figure 16 shows a much more active period of the flow. Here, an ejection event is tracked in an already highly turbulent period. In this case high-amplitude waves are present in the channel core before the ejection event, and turbulent ejections appears to act in phase with the observed internal waves. The density gradient in the core is sharpened during highly turbulent activity, but once events have dissipated the contours relax and the density field smooths out (as is apparent at $t=0.35$). The event of figure 16 appears to be a local source of high dissipation, particularly at $t=0.1$, in contrast to the clear ejection of fluid in figure 15. However, this is due to the purely planar visualisation which cannot capture the three-dimensionality of the ejection, evidenced through the appearance of the ejection ‘tail’ at $t=0.3$ in figure 16. All observed ‘wave-breaking’ events are consistent and appear to be induced by ejection events rather than linear instability.

Contours of $v'$ are also plotted for the turbulent event of figure 16 in figure 17. The turbulent ejection clearly acts in phase with the vertical velocity of the internal wave. The turbulent structure is fairly consistent in height during the first three snapshots ($t = 0$ to $t=0.1$), with $y \approx 0.9$, due to the large region of negative vertical velocity above it, suppressing upward vertical flow. This changes during $t=0.15$ to $t=0.2$, where the vertical velocity in the channel core is positive and the turbulent structure rapidly moves upwards. At $t=0.25$ to $t=0.3$ the structure is suppressed by the following wave crest and dissipates its remaining TKE. The opposite is true for the upper region of the channel, where turbulent events are typically ejected when the vertical velocity in the core is locally negative. As such, ejections on either side of the channel core rarely interact since they are typically out of phase. This process suggests that when high-amplitude waves are present in the channel core, turbulent ejections are enhanced by the wave motion which draws in turbulent structures from the outer region.

Figure 17. Time series of $z$-slice data focused on a turbulent ejection event in a Lagrangian frame of reference, where $t=0$ represents the start of the event. The streamwise coordinate is given by $x' = x - U_c t$ where $U_c$ is the average velocity in the given window. Filled contours are vertical momentum sources due to fluctuating pressure and density fields, overlaid by iso-contours of constant density with a difference between them of $\Delta \rho = 0.1$. The thick solid contour represents $\rho = 0$. Negative contours are dashed lines, positive contours are solid lines. The grey regions indicate turbulent events ($\varepsilon > 5$). The dashed box focuses on a single ejection event.

The wave-breaking mechanisms can be analysed by assessing the budgets of spanwise vorticity, $\varOmega _z$, where $\boldsymbol {\varOmega } = \boldsymbol {\nabla } \times \boldsymbol {U}$ is the three-dimensional vorticity field. Taking the curl of the momentum equation leads to

(3.16)\begin{equation} \frac{\textrm{D} \boldsymbol{\varOmega}}{\textrm{D}t} = (\boldsymbol{\varOmega} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{U} - \textit{Ri}_\tau \boldsymbol{\nabla} \rho \times \hat{\boldsymbol{y}} + \frac{1}{\textit{Re}_\tau} \nabla^{2} \boldsymbol{\varOmega}, \end{equation}

demonstrating that the vorticity can change due to vortex stretching and baroclinic torques as well as due to viscous effects. Tracking the transport of vorticity yields insight into how wave breaking occurs in the system. It is particularly insightful to identify the mechanisms by which $\varOmega _z$ is generated at the edges of the channel core. Contours of the spanwise vorticity field and corresponding budgets are presented in figure 18 for a single snapshot. The ejection event in the upper region of the channel core is characterised by strong positive spanwise vorticity at the interface between the turbulent structure and the channel core region; it is this positive vorticity that leads to wave breaking at the core edge. Wave-breaking events are associated with negative vorticity for $y<1$ and positive vorticity for $y > 1$. As shown in figure 18 spanwise vorticity generation is dominated by vortex stretching, particularly the term involving the spanwise vorticity $\varOmega _z \partial _z U_z$, indicating that the turbulent structure is highly three-dimensional and has a strong spanwise vorticity field. The contribution from the baroclinic vorticity term, $-(\textit {Ri}_\tau \boldsymbol {\nabla } \rho \times \hat {\boldsymbol {y}})_z = -\textit {Ri}_\tau \partial _x \rho$ has a leading-order effect and is of the same sign as the vortex stretching, but is approximately 2–3 times smaller in magnitude. Both the spanwise vortex stretching and the baroclinic source act to generate positive vorticity on the downstream edge and the core-facing edge of the ejection, and negative vorticity at the upstream edge of the ejection, subsequently suppressing downward flow and overturning the local density contours. The vortex stretching term involving the wall-normal vorticity $\varOmega _y \partial _y U_z$ is also significant at leading order, acting to oppose spanwise vorticity generation.

Figure 18. Vorticity field and budgets on a $z$-normal slice at $t=0.15$ of figure 15. Solid iso-contours are density with a difference between them of $0.1$. The thick solid contour represents $\rho = 0$. Negative contours are dashed lines, positive contours are solid lines. The dashed box focuses on a single ejection event.

The interface between the channel core and outer region is clearly intermittent with occasional high dissipation events that trigger high-amplitude waves. The intermittency of Case S is quantified in figure 19. The time dependences of several statistics are presented in the window $x' \in [0,{\rm \pi} ]$ for the lower portion of the channel core, $y \in [0.8,1]$, and the upper portion of the outer region, $y \in [0.5,0.8]$. Statistics are calculated from the $z$-normal slice data, snapshots of which are plotted in figures 15 to 18. The length $l_c$ is defined as the length of the density contour $\rho = 0.2$ in the interval $x' \in [0,{\rm \pi} ]$, normalised by its minimum possible value of ${\rm \pi}$, such that $l_c=1$ corresponds to a perfectly flat density contour. High values of $l_c$ correspond to wave-breaking events where the contour is stretched, overturned and made multi-connected by the turbulent events. Spatially averaged TKE and buoyancy dissipation rates are also plotted, defined by $\langle \varepsilon \rangle = \textit {Re}_\tau ^{-1} \langle \partial _i u_j'\partial _i u_j' \rangle$ and $\langle \mathcal {X} \rangle = - \textit {Ri}_\tau \textit {Pr}^{-1}\textit {Re}_\tau ^{-1} \langle \partial _i \rho ' \partial _i \rho ' \rangle / \partial _y \bar {\rho }$, where the buoyancy dissipation rate is (appropriately) normalised by the square of the characteristic buoyancy frequency (Caulfield Reference Caulfield2021). The angled brackets denote local spatial averaging in the interval $x' \in [0,{\rm \pi} ]$ and either $y \in [0.8,1]$ or $y \in [0.5,0.8]$. Grey regions are based upon strong peaks in $l_c$, indicative of vigorous perturbation of the central core region. The intermittency of the flow is clear in the time series of $l_c$, $\langle \varepsilon \rangle$ and $\langle \mathcal {X} \rangle$, averaged in the interval $y \in [0.8,1]$. The large distinct peaks in these signals correspond to events similar to those displayed in figures 15 and 16, where turbulent ejections cause high dissipation and wave breaking. There are quiet periods of low dissipation, short periods of high dissipation, and longer periods where the flow is quite active. Note that quiet periods have an ${O}(\delta /u_\tau )$ time scale (or $\Delta t \approx {O}(1)$), necessitating the long integration times adopted in the present study. Peaks in $l_c$ clearly correspond to peaks in $\langle \varepsilon \rangle$, and $\langle \mathcal {X} \rangle$ in $y \in [0.8,1]$. Further, there is a strong correlation between activity in the core $y \in [0.8,1]$ and activity in the outer region $y \in [0.5,0.8]$. Specifically, high dissipation in the channel core corresponds to strong turbulent activity in the outer region of the flow. There also appears to be a lag between events in the core and corresponding peaks in the outer region, where activity in the core ‘leads’ activity in the outer layer. This seemingly contradictory behaviour arises due to the Lagrangian frame of reference where flow in the outer region is convected more slowly than flow in the core. Ejection events therefore appear in the core slightly before they appear in the outer region (see e.g. $t\approx 14.25$ of figure 19).

Figure 19. Intermittency of the spatially averaged dissipation rates of TKE ($\varepsilon$) and buoyancy ($\mathcal {X}$) in the lower portion of the channel core ($y \in 0.8,1.0$), and the near-core portion of outer region ($y\in [0.5,0.8]$), in the interval $x' \in [0,{\rm \pi} ]$; $l_c$ is the length of the contour $\rho =0.2$ normalised by ${\rm \pi}$, such that $l_c=1$ corresponds to a perfectly flat contour. Grey regions correspond to peaks in $l_c$, indicative of wave-breaking ejection events.

Time-series data therefore suggest that dissipative events in the channel core are dependent on turbulent activity in the outer region of the flow. These time series also show that the dissipation rates of TKE and buoyancy are well correlated and of a similar order, with $\langle \mathcal {X}\rangle$ typically half the value of $\langle \varepsilon \rangle$. Correlations between turbulent events and locally averaged density and velocity gradients ($\langle \partial _y \rho \rangle$ and $\langle \partial _y U \rangle$) are less clear, and subsequently omitted from figure 19.

The three-dimensionality of the turbulent structures is visualised by iso-contours of swirl strength $\lambda _{ci}$ in figure 20, where $\lambda _{ci}$ is the imaginary part of complex eigenvalues for the velocity gradient tensor. Hairpin vortices dominate at the interface between the outer and core regions of the channel, which are not observed in such abundance for moderate to high Reynolds number developed turbulent flows (Schlatter et al. Reference Schlatter, Li, Örlü, Hussain and Henningson2014). Hairpins dominate the flow, with some developing in isolation and others grouping together as hairpin packets. Figures 15 and 16 are therefore associated with ejections of hairpin structures in isolation, and as hairpin packets, respectively. These structures are not a dominant feature of high Reynolds number channel flow but have been detected by Hack & Schmidt (Reference Hack and Schmidt2021) at $\textit {Re}_\tau = 2000$ by conditionally averaging the most dissipative events (99.9th percentile). However, they clearly dominate in this weakly stratified regime. An isolated hairpin vortex is visualised in figure 21, with iso-surfaces of positive and negative $u'$ and $v'$. The hairpin structure is consistent with visualisations in previous (unstratified) simulations and experiments (see e.g. Dennis & Nickels Reference Dennis and Nickels2011; Hack & Moin Reference Hack and Moin2018). Negative $u'$ and positive $v'$ are observed in the vortex core (Q2 event, referring to an event in the second quadrant of the $u' - v'$ plane), and positive $u'$ and negative $v'$ at the edges of the vortex (Q4 event, referring to an event in the fourth quadrant of the $u' - v'$ plane). The low speed upward moving fluid in the centre of the hairpins is a Q2 event, as per typical hairpin structures that dominate low Reynolds number boundary layer flows.

Figure 20. Iso-surfaces of swirl strength, $\lambda _{ci}$, in the lower half of the channel for Case S, coloured by vertical coordinate $y$. The colour map is centred at $y=0.8$, such that ‘red’ regions highlight structures in the core region of the channel. The flow direction is approximately ‘north–west’. The $\lambda _{ci}$ contour value is 0.5 % of its maximum value.

Figure 21. An isolated hairpin vortex visualised by an iso-surface of swirl strength, $\lambda _{ci}$, in the lower half of channel for Case S (grey). Red and blue iso-surfaces represent $u' = 1$ and $u'=-1$, respectively (a), and $v' = 1$ and $v'=-1$, respectively (b). The $\lambda _{ci}$ contour value is 0.5 % of its maximum value.

The Case S hairpins are quantitatively different to those of Case U, as shown in figure 22. Here hairpins have been visualised by $\lambda _{ci}$ along with contours of $u' = -1$, visualising low speed Q2 streaks. Long and wide streaks are associated with Case U which extend from the wall up to the channel centre. Hairpins are not a dominant feature, but some vortical structures can be observed wrapped around the low speed streaks which originate near the wall. Hairpins are the dominant coherent structure of Case S, and are fundamentally different to those of Case U. Case S streaks are much narrower and shorter, and it is unclear if they originate from near the wall. This is consistent with the spatial spectra of figures 4 to 6, which show that stratification suppresses long wide structures, or global modes. Further, Garcia-Villalba & del Alamo (Reference Garcia-Villalba and del Alamo2011) have shown that outer region turbulence scales like the local Obukhov length (3.2). It is therefore unlikely that hairpins originate from the inner region of the flow, and are instead generated in the outer region.

Figure 22. Coherent hairpin structures visualised by $\lambda _{ci}$ with iso-surfaces of $u'=-1$ (green). Cases are U (a) and S (b). The $\lambda _{ci}$ contour value is 0.5 % of its maximum value. Both sets of iso-surfaces, $\lambda _{ci}$ and $u'$, are coloured by the vertical coordinate $y$.

Similar hairpin structures have recently been observed in stratified shear flow without walls (Pham & Sarkar Reference Pham and Sarkar2010, Reference Pham and Sarkar2011; Pham, Sarkar & Winters Reference Pham, Sarkar and Winters2012; Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019). Watanabe et al. (Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019) observed remarkably similar hairpins at the edge of a stably stratified shear layer. The abundance of hairpin structures in the stratified shear layer was explained by reporting the integral shear parameter (Jiménez Reference Jiménez2018), $S^{*}$, defined as the ratio between the eddy turnover time and the time scale of the mean shear

(3.17)\begin{equation} S^{*} = \frac{\overline{u_i'u_i'}}{\overline{\varepsilon}} \frac{\partial {\overline{U}}}{\partial {y}}. \end{equation}

A Value of $S^{*} \gg 1$ implies the evolution of turbulent structures is strongly dependent on mean shear, typically peaking in the buffer region of channel flows and boundary layers before decaying to $S^{*} \approx 10$ in the logarithmic and outer regions of the flow. Watanabe et al. (Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019) noted that hairpins were observed when $S^{*} \gtrsim 10$ at the centre of the stably stratified shear layer. The $y$-dependence of the integral shear parameter $S^{*}$ for cases U and S is plotted in figure 23. Under stable stratification $S^{*}$ peaks at a higher value in the channel core than in the inner region, in contrast to the unstratified case where $S^{*}$ tends to 0 in the region $y \gtrsim 0.4 - 1.0$. The second peak in Case S is due to the local maximum of $\partial _y \overline{U}$ (figure 2) and the low TKE dissipation rate. Further, Case S exhibits a consistently higher $S^{*}$ in the outer region of the flow than Case U, particularly in the region $y > 0.5$, where hairpins dominate. The high value of $S^{*}$ implies coherent turbulent structures are strongly dependent on shear, and can therefore be considered to be quasi-linear (Jiménez Reference Jiménez2018).

Figure 23. Variation with $y$ of the integral shear parameter $S^{*} = \overline {u_i'u_i'} \partial _y \overline{U} /\varepsilon$ for Case U (solid line) and Case S (dashed line).

Hairpins were also observed in stratified shear flow by Pham & Sarkar (Reference Pham and Sarkar2010) and Pham et al. (Reference Pham, Sarkar and Winters2012). They found hairpin vortices (which they referred to as ‘horseshoe’ vortices) appeared to originate as spanwise instabilities of Kelvin–Helmholtz rollers in an unstable shear layer. Spanwise vortices were subsequently drawn into the stable shear layer beneath and stretched into hairpins. However, we did not detect any evidence of Kelvin–Helmholtz instabilities in the present LES. The study of Pham & Sarkar (Reference Pham and Sarkar2011) observed similar hairpins to Pham & Sarkar (Reference Pham and Sarkar2010) and Pham et al. (Reference Pham, Sarkar and Winters2012) but without the presence of Kelvin–Helmholtz rollers. No insight into their origin was given, but it is entirely reasonable that hairpins in stratified channel flow are governed by an analogous process. In particular, the well-mixed outer region can be interpreted as an unstable shear layer, neighbouring the stably stratified channel core. One hypothesis is that the initial spanwise vortices originate in the outer region due to the internal waves in the channel core, noting that the waves can be detected in spatial spectra at $y=0.75$ (figure 5). Alternatively (or additionally) spanwise vortices may be generated locally in the outer region through quasi-linear processes, due to high $S^{*}$. Regardless of their origin, outer layer spanwise vortices are stretched and lifted by background shear and ejected into the channel core where they dissipate and perturb internal waves.