The coupled dynamics of internal waves and hairpin vortices in stratified plane Poiseuille flow

Abstract

A simulation of stably stratified plane Poiseuille flow at a moderate Reynolds number ($\textit {Re}_\tau = 550$) and Richardson number ($\textit {Ri}_\tau = 480$) is presented. For the first time, the dynamics in the channel core are shown to be described as a series of internal waves that approximately obey a linear wave dispersion relationship. For a given streamwise wavenumber $k_x$ there are two internal wave solutions, a dominant low frequency mode and a weaker-amplitude high-frequency mode, respectively corresponding to ‘backward’ and ‘forward’ propagating internal waves relative to the mean flow. Analysis of linearised equations shows that the dominant low-frequency mode appears to arise due to a particularly sensitive response of the mean flow profiles to incoherent forcing. Instantaneous visualisations reveal that hairpin vortices dominate the outer region of the channel flow, neighbouring the buoyancy dominated channel core. These hairpins are fundamentally different from those observed in canonical unstratified boundary layer flows, as they arise via quasi-linear local processes far from the wall, governed by background shear. Outer region ejection events are common and can be induced by high amplitude waves. Ejected hairpins are transported into the channel core, in turn ‘ringing’ the prevailing strong buoyancy gradient and thus generating high-amplitude internal waves, high dissipation and wave breaking, induced by spanwise vortex stretching and baroclinic vorticity generation. Such spontaneous and sustained generation of quasi-linear internal waves by wall-bounded sheared turbulence may provide novel idealised solutions for, and insight into, large-scale turbulent mixing in a wide range of environmental and industrial flows.

1. Introduction

Stably stratified shear flows are commonplace in nature, dominating fluid dynamical processes in the ocean (Wells & Dorrell 2021) and atmosphere (Mahrt 2014). The restoring forces of buoyancy play a key role in mixing processes and have been subject to considerable research over the last few decades (Caulfield 2021). Dynamics are rich even under weak (stable) stratification and can induce strong anisotropy, intermittency, layering, and internal waves (Caulfield 2021). These dynamics are crucial for the understanding, and parameterisation in larger scale models, of processes fundamental to natural flows, including: transport of scalar properties, e.g. temperature and salinity (Garaud 2018) and particulate concentration (Hung, Niu & Chou 2020); mixing and entrainment of ambient fluids (Wells, Cenedese & Caulfield 2010); and energy transport (Winters et al. 1995).

To understand the complex dynamics of stratified turbulent flow this paper focuses on a particular idealised wall-bounded flow, namely stratified plane Poiseuille flow, also referred to as channel flow. Simulations of stratified channel flow have been the focus of much research over the last two decades (Garg et al. 2000; Armenio & Sarkar 2002; Iida, Kasagi & Nagano 2002; Moestam & Davidson 2005; Garcia-Villalba & del Alamo 2011; Zonta, Onorato & Soldati 2012; Zonta 2013; Zonta & Soldati 2018). Flows have typically been modelled using direct numerical simulation (DNS) or large-eddy simulation (LES). However, there has been limited quantitative analysis of the spatio-temporal structure of these flows.

Stratified channel flow is governed by three dimensionless quantities: the friction Reynolds number ($\textit {Re}_\tau$); the friction Richardson number ($\textit {Ri}_\tau$); and the Prandtl number ($\textit {Pr}$), defined as

(1.1a–c)\begin{equation} \textit{Re}_\tau = \frac{u_\tau \delta}{\nu}, \quad \textit{Ri}_\tau = \frac{\delta g \Delta \rho}{ \rho_0 u_\tau^{2}}, \quad \textit{Pr} = \frac{\nu}{\kappa}, \end{equation}

