Numerical Simulations of Scalar Transport on Rough Surfaces

Abstract

Numerical simulations provide unfettered access to details of the flow where experimental measurements are difficult to obtain. This paper summarises the progress achieved in the study of passive scalars in flows over rough surfaces thanks to recent numerical simulations. Townsend’s similarity applies to various scalar statistics, implying the differences due to roughness are limited to the roughness sublayer (RSL). The scalar field exhibits a diffusive sublayer that increasingly conforms to the roughness surface as $k_s^{+}$ or $Pr$ increase. The scalar wall flux is enhanced on the windward slopes of the roughness, where the analogy between momentum and scalar holds well; the momentum and scalar fields, however, have very different behaviours downwind of the roughness elements, due to recirculation, which reduces the scalar wall flux. Roughness causes breakdown of the Reynolds analogy: any increase in $St$ is accompanied by a larger increase in $c_f$. A flattening trend for the scalar roughness function, $\Delta {\Theta }^{+}$, is observed as $k_s^{+}$ increases, suggesting the possibility of a scalar fully rough regime, different from the velocity one. The form-induced (FI) production of scalar fluctuations becomes dominant inside the RSL and is significantly different from the FI production of turbulent kinetic energy, resulting in notable differences between the scalar and velocity fluctuations. Several key questions remain open, in particular regarding the existence of a fully rough scalar regime and its characteristics. With the increase in $Re$ and $Pr$, various quantities such as scalar roughness function, the dispersive fluxes, FI wall flux, etc., appear to trend towards saturation. However, the limited range of $Re$ and $Pr$ achieved by numerical simulations only allows us to speculate regarding such asymptotic behaviour. Beyond extending the range of $Re$ and $Pr$, systematic coverage of different roughness types and topologies is needed, as the scalar appears to remain sensitive to the geometrical details.

1. Introduction

1.1. Scalar Transport

The transport of scalar quantities in turbulent flows is of great importance in many engineering applications, as well as a central phenomenon in nature. The transfer of internal energy or enthalpy affects applications such as heat exchangers, gas turbines, and nuclear reactors, while mass transport governs the dispersion of particles in the atmosphere, or hazardous material release. Scalar transport adds additional transport equations with a new parameter, the Prandtl number $Pr=\nu /\alpha$ or Schmidt number $Sc=\nu /D$, where $\nu$ is the fluid kinematic viscosity, and $\alpha$ or D are the conductivity or diffusivity of the scalar, respectively. Scalars are defined as “passive” when their effect on the fluid or flow is negligible, that is, they do not change significantly the fluid viscosity, density, etc.

1.2. Stanton Number and the Reynolds Analogy

Scalar transport (especially for heat transfer applications) has been studied for over a century. Reynolds [1] suggested that momentum and heat in a fluid are transferred in the same way and, therefore, proportionality must exist in geometrically similar systems; this relation is commonly referred to as the “Reynolds analogy”. Generally, it applies to gases where the fluid viscosity and scalar diffusivity are similar, and is based on the notion that the same turbulent eddies transport both scalar quantities and momentum [2,3,4]. The Reynolds analogy is based on the formal similarity between the transport equations for momentum and scalar quantities, as long as the pressure gradient is insignificant and no additional forces are present [1,2]. The skin-friction coefficient and the scalar-flux coefficient, $St$, are defined as

$$C_f=\frac{2{\tau }_w}{\rho U_b^2};St=\frac{q_w}{\rho c_p{\Theta }_m{U}_b}=\frac{1}{{\Theta }_m^{+}{U}_b^{+}},$$

(1)

where ${\tau }_w$ is the wall stress, $q_w$ the wall flux, $\rho$ the fluid density, $c_p$ the specific heat capacity, ${\Theta }_m$ is the mixed-mean scalar, and $U_b$ is the average velocity. The Reynolds analogy implies that the ratio $2C_f/St$, known as the Reynolds analogy factor, $RA=U_b^{+}/{\Theta }_m^{+}$, is constant and equal to one for $Pr$ near unity [2,3,4,5]. While this analogy also applies to mass transfer (where the transfer coefficient is the Sherwood number, $Sh$) or any other type of scalar transport, the rest of the paper will refer to heat transfer and temperature. It is to be understood that other scalars can be treated in the same way.

The analogy may also hold when $Pr\ne 1$; in this case, $RA\ne 1$ (but remains constant) [4,6]. Classic approximations of $RA$ for $Pr\le 5$ are power laws of the form $RA\sim {Pr}^m$, and various values of m have been suggested [3,7,8]. For larger $Pr$, a more complicated formulation is necessary [4,8,9].

1.3. Effects of Roughness on Scalar Transport

The Reynolds analogy has been extensively verified for smooth surfaces [10,11,12,13,14]. However, hydraulically smooth behaviour can only be achieved at sufficiently low $Re$ [15,16]. At higher $Re$, when the roughness scale becomes comparable to the viscous sublayer, any surface becomes rough. Roughness has many effects on momentum transfer; the main ones are (1) an increase in the drag [17,18]; (2) the amplification of the wall-normal and spanwise fluctuations, at the expense of the streamwise ones, resulting in decreased flow anisotropy [19,20,21,22]; and (3) the breakup of the near-wall structures and modifications of the near-wall turbulence generation cycle [23,24]. The increased wall stress results in a downward shift of the logarithmic layer, the “roughness function”, $\Delta U^{+}$ (here, + denotes wall units, i.e., quantities normalised using $u_{\tau }$ and $\nu$, where $u_{\tau }=\sqrt{{\tau }_w/\rho }$ is the friction velocity). Townsend [25] proposed a similarity hypothesis stating that, at sufficiently high Reynolds numbers, the only role of the region below the roughness crest is to set the length and velocity scales for the outer flow. Approximately 3–5 roughness heights above the crest, the turbulence statistics thus normalised collapse with the corresponding smooth-wall case [15,16,26,27]. This region is known as the roughness sublayer, RSL. Reviews of this topic can be found in Flack and Schultz [15], Chung et al. [16], Raupach et al. [23], Jiménez [28].

There are many similarities between the effects of roughness on momentum and scalar transport. $St$ increases, and the mean scalar profile, which for smooth wall presents a log-law behaviour, is also shifted down by a scalar roughness function, $\Delta {\Theta }^{+}$. Townsend’s similarity hypothesis was found to apply to the scalar variance and turbulent fluxes [29,30,31,32,33], to the terms in their budget [29,31,34], to skewness and flatness [29], and to some of the higher-order moments [29]. However, experimental studies [9,35] showed that when roughness is introduced $St$ remains sensitive to scalar diffusivity and decreases with $Re$, while $c_f$ approaches a constant value in the fully rough regime [9,15,17,35]. This implies the breakdown of the Reynolds analogy when the surface is sufficiently rough. Since the roughness elements create stagnation points and recirculation regions, significant pressure gradients occur that break the formal similarity between momentum and scalar transport equations. Away from the wall, on the other hand, the scalar transport is dominated by the turbulent transport, so that Townsend’s similarity hypothesis [25] also applies to the scalar.

1.4. Numerical Simulations of Rough-Wall Scalar Transport

Numerical simulations have recently begun to contribute to our understanding of scalar-transport physics over rough surfaces. Direct numerical simulations (DNSs) (which requires that all the scales of motion are resolved) and large-eddy simulations (LESs) (in which only the large, energy-carrying eddies are resolved) have been used to highlight the flow physics and validate turbulence models for the Reynolds-averaged Navier–Stokes (RANS) equations. Early simulations were limited to simple roughness shapes, such as two-dimensional (2D) square bars, ribs with various orientations, circular rods, wavy walls, or ridges [36,37,38,39,40,41,42]. Direct and large-eddy simulations have also been used to study three-dimensional roughness [5,30,31,32,34,41,43,44,45,46,47,48,49,50,51,52,53,54]. These studies considered $Pr=0.2–6$ and $k_s^{+}$ up to 700, where $k_s^{+}$ is the equivalent sand-grain size (that is, the sand-grain height from Nikuradse’s experiments [17] for which the same roughness function, $\Delta U^{+}$, value as the present surface is achieved); however, most of them were focused on $Pr$ near unity (in particular $0.7$ and $1.0$), and $k_s^{+}\le 400$.

Table 1. Summary of recent numerical simulations of wall-resolved passive scalar transport with roughness.

Study	Max ${\mathit{k}}_{\mathit{s}}^{+}$ I	$\mathit{Pr}$ or $\mathit{Sc}$	Roughness Type
Doosttalab et al. [29]	11	$0.7$	Grit blasted
Stroh et al. [42]	12 I	$0.7$	Streamwise-aligned triangular ridges
Peeters [48]	52 II	2.0–6.0	Grit blasted
Peeters and Sandham [46]	102	$1.0$	Grit blasted
Hantsis and Piomelli [31]	110	0.7–1.4	Sand grain
Yang et al. [55]	120	$0.7$	Ice accretion, sandpaper, combustion chamber deposits
Miyake et al. [36]	130 I	$0.7$	Sand grain III, 2D ribs
Jelly et al. [5]	130	$0.7$	X-ray computed tomography
Orlandi et al. [43]	200 I	$1.0$	Regularly arranged protrusions
Kuwata [52]	200	$1.0$	Distributed hyperbolic elements
Nagura et al. [56]	230	$0.7$	Sinusoidal wall
Forooghi et al. [51]	300	$0.7$	Real (scanned) surfaces, various distributed protrusions
Zhong et al. [45]	300	0.5–2.0	3D sinusoidal
Rowin et al. [57]	340	$0.7$	3D sinusoidal
MacDonald et al. [30,44]	380 IV	$0.7$	3D sinusoidal
Hantsis [32]	485	0.2–2.0	Sand grain
Leonardi et al. [39]	600 I	$1.0$	2D circular and square bars V
Zhang et al. [58,59]	600 I	$0.7$	2D sinusoidal wall
Garg et al. [54]	633	$0.7$	Additive manufacturing scanned surfaces
Kuwata [60]	700	$1.0$	Distributed hyperbolic elements
Zhang et al. [58]	700	$0.7$	Distributed hyperbolic elements

^I If not explicitly reported, value is derived from the given values of $\Delta U^{+}$. ^II The DNS conducted in the study was limited; however, a wider range of $k_s^{+}$ was considered by using results from other studies. ^III Not explicitly solved but represented by a line force emanating from the wall. ^IV This value was later revised to 233 by [45]. ^V Only on one side; other wall is smooth.

contains a non-exhaustive list of key studies.

1.5. Scope of This Paper

The focus of this paper is on studies that resolve the near-wall region and provide physical understanding of the flow behaviour. First, we will present the formulation of the models discussed here, and then we will discuss some of the novel findings obtained from numerical simulations. A list of open questions and some concluding remarks will follow.

2. Problem Description

The problem at hand can be summarised as follows: for a given roughness topology, knowing the velocity and pressure fields, we would like to predict the behaviour of the scalar field and its statistics. A particular question is whether a regime equivalent to the “fully rough” regime exists for the scalar quantity and, if so, how it manifests itself. Since Townsend’s similarity hypothesis extends to scalar transport [29,31,44,46,52], the focus of this review will be the roughness sublayer (RSL), which extends a few roughness heights, k (typically 3–5) above the roughness crest. Most of the studies considered focus on canonical flows: plane channels or flat-plate boundary layers. A passive scalar is then transported by the flow, with the rough boundaries either kept at a constant scalar value, or providing a constant and uniform flux.

2.1. Governing Equations

In incompressible flows, conservation of mass, momentum, and passive scalar can be written as

$$\begin{array}{cc}\hfill \frac{\partial u_k}{\partial x_k}& =0\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\partial u_i}{\partial t}+\frac{\partial \left(u_i{u}_k\right)}{\partial x_k}& =-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_k}\left(2\nu S_{ik}\right)\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{\partial \theta }{\partial t}+\frac{\partial \left(\theta u_k\right)}{\partial x_k}& =\frac{\partial }{\partial x_k}\left(\alpha \frac{\partial \theta }{\partial x_k}\right)+Q\hfill \end{array}$$

