2. Computational campaign

We perform DNSs of stably stratified turbulent channel flows to develop a framework for describing their mean velocity profile, which is subsequently validated using the simulation data. The Navier–Stokes equations are solved in the low-Mach-number regime (Majda & Sethian Reference Majda and Sethian1985; Cook & Riley Reference Cook and Riley1996), with a prescribed temperature difference of $1\,\%$ . This set-up ensures minimal variation in thermodynamic properties, enabling a reasonable comparison with cases in the literature that use the Oberbeck–Boussinesq approximation (see, e.g. Garcia-Villalba & Del Alamo Reference Garcia-Villalba and Del Alamo2011). The governing equations are

(2.1) \begin{align} \frac {\partial \rho }{\partial t} + \frac {\partial \rho u_j}{\partial x_j} = 0, \\[-28pt] \nonumber \end{align}

(2.2) \begin{align} \rho \frac {\partial u_i}{\partial t} + \rho u_j\frac {\partial u_i}{\partial x_j} = -\frac {\partial p}{\partial x_i} + \frac {\partial \tau _{\textit{ij}}}{\partial x_j}+\rho g_i+f_i\delta _{i,1}, \\[-28pt] \nonumber \end{align}

(2.3) \begin{align} \rho C_p\frac {\partial T}{\partial t} + \rho C_p u_j\frac {\partial T}{\partial x_j} = \frac {\partial }{\partial x_j}\left (\lambda \frac {\partial T} {\partial x_j}\right) + \frac {\mathrm{d} p_0}{\mathrm{d} t}, \\[-2pt] \nonumber \end{align}

where $\tau _{\textit{ij}} = \mu (\partial u_i/ \partial x_j + \partial u_j/ \partial x_i - 2/3\ \partial u_k/\partial x_k\ \delta _{i,j})$ . Here, $p_0$ denotes the spatially invariant thermodynamic pressure, while $p$ represents the hydrodynamic pressure. The walls are maintained at fixed temperatures to enforce stable stratification. The arithmetic mean of the wall temperatures is chosen as the reference temperature at which the reference thermodynamic properties such as thermal conductivity ( $\lambda$ ), heat capacity ( $C_p$ ), dynamic viscosity ( $\mu$ ) and density ( $\rho$ ) are evaluated. The density is related to temperature ( $T$ ) via the ideal gas law ( $\rho = p_0/ R T$ ), where $R$ is the universal gas constant. The gravitational acceleration vector $g_i$ acts vertically downward, and $f_1$ represents the constant pressure difference driving the flow along the streamwise direction. Periodic boundary conditions are imposed in the streamwise and spanwise directions, while the no-slip condition is enforced at the walls. The equations were solved on staggered Cartesian grids, with velocities at cell faces and scalars at the cell centres. The convective and diffusive terms were discretised in space using a second-order, finite volume scheme (Costa Reference Costa2018), and advanced in time using Wray’s low-storage third-order Runge–Kutta method, with implicit treatment for wall-normal diffusion.

This set-up is characterised by a friction Reynolds number $\textit{Re}_{\tau } = \rho _0 u_{\tau } h/\mu _0$ , and a friction Richardson number, $Ri_{\tau } = \Delta \rho g h/\rho _0 u_{\tau }^2$ (Garg et al. Reference Garg, Ferziger, Monismith and Koseff2000); $\rho _0$ and $\mu _0$ denote the fluid density and dynamic viscosity at a reference temperature $\theta _0$ , $h$ is the half-channel height and $\Delta \rho$ the bottom-to-top wall density difference. The DNSs were performed at $\textit{Re}_{\tau } = 395$ and $550$ . Additionally, the data at $\textit{Re}_{\tau } = 1000$ were obtained from Zonta et al. (Reference Zonta, Sichani and Soldati2022). For each $\textit{Re}_{\tau }$ , a range of $Ri_{\tau }$ was explored ( $0 \leqslant Ri_{\tau } \leqslant 900$ ), at fixed Prandtl number $\textit{Pr} = 0.71$ . Figure 2 summarises the cases.

All simulations were conducted in a domain of dimensions, $L_x \times L_y \times L_z = 6\pi h \times 2\pi h \times 2h$ , in the streamwise, spanwise and wall-normal directions, respectively. At $\textit{Re}_{\tau } = 395$ , nine different $Ri_{\tau }$ were simulated using a grid of $ N_x \times N_y \times N_z = 1024 \times 512 \times 300$ points. On the other hand, eight different $Ri_{\tau }$ were simulated at $\textit{Re}_{\tau } = 550$ using a grid of $ N_x \times N_y \times N_z = 1536 \times 768 \times 480$ points. Uniform grid spacing was used in the periodic directions, while wall-normal grid stretching was applied to resolve the Kolmogorov length scale near the wall. Previous studies (e.g. Garcia-Villalba & Del Alamo (Reference Garcia-Villalba and Del Alamo2011)) have indicated that increasing stratification does not affect the smallest turbulent scales, so grid resolution requirements remain unchanged with increasing $Ri_{\tau }$ . Simulations at higher $Ri_{\tau }$ were initialised using steady-state fields from the preceding lower $Ri_{\tau }$ case, following the strategy employed in earlier studies. Grid details of cases at $\textit{Re}_{\tau } = 1000$ are available in Zonta et al. (Reference Zonta, Sichani and Soldati2022).