where $u_\tau = \sqrt {\tau _w/\rho _0}$ is the friction velocity, $\tau _w$ the wall shear stress, $\rho _0$ a reference density, $\delta$ the channel half-height, $g$ the gravitational acceleration, $\Delta \rho$ the density difference between the upper and lower walls, $\nu$ the kinematic viscosity and $\kappa$ the diffusivity of the density field. For constant Reynolds and Prandtl numbers, Garcia-Villalba & del Alamo (2011) demonstrated the existence of three flow regimes, which they referred to as strongly stratified laminar flow, strongly stratified turbulent flow and weakly stratified turbulent flow. Each regime describes the state of turbulence, which is still stratified with depth (Zonta & Soldati 2018). At sufficiently high Richardson numbers, turbulence is entirely suppressed resulting in the strongly stratified laminar flow regime. As the Richardson number reduces, local regions of intermittent turbulence can be sustained near the walls, leading to the strongly stratified turbulent flow regime. Critical values of $\textit {Ri}_\tau$ for the transition from laminar flow can be well predicted through linear stability analysis (Gage & Reid 1968). The transition to weakly stratified turbulence is (inevitably) less well defined, but is quantified when near-wall turbulence is no longer intermittent. Here, the core of the channel is strongly affected by buoyancy while the near-wall dynamical behaviour is dominated by turbulence. In this regime there is strong spatial variability, which as we discuss below, leads to subtle interplay between (relatively) localised internal waves and disordered turbulent motions.

This paper is concerned with this interesting case of weakly stratified flow, where both turbulence and buoyancy control the flow dynamics in complex and interconnected ways. Heretofore, the most extensive data set for weakly stratified turbulent channel flow was created by Garcia-Villalba & del Alamo (2011), using DNS at $\textit {Pr} = 0.7$, $\textit {Re}_\tau = 180$ and $\textit {Re}_\tau = 550$, and a Richardson number range of up to $\textit {Ri}_\tau = 1920$. The high Reynolds number allowed Garcia-Villalba & del Alamo (2011) to simulate a wider range of weakly stratified cases compared with previous authors who were limited by $\textit {Re}_\tau =180$, where the strongly stratified regime onsets at lower $\textit {Ri}_\tau \approx 60$ (Garcia-Villalba & del Alamo 2011). In contrast, simulations at $\textit {Re}_\tau = 550$ and $\textit {Ri}_\tau = 960$ were still found to be ‘weakly stratified’, in this specific sense.

In the weakly stratified flow regime the dynamical influence of buoyancy varies qualitatively in the vertical (wall-normal) direction. Near the wall, effects of buoyancy are relatively small and the flow is turbulent, while at the centre of the channel there is evidence that relatively coherent internal waves are present (Armenio & Sarkar 2002; Iida et al. 2002; Moestam & Davidson 2005; Garcia-Villalba & del Alamo 2011). Moestam & Davidson (2005) related periodic variations in stratified shear flow to the (local) dispersion relation of linear internal waves, Doppler shifted in a spatially varying streamwise flow

(1.2)\begin{equation} (\omega - U_c k_x)^{2} = N^{2} \left(\frac{k_x + k_z}{k_x + k_y + k_z}\right)^{2}, \end{equation}

where $\omega$ is the temporal frequency, $U_c$ the (time- and planar-averaged) velocity of the streamwise background flow and $k_i$ the three-dimensional wavenumber. Here, $N$ is the buoyancy frequency associated with the time- and planar-averaged density, which, like $U_c$, varies in the vertical direction. However, Moestam & Davidson (2005) neglected the effect of Doppler shift ($U_c k_x$), assumed (in a WKB fashion) that the buoyancy frequency varies slowly relative to the vertical wavelength of the internal waves, and approximated spatial wavenumbers using the domain size. These are sweeping assumptions, particularly in the channel core, where internal waves have been observed over a range of wavenumbers (Iida et al. 2002; Garcia-Villalba & del Alamo 2011), the mean velocity is large, and the buoyancy frequency varies rapidly.

Iida et al. (2002) and Garcia-Villalba & del Alamo (2011) characterised internal waves by investigating spatial energy spectra. Both studies found peaks in vertical velocity energy spectra over a range of wavelengths $k_x \approx 1$ to 4. In addition, both studies reported a phase shift of ${\rm \pi} /2$ between the density and vertical velocity signals, consistent with solutions to the linearised energy equation for internal waves. The higher Reynolds number simulations of Garcia-Villalba & del Alamo (2011) allowed separation of scales of motion in spatial spectra; providing clear evidence that internal waves were confined to the channel core at $\textit {Re}_\tau = 550$, while Garcia-Villalba & del Alamo (2011) reported that confinement of internal waves at $\textit {Re}_\tau = 180$ was difficult to determine due to the smaller scale separation between the outer and inner regions.