(4)

where $x_i$, for $i=1-3$ (or x, y, and z), are the Cartesian coordinates in the streamwise, wall-normal and spanwise directions, respectively. $u_i$ (or u, v, and w) are the velocity components in the Cartesian directions. $p=P/\rho$ is the hydrodynamic pressure divided by the density, $\nu$ is the kinematic molecular viscosity (assumed constant), and $S_{ij}$ is the strain-rate tensor:

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

(5)

In Equation (4), $\alpha$ is the molecular diffusivity and $\theta =T-T_w$, where T is the transported scalar and $T_w$ is its value at the wall. Q is a source term used to facilitate the use of periodic boundary conditions. For a constant-flux condition,

$$Q=u\frac{d{T}_w}{dx}=\frac{q_w}{\rho c_p\delta }\frac{u}{U_b},$$

(6)

where $c_p$ is the specific heat, $\delta$ is the boundary layer thickness (or channel half-width), and $q_w$ is the prescribed scalar wall flux. For the constant temperature (Dirichlet) boundary condition, $Q=0$. A discussion of this approach and its limitations can be found in [61].

2.2. Averaging Process and Mean Equations

Roughness causes the time-averaged fields to be spatially inhomogeneous. Consequently, three types of averages can be defined: the usual time average, and two types of spatial averages, the superficial average, $\langle f{\rangle }_s$, and the intrinsic one, $\langle f\rangle$. The intrinsic average is only performed over the fluid portion of the averaging region, while the superficial one is performed over the averaging region, including both fluid and solid [23,62,63]. If the averaging region does not contain any solid portion, both averages are identical and correspond to the spatial mean of the averaged quantity; this framework allows for consistent comparison between averages of smooth- and rough-wall cases. Using the intrinsic and time averages, a triple decomposition can be defined. Any variable can be written as

$$f(x_i,t)=\overline{f}\left(x_i\right)+f^{\prime }(x_i,t)=\langle \overline{f}\rangle +\tilde{f}\left(x_i\right)+f^{\prime }(x_i,t).$$

(7)

Here, $\langle \overline{f}\rangle$ is the double-averaged value of f, $\tilde{f}$ represents the spatial inhomogeneity of the time-averaged f, and is known as the dispersive or wake component, and $f^{\prime }$ is the stochastic part. The differentiation of the averaged quantities introduces additional terms in the governing equations. The double-averaged momentum and passive-scalar equations are [23,63]

$$\begin{array}{cc}\hfill -\langle \overline{\Pi }{\rangle }_s=& -\langle \overline{u_j}{\rangle }_s\frac{\partial \langle \overline{u_i}\rangle }{\partial x_j}-\frac{\partial }{\partial x_j}\langle \overline{u_i^{\prime }{u}_j^{\prime }}{\rangle }_s-\frac{\partial }{\partial x_j}\langle {\tilde{u}}_i{\tilde{u}}_j{\rangle }_s+\frac{\partial }{\partial x_j}\left(\nu \frac{\partial \langle {\overline{u}}_i{\rangle }_s}{\partial x_j}\right)\hfill \\ \hfill & +\langle \frac{\partial }{\partial x_j}\left(\nu \frac{\partial {\tilde{u}}_i}{\partial x_j}\right){\rangle }_s-\langle \frac{\partial \tilde{p}}{\partial x_i}{\rangle }_s\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill -\langle \overline{Q}{\rangle }_s=& -\langle \overline{u_j}{\rangle }_s\frac{\partial \langle \overline{\theta }\rangle }{\partial x_j}-\frac{\partial }{\partial x_j}\langle \overline{{\theta }^{\prime }{u}_j^{\prime }}{\rangle }_s-\frac{\partial }{\partial x_j}\langle \tilde{\theta }{\tilde{u}}_j{\rangle }_s\hfill \\ \hfill & +\frac{\partial }{\partial x_j}\left(\alpha \frac{\partial \langle \overline{\theta }{\rangle }_s}{\partial x_j}\right)+\langle \frac{\partial }{\partial x_j}\left(\alpha \frac{\partial \tilde{\theta }}{\partial x_j}\right){\rangle }_s.\hfill \end{array}$$

(9)

The third term on the RHS of Equation (8) is the dispersive stress, while the last two terms are the viscous and pressure fluxes; integrating them over the RSL yields the viscous and form drag, respectively. In Equation (9), conversely, the third and fifth terms on the RHS are the dispersive flux and the form-induced (FI) flux, sometimes referred to as “wall-interaction term” [33,60].

3. Achievements

3.1. Instantaneous Scalar Field

Numerical simulations provide unfettered access to details of the flow field, both for momentum and scalar. Examining the instantaneous and time-averaged scalar field is the most direct way to understand the effects of roughness on scalar quantities.

The key differences between the scalar and momentum fields occur near and around the roughness elements, which frequently have stagnation points in front and recirculation regions behind them, leading to significant pressure gradients. As is well known, in the fully rough regime closely packed roughness elements “shelter” each other [62], and the pressure drag (which is weakly dependent on $Re$) is much larger than the $Re$-dependent viscous one. This behaviour has been observed in many investigations [30,31,32,39,41,44,46,52].

Figure 1 shows the instantaneous velocity (a, e) and scalar contours (b–d, f–h) for a range of $Pr$ and $k_s^{+}$. It is apparent that, as $k_s^{+}$ increases, the scalar gradients are confined to a thinner and thinner layer that closely follows the geometry [30,44]. This is demonstrated by the thickness of the regions covered by $u^{+}\le 5$ or ${\theta }^{+}\le 5$, respectively, which indicate the extent of viscous and diffusive sublayers.

The same effect is observed when $Pr$ is increased, and the low diffusivity also acts to confine the gradients to a thin layer. For the lower values of $Pr$ and $k_s^{+}$, however, we still observe some effect of the sheltering and of the recirculation regions. When the convection behind the roughness elements is less vigorous, the scalar forms regions where it is well mixed, where the small gradient limits the local transfer rate of the scalar quantity [30,46,51]. The much thicker recirculation regions behind the roughness protrusions act as “resistance” to the scalar wall flux [30,32,45,51].

In a smooth-wall case, an increase in $Re$ affects the thickness of both viscous and diffusive sublayers (and the boundary layer thickness itself). For a fully rough case, on the other hand, the increase in $Re$ does not affect the viscous sublayer while still thinning the diffusive sublayer; therefore, the mechanisms that govern the viscous sublayer in the smooth-wall case do not apply directly to the rough-wall case. It is then clear that the effects of $Re$ and $Pr$, while similar, work separately from each other—a Péclet-number scaling $P{e}_{\tau }=R{e}_{\tau }Pr$, commonly used in the smooth-wall case, is not directly applicable in the presence of roughness [34,45,57].

Repeated impingement of bulk fluid on the exposed (i.e., not sheltered by upstream elements) windward side of the roughness produces an attached thin diffusive sublayer over the windward surface [45,57]. Zhong et al. [45] found that the ratio between the local diffusive and viscous sublayer thicknesses ${\delta }_{\alpha }/{\delta }_{\nu }\propto {Pr}^{-1/3}$ at the crests, supporting an argument originally put forward by Owen and Thomson [35] and Yaglom and Kader [64] that “smooth-wall-like“ boundary-layer behaviour is seen at the local roughness crests. The diffusive sublayer is observed to be thinnest on the exposed windward side and around the crest, resulting in sharp scalar gradients that enhance the scalar wall flux [45,46]. Both viscous and diffusive sublayers at the roughness crests were shown to be thinner compared with an equivalent smooth-wall case (under the same $R{e}_{\tau }$ and $Pr$) [45,57]; however, as $k_s^{+}$ increased, their thickness approached that of the smooth-wall case (this was shown for $Pr=0.5–2.0$; it is assumed that this behaviour is maintained at higher $Pr$) [30,32,45].

Conversely, on the leeward side, the diffusive sublayer is much thicker than the equivalent smooth-wall one [45,57]. This can be seen best on the leeward side right behind the crest in Figure 1, where the diffusive sublayer can sometimes be the same thickness as the roughness elements themselves, and this is further supported by the effective reduction in the scalar wall flux in the recirculation regions. On the windward side, the rough-wall diffusive sublayer is analogous to a contorted, smooth-wall diffusive sublayer (at sufficiently high $k_s^{+}$) [44,45]; a similar analogy cannot be made on the leeward side. The variation of the diffusive-sublayer thickness (compared with an equivalent smooth-wall case) results in enhanced scalar transport in some regions but attenuation in others; the integral of all of these results in a net increase or decrease in scalar transport relative to the smooth-wall case.

The recirculation affects the fluid that impinges on the windward slope, drawing it down and backwards (opposite to the direction of the bulk flow) towards the base of the roughness [44,46]. The lack of a pressure-equivalent mechanism for the scalar results in reduced advection of the scalar away from the troughs and, therefore, results in some accumulation of the scalar in “pockets” [45,46]. These “pockets” extend into sheltered regions where the velocity is much smaller, or even negative [44,46]; their extent seems to be well correlated with the location where $\overline{u}$ becomes zero [45]. This mechanism generates larger scalar gradients locally, which partially compensates the reduced scalar transfer in the recirculation regions (more details in Section 3.3).

3.2. Mean Profiles and the Scalar Roughness Function

Mean velocity and scalar profiles have been extensively measured experimentally ([11,12,13,17,23,26,27,28,35,62,65,66,67,68,69] and others). These studies produced empirical correlations and phenomenological models that, in turn, led to predictions of the scalar roughness function, $\Delta {\Theta }^{+}$, in the fully rough regime.

Disagreement amongst these models resulted in very different predicted behaviours of $\Delta {\Theta }^{+}$ when $k_s^{+}\to \infty$. The underlying physical assumptions that underpin these models were difficult to evaluate, due to a lack of high-fidelity data, in particular for dense roughness. Numerical simulations provide a tool to test the underlying physical assumptions, and to sweep systematically across a range of parameters (geometry, Prandtl number, Reynolds number, etc.), which may be difficult to achieve experimentally. In particular, the minimal-channel approach facilitates a relatively low-cost way to achieve such parameter sweeps.

Figure 2 show the profile of the mean scalar as $Re$ and $Pr$ are varied. Similar to the mean velocity, the mean scalar also shows a logarithmic region, which can be expressed as:

$${\Theta }^{+}=\langle \overline{\theta }{\rangle }^{+}=\frac{1}{{\kappa }_{\theta }}log\left[{\left(y-d\right)}^{+}\right]+A_{\theta }-\Delta {\Theta }^{+}$$

(10)

where ${\kappa }_{\theta }\approx 0.46$ [31,70] and $A_{\theta }\left(Pr\right)$ is the smooth-wall intercept. Zhong et al. [45] report other equivalent forms of ${\Theta }^{+}$. $A_{\theta }$ increases with $Pr$ [3,70], resulting in an upwards shift of the logarithmic region. In contrast, the increase in $Re$ (and therefore $k_s^{+}$) results in a downward shift of the logarithmic region.

The downwards shift of the logarithmic region between the smooth- and rough-wall cases is quantified by the scalar roughness function, $\Delta {\Theta }^{+}$, which is equivalent to the Clauser–Hama roughness function, $\Delta U^{+}$. $\Delta {\Theta }^{+}$ depends on $Pr$ and various geometrical roughness parameters, as illustrated in Figure 2. In Figure 2a, $Pr=1.0$ is fixed and $k^{+}$ is increased (by varying $Re$); $\Delta {\Theta }^{+}$ increases with $k^{+}$, but seems to approach an asymptotic state where all curves collapse at high $k^{+}$. This trend was also noticed in other types of roughness [5,30,34,45,57]. In Figure 2b, the roughness height is fixed and $Pr$ is varied. Again, there is an upwards shift of the logarithmic layer with $Pr$, which is much larger in the smooth-wall case [34,64]. As a consequence, the net shift in $\Delta {\Theta }^{+}$ between the smooth- and rough-wall cases (for each $Pr$) increases with $Pr$ [34,48,57].