Figure 2. Summary of the parameter space sampled in the DNS database, with data at $\textit{Re}_{\tau } = 1000$ from Zonta et al. (Reference Zonta, Sichani and Soldati2022). Data for each case will be presented consistently with these marker colours.

3. Suitability of MOST for pressure-driven channels

While MOST was originally formulated for the ABL, we investigate its suitability in capturing the effects of stratification in turbulent channel flow. The applicability is assessed by examining the dependence of normalised velocity gradients on the wall-based $(z/L)$ and local $(z/\varLambda)$ stability parameters, as introduced in (1.1) and (1.2).

The function $\phi$ , commonly referred to as the stability correction function, is typically obtained empirically by fitting experimental data. The linear relations proposed by Webb (Reference Webb1970), Businger & Yaglom (Reference Businger and Yaglom1971)and Dyer (Reference Dyer1974) provide a good fit for observations within the dry ABL, particularly under weak stratification. However, the functional form is not universal, with these studies reporting slightly different slopes for the linear relation, typically ranging between $4.5$ and $5.5$ . Moreover, several alternative formulations for  $\phi$ have also been proposed, including higher-order functions that reduce to a linear form under weak stratification but approach a constant value under strong stratification (Beljaars & Holtslag Reference Beljaars and Holtslag1991; Chenge & Brutsaert Reference Chenge and Brutsaert2005). Through their linear behaviour under weak and nonlinear corrections under strong stratification, these formulations recover the well-known ‘ $z$ -less’ region observed in the ABL (Nieuwstadt Reference Nieuwstadt1984). However, this region has not been observed in channel flows (Armenio & Sarkar Reference Armenio and Sarkar2002; Garcia-Villalba & Del Alamo Reference Garcia-Villalba and Del Alamo2011; Zonta & Soldati Reference Zonta and Soldati2018).

In turbulent channel flows, the mean velocity (and, thus, its gradient) has been reported to increase with stratification. This is attributed to the reduction in turbulent shear stress under stable stratification. To represent this behaviour, we adopt the linear formulation proposed by Businger & Yaglom (Reference Businger and Yaglom1971), which captures the monotonic increase in velocity gradients with increasing stability. Figure 3 presents the dimensionless velocity gradients, $(\kappa z/u_{\tau } \ \mathrm{d} u/\mathrm{d}z)$ , obtained from DNSs across a range of stratification levels as a function of the stability parameter: $z/L$ in panel (a), and $z/\varLambda$ in panel (b). The figure also shows the Businger relation $1 + 4.7\zeta$ , where $\zeta$ is the stability parameter. Clearly, the local scaling captures the variation in velocity gradients more consistently than the wall-based scaling, particularly in representing the increase in gradient magnitude. Interestingly, although the Businger relation has been primarily validated within the ‘constant-flux’ layer, we find that it captures the velocity gradients reasonably well as long as turbulent shear remains the dominant mechanism and stratification effects are relatively weak ( $z/\varLambda \lt 1$ ). For $z/\varLambda \gt 1$ , deviations from the linear relation become evident, with the Businger relation tending to overestimate the velocity gradients (see short dash-dotted lines in figure 3).

To improve the similarity hypothesis, in the context of ABL, Yokoyama, Gamo & Yamamoto (Reference Yokoyama, Gamo and Yamamoto1979) proposed using ‘ $z$ ’-dependent surface parameters (e.g. $u_{\tau }(z) = u_{\tau }(1 - z/h)^{\gamma }$ ) that reliably estimate local fluxes. This method was subsequently adopted by Sorbjan (Reference Sorbjan1986) and more recently by Gryning et al. (Reference Gryning, Batchvarova, Brümmer, Jørgensen and Larsen2007) and Shen et al. (Reference Shen, Liu, Lu and Stevens2025). The exponent $\gamma$ is determined empirically and varies across different studies. Following a similar approach, we apply a correction of the form $(1 - z/h)^{\gamma }$ to the Businger relation, with $\gamma = 1/2$ . As shown in figure 3(b), the corrected relation captures the trend of the velocity gradients reasonably well for $z/\varLambda \gt 1$ , particularly in high $Ri_{\tau }$ cases, where this region is more distinct. Furthermore, the value of $\gamma =1/2$ also aids analytical derivation through simplification of equations which will be discussed in the following sections.

