3.1. Qualitative behaviour of the flow structure

We look first at the qualitative structure of the flow, focusing in particular on the instantaneous temperature distribution $\theta$ on a ($y-z$) cross-section of the channel located at $x=L_x/2$. Results, which are shown in figure 2 for $Ri_{\tau }=0$ (panel a, neutrally buoyant case) and for $Ri_{\tau }=300$ (panel b, stably stratified case), will be conveniently discussed by keeping the neutrally buoyant case ($Ri_{\tau }=0$, figure 2a) as reference. For $Ri_{\tau }=0$, temperature is a passive scalar and, as such, it is purely transported by velocity. Vortical structures rising from the boundaries are therefore free to travel over long distances, since they are only bounded by the physical constraint imposed by the walls. Under the action of these vortical structures, a fluid particle with a given temperature is brought to a region with a different temperature, where it thermalizes by diffusion (see figure 2a). This naturally gives a high degree of mixing. For $Ri_{\tau }=300$, the situation is different. While vortical structures are still dominating the near-wall region, their influence away from the boundary appears limited. The reason is that vortical structures are in this case subject to an additional vertical (i.e. in the wall-normal direction) constraint imposed by gravity, which reduces their range of influence to approximately half of the channel height (this is rather clear upon comparison of figure 2a,b). As a consequence, the channel results divided into a top, hot part, and a bottom, cold part (figure 2b). These two parts behave almost independently each other, and are separated by IGWs (streaky structure developing at the channel centre). The physics of IGWs is simple: because of the background density profile – with density decreasing with height – a fluid particle that is displaced in the wall-normal direction by velocity fluctuations is subject to a restoring buoyancy force that tends to bring it back to its initial position. The fluid particle trespasses its initial equilibrium position and overshoots inertially, giving rise to an oscillation that constitutes the essence of IGWs. We anticipate here, but it will become clear by looking at the fluid statistics in the next paragraphs, that IGWs appear inside a thermocline – i.e. a region where temperature changes much more than it does above and below it, hence representing a sort of thick interface that hinders the wall-normal transfer of momentum and heat. Due to their importance in the dynamics of stably stratified flows, IGWs have been analysed in detail in a number of previous studies (we refer the reader to Staquet & Sommeria (Reference Staquet and Sommeria2002), for a comprehensive review on the topic). IGWs will be further characterized from a more quantitative and qualitative viewpoint in § 3.4.

Figure 2. Contour maps of temperature, $\theta$, on a $(y\unicode{x2013}z)$ cross-section located at $x=L_x/2$. (a) Refers to the neutrally buoyant case, $Ri_{\tau }=0$ and (b) refers to the stably stratified case at $Ri_{\tau }=300$.

To appreciate the different flow structure induced by the stable stratification, we now turn our attention to the distribution of temperature $\theta$ and axial velocity $u_x$ on a longitudinal $(x\unicode{x2013}z)$ plane located at $y=L_y/2$. Results are shown in figure 3: panel (a) and panel (b) refer to $Ri_{\tau }=0$ (temperature and axial velocity, respectively), while panel (c) and panel (d) refer to $Ri_{\tau }=300$ (temperature and axial velocity, respectively). The flow moves from left to right. As expected, for $Ri_{\tau }=0$, temperature (figure 3a) – which is a passive scalar – is efficiently mixed throughout the entire height of the channel by the dominant vortical structures. The flow (figure 3b) appears organized into taller vortices which are emitted from the wall and contain ensembles of smaller-scale vortices. The wall-normal extension of these taller vortices scales with the channel height ${\sim }2h$. For $Ri_{\tau }=300$, on the other hand, the situation is controlled by the presence of IGWs, which – as clearly visible in both the temperature and velocity maps – dominate the central region of the channel. The presence of these structures naturally modifies the entire dynamics of the flow, inducing an extra confinement to the wall-normal development of vortices, and reducing at the same time their capability of effectively mixing the flow. As a side observation, we note that the temperature field in the proximity of the channel centre appears stretched and tilted at an angle with respect to the horizontal direction. This is due to the presence of a strong vertical shear in that region.

Figure 3. Contour maps of temperature $\theta$ (a,c) and streamwise velocity $u_x$ (b,d) on a $(x\unicode{x2013}z)$ longitudinal section located at $y=L_y/2$. (a,b) Refer to the neutrally buoyant case, $Ri_{\tau }=0$ and (c,d) refer to the stably stratified case at $Ri_{\tau }=300$. The temperature tilting induced by the vertical shear at the channel centre is also explicitly indicated (c).

3.2. Velocity and temperature statistics

We now characterize the flow from a statistical viewpoint. Unless differently stated, all results will be presented in wall units (denoted by superscript $+$), obtained normalizing velocities by $u_{\tau }$, lengths by $l_{\tau }=\nu /u_{\tau }$, times by $t_{\tau }=\nu /u_{\tau }^{2}$ and temperatures by $\theta _{\tau }=q_w/u_{\tau }$. In figure 4 we show the behaviour of the mean streamwise velocity $\langle u_x^{+} \rangle$ as a function of the wall-normal coordinate, in linear (figure 4a) and semilog (figure 4b) scale. Brackets indicate time and space average over the homogeneous directions. Results are rendered according to the following colour code: blue (up-triangle) refers to the neutrally buoyant case ($Ri_{\tau }=0$), yellow (down-triangle) refers to $Ri_{\tau }=50$, green (square) refers to $Ri_{\tau }=100$, purple (diamond) refers to $Ri_{\tau }=200$ and red (circle) refers to $Ri_{\tau }=300$. The law of the wall, $\langle u_x^{+} \rangle =z^{+}$, and $\langle u_x^{+} \rangle =(1/\kappa ) \log (z^{+})+5.5$, with $\kappa$ the von Kármán constant, is also shown in figure 4(b) by a solid line. As expected (see in particular figure 4b), in the neutrally buoyant case the mean velocity closely follows the law of the wall, since temperature is a passive scalar that does not influence the velocity field. In the stably stratified cases, we observe an increase of the mean axial velocity, which is particularly pronounced in the core part of the channel (see both figure 4a,b), as the shape of the velocity profile deviates significantly from the classical logarithmic behaviour and approaches a nearly laminar, parabolic shape (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011). This tendency towards a laminarization in the core part of channel is due to the conversion of turbulent kinetic into potential energy, which occurs when a fluid particle moves in the wall-normal direction within the flow. Note also that, since the mean pressure gradient is kept constant among the different simulations, the mean wall stress remains constant, and so does the slope of the velocity at the wall, with all profiles collapsing onto that of the neutrally buoyant case ($z^{+}/h^{+}<0.1$).