Figure 2. Variation of the mean (double-averaged) scalar profiles with $R{e}_{\tau }$ and $Pr$. (a) $Pr=1.0$ at different $R{e}_{\tau }$; (b) $Pr$ between $0.2$ and $2.0$ at $R{e}_{\tau }=4200$ ($k_s^{+}=265$). Dashed line in (a) indicates the smooth-wall curve for $Pr=1.0$. Reprinted with permission from Hantsis [32].

Figure 2. Variation of the mean (double-averaged) scalar profiles with $R{e}_{\tau }$ and $Pr$. (a) $Pr=1.0$ at different $R{e}_{\tau }$; (b) $Pr$ between $0.2$ and $2.0$ at $R{e}_{\tau }=4200$ ($k_s^{+}=265$). Dashed line in (a) indicates the smooth-wall curve for $Pr=1.0$. Reprinted with permission from Hantsis [32].

Several empirical models have been proposed to relate $\Delta {\Theta }^{+}$ to the equivalent sand-grain roughness, Prandtl number, and scalar von Kármán constant [3,9,35,64,69]. These models result in disparate predictions of $\Delta {\Theta }^{+}$, generally suggesting a monotonic decrease at higher $k_s^{+}$. These studies also cannot agree on the value of $k_s^{+}$ at which this decrease starts or on the decrease rate.

A different trend was observed by MacDonald et al. [30,44], who suggested a flattening of $\Delta {\Theta }^{+}$ and the possibility of an asymptotic value $\Delta {\Theta }_{\mathrm{FR}}^{+}$ for $k_s^{+}\to \infty$. This behaviour is consistent with the observation of a thin diffusive sublayer that conforms to the topography in the fully rough regime. Other studies [5,31,41,46,47,52,53,57,71], covering various types of roughness and a wider range of $Pr$, support this trend (Figure 3).

Zhong et al. [45], however, found that the scalar roughness function is largest at $k^{+}\approx 40$, and then gradually decreases. They suggest that this trend would lead to a decrease of the wall flux, which contradicts the current understanding of rough-wall scalar transfer [72,73]. This discrepancy may result from the use of the minimal-channel approach, in which the mean profiles away from the wall may be artificially altered.

Figure 3. Scalar roughness function, $\Delta {\Theta }^{+}$, vs. $k_s^{+}$ from selected studies; (a) shows the effect of $Pr$ and (b) the effect of the frontal solidity, ${\lambda }_f$. For (a): Colours and lines indicate the Prandtl number: $Pr=0.2$; $Pr=0.5$; $Pr=0.7$; $Pr=1.0$; $Pr=2.0$. Markers indicate: ● Hantsis [32]; ▼ Zhong et al. [45]; ▲ Rowin et al. [57]; ◄ MacDonald et al. [44]; ■ Jelly et al. [5]; ◆ Peeters and Sandham [46], Peeters [48]; ★ De Maio [74]; 🟌 Kuwata [60]. Grey dash–dot line: The velocity roughness function, $\Delta U^{+}$, from Nikuradse’s sandpaper experiment [17]. For (b): $Pr=0.7$, data from Rowin et al. [57] for $\Lambda =0.09$; $\Lambda =0.18$; $\Lambda =0.36$. Reprinted with permission from [57].

Figure 3. Scalar roughness function, $\Delta {\Theta }^{+}$, vs. $k_s^{+}$ from selected studies; (a) shows the effect of $Pr$ and (b) the effect of the frontal solidity, ${\lambda }_f$. For (a): Colours and lines indicate the Prandtl number: $Pr=0.2$; $Pr=0.5$; $Pr=0.7$; $Pr=1.0$; $Pr=2.0$. Markers indicate: ● Hantsis [32]; ▼ Zhong et al. [45]; ▲ Rowin et al. [57]; ◄ MacDonald et al. [44]; ■ Jelly et al. [5]; ◆ Peeters and Sandham [46], Peeters [48]; ★ De Maio [74]; 🟌 Kuwata [60]. Grey dash–dot line: The velocity roughness function, $\Delta U^{+}$, from Nikuradse’s sandpaper experiment [17]. For (b): $Pr=0.7$, data from Rowin et al. [57] for $\Lambda =0.09$; $\Lambda =0.18$; $\Lambda =0.36$. Reprinted with permission from [57].

Figure 3a shows the variation of $\Delta {\Theta }^{+}$ with $k_s^{+}$ at different $Pr$. One commonly observed feature is $\Delta {\Theta }^{+}<\Delta U^{+}$ at sufficiently high $k_s^{+}$. For $Pr>1.0$, there exists a finite range of $k_s^{+}$ for which the opposite is true; the value of $k_s^{+}$ for which $\Delta {\Theta }^{+}$ intersects the velocity asymptote increases with $Pr$. Since higher values of $\Delta U^{+}$ and $\Delta {\Theta }^{+}$ indicate enhancement of momentum and scalar transport, respectively, the range in which $\Delta U^{+}<\Delta {\Theta }^{+}$ is of significant interest.

The implication of an asymptotic value $\Delta {\Theta }_{\mathrm{FR}}^{+}$ for all $Pr$ is the possibility that an “ultimate regime”, a scalar equivalent of the fully rough regime, exists, characterised by $\Delta {\Theta }^{+}\to \Delta {\Theta }_{\mathrm{FR}}^{+}$. Peeters [47,48] suggested a model for $\Delta {\Theta }_{\mathrm{FR}}^{+}$, noting that the flattening of $\Delta {\Theta }^{+}$ at sufficiently high $k_s^{+}$ was reproduced when only the viscous–convective and diffusive ranges were considered; otherwise an increase in $\Delta {\Theta }^{+}$ would be observed. Peeters [48] concluded that the flattening must be the result of the interaction between the viscous–convective and diffusive scales with the rough wall.

$\Delta {\Theta }_{\mathrm{FR}}^{+}$ increases with $Pr$ for a fixed roughness geometry [31,32] (Figure 3a). For a fixed $Pr$, it is sensitive to topological details such as the frontal solidity, ${\lambda }_f$ (the ratio between the frontal area of the roughness and the total plan area), as shown in Figure 3b. $\Delta {\Theta }_{\mathrm{FR}}^{+}$ also increases with ${\lambda }_f$, unlike $\Delta U^{+}$, which saturates. Figure 3a also raises the possibility of $\Delta {\Theta }_{\mathrm{FR}}^{+}$ clustering with $Pr$, even for different roughness types and topologies, suggesting stronger dependence on $Pr$ and weaker dependence on the roughness details.

Because of the limited range of $Pr$ and $k_s^{+}$ achieved by numerical studies, we can only speculate on the behaviour of $\Delta {\Theta }^{+}$ at higher $k_s^{+}$, and further investigations are required to clarify the issue. However, it is expected that at the ultimate regime the enhancement of the scalar transport will saturate and will be accompanied by increasing drag.

MacDonald et al. [30] and Rowin et al. [57] analytically derived a relation between $\Delta {\Theta }^{+}$ and the Stanton number, $St$:

$$\begin{array}{cc}\hfill \Delta {\Theta }^{+}& =\sqrt{\frac{2}{c_{f0}}}\left(\frac{1}{{St}_0}-\frac{1}{{\kappa }_{\theta }\kappa }\right)-\sqrt{\frac{2}{c_f}}\left(\frac{1}{St}-\frac{1}{{\kappa }_{\theta }\kappa }\right)\hfill \\ \hfill & =\underset{\mathrm{Term}1}{\underbrace{A_{\theta }+\frac{\kappa }{{\kappa }_{\theta }}\left(\Delta U^{+}-A\right)}}-\underset{\mathrm{Term}2}{\underbrace{\sqrt{\frac{2}{c_f}}\left(\frac{c_f}{2St}-\frac{\kappa }{{\kappa }_{\theta }}\right)}},\hfill \end{array}$$

(11)

where the subscript 0 indicates the equivalent smooth-wall case value; and A and $A_{\theta }\left(Pr\right)$ are the previously mentioned velocity and scalar log intercepts for the smooth-wall case. This expression identifies two types of contributions: Term 1 depends only on inner-scale quantities that are invariant to bulk (outer-scale) $Re$, while Term 2 only involves bulk quantities (such as the skin-friction coefficient and Stanton number), which inherently depend on the bulk $Re$; note that the fraction in Term 2 is the inverse of the Reynolds analogy factor. The Reynolds number does not appear explicitly, implying that $\Delta {\Theta }^{+}$ does not scale with $Re$. Furthermore, in the fully rough regime, only $\Delta U^{+}$ and $St$ change for a fixed Prandtl number, with $\Delta U^{+}=log\left(k_s^{+}\right)/\kappa -3.5$ [15,16,17]. Thus, we can then rewrite Equation (11) as:

$$\Delta {\Theta }^{+}=\left[A_{\theta }+\frac{\kappa }{{\kappa }_{\theta }}\left(\sqrt{\frac{2}{c_f}}-A\right)\right]+\left[\frac{\kappa }{{\kappa }_{\theta }}\left(\frac{1}{\kappa }log\left(k_s^{+}\right)-3.5\right)-\sqrt{\frac{c_f}{2}}\frac{1}{St}\right]$$

(12)

where the expression in the first bracket is constant for fixed $Pr$ and given roughness geometry. An asymptotic value $\Delta {\Theta }_{\mathrm{FR}}^{+}$ implies that the terms in the second bracket must also be constant and that $St\sim {log}^{-1}\left(k_s^{+}\right)$. This provides an alternative indicator to verify the behaviour of $\Delta {\Theta }^{+}$ at $k_s^{+}\to \infty$ using a bulk quantity ($St$), which is easier to measure experimentally.

3.3. Scalar-Flux Coefficient

The scalar-flux coefficient, known, in normalised form, as the Stanton number, $St$, is the scalar equivalent of the skin-friction coefficient, $c_f$. The Reynolds analogy factor, $RA$, is the ratio between $St$ and $c_f$. Both $St$ and $RA$ are commonly used to quantifying the effects of roughness, by comparing them to the equivalent smooth-wall case. Hereinafter, we use the subscript 0 to denote the values from the equivalent smooth-wall case (same $R{e}_{\tau }$ and $Pr$), i.e., ${St}_0$, $c_{f0}$, and ${RA}_0$.

Many experimental studies (such as [3,9,12,35,64,65,75] and others) focus on the bulk values of $c_f$ and $R$.

A common objective is to find a predictive empirical correlation for $St$ in various roughness conditions which, in conjunction with known $c_f$ correlations, can be used to calculate $RA$. Numerical simulations, on the other hand, have allowed us to examine the local distributions of their quantities, and to access the flow properties below the roughness crest.

Local values of both $St$ and $c_f$ (and therefore $RA$) can be defined using either the instantaneous or time-averaged scalar flux, $q_w$, shear velocity, $u_{\tau }$, and shear scalar, ${\theta }_{\tau }$ [30,44,45,46]. MacDonald et al. [44] separated the friction- and pressure-drag components of $c_f$ (which is difficult or impossible in experimental studies). The viscous component, $c_{f\nu }$, was then compared with $St$ locally. They are shown over a section of the roughness surface (which was a 3D sinusoid) in Figure 4. For $Pr=0.7$, the behaviours of $St$ and $c_{f\nu }$ are overall similar; the most distinct differences occur in the troughs between roughness elements. Unlike $St$, which is a non-negative quantity, the viscous component, $c_{f\nu }$, is negative in the recirculation regions on the leeward side of the roughness elements. For this sinusoidal roughness, the negative region covers approximately $40\%$ of the plan area, both in the transitionally and fully rough regimes [44]. Thus, $RA$ is more likely to be close to one on the windward side of the roughness elements, or near the crests [44].