In the region where $z/\varLambda \gt 1$ , note that the current approach is equivalent to using a local friction velocity, $u_{\tau }(z) = u_{\tau }(1 - z/h)^{1/2}$ , to normalise the velocity in (1.1) $(\kappa z\ \mathrm{d} u^+/\mathrm{d}z)$ , while maintaining the linear Businger relation. However, for the remainder of the article, we interpret this factor as a correction to the Businger relation rather than as an estimation of the local friction velocity, and use a mixing length model that consistently accounts for these variations.

Figure 3. The dimensionless velocity gradients, $(\kappa z\ \mathrm{d} u^+/\mathrm{d}z)$ , obtained from DNSs across a range of stratification levels as a function of $z/L$ in panel (a) and $z/\varLambda$ in panel (b). Increasing shades of blue correspond to $Ri_{\tau } = 60,\ 240,\ 720$ at $\textit{Re}_{\tau } = 550$ , the green line corresponds to $Ri_{\tau } = 720$ at $\textit{Re}_{\tau } = 395$ and the red one to $Ri_{\tau } = 600$ at $\textit{Re}_{\tau } = 1000$ . These gradients are compared against the Businger relation, $1+4.7\ \zeta$ , where $\zeta$ is the stability parameter.

4. Heat and momentum flux balance

Following Donda et al. (Reference Donda, van Hooijdonk, Moene, Jonker, van Heijst, Clercx and van de Wiel2015), $\varLambda$ can be expressed in terms of the wall-normal coordinate using the heat and momentum flux distributions. For pressure-driven turbulent channel flows, with walls at fixed temperatures, the mean total heat flux ( $q$ ) is constant across the channel, while the total shear stress ( $\tau$ ) varies linearly

(4.1a) \begin{align} q(z) &= {\alpha }\frac {\mathrm{d}\overline {\theta }}{\mathrm{d}z} - \ \ \ \ \overline {w^{\prime }\theta ^{\prime }} = q_{{w}}, \\[-9pt] \nonumber \end{align}

(4.1b) \begin{align} \tau (z) &= \mu \frac {\mathrm{d} \overline {u}}{\mathrm{d} z} - {\rho _0}{\overline {u^{\prime }w^{\prime }}} = {\tau _{{w}}} \left (1 - \frac {z}{h}\right)\!, \\[9pt] \nonumber \end{align}

