Contrasts Between Momentum and Scalar Transport Over Very Rough Surfaces

TL;DR

This study uses large-eddy simulations to compare momentum and scalar transport over rough surfaces, revealing significant differences in dispersive fluxes and the influence of surface geometry on flow characteristics.

Contribution

It provides new insights into the contrasting behaviors of momentum and scalar transport, especially dispersive fluxes, over very rough surfaces in turbulent flows.

Findings

Dispersive fluxes are a significant part of total fluxes within roughness sublayers.

Turbulence transports momentum more efficiently than scalars, but dispersive fluxes show the opposite trend.

Surface geometry, especially frontal density, strongly affects fluxes and flow regimes.

Abstract

Large-eddy simulations are conducted to contrast momentum and passive scalar transport over large, three-dimensional roughness elements in a turbulent channel flow. Special attention is given to the dispersive fluxes, which are shown to be a significant fraction of the total fluxes within the roughness sublayers. Based on point-wise quadrant analysis, the turbulent components of the transport of momentum and scalars are found to be similar in general, albeit with increasing dissimilarity for roughnesses with low frontal blockage. However, strong dissimilarity is noted between the dispersive momentum and scalar fluxes, especially below the top of the roughness elements. In general, turbulence is found to transport momentum more efficiently than scalars, while the reverse applies for the dispersive contributions. The effects of varying surface geometries, measured by the frontal density, can be pronounced on both turbulent and dispersive fluxes. Increasing frontal density induces a general transition in the flow from a rough boundary-layer type to a mixed-layer-like type. This transition results in an increase in the efficiency of turbulent momentum transport, but the reverse occurs for scalars due to reduced contributions from large scale motions in the roughness sublayers. This study highlights the need for distinct parameterizations of the dispersive scalar fluxes and the importance of considering the contrasts between momentum and scalar transport for turbulent flows over very rough surfaces.

Figures

Tables

Table 1. Table 1: Summary of simulation parameters. N i b superscript subscript 𝑁 𝑖 𝑏 N_{i}^{b} is the number of grid points resolving one obstacle in the i 𝑖 i direction, while N i subscript 𝑁 𝑖 N_{i} is the total number of computational points in that direction. L x subscript 𝐿 𝑥 L_{x} and L y subscript 𝐿 𝑦 L_{y} are the domain dimensions in the streamwise and cross-stream directions. The two digits at the end of each case name are λ f × 100 subscript 𝜆 𝑓 100 \lambda_{f}\times 100 for that case. Sf stands for ‘Stagger frontal’.

Case	${\lambda }_p$	${\lambda }_f$	$N_x^b,N_y^b,N_z^b$	$N_x,N_y,N_z$	$L_x/\delta$	$L_y/\delta$
Cube25	0.25	0.25	08,08,08	192,096,064	3	1.5
Slender06	0.12	0.06	16,03,08	200,0100,064	3.125	1.5625
Sf08	0.12	0.08	12,04,08	200,0100,064	3.125	1.5625
Sf12	0.12	0.12	8,06,08	200,0100,064	3.125	1.5625
Sf16	0.12	0.16	6,08,08	200,0100,064	3.125	1.5625
Sf24	0.12	0.24	4,012,08	200,0100,064	3.125	1.5625
Wide32	0.12	0.32	03,016,08	200,0100,064	3.125	1.5625

Equations

$$\frac{\text{p}{u}_i}{\text{p}{x}_i}=0,$$

$$\frac{\text{p}{u}_i}{\text{p}t}+u_j\left(\frac{\text{p}{u}_i}{\text{p}{u}_j}-\frac{\text{p}{u}_j}{\text{p}{u}_i}\right)=-\frac{\text{p}p}{\text{p}{x}_i}-\frac{\text{p}{\tau }_{ij}}{\text{p}{x}_j}+F_i+B_i,$$

$$\frac{\text{p}\theta }{\text{p}t}+u_i\frac{\text{p}\theta }{\text{p}{x}_i}=-\frac{\text{p}{q}_i^s}{\text{p}{x}_i}.$$

Full text

1 Introduction

The dynamics of turbulent shear flows over rough walls has been an active area of research because of its relevance in the design of engineering systems and in environmental fluid mechanics. Momentum and scalar exchanges between the wall and the fluid in such flows are of interest in a wide range of disciplines (Belcher et al., 2012). Previous field experiments over natural vegetation (Poggi et al., 2004; Katul et al., 1997b), wind tunnel studies over obstacles of regular shapes (Macdonald & Griffiths, 1998; Castro et al., 2006), and numerical simulations over three-dimensional roughness elements (Coceal et al., 2007a; Finnigan et al., 2009; Leonardi et al., 2015) have advanced our understanding of this problem significantly. For example, these past studies have underlined the importance of dispersive fluxes inside and close to the roughness elements (Poggi & Katul, 2008; Poggi et al., 2004; Jelly & Busse, 2018). Nevertheless, there remains significant gaps in our knowledge, particularly concerning the transport of scalars and how it compares to that of momentum, over surfaces that consist of “large” three-dimensional bluff-body-type roughness elements. Large here implies that the roughness protrudes significantly into the inertial layer and that the details of the flow below the top of the roughness elements, often called the canopy sublayer, are important for the application. In particular, Jiménez (2004) in his review paper limited the discussions of rough-wall boundary layers to $H/\delta<0.025$, where $H$ is the roughness element height and $\delta$ is depth of the boundary layer. In many natural settings and engineering applications, $H/\delta$ often exceeds 0.1, a regime sometimes termed the very-rough surface (Castro et al., 2006).