A different approach was employed by Zhong et al. [45], who used a compensated viscous-to-diffusive sublayer ratio, ${\mathcal{R}}_{\nu \alpha }=\left({\delta }_{\nu }/{\delta }_{\alpha }\right){Pr}^{1/3}$, to identify where a Reynolds analogy-type behaviour exists. The Reynolds analogy can be assumed to be valid when this ratio is approximately constant (with a value depending on $Pr$). Departures from a constant value indicate weakening of the analogy. It was found that, at the shear-driven windward side of the roughness elements, the Reynolds analogy usually holds locally, as the ratio is almost constant, in line with the findings of [30]. The thicknesses of the viscous and diffusive sublayers on the windward side are comparable to those in the equivalent smooth-wall case. At the crests, a mild departure from a constant is reported, suggesting a partial analogy [45,58]. It is in the detached regions behind the roughness elements where the analogy clearly breaks down, due to the prominence of reversed flow and weak shear. Zhong et al. [45] noted that, while a Reynolds analogy-type mechanism is dominant on the windward side, behind the roughness elements an ensemble of different local behaviours determine the attenuation of the scalar wall flux. The degree of dissimilarity also depends on the relaxation of the detached region, and whether reattachment occurs before the next roughness element (i.e., on the level of sheltering).

This is also confirmed by examining the probability density functions (PDFs) of the instantaneous $c_f$ and $St$, as well as their joint PDF (jPDF), conditioned on flow reversal. Indeed, the conditioning on $u>0$ (attached flow) resulted in similar PDFs for $c_f$ and $St$. In contrast, the condition $u<0$ (reversed-flow region) produced significantly different PDFs for $c_f$ and $St$ [45,46]. Similarly, both the shape and trend of the jPDFs between $c_f$ and $St{Pr}^{2/3}$ (the components of ${\mathcal{R}}_{\nu \alpha }$) on the windward side were very similar to a smooth-wall case; in contrast, the trough regions resulted in very different-looking jPDFs for all the cases examined [45]. These results were consistent over a wide range of roughness heights.

Another qualitative observation of [44] was quantitatively confirmed by the PDFs: higher values of $St$ are associated with the larger positive slopes on the windward side, while lower $St$ values are associated with the negative slopes where the flow is reversed in the troughs, which act as local resistance [40,46]. That means that taller roughness elements are likely to have higher-than-average $St$ on the windward side, and lower-than-average $St$ on the leeward side [30,46,56]. This indicates that describing the rough wall as a “contorted smooth wall” with increased effective wall area (compared with the plan area) is only reasonable in regions where Reynolds analogy-type behaviour exists.

Comparing ${\delta }_{\nu }$ and ${\delta }_{\alpha }$ on the rough surface to the equivalent smooth-wall value is a measure of the local augmentation or attenuation of the scalar wall flux [45]. The deviation of the rough-wall values from the smooth-wall ones can be integrated over the entire wall to estimate the bulk effect of the roughness.

On the windward side of the roughness elements, where the flow is attached and exposed to high shear, a smooth-like boundary layer is observed, thickest at the crest [35,45,64]. Here, both ${\delta }_{\nu }$ and ${\delta }_{\alpha }$ are smaller than the equivalent smooth-wall values, and ${\delta }_{\nu }$ is attenuated more significantly than ${\delta }_{\alpha }$ at $Pr=1.0$ [45]. As $k_s^{+}$ increases, however, both ${\delta }_{\nu }$ and ${\delta }_{\alpha }$ approach their respective smooth-wall values, albeit at a decreasing rate. This trend appears to become insensitive to $Pr$ as $k_s^{+}$ increases; the slowdown of the growth rate is much more drastic for ${\delta }_{\alpha }$ than for ${\delta }_{\nu }$. This implies a local Reynolds analogy-type behaviour at higher $k_s^{+}$, and an increase in $St$ relative to the smooth-wall case on the windward side. Furthermore, for a wavy rough wall, [45] suggested the following correlations for $c_f$ and $St$ at the crests:

$$\begin{array}{cc}\hfill & {c_f}_{\mathrm{crest}}=2\frac{u_{{\delta }_{\nu }}}{U_k}\frac{1}{U_k^{+}{\delta }_{\nu }^{+}}\approx 0.46{\left(k^{+}\right)}^{-0.5}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & {St}_{\mathrm{crest}}=\frac{{\theta }_{{\delta }_{\alpha }}}{{\Theta }_k}\frac{1}{U_k^{+}{\delta }_{\alpha }^{+}}\approx 0.16{\left(k^{+}\right)}^{-0.42}{Pr}^{-0.72}\hfill \end{array}$$

(14)

where $u_{{\delta }_{\nu }}$ and ${\theta }_{{\delta }_{\alpha }}$ are the velocity and scalar at the edge of the viscous or conductive sublayers, and $U_k=\langle \overline{u}{\rangle }_{y=y_{\mathrm{crest}}}$ and ${\Theta }_k=\langle \overline{\theta }{\rangle }_{y=y_{\mathrm{crest}}}$ are the averaged crest velocity and scalar, respectively. In this context, the roughness height, k, is the wavy-wall amplitude, related to the equivalent sand-grain size by $k_s^{+}=2.7k^{+}$.

Numerical simulations also help understand the effect of $Pr$ on $St$ in the different regions. In general, $St\sim {Pr}^{-m}$, where for a smooth wall $m=0.48–0.5$ [3,7,67]. On the windward side and crests, $m\ge 0.5$, while on the leeward side and in recirculation bubbles, $m\le 0.5$. In the wavy-wall example mentioned earlier, $St\sim {Pr}^{-0.38}$ in the troughs, while $St\sim {Pr}^{-0.72}$ at the crests [45]. Rowin et al. [57] further refined this observation, noting that, when the windward slope is partially sheltered, $St\sim {Pr}^{-2/3}$ (corresponding to $RA=1$) on the exposed side of the slope, while $St\sim {Pr}^{-0.45}$ on the sheltered side. They used correlations similar to Equations (13) and (14) to develop a model predicting the scalar transport for rough walls, accounting for the “sheltered” and “exposed” regions of the roughness. The model gives good results for the cases tested (3D sinusoidal roughness with high frontal solidity, and $Pr=0.5–2.0$). Further improvements of the model for lower solidity are, however, desirable.

The local contributions can be integrated over the roughness surface to determine the net effect of roughness on the bulk scalar, which is a more relevant measure for practical engineering applications. The increased contribution of the windward side overcomes the resistance due to the recirculating regions, resulting in a net increase in the scalar transport compared with the smooth-wall case, that is, $St>{St}_0$ [9,35,44,46].

The smooth-wall Stanton number and skin-friction coefficients can be obtained from correlations of the experimental data [3,76]:

$${St}_0=0.021{Re}^{-1/5}{Pr}^{-1/2};c_{f0}=0.073{Re}^{-1/4}.$$

(15)

Both quantities decrease with $Re$ at different rates. In the rough-wall case, $St$ increases with equivalent sand-grain size $k_s^{+}$, reaching a peak in the transitionally rough regime, followed by a monotonic decrease in the fully rough regime [35,39,44,46,48,51]. In contrast, the bulk $c_f$ approaches an asymptotic value in the fully rough regime [17].

Several numerical studies [40,44,45,51,56,57,60,77] examined various types of roughness, systematically changing some topological parameters while maintaining the others fixed. They included variation of the frontal solidity, spacing, density, orientation, slope, etc. They showed that $St$ is both sensitive to the type of roughness and to its topology. The rate of decrease for $St$ seems to be specific to the type of roughness and details of the topology [5,51,52] and, perhaps, also to the value of $Pr$.

This corresponds to the scalar roughness function behaviour discussed in Section 3.2. Some authors (e.g., [5,30]) reported that $St$ decreases at a slower rate than ${St}_0$, indicating that the $St$ increase by roughness becomes larger with $Re$ (although at a decreasing rate), while others [45,58] reported the opposite. This suggests that, at high enough $k_s^{+}$, $St$ may be smaller than in the corresponding smooth-wall case at the same values of $Re$ and $Pr$. It would be worthwhile to explore this issue in a systematic manner.

The frontal solidity, ${\lambda }_f$, is a key parameter associated with impingement and recirculation. The effective slope, ES, plays the same role as ${\lambda }_f$, as the two are directly related by $ES=2{\lambda }_f$ [78]. Therefore, low values of ${\lambda }_f$ correspond to sparsely packed shallow roughness, while high values of ${\lambda }_f$ correspond to densely packed steep roughness [16]. Depending on the roughness geometry, a peak of $St$ can be seen in $0.1<{\lambda }_f<0.3$, usually corresponding to maximum $c_f$ [39,45,51].

Forooghi et al. [51] used distribution of truncated cones, and varied independently the height, effective diameter (i.e., the diameter of a cylinder with the same frontal area), and the spacing between the roughness elements. They found that as the spacing decreases (while maintaining the other parameters, i.e., the roughness becomes more “dense”) $St$ reaches a maximum and then decreases. It was suggested that this behaviour occurs because more of the upwind portions becomes sheltered (partially or fully), which limits the amount of unmixed fluid that could impinge on the surface [45,51]. In the limit of zero spacing (or infinite density), a smooth-wall equivalent is expected with $St\to {St}_0$, consistent with this trend. In contrast, increasing in the slope by scaling the geometry while maintaining the elevation resulted in a monotonic increase in $St$.

As discussed, $c_f$ and $St$ behave differently at increasing $Re$ and for different roughness geometries, and, therefore, $RA$ depends on the morphological parameters for both regular and irregular roughness [5,44,51]; this is a result of the large differences between momentum and scalar in the recirculation region downstream of roughness elements.

An expression for $RA$ can be directly derived from a semi-empirical expression for the scalar roughness function proposed by Kays and Crawford [3]. Kuwata [60] modified this expression based on the phenomenological behaviour of $\Delta {\Theta }^{+}$ discussed in Section 3.2 and on their DNS data; their expression predicted all the results within a $4.2\%$ error. Its universality, however, was not established and the authors suggest the possibility of adding a coefficient dependent on roughness type.

Figure 5 compares $c_f/c_{f0}$ and $St/{St}_0$. In most cases, the skin friction increases faster than the Stanton number, implying $RA<1$; in Ref. [41], however, the opposite behaviour was observed. It should be pointed out that in this paper high-aspect-ratio streamwise ribs were used, whose behaviour differs significantly from that of distributed roughness because of the rollers generated by the Kelvin–Helmoltz instability, which suppresses sweeps and ejections, resulting in lower Reynolds stress [41,77]. A similar effect can be achieved using blowing or suction at the wall [77]. It is important to note that $RA>1$ can only be attained for a limited range of $Re$; once the rib spacing, $s^{+}$, becomes too large the increased drag overcomes the augmentation of scalar transfer, leading to $RA<1$; this is true for various types of longitudinal ribs and riblets [41,77] where the range of $s^{+}$ at which $RA>1$ depends on the geometry.

It is also possible to examine how $RA$ compares to the equivalent smooth-wall value ${RA}_0$ as function of $k_s^{+}$. The rough-to-smooth ratio of the Reynolds analogy factor, $\eta$,

$$\eta =\frac{RA}{{RA}_0}=\frac{St}{{St}_0}\frac{c_{f0}}{c_f}.$$

(16)

is an efficiency parameter commonly known as the “aero-thermal efficiency factor”. We will use the same terminology here, understanding it is applicable to the various scalar quantities. In most cases, $\eta <1$ (i.e., $\Delta U^{+}>\Delta {\Theta }^{+}$); riblets or other high-aspect-ratio streamwise longitudinal ribs have a different behaviour [41,77]. Using the discrete element approach (in which the flow is spatially averaged over the roughness elements, modifying the mean flow equations to account for blockage effects, as well as the drag and scalar flux), Aupoix [80] suggested:

