1. Introduction
Closely packed roughness is common in urban boundary layer flows (Barlow & Coceal Reference Barlow and Coceal2008). A prominent example is closely packed buildings in a metropolitan area (Giometto et al. Reference Giometto, Christen, Egli, Schmid, Tooke, Coops and Parlange2017; Krayenhoff et al. Reference Krayenhoff2020). In this work, we study flow over closely packed cubes.
Flows over cube arrays have been extensively studied experimentally and using scale-resolving computational tools (Cheng et al. Reference Cheng, Hayden, Robins and Castro2007; Leonardi & Castro Reference Leonardi and Castro2010; Yang Reference Yang2016; Li & Bou-Zeid Reference Li and Bou-Zeid2019). Depending on the packing density, wall-mounted cubes can be generally categorized into: isolated roughness, k-type roughness and d-type roughness (Oke Reference Oke1988; Jiménez Reference Jiménez2004). In the isolated regime, the roughness elements are sparsely packed and each roughness element acts as an isolated roughness element (see Yang et al. (Reference Yang, Xu, Huang and Ge2019) for a more detailed discussion). In the k-type roughness regime, the packing density/surface coverage density is such that the wakes of upstream roughness elements reaches downstream roughness element. In this regime, the flow in the outer layer actively exchanges momentum with the flow in the roughness-occupied layer, and the equivalent roughness height 
$z_o$ is an increasing function of the packing density 
$\lambda$. In the d-type roughness regime, the packing density is high and the flow in the outer layer does not actively exchange momentum with the flow in the roughness-occupied layer. In this regime, the equivalent roughness height 
$z_o$ is a decreasing function of the packing density. The common expectation is that roughness morphologies in the same roughness regime will have similar hydro/aerodynamic properties and flow features. For example, with closely packed cubes (
$\lambda >0.4$, d-type roughness), the common expectation is that cubes with higher packing densities should have similar properties. Owing to this common expectation, few have studied cubes with packing densities above 
$\lambda =O(0.4)$. The objective of this work is to study these higher coverage densities more extensively than has appeared to date. We study the behaviour of the effectiveness roughness height 
$z_o$ and the displacement height 
$d$ as a function of the surface coverage density for closely packed cubes. We would show that the displacement height approaches the cube height as the surface coverage density approaches one. Another motivation of this work is to study flow in the ‘street canyons’, that is, narrow streets between closely packed buildings. Closely packed buildings block momentum exchange between street canyons and the atmosphere, leading to low wind speed inside these street canyons and poor air quality at the pedestrian level (Li, Liu & Leung Reference Li, Liu and Leung2008, Reference Li, Liu and Leung2005; Huang et al. Reference Huang, Akutsu, Arai and Tamura2000).
In the following, we summarize a few relevant studies on flow over cubical roughness. Cheng & Porté-Agel (Reference Cheng and Porté-Agel2015) and Yang (Reference Yang2016) found that a flat plate turbulent boundary layer (TBL) becomes fully developed shortly after the transition from a smooth wall to a rough wall. When fully developed, roughness directly affects the flow in the roughness sublayer (RSL), above which the flow is nearly homogeneous in the horizontal directions. For closely packed cubical roughness, the height of the RSL is usually thin with 
$\delta _{RSL}<2 h$ (Cheng et al. Reference Cheng, Hayden, Robins and Castro2007; Leonardi & Castro Reference Leonardi and Castro2010), where 
$\delta _{RSL}$ is the height of the RSL and 
$h$ is the cube height. A thin RSL is also observed for other closely packed roughness morphologies, see, e.g., Chan et al. (Reference Chan, MacDonald, Chung, Hutchins and Ooi2015) where the authors conducted DNSs of pipe flow with closely packed sinusoidal bumps on the pipe surfaces. A thin RSL justifies the use of a small vertical domain (Coceal et al. Reference Coceal, Thomas, Castro and Belcher2006; Jiang et al. Reference Jiang, Jiang, Liu and Sun2008). In particular, Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) studied the effects of domain size and concluded that a domain of size 
$4h\times 4h\times 4h$ is sufficient for capturing the mean flow. Blackman, Perret & Calmet (Reference Blackman, Perret and Calmet2016) and Blackman et al. (Reference Blackman, Perret, Calmet and Rivet2017) conducted in-depth analysis of the turbulent kinetic energy (TKE) budget in the RSL. The authors found that the production term in the TKE equation is not necessarily positive. A negative production term transfers energy from the turbulence to the mean flow, and is responsible for the generation of secondary turbulent motions. From a phenomenological standpoint, secondary turbulent motions are spatial variations of the mean flow. For example, spanwise heterogeneity in surface roughness, e.g. spanwise-alternating regions of low and high roughness, gives rise to spanwise-alternating regions of momentum deficit and surplus in the RSL. These regions of momentum deficit and surplus are known as low- and high-momentum pathways (LMPs and HMPs) (Barros & Christensen Reference Barros and Christensen2014). LMPs and HMPs have received significant attention in the recent literature. For example, Mejia-Alvarez & Christensen (Reference Mejia-Alvarez and Christensen2013), Nugroho, Hutchins & Monty (Reference Nugroho, Hutchins and Monty2013) and Barros & Christensen (Reference Barros and Christensen2014) established that any spanwise heterogeneous roughness can give raise to LMPs and HMPs in the mean flow, and Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015) showed that the HMPs and LMPs observed in previous works (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2013; Barros & Christensen Reference Barros and Christensen2014; Willingham et al. Reference Willingham, Anderson, Christensen and Barros2014) are of Prandtl's second kind, that is, they arise due to gradients of the Reynolds stress components.
In this work, we study closely packed cubes in an aligned arrangement. We vary the surface coverage densities from 
$0.25$ to 1, where 
$\lambda =1$ corresponds to a plane channel. We report essential flow statistics including mean velocity profiles, Reynolds stresses, dispersive stresses and roughness properties. Mean velocity profiles and roughness properties have been reported extensively for cubes with low-to-moderate packing densities (Macdonald, Griffiths & Hall Reference Macdonald, Griffiths and Hall1998; Leonardi & Castro Reference Leonardi and Castro2010; Lee, Sung & Krogstad Reference Lee, Sung and Krogstad2011; Volino, Schultz & Flack Reference Volino, Schultz and Flack2011; Yang & Meneveau Reference Yang and Meneveau2016, among others). Here, we augment the literature by studying cubes with high packing densities. It is worth noting that because the width of the thin slots between the cubes are comparable to the viscous length scale as 
$\lambda \rightarrow 1$, our results are inevitably affected by the Reynolds number. The same is true for fractal roughness where the multi-scale nature of the surface roughness dictates that there cannot be a clear scale separation between the roughness length scale and the viscous length scale (Busse, Thakkar & Sandham Reference Busse, Thakkar and Sandham2017). In addition to the basic flow statistics, we put an emphasis here on secondary turbulent motions. We analyse the dispersive kinetic energy (DKE) budget and study the generation and destruction of the secondary turbulent motions. We also determine whether the physics is such that the correct prediction of the secondary turbulent motions depends critically on the correct prediction of the small-scale turbulence.
Before proceeding with the computational details, we define the friction velocity used for normalizing the mean velocity. For a rough wall TBL, there is more than one way to define the friction velocity. Consider a half channel with a rough wall and a free slip top boundary, a common definition of the friction velocity is based on the momentum balance:
\begin{equation} u_\tau=\sqrt{\frac{\text{body force}\times \text{volume occupied by the fluid}}{\text{fluid density}\times \text{planar area}}}=\sqrt{\frac{f_b}{\rho} \left(L_z-\lambda h\right)}, \end{equation}
where 
$L_x$, 
$L_y$ and 
$L_z$ are the extents of the half channel in the streamwise (
$x$), spanwise (
$y$) and wall-normal (
$z$) directions, respectively, 
$f_b$ is the body force and 
$\rho = \textrm {unity}$ is the fluid density. It follows that the definition of the friction Reynolds number is
\begin{equation} Re_\tau=\frac{u_\tau L_z}{\nu}. \end{equation}A second definition leads to the so-called nominal friction velocity
\begin{equation} u_{\tau,N}=\sqrt{\frac{\text{body force}\times \text{volume of the channel}}{\text{fluid density}\times \text{planar area}}}=\sqrt{\frac{f_b}{\rho}L_z}. \end{equation}
A third definition is by Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006). The authors argued that to scale the velocity in the logarithmic layer, the friction velocity must only account for the flow above the virtual ground, that is, above the zero-plane displacement height 
$d$, i.e.
\begin{equation} u_{\tau,C}=u_{\tau,N}\sqrt{1-d/L_z}. \end{equation}
Chan-Braun, García-Villalba & Uhlmann (Reference Chan-Braun, García-Villalba and Uhlmann2011) and Forooghi et al. (Reference Forooghi, Stroh, Schlatter and Frohnapfel2018) proposed to define the friction velocity by extrapolating the total stress from the outer region to the wall, which leads to a fourth definition for 
$u_\tau$. Lastly, we define
\begin{equation} u_{\tau,T}=\sqrt{\frac{f_b}{\rho}(L_z-h)}, \end{equation}
which accounts for the flow above the roughness only, and provides a lower bound for the friction velocity. The ambiguities in the definition of the friction velocity are often overlooked. Although these definitions may end up returning similar values, we think it is important that we explicitly mention the ambiguities involved in the definition of 
$u_\tau$.
Having defined the friction velocity, the mean flow in the inertial layer conforms to
\begin{equation} U^+{=}\frac{1}{\kappa}\ln{z^+}+C-\Delta U^+, \end{equation}
where 
$U^+=U/u_\tau$, 
$z^+=zu_\tau /\nu$, 
$\kappa$ is the von Kármán constant, 
$C$ is a constant and 
$\Delta U$ is the roughness function. Unless otherwise noted, in the following, all velocities are normalized using plus units. An alternative expression of the above scaling reads
\begin{equation} \frac{U}{u_\tau}=\frac{1}{\kappa}\log\frac{z-d}{z_0}. \end{equation}Equating (1.6) and (1.7), we have
\begin{equation} \frac{1}{\kappa}\log\frac{z-d}{z_0}=\frac{1}{\kappa}\log\left[({z-d})^+\right]+C-\frac{1}{\kappa}\log\left[z_0^+\exp(\kappa C)\right], \end{equation}and
\begin{equation} \Delta U^+{\approx} {1}/{\kappa}\log\left[{z_0^+}{\exp(\kappa C)}\right]. \end{equation}
That is, 
$\Delta U$ and 
$z_o$ are correlated. In addition to the flow in the inertial layer, we report the mean flow 
$U$ in the roughness-occupied layer. The definition of 
$U$ involves horizontal averaging in the roughness-occupied layer. There are two ways one can go approximately horizontal averaging: comprehensive spatial averaging (CSA) and intrinsic spatial averaging (ISA). When conducting CSA, the fluid velocity within a roughness element is treated as zero, and the roughness occupied space is not excluded from spatial averaging. On the other hand, when conducting ISA, one excludes the solid volume entirely. For cubical roughness, the fluid occupied horizontal area is 
$A_s(1-\lambda )$ for 
$z< h$ and 
$A_s$ for 
$z>h$, where 
$A_s$ is the planar area. Owing to this discontinuity in the fluid occupied area at 
$h$, ISA leads to a discontinuity in the mean velocity at 
$z=h$ (Xie & Fuka Reference Xie and Fuka2018), which is rather undesirable. Accordingly, in this work, unless otherwise specified, horizontal averaging refers to CSA.
The rest of paper is organized as follows. The computational details are presented in § 2. We report the most essential flow quantities including the mean velocity profiles, the Reynolds stresses, the dispersive stresses and the roughness properties in § 3. It is observed, as hypothesized, that at a surface coverage density of 
$\lambda >0.4$, the roughness fields are d-type and behave qualitatively similarly as aligned cubes with a surface coverage density of 
$\lambda =0.25$. We devote § 4 to secondary turbulent motions. We observe that these motions exhibit significant variation within the d-type regime. Conclusions are given in § 5, and we present further discussion approximately the domain size and the grid resolution in Appendix A.
2. Computational details
Figure 1 shows a sketch of the flow configuration. The computational domain is a half channel. Periodic boundary conditions are imposed in the streamwise and the spanwise directions. The flow is driven by a constant body force 
$f_b$, which is set such that the nominal friction Reynolds number
\begin{equation} Re_{\tau,N}=\frac{u_{\tau,N} L_z}{\nu} \end{equation}
is 500. We note that imposing a constant 
$U_b$ is not the common practice in DNSs of channel flow. The standard practice is to directly impose a constant body force, or to use a constant 
$U_b$ until the flow is fully developed and then switch to a constant body force (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Hoyas & Jiménez Reference Hoyas and Jiménez2006; Graham et al. Reference Graham2016; Yamamoto & Tsuji Reference Yamamoto and Tsuji2018). Details of our DNS simulations are summarized in table 1. We vary 
$\lambda$ from 0.25 to 1. For R100, i.e. for 
$\lambda =1$, the top surfaces of the cubes form a continuous wall at 
$z=h$, and the flow reduces to a plane channel with half channel height 
$L_z=3h$.The half channel height is otherwise 
$L_z=4h$ following Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006). For reference purposes, the computational domain size in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) is 
$L_x\times L_y\times L_z=4h\times 4h \times 4h$, and 
$L_x\times L_y\times L_z=8h\times 6h \times 6h$ in Leonardi & Castro (Reference Leonardi and Castro2010). In § 3, we further justify our use of a somewhat narrower channel. Figure 2(a–c) shows the top view of the rough walls for R25, R50 and R90. The cubes are aligned in both the streamwise and spanwise directions. Figure 2(d–f) show the mesh for R25, R50 and R90. The grid resolution is such that 
$\Delta _n^+< \Delta ^+_{\nu,l}$ at the wall and 
$\Delta z^+\approx 5 \Delta ^+_{\nu }$ in the bulk. Here, 
$\Delta ^+_n$ is the wall-normal grid spacing, 
$\Delta _{\nu, l}^+$ is the locally defined viscous unit, 
$\Delta ^+_{x/y/z}$ is the grid spacing in the 
$x$, 
$y$ and 
$z$ directions, and 
$\Delta ^+_{\nu }=\nu /u_\tau$ is the viscous unit. Table 2 lists the detailed grid information. For reference purposes, the streamwise and spanwise grid resolution is 
$\Delta _{x/y}^+=18.75$ in Leonardi & Castro (Reference Leonardi and Castro2010), 
$\Delta _{x/y/z}^+=7.8$ to 
$15.6$ in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) and 
$\Delta _z^+\approx 8$ to 
$13$ in the recent work by MacDonald et al. (Reference MacDonald, Ooi, García-Mayoral, Hutchins and Chung2018). The resolution of a DNS should be proportional to the Kolmogorov length scale:
\begin{equation} \eta(z)=\left[{\nu^3}/{\epsilon}\right]^{1/4}, \end{equation}
where 
$\epsilon$ is the viscous dissipation. Figure 3(a,b) show our grid spacing in terms of the Kolmogorov length. For reference, the grid spacing is approximately 
$5$ to 
$10\eta$ in a typical channel DNS. Two coarse-grid DNSs are conducted at surface coverage densities 0.5 and 0.8. We compare the two coarse-grid DNSs with their fine-grid counterparts in § A.
Figure 1. A sketch of the flow configuration. The half channel height is 
$L_z$. Here 
$\Delta l$ is the inter-cube distance and 
$x$, 
$y$ and 
$z$ are the streamwise, spanwise and wall-normal directions.
Figure 2. (a–c) Top view of cube roughness topology for R25, R50 and R90. (d–e) Mesh cross-section in 
$x$–
$z$ plane for R25, R50 and R90.
Figure 3. (a) The maximum grid spacing in the 
$x$ and 
$y$ directions, normalized by the Kolmogorov length scale, as a function of 
$z$. (b) Grid spacing in the 
$z$ direction normalized by the Kolmogorov length scale, as a function of 
$z$.
Table 1. DNS details. ‘R’ stands for ‘regular domain size and regular grid resolution’, ‘C’ stands for ‘coarse grid resolution’ and ‘L’ stands for ‘large domain”. The nominal friction Reynolds number, i.e. 
$Re_{\tau,N}=u_{\tau,N}L_z/\nu =500$ is held constant. Here 
$Re_{bulk}=u_{b}L_z/\nu$ is the bulk Reynolds number, 
$u_{b}$ is the bulk velocity, ‘Mesh’ is the total number of grids in million cells (M) and ‘
$N$’ is the number of wall-mounted cubes in 
$n_{{cube},x} \times n_{{cube},y}$, where 
$n_{{cube},x}$ and 
$n_{{cube},y}$ are the numbers of cubes in 
$x$ and 
$y$ directions, respectively.
Table 2. Details of the DNS grids. Here 
$n_{x,y,z}$ is the grid number in 
$x$, 
$y$ and 
$z$ directions, respectively. There is no grid within the wall-mounted cubes.
The DNS solver is the in-house finite-volume code CharLES. It uses a fourth-order accurate spatial discretization and a third-order accurate temporal discretization (Mahesh, Constantinescu & Moin Reference Mahesh, Constantinescu and Moin2004). This code has been extensively used for wall-bounded flow calculations (Ma, Yang & Ihme Reference Ma, Yang and Ihme2018; Yang et al. Reference Yang, Xu, Huang and Ge2019). Further details of the code can be found in Khalighi et al. (Reference Khalighi, Ham, Nichols, Lele and Moin2011) and Bermejo-Moreno et al. (Reference Bermejo-Moreno, Campo, Larsson, Bodart, Helmer and Eaton2014). Time averages are taken over approximately 200 flow-through times after the flow reaches a statistically stationary state. We verify the adequacy of the time averaging by comparing the sum of the Reynolds stress 
$-\langle \overline {u^\prime w^\prime }\rangle ^+$, the dispersive stress 
$-\langle {\bar {u}''} {\bar {w}''}\rangle ^+$ and the viscous stress 
$\langle {\textrm {d}\bar {u}^+}/{\textrm {d}z^+}\rangle$ to 
$1-{z}/{L_z}$, where 
$u$, 
$v$ and 
$w$ are the velocity in the three Cartesian directions, 
$'$ denotes temporal fluctuation, 
$\langle \cdot \rangle$ denotes horizontal averaging, 
$\bar {\cdot }$ denotes time averaging and 
$''$ denotes spatial variation. We also use 
$D_{ij}$ to denote the dispersive stress tensor and 
$R_{ij}$ to denote the Reynolds stress tensor. Figure 4(a–f) show the results. We see that the total stresses follow 
$1-z/L_z$ closely above 
$z=h$. Hence, we can safely conclude that our time averaging is adequate (for first- and second-order statistics).
Figure 4. Reynolds stress 
$\langle \overline {u^\prime w^\prime }\rangle$, dispersive stress 
$\langle \bar {u}' \bar {w}''\rangle$ and viscous stress in (a) R25, 
$(\textit {b})$ R50, 
$(\textit {c})$ R60, 
$(\textit {d})$ R70, 
$(\textit {e})$ R80 and 
$(\textit {f})$ R90. Total stress is the sum of dispersive, Reynolds and viscous stresses. Normalization by the friction velocity 
$u_{\tau,N}$ in (1.3). Bold solid line is 
$1-z/L_z$.
3. Basic flow phenomenology
In this section, we report the most essential flow quantities including the mean velocity profiles 
$\langle \bar {u}\rangle$, the Reynolds stresses 
$R_{ij}$, the dispersive stresses 
$D_{ij}$, the equivalent roughness heights 
$z_o$ and the displacement heights 
$d$. The roughness including R25 and certainly including R50 to R90 are d-type. (Note that staggered cubes with 
$\lambda =0.25$ would belong to k-type roughness.) The common expectation is that the roughness in the same roughness regime would have similar roughness properties. Specifically, the expectation is that 
$\langle \bar {u}\rangle$, 
$R_{ij}$, 
$D_{ij}$, 
$z_o$ and 
$d$ are monotonic functions of 
$\lambda$ within the d-type regime.
3.1. Mean velocity profiles
Figures 5(a) and 5(b) show the time and horizontally averaged velocity profiles above and below 
$z=h$, respectively. The mean velocities above 
$z=h$ follow approximately a logarithmic scaling from 
$(z-h)^+\approx 30$. Below 
$z=h$, the flow is impeded by the closely packed roughness, and the velocity decreases quickly towards the bottom wall. Quantitatively, the mean flow follows the exponential scaling 
$U/U_h=\exp [a(b-z)/h]$ in R25, R50, R60 and R70, where 
$a$ is the attenuation factor, 
$b$ is a constant and 
$U_h=\langle \bar {u}\rangle (z=h)$ (Cionco Reference Cionco1965; Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016).
Figure 5. (a) Time and horizontally averaged velocity profiles above the cube crests. Normalization is by 
$u_{\tau,T}$, which is defined in (1.5). LoW corresponds to 
$\langle \bar {u}\rangle ^+=1/\kappa \log ((z-h)^+)+C$, where 
$\kappa =0.384$ and 
$C=4.6$. 
$(\textit {b})$ Time and horizontally averaged velocity profiles beneath the cube crests. The dashed lines are best fits to the exponential law.
3.2. Roughness properties
An important objective of rough-wall boundary-layer flow modelling is to correlate roughness properties with roughness morphology. Historically, this goes hand in hand with drag partition analysis. In this section, we report the partition of the drag forces and our measurements of rough wall properties. Figure 6(a) shows 
$\tau _S/\tau _R$ in our DNSs, where 
$\tau _S$ is the friction force on the bottom wall and 
$\tau _R$ is the force on the wall-mounted cubes. The roughness are d-type, and the flow skims over the roughness without much interaction with the bottom wall, leading to a small 
$\tau _S/\tau _R$ ratio. A more detailed drag partition is shown in figure 6(b), where we show the breakdown of the drag force into the pressure, the friction on the top surfaces, the friction on the side surfaces and the friction on the bottom wall. The pressure force dominates at 
$\lambda =0.25$. For 
$\lambda >0.6$, the top surface friction dominates, and the flow transits from being fully rough to transitionally rough.
Figure 6. (a) The ratio of the friction force on the bottom wall, 
$\tau _S$, and the drag due to the wall-mounted cubes, 
$\tau _R$, as a function of 
$\lambda$. (b) A breakdown of the drag forces into the pressure force, the friction on the top surfaces, the friction on the bottom wall and the friction on the side surfaces.
Next, we measure the effective roughness height 
$z_0$ and the displacement height 
$d$. We use a von Kármán constant of 
$\kappa =0.384$ following Marusic et al. (Reference Marusic, Monty, Hultmark and Smits2013). The friction velocity is 
$u_{\tau,C}$ following Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006). Jackson (Reference Jackson1981) argued that the virtual ground must be located at the height where the resultant drag force acts
\begin{equation} \textrm{d}=\frac{\displaystyle\int{D
z\,\textrm{d}A}}{\displaystyle\int{D\,\textrm{d}A}}.
\end{equation}
Here, 
$D$ is the drag force, and the integration is from the bottom wall to 
$z=h$. For sparsely packed roughness, where the boundary layer interacts with the flow near the bottom, the argument of Jackson (Reference Jackson1981) is valid, but for a boundary layer that interacts with only the top part of the roughness, (3.1) underestimates the height of the virtual ground. To see that more clearly, we consider flow over an array of slender roughness elements, as sketched in figure 7. The flow is driven by an imposed body force, the same as in our DNSs. A boundary layer forms above the roughness crests, and it interacts with only the top part of the roughness, i.e. the shaded part in figure 7. The flow below the dashed line is similar to that in a porous media, where the body force drives the flow through an array of obstacles. The drag force on the shaded part balances the body force imposed on the fluid in the boundary layer. Hence, the virtual ground of the boundary layer should be located between the dashed line and the crests of the roughness elements, i.e. at 
$d_1$. However, because the fluid near the bottom wall is also subjected to the imposed body force, which is balanced by the drag on the unshaded roughness elements, (3.1) yields a virtual ground at 
$d_2$, which is an underestimate. That is, when spacing between roughness elements is small, the virtual ground should approach the cubes’ top surface (MacDonald et al. Reference MacDonald, Chan, Chung, Hutchins and Ooi2016; Sharma & García-Mayoral Reference Sharma and García-Mayoral2020; Chung et al. Reference Chung, Hutchins, Schultz and Flack2021). That being said, because detailed measurements of the drag force are often unavailable, (3.1) is not widely used, and the zero-plane displacement height is usually such that it yields a log layer in the mean velocity profile (Zhu et al. Reference Zhu, Iungo, Leonardi and Anderson2017; Volino & Schultz Reference Volino and Schultz2018; Womack, Meneveau & Schultz Reference Womack, Meneveau and Schultz2019).
Figure 7. A sketch of flow over an array of slender roughness elements. The flow is driven by a body force. A boundary layer forms above the crests of the roughness elements. Here 
$d_1$ is the height of the virtual ground and 
$d_2$ is the height where the drag force acts.
Figure 8(a) shows the displacement heights computed according to (3.1), and figure 8(b) shows the displacement heights such that they yield a log layer in the mean velocity profile. In the fitting, the log layer is from 
$z/h=1.5$ to approximately 
$(z-d)/h=1$, within which layer the dispersive stress is essentially zero, and the flow is not susceptible to the effects of the top boundary. Compared with 
$d$ that yields a log layer in the mean profile, (3.1) underestimates the zero-plane displacement height 
$d$. Next, we fit for the effective roughness height 
$z_0$. The results are shown in figure 9(a,b). Although the value of 
$d$ depends very much on how it is measured, the effective roughness heights seem to be quite robust. In fact, the displacement heights in figure 8(a,b) lead to essentially the same effective roughness heights. In figures 8 and 9, we have also included measurements in Cheng et al. (Reference Cheng, Hayden, Robins and Castro2007) of aligned cubes at a surface coverage density of 0.25, which gives a rough idea of the uncertainty in the data.
Figure 8. (a) Zero-plane displacement 
$d$ computed according to (3.1). Note that the 
$y$ axis does not start from zero. (b) Zero-plane displacement 
$d$ such that it yields the best log layer. The cross symbols are measurements reported in Cheng et al. (Reference Cheng, Hayden, Robins and Castro2007).
Figure 9. (a,b) The effective roughness heights that best fit the velocity profiles with the zero-plane displacements in figure 8(a,b).
3.3. Reynolds and dispersive stresses
Figure 10 shows the diagonal components of the Reynolds and the dispersive *stress tensors. The streamwise components 
$R_{11}$ and 
$D_{11}$ dominate. The roughness suppresses turbulence for 
$z< h$. As a result, 
$R_{11}$ decreases as 
$\lambda$ increases below 
$z=h$. Above the cube-occupied layer, the peak magnitude in 
$R_{11}$ increases as 
$\lambda$ increases. The 
$R_{11}$ profiles do not depend significantly on the surface coverage density above 
$z/h\approx 1.6$, suggesting a thin RSL. The dispersive component 
$D_{11}$ peaks at 
$z\approx h$, and decreases with 
$\lambda$ in the roughness-occupied layer. Above 
$z=h$, the dispersive stress decays rapidly to zero, again, suggesting a thin RSL.
Figure 10. (a) The normal components of the Reynolds stress tensor 
$R_{ij}=\langle \overline {u'_iu'_j}\rangle ^+$ and the dispersive stress tensor 
$D_{ij}=\langle \bar {u}''_i\bar {u}''_j\rangle ^+$ in (a) R25, (b) R50, (c) R60, (d) R70, (e) R80 and (f) R90. Normalization is by the friction velocity 
$u_{\tau }$ in (1.1). The vertical line is at 
$z/h=1$.
4. Secondary turbulent motions
We pay special attention to the secondary turbulent motions. In this section, we report the basic phenomenology of these motions and determine the physical processes responsible for their generation and destruction.
4.1. Secondary flow phenomenology
Figure 11(a–f) show the contours of the mean streamwise velocity at a constant 
$x$ plane through the centre plane of the wall-mounted cubes. The highlighted contour lines go through 
$y=0$, 
$z=1.2$. The contour lines are relatively flat in R50, R80 and R90, suggesting spanwise homogeneity. The lines are concave in R25, and convex in R60 and R70, suggesting the presence of LMPs at the cube location in R25, and HMPs at the cube locations in R60 and R70. The location of the HMPs and LMPs are more clearly visualized in figure 12, where we show the spatial variations of the time-averaged streamwise velocity, i.e. 
$\bar {u}^{+\prime \prime }=\bar {u}^+-\langle \bar {u}\rangle ^+$, at a distance 
$0.2h$ above the cubes. In R25, a LMP is found at the cube location, and in R60 and R70, a HMP is found at the cube location, consistent with figure 11. Evidently, the positions of the LMPs and HMPs depend on the roughness geometry, and they do not necessarily bring high-momentum fluid to the roughness (Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015). In addition, we see from figure 12 that the flow is (mostly) homogeneous in the spanwise direction in R50, R80 and R90 (figure 12b,e,f). The presence of LMPs and HMPs is usually accompanied by streamwise vorticity. Here, we choose to examine 
$\bar {u}''$ instead of the streamwise vorticity, because the connection between LMP, HMP and the streamwise vorticity is not as direct as that between LMP, HMP and 
$\bar {u}''$.
Figure 11. (a) Contours of the mean streamwise velocity at a constant 
$x$ location through the centre of a wall-mounted cube in (a) R25, (b) R50, (c) R60, (d) R70, (e) R80 and (f) R90. Normalization is by the friction velocity 
$u_\tau$ in (1.1).
Figure 12. (a) Contours of the spatial variation of the mean streamwise velocity at a wall-normal height 
$z=1.2h$, i.e. 
$0.2 h$ above the cubes, in (a) R25, (b) R50, (c) R60, (d) R70, (e) R80 and (f) R90. The dashed line indicates the cube location. Normalization is by the friction velocity 
$u_\tau$ in (1.1).
Figure 13(a) shows a sketch of the principal vortical structures that arise in a cube-roughened boundary layer flow (Martinuzzi & Tropea Reference Martinuzzi and Tropea1993; Hussein & Martinuzzi Reference Hussein and Martinuzzi1996), and figure 13(b) show the visualization of these principal vortical structures in DNS via the 
$Q$-criterion at the surface coverage density of 
$\lambda =0.5$. As the 
$Q$-criterion only gives a qualitative description of the vortical structures in the flow field, we more closely examine the induced flow due to the principal vortical structures in figure 14. The tip vortices give rise to LMPs and HMPs. The arch vortex is confined between two neighbouring cubes and gives rise to vortical motions as seen in the DNS results in figure 14. Specifically, figure 14(a,c) shows the in-plane streamlines and the in-plane velocity magnitude contours at a constant 
$y$ location through the centre of a wall-mounted cube in R50 and R70; figure 14(b,d) shows the in-plane streamlines and the in-plane velocity magnitude contours at a constant 
$z=0.8h$ location in R50 and R70. The head of the arch vortex gives rise to the vortex close to 
$z/h=1$ in figure 14(a,c), and the two legs of the arch vortex give rise to the vortices close to 
$y=\pm 0.5h$ in figure 14(b,d) The vortices are also confined by the back surface of the upstream cube and the front surface of the downstream cube, and their sizes decrease as the cubes are packed more closely. In addition, the roughnesses are obviously d-type, i.e. the flow above the cubes does not actively interact with the flow in the roughness-occupied layer.
Figure 13. (a) A sketch of the vortical structures when the roughness elements are closely packed. (b) Visualization of the vortical structures in R50 via the 
$Q$-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988).
Figure 14. (a) Contours of the mean streamwise velocity and the in-plane streamlines at a constant 
$y$ location through the centre of a row of wall-mounted cubes in R50. (b) Contours of the mean streamwise velocity and the in-plane streamlines at a constant 
$z$ location 
$z=0.8h$ in R50. (c) Same as (a) but for R70. (d) Same as (b) but for R70. Normalization is by the friction velocity 
$u_\tau$ in (1.1).
These secondary motions lead to a redistribution of the momentum/momentum flux in the RSL. Figure 15(a–d) shows contours of 
$\overline {u'w'}$ at a constant 
$y$ location through the centre location between two rows of cubes, and figure 15(e–h) show 
$\overline {u'w'}$ at a constant 
$y$ location through the centre of a row of cubes. R80 is similar to R90 and is not included for brevity. R60 is similar to R70 and is not shown for brevity as well. Form drag leads to momentum extraction at the front surfaces of the wall-mounted cubes, which in turn results in a large momentum flux from the bulk to the roughness just upstream of the cubes in R25 and R50 (figure 15e,f), as expected. Interestingly, if we compare figure 15(b,f) or figure 15(c,g), we see that, in R50 and R70, turbulence is actively transporting momentum to the thin slit between two rows of cubes, which, in turn, leads to a large friction force on the side surfaces as shown in figure 6 (b).
Figure 15. (a–d) Contours of 
$\overline {u'w'}$ through the centre location of two columns of cubes in (a) R25, (b) R50, (c) R70 and (d) R90. (e–h) Contours of 
$\overline {u'w'}$ through the centre of a column of cubes in (e) R25, (f) R50, (g) R70 and (h) R90. The normalization is by the nominal friction velocity 
$u_{\tau,N}$. The solid lines go through 
$x/h=1$, 
$z/h=2$ and 
$x/h=1$, 
$z/h=1.2$. The contour levels are kept unchanged.
4.2. Budget
Next, we examine the physical processes responsible for the generation and the destruction of the secondary motions by analysing the budget of 
$\langle \bar {u}''_i\bar {u}_i''\rangle$ and 
$\langle u'_iu_i'\rangle$. In this work, turbulent secondary motions refers to LMPs and HMPs, which manifests as 
$\bar {u}''$, so studying the generation and destruction of 
$\bar {u}''$ is to study the fate of secondary turbulent motions. This is like 
$u'$ and the TKE budget. The TKE budget governs the generation and the destruction of 
$u'$, and we study the TKE budget to understand the fate of 
$u'$.
In general, the total kinetic energy is (note that the prefactor 1/2 is not included here in the definition)
\begin{equation} \left\langle \overline{u_i u_i}\right\rangle=\left\langle \bar{u}_i\bar{u}_i\right\rangle+\left\langle \overline{u'_iu'_i}\right\rangle+\left\langle \bar{u}''_i\bar{u}''_i\right\rangle, \end{equation}where the first term on the right-hand side represents the kinetic energy in the mean flow (MKE), the second term represents the kinetic energy in the fluctuating turbulence (temporal variation of the fluid velocity) and the third term is the kinetic energy in the secondary motions (spatial variation of the time-averaged velocity). In anticipation of the results in the next section, we write the budget for the TKE in a periodic channel:
$$\begin{gather} 0={-}\left\langle \overline{u_1^\prime u_3^\prime}\right\rangle\frac{\partial\left\langle \bar{u}_1\right\rangle}{\partial {x_3}} -\left\langle \overline{u_i^\prime u_j^\prime}^{\prime \prime} \frac{\partial \overline{u_i}^{\prime \prime}}{\partial x_j}\right\rangle -\frac{\partial}{\partial { x_3}}\left\langle \frac{\overline{{u_3^\prime} {u'_iu'_i}}}{2} \right\rangle -\frac{\partial}{\partial { x_3}}\left\langle \frac{\overline{{u_3}}^{\prime \prime}\overline{u'_iu'_i}^{\prime \prime}}{2}\right\rangle\nonumber\\ -\,\frac{\partial}{\partial {x_3}}\left\langle {u_3^\prime} p^\prime \right\rangle +\nu\frac{\partial^2}{\partial {x_3^2}}\left\langle \frac{\overline{u'_iu'_i}}{2}\right\rangle -\nu \left\langle \overline{\frac{\partial u_i^\prime}{\partial x_j}\frac{\partial u_i^\prime}{\partial x_j}} \right\rangle, \end{gather}$$
where we have neglected buoyancy and assumed periodicity in both the streamwise and the spanwise directions. The flow is at a statistically stationary state, hence the 0 on the left-hand side. The terms on the right-hand side are the production terms (the first term, 
$P_t$, is the turbulent production, and the second term, 
$P_d$, is the wake production), the transport terms (the third term, 
$\varPi _t$, is the turbulent transport and the fourth term, 
$\varPi _d$, is the wake transport), the pressure work (the fifth term 
$W$), the viscous diffusion (the sixth term 
$D$), and the dissipation term (the last term 
$\epsilon$). In addition, the budget of the DKE reads (Raupach & Shaw Reference Raupach and Shaw1982)
$$\begin{gather} 0={-}\left\langle {\bar{u}_1'' \bar{u}_3''}\right\rangle\frac{\partial\left\langle {\bar{u}_1}\right\rangle}{\partial {x_3}} +\left\langle \overline{u_i^\prime u_j^\prime}^{\prime \prime} \frac{\partial \overline{u_i}^{\prime \prime}}{\partial x_j}\right\rangle -\frac{\partial}{\partial {x_3}}\left\langle \frac{{\bar{u}_3''} \bar{u}''_i\bar{u}''_i}{2} \right\rangle -\frac{\partial}{\partial { x_3}}\left\langle {\bar{u}_i''\overline{u'_i{u_3'}}^{\prime \prime}}\right\rangle \nonumber\\ -\frac{1}{\rho}\frac{\partial}{\partial { x_3}}\left\langle {\bar{u}_3''} \bar{p}'' \right\rangle +\nu\frac{\partial^2}{\partial { x_3^2}}\left\langle \frac{\overline{u''_i} \overline{u''_i}}{2}\right\rangle -\nu \left\langle \frac{\partial \overline{u_i''}}{\partial x_j}\frac{\partial \overline{u_i''}}{\partial x_j} \right\rangle +\frac{1}{\rho}\left\langle \bar{u}_i\right\rangle\left\langle \frac{\partial \bar{p}''}{\partial x_i}\right\rangle. \end{gather}$$The left-hand side is, again, 0. The first term on the right-hand side is the DKE-specific production term. The second term is the wake production term. This term shows up in the TKE budget as well, and it transfers energy from/to turbulence eddies to/from secondary eddies. The third, fourth and fifth terms are DKE-specific transport terms. The sixth term is the DKE-specific viscous diffusion term. The seventh term is the DKE-specific viscous dissipation term. The last term is a source term and in a periodic channel is zero.
Figure 16 shows the dominant terms in (4.2). All terms are very small in the roughness-occupied layer and we show results for 
$z>h$ only. The profiles collapse above approximately 
$(z-h)^+=60$, again suggesting thin RSLs in our DNSs. Figure 16(a) shows the production terms, i.e. the first and the second terms on the right-hand side of (4.2). The wake production is negative above the roughness occupied layer. The magnitude of the wake production increases as a function of 
$\lambda$ (from R25 to R50), stays constant (from R50 to R60) and then decreases to zero (from R60 to R90). Figure 16(b) shows the viscous dissipation. Its magnitude is similar to the production term. Figure 16(c) shows the viscous diffusion. This term transports energy from the upper part of the RSL to the lower part of the RSL. The transport terms and the pressure term are small and are not shown here for brevity. Figure 17 shows the significant terms in (4.3), i.e. the DKE-specific production term, the DKE-specific dissipation term, the DKE-specific transport term and the DKE-specific diffusion term. All the terms peak at 
$z=h$, are decreasing functions of 
$\lambda$ and vanish as 
$\lambda$ approaches one (at 
$100\,\%$ surface coverage density the dispersive stresses are zero).
Figure 16. (a) The production term. The bold lines represent turbulent production. The thin lines represent dispersive production. (b) The dissipation term. (c) The viscous diffusion term. The friction velocity 
$u_{\tau,T}$ in (1.5) is used for normalization.
Figure 17. (a) Production term, i.e. the first term on the right-hand side of (4.3). (b) The dissipation term, i.e. the seventh term on the right-hand side of (4.3). (c) The transport terms, i.e. the sum of the third and the fourth terms on the right-hand side of (4.3). (d) The diffusion term, i.e. the sixth term on the right-hand side of (4.3). Normalization is by the friction velocity 
$u_{\tau,N}$ in (1.3).
The data in figures 16 and 17 leads us to the following conclusions. TKE receives energy from MKE in both the lower and the upper part of the RSL (figure 16(a), bold lines). Viscous diffusion redistributes TKE in the RSL, transferring energy from the upper part of the RSL to the lower part of the RSL (figure 16d). In the lower part of the RSL, the energy from MKE and viscous diffusion are transferred to viscous dissipation (figure 16d). Wake production (figure 16(a), thin lines) transfers energy from TKE to DKE and is the dominant energy flux to DKE. DKE also gets energy from MKE (figure 17a) and viscous diffusion (figure 17d). However, both terms are smaller than the wake production term. The energy influx is balanced by the DKE-specific dissipation (figure 17b) and the DKE-specific transport. Because small-scale turbulence does not contribute significantly to the two most dominant terms in the DKE-budget equation, i.e. the wake production term and the DKE-specific dissipation term (note, the DKE specific dissipation term involves only the mean flow), the results here suggest that the accurate prediction of secondary turbulent flows’ presence does not critically depend on accurate modelling of the small-scale turbulence. This is consistent with, e.g., Anderson et al. (Reference Anderson, Barros, Christensen and Awasthi2015), where large-eddy simulation (LES) results agree well with the experiments although the small scales are modelled in LES.
5. Concluding remarks
In this work, we report DNS results of flow over aligned cubes. The wall-mounted cubes are closely packed with surface coverage densities of 0.25, 0.5, 0.6, 0.7, 0.8, 0.9 and 1. The nominal friction Reynolds number is kept constant at 
$Re_{\tau,N}=500$. The flow is in the fully rough regime at the surface coverage 0.25; at the surface coverage density 1, the top surfaces of the wall-mounted cubes form a new wall, and the flow reduces to plane channel flow. As few have studied roughness at such high surface coverage densities, this work fills a gap in the literature.
We have reported the essential flow quantities including the mean velocity profiles, the Reynolds stresses, the dispersive stresses and the roughness properties. The dispersive stress decays to zero at approximately 
$z=1.5h$,suggesting a thin RSL in all the DNSs. The equivalent roughness height 
$z_o$ is a decreasing function of the surface coverage density and the zero-plane displacement height 
$d$ is an increasing function of 
$\lambda$. We show that the method of Jackson (Reference Jackson1981) for determining 
$d$, which works well for cubes with low-to-moderate surface coverage densities, leads to underestimation of the zero-plane displacement for closely packed cubes.
Special attention has been given to secondary turbulent motions in the RSL. The spanwise alternating pattern of slits and cube surfaces gives rise to LMPs and HMPs above the cube crests. The strengths and the spanwise positions of these LMPs and HMPs depend on the surface coverage density for roughness within the d-type roughness regime. LMPs and HMPs are found in R25, R60 and R80 only. In R25, LMPs lie above cubes, whereas in R60 and R70, LMPs lie above the thin slits between two cubes. This result calls for more detailed roughness categorization (than the simple k-type and d-type categorization). More specifically, more data is needed before we can start to understand the positioning of LMPs and HMPs in the d-type roughness regime. To determine the physical processes responsible for the generation and destruction of the secondary flows, we analyse the DKE budget. The analysis shows that the secondary motions get energy from the mean flow and the fluctuating turbulence, and lose energy to DKE-specific dissipation.
Acknowledgements
Inputs from Dr. Stephen Lynch are gratefully acknowledged.
Funding
X.Y. and R.K. thank ONR for financial support, contract N000142012315, with Dr. Peter Chang as Technical Monitor.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Further discussion
We discuss the effects of domain size and grid resolution.
A.1. Domain size
In general, the size of the computational domain would have no/little effect on the flow statistics as long as the velocity correlation drops to zero within the domain. Figure 18 shows the correlation of the streamwise velocity fluctuation 
$u'$ in the 
$x$ and 
$y$ directions. Although the spanwise velocity correlation is zero at a distance of approximately 
$r_y=\delta$, the streamwise velocity correlation 
$\langle u'(x+r,y,z)u'(x,y,z)\rangle$ is non-zero in our R cases as well as the L80 case and would remain non-zero at the distance 
$r=20\delta$ (Sillero, Jiménez & Moser Reference Sillero, Jiménez and Moser2014), where 
$\delta$ is the half channel height. This puts a strict requirement on the domain size.
Figure 18. Correlation of the streamwise velocity fluctuation 
$u'$ in the streamwise and the spanwise directions in R25, R60 and R/L80 at 
$z=1.2h$. The bold lines are the correlations in the spanwise direction, and the thin lines are the correlations in the streamwise direction. The thin black solid line is at 0. The undulations in the correlations are due to turbulent dispersion (Borgas, Flesch & Sawford Reference Borgas, Flesch and Sawford1997). The reader is directed to Jimenez (Reference Jimenez1983) and Kitoh & Umeki (Reference Kitoh and Umeki2008) for more details.
Lozano-Durán & Jiménez (Reference Lozano-Durán and Jiménez2014) discussed the effect of domain size and concluded that the domain size has no effect on the first- and second-order statistics if 
$L_x>6 \delta$, 
$L_y>3\delta$. This domain size, i.e. 
$L_x=6\delta$, 
$L_z=3\delta$, is often referred to as the ‘minimal (plane) channel’. The effect of domain size on flow statistics in a rough wall TBL was discussed in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006), where the authors concluded that the computational domain has no effect on the mean flow as long as the domain is larger than 
$L_x\times L_y=4h\times 4h$. Our domain size of the R25-90 is 
$L_x\times L_y>8h\times 8h$ (in terms of half channel height, 
$L_x\times L_y>1.67\delta \times 1.67\delta$, 
$\delta =4h-h=3h$ is the half channel height), which is more than twice the suggested domain size in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006), and should be sufficient.
Nonetheless, to re-confirm the adequacy of our domain, we compare R80 and L80. The streamwise domain size of L80 is 
$L_x=19h$ (and 
$6.3\delta$, in terms of half channel height), which is more than twice the streamwise domain size of R80. We have kept the spanwise domain size the same between R80 and L80 considering that the velocity correlation drops to zero in the spanwise direction. In figures 19(a) and 19(b), we compare the mean velocity and the Reynolds stresses in R80 and L80, and there is barely any difference between the two DNSs. Hence, we re-confirm the conclusions in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006). It is worth noting that the fact that we are using a larger computational domain than Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) does not mean that the minimum (rough wall) channel suggested in Coceal et al. (Reference Coceal, Thomas, Castro and Belcher2006) is not adequate. The discussion here only shows that our domain has no effect on the resolved part of the flow.
Figure 19. (a) Time and horizontally averaged streamwise velocity in R80 and L80. Normalization is by 
$u_{\tau,T}$. (b) TKE in R80 and L80. Normalization is by 
$u_{\tau }$.
In addition to the mean velocity profiles and the TKE, we also examine how domain size affects the secondary flow structures. Figure 20(a) shows the streamwise velocity contour at a constant 
$x$ plane through the centre location of a wall-mounted cube for L80, and figures 20(b) and 20(c) show the spatial variation of mean streamwise velocity contour at 
$z=1.2h$, and there is barely any difference between L80 and R80 (see also figure 11e).
Figure 20. (a) Contours of the spatial variation of the streamwise velocity on a constant 
$x$ plane through the centre of a wall-mounted cube for L80. The contour line that go through 
$y=0$, 
$z=1.2h$ are highlighted. (b,c) Contours of spatial variation of the mean streamwise velocity at wall normal height 
$z=1.2h$, i.e. 
$0.2h$ above the cubes for (b) L80 and (c) R80.
A.2. Grid resolution
To reconfirm the adequacy/necessity of our somewhat excessive grid resolution, we compare two coarse-grid DNSs, i.e. C50 and C80, with their fine-grid counterparts, i.e. R50 and R80.
Figures 21(a) and 21(b) show the mean velocity profiles in R/C50 and R/C80. Figures 22(a) and 22(b) show the Reynolds stresses. The coarse-grid calculations yield practically the same mean flows and the Reynolds stresses, except for 
$R_{11}$ in R/C80. The coarse grid calculation C80 leads to a slightly larger 
$R_{11}$ than its fine-grid counterpart R80. This difference could be attributed to the less well-resolved small scales and the resulting less well-resolved turbulent dissipation. In all, our grid resolution, which is approximately twice as fine than the typical DNS resolution of DNS, helps us better resolve the small scales and produce more accurate results.
Figure 21. (a) Mean flow in R50 and C50. The profile above the roughness layer is shown. Normalization is by the friction velocity 
$u_{\tau,T}$ in (1.5). (b) Same as (a) but for R80 and C80.
Figure 22. Reynolds stresses in (a) R/C50 and (b) R/C80. Normalization is by the friction velocity 
$u_\tau$. Symbols are for C50/80 and lines are for R50/80.
Open access