Various previous numerical studies (e.g. Kanda et al. (2004); Castro et al. (2006); Coceal et al. (2007a); Orlandi & Leonardi (2008); Anderson et al. (2015); Giometto et al. (2016); Anderson (2016); Li et al. (2016a, b)) have probed the details of these very-rough-wall flows such as the morphology of coherent structures and the effects of the roughness in such regimes. Moreover, recent large eddy simulations (LES) and direct numerical simulations (DNS) (Finnigan et al., 2009; Boppana et al., 2010, 2012; Philips et al., 2013; Park & Baik, 2013; Leonardi et al., 2015; Li et al., 2016b) investigated the transport of scalars. However, compared to the extensive number of previous studies focusing on the flow and the momentum transport, research on scalars and their transport in flows over very rough walls remains quite limited. In addition, a detailed comparative analysis of momentum and scalar transport dynamics has not been performed before. These open gaps motivate this present paper: we investigate both momentum and scalars at very-high Reynolds numbers over three-dimensional, large roughness using the LES technique. The region within the roughness elements is defined as the canopy sublayer; the region just above the canopy where the flow is horizontally inhomogeneous is called the roughness sublayer. This paper focuses on these two regions. A logarithmic layer will exist further aloft if the roughness sublayer does not extend all the way to the top of the inertial layer. In particular, we focus on two aspects of the problem, namely effects of different roughness geometries and comparisons between momentum and scalar transport. After descriptions of the numerical setup of the problem in Section two, Section three analyzes the changes in turbulent flow characteristics for different surface geometries; the spatially-coherent dispersive fluxes are then investigated with a focus on the differences between momentum and scalar transport. Section four concludes with a summary and discussion.

2 Numerical setup

The LES code uses the immersed boundary method (IBM) to account for the presence of the roughness elements, in which a discrete-in-time momentum forcing is used to simulate the immersed boundary force (Tseng et al., 2006; Chester et al., 2007; Li et al., 2016a). The non-dimensional filtered incompressible continuity (1), Navier-Stokes (2), and scalar conservation (3) equations are solved assuming hydrostatic equilibrium of the mean flow:

[TABLE]

[TABLE]

[TABLE]

All the variables we will discuss are filtered components, so the usual tilde above the symbols is omitted for simplicity. The density $\rho$ is taken to be unity and is uniform: buoyancy forces are not considered. $x$, $y$ and $z$ denote the streamwise, cross-stream and wall-normal directions respectively, and $u$, $v$ and $w$ are the velocity components in these respective directions. $t$ denotes time; $u_{i}$ is the resolved velocity vector; $x_{i}$ is the position vector; $p$ is a modified pressure (Bou-Zeid et al., 2005); $\tau_{ij}$ is the deviatoric part of the subgrid stress tensor; $B_{i}$ is the immersed boundary force representing the action of the obstacles on the fluid; and $F_{i}$ is the body force driving the flow (here simply a homogeneous steady horizontal pressure gradient along the $x$ direction with magnitude $u_{*}^{2}/\rho\delta$). A friction velocity of $u_{*}=1$ m s*-1*, and a half channel width of $\delta$=100 m, are used to non-dimensionalize $t$ and all outputs of the code; their imposed numerical values are thus inconsequential here. The top boundary is impermeable with zero stress; the simulations therefore are similar to half a channel. In Eq. (3), $\theta$ denotes a passive scalar quantity, which for illustration is considered to be temperature in Kelvin in the current simulation, and $q_{i}^{s}$ is the $i$th component of the subgrid scale scalar flux. Further numerical details on the code and the subgrid-scale model can be found in Bou-Zeid et al. (2005), while detailed validations for the flow and scalar transport can be found in Li et al. (2016a) and Li et al. (2016b), respectively. The LES uses a wall-model for momentum and scalars that has been developed for a hydrodynamically smooth walls (here each facet of a building/cube is such a smooth wall) at high Reynolds numbers based on Kader & Yaglom (1972). A constant surface temperature scalar boundary condition is used. Wall modelling for complex walls remains an on going area of research (Yang et al., 2015) and an open challenge (Bose & Park, 2018); however, the performance of the current approach has been evaluated and shown to be quite satisfactory in Li et al. (2016b).

To explore the effects of the geometry, we conducted simulations of different cases as summarized in table 1 and illustrated in figure 1. $\lambda_{f}$ and $\lambda_{p}$ are the frontal area ratio and plan area ratio, respectively. $\lambda_{f}$ is defined as the total projected frontal (mean-flow normal) area of the roughness elements per unit wall-parallel area (i.e. land area); $\lambda_{p}$ is the ratio between the crest plan area (i.e. roof area, shaded in black figure 1) and the wall-parallel area. The first case is the classic cuboid obstacles. In the remaining cases, we gradually change the horizontal aspect ratio of the obstacles, while maintaining the same height and area ($\lambda_{p}$). Figure 1 shows the cubic case and the two extreme aspect ratios of the other cases. All cases shown in Table 1 were run for 50 eddy turn-over times ($T\approx\delta/u_{*}$) and averaged for the last 25 $T$. The averaging time is comparable to that in Leonardi & Castro (2010), where averaging time is 2000$T_{b}$ for $T_{b}=H/U_{b}$; $U_{b}$ is the bulk streamwise velocity and $H$ is the obstacle height. Throughout the paper, $\left\langle X\right\rangle$ defines a volume average of $X$ over $L_{x}L_{y}\Delta z$, for $\Delta z$ being the vertical grid size; $\overline{X}$ defines temporal average of $X$.