$$\eta =\frac{1-\Delta U^{+}\sqrt{c_{f0}/2}}{1-{RA}_0\Delta {\Theta }^{+}\sqrt{c_{f0}/2}}$$

(17)

which was found to give good predictions for various types of roughness [46,51].

Forooghi et al. [51] observed that $\eta$ is a convex function of the roughness density and the frontal solidity, ${\lambda }_f$; predicting the ${\lambda }_f$ value of the minimum for different types of roughness is still an open question. It was also shown [5,51,56] that $\eta$ decreases both with an increase in $k_s^{+}$ or in the effective slope, $ES$. This decrease eventually saturates, Figure 5b. Forooghi et al. [51] derived an expression for $\eta$ based on a systematic DNS study, using various types of roughness for $Pr=0.7$:

$$\eta =0.55+0.45e^{-k_s^{+}/130}$$

(18)

similar to the proposal by Bons [72], but using $k_s^{+}$ instead of the outer-scaled $k_s$. This is meant to be a simplified single-parameter correlation rather than a general expression that accounts for all possible effects; nevertheless, although only requiring the knowledge of $k_s^{+}$, Equation (18) predicts all the roughness morphologies considered by Forooghi et al. [51] (with a single exception) within a $10\%$ error, as well as the results of Jelly et al. [5], Nagura et al. [56], Aupoix [80]. Interestingly enough, Equation (18) also predicts well the data of Peeters and Sandham [46] (grit-blasted surfaces) and Kuwata and Kawaguchi [33] (irregular roughness) for $Pr=1.0$. One exception is the additive-manufacturing (AM) roughness of [54], which follows the same trend but falls below all the other results. This may indicate some significant difference in the properties of AM compared with other types of roughness.

3.4. Dispersive Flux

The roughness geometry generates a stationary inhomogeneity in the velocity and scalar fields, referred to as a wake or dispersive field. The wake velocity, $\tilde{u}$, and scalar, $\tilde{\theta }$, fields are compared in Figure 6 for $Pr=1$, showing overall similarity throughout the RSL. A difference can be observed in the troughs between roughness elements, where $\tilde{u}$ and $\tilde{\theta }$ tend to have opposite signs [31,58], the result of the accumulation of scalar near the wall due to the lack of a pressure-like mechanism, discussed in Section 3.1. For $Pr\ge 1.0$, it is expected that the $\tilde{\theta }>0$ regions would be more extensive than the $\tilde{u}>0$ regions.

The normal dispersive stresses, ${\tilde{u}}_k{\tilde{u}}_k$, tend to be largest at the crest of the elements. Their streamwise component, $\tilde{u}\tilde{u}$, is the largest in magnitude, and has been found to be insensitive to $k_s^{+}$ [58,60]. A secondary peak of high ${\tilde{u}}_k{\tilde{u}}_k$ can be found around the troughs, usually at the edge of the recirculation bubble. These locations generally correspond to the points where the stochastic normal Reynolds stresses, $\langle \overline{u_k^{\prime }{u}^{\prime }k}\rangle$, are smaller [32,58].

Contours of the shear components that contribute to the momentum and scalar transport, $\langle \tilde{u}\tilde{v}\rangle$ and $\langle \tilde{\theta }\tilde{v}\rangle$, from Zhang et al. [58], are shown in Figure 7. The wall-normal wake velocity, $\tilde{v}$, is particularly large in magnitude around the crests and close to the trough, and is directed towards the boundary [58,60,81]. The difference between the dispersive stresses and fluxes is due to the differences in size and magnitude between $\tilde{\theta }$ and $\tilde{u}$, as well as their correlation with locations of intense $\tilde{v}$ (note that in a channel $\langle \overline{v}\rangle =\langle \overline{w}\rangle =0$, so that $\tilde{v}=\overline{v}$ and $\tilde{w}=\overline{w}$). Both $\tilde{u}\tilde{v}$ and $\tilde{\theta }\tilde{v}$ are large and negative around the crest but positive on the leeward slope, where $\tilde{\theta }\tilde{v}$ is more dominant. In the secondary peak near the trough, $\tilde{u}\tilde{v}$ is positive, while $\tilde{\theta }\tilde{v}$ is negative, and peaks slightly further downstream, close to the reattachment point [58]. This is due to the positive region of $\tilde{\theta }$ extending towards the recirculation region and into regions of negative $\tilde{v}$, as discussed before.

Zhang et al. [58] noted that the maximum value of the various dispersive quantities, in outer units, increased linearly as the slope became steeper. The dispersive normal stresses, shear stress, and vertical scalar dispersive flux all showed similar trends, while the streamwise dispersive scalar flux showed a more moderate increase.

From an average perspective, the streamwise component is the the dominant one for both dispersive fluxes and stresses [60]. The flux component, $\langle \tilde{\theta }\tilde{v}{\rangle }_s$, was reported to vary between $25\%$ [32] and $50\%$ [82] of the total scalar flux inside the RSL. At sufficiently large $k_s^{+}$, the mean dispersive stresses and fluxes become independent of $k_s^{+}$ and $Re$ [32,40,52,60]; a similar trend was found for $\langle \tilde{\theta }{\tilde{u}}_j{\rangle }_s$ with increasing $Pr$ [32]. Both the dispersive stresses and fluxes remain dependent on the roughness topological details.

Overall, the dispersive fluxes are not well studied and may be a very useful avenue of numerical exploration, since they are very difficult to measure experimentally. They seem to be less sensitive to $k_s^{+}$ [60] and even reach saturation [32,58,60]. They are also important for the production of turbulent fluxes [31,32] and in the context of modelling, especially if the double-averaged NS equations (DANS) are used [80,82,83].

3.5. Form-Induced Scalar Flux

A key difference between the transport of momentum and scalar, in the context of roughness, is the mechanism by which the flux across the solid boundary takes place. The transfer of momentum becomes progressively dominated by the pressure drag, $F_p$, as $k_s^{+}$ increases, while the viscous drag, $F_{\nu }$, decreases. The scalar wall transfer, however, is only due to diffusive mechanisms [3,9,35], and is denoted here as“form-induced (FI) transfer”, $F_{\alpha }$ (sometimes referred to as “wall-interaction transfer term” [52,56,60]). The spatial distributions of $F_{\nu }$, $F_p$, and $F_{\alpha }$ are difficult to obtain experimentally but can be easily calculated numerically [23,62,84].

The FI scalar transfer, $F_{\alpha }$, is the spatial integral of the local FI flux $f_{\alpha }$, found on the RHS of Equation (9):

$$f_{\alpha }=\langle \frac{\partial }{\partial x_j}\left(\alpha \frac{\partial \tilde{\theta }}{\partial x_j}\right){\rangle }_s=-\frac{1}{A_0}\underset{\partial S_f}{\oint }\alpha \frac{\partial \tilde{\theta }}{\partial x_k}\left(n_{k}dl\right)$$

and the viscous drag, $F_{\nu }$, and pressure drag, $F_p$, are the integrals of

$$\begin{array}{cc}\hfill & f_{\nu }=\langle \frac{\partial }{\partial x_j}\left(\nu \frac{\partial {\tilde{u}}_i}{\partial x_j}\right){\rangle }_s=-\frac{1}{A_0}\underset{\partial S_f}{\oint }\nu \frac{\partial {\tilde{u}}_i}{\partial x_k}\left(n_{k}dl\right)\hfill \\ \hfill & f_p=-\langle \frac{\partial \tilde{p}}{\partial x_i}{\rangle }_s=-\frac{1}{A_0}\underset{\partial S_f}{\oint }\tilde{p}\left(n_{k}dl\right)\hfill \end{array}$$

respectively, found on the RHS of Equation (8). The integration domain, $\partial S_f$, is the part of the roughness surface that intersects the averaging domain of size $A_0$ (i.e., the perimeter of the roughness surface intersecting a horizontal plane of size $A_0$ for a channel) and $n_k$ is the k-th component of the outward-pointing unit vector normal to the roughness surface.

As the bottom of the RSL is approached, $F_p$ and $F_{\nu }$ dominate the mean momentum transfer (with $F_p$ being more significant at increasing $k_s^{+}$), while $F_{\alpha }$ dominates the mean scalar transfer [31,32,33].

The total FI momentum flux, $f_{\nu }^{+}+f_p^{+}$, is generally larger than the scalar FI one, $f_{\alpha }^{+}$, [32,33], which results in the total drag ($F_{\nu }^{+}+F_p^{+}$) being larger than the scalar wall transfer ($F_{\alpha }^{+}$). Kuwata and Kawaguchi [33] noted that, in the region closest to the bottom of the roughness, the opposite ($f_{\nu }^{+}+f_p^{+}<f_{\alpha }^{+}$) can occur locally due to the negative local shear stress.

Kuwata [52] showed that the total momentum FI, $f_{\nu }^{+}+f_p^{+}$, and scalar FI flux, $f_{\alpha }^{+}$, have similar behaviours when scaled in inner units and plotted against non-dimensionalised wall distance. Furthermore, the increase in the effective slope, $ES$, shifts the main contribution of $f_{\alpha }^{+}$ away from the base of the roughness. A possible explanation for this is the increase in surface area away from the base, as well as the presence of the well-mixed fluid between the roughness elements (see Section 3.1), referred to as an “insulation layer”, near the base of the roughness. Consistent with the discussion of the thinning of the diffusive layer, which tends to conform to the roughness shape when either $Re$ or $Pr$ increase, $F_{\alpha }^{+}$ becomes proportional to the solid fraction [31,32], as shown in Figure 8.

3.6. Scalar Fluctuations and Turbulent Structures

Roughness has a profound impact on both velocity and scalar fluctuations, and numerical simulations are able to access them even below the roughness crest. Several studies compared stochastic stresses and scalar correlations inside the RSL, contributing to a better understanding of the dynamics of momentum and scalar transport.

Roughness breaks up both velocity and scalar streaks, resulting in more isotropic turbulence, evidenced by both stochastic stresses, $\langle \overline{u_i^{\prime }{u}_j^{\prime }}{\rangle }_s$, and scalar correlations, $\langle \overline{u_i^{\prime }{\theta }^{\prime }}{\rangle }_s$. Numerical studies show that both the magnitude and peak location of $\langle \overline{{\theta }^{\prime }{\theta }^{\prime }}{\rangle }_s$ are affected: for $Pr\ge 1$, the peak is lower in the rough-wall case [32,46,52], while for lower values of $Pr$ it increases [32]. Roughness attenuates the mean streamwise velocity fluctuations (normalised in wall units), $\langle \overline{u^{\prime }{u}^{\prime }}{\rangle }_s$, more than the scalar ones, $\langle \overline{{\theta }^{\prime }{\theta }^{\prime }}{\rangle }_s$, for $Pr$ near unity [29,31,46,52]. It could be expected that the mean wall-normal stochastic flux, $\langle \overline{{\theta }^{\prime }{v}^{\prime }}{\rangle }_s$, in wall units, is less affected by roughness than $\langle \overline{u^{\prime }{v}^{\prime }}{\rangle }_s$. For $Pr=1.0$, however, they are similar at low-to-moderate $k_s^{+}$. However, the differences between the two inside the RSL increase with roughness size [34,46,60]; outside the RSL, $\langle \overline{{\theta }^{\prime }{v}^{\prime }}{\rangle }_s$ obey Townsend’s outer-layer similarity [25,32].

Various studies [29,31,32,39,46,52,60] examined the time-averaged field, $\overline{{\theta }^{\prime }{v}^{\prime }}$. As mentioned in Section 3.1, the main differences between the scalar and velocity fields occur in recirculation regions behind roughness elements and in the sheltered portion of the windward slopes. This is also evident for the turbulent scalar fluxes, Figure 9, where the time-averaged wall-normal Reynolds stress, $-{\overline{u’v’}}^{+}$, and scalar flux, $-{\overline{{\theta }^{\prime }{v}^{\prime }}}^{+}$, are shown in a vertical plane, for $k_s^{+}\approx 200$ and $Pr=1.0$. For the most part, the two fields are visually similar, which is expected due to the fact that both are driven by the same turbulent eddies and associated with $v^{\prime }$. There are noticeable differences in the recirculation regions downstream of larger roughness elements, where the scalar flux is more intense than the stress (both in positive and negative values). Quasi-streamwise elongated vortices develop in narrow gaps in the roughness, which strengthen both the turbulent stress and flux close to the rough surface [60]. For $Pr>1$, this enhances the scalar turbulent flux more than the stress.