Figure 4. Mean fluid streamwise velocity $\langle u^{+}_x \rangle$ as a function of the wall-normal direction, $z^{+}$, in linear (a) and semilog scale (b) for the different cases considered in the present study. Comparison between the reference case of unstratified turbulence ($Ri_{\tau }=0$), and the stratified turbulence at $Ri_{\tau }=50$, $Ri_{\tau }=100$, $Ri_{\tau }=200$ and $Ri_{\tau }=300$ (filled symbols). The classical law of the wall $\langle u_x^{+} \rangle =z^{+}$ and $\langle u_x^{+} \rangle =(1/\kappa ) \log (z^{+})+5.5$, with $\kappa$ the von Kármán constant, is also shown for comparison in (b) (solid line).

As shown in figure 5, stratification modifies the behaviour of the mean temperature field $\langle \theta \rangle$, which takes a layered structure formed by a near-wall layer ($0< z^{+}/h^{+}<0.1$), a transition layer ($0.1< z^{+}/h^{+}<0.8$) and a core layer ($0.8< z^{+}/h^{+}<1$). Note that temperature in figure 5(a) is shown in outer units, i.e. not rescaled in wall units. Compared with the neutrally buoyant case, stratification reduces the mean temperature gradient in the near-wall layer (i.e. it reduces the Nusselt number), while at the same time increasing it in the core layer. This latter observation is of particular importance, since it indicates the tendency for turbulent stratified channel flows to develop a kind of thick interface – the thermocline – inside which the temperature changes more vigorously than it does immediately above and below. The thermocline, which forms right where the mean shear vanishes, constitutes a barrier for wall-normal momentum and heat transport. Reportedly (Ferziger, Koseff & Monismith Reference Ferziger, Koseff and Monismith2002; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Zonta et al. Reference Zonta, Onorato and Soldati2012b; Zonta & Soldati Reference Zonta and Soldati2018), there is a strong connection between the presence of a thermocline and the presence of IGWs (which, as discussed in § 3.1, occur at the channel centre), in the sense that IGWs develop where a thermocline exists. Interestingly, the temperature gradient in the core region of the channel does not increase monotonically for increasing $Ri_\tau$: it first increases (going from $Ri_\tau =0$ to $Ri_\tau =100$), and then reduces (going from $Ri_\tau =100$ to $Ri_\tau =300$). We anticipate here, but we will come back to it later, that this non-monotonic trend of the temperature profile at the channel centre can be explained by looking at the intertwined behaviour of the turbulent (or buoyancy) and diffusive fluxes. In the transition layer, between the near-wall and the core layers, the temperature gradient remains small, as the flow is characterized by a higher degree of mixing. Not surprisingly, when rescaled in wall units, i.e. by keeping the friction temperature $\theta _\tau =q_w/(\rho c_p u_\tau )$ as reference, the mean temperature profile $\varTheta ^{+}=(\theta -\theta _w)/\theta _\tau$ recovers a monotonic behaviour (figure 5b). This is particularly visible in the inset of figure 5(b), where a close-up view of the mean temperature profile in the core region of the channel is given. Note also that, for the neutrally buoyant case, and similarly to what happens for the streamwise velocity $\langle u^{+}_x \rangle$, the behaviour of $\varTheta ^{+}$ can be parameterized by $\varTheta ^{+}=Pr \cdot z^{+}$, in the boundary layer, and by $\varTheta ^{+}=2.3 \log (z^{+})+3.11$, in the core region (see Alcántara-Ávila et al. Reference Alcántara-Ávila, Hoyas and Jezabel Pérez-Quiles2021).

Figure 5. Mean fluid temperature $\langle \theta \rangle$ in linear (a) and log (b) scale as a function of the wall-normal distance expressed in wall units, $z^{+}$. Comparison between the reference case of unstratified turbulence ($Ri_{\tau }=0$), and the stratified turbulence at $Ri_{\tau }=50$, $Ri_{\tau }=100$, $Ri_{\tau }=200$ and $Ri_{\tau }=300$ (filled symbols). The correlation $\langle \varTheta ^{+} \rangle =Pr \cdot z^{+}$ and $\langle \varTheta ^{+} \rangle =(1/\kappa ) \log (z^{+})+B_{\theta }$, with coefficients $\kappa _{\theta }=0.436$ and $B_{\theta }=3.11$ taken from Alcántara-Ávila, Hoyas & Jezabel Pérez-Quiles (Reference Alcántara-Ávila, Hoyas and Jezabel Pérez-Quiles2021), is also shown for comparison (solid line, b). A close-up view of the mean temperature in the core region of the channel is offered in (b) for clarity.