A sensitivity test was carried out to examine the effect of the number of points used to represent the obstacles and the domain length, particularly for cases where one dimension of the cubes is only spanned by three grid points. The comparisons for case Wide32 between the original simulation, a doubled resolution run and a run with doubled $L_{x}$ are shown in the appendix. Both simulations show converged results in the canopy sublayer and roughness sublayer. We thus conclude that the numerical resolution is sufficient for the present purposes.

3 Results and discussion

3.1 Turbulent transport

3.1.1 Effects of surface geometry

We first focus on the impact of surface geometry on the general turbulent flow characteristics. The parameter space that characterizes the surface geometry is large. A non-exhaustive list of studies includes previous works that have investigated height variations (Yang et al., 2016), geometric shapes (Leonardi et al., 2015; Yang et al., 2016; Llaguno-Munitxa et al., 2017; Llaguno-Munitxa & Bou-Zeid, 2018) of regular (Kanda et al., 2004; Placidi & Ganapathisubramani, 2015) or irregular (Chester et al., 2007; Yuan & Jouybari, 2018) surface roughness elements, as well as statistical moments of roughness elements (Zhu et al., 2017). We do not aim to comprehensively examine the parameter space in this paper, but we are more interested in the general transition of the flow as the roughness, conserving the same area density $\lambda_{p}$, changes from slender elements with low $\lambda_{f}$ to wide ones with high $\lambda_{f}$. Ghisalberti (2009) analyzed the dynamics of flows over many different types of canopies and used the term ‘obstructed shear flow’ to characterize them. It is found that this type of canopy flow can be characterized by the penetration depth ($\delta_{e}$) of the vortices into the canopy sublayer. Figure 2a shows the shear scale $L_{s}$ as a function of the frontal blockage ratio $\lambda_{f}$. $L_{s}$ is a basic length scale in canopy flows similar to the vorticity thickness defined in a plane mixing layer; it is defined as $L_{s}=\frac{\langle\bar{u}\rangle}{\langle\bar{u}/dz\rangle}$, here computed at $z=H$. $L_{s}$, which is analogous to the penetration depth $\delta_{e}$, decreases with increasing $\lambda_{f}$ (figure 2a), which signifies a larger shear strength and hence more deviation from the classical surface layer profile as the blockage increases. $L_{s}$ for Slender 06 and Wide32 differ significantly and additional intermediate cases show a monotonic decrease with $\lambda_{f}$.

The correlation coefficient $r_{uw}$ computed as $\left\langle\frac{\overline{u^{\prime}w^{\prime}}}{\sigma_{u}+\sigma_{w}}\right\rangle$ at $z=H$ is also depicted in figure 2a. Its magnitude can be interpreted as a vertical turbulent momentum transport efficiency (Li & Bou-Zeid, 2011). As $\lambda_{f}$ increases, $-r_{uw}$ increases approximately from 0.2 to 0.46, the latter being consistent with what is typically found in canopy flow and mixing layer flow ($-r_{uw}\approx 0.5$) (Finnigan, 2003). It is worth noting that obstacles in all the cases presented in figure 2 have the same plan area density (i.e. $\lambda_{p}=0.12$); but by changing $\lambda_{f}$, we are observing a transition from a canonical boundary-layer flow to a mixing-layer one where momentum turbulent transport is more efficient. As this transition occurs, and unlike for its momentum counterpart, the vertical transport efficiency for scalar exhibits a non-monotonic behaviour (figure 2b), reaching a peak of approximately 0.23 at around $\lambda_{f}=0.12$. For large $\lambda_{f}$, while momentum $r_{uw}$ is more than twice larger than its value at $\lambda_{f}=0.06$, the scalar $r_{\theta w}$ changes less significantly with changing geometry. Overall, the results show that as the flow becomes more mixing-layer like, the momentum and scalar transport efficiencies increasingly diverge. The ratio of their correlation coefficients become appreciably larger than 1 (figure 2b), making the Reynolds analogy (which postulates that momentum and scalars are transported similarly) less applicable. Note that we do not imply that the very rough canopy or vegetation canopy completely resembles the canonical mixing layer flows. Rather, the examples shown in the present study demonstrate that the rough wall-bounded flows can exhibit mixed properties of wall-bounded and mixing-layer flows (Kanda, 2006).

To gain a better understanding of the turbulent statistics at representative points with respect to the roughness elements, we sampled data at four different points $x_{1}$-$x_{4}$ (see figure 1) for all heights. The point-wise time series of relevant quantities are recorded at a frequency of $1/(2500T)$, where $T$ is the eddy turn-over time defined previously. Data were sampled at each time step for a total period of 2$T$. This time averaging alone is not sufficiently long to guarantee complete statistical convergence of the results. However, the statistics were also averaged for all repeating units to improve statistical convergence and the primary characteristics of the flow that we will focus on are already very clear with this limited averaging time. We first examine the skewness of the fluctuating components shown in figure 3 for point $x_{3}$ only for $u/u_{*}$, $w/u_{*}$ and $\theta/\theta_{1}$, where $\theta_{1}=\theta(z=\delta)$. We only illustrate three cases with the same $\lambda_{p}$, spanning the full range of $\lambda_{f}$. Slender06 has close to zero skewness for $u^{\prime}$ and $w^{\prime}$ (figures 3a and b) in the roughness sublayer but higher values are observed as $\lambda_{f}$ increases in cases Sf16 and Wide32, more typical of canopy flows from experimental measurements over vegetation canopies (Raupach, 1981; Rotach, 1993; Finnigan, 2003). The variation in skewness presents further evidence of the role of frontal area, and hence pressure drag, in dictating the characteristics of the roughness sublayer dynamics. The negative skewness for $\theta^{\prime}$ in figure 3c signifies dominance of downward events that tend to bring down cooler fluid regardless of the underlying surface geometry (this is also consistent with the negative skewness of $w^{\prime}$ in figure 3b). This is similar to findings for temperature, as a passive scalar, over rough surfaces in the DNS by Orlandi & Leonardi (2008).