Quadrant analysis (QA) [85] can be used to examine the instantaneous behaviour through the probabilities of turbulent events, as well as their relative contribution to the Reynolds stresses and scalar fluxes. The quadrants are defined as Q1: $u^{\prime }>0,v^{\prime }>0$; Q2: $u^{\prime }<0,v^{\prime }>0$; Q3: $u^{\prime }<0v^{\prime }<0$; and Q4: $u^{\prime }>0,v^{\prime }<0$. Q2 and Q4 events, in flat-plate boundary layers and plane channels, correspond to ejections (slow fluid moving away from the wall) and sweeps (fast fluid moving towards the wall), respectively. They are the main contributors to the Reynolds shear stress, $\overline{u^{\prime }{v}^{\prime }}$. Replacing $u^{\prime }$ with ${\theta }^{\prime }$ provides an additional way to examine the relationship between the flow topology and the scalar fluxes. Note that in the classical quadrant analysis the fluctuations normal to the wall are used; here, instead, $u^{\prime }$ and $v^{\prime }$ are fluctuations parallel and normal to the surface at the bottom of the roughness, and not to the roughness surface. Thus, Q2 events represent fluid going away from the roughness layer towards the outer flow, and, conversely, Q4 events represent outer-layer fluid advected into the troughs.

Figure 10 shows the distribution of events, sorted by quadrants, in the near-wall region. Above the roughness crest, the smooth-wall distribution is found both in the $u^{\prime }-v^{\prime }$ and ${\theta }^{\prime }-v^{\prime }$ planes, with Q2 and Q4 events dominating. Inside the RSL, on the other hand, the probability of the four quadrants becomes roughly equal for the stress, tending to approximately $25\%$ as the base of the roughness is approached. For the scalar flux (at $Pr=1.0$), however, the probability of Q4 motions, where ${\theta }^{\prime }>0$ & $v^{\prime }<0$, is meaningfully diminished, while the probability of Q3 events, where ${\theta }^{\prime }<0$ & $v^{\prime }>0$, is significantly enhanced (Figure 10b) [39,46]. The difference in event probability, however, does not translate into a contribution to the mean flux and stress [46,82]. Examining the relative contributions to the mean (double-averaged) scalar flux, ${\langle \overline{{\theta }^{\prime }{v}^{\prime }}\rangle }_i/\langle \overline{{\theta }^{\prime }{v}^{\prime }}\rangle$ shows that the relative contribution of Q4 events remains larger than that of Q2 events (Figure 10d, similar to the counterpart quadrants ${\langle \overline{u^{\prime }{v}^{\prime }}\rangle }_i$. This implies that fewer Q4 motions occur (e.g., motions of hot fluid towards the troughs) but they are generally higher in intensity compared with the smooth-wall case [29,39,46]. These results are consistent with the observations of Doosttalab et al. [29] that roughness enhances the transport of both ${\overline{v^{\prime }{u}^{\prime }{v}^{\prime }}}^{+}$ and ${\overline{v^{\prime }{\theta }^{\prime }{v}^{\prime }}}^{+}$ towards the troughs, as well as the increasing the skewness of $v^{\prime }$ and ${\theta }^{\prime }$, mentioned earlier.

Despite the conclusion that the relative contributions of Q1–Q4 to the mean are similar between momentum and scalar, there is meaningful dissimilarity in the type of motions experienced by the momentum and scalar, implying a breakdown of analogy locally between the two in the near-wall region [46], which corresponds to $y/k<1/3$ in the case presented in Figure 10.

3.7. Scalar-Fluctuation Budget

Roughness was shown to affect the scalar fluctuations and the TKE differently [31,32,46,60]. These differences are reflected in the budgets for the double-averaged scalar fluctuation variance ${\mathcal{K}}_{\theta }\equiv \langle \overline{{\theta }^{\prime }{\theta }^{\prime }}{\rangle }_s/2$ and TKE. Figure 11 shows the various budget terms in the fully rough regime for $Pr=1$. A key reason for the difference between velocity and scalar fluctuations lies in the production of the TKE and scalar variance near the roughness crests. In addition to the shear production, ${\mathcal{P}}^s$, which also exists in a smooth-wall case, a second “bypass” production mechanism exists, which generates velocity and scalar fluctuations at a scale comparable to the roughness height [23,31,62,63]. This form-induced (FI) production, ${\mathcal{P}}^{fi}$, stems from the gradients of the wake fields of the velocity and scalar:

$$\begin{array}{cc}\hfill {\mathcal{P}}_{ij}^{fi}& =-\langle \overline{u_i^{\prime }{u}_m^{\prime }}\frac{\partial {\tilde{u}}_j}{\partial x_m}+\overline{u_j^{\prime }{u}_m^{\prime }}\frac{\partial {\tilde{u}}_i}{\partial x_m}{\rangle }_s\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathcal{P}}_{\theta }^{fi}& =-2\langle \overline{\theta u_m^{\prime }}\frac{\partial \tilde{\theta }}{\partial x_m}{\rangle }_s.\hfill \end{array}$$

(20)

The total production is $\mathcal{P}={\mathcal{P}}^s+{\mathcal{P}}^{fi}$.

If the wall is rough, the total production (in wall units) of both velocity and scalar fluctuations is reduced, compared with the smooth-wall case; this reduction depends on $k_s^{+}$ [29,31,32]. The reduction in the velocity fluctuation production is more significant than that of the scalar ones: for transitional roughness, Doosttalab et al. [29] reported a reduction of ∼5% in the production of ${\mathcal{K}}_{\theta }$ compared with ∼17% for TKE production, while for a fully rough regime Hantsis and Piomelli [31] observed a reduction closer to $50\%$ and $70\%$ for ${\mathcal{K}}_{\theta }$ and TKE, respectively.

While the scalar shear production, ${\mathcal{P}}_{\theta }^s$, increases with $Pr$, it quickly reaches saturation near $Pr\simeq 1$, where the curves for increasing $Pr$ collapse [32]; this is consistent with the observation that the bulk contribution of the shear production is in the upper portion of the roughness due to strong shear layers emanating from the crests, where the motions are generally insensitive to $Pr$. The shear production term for both TKE and ${\mathcal{K}}_{\theta }$ decreases with $k_s^{+}$; the two remain similar for $Pr\ge 1$ across the range of $k_s^{+}$ [31,32]. It is not known if the shear production continues to decrease with $k_s^{+}$ or eventually reaches saturation.

In contrast, the scalar FI production, ${\mathcal{P}}_{\theta }^{fi}$, becomes increasingly dominant inside the RSL as $k_s^{+}$ increases, compared with the TKE counterpart, ${\mathcal{P}}_{kk}^{fi}$; the magnitude of ${\mathcal{P}}_{\theta }^{fi}$ increased with $Pr$ and does not seem to saturate [31,32]. Furthermore, for a given $Pr$, the ${\mathcal{P}}_{\theta }^{fi}$ curves collapse as $k_s^{+}$ increases [32], as seen in Figure 11b, which corresponds to the geometry-conforming behaviour of the scalar field discussed in Section 3.1. Thus, it is the FI portion that is responsible for the increased scalar fluctuations’ production. The bulk contribution of ${\mathcal{P}}_{\theta }^{fi}$ occurs in the lower portion of the roughness, and is small around the crest. As $k_s^{+}$ increases, the scalar variance budget becomes dominated by the FI production, ${\mathcal{P}}_{\theta }^{fi}$, diffusion, and dissipation; both diffusion and FI production tend to saturate with increasing $k_s^{+}$ [32]. A saturation in the shear production then implies saturation in the dissipation as well, resulting in a budget that is insensitive to $k_s^{+}$. A more systematic study is required to ascertain these trends at higher $k_s^{+}$, for a wider range of $Pr$ and various types of roughness. Asymptotic trends such as these are useful, in particular in the context of modelling.

3.8. Turbulent Prandtl Number

The turbulent Prandtl number is the ratio between the eddy viscosity and the scalar eddy diffusivity, ${Pr}_T={\nu }_T/{\alpha }_T$. It is used to model scalar transport [3] in analogy to momentum transport. A direct consequence of the Reynolds analogy is that the turbulent Prandtl number, ${Pr}_T$, needs to be constant and close to unity for $Pr\approx 1.0$. Although many experimental studies examined ${Pr}_T$, to properly calculate it, an accurate and simultaneous measurement of at least four quantities (stochastic stress and flux, mean velocity, and scalar gradients) is required, at multiple locations—which is extremely difficult experimentally [3,86]. This difficulty resulted in a relatively large scatter between experimental results [3,29]; numerical simulations can compute these quantities more accurately.

While ${Pr}_t$ is generally obtained, in thin shear-layer flows, from the only non-zero Reynolds shear stress (or scalar flux) and velocity (or scalar) gradient, a more general form is obtained through a contraction with the strain rate or scalar gradient [87]:

$${Pr}_T=\frac{{\nu }_T}{{\alpha }_T}=\frac{\langle \overline{u_i^{\prime }{u}_j^{\prime }}\rangle \langle {\overline{S}}_{ij}\rangle }{2\langle {\overline{S}}_{lm}\rangle \langle {\overline{S}}_{lm}\rangle }\times \frac{\langle \frac{\partial \overline{\theta }}{\partial x_k}\rangle \langle \frac{\partial \overline{\theta }}{\partial x_k}\rangle }{\langle \overline{u_n^{\prime }{\theta }^{\prime }}\rangle \langle \frac{\partial \overline{\theta }}{\partial x_n}\rangle }.$$

(21)

${Pr}_T$ decreases with wall-normal distance, starting from a value near one at the wall (e.g., ${Pr}_T=1.1$ for $Pr=0.7$ [88]) to 0.6–0.8 away from the wall [89]. It tends to a roughly constant value in the logarithmic layer, with reported values falling within the range of 0.8–1.0, where the most common values are around 0.85–0.9 (see [3] for comprehensive discussion). In the context of modelling, however, it is often taken as a constant, usually the logarithmic layer value.

All the terms comprising ${Pr}_T$ are affected by roughness within the RSL, while away from the wall, corresponding with Townsend’s outer-layer similarity hypothesis [25], the various quantities collapse on the smooth-wall curves; the differences in ${Pr}_T$ between smooth- and rough-wall cases are limited to the RSL, where numerical simulations can provide more information than experimental measurements.

In the presence of roughness, ${Pr}_T$ is ill-behaved close to the bottom of the roughness, as the denominator becomes small [29,34,46,52]; the magnitude of these changes increases with $k_s^{+}$ and seems to depend on the roughness geometry [52,60]. For $Pr=1$, the rough-wall ${Pr}_T$ can be much greater than one; for example, [46] reported peak values of 5. A higher peak indicates a decrease in the effective scalar diffusivity, resulting in a decreased scalar flux/transfer (i.e., in thermal resistance). The erratic behaviour near the bottom can be related to the difference in the behaviour of the stochastic stresses and fluxes in the recirculation regions [33,34,46], as is demonstrated in Figure 9. In fact, it may be argued that in the vicinity of the wall, where the Reynolds stress and turbulent scalar flux are nearly damped out, ${Pr}_T$ does not have physically meaningful values [33,46].

As the roughness crest is approached, a more predictable behaviour is observed, with ${Pr}_T$ peaking and then decreasing towards the smooth-wall value [29,34,46]. The peak magnitude increases with $k_s^{+}$ and seems to depend on the details of the geometry as well [46,52,60]. The location of the peak may depend on the roughness geometry [29,34,46,52].