To evaluate the influence of stratification on turbulence, we now look at the root mean square of the streamwise $\langle u'^{+}_{x,rms}\rangle$, spanwise $\langle u'^{+}_{y,rms}\rangle$ and wall-normal $\langle u'^{+}_{z,rms}\rangle$ velocity fluctuations, and at the root mean square of temperature fluctuations $\langle \varTheta '^{+}_{rms}\rangle$. Results, which are presented in figure 6, clearly show that all velocity and temperature fluctuations are essentially unaffected by stratification in the near-wall region ($z^{+}/h^{+} < 0.1$), where they recover the behaviour of canonical near-wall turbulence. Farther from the wall, in the region $0.3< z^{+}/h^{+}< 0.8$, $\langle u'^{+}_{x,rms}\rangle$ decreases first and increases later compared with the neutrally buoyant case, while the opposite behaviour – increasing first and decreasing later – is observed for $\langle u'^{+}_{z,rms}\rangle$. Note that the cross-over between the profiles of stratified and neutrally buoyant turbulence occurs around $z^{+}/h^{+}\simeq 0.5$, and is somehow influenced by $Ri_{\tau }$. In contrast, $\langle u'^{+}_{y,rms}\rangle$ is characterized by a clear increase that – provided $Ri_{\tau }>0$ – seems independent of the value of $Ri_{\tau }$. This behaviour, and in particular the increase of $\langle u'^{+}_{x,rms}\rangle$ and $\langle u'^{+}_{y,rms}\rangle$, is associated with the increase of velocity gradient (see also figure 7(a) and the discussion therein) in that region, which enhances turbulence production, while the corresponding decrease of $\langle u'^{+}_{z,rms}\rangle$ is due to the conversion of turbulent kinetic energy (in the vertical direction) into potential energy. In the core region of the channel, $z^{+}/h^{+}>0.8$, $\langle u'^{+}_{x,rms}\rangle$ strongly decreases while $\langle u'^{+}_{z,rms}\rangle$ develops a peak that is not visible in neutrally buoyant turbulence and is due to the presence of IGWs. Correspondingly, a marked peak at the channel centre is also observed for the temperature fluctuations, figure 6(d). Note that, for $\langle \varTheta '^{+}_{rms}\rangle$, there is a clear cross-over between the different cases of stratified turbulence: for increasing $Ri_{\tau }$, temperature fluctuations tend to increase in the transition layer, $0.1< z^{+}/h^{+}<0.8$, and decrease in the core layer $0.8< z^{+}/h^{+}<1$ (although they remain much larger than the neutrally buoyant case). To summarize, previous observations indicate that the structure of turbulence is not influenced near the wall, since velocity fluctuations, and also the ratio between the wall-normal and the streamwise velocity fluctuations, $\langle u'^{+}_{z,rms}\rangle /\langle u'^{+}_{x,rms}\rangle$, remains almost constant among the different cases. At the same time, in the buffer region there is a decrease of $\langle u'^{+}_{z,rms}\rangle /\langle u'^{+}_{x,rms}\rangle$, which indicates that stratification hinders the energy transfer from the streamwise to the wall-normal component. Then, at the channel centre, velocity and temperature fluctuations are essentially induced by IGWs.

Figure 6. Wall-normal behaviour of the root mean square of the velocity fluctuations in the streamwise direction ($\langle u'^{+}_{x,rms} \rangle$, panel a), in the spanwise direction ($\langle u'^{+}_{y,rms} \rangle$, panel b) and in the wall-normal direction ($\langle u'^{+}_{z,rms} \rangle$, panel c). The wall-normal behaviour of the temperature fluctuations is also shown ($\langle \varTheta '^{+}_{rms} \rangle$, panel d). Comparison between the reference case of unstratified turbulence ($Ri_{\tau }=0$), and the stratified turbulence at $Ri_{\tau }=50$, $Ri_{\tau }=100$, $Ri_{\tau }=200$ and $Ri_{\tau }=300$ (filled symbols).

Figure 7. (a) Wall-normal behaviour of the viscous shear stress, $\tau ^{v}_{xy}=\partial \langle u^{+}_x\rangle / \partial z^{+}$ and turbulent shear stress, $\tau ^{t}_{xy}=\langle u'^{+}_x u'^{+}_z\rangle$. The linear behaviour of the total shear stress, $\tau _{tot}/(\rho u_{\tau }^{2})$ is also shown by the black solid line. Comparison between the reference case of unstratified turbulence ($Ri_{\tau }=0$), and the cases of stratified turbulence at $Ri_{\tau }=50$, $Ri_{\tau }=100$, $Ri_{\tau }=200$ and $Ri_{\tau }=300$ (filled symbols). (b) Wall-normal behaviour of the diffusive heat flux, $q_{d}=Pr^{-1} \partial \langle \varTheta ^{+} \rangle / \partial z^{+}$ and turbulent heat flux (buoyancy flux), $q_t= \langle u'^{+}_z \varTheta '^{+} \rangle$. Comparison between the reference case of unstratified turbulence ($Ri_{\tau }=0$), and the cases of stratified turbulence at $Ri_{\tau }=50$, $Ri_{\tau }=100$, $Ri_{\tau }=200$ and $Ri_{\tau }=300$ (filled symbols).

3.3. Momentum and heat fluxes

We focus here on the wall-normal behaviour of momentum and heat fluxes, two key quantities in turbulent transport phenomena. The momentum flux can be obtained from the Reynolds-averaged streamwise momentum equation as

(3.1)\begin{equation} \tau_{tot}=1-\dfrac{2z^{+}}{Re_{\tau}}=\underbrace{\dfrac{\partial \langle u^{+}_x \rangle}{\partial z^{+}}}_{\tau^{v}_{xz}}+\underbrace{\langle u'^{+}_x u'^{+}_z\rangle}_{\tau^{t}_{xz}}, \end{equation}