Since in physical space we observed that changes in surface roughness modify the characteristic length scale and turbulence statistics, we now examine in spectral space if that also implies a modification in the range of length scales that contributes the most to the turbulent transport. In spectral space, the correlation spectrum is defined as $R_{XY}=Co_{XY}/\left(\Gamma_{X}\Gamma_{Y}\right)$, where $Co_{XY}$ is the cospectrum of time series $X$ and $Y$ and $\Gamma_{X}$ and $\Gamma_{Y}$ are the power spectra of $X$ and $Y$, respectively.

For $Y$ being vertical velocity fluctuation $w^{\prime}$ and $X$ being $u^{\prime}$ or $\theta^{\prime}$, $R_{XY}$ is interpreted as a wavelength specific transport efficiency in the spectral space. The total length of the time series defined in time unit $L_{z}/U_{0}$ is about 60, where $U_{0}$ is the free stream velocity. Using Taylor’s frozen turbulence hypothesis to convert time to length scale with $U_{0}$, and performing ensemble averaging over 20 shorter time series (obtained by sectioning the total time series into 20 sub-series) at each representative point $x_{i}$, the correlation spectra are computed. The results are shown in figure 4, where only results at $x_{3}$ (in front of the obstacle) are presented; similar conclusions can be obtained from the results at other points (not shown here). Eddies of large scales (close to $\delta=L_{z}$) contribute to the efficient transport of both momentum, and scalars, but a wider range of scales contribute to $R_{\theta w}$. Higher correlation spectra values can be noted for scalars, compared to momentum, for all three cases between $H/\delta$ and $3H/\delta$. As the roughness frontal blockage increases (e.g. for case Wide32), the large scale motions ($O(\lambda_{x}/\delta)=1$) can penetrate into the canopy more easily and they become more efficient in transporting momentum below the canopy top (figures 4a to c). For scalars, however, it seems they are most efficient at transport deep inside the canopy at intermediate densities $\lambda_{f}$ (case Sf16), becoming less efficient in case Wide32.

In summary, figure 4 shows the change in length scales of the most-efficient transport eddies across different geometries. As the transition to a mixed-layer-type flow indicated in figure 2 occurs, larger scale motions contribute more effectively to momentum transport in the canopy and roughness sublayer. However, for scalar transport, there seems to be an optimal configuration for the large scalar motions to interact effectively with the lower part of the canopy.

3.1.2 Flow and passive scalar structures

To understand how turbulent momentum and scalar transport behave as the flow regime transition illustrated in the previous sub-section occurs, we analyze the spatial correlations of turbulent quantities. The two-point streamwise correlation for a quantify $s$ at reference point $x$ and for a separation distance $X_{0}$ can be computed as \overline{\rho_{ss}}\big{(}X_{0},x\big{)}=\overline{s^{\prime}\big{(}X_{0}\big{)}s^{\prime}\big{(}x\big{)}}/\big{(}s_{rms}(X_{0})s_{rms}(x)\big{)}, where rms denotes the root mean square of the turbulent fluctuations. These correlations are shown in figure 5 for cases Slender06 and Wide32. The angle of inclination of the contours was calculated following the approach of Leonardi et al. (2015), where $\alpha$ is estimated as the angle between the horizontal direction and the segment $X_{0}X_{1}$, indicated by the green solid line, the point $X_{1}$ being the furthest from $X_{0}$ on the contour $R_{uu}$ =0.3. The increased interactions between the canopy top and the roughness sublayer are indicated by larger $\alpha$, especially $R_{uu}(x_{2})$ in figure 5f compared to figure 5b. Compared to the two-dimensional bars studied in Leonardi et al. (2015), the three-dimensional roughness geometry and different configurations also introduce spanwise variation in the correlation contours. The mixing-layer like flow regime results in increased interactions between the surface and the roughness sublayer for both $u^{\prime}$ and $\theta^{\prime}$. However, the increase is more pronounced for the velocity, indicating that the correlation contours $R_{\theta\theta}$ are less sensitive to the increase of $\lambda_{f}$ than their velocity counterparts. This implies that even for a surface-layer like flow regime, where the coupling between canopy sublayer top and the roughness sublayer aloft is not as strong as that for the mixing-layer like case, the effect of the roughness geometry on passive scalar can extend further up. This result is also consistent with the larger scale of the motions contributing to $R_{\theta w}$ than to $R_{uw}$ in figure 4.