While ${Pr}_T$ is a key element for modelling, it may be too sensitive to the details of the roughness to be used effectively. The existence of additional dispersive stresses and fluxes (discussed in Section 3.4) in the mean transport equations, which are neglected in the classical definition Equation (21), motivated the definition of an “effective” Prandtl number [34,52]:

$${Pr}_{T,\mathrm{eff}}=\frac{{\nu }_T+{\tilde{\nu }}_T}{{\alpha }_T+{\tilde{\alpha }}_T}=\frac{\left(\langle \overline{u_i^{\prime }{u}_j^{\prime }}\rangle +\langle \widetilde{u_i}\widetilde{u_j}\rangle \right)\langle {\overline{S}}_{ij}\rangle }{2\langle {\overline{S}}_{lm}\rangle \langle {\overline{S}}_{lm}\rangle }\times \frac{\langle \frac{\partial \overline{\theta }}{\partial x_k}\rangle \langle \frac{\partial \overline{\theta }}{\partial x_k}\rangle }{\left(\langle \overline{u_i^{\prime }{\theta }^{\prime }}\rangle +\langle \widetilde{u_i}\tilde{\theta }\rangle \right)\langle \frac{\partial \overline{\theta }}{\partial x_i}\rangle }.$$

(22)

The expectation is that the inclusion of the dispersive stresses and fluxes will make ${Pr}_{T,\mathrm{eff}}$ better behaved, even where the stochastic stresses and fluxes go to zero near the wall. Furthermore, in the context of modelling, the dispersive stresses and fluxes are required for the closure relations [83].

Hantsis and Piomelli [34] compared ${Pr}_{T,\mathrm{eff}}$ and ${Pr}_T$, and found them to be generally similar in the $RSL$, collapsing above the roughness crest, where the dispersive stress and fluxes vanish. ${Pr}_{T,\mathrm{eff}}$ is also ill-behaved near the bottom of the rough surface. For $Pr\simeq 1$, ref. [34] found that ${Pr}_{T,\mathrm{eff}}$ is larger in magnitude than ${Pr}_T$. A peak, larger than that observed in ${Pr}_T$, was also observed for ${Pr}_{T,\mathrm{eff}}$, but occurred closer to the roughness crest. Due to the scarcity of studies considering ${Pr}_{T,\mathrm{eff}}$, and comparing ${Pr}_{T,\mathrm{eff}}$ to ${Pr}_T$ under comparable conditions, it is unclear if this behaviour is particular to this setup or has more general validity.

The dependence on the details of the roughness is illustrated in Figure 12, showing ${Pr}_{T,\mathrm{eff}}$ for different types of roughness and topologies, for $Pr$ of unity. Four representative roughness surfaces from [52] show the behaviour at different skewness ($S_k$) and effective-slope $\left(ES\right)$ values for an irregular roughness, while the sand-grain roughness of [34] shows the change of roughness type. A common pattern of positive and negative peaks can be noticed as the wall is approached. The magnitudes and locations of the peak values (where the jumps occur), however, differ significantly between the various cases.

Neither ${Pr}_T$ nor ${Pr}_{T,\mathrm{eff}}$ seem to be effective modelling tools in generalised rough-wall flows, and an alternative might be required for modelling. In the context of smooth-wall flows, it was found that quantities such as the time-scale ratio, $\mathcal{R}$ (between the scalar variance and the turbulent kinetic energy) are more useful in modelling, and can replace the need for ${Pr}_T$ Abe and Antonia [90]. While a similar observation was made for a rough-wall case [34], only one case was considered, and further study would help in examining the viability of this alternative.

4. Open Issues

4.1. Boundary Conditions

The mixed-boundary condition proposed by Kasagi et al. [61] is commonly applied for the uniform scalar isoflux condition, when periodic domains are used. Kasagi et al. [61] derived this condition for a smooth-wall case, where, for a limited range of $Pr$ around unity, the fluctuations at the wall were negligible, and the wall could be assumed to be locally isothermal [91]. As was discussed in Section 3.6, roughness tends to reduce scalar fluctuations near the wall, which suggests the assumption is still valid for a rough-wall case; indeed, many studies implemented this approach [5,30,31,32,34,44,45,46,51,54,55,57,59,60,92].

Two key issues arise from this setup: (1) the use of this condition is limited to a narrow range of $Pr$, and (2) a proper validation of the locally isothermal assumption was not conducted in the context of rough-wall flows (see [93] for a smooth-wall case comparison). Thus, a boundary condition that allows for fluctuations at the wall is desirable.

Another issue arises from the use of the mean wall-scalar, $T_w$, in the case of a wall flux condition: for a smooth wall, $T_w$ has a constant streamwise gradient, equal to that of the bulk scalar, $T_b$. For the rough-wall case, the spatial inhomogeneity of the roughness implies that the streamwise gradient of $T_w$ depends on the surface area of the roughness [52], and a local spanwise gradient may also exist. Kuwata [52] suggested a correction by considering the local time-averaged scalar budget at the roughness surface and relating $T_w$ to $T_b$; this yields

$$\frac{d{T}_b}{dx}=\frac{q_w}{\rho c_p\dot{V}}{\mathcal{l}}_{yz}$$

(23)

where $\dot{V}$ is the flow rate and ${\mathcal{l}}_{yz}$ is the local roughness length scale. The source term Q then becomes

$$Q=\frac{q_w}{\rho c_p\dot{V}}\left(u{\mathcal{l}}_{yz}-\alpha \frac{d{\mathcal{l}}_{yz}}{dx}\right)$$

(24)

For a sufficiently large roughness surface of a homogeneous and isotropic nature, ${\mathcal{l}}_{yz}$ can be approximated as ${\mathcal{l}}_{yz}\approx A_r/L_x$, where $A_r$ is the total roughness wetted area and $L_x$ is the streamwise domain length (assuming roughness covers the entire boundary). The source term Q then can be approximated as

$$Q\approx \frac{A_r}{A}\frac{q_w}{\rho c_p{\delta }_e}\frac{u}{U_b}$$

(25)

where ${\delta }_e$ is the effective half-channel height (${\delta }_e=\delta -d$). While this formulation better accounts for the local wall flux, there has not been a comparative study examining the significance of this correction.

4.2. Limited Coverage of the Parameter Space

The study of flow over rough surfaces is a complex problem with a large parameter space [16]. As discussed in earlier sections, the parameter space associated with scalar transport is larger than that of momentum transport, at the very least by the addition of a Prandtl number (or equivalent). The need to resolve all (for DNS) or most (for LES) of the scales of motion becomes prohibitive as the Reynolds number increases. For $Pr>1$, the resolution requirements for the scalar are stricter than those for momentum, making the calculations even more expensive.

The presence of roughness increases the parameter space further. Jiménez [28] recommends $\delta /k>40$ and $k^{+}>100$ to achieve the fully rough regime with negligible blockage effects. Most present studies relax one or both of these requirements because of the costs associated with them.

Some studies [32,45] indicated the possibility that, with proper scaling at sufficiently high $k_s^{+}$, various scalar statistics may become independent of $Pr$. A much wider range of $Pr$ is required to validate this conjecture, which could have repercussions on modelling of these problems.

The roughness geometry is arguably the parameter space that is most difficult to cover. Beyond covering an immense spectrum of roughness types (i.e., regular, irregular, k-type, d-type, etc.), it also includes different morphological characteristics of those roughness types. To examine if any of the findings and trends have some level of universality, a wider portion of this parameter space needs to be covered systematically. In most cases, a single roughness type is considered, with varying the roughness height, usually by increasing $Re$, and usually for a single value of $Pr$.

Since only near-wall phenomena are of interest, the minimal channel approach, proposed by Jiménez and Moin [94], can be used to decrease the computational expense of channel-flow calculations and facilitate sweeping across the parameter space. Chung et al. [95] demonstrated the applicability of the approach to rough-wall flows; MacDonald et al. [44] extended it to also include scalar transport. Limiting the width or length of the domain (see [95,96] for more details) removes the larger motions while retaining many of the dynamically relevant scales within the RSL [97], in particular at higher $Re$. This approach was used successfully to sweep across the topological space of 3D sinusoidal roughness and $k_s^{+}$ [45,57,84,95,96,97,98].

Relatively few studies ([51,55,80]) systematically cover a range of roughness types. Thus, it is difficult to form conclusions on specific quantities across the parameter space by mapping out results from various studies, and a more systematic approach is desirable. Forooghi et al. [51] performed a systematic study of several roughness types, and several investigations consider variations of a single roughness type [30,44,45,57,92]. The importance of the geometrical detail, however, remains an open question.

Because of the wide parameter space and its associated computational cost, only a limited number of the studies mentioned here reached $k_s^{+}$ above 500 (see Table 1 for a non-exhaustive list), while some remained in the transitional regime. Most studies tend to be in the range $k_s^{+}=200–400$. While the fully rough regime starts at $k_s^{+}\approx 100$ [28], the onset of the scalar “ultimate regime” seems to depend on both the roughness geometry and $Pr$. For $Pr=0.7$, MacDonald et al. [44] suggested the onset occurring beyond $k_s^{+}=250$. Hantsis [32] observed that as $Pr$ the onset of the ultimate regime occurs at higher $k_s^{+}$ (see Figure 3a).

As such, much of the discussion around the ultimate regime (e.g., Section 3.2) and the behaviour of various quantities such as $\Delta {\Theta }^{+}$ requires a wider range of $k_s^{+}$ to determine if the trends and findings extend in their applicability beyond the limited scope of $k_s^{+}$ examined by these studies, leaving much room for speculation. This is particularly true for higher Prandtl numbers.

The increase in Prandtl numbers also results in improved scalar transfer for both smooth- and rough-wall cases, which is important for various applications. For example, in the context of nuclear reactor cooling, Gowen and Smith [65] considered a Prandtl number range of $0.7–14.3$. Most wall-resolved numerical studies remained in the vicinity of $Pr=0.7–1.0$, with $Pr=0.7$ being of particular interest. In some cases it was due to the limitation computational cost, while in other cases it was either due to the relevance to practical applications (e.g., $Pr\approx 0.7$ for air) or existing experimental data to compare against, or making the most direct comparison between the velocity and scalar (with $Pr=1.0$ providing this case). Furthermore, most studies only considered a single $Pr$, while the few studies that examined a range of $Pr$ were usually limited to a small range around unity, say $Pr=0.5–2.0$ [31,32,45,48]. While the latter can be due to the increased computational costs associated with higher $Pr$ values, it is also (at least in part) due to the limitation on the boundary conditions, as discussed in Section 4.1. Some studies [32,45] indicated the possibility that at sufficiently high $k_s^{+}$ various scalar statistics may become independent of $Pr$ when scaled accordingly. A much wider range of $Pr$ is required, in particular $Pr\gg 1.0$, to see if the various findings and trends, identified for the much more limited range of $Pr$, are valid and applicable at much higher $Pr$.

4.3. Effect of Blockage

Roughness generates drag and scalar flux, along with blockage of the flow [2]. As mentioned above, Jiménez [28] recommends roughness no higher than $1/40$ of the channel half-height or boundary-layer thickness.

This limitation seems to be less strict in the transitionally rough regime; Thakkar et al. [78] demonstrated that $k/\delta \approx 1/6$ hardly altered the flow in the outer layer by using 17 different irregular surfaces. Even if the outer flow remains mostly unaffected, a large blockage ratio will modify the effective boundary-layer height (${\delta }_{\mathrm{eff}}<\delta$), impact $k_s^{+}/k^{+}$, and/or modify $u_{\tau }$, etc. [28,57].

For a given roughness topology, increasing $k_s^{+}$ can be achieved by either increasing $Re$ or by increasing the physical roughness height, k. The former is associated with higher computational costs, leading authors to opt for the latter option. The latter option, however, increases the blockage; if the flow is significantly constricted, it should be considered flow over obstacles rather than roughness.