where the turbulent $(\tau ^{t}_{xy})$ and viscous $(\tau ^{v}_{xy} )$ counterparts to the overall stress are explicitly indicated. The behaviour of $\tau ^{t}_{xy}$ and $\tau ^{v}_{xy}$ (symbols) is shown in figure 7(a), together with the behaviour of $\tau _{tot}$ (solid black line). Increasing $Ri_{\tau }$, we note a general reduction of $\tau ^{t}_{xy}$, in particular in the core region, where $\langle u'^{+}_x u'^{+}_z\rangle \simeq 0$ for $Ri_{\tau }\geqslant 200$. Accordingly, a corresponding increase of $\tau ^{v}_{xy}$ – and hence of the velocity gradient $\partial u_x/\partial z$ – is observed, so that the overall linear behaviour of the total shear stress is recovered (3.1). This is clearly visualized in the inset of figure 7(a), where a close-up view of the behaviour of $\tau ^{t}_{xy}$ and $\tau ^{v}_{xy}$ at the channel centre, $0.8< z^{+}/h^{+}<1$, is given. In view of the present results, it is clear that the velocity increase observed at the channel centre (see figure 4) is due to the reduction of turbulent momentum flux in the wall-normal direction (i.e. reduction of $u'_x u'_z$), and to the corresponding increase of the relative importance of $\tau ^{v}_{xy}$ therein (Armenio & Sarkar Reference Armenio and Sarkar2002; Yeo et al. Reference Yeo, Kim and Lee2009).

Linked to the previous analysis of the wall-normal momentum flux, we now consider the wall-normal heat flux, whose behaviour can be obtained from the Reynolds-averaged energy balance equation as

(3.2)\begin{equation} \langle u'_z \theta' \rangle - \alpha \frac{\partial \langle \theta \rangle}{\partial z}={-}\alpha \left[ \frac{\partial \langle \theta \rangle }{\partial z}\right]_w=q_w, \end{equation}

where $\alpha =\nu /Pr$ is the thermal diffusivity. Normalizing (3.2) by $q_w=-\alpha [{\partial \theta }/{\partial z}]_w$, and recalling that $q_w=\theta _{\tau } u_{\tau }$, we finally obtain (Armenio & Sarkar Reference Armenio and Sarkar2002; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011)

(3.3)\begin{equation} \frac{ \langle u'_z \theta' \rangle }{q_w} - \frac{\alpha \dfrac{\partial \theta}{\partial z}}{q_w}=\underbrace{ \langle u'^{+}_z \varTheta'^{+} \rangle }_{{q_t}}- \underbrace{ \frac{1}{Pr} \frac{\partial \langle \varTheta^{+}\rangle}{\partial z^{+}}}_{{q_d}}=1. \end{equation}

The two terms $q_t$ and $q_d$, explicitly indicated in (3.3), represent the turbulent (usually referred to as buoyancy flux) and the diffusive counterparts to the total heat flux, and their behaviour is shown in figure 7(b). By looking at the profile of $q_t$, it is apparent that, while moving away from the wall – and no matter the value of $Ri_{\tau }$ – $q_t$ increases sharply until it reaches a maximum value of approximately $q_t\simeq 0.95$ around $z^{+}/h^{+}\simeq 0.1$. In neutrally buoyant conditions, this peak value is kept almost unaltered throughout the entire channel, clearly corresponding to the constant flux hypothesis customarily assumed in neutral boundary layers (Tennekes & Lumley Reference Tennekes and Lumley1972; Ortiz-Suslow et al. Reference Ortiz-Suslow, Kalogiros, Yamaguchi and Wang2021). At larger $Ri_{\tau }$, we observe a significant decrease of $q_t$ in the core region of the channel. This decrease is so important that, for $Ri_{\tau }>200$, $q_t\simeq 0$. Interestingly, and in agreement with previous observations (Ohya et al. Reference Ohya, Neff and Meroney1997; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011), we cannot find evidence of the speculated mean counter-gradient flux (or, in other words, negative $q_t$) at the centre of the channel (Komori et al. Reference Komori, Ueda, Ogino and Mizushina1983; Armenio & Sarkar Reference Armenio and Sarkar2002).

We further observe that previous studies discussing the possibility of a mean counter-gradient flux considered low-$Re_{\tau }$ flows. Here, we focus on flows at much higher $Re_{\tau }$, for which the energy containing boundary layer structures cannot reach very far from the wall. Since we may assume that the mean counter-gradient flux is the outcome of a number of local counter-gradient events (uncorrelated thermal and flow structures) driven by the interaction between near-wall structures and IGW, it is unlikely that a mean counter-gradient flux can be observed at high $Re_{\tau }$. However, we cannot exclude the presence of a mean counter-gradient flux in strongly stratified conditions, when the flow becomes largely inhomogeneous and intermittent. Between the near-wall and the core regions of the channel there is a region in which $q_t$ remains almost constant and close to unity for all cases considered here, with only a slight decrease for increasing $Ri_{\tau }$. The diffusive heat flux $q_d$ has a mirror-like behaviour compared with $q_t$, since the total heat flux is constant across the channel (see (3.3)): $q_d$ decreases sharply while moving away from the wall and it subsequently increases – with the only exception of $Ri_{\tau }=0$, for which it remains uniform and very low – in the core region of the channel. This trend of $q_d$ is important to explain the non-monotonic behaviour of the temperature profile observed in figure 5 ($q_d$ is, by definition, proportional to the mean temperature gradient). Note indeed that the mean temperature gradient along the wall-normal direction can be conveniently expressed as (Armenio & Sarkar Reference Armenio and Sarkar2002)

(3.4)\begin{equation} \frac{\partial \langle \theta \rangle}{\partial z}=\left[\frac{\partial \langle \theta \rangle}{\partial z}\right ]_w \left[ 1-\frac{\langle u'_z \theta' \rangle}{u_{\tau}\theta_{\tau}} \right], \end{equation}

where $[\partial \langle \theta \rangle /\partial z ]_w$ is the dimensionless mean temperature gradient at the wall, i.e. the Nusselt number. Equation (3.4), together with the observation that the buoyancy flux $q_t=\langle u'_z \theta '\rangle /u_{\tau }\theta _{\tau }= \langle u_z'^{+} \varTheta '^{+} \rangle$ becomes almost zero for $Ri_\tau >200$ (see figure 7b), indicates that for large stratification the temperature gradient at the channel centre perfectly matches the temperature gradient at the wall. Since the temperature gradient at the wall (i.e. the Nusselt number) decreases for increasing $Ri_{\tau }$, the same does the temperature gradient at the channel centre (but only once $Ri_{\tau }$ is large enough for $q_t$ to be approximately zero). To summarize, the temperature gradient at the channel centre initially increases for increasing stratification, until the stratification becomes so strong to completely damp the turbulent heat transfer (buoyancy flux), i.e. $q_{t}=\langle u'^{+}_z \varTheta '^{+} \rangle \simeq 0$, at the channel centre. When it happens, the temperature gradient at the channel centre – which becomes equal to the temperature gradient at the wall (Nusselt number) – decreases for increasing $Ri_{\tau }$.