The iso-surfaces of high temperature fluctuation, $\theta^{\prime}_{high}$, are shown in figure 6. $\theta^{\prime}_{high}$ is defined here as $1.5\theta^{\prime}_{rms}$ at any instant. However, here we colour these iso-surfaces by the normalized Reynolds stress and wall-normal scalar flux (figures 6b, d and figures 6 e, f, respectively). First, although both Slender06 and Wide32 demonstrate streamwise thermal ‘streaks’ (Hetsroni et al., 1999; Leonardi et al., 2015), Wide32 is characterized by more patchy structures in most of the domain; whereas the temperature iso-surfaces in Slender06 are more elongated in the streamwise direction and are aligned in between the staggered roughness elements. Stronger flow disturbances by the obstacles are observed in the mixing-layer-like flow where the more obstructive three-dimensional geometry is found to significantly disrupt the structures, in agreement with the experimental observations in Hetsroni et al. (1999). Second, there is a general spatial correspondence between $w^{\prime}\theta^{\prime}$ and $-w^{\prime}u^{\prime}$ “events”, especially in the roughness sublayers (see figures 6a and b, $x/H$=10-15; figures 6c and d). Nevertheless, for instantaneous flow structures, there can also be discernible dissimilarities between them in the canopy sublayer, suggesting that some eddies might be carrying momentum but not scalars or vice-versa. Thus, we conduct quadrant analysis in next section to further investigate the (dis)similarity of momentum and scalar.

3.1.3 Quadrant analysis of turbulent momentum and scalar transport

Quadrant analysis is a useful and widely used technique for probing how turbulent motions evolve and transport momentum and scalars in the wall-normal direction (Wallace, 2016). Thus, we apply this technique here to compare the momentum and scalar turbulent transports over the canopy and roughness sublayers. The definition of each quadrant for momentum flux follows previous studies (Katul et al., 1997a, b; Li & Bou-Zeid, 2011). Q1 events are classified as $s^{\prime}>0$ and $w^{\prime}>0$; Q2 as $s^{\prime}<0$ and $w^{\prime}>0$; Q3 as $s^{\prime}<0$ and $w^{\prime}<0$; Q4 as $s^{\prime}>0$ and $w^{\prime}<0$, where $s$ is $u$ or $\theta$. Here a prime denotes the turbulent perturbation of an instantaneous value from its Reynolds (time) average denoted by an overbar. We applied quadrant analysis to time series collected at four representative horizontal locations ($x_{1}$ to $x_{4}$) indicated in figure 1 and at every height, but we will only show results for $x_{2}$ and $x_{3}$ since the other locations convey the same information, which are shown in the appendix. The same averaging procedure detailed for the results of the skewness is used here. For the momentum flux, Q2 and Q4 are termed ejections and sweeps respectively; for scalar flux, Q1 (motion that transports higher concentrations of the scalar upward, the wall being a source here) and Q3 (motion that brings lower concentrations downward) are termed ejections and sweeps, respectively. Compared to the other two quadrants, the ejections and sweeps are the dominant events in transporting momentum and scalars; they contribute to down-gradient transport.

Figure 7 shows the contributions from various quadrants (momentum in (a) and (c); scalar in (b) and (d)) for case Cube25, where the contribution of quadrant $i$ is defined as $\lvert\overline{s^{\prime}w^{\prime}}\rvert_{i}/\ \sum_{i}\lvert\overline{s^{\prime}w^{\prime}}\rvert$. The contribution of each quadrant to scalar fluxes is broadly similar to that of momentum across all points, though slight differences between momentum and scalars are seen inside the canopy sublayer below $z/H=1$. The spatial variation across the four different horizontal location in the canopy sublayer is on the other hand significant, underlining the complexity of the flow within these three dimensional roughness arrays. These differences blend away above the canopy. The number of ejection events exceeds that of sweeps below $z/H=1.25$ across all four points (not all shown here), indicating that below this height sweeping events are stronger compared to the more frequent and less intense ejections. The crossover point at $z/H=1.25$ is consistent with results for the Reynolds stress from DNS reported by Coceal et al. (2007a) for the same underlying roughness geometry of staggered cubes, although they averaged over all points at a given $z$. It has been previously observed that for momentum transport, the sweeps’ contribution to total stress dominates over ejections, under near-neutral static stability and in the roughness sublayers, in vegetation canopies (Dupont & Patton, 2012) and in a realistic urban area characterized by large roughness elements (Rotach, 1993; Christen et al., 2007). Dominance of sweeps in momentum transport over very rough walls, compared to the dominance of ejections over smooth and less-rough counterparts, is in general agreement with the picture proposed and discussed in detail by Raupach et al. (1996) who adopt a mixing-layer analogy for these very-rough walls to account for the differences compared to a typical smooth-wall surface layer. In the outer layer, the dominance of ejections is restored as typical in smooth-wall boundary layers (Adrian, 2007).

Our results for case Cube25 demonstrate that the turbulent transport of a passive scalar exhibits similar behaviour as for momentum. Similarity or dissimilarity between momentum and passive scalar transport over wall-bounded turbulent flows has also been discussed in the literature. For instance, similarity in the transport is found in laboratory results (Perry & Hoffmann, 1976; Nagano & Tagawa, 1988) over smooth surfaces and field measurements in the atmospheric boundary layer over vegetation canopy and in the roughness sublayer over an urban area (Dupont & Patton, 2012; Wang et al., 2014), and in DNS over two-dimensional square-shaped obstacles (Leonardi et al., 2015); whereas dissimilarity is not as extensively reported (Christen et al., 2007). The generalizability of the similarity between turbulent transport of momentum and scalar observed above for case Cube25 to other cases is assessed through quadrant-analysis based turbulent transport efficiencies, defined as $\eta=F_{total}/\left(F_{ejection}+F_{sweep}\right)$ (Wyngaard & Moeng, 1992; Li & Bou-Zeid, 2011), where $F$ is the (total or from a single quadrant) flux of momentum or scalars. The results for cases Cube25, Slender06 and Wide32 are shown in figures 8-10. Figures 8-10(a) and (b) show the vertical profiles of efficiencies at each horizontal location, and an average over the four points $x_{1}$ to $x_{4}$, while figures 8-10(c) shows the ratio of the $\eta_{m}$ to $\eta_{s}$, averaged over the four horizontal locations.