The limits on the blockage ratio are acknowledged by the various studies but it is not uncommon for numerical simulations to reach a blockage factor, $k/\delta$, of $0.2$ or more, which is significant. In contrast, the work performed on minimal channels by MacDonald et al. [30], Zhong et al. [45], Rowin et al. [57] achieves better blockage ratios, with $k/\delta$ between $0.015$ and $0.05$. The sand-grain roughness of Hantsis and Piomelli [31], Hantsis [32], Hantsis and Piomelli [34] also maintained $k/\delta =0.04$.

4.4. Virtual Origin

In a smooth-wall case, the outer turbulent flow “perceives” its origin to be at $y=0$, where the no-slip condition is applied. This is not the case when the wall is rough, and the “virtual origin” of the flow is located somewhere between the base and the crest of the roughness [23,99]. The virtual origin is sometimes referred to as “wall offset” or “zero-plane displacement”, and is commonly denoted by d or $h_m$. This is equivalent to the “canopy penetration depth” concept used in atmospheric flows and dense vegetation flows.

To obtain robust measures of the logarithmic intercepts, the wall distance need to be measured relative to this virtual origin [23,28,99]. The wall distance thus modified, $y_e$, yields the required collapse of various rough-wall statistics on the corresponding smooth-wall ones in the outer region, satisfying Townsend’s outer-layer similarity [45].

While the concept is widely recognised, there is debate over the proper choice of the virtual origin, which manifests itself in various definitions of d. Two of the most common ones are based either on the roughness geometry or on the center of drag. The former commonly defines d as the geometrical mean of the roughness, and is insensitive to the conditions of the flow. The latter usually defines the virtual origin as the displacement that situates the centre of drag acting on the rough surface [99]. Both definitions are common throughout the literature and generally result in similar values. However, this creates uncertainty, which propagates into the results and statistics, in particular when comparing studies that use different definitions for d.

One limitation of the centroid method suggested by Jackson [99] is that it assumes a zero pressure gradient and fully developed flow. This was partially amended by [100], who suggested an extension to include pressure gradients, which seems promising.

Zhong et al. [45] studied the variation of $d/k$, and concluded that the influence of the uncertainty in d is largely inconsequential when measuring the logarithmic intercept, which translates into small variations in the roughness functions for velocity, $\Delta U^{+}$, and scalar, $\Delta {\Theta }^{+}$. Thus, while there is a level of unreliability when comparing between studies that use different definitions of d, it is not significant enough to affect the conclusions.

Another geometry-based approach to define the centroid-shifted coordinate was suggested by [33] and used by [52,60,101]. The effective wall-normal distance from the rough surface is introduced as

$$y_e={\int }_0^y\varphi dy$$

(26)

where the plane porosity $\varphi =A_f/A_0$, and $A_f$ is the plane area occupied by the fluid. While this formulation produces the same shifted coordinate, $y-d$, outside the rough wall, its origin remains at zero. This is in contrast to the previously mentioned simple shifting, whose origin becomes negative (i.e., $y_e=y-d<0$).

The uncertainty of d also influences the ratio of $ks/\delta$ [57], which may cause difficulties properly determining $k_s^{+}$. Not only does the choice of virtual origin tend to have more significant impact especially for larger $k_s^{+}$ [16,45], but also there is no reason for the scalar virtual origin, $d_{\theta }$, to correspond with the momentum one [57] (despite it being common practice).

While these differences may not change trends and conclusions, they can affect the resulting values. For example, both MacDonald et al. [30] and Rowin et al. [57] used the same 3D sinusoidal roughness at $Pr=0.71$ and identified a flattening of $\Delta {\Theta }^{+}$ at the ultimate regime. Using different approaches to determine the virtual origin (assuming $d_{\theta }=d$ in both studies), [30] reported $\Delta {\Theta }_{\mathrm{FR}}^{+}\approx 4.4$, while [57] reported $\Delta {\Theta }_{\mathrm{FR}}^{+}\approx 3.8$—a difference of approximately $13\%$, which is not insignificant. This may cause discrepancies between studies that are not physical, but rather due to a methodological choice.

A similar problem also exists with respect to experimental studies. Therefore, a more universal and consistent method of determining d and $d_{\theta }$ is desirable.

5. Conclusions

Roughness is known to have significant effects on the transport of momentum and scalar quantities. Despite significant research, our understanding of the latter lags behind that of the former. Wall-resolved numerical simulations provide unfettered access to details of the flow, both for momentum and scalar; this is particularly important when experimental measurements are difficult to obtain (e.g., dense roughness). This provides a tool to study the physics of the flow inside the roughness sublayer (RSL) and examine the physical assumptions taken by various phenomenological modes. Moreover, numerical simulations provide an effective technique for systematic studies and sweeping across the parameter spectrum. Here we summarise the progress in the study of passive scalars in flows over rough surfaces, thanks to recent numerical simulations.

Wall-resolved numerical simulations provided a high-fidelity description of the instantaneous scalar field, showing that the scalar field exhibits a diffusive sublayer that increasingly conforms to the roughness surface as $k_s^{+}$ or $Pr$ increase. The diffusive sublayer is thinnest on the windward roughness surface and around the crest, resulting in sharp gradients that enhance the scalar wall flux. The Reynolds analogy was found to hold locally on the exposed portion of the windward slope; it breaks on the sheltered portion and in the recirculation regions behind the roughness elements, where the scalar wall flux is attenuated. The differences between the velocity and scalar fields are focused on the leeward side of the roughness elements and in the recirculation region behind them. On the windward side, the rough-wall diffusive sublayer is analogous to a contorted, smooth-wall diffusive sublayer (at sufficiently high $k_s^{+}$); a similar analogy cannot be made on the leeward side.

The downwards shift of the logarithmic region between the smooth- and rough-wall profiles of the mean scalar is quantified by the scalar roughness function, $\Delta {\Theta }^{+}$, which is equivalent to the Clauser–Hama roughness function, $\Delta U^{+}$. It was observed that $\Delta U^{+}<\Delta {\Theta }^{+}$ for sufficiently large $k_s^{+}$, regardless of the Prandtl number, $Pr$. For $Pr>1.0$, there exists a finite range of $k_s^{+}$ for which the opposite is true; the value of $k_s^{+}$ for which $\Delta {\Theta }^{+}$ intersects the velocity asymptote increases with $Pr$. Various numerical studies observed a flattening of $\Delta {\Theta }^{+}$ as $k_s^{+}\to \infty$, suggesting the possibility that an “ultimate regime”, a scalar equivalent of the fully rough regime, exists. This regime is characterised by a finite asymptotic value $\Delta {\Theta }_{\mathrm{FR}}^{+}$, which depends strongly on the Prandtl number, $Pr$, and weakly on the roughness geometry. While this behaviour deviates from experimental findings, it is consistent with the observed topography-conforming behaviour for the scalar field. Further investigations are required to clarify the issue.

The scalar roughness function, $\Delta {\Theta }^{+}$, is related to the total scalar wall flux and the Stanton number, $St$. As such, $St$ also depends on $Pr$ and the roughness parameters such as the frontal solidity, ${\lambda }_f$. A peak in $St$ was found for $0.1<{\lambda }_f<0.3$, usually corresponding to the maximum of the skin-friction coefficient, $c_f$. This is a result of an interplay between the impingement of bulk fluid and the sheltering by upstream roughness elements; the earlier enhances $St$, while the latter attenuates it. The rough-wall $St$ is greater than the smooth-wall case, ${St}_0$, and approaches it in the limit of infinite roughness density. The existence of an ultimate regime implies that the enhancement of the scalar transport will saturate and will be accompanied by increasing drag, resulting in the Reynolds analogy factor $RA<1$.

The rough-to-smooth ratio of the Reynolds analogy factor, $\eta$, is an efficiency parameter indicating the augmentation to $St$ relative to the augmentation to $c_f$. Roughness results in $\eta <1$, which tends to decrease with $k_s^{+}$ until reaching saturation; it was shown to be a convex function of the roughness density and ${\lambda }_f$. Predicting the ${\lambda }_f$ value of the minimum for different types of roughness is still an open question.

The wake velocity, $\tilde{u}$, and scalar, $\tilde{\theta }$, fields are overall similar throughout the RSL, with differences observed in the troughs between roughness elements, where $\tilde{u}$ and $\tilde{\theta }$ tend to have opposite signs. The mean vertical scalar dispersive flux, $\langle \tilde{\theta }\tilde{v}{\rangle }_s$, is larger than $\langle \tilde{u}\tilde{v}{\rangle }_s$ for $Pr\ge 1.0$. It can account for a significant portion of the mean vertical scalar flux inside the RSL and is also important for the production of turbulent fluxes. Overall, the dispersive stresses, $\langle {\tilde{u}}_i{\tilde{u}}_{k}j{\rangle }_s$, fluxes, $\langle \tilde{\theta }{\tilde{u}}_j{\rangle }_s$, are not well studied and may be a very useful avenue of numerical exploration, since they are very difficult to measure experimentally. Additionally, they seem to be less sensitive to $k_s^{+}$ and even reach saturation with either $k_s^{+}$ or $Pr$. Dispersive fluxes are also needed for closing models that are based on the discrete-element method or the double-averaged NS equations (DANS).

When scaled in inner units, both the velocity and scalar fluctuations are attenuated by roughness; the former is attenuated more significantly than the latter. This may be due to the meaningful difference in the production of the velocity and scalar fluctuations, particularly in the form-induced (FI) production term, which produces fluctuations at a scale comparable to the roughness scale, k. As $k_s^{+}$ increases, the scalar FI production, ${\mathcal{P}}_{\theta }^{fi}$, becomes dominant in the RSL and is much more significant than the corresponding FI production of TKE, ${\mathcal{P}}_{kk}^{fi}$. ${\mathcal{P}}_{\theta }^{fi}$ increases with $Pr$ and does not seem to saturate; however, for a fixed $Pr$, the ${\mathcal{P}}_{\theta }^{fi}$ curves collapse as $k_s^{+}$ increases. The shear production component remains similar between the TKE and scalar variance; it continually decreases with $k_s^{+}$ and may reach saturation. A possible result is that for a fixed $Pr$ the scalar variance budget becomes insensitive to $k_s^{+}$. Asymptotic trends such as these are useful, in particular in the context of modelling. A more systematic study is required to ascertain these trends at higher $k_s^{+}$, and for a wider range of $Pr$ and various types of roughness.

Both the turbulent Prandtl number, ${Pr}_T$, and the effective turbulent Prandtl number, ${Pr}_{T,\mathrm{eff}}$, are ill-behaved near the bottom of the rough surface. Neither seem to be an effective modelling tool in generalised rough-wall flows and an alternative might be required for modelling.

Several key issues, however, remain open. Despite providing an effective technique for systematic studies and sweeping across the parameter spectrum, the computational cost associated with higher Reynolds and Prandtl numbers has limited the parameter range covered by numerical simulations. One way to circumvent this limitation and reach higher $k_s^{+}$ is to increase the roughness height, k, at relatively low $Re$. This approach, however, results in excessive blockage, which can have a detrimental effect on the simulation, turning it effectively to a flow over obstacles with significant acceleration.

The boundary conditions commonly used to implement constant scalar flux at the wall assume negligible fluctuations at the wall surface. These assumptions are only valid for relatively low Prandtl numbers, say $Pr<5$; these may also be invalid at sufficiently high Reynolds numbers.

Finally, a common problem with numerical and experimental studies is the definition of the virtual origin, d, or the location of the zero-plane displacement, which is where the outer turbulent flow “perceives” its origin to be. Various definitions of d are employed, such as the geometric average or the location where the center of drag is acting on the rough surface; the latter is sensitive to the flow conditions, while the former is not. Both definitions are common throughout the literature and all the statistical quantities are measured with respect to the displaced coordinate. This creates uncertainty, which propagates into the results and statistics, in particular when comparing studies that use different definitions for d.