There are, however, several important unanswered questions regarding the nature of the internal wave field in weakly stratified channel flow. While there is some evidence that these waves are internal waves, it is fair to say that an appropriate form for their dispersion relation has yet to be clearly determined, and in particular the relevance of linearised descriptions to what are clearly finite-amplitude structures is uncertain. In addition, the mechanisms both generating and sustaining waves in stratified channel flow have not yet been quantified, nor has their interaction with near-wall turbulence.

This paper aims to answer these questions by presenting a simulation of the weakly stratified channel flow at moderate Reynolds and Richardson numbers ($\textit {Re}_\tau = 550$, $\textit {Ri}_\tau = 480$). Numerical methodology is detailed in § 2 before analysing time- and planar-averaged statistics, multi-dimensional spectra and dynamic mode decomposition (DMD) in §§ 3.1 to 3.3, where highly coherent internal waves are found to dominate the core of the channel in the region $0.8 < y < 1.2$, where $y$ is the vertical wall-normal coordinate scaled by the channel half-height. The dominant mode identified is found to be a ‘backward’ travelling (relative to the spatially and temporally averaged flow velocity in the channel core) internal wave approximately satisfying a linear dispersion relation, $\omega = \overline{U}_{mean} k_x - N_{mean}$, where $\overline{U}$ is the time- and planar-averaged streamwise velocity and subscript ‘mean’ denotes spatial averaging over the full channel core region ($0.8 < y < 1.2$). Solutions to the linearised governing equations are also sought and tested for relevance, formulated as both a stochastically forced initial value problem and a differential eigenvalue problem in §§ 3.4 and 3.5. Internal waves are found to emerge as a sensitive filtered response to turbulent perturbations at the core edge over a range of wavenumbers and frequencies. Analysis of the instantaneous structure of plane Poiseuille flow (§ 3.6) reveals that the largest turbulent perturbations to the core edge are coherent hairpin vortices. Crucially, these vortices are generated via local processes in the buoyancy-affected outer region of the channel. These dynamics are clearly qualitatively different from canonical boundary layer flows, where hairpin vortices are generated in the inner region much closer to the wall. These (stratified) hairpins are ejected from the turbulent region of the channel and into the buoyancy-dominated channel core, ‘ringing’ the relatively strong density gradient and thus both inducing periodic internal waves and triggering local wave breaking. Results are discussed in § 4 and placed in context, before we conclude with some possible future avenues of research inspired by these findings.

2. Methodology

A fully developed stably stratified channel flow is numerically simulated. The flow is driven past no-slip walls by a constant pressure gradient (constant $\textit {Re}_\tau$), and stable stratification is maintained by imposing constant upper and lower wall densities. The dimensionless continuity, momentum and density scalar transport equations are solved:

(2.1)$$\begin{gather} \frac{\partial {U_i}}{\partial {x_i}} = 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\textrm{D} {U_i}}{\textrm{D} t} ={-} \frac{\partial {P}}{\partial {x_i}} + \frac{1}{\textit{Re}_\tau} \frac{\partial^{2} {U_i}}{\partial {x_j}^{2}} + f_i - \textit{Ri}_\tau \rho' \,\hat{y}_i, \end{gather}$$

and

(2.3)\begin{equation} \frac{\textrm{D} {\rho}}{\textrm{D} t} = \frac{1}{\textit{Re}_\tau \textit{Pr}} \frac{\partial^{2} {\rho}}{\partial {x_j}^{2}}, \end{equation}