Despite the large fluctuations of temperature and wall-normal velocity observed around the channel centre for increasing $Ri_{\tau }$ (figure 6c,d), the buoyancy flux $\langle u'^{+}_z \varTheta '^{+} \rangle$ reduces remarkably down to the point at which, for $Ri_{\tau }\geqslant 200$, it completely vanishes. To understand the reason behind this behaviour, we focus on the distribution of temperature, $\varTheta '^{+}$, and wall-normal velocity fluctuations, $u'^{+}_z$, on a wall-parallel plane $(x\unicode{x2013}y)$ located at the channel centre. This is shown in figure 8 for $Ri_{\tau }=0$ and $Ri_{\tau }=300$. While for $Ri_{\tau }=0$ (figure 8a,b), the picture displays the typical features of fully developed turbulence, with spots of positive and negative temperature and wall-normal velocity randomly distributed over the considered plane, for $Ri_{\tau }=300$ (figure 8c,d) the situation is different, and there seems to be a connection between $\varTheta '^{+}$ and $u'^{+}_z$, which are both organized into stripes extending along the spanwise direction $y$. However, stripes of $u'^{+}_z$ appear to lag behind those of $\varTheta '^{+}$. The phase lag between $u_z'^{+}$ and $\varTheta '^{+}$ distributions can be estimated by looking at the coherency spectrum

(3.5)\begin{equation} C_{u'^{+}_z,\varTheta'^{+}}(\kappa_w)=\int_{-\infty}^{\infty} {R_{u'^{+}_z,\varTheta'^{+}}(s) \exp \left({-}2{\rm \pi} j s \kappa_w \right) \mbox{d} s }, \end{equation}

with $\kappa _w=kh$ the dimensionless wavenumber and $R_{u'^{+}_z,\varTheta '^{+}}$ the correlation coefficient between $u'^{+}_z$ and $\varTheta '^{+}$. In particular, it is intriguing to focus on the phase of $C_{u'^{+}_z,\varTheta '^{+}}(\kappa _w)$. As shown in figure 8( f), we clearly notice that, for $Ri_{\tau }=300$, $u'^{+}_z$ and $\varTheta '^{+}$ are shifted by ${\rm \pi} /2$ (see Iida et al. Reference Iida, Kasagi and Nagano2002). Such a phase delay is instead not observed at $Ri_{\tau }=0$ (see figure 8e). It is important to note that the presence of this phase shift explains why, although temperature and wall-normal fluctuations are both very large at the channel centre, their correlation – that is, the buoyancy flux $\langle u'^{+}_z \varTheta '^{+} \rangle$ – is almost zero: where $\varTheta '^{+}$ is maximum in magnitude, $u'^{+}_z \simeq 0$ and vice versa, consistently with the presence of a wavy motion (IGW).

Figure 8. Temperature and velocity fluctuations on a wall-parallel $(x\unicode{x2013}y)$ plane located at the channel centre for $Ri_{\tau }=0$ (a,b) and for $Ri_{\tau }=300$ (c,d). Phase of the cross-spectrum, $\phi (C_{u'^{+}_z, \varTheta '^{+}} )$, as a function of the wavenumber $\kappa _w=kh$ for $Ri_{\tau }=0$ (e) and for $Ri_{\tau }=300$ ( f).

Not only the average value of the buoyancy flux $q_t$ is important, but also its distribution in space and time. To understand it, we focus on the two limiting cases $Ri_{\tau }=0$ and $Ri_{\tau }=300$ and we look at the behaviour of $q_t$ on a wall-parallel plane $(x\unicode{x2013}y)$ located at the channel centre (figure 9a,b), and normalized by the corresponding maximum value $q_{t,max}$ observed on that plane. Then, we compute the probability density function $\varPi (q_t/q_{t,max})$ (figure 9c). For $Ri_{\tau }=0$ (blue line in figure 9c), $\varPi (q_t/q_{t,max})$ is highly asymmetric, with the most probable value occurring for $q_t/q_{t,max}=0$ and with larger positive fluctuations compared with negative ones. This suggests that, although the mean temperature difference between the walls induces a net positive wall-normal energy flux (i.e. the mean value $\langle u'^{+}_z \varTheta '^{+} \rangle > 0$), $u'^{+}_z \varTheta '^{+}$ can often be negative, indicating the presence of regions characterized by local counter-gradient heat fluxes. The occurrence of localized counter-gradient heat fluxes is an extremely important phenomenon that has been observed also in other situations (Huisman et al. Reference Huisman, van Gils, Grossmann, Sun and Lohse2012; Zonta & Chibbaro Reference Zonta and Chibbaro2016; Hadi Sichani et al. Reference Hadi Sichani, Marchioli, Zonta and Soldati2020). From a physical point of view, small positive and negative values of $\varTheta '^{+} u'^{+}_z$ are due to turbulence, which is characterized by uncorrelated velocity and temperature fluctuations. These small positive and negative values of $\varTheta '^{+} u'^{+}_z$, which are equally likely events, balance each other and do not contribute to the net heat transport (Shang, Tong & Xia Reference Shang, Tong and Xia2005). Only large velocity and temperature fluctuations produced by larger coherent structures are correlated and contribute to the net heat flux. And clearly, large positive fluctuations are larger than negative ones, so to produce a positive net heat flux.