Above $z/H=1.25$, $\eta_{m}-\eta_{s}$ approaches 0 and spatial variability is blended out at all points and across all cases shown in figures 8-10. Below that blending height, case Cube25 (figure 8) shows the most significant spatial variation of efficiencies across four points: $\eta_{m}$ and $\eta_{s}$ for Slender06 (figure 9) vary much less while for Wide32 only $\eta$ at $x_{4}$ deviates appreciably from the other points. Although there are some discrepancies between $\eta_{s}$ and $\eta_{m}$ at some locations, such as at $x_{2}$ and $x_{4}$ within the canopy sublayer of Cube25 and Wide32, the general trends for momentum and scalar efficiency remain similar for those cases. On the other hand, Slender06 displays more dissimilarity between momentum and scalar in the canopy sublayer as shown by both $\eta$ values (figure 9), with a higher transport efficiency for the scalar. While sweeps dominate the transport of the scalar in Slender06 (figure B.2), both sweeps and ejections contribute about equally to the total momentum flux (figure B.1), especially at points $x_{3}$ and $x_{4}$. A physical explanation for this dissimilarity is that the penetration of high horizontal momentum sweeps downward is hindered by the blockage effect of the slender-shaped roughness element, which induces negative fluctuations in $u^{\prime}$ and thus increased inward interactions and reduced $\eta_{m}$. Nevertheless, sweeps that are associated with cold (low $\theta^{\prime}$) fluid but that do not have a large positive $u^{\prime}$ can still efficiently penetrate the canopy and transport the passive scalar. This Slender06 case is characterized by higher $\eta_{s}$ than $\eta_{m}$. Overall, turbulent momentum and scalar transport are broadly similar for Cube 25 and Wide 32, but not for Slender06. It would be instructive to analyze further why the slender blocks induce more dissimilarity than the wider ones, but this is not at the core of the aim of this paper and cannot be pursued here. What we can conclude is that, as the general characteristics of turbulent flows transitions from surface-layer-like to mixing-layer-like regimes, $u^{\prime}$ and $\theta^{\prime}$ exhibit more similar correlations with $w^{\prime}$. For these mixing-layer-like flows, turbulent momentum and scalar transports are broadly similar but with a slightly more efficient momentum transport ($\eta_{m}>\eta_{s}$). However, this similarity is reduced for surface-layer-like turbulence similar and the more efficient exchange is now for scalars with $\eta_{m}<\eta_{s}$). For all geometries, the results suggest a turbulent Schmidt or Prandtl number $\neq 1$.

3.2 Dispersive transport

An important aspect that sets very rough walls apart is that dispersive stresses or fluxes, which arise from the spatial inhomogeneity of the time-averaged flow field, can be important contributors to total transport. In the multiply connected space inside the roughness arrays, the spatial averaging and differentiation operations do not commute (Finnigan, 1985, 2003). Any mean (time-averaged) quantity $\phi$ can be decomposed into $\phi=\left\langle\phi\right\rangle+\phi^{\prime\prime}$, where the double prime represents the spatial deviation of the time-averaged variable from that of the spatial average. In the present paper, we consider $\left\langle\phi\right\rangle$ representing the planar average over an $x-y$ plane at a given height. The dispersive fluxes then arise from the spatial averaging of Eqs. (2) and (3). Therefore, the spatially-local, time averaged dispersive momentum stress is $\overline{u_{i}^{\prime\prime}u_{j}^{\prime\prime}}$ and the dispersive scalar flux is $\overline{u_{i}^{\prime\prime}\theta^{\prime\prime}}$, which can then also be spatially averaged over the plane.

Figure 11 shows the ratios $F_{dis}/F_{total}=\left\langle\overline{u^{\prime\prime}w^{\prime\prime}}\right\rangle/\left(\left\langle\overline{u^{\prime\prime}w^{\prime\prime}}\right\rangle+\left\langle\overline{u^{\prime}w^{\prime}}\right\rangle\right)$ and $F_{turb}/F_{total}=\left\langle\overline{u^{\prime}w^{\prime}}\right\rangle/\left(\left\langle\overline{u^{\prime\prime}w^{\prime\prime}}\right\rangle+\left\langle\overline{u^{\prime}w^{\prime}}\right\rangle\right)$ in (a) and their counterparts $\left\langle\overline{w^{\prime\prime}\theta^{\prime\prime}}\right\rangle/\left(\left\langle\overline{w^{\prime\prime}\theta^{\prime\prime}}\right\rangle+\left\langle\overline{w^{\prime}\theta^{\prime}}\right\rangle\right)$ and $\left\langle\overline{w^{\prime}\theta^{\prime}}\right\rangle/\left(\left\langle\overline{w^{\prime\prime}\theta^{\prime\prime}}\right\rangle+\left\langle\overline{w^{\prime}\theta^{\prime}}\right\rangle\right)$ in (b) for all the cases with $\lambda_{p}=0.12$ and for increasing $\lambda_{f}$ (Slender06 to Wide32).