where $U_i = (U,V,W)'$ represents velocity, $x_i = (x,y,z)'$ represents Cartesian spatial coordinates, $t$ represents time, $P$ represents kinematic pressure, $\rho$ is the density field and $\hat {y}_i$ is the $y$-direction unit vector. The flow is forced by a constant streamwise (negative) pressure gradient, $f_i = (1,0,0)'$. Flow variables are made dimensionless by the channel half-height $\delta$, friction velocity $u_\tau$ and the density difference between the upper and lower walls, $\Delta \rho$. Resultant dimensionless variables are the friction Reynolds number ($\textit {Re}_\tau$), the friction Richardson number ($\textit {Ri}_\tau$) and the Prandtl number ($\textit {Pr}$). The buoyancy force $\textit {Ri}_\tau \rho ' \,\hat {y}_i$ is dependent on the density fluctuation, $\rho ' = \rho - \bar {\rho }$ where $\bar {\rho }$ is the time (averaged over time $T$ in table 1) and planar-averaged (averaged in homogeneous directions $x$ and $z$) density field. In this case the purely wall-normal $y$-dependent density field ($\bar {\rho }$) is absorbed (hydrostatically) into the pressure term. In practice $\bar {\rho }$ is approximated at each timestep by the instantaneous planar-averaged density, consistent with previous work (e.g. Armenio & Sarkar 2002; Garcia-Villalba & del Alamo 2011). However, it should be noted that preliminary simulations carried out with a buoyancy term based upon the full density field $\rho$ led to no difference in turbulent statistics, aside from the pressure term.

Table 1. Cases with time-averaged parameters and case-specific quantities. The bulk Reynolds and Richardson numbers are given by $\textit {Re}_b = U_b \delta /\nu$ and $\textit {Ri}_b = \Delta \rho g \delta /2 \rho _0 U_b^{2}$ where $U_b$ is the bulk velocity. $Nu = -2 \partial _y \bar {\rho }|_w$ is the Nusselt number and $\varLambda = 2 \textit {Re}_\tau \textit {Pr} / \varkappa \textit {Ri}_\tau Nu$ is the Obukhov length normalised by the channel half-height $\delta$ with the Kármán constant $\varkappa = 0.41$; $\varDelta ^{+}_i$ quantifies the element sizes in wall units (normalisation by $\delta _\nu =\nu /u_\tau$). Note that each element is further discretised by $8^{3}$ GLL points.

The domain is presented in figure 1. The reference density is taken as the mean density in the channel, equal to the density at the channel centreline, such that the (dimensionless) density boundary conditions at the walls are $\rho (y=0) = \tfrac {1}{2}$ and $\rho (y=2) = -\tfrac {1}{2}$. Velocity boundary conditions are no slip at the walls, and periodic conditions are applied at the $x$- and $z$-normal boundaries.

Figure 1. Numerical domain. Coloured regions represent the inner region at $y\leqslant 0.2$ and $y\geqslant 1.8$ (blue), the outer region at $0.2 < y \leqslant 0.8$ and $1.2 \leqslant y < 1.8$ (pink) and the channel core at $0.8 < y < 1.2$ (red). The dashed line represents the channel centreline at $y=1$.

Equations are discretised and solved using Nek5000, a spectral element code developed by Fischer et al. (2008). The dimensionless domain size is $L_x \times L_y \times L_z = 8 {\rm \pi}\times 2 \times 3 {\rm \pi}$, consistent with Garcia-Villalba & del Alamo (2011). The computational domain is decomposed into $N_x \times N_y \times N_z$ hexahedral spectral elements on which the governing equations are discretised using a Galerkin method. Vertical element sizes are determined by a hyperbolic stretching function, with a maximum element size at the centre of the channel. Equations are solved by means of local approximations based on high-order tensor-product polynomial basis on Gauss-Lobatto-Legendre (GLL) nodes. Simulations adopt $8^{3}$ GLL nodes in each element (seventh-order polynomials) and equations are solved using the $\mathbb {P}_n\text {-}\mathbb {P}_n$ scheme, where all variables are represented by the same polynomial order (Tomboulides, Lee & Orszag 1997). Temporal discretisation uses third-order backward differencing with a dimensionless timestep of $1\times 10^{-4}$, and nonlinear terms are dealiased using the $3/2$ rule (Orszag 1979; Canuto et al. 2012). Grid resolution is detailed in table 1.