Figure 9. Contour map of the turbulent heat flux $q_t=\langle \varTheta ^{+} u'^{+}_z\rangle$ on a wall-parallel $(x\unicode{x2013}y)$ plane located at the channel centre for $Ri_{\tau }=0$ (a), for $Ri_{\tau }=300$ (b) and corresponding probability density function $\varPi (q_t/q_{t, max})$ (c).

The asymmetry of $\varPi (q_t/q_{t,max})$ is almost completely absent at $Ri_{\tau }=300$. In particular, while the probability of observing positive fluctuations is globally reduced (i.e. there is a reduction of large correlated velocity and temperature fluctuations), negative fluctuations are rather persistent, so that the average value of $q_t/q_{t,max}$ is close to zero. This is nicely rendered by the contour maps of $q_t/q_{t,max}$ shown in figure 9(b): while most of the plane is characterized by a correlation $q_t/q_{t,max}\approx 0$ (and corresponding to the most probable value of $\varPi (q_t/q_{t, max})$ in figure 9), small and rare patches of large positive $q_t/q_{t,max}$ coexist with small and rare patches of large negative $q_t/q_{t,max}$.

3.4. Internal gravity waves

In the previous chapter, we described IGWs on a qualitative basis (see in particular figures 2–3 and comments therein). To characterize such structures on a more quantitative basis, we introduce here the buoyancy frequency $N$ (or Brunt–Väisälä frequency), which represents the frequency of oscillation – in a stratified medium – of a lump of fluid that is displaced by velocity fluctuations in the wall-normal direction. In dimensionless form, it is defined as

(3.6)\begin{equation} {\mathcal{N}} =N h / u_{\tau}= \sqrt{- g \beta \frac{\partial \langle \theta \rangle}{\partial z}}\frac{h}{u_{\tau}}, \end{equation}

and can naturally be used to evaluate the characteristic frequency of oscillation of IGWs. The behaviour of ${\mathcal {N}}$ as a function of the wall-normal distance is shown in figure 10(a) for the different values of $Ri_{\tau }$. We observe that ${\mathcal {N}}$ is large at the centre of the channel – and peaks at the centreline – due to the sharpening of the mean temperature profile therein (figure 5), while it is rather small in the remaining part of the channel. Note that the behaviour of ${\mathcal {N}}$ in the proximity of the wall – shaded area in figure 10(b) – is not particularly important for the present discussion, since the temperature gradient in that region is mostly influenced by the presence of the boundary rather than by buoyancy. In the framework of the present discussion, it is interesting to look at the time-averaged premultiplied streamwise energy spectra of wall-normal velocity fluctuations, $k_x E_{u'^{2}_z} (k_x)$, at the channel centreline (Iida et al. Reference Iida, Kasagi and Nagano2002; Yeo et al. Reference Yeo, Kim and Lee2009; Zonta et al. Reference Zonta, Onorato and Soldati2012b). Results are presented in figure 10(b) for the different values of $Ri_{\tau }$. We consider in particular the case $Ri_{\tau }=300$, for which stratification is more important. We note that $k_x E_{u'^{2}_z} (k_x)$ develops a peak – which is not observed for $Ri_{\tau }=0$ – at the dimensionless wavenumber $\hat {k}_x=k_{x} h \approx 5$ (where $k_{x}$ is the wavenumber in dimensional form). We can hypothesize that this peak is due to the presence of IGWs, which induce vertical velocity fluctuations at the frequency ${\mathcal {N}}$ (Zonta et al. Reference Zonta, Onorato and Soldati2012b). To demonstrate this relationship, we assume that IGWs move in the streamwise direction with the mean centreline velocity, $u_{IGW}$. For $Ri_{\tau }=300$, $u^{+}_{IGW} \simeq 47$ (or, in dimensional form, $u_{IGW} = u_{\tau } u^{+}_{IGW} \simeq 1 \ {\rm m}\ {\rm s}^{-1}$). The corresponding frequency of such waves is $\omega _{IGW} = u_{IGW} \times k_{x} / 2 {\rm \pi}\simeq 1 \ {\rm s}^{-1}$. Considering now the value of the buoyancy frequency at the channel centreline, we get ${\mathcal {N}}\simeq 38$ or, in dimensional form, $N = {\mathcal {N}} u_{\tau } / h \simeq 1 \ {\rm s}^{-1}$, in close agreement with the previous estimate. Also shown in figure 10(a) is the range of frequency peaks, which correlates well with the range of wavenumber peaks shown in figure 10(b). The previous estimate is based on an oversimplified scenario, in which IGWs are monochromatic waves that travel at the advection velocity and are characterized by the frequency $N$. However, it is known (Staquet & Sommeria Reference Staquet and Sommeria2002) that the frequency of IGWs obeys the dispersion relation

(3.7)\begin{equation} \omega={\pm} N \cos \phi={\pm} N \frac{\sqrt{k_x^{2}+k_y^{2}}}{\sqrt{k_x^{2}+k_y^{2}+k_z^{2}}}, \end{equation}

in which the wave vector $\boldsymbol {k}=(k_x,k_y,k_z)$ is not necessarily directed along $x$, and can have an angle $\phi$ to the horizontal plane. From (3.7), it clearly emerges that IGWs are characterized by frequencies in the range $0< \omega \leqslant N$. In addition, due to the advection velocity $u_{IGW}$ at the channel centre, the observed wave frequency is modified by the Doppler shift, which takes the form $\omega _{DS}=u_{IGW} \hat {k}_x$. Therefore, the frequency of IGWs can be expressed by the relation

(3.8)\begin{equation} \omega=\omega_{DS}\pm N \cos \phi = u_{IGW} \hat{k}_x \pm N \cos \phi. \end{equation}

From (3.8), and under the assumption of $\cos \phi =1$, it emerges that the observed peak wavenumber $\hat {k}_x$ in figure 10(b) corresponds to the Doppler shift $u_{IGW} \hat {k}_x=N$. This result seems consistent with the recent analysis of Maffioli, Delache & Godeferd (Reference Maffioli, Delache and Godeferd2020), in which – employing spatio-temporal analysis of an unbounded stratified flow – it was shown that the dominant contributions to the wave energy in stratified turbulence come from frequencies close to $N$.