where $\alpha$ is the thermal diffusivity, and the subscript ` ${w}$ ’ denotes wall values. Substituting (4.1a ) and (4.1b ) into the definition of $\varLambda$ (1.2) yields

(4.2) \begin{align} \varLambda (z) = L \left (1 - \frac {z}{h}\right)^{3/2}, \end{align}

with $L$ the Obukhov length scale as given in (1.1). Approximating the turbulent flux in (4.1b ) using the mixing length, we can write the stress balance in the well-known form

(4.3) \begin{align} \nu \frac {\mathrm{d} \overline {u}}{\mathrm{d} z} + \ell _m^2 \left |\frac {\mathrm{d}\overline {u}}{\mathrm{d}z}\right | \frac {\mathrm{d} \overline {u}}{\mathrm{d} z} = u_{\tau }^2 \left (1 - \frac {z}{h}\right)\!, \end{align}

where $u_\tau = \sqrt {\tau _w/\rho _0}$ and $\nu = \mu /\rho _0$ . This classical zero-equation model forms the basis for the proposed scaling laws. The mean velocity profile is naturally obtained by integrating (4.3) across the channel.

The mixing length $(\ell _m)$ is governed by the relative influence of buoyancy forces and turbulent shear. It is defined by scaling its value for pressure-driven turbulent channel flows in the neutral limit $(\ell _m^{{N}})$ , with the stability function $\phi$ that accounts for effects of stable stratification

(4.4) \begin{align} \ell _m = \frac {\ell _m^{{N}}}{\phi }, \end{align}

with $\phi$ being a function of $z/\varLambda$ (recall the discussion in § 2). For $\ell _m^{{N}}$ , we adopt the mixing length reported by Pirozzoli (Reference Pirozzoli2014) for channel flows

(4.5) \begin{align} \ell _m^{{N}} = \kappa z\left (1 - \frac {z}{h}\right)^{1/2}. \end{align}

Note that, in (4.4), $\ell _m^{{N}}$ does not depend on the strength of stratification, and (4.5) is used across the entire channel. In contrast, the function $\phi$ accounts for stratification, and the stability correction should be used depending on the value of $z/\varLambda$ , as discussed in § 2. Hence, based on the local balance between viscous and turbulent stresses, and the competing influence of buoyancy and shear, (4.3) can be simplified locally to identify distinct flow regions, as discussed in the next section.

5. Proposed mean velocity scaling

Figure 4 shows the typical distribution of viscous and turbulent shear stresses in a pressure-driven turbulent channel flow. Our analysis highlights five characteristic regions, as illustrated in figure 4(a), based on the distribution of the local stability parameter ( $z/\varLambda$ ) and the relative magnitudes of viscous and turbulent stresses. For convenience of reporting, the regions are hereafter given descriptive names.

Figure 4. (a) Schematic illustrating the distinct regions in a stably stratified turbulent channel flow, classified based on the influence of stratification (quantified by $z/\varLambda$ ) and the relative contributions of viscous stress ( $\tau _{{v}}$ ) and turbulent stress ( $\tau _{{t}}$ ) to the total stress ( $\tau$ ). The identified regions are: viscous sublayer (I), shear-dominated sublayer (II), stratified outer region (III), turbulent-viscous transition layer (IV) and viscous core (V). The illustration corresponds to the case with $Ri_{\tau } = 720$ at $\textit{Re}_{\tau } = 550$ . (b) Stress profiles for various cases. Blue curves: $\textit{Re}_{\tau } = 550$ , $Ri_{\tau } = 0,\ 60,\ 240,\ 720$ (increasing darkness). Green: $\textit{Re}_{\tau } = 395$ , $Ri_{\tau } = 720$ . Red: $\textit{Re}_{\tau } = 1000$ , $Ri_{\tau } = 600$ (Zonta et al. Reference Zonta, Sichani and Soldati2022). (c) Close-up of the profiles close in the near-wall region to illustrate the slow rise of turbulent stresses with increasing stratification. (d) Close-up near the centre of the channel, highlighting the drop-off of turbulent stresses to zero.

Closest to the wall, the classical viscous sublayer (I) (and the buffer layer) is unaffected by buoyancy forces as the strong shear dominates the stratification effects in this region. This is followed by a shear-dominated sublayer (II), where turbulent stresses dominate with minimal stratification effects. Further away from the wall, buoyancy effectively suppresses the larger eddies while the turbulent stresses are still much greater than the viscous stresses – a region termed here as the stratified outer region (III). Beyond lies the turbulent-viscous transition layer (IV), where viscous and turbulent stresses become comparable. Finally, from the channel centre, turbulence may vanish and lead to a laminar region, forming a viscous core layer (V) governed by viscous stresses. We now derive the mean velocity profile in these regions, excluding the well-established $u^+ = z^+$ scaling in the viscous sublayer.

5.1. Shear-dominated sublayer (II) – classical MOST

Following the viscous sublayer, viscous stresses rapidly decrease while turbulent stresses increase, reaching a peak at approximately $z^+ \approx 30$ wall units (seen in figure 4 c). In the shear-dominated sublayer, the local stability parameter $z/\varLambda$ is less than one, indicating that turbulent shear dominates over buoyancy effects. This is typically classified as the ‘surface layer’ in the ABL context, where the turbulent fluxes remain approximately constant and the classical MOST is applicable. In this region, (4.3) simplifies to

(5.1) \begin{align} \ell _m \frac {\mathrm{d} \overline u}{\mathrm{d} z} = u_{\tau }\left (1 - \frac {z}{h}\right)^{1/2}. \end{align}

The mixing length in the shear-dominated sublayer $(\textrm{II})$ , denoted as $\ell _m^{\textrm{II}}$ , follows the formulation in (4.4). In this region, the effect of stable stratification is incorporated through the stability function represented by the linear Businger relation, $\phi = 1 + 4.7 \zeta$ , where $\zeta = z/\varLambda$ corresponds to the local stability parameter (recall figure 3 and the discussion thereabout). Hence, the mixing length is given by

(5.2) \begin{align} \ell _m^{\textrm{II}} = \frac {\ell _m^{{N}}}{\phi } =\frac {\ell _m^{{N}}}{1 + 4.7 \ \zeta } = \frac {\kappa z}{1 + 4.7 \ \zeta }\left (1 - \frac {z}{h}\right)^{1/2}\text{;} \quad \ \ \ \zeta = \frac {z}{\varLambda }. \end{align}

Figure 5 presents a comparison between the mixing length measured from DNS data, $(-\overline {u^{\prime}w^{\prime}})^{1/2}/(\mathrm{d}\overline {u}/\mathrm{d}z)$ , and the expression given by (5.2). A good agreement is observed in the region where turbulent shear dominates (i.e. $z/\varLambda \lt 1$ , see the fine dashed lines). The expression begins to deviate from the DNS profile when buoyancy effects become prominent, indicating a transition to a different regime where an alternative mixing length formulation is required.

Figure 5. Comparison of the mixing length obtained from DNS ( $\ell _m^{\textit{DNS}} = (-\overline {u^{\prime}w^{\prime}})^{1/2}/(\mathrm{d}\overline {u}/\mathrm{d}z)$ , solid lines) with the mixing lengths used in this study. The short dashed lines correspond to the mixing length expression for the shear-dominated sublayer, $\ell _m^{\textrm{II}}$ (5.2), while the long-dashed lines correspond to the mixing length formulation in the stratified outer region, $\ell _m^{\textrm{III}}$ (5.6). Line colours indicate different cases: green for $\textit{Re}_{\tau } = 395$ , $Ri_{\tau } = 720$ ; blue for $\textit{Re}_{\tau } = 550$ , $Ri_{\tau } = 720$ ; and red for $\textit{Re}_{\tau } = 1000$ , $Ri_{\tau } = 600$ .

Substituting the mixing length expression from (5.2) into (5.1) yields the velocitygradient for this region

(5.3) \begin{align} \frac {\mathrm{d} u^+}{\mathrm{d} z^+} = \frac {1}{\kappa z^+}\left (1 + 4.7 \frac {z}{\varLambda }\right)\!, \end{align}

which we have shown in figure 3 (see short dash-dotted lines, $z/\varLambda \lt 1$ ).

Further, upon substituting $\varLambda$ from (4.2) in (5.3) and integrating, we obtain the velocity profile

(5.4) \begin{align} u^+ = \frac {1}{\kappa } \ln (z^+) + \frac {4.7}{\kappa } \frac {2h}{L} \left (1 - \frac {z}{h}\right)^{-1/2} + C_1, \end{align}

where $C_1$ is the integration constant, obtained by matching (5.4) to the viscous sublayer profile. In neutral conditions, this matching typically occurs at approximately $z^+ \approx 11$ where the viscous and turbulent stresses intersect. With increasing $Ri_{\tau }$ , this point shifts farther away from the wall, as the rate at which turbulent stresses rise is reduced, indicating deeper penetration of buoyancy effects in the near-wall region. This is evident in figure 4(c), where the blue curves ( $\textit{Re}_{\tau } = 550$ ) show that the darker shades, representing stronger stratification, rise progressively more slowly. Furthermore, for the same $Ri_{\tau } = 720$ , the case at $\textit{Re}_{\tau } = 395$ (green) shows an even slower rise compared with $\textit{Re}_{\tau } = 550$ . This is because the stronger shear tends to resist the suppression of turbulence by buoyancy forces at higher values of $\textit{Re}_{\tau }$ . The constant used in the present work is described in Appendix A.

Interestingly, the half-channel height, $h$ , appears in (5.4), indicating the influence of an outer length scale on the mean velocity profile near the wall. In the neutral limit ( $L \rightarrow \infty$ ), however, this dependence vanishes, and the classical logarithmic law is recovered. With the introduction of stratification, a new length scale $\varLambda$ emerges, determined by the local stress distribution. This distribution is constrained by the stress boundary conditions at the wall and at the channel centreline, thereby introducing a dependence on $h$ in the formulation.

The shear-dominated sublayer extends from $z^+ \approx 30$ to where stratification effects become prominent relative to the turbulent shear. This transition point can be determined analytically by identifying the location where $z/\varLambda \approx 1$ . Substituting $\varLambda = z$ into (4.2) yields

(5.5) \begin{align} \left (\frac {z_t}{h}\right)^3 + \left (\frac {h^2}{L^2} - 3\right)\ \left (\frac {z_t}{h}\right)^2 + 3 \left (\frac {z_t}{h}\right) -1 = 0, \end{align}

whose real root gives the normalised transition location, $z_t/h$ .

For neutral flows, the transition point $(z_t)$ is located at the channel centreline, indicating the absence of buoyancy effects and the dominance of shear throughout the channel. Consequently, the mixing length given by (5.2) is applied up to the channel centre. This assumes a logarithmic law of the wall across the remainder of the neutral velocity profile, which is a reasonable approximation given the relatively weak wake contribution observed in channel flows. With increasing stratification, the transition point shifts closer to the wall, marking the onset of the stratified outer region.

5.2. Stratified outer region (III)

This region begins approximately at the transition point ( $z = z_t$ ), beyond which turbulent stresses keep decreasing while still dominating the stress budget (figure 4). Yet, the influence of buoyancy forces becomes increasingly significant $(z/\varLambda \gt 1)$ , with large eddying motion being more effectively suppressed. In this region, (5.1) remains valid since turbulent stresses remain much larger than viscous stresses.

However, the mixing length in this region ( $\ell _m^{\textrm{III}}$ ) differs from that in the shear-dominated sublayer, as stratification effects become significant ( $z/\varLambda \gt 1$ ). In this case, $\ell _m^{{N}}$ is normalised using the Businger relation that is corrected by a factor $(1 - z/h)^{1/2}$ . As discussed in § 2, similar corrections with varying exponents have been used in previous studies (e.g. Yokoyama et al. (Reference Yokoyama, Gamo and Yamamoto1979), Sorbjan (Reference Sorbjan1986)) to extend (1.1) across the ABL. Here, the correction is applied directly to the Businger relation and is motivated by the trends observed in figure 3. Thus,

(5.6) \begin{align} \ell _m^{\textrm{III}} = \frac {\ell _m^{{N}}}{1 + 4.7\zeta }\left (1 - \frac {z}{h}\right)^{-1/2} = \frac {\kappa z}{1 + 4.7\zeta }. \end{align}

This mixing length model is also depicted in figure 5 (long dashed lines), agreeing well with the DNS data. The formulation in (5.6) reasonably captures the DNS trend in the region dominated by buoyancy effects.

Substituting (5.6) in (5.1) gives the velocity gradient in the stratified outer region

(5.7) \begin{align} \frac {\mathrm{d}u^+}{\mathrm{d}z} = \frac {1}{\kappa z}\left (1 + 4.7 \frac {z}{\varLambda }\right)\left (1 - \frac {z}{h}\right)^{1/2}, \end{align}

which, after substituting the definition of $\varLambda$ from (4.2) and integrating, yields the velocity profile

(5.8) \begin{align} u^+ = \frac {1}{\kappa }\ln {\frac {z}{h}}-\frac {4.7}{\kappa } \frac {h}{L} \ln {\left (1 - \frac {z}{h}\right)} + \frac {2}{\kappa } \sqrt {1 - \frac {z}{h}} - \frac {2}{\kappa } \ln {\left (1 + \sqrt {1 - \frac {z}{h}}\right)} + C_2, \end{align}

with $C_2$ , being an integration constant. Note that a similar expression was derived for open-channel flows by van de Wiel et al. (Reference van de Wiel, Moene and Jonker2012), Donda et al. (Reference Donda, van Hooijdonk, Moene, Jonker, van Heijst, Clercx and van de Wiel2015) in the context of ABL.

Furthermore, we observe that $C_2$ approximately corresponds to the centreline velocity of neutral flows at a given $\textit{Re}_{\tau }$ , i.e. $C_2 \approx {1}/{\kappa }\ln {\textit{Re}_{\tau }} + 5.2$ . This can be reconciled by noting that, in neutral limit ( $L \rightarrow \infty$ ), the shear-dominated sublayer extends up to the centreline. Consequently, the stratified outer region reduces to a single point at $z = h$ . In this limit, (5.8) yields the centreline velocity of the neutral flow. Note also that, in this case, the buoyancy-related term in (5.8) vanishes, and in the overlap region ( $0 \ll z \ll h$ ) the expression also recovers the classical velocity defect law (Pope Reference Pope2000), demonstrating the consistency of the approach.

The profile described by (5.8) results from integrating only the turbulent stress contribution (5.1). However, farther from the walls, turbulent stresses gradually become comparable to viscous stresses, and the latter must also be taken into account (see also figure 4). Based on observations across a wide range of cases, this transition typically occurs at a location where $z/\varLambda \approx 8$ . This location also marks the upper bound beyond which $z/\varLambda$ ceases to be a relevant parameter, as the turbulent momentum flux is no longer accurately represented by the total stress in the definition of $\varLambda$ . Beyond this, (5.1) is no longer valid, marking the onset of the transition layer. Within the turbulent–viscous transition layer, the viscous and turbulent diffusion play comparable roles. Due to the limited extent of this region, a simple linear blending of the adjacent layers provides a reasonable estimate of the flow velocity here. While this blending approach is not explicitly discussed in the main text, its implementation and resulting composite profile are illustrated in Appendix D and also the accompanying Jupyter notebook.

5.3. Viscous core (V)

Finally, towards the centre of the channel, stratification may become so influential that the eddies are fully suppressed, turbulence ceases to exist and the local flow is laminar. Since only viscous stresses contribute to the momentum balance in this region, (4.3) simplifies and can be integrated, with the centreline velocity ( $u_{{cl}}^+$ ) appearing as the constant of integration

(5.9) \begin{align} \frac {h}{\textit{Re}_{\tau }} \frac {\mathrm{d}u^+}{\mathrm{d}z} = 1 - \frac {z}{h}, \quad u_{{cl}}^+ - u^+ = \frac {\textit{Re}_{\tau }}{2}\left (1 - \frac {z}{h}\right)^2. \end{align}

The distance from the centreline where the viscous and turbulent stresses become comparable can be estimated as the location up to which the mean velocity follows a parabolic profile. A reasonable estimate for this distance, from the channel centre, is where $h/\textit{Re}_{\tau } \sim 2 \ell _m \sqrt {\tau /\tau _w}$ (see Appendix B). This point ( $z^l_{\textrm{V}}$ ) marks the beginning of the viscous core, which extends up to the centreline. If the flow is strongly stratified, then the turbulent stresses vanish earlier, as can be seen in figure 4(d). The estimated distance is given as

(5.10) \begin{align} 1 -\frac {z^l_{\textrm{V}}}{h} = \frac {1}{\sqrt {\textit{Re}_L}}\ \sqrt {\frac {4.7}{2\kappa }}, \end{align}

where $\textit{Re}_{L} = u_{\tau } L/\nu$ is the friction Reynolds number based on the Obukhov length scale.

5.4. Important parameterisations for predicting the mean velocity profile

The velocity profiles ((5.4), (5.8), (5.9)) depend on the governing input parameters $\textit{Re}_{\tau }$ and $Ri_{\tau }$ . Yet, an important parameter naturally emerging throughout the analysis is a Richardson number based on wall fluxes, $Ri_{{w}}=h/L$ (recall 1.1). Hence, to obtain a closed form of the velocity profiles, $Ri_{{w}}$ must be estimated a priori. To this end, we derive a relation (see Appendix C) for it in terms of $\textit{Re}_{\tau }$ and $Ri_{\tau }$ , relying on the Nusselt number scaling proposed by Zonta et al. (Reference Zonta, Sichani and Soldati2022) that was obtained for a range of cases (either simulated at $\textit{Pr} = 0.71$ or rescaled accordingly). The proposed scaling is given by

(5.11) \begin{align} Ri_{{w}} = \frac {h}{L} = \frac {\kappa \ (g/\theta _0)\ q_w}{u_{\tau }^3/h} \sim \ \left (\frac {Ri_{\tau }}{\sqrt {\textit{Re}_{\tau }}}\right)^{2/3}. \end{align}

Figure 6. (a) The proposed scaling to estimate a priori the parameter $Ri_{{w}} = h/L$ . (b) Parametrisation of the deviation of the centreline velocity from the neutral case at same $\textit{Re}_{\tau }$ . Increasing darkness of the marker colours denotes increasing Richardson numbers. The symbols and colours are as indicated in figure 2.

Figure 6(a) shows the agreement of this scaling for the different cases considered in this study. Additionally, the viscous core solution requires a closure for stratification effects on the centreline velocity, $u_{{cl}}^+$ . This is done by modelling its deviation from the neutral case at the same $\textit{Re}_{\tau }$ , given by $u_{{cl,N}}^+ \approx {1}/{\kappa } \ln \textit{Re}_{\tau } + 5.2$ . With stratification, turbulence suppression near the centreline leads to an increased $u_{{cl}}^+$ that we found to be well captured by

(5.12) \begin{align} \Delta u_{{cl}}^+ = u_{{cl}}^+ - u_{{cl,N}}^+ = Ri_{\tau }^{0.63} \textit{Re}_{\tau }^{-0.12}, \end{align}

as shown in figure 6(b). Indeed, $\Delta u_{{cl}}^+$ increases with $Ri_{\tau }$ as expected due to stronger stratification effects, and depends weakly on $\textit{Re}_{\tau }$ , as most of its influence is already accountedfor in $u_{{cl,N}}^+$ .

These parametric relations and ((5.4), (5.8), (5.9)) are the ingredients needed for reconstructing the mean velocity profile across the entire channel height. In the following section, we compare the predictions against our DNS dataset to evaluate the applicability of MOST in this confined flow configuration.

6. Assessment of the proposed scaling with DNS

Figure 7 compares the composite mean velocity profile predicted by the proposed approach against DNS results across several $\textit{Re}_{\tau }$ and $Ri_{\tau }$ cases. The profiles, plotted in linear coordinates from the wall to the channel centreline, provide a global assessment of the model performance. The individual contributions from ((5.4), (5.8), (5.9)) are shown as dashed lines, with the highlighted segments indicating the region over which each expression is valid. The symbols are used to distinguish the different flow regions discussed earlier.

Figure 7. The composite mean velocity profiles obtained from the proposed formulation (dashed lines) are compared with DNS data (solid lines). The markers indicate the boundaries of the different regions, as defined in figure 4. Results are shown for $\textit{Re}_{\tau } = 395$ (greens), $\textit{Re}_{\tau } = 550$ (blues) and $\textit{Re}_{\tau } = 1000$ (reds). The implementation of the mean velocity profile reconstruction and the estimation of the skin-friction coefficient can be found in this notebook.

In addition to this overall comparison, figure 8 offers a closer picture of the predicted and DNS velocity profiles, focusing on the different flow regions. To examine the influence of stratification, results are presented for flows with $Ri_{\tau } = 0, \ 60, \ 240 \ \text{and} \ 720$ at $\textit{Re}_{\tau } = 550$ (shown in increasingly darker shades of blue). Moreover, to compare the Reynolds number effects, results for flows with $Ri_{\tau } = 720$ are presented at $\textit{Re}_{\tau } = 395,\ 550$ (green, dark blue). Since the case at $Ri_{\tau } = 720$ was unavailable for $\textit{Re}_{\tau } = 1000$ , the result of the closest available $Ri_{\tau } = 600$ is shown in dark red.

Figure 8. Mean streamwise velocity profiles in stably stratified turbulent channel flows. The solid lines correspond to DNS data which are compared against (a) (5.4) in the shear-dominated sublayer, (b) (5.8) in the stratified outer region and (c) (5.9) in the viscous core. The predictions are indicated with dashed lines. Blue curves: $\textit{Re}_{\tau } = 550$ , $Ri_{\tau } = 0,\ 60,\ 240,\ 720$ (increasing darkness). Green: $\textit{Re}_{\tau } = 395$ , $Ri_{\tau } = 720$ . Red: $\textit{Re}_{\tau } = 1000$ , $Ri_{\tau } = 600$ . Symbols , and indicate the bounds of different regions. (d) Comparison of skin-friction coefficient, $C_f$ , predicted by the proposed approach (coloured dashed lines) with DNS data (symbols) as a function of $Ri_{\tau }$ . The grey dashed line indicates the $C_f \sim Ri_{\tau }^{-1/3}$ scaling shown in Zonta et al. (Reference Zonta, Sichani and Soldati2022). The inset shows the percentage error ( $\epsilon$ ) between the model and DNSs at different $Ri_{\tau }$ across all cases considered. The shaded region indicates a $\pm 2 \,\%$ error band. The implementation of the mean velocity profile reconstruction and the estimation of the skin-friction coefficient can be found in this notebook.

In the shear-dominated sublayer (II), it can be observed that (5.4) accurately predicts the mean velocity profile with discrepancies appearing only for strongly stratified flows; for instance, $Ri_{\tau } = 720$ at $\textit{Re}_{\tau } = 395$ (figure 8 a). This is due to insufficient scale separation in the flow, particularly at the imposed stratification. When $\textit{Re}_{\tau }$ is not sufficiently large and stratification is strong, the smallest eddies in the flow are not small enough to remain unaffected by buoyancy. As a result, a substantial portion of the turbulence spectrum is suppressed. The estimated extent of the shear-dominated sublayer ( $z_t$ ) is indicated for different cases with $\blacktriangle$ . If $z_t^+ \lt 150$ , then the buoyancy effects are perceived too close to the wall, disrupting the near-wall region of the flow and thus violating the assumption of a logarithmic velocity profile in this region.

In the stratified outer region (III), predictions with (5.8), with $C_2 \approx 1/\kappa \ln {\textit{Re}_{\tau }} + 5.2$ , agree well with all the cases except for the case at $\textit{Re}_{\tau } = 395$ (figure 8 b), for the same reason as described above. While the equation effectively reproduces the general shape of the velocity profile, we observe a shift for low $\textit{Re}_{\tau }$ cases. The point marked with $\blacktriangledown$ ( $z_{\textrm{III}}^u$ ) in figure 8(b) indicates the end of the stratified outer region and the onset of the turbulent-viscous transition layer, where the influence of viscous stresses becomes comparable to that of turbulent stresses. Finally, figure 8(c) also confirms the agreement of the parabolic velocity profile with DNS data in the viscous core (V). The velocity profile in this region follows a parabolic form parametrised solely by $\textit{Re}_{\tau }$ . Hence, a distinct parabolic curve is obtained for each $\textit{Re}_{\tau } = 395,\ 550,\ 1000$ . The extent to which the DNS profiles follow the parabolic form, measured from the channel centre, increases with stratification. This is evident for $Ri_{\tau } = 0,\ 60,\ 240,\ 720$ (increasing shades of blue) at $\textit{Re}_{\tau } = 550$ . Specifically, the $Ri_{\tau } = 720$ case follows the parabolic profile over a larger portion of the core compared with $Ri_{\tau } = 60\ \text{and}\ 240$ , while $Ri_{\tau } = 0$ does not exhibit this behaviour in the centre. The boundary marking the approximate location beyond which the velocity profiles begin to diverge from the parabola is indicated by the symbol . We observe that, for $Ri_{\tau } =720$ at $\textit{Re}_{\tau } = 395,\ 550$ and $Ri_{\tau } = 600$ at $\textit{Re}_{\tau } = 1000$ , this is estimated reasonably well. However, at lower $Ri_{\tau }$ values, corresponding to $\textit{Re}_{\tau } = 550$ , this boundary is estimated to be at a slightly larger distance from the channel centre than the actual location where the velocity profiles begin to diverge from the parabola. Nevertheless, this difference is visually exaggerated due to the logarithmic horizontal axis, and the estimate is acceptable for practical purposes.

The mean velocity profile obtained from the proposed approach can be integrated to obtain bulk velocity ( $u_b$ ), and consequently estimate the skin-friction coefficient ( $C_f = 2\tau _w/\rho u_b^2 = 2/{u_b^+}^2$ ). Figure 8(d) compares these estimates with values obtained from DNS, with the scaling $C_f \sim Ri_{\tau }^{-1/3}$ reported by Zonta et al. (Reference Zonta, Sichani and Soldati2022) also shown for reference. While this relation may hold in a narrower parameter range, it tends to deviate from DNS results over the wider parameter space explored in this study. Integrating the velocity profile obtained from the proposed approach yields consistently accurate results, illustrating the robustness of the present framework. The relative errors, defined as $\epsilon = (C_f^{\textit{DNS}} - C_f)/C_f^{\textit{DNS}}$ , are generally within a $2\,\%$ margin, with the exception of the high $Ri_{\tau }$ cases at the lowest $\textit{Re}_{\tau } = 395$ , where deviations up to $4\,\%$ are observed. This slightly larger error is, once more, an expected consequence of the lack of scale separation in these cases.