The effects of sub-grid-scale dissipation are accounted for via modal based explicit filtering, equivalent to deconvolution LES or hyper-viscosity (Stolz, Schlatter & Kleiser 2005). This method solves the unfiltered governing equations (2.1) to (2.3) and applies an explicit low-pass filter at the end of each time step to dissipate energy artificially at the highest polynomial modes. These simulations adopt an attenuation amplitude of 0.05 following previous work modelling wall-bounded turbulent flows (Lai, Merzari & Hassan 2019; Merzari et al. 2020a,b; Jin et al. 2021). For further details see Fischer & Mullen (2001) and Chatterjee & Peet (2018). While simulations are coarser than DNS in the spanwise and streamwise directions, the stratified flow case is fully resolved in the vertical direction (see e.g. Garcia-Villalba & del Alamo 2011) such that artificial dissipation is primarily felt in the homogeneous directions. The explicit filtering methodology is validated against the stratified channel flow DNS of Garcia-Villalba & del Alamo (2011) in Appendix A, where a grid sensitivity study is also presented. Despite the close agreement between the LES herein and the DNS of Garcia-Villalba & del Alamo (2011) it should be noted that all presented data are of filtered quantities and therefore should be treated as approximate values, particularly variables derived from the dissipation rates.

Two simulations are presented in this paper, detailed in table 1 ($\textit {Pr} = 1$ for both simulations). Case U, of an unstratified flow with a passive scalar field, is initialised with a perturbed logarithmic velocity profile while Case S, of a stratified flow, is initialised from pseudo-steady U data, effectively by instantaneously ‘turning on’ gravity, and so making the passive scalar field (density) dynamically significant which subsequently feeds back on the velocity field via the buoyancy force. Before data collection, all cases are advanced to a statistically steady state in two steps. First a coarse simulation is run with a polynomial order of 5 ($6^{3}$ GLL points per element). Secondly, the simulation is run with the desired polynomial order of 7 ($8^{3}$ GLL points). Convergence of both steps is assessed by monitoring the friction Reynolds numbers at the walls and the bulk Richardson number. This process reduces overall computational cost while ensuring the flow is fully developed before data collection. Subsequent data collection is carried out over a dimensionless time interval of $T$, which varied between the two cases according to table 1.

Over integration period $T$ time-averaged statistics are calculated at all GLL points which are subsequently plane averaged in the homogeneous directions ($x$ and $z$). In order to assess temporal dependence of flow variables time-series data are collected at several slices through the domain, several normal to $y$, and a slice normal to $z$ at $z=L_z/2$. Slice data are collected on uniform grids at a resolution equal to the total number of degrees of freedom in each direction, obtained using spectral interpolation. Slice data containing instantaneous velocity ($U,V,W$), pressure ($P$), density ($\rho$), velocity gradients ($\partial U_i /\partial x_j$) and density gradients ($\partial \rho / \partial x_j$) are collected every 20 timesteps (every $2\times 10^{-3}$ dimensionless time units) over duration $T$, resulting in a time series of at least 20 000 slices. The integration time ($T$) is over 4 times longer than previous studies (e.g. Garcia-Villalba & del Alamo 2011), enabling convergence of temporal spectra and DMD data. This paper will demonstrate that intermittent ejection events near the channel core occur over ${O}(\delta /u_\tau )$ time scales, necessitating the longer integration time used herein.

3. Results

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 2000), 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 (2011).

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 2000). 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 (2011). 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 (2011). 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 2003). 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 (2011) 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 (2011) 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 (2011) 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 1984)

(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 (2011) 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 (2011).

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 (2003), where waves at such an angle were observed propagating away from a turbulent mixed region forced by an oscillating grid. Dohan & Sutherland (2003) 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 2014) 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 2014) ‘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 2010) 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 (2010): 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 (2011).

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 (2003), 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 1993; Jovanović & Bamieh 2005; Hwang & Cossu 2010; Jovanović 2021), with previous work often adopting an approach based on linear dynamical systems analysis (see 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 (2012) 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. 2020). 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. 2012; Lian, Smyth & Liu 2020) 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. (2020), 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. (2020) 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 2021). 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. 2014). 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 (2021) 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 2011; Hack & Moin 2018). 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 (2011) 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 2010, 2011; Pham, Sarkar & Winters 2012; Watanabe et al. 2019). Watanabe et al. (2019) 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 2018), $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. (2019) 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 2018).

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 (2010) and Pham et al. (2012). 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 (2011) observed similar hairpins to Pham & Sarkar (2010) and Pham et al. (2012) 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.

4. Discussion