Figure 10. (a) Wall-normal behaviour of the mean Brunt–Väisälä frequency $N$ for the different cases of stratified turbulence ($Ri_{\tau }=50$, $Ri_{\tau }=100$, $Ri_{\tau }=200$ and $Ri_{\tau }=300$). (b) Time-averaged premultiplied streamwise energy spectra of wall-normal velocity fluctuations $k_xE(u'_z u'_z)$ evaluated on a wall-parallel $(x\unicode{x2013}y)$ plane located at the channel centre.

Using the buoyancy frequency $N$, we can calculate the gradient Richardson number as

(3.9)\begin{equation} Ri_g = \frac{N^{2}}{S^{2}} = \frac{- g \beta \dfrac{\partial \langle \theta \rangle}{\partial z}}{\left( \dfrac{\partial \langle u_{x} \rangle}{\partial z} \right)^{2}}, \end{equation}

where $S$ is the mean shear rate. The gradient Richardson number is a measure of the relative importance of the buoyant consumption of turbulence, $N^{2}$, compared with the shear production of turbulence, $S^{2}$. In the literature (Miles Reference Miles1961; Rohr et al. Reference Rohr, Itsweire, Helland and Van Atta1988; Holt, Koseff & Ferziger Reference Holt, Koseff and Ferziger1992; Piccirillo & Van Atta Reference Piccirillo and Van Atta1997), $Ri_g \simeq 0.25$ has been identified as a critical threshold to evaluate the flow structure in unbounded stratified flows, where $Ri_g$ is uniform in space. For $Ri_g \simeq 0.25$, turbulence does not grow or decay. At larger $Ri_g$, turbulence decays, while at smaller $Ri_g$ turbulence grows. Naturally, in wall-bounded flows, $Ri_{g}$ varies in the wall-normal direction from small values close to the wall, to large values close to the centre of the channel. This is represented in figure 11, where we plot $Ri_{g}$ in logarithmic scale as a function of the dimensionless wall-normal distance, for different $Ri_{\tau }$. Current results are in good agreement with previous observations (Armenio & Sarkar Reference Armenio and Sarkar2002; Taylor, Sarkar & Armenio Reference Taylor, Sarkar and Armenio2005; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Zonta et al. Reference Zonta, Onorato and Soldati2012b), showing in particular the sharp increase of $Ri_g$ around $Ri_g \simeq 0.2$. This sharp increase, which occurs at the channel centre (where $\partial \langle \theta \rangle / \partial z$ is large while $\partial \langle u_x \rangle / \partial z \rightarrow 0$) identifies the region of the channel in which buoyant consumption of turbulence dominates over shear production of turbulence, and sets the boundaries where IGWs are observed ($z^{+}/h^{+}>0.8$). Based on the previous observations, it seems reasonable to set the value $Ri_g=0.2$ as the critical threshold above which turbulence is significantly damped by buoyancy in wall-bounded stratified turbulence.

Figure 11. Wall-normal behaviour of the gradient Richardson number $Ri_g$ (a) and of the turbulent Prandtl number $Pr_t$ (b) for the different values of $Ri_{\tau }$. The threshold value $Ri_g=0.2$ is also explicitly shown in panel (a) (dashed line).

A further crucial aspect of stratified turbulence is the evaluation of irreversible mixing of momentum and energy (Fernando Reference Fernando1991; Ivey, Winters & Koseff Reference Ivey, Winters and Koseff2008). A common approach to evaluating mixing in stratified turbulence is to rely on the concept of turbulent eddy viscosity, ${\mathcal {K}}_m$, and diffusivity, ${\mathcal {K}}_\theta$, which are defined as (Pope Reference Pope2000)

(3.10a,b)\begin{equation} {\mathcal{K}}_m = \frac{\langle u^{\prime} w^{\prime} \rangle} {\dfrac{\partial \langle u \rangle}{\partial z}},\quad {\mathcal{K}}_{\theta} = \frac{\langle \theta^{\prime} w^{\prime} \rangle} {\dfrac{ \partial \langle \theta \rangle}{\partial z}}. \end{equation}

The ratio between ${\mathcal {K}}_m$ and ${\mathcal {K}}_\theta$ is usually referred to as turbulent Prandtl number (Launder & Spalding Reference Launder and Spalding1972), and can be written in dimensionless form as

(3.11)\begin{equation} Pr_T = \frac{\langle u^{\prime} w^{\prime} \rangle \dfrac{\partial \langle \theta \rangle}{\partial z}}{\langle w^{\prime} \theta^{\prime} \rangle \dfrac{\partial \langle u \rangle}{\partial z}}. \end{equation}

The evaluation of $Pr_T$ as a function of the distance from the wall and for different stratification levels (i.e. different $Ri_\tau$), is a major aspect in the field of turbulence modelling (Grachev et al. Reference Grachev, Andreas, Fairall, Guest and Persson2007; Venayagamoorthy & Stretch Reference Venayagamoorthy and Stretch2010). In figure 11(b) we plot $Pr_T$ as a function of $z^{+}/h^{+}$. Near the boundary, and regardless of $Ri_\tau$, we observe that $Pr_T \simeq 1$. Above $z^{+}/h^{+} \simeq 0.2$, there is slight dependence of $Pr_T$ on $Ri_\tau$, but $Pr_T$ remains bounded in the range $0.8< Pr_T<1.2$. At the channel centre, where stratification effects are stronger, the behaviour of $Pr_T$ does not have a clear trend for increasing $Ri_\tau$, although $Pr_T$ does not depart much from unity. We therefore conclude that the assumption $Pr_{T} \approx 1$, customarily taken in turbulence modelling (Laskowski et al. Reference Laskowski, Kearney, Evans and Greif2007; Zonta Reference Zonta2013), seems a reasonable choice for mild to moderate stratification levels.