Examples of previous studies on the momentum dispersive fluxes (Poggi et al., 2004; Christen & Vogt, 2004; Coceal et al., 2007b; Martilli & Santiago, 2007; Poggi & Katul, 2008; Leonardi et al., 2015; Giometto et al., 2016) and a few on dispersive scalar fluxes (Christen & Vogt, 2004; Leonardi et al., 2015) have demonstrated their importance within the canopy sublayer, as well as in the roughness sublayer. Similar to previous findings, our simulations also indicate that the dispersive fluxes are significant within the roughness arrays. The dispersive flux can contribute $\approx$ 50$\%$ of the total momentum or scalar flux (where the contribution of the sub-grid scale fluxes is included in the total) below the roughness elements height for some cases. Interestingly, there are no monotonous trends of how the dispersive fluxes vary with $\lambda_{f}$. The fractions of dispersive fluxes, in general, are the highest for the most eccentric geometries Slider06 and Wide32.

Note that Leonardi et al. (2015) and Coceal et al. (2007b) commented that the dispersive fluxes are only important on the intermediate time scales. We performed our analysis for the time averaged quantities over time spans of approximately $H/U_{b}$ = 1200 ($U_{b}$ is the bulk velocity) and longer than the $H/U_{b}$=600 suggested in Leonardi et al. (2015). Still, temporal averaging for large $H/U_{b}$ gives converging results of the dispersive fluxes that continue to show that they are significant for the canopy and roughness sublayers. Long-lived streamwise rolls can develop and persist in periodic numerical simulations due a “locking” of these rolls as they begin to interact with themselves across the periodic boundaries. Such rolls would increase the estimated magnitude of the dispersive fluxes unless very long time averaging is performed. This might have been the argument of previous studies, but in our simulations the large dispersive fluxes are not a numerical artifact as will be further illustrated below.

Both dispersive stress and scalar flux are important portions of the total fluxes especially within the canopy sublayer, and their differences are now studied in detail. Figure 12 shows $u^{\prime\prime}_{N}=u^{\prime\prime}/\langle u(z=H)\rangle$ and $\theta^{\prime\prime}_{N}=\theta^{\prime\prime}/\Delta\theta$, where $\Delta\theta=\langle\theta(z=H)-\theta(z=0.125H)\rangle$ for an $x$-$z$ cross session indicated by the blue line in figure 1. The pseudocolor plots and streamlines along the horizontal line shown in figure 1 are spatially averaged for all repeating units in the domain. Negative $u^{\prime\prime}_{N}$ in figure 12(a) is due to loss of horizontal momentum $u$ in the wakes produced by the obstacle. While figure 12(d) shows positive $\theta^{\prime\prime}_{N}$ as a result of the surfaces being kept at a higher temperature. As $\lambda_{f}$ increases, in figures 12(b) and (e) and figures 12(c) and (f), the most pronounced momentum-scale distinction can be found upstream of the obstacle where the mean recirculation pattern results in “counter-gradient” dispersive momentum transport (fluid slowed by the pressure field generated by the obstacle being advected downwards) but “down-gradient” scalar transport (cooler or lower concentration fluid transported downwards). This difference emerges from the non-local action of pressure on momentum (the fluid streamwise velocity has to decrease as it approaches the windward face even before it contacts that face), but its absence from the dynamics of scalars (the fluid has to “touch” the surface to uptake scalars). The sign of the dispersive fluxes is thus the same as the total fluxes, while there is partial cancellation of dispersive fluxes for momentum.

The quadrant analysis approach can also be applied to dispersive fluxes to compute a transport efficiency $\eta=F_{total}/\langle F_{ejection}+F_{sweep}\rangle$ for momentum ($\eta^{d}_{m}$) and scalars ($\eta^{d}_{\theta}$). For example, when $u^{\prime\prime}<0$ and $w^{\prime\prime}>0$, we denote it as a dispersive sweep (this is no longer an event but a persistent spatial feature). Figures 13a-c indicate that there is distinct dissimilarity between the vertical distribution of $\eta^{d}_{\theta}$ and $\eta^{d}_{m}$, related to the physical differences that arise from the role of pressure discussed above. In general, it is observed that $\eta_{m}^{d}<\eta_{s}^{d}$, except for case Slender06 shown in figure 13d. This is in stark contrast to the turbulent transport efficiencies shown in figure 8 where the slender case was the one showing a lower turbulent momentum transport efficiency. As $\lambda_{f}$ increases, a non-monotonic trend is observed for the difference between “efficiency” of dispersive transport of momentum and passive scalar, with the most eccentric geometries again showing the highest transport efficiencies. This is in agreement with the observation above that these eccentric cases also had the highest fraction of the fluxes carried by the dispersive part.

4 Conclusion

We use LES to compare and contrast the transport of momentum and passive scalars over very rough surfaces consisting of three-dimensional cuboid roughness elements. The paper focuses on the comparisons between momentum and scalar exchanges, and between turbulent and dispersive contributions, over rough surfaces in the canopy and roughness sublayers. The effect of the frontal blockage of the roughness geometry on the fluxes is also investigated.