It has been shown that stably stratified channel flow can be split into three regions: the inner region ($y < 0.2$); the outer region ($0.2 < y < 0.8$); and the channel core ($y > 0.8$), all symmetric about the channel centreline, $y=1$. These definitions are consistent with canonical unstratified channel flow and boundary layer definitions, apart from the buoyancy-affected channel core. This paper has demonstrated that a series of internal waves dominate the channel core, approximately following a linear dispersion relation with constant flow velocity and buoyancy frequency $\omega = \pm {N}_{mean} + \overline{U}_{mean} k_x$, where subscript mean denotes averaging of variables over the full channel core height. Waves are generated over a wide range of wavenumbers and frequencies, with a dominant low-frequency ‘backward’ propagating mode, and a low-amplitude high-frequency ‘forward’ propagating mode, relative to the background mean flow. Spontaneous development of dominant (at least apparently) linear waves in stratified flows may offer insight into the large-scale mixing dynamics of environmental flows (Caulfield 2021; Wells & Dorrell 2021). Analysis of linearised equations have shown that the dominant relatively low-frequency mode can be interpreted as arising either due to a particularly sensitive response of the mean flow profiles to incoherent forcing, or alternatively due to the exponential growth of instabilities of the mean flow profiles themselves. Either (linear) interpretation implies that perturbations can ring the system and generate highly coherent structures spanning the full height of the channel core. However, inherently linear mechanisms do not appear to reproduce the higher-frequency, yet lower-amplitude, internal waves which are also observed in the full nonlinear numerical simulations.

Instantaneous (three-dimensional) data reveal intermittent spontaneous highly dissipative events, where turbulent structures are ejected into the channel core from the turbulent outer region. These extreme structures generate high-amplitude waves and drive local wave breaking as they pass from the outer region and into the core. These structures are largely coherent hairpin vortices which dominate the outer region, arising as isolated structures or hairpin packets, and are fundamentally different to those typically observed in canonical boundary layer and channel flows. In particular, these stratified hairpins arise via local processes rather than originating near the wall. The low speed streaks at the centre of these hairpins are much thinner and shorter than in unstratified channel flow, and Garcia-Villalba & del Alamo (2011) have shown that spectra in this region scales like the local Obukhov length (Nieuwstadt 1984). Similar structures have been observed in stratified shear flow (e.g. Pham et al. 2012; Watanabe et al. 2019), in the absence of walls, and can be partially explained through consideration of the shear parameter $S^{*}$, as defined in (3.17), and as the ratio of the eddy turnover time to the time scale of mean shear. $S^{*}$ has a peak at the edge of the channel core, larger than its peak in the buffer layer for Case S, indicating that outer region turbulent structures are quasi-linear and strongly affected by mean shear, consistent with typical hairpin vortices in the near-wall region of unstratified flow. It is hypothesised that these vortices evolve due to similar processes to those observed by Pham & Sarkar (2011), where spanwise vortex tubes originate in the outer region, analogously to in an unstable shear layer. The origin of these spanwise vortices is hypothesised to be due to local quasi-linear processes, or alternatively as a result of internal waves in the channel core. The evolution of these structures is governed by the background shear which stretches and rolls vortices into hairpins. Hairpins are ejected from the outer region into the channel core, causing large-scale coherent internal waves. As hairpins pass through the critical layer they cause wave breaking, due to spanwise vortex stretching and baroclinic vorticity generation.

Analysis of the linearised system has shown that the high-frequency low-amplitude modes observed in the 3-D simulations must arise from either deviations from the mean flow profiles or nonlinearity. A reasonable hypothesis is that hairpin ejections cause strong deviations to the mean flow profiles that may allow the high-frequency mode to propagate. However, conditional averaging over particularly high $\varepsilon$ events in the channel core does not produce deviations large enough to alter the dispersion relation obtained from the viscous Taylor–Goldstein equations, and only has a minor distorting effect on the modal spatial structure (not shown). It is therefore concluded that the high-frequency mode must arise from the strong and nonlinear ejections of hairpin vortices.