3.5. Macroscopic characterization of the flow: $Ri_b$, $C_f$, $Nu$

Although the shear Richardson number $Ri_{\tau }$ is customarily used to characterize the flow in numerical simulations of wall-bounded stratified turbulence (Zonta & Soldati Reference Zonta and Soldati2018), its use in experiments/field measurements is much more limited. The reason is the difficult experimental evaluation of the shear velocity, which in turn requires the measurement of the shear stress at the wall. As a consequence, in many experiments the bulk Richardson number $Ri_b=\beta (\Delta T/2) g h/(2u^{2}_{b})$ is usually preferred, since the bulk velocity $u_b$ is an easily accessible quantity. To draw a link between simulations and experiments/field measurements, it is interesting to evaluate the behaviour of the bulk Richardson number $Ri_{b}$, as a function of the shear Richardson number $Ri_{\tau }$. This behaviour is shown in figure 12. Present results (filled symbols, $\bullet$, red), which are plotted together with literature results obtained at lower $Re_{\tau }$ (open symbols), confirm previous indications that $Ri_b\propto Ri_{\tau }^{2/3}$. In view of its independence on $Re_{\tau }$ – at least over approximately one decade, from $Re_{\tau }=180$ up to $Re_{\tau }=1000$ – the proposed $Ri_b$ scaling seems a general property of the flow that can be used as a robust parametrization tool, as will be shown below.

Figure 12. Bulk Richardson number, $Ri_{b}$, as a function of the shear Richardson number, $Ri_{\tau }$, for the present simulations at $Re_{\tau }=1000$, and for previous simulations at $Re_{\tau }=180$ and $Re_{\tau }=550$ (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011). The proposed scaling $Re_{b}\sim Ri_{\tau }^{2/3}$ is also explicitly indicated.

Perhaps the most important quantities to be monitored in wall-bounded stratified turbulence are the overall momentum and heat transfer rates, which are commonly quantified in terms of the friction factor $C_f$ – shear stress to kinetic energy ratio – and by the Nusselt number $Nu$ – convective to conductive heat transfer ratio – as

(3.12)\begin{equation} \left.\begin{array}{l@{}} C_f=\dfrac{2 \tau_w}{\rho u_b^{2}},\\ Nu=\dfrac{2 q_w h}{\lambda \Delta T }. \end{array} \right\} \end{equation}

Results obtained from the present simulations are shown by filled symbols ($\bullet$, red) in figures 13–14, together with the results obtained in previous studies (Garg et al. Reference Garg, Ferziger, Monismith and Koseff2000; Armenio & Sarkar Reference Armenio and Sarkar2002; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Zonta et al. Reference Zonta, Onorato and Soldati2012b; Zonta Reference Zonta2013). Focusing on $C_f$ (figure 13), it is clear that an increase of stratification (increase of $Ri_{\tau }$) reduces the wall-normal momentum transfer, via the reduction of wall-normal turbulent transport and via the corresponding increase of the volume flowrate (keeping the imposed pressure gradient constant). Interestingly, recalling that $Ri_b\propto Ri_{\tau }^{2}/3$ (see figure 12), and considering that $C_f/4=Ri_b/Ri_{\tau }$, we obtain $C_f \sim Ri_{\tau }^{-1/3}$ (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011). This scaling law, which is shown in figure 13 together with the reported numerical results, predicts fairly well the behaviour of $C_f$ for a broad range of $Ri_{\tau }$, with some departure observed only at very large stratification levels (when local flow laminarization is likely to appear).

Figure 13. Friction factor $C_f$ as a function of the shear Richardson number $Ri_{\tau }$. Results of present study (filled symbols) are shown together with results obtained in previous studies (Garg et al. Reference Garg, Ferziger, Monismith and Koseff2000; Armenio & Sarkar Reference Armenio and Sarkar2002; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Zonta et al. Reference Zonta, Onorato and Soldati2012b; Zonta Reference Zonta2013). The proposed scaling $C_{f}\sim Ri_{\tau }^{-1/3}$ is also explicitly indicated.

Figure 14. Rescaled Nusselt number $Nu \times Re_{\tau }^{-2/3}$ as a function of the shear Richardson number $Ri_{\tau }$. Results of present study (filled symbols) are shown together with results obtained in previous studies (Garg et al. Reference Garg, Ferziger, Monismith and Koseff2000; Armenio & Sarkar Reference Armenio and Sarkar2002; García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Zonta et al. Reference Zonta, Onorato and Soldati2012b; Zonta Reference Zonta2013). Note that the values of $Nu$ in Zonta et al. (Reference Zonta, Onorato and Soldati2012b) and Zonta (Reference Zonta2013), obtained at $Pr=3$, have been rescaled by $Pr^{1/3}$. The proposed scaling $Nu \times Re_{\tau }^{-2/3} \sim Ri_{\tau }^{-1/3}$ is also explicitly indicated.

The fair collapse observed for $C_f$ is not recovered for the Nusselt Number $Nu$ (not shown in figure 14). Interestingly, when rescaled by $Re_{\tau }^{-2/3}$, the collected results collapse nicely (figure 14) and scale as $Nu \times Re_{\tau }^{-2/3} \sim Ri_{\tau }^{-1/3}$ (therefore giving $Nu \times Re_{\tau }^{-2/3} \sim C_f$). The proposed scaling differs from classical analogies widely used in the literature – such as the Chilton and Colburn one $C_f/2=Nu \times Re^{-1} \times Pr^{-1/3}$ (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot2007) – to relate heat, mass and momentum transfer coefficients in forced convection phenomena. It is reasonable to expect that the interaction between buoyancy and momentum modifies the main transport processes in a non-trivial way that is hard to predict by simplified assumptions. An in-depth analysis and a corresponding accurate parametrization of the actual flow field are therefore required to shed some light on the proposed scaling, which, however, appears rather robust.