The influence of different three-dimensional roughness elements can be seen to alter the general turbulent flow characteristics. We observe a general transitioning behaviour from a surface-layer-like to a typical mixing-layer-like flow or ‘obstructed shear flow’ when the frontal area ratio increases; this ratio $\lambda_{f}$ encodes the flow blockage by the roughness elements (the other widely studied geometric parameter, the planar density of the roughness $\lambda_{p}$, is maintained constant in our study). This is particularly illustrated by the decrease in the shear length scale and the dampened penetration of large eddies into the canopy layer as $\lambda_{f}$ increases. This transition causes a monotonic and significant increase in the momentum transport efficiency as the flow becomes more mixing-layer-like, while on the contrary these transport efficiencies for scalars measured by the correlation coefficient change more mildly and non-monotonically with increasing $\lambda_{f}$. This is mainly a result of the distinct contributions from large-scale motions in the roughness sublayer, which are able to increase the interaction between the canopy sublayer and the roughness sublayer more effectively for momentum than for scalars as the geometry is modified to increase blockage ratio.

Turbulent transports of momentum and passive scalar are found to be similar in general, as evidenced by a quadrant analysis, with higher averaged turbulent momentum transfer efficiency. The unique exception is case Slender06, where the turbulent momentum and scalar transports difference becomes more substantial. However, the dispersive momentum and scalar fluxes show more pronounced differences than their turbulent counterparts (Slender06 again being the exception and displaying the strongest scalar-momentum similarity for the dispersive part). These dispersive fluxes are significant, particularly for the most eccentric geometries (Slender06 and Wide32) where they amount to $\approx 50\%$ of the total fluxes. The least eccentric geometries on the other hand show the most substantial discrepancy between scalar and momentum exchange. The differences between momentum and scalar dispersive fluxes are traced back to the non-local action of pressure in the momentum dispersive contributions, and the absence of a pressure-counterpart for scalar transport. This results in stronger interactions between roughness elements wakes as the surface geometry changes, giving rise to the characteristic recirculation patterns and leading to consistent differences between the dispersive transport efficiencies of momentum and scalar.

Although we only simulated flows with constant scalar boundary conditions, the results about dis(similarity) between turbulent and dispersive fluxes and between momentum and scalars should hold regardless of the scalar boundary condition. On the other hand, if the scalar concentration influences buoyancy (e.g. with large temperature fluctuations), the active role of the scalars will then have a strong impact on the dynamics (Li & Bou-Zeid, 2011) and the present results will be altered. This study demonstrates that the (dis)similarity between momentum and scalar transport in three-dimensional very rough surfaces can be complicated by the the spatial variability of the roughness elements and the specific topology of the underlying geometry. However, the present findings, especially the significance of the scalar dispersive flux contributions over dense canopy, can inform the interpretation of experimental measurements, which are often point-wise data where dispersive fluxes cannot be estimated: such measurement in the canopy and roughness sublayers are missing nearly half of the total fluxes. In addition, in the dense canopy sublayer, the mechanistic difference (arising from the pressure term) in generating dispersive momentum and scalar fluxes need to be incorporated in model development. Current parameterizations either neglect dispersive fluxes or do not distinguish them from turbulent fluxes, but as our analyses show their physics are quite different and this difference is not the same for momentum and scalars.

5 Acknowledgment

This work was supported by the U.S. National Science Foundation’s Sustainability Research Network Cooperative Agreement 1444758 and Grant 1664091. The simulations were performed on the supercomputing clusters of the National Center for Atmospheric Research through projects P36861020 and UPRI0007.

Appendix A Resolution test

A sensitivity test is carried out to examine the effect of the number of points used in representing the obstacles. It was found previously by Tseng et al. (2006) that at least 6 points per dimension of the obstacle are required to represent a solid in LES (they used a similar approach and LES code to the ones we are using here). Therefore, we doubled the resolutions for the case Wide32 listed in table 1. In addition, a test case with all other parameters being the same as Wide32 except that $L_{x}/\delta$ is doubled to 6.25 is also conducted. Figure A.1 shows the comparison between the original simulation, the higher resolution case and the case with longer domain length, where $\theta_{1}$ is $\langle\theta\rangle$ at $z=\delta$ and $\theta_{0}$ is $\langle\theta\rangle$ at $z=1/8H$. Both test cases were run for about 30 eddy turn over times and are averaged for the last 15 eddy turn over times. We compute the height-averaged relative mean-square-error defined as r_{e}^{n}=\bigg{\langle}\frac{\sqrt{(q_{n}-q_{0})^{2}}}{q_{0}}\bigg{\rangle}_{z}, where $q_{n}$ is some quantity for comparison in case $n$ ($n=1$ being doubled $L_{x}$ and $n=2$ being the case of double resolution in all directions) and $q_{0}$ is that same quantity from the initial setup of case Wide32. $r_{e}^{1}$ in percentage for figures (a)-(d) are 1.7, 1.5, 15 and 27 $\%$; $r_{e}^{2}$ in percentage for figures (a)-(d) are 1.7, 2.0, 5.5 and 4.2 $\%$. $r_{e}^{2}$ are all within 6 $\%$ compared to case Wide32, which shows that we achieve good grid convergence. The larger deviation seen in case of doubled $L_{x}$ in the second-order moments is due to the fact that in the larger domain the streamwise roll are less persistent and the convergence to a linearly decreasing stress (as expected) is faster. However, for $z/H<2$, the deviation in the case of doubled $L_{x}$ is small, likely because of the dominance of the small-scale wall-attached eddies.

Figure A.2 shows the percentage of dispersive momentum and scalar flux to total flux for these three cases. Good agreements between the three cases confirm that the differences observed between dispersive momentum and scalar fluxes are robust, especially for the roughness sublayer, which is the focus of this study.

Appendix B Quadrant analysis at different points

Quadrant analysis similar to results in figure 7 are presented here for Slender06 in figures B.1 and B.2; Wide32 in figures B.3 and B.4.