The instantaneous vertical velocity profiles of figure 17 actually demonstrate an intricate coupling between the hairpin ejections and the channel core, where ejections typically occur in phase with the internal waves. High-amplitude internal waves therefore enhance ejection events, allowing them to penetrate deeper into the core. It is plausible that a self-sustaining process is present, where hairpin ejections lead to high-amplitude waves which subsequently draw in further hairpins from the outer region. In addition, high-amplitude internal waves may sustain further hairpins by inducing spanwise vorticity in the outer layer, subsequently leading to the rolling and stretching of vortices into hairpins. This cycle is not necessary for waves to persist, but does explain the occasional, highly active periods in the flow (e.g. figure 19 at $t \approx 15$).

An immediate question that arises from this study is how these mechanisms and structures scale with Reynolds, Richardson and Prandtl numbers. Some insight may be gained from the study of Garcia-Villalba & del Alamo (2011), who simulated the stratified channel flow at $\textit {Ri}_\tau = 960$. They obtained very similar spatial spectra and density/velocity profiles as in their case with $\textit {Ri}_\tau = 480$ and the present study. Given the primary backward travelling waves appear to arise due to instability of the mean flow profiles, the waves reported by Garcia-Villalba & del Alamo (2011) are likely governed by the same processes identified here.

The influence of $\textit {Ri}_\tau$ on hairpins and subsequent forward travelling modes is less clear. As $\textit {Ri}_\tau$ increases, the velocity profile becomes more jet like (Garcia-Villalba & del Alamo 2011). Subsequently, one may expect the outer region to sustain more hairpin vortices due to an increase in $S^{*}$, which will ring the sharp density profile more frequently. However, an increase in $\textit {Ri}_\tau$ will also act to stabilise the channel core and may therefore suppress strong ejection events. It is also unclear if the same structures and mechanisms are present in the strongly stratified regime, where near-wall turbulence is intermittent (Garcia-Villalba & del Alamo 2011). Clearly these dynamics are complex and should be considered for future research.

5. Conclusions

In conclusion, the present investigation uses numerical simulations to quantify the internal waves at the core ($0.8 < y < 1.2$) of relatively weakly stratified plane Poiseuille flow. Spatio-temporal spectra and DMD have demonstrated, for the first time, that coherent structures spanning the full channel core height largely follow the linear dispersion relation for internal waves (with constant buoyancy frequency and background streamwise velocity) over a wide range of wavenumbers and frequencies. This particular demonstration of the relevance of linear stability to turbulent stratified shear flows opens up pathways of future research of natural flows. The dominant mode is a relatively low-frequency internal wave-like structure, ‘backwards’ propagating relative to the mean flow, which can be interpreted as arising from a sensitive forced response or from a linear instability inherent in the mean flow profiles, demonstrated by solving the stochastically forced linearised governing equations and the viscous Taylor–Goldstein equations respectively. Heuristically, waves may be thought of as being generated via ‘ringing’ of the enhanced density gradient at the edge of the central core region due to perturbations in the outer, turbulent, regions of the flow.

Flow visualisation demonstrates that the resulting high-amplitude internal waves are a result of intermittent highly dissipative ejections from the outer regions and into the channel core. These ejections are associated with coherent hairpin vortices that dominate the boundary between the buoyancy-affected core region and the turbulent outer region. These stratified hairpins are fundamentally different from the attached eddies of canonical boundary layer flows, and arise from local quasi-linear processes where spanwise vortices are stretched and lifted by background shear. Indeed, it appears that the dynamics of the internal wave-like structures and the stratified hairpins are strongly coupled. Specifically, if internal waves are present they can induce further ejections which act in turn in phase with the waves. As hairpins pass through the boundary of the central core region they cause local wave-breaking due to vortex stretching and baroclinic vorticity generation.

This paper quantifies the interaction between outer region turbulence and internal waves in stably stratified plane Poiseuille flow. However, it raises several unanswered questions. For example, it is still unclear how the observed coherent and robust structures scale with key dimensionless parameters such as $\textit {Re}_\tau$, $\textit {Ri}_\tau$ and $\textit {Pr}$. Perhaps most significantly, what role do coherent hairpins and internal waves play in the irreversible turbulent mixing within this flow, and how does their subtle interplay affect the energetic partitioning (i.e. the ‘efficiency’) of the mixing in wall-bounded stratified turbulent flow.