Flow over closely packed cubical roughness

Abstract Cube arrays are one of the most extensively studied types of surface roughness, and there has been much research on cubical roughness with low-to-moderate surface coverage densities. In order to help populate the literature of flow over cube arrays with high surface coverage densities, we conduct direct numerical simulations (DNSs) of flow over aligned cube arrays with coverage densities $\lambda =0.25$ (for validation and comparison purposes), $0.5$, $0.6$, $0.7$, $0.8$ and $0.9$. The roughness are in the d-type roughness regime. Essential flow quantities, including the mean velocity profiles, Reynolds stresses, dispersive stresses and roughness properties, are reported. Special attention is given to secondary turbulent motions in the roughness sublayer. The spanwise-alternating pattern of the thin slots between two neighbouring cubes gives rise to spanwise-alternating regions of low- and high-momentum pathways above the cube crests. We show that the strength and spanwise location of these low- and high-momentum pathways depend on the surface coverage density, and that the high-momentum pathways are not necessarily located directly above the roughness elements. In order to determine the physical processes responsible for the generation and the destruction of these secondary turbulent motions, we analyse the dispersive kinetic energy (DKE) budget. The data shows that the secondary motions get their energy from the DKE-specific production term and the wake production term, and lose energy to the DKE-specific dissipation term.


Introduction
Closely packed roughness is common in urban boundary layer flows (Barlow & Coceal 2008). A prominent example is closely packed buildings in a metropolitan area (Giometto et al. 2017;Krayenhoff et al. 2020). 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. 2007;Leonardi & Castro 2010;Yang 2016;Li & Bou-Zeid 2019). Depending on the packing density, wall-mounted cubes can † Email address for correspondence: xzy48@psu.edu be generally categorized into: isolated roughness, k-type roughness and d-type roughness (Oke 1988;Jiménez 2004). In the isolated regime, the roughness elements are sparsely packed and each roughness element acts as an isolated roughness element (see Yang et al. (2019) 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 λ. 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 (λ > 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 λ = 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 2008, 2005Huang et al. 2000).
In the following, we summarize a few relevant studies on flow over cubical roughness. Cheng & Porté-Agel (2015) and Yang (2016) 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 δ RSL < 2h (Cheng et al. 2007;Leonardi & Castro 2010), where δ 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. (2015) 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. 2006;Jiang et al. 2008). In particular, Coceal et al. (2006) studied the effects of domain size and concluded that a domain of size 4h × 4h × 4h is sufficient for capturing the mean flow. Blackman, Perret & Calmet (2016) and Blackman et al. (2017) 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) . LMPs and HMPs have received significant attention in the recent literature. For example, Mejia-Alvarez & Christensen (2013), Nugroho, Hutchins & Monty (2013) and  established that any spanwise heterogeneous roughness can give raise to LMPs and HMPs in the mean flow, and Anderson et al. (2015) showed that the HMPs and LMPs observed in previous works (Mejia-Alvarez & Christensen 2013; Willingham et al. 2014) 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 λ = 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 1998;Leonardi & Castro 2010;Lee, Sung & Krogstad 2011;Volino, Schultz & Flack 2011;Yang & Meneveau 2016, 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 λ → 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 2017). 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: 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 ρ = unity is the fluid density. It follows that the definition of the friction Reynolds number is A second definition leads to the so-called nominal friction velocity A third definition is by Coceal et al. (2006). 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.
(1.4) Chan-Braun, García-Villalba & Uhlmann (2011) and Forooghi et al. (2018) 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 τ . Lastly, we define 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 τ . Having defined the friction velocity, the mean flow in the inertial layer conforms to where U + = U/u τ , z + = zu τ /ν, κ is the von Kármán constant, C is a constant and 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 (1.7) Equating (1.6) and (1.7), we have and U + ≈1/κ log z + 0 exp(κC) . (1.9) That is, 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 − λ) 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 2018), 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 λ > 0.4, the roughness fields are d-type and behave qualitatively similarly as aligned cubes with a surface coverage density of λ = 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.  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

Computational details
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 1987;Hoyas & Jiménez 2006;Graham et al. 2016;Yamamoto & Tsuji 2018). Details of our DNS simulations are summarized in table 1. We vary λ from 0.25 to 1. For R100, i.e. for λ = 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. (2006). For reference purposes, the computational domain size in Coceal et al. (2006) is Leonardi & Castro (2010). 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 + n < + ν,l at the wall and z + ≈ 5 + ν in the bulk. Here, + n is the wall-normal grid spacing, + ν,l is the locally defined viscous unit, + x/y/z is the grid spacing in the x, y and z directions, and + ν = ν/u τ is the viscous unit. Table 2 lists the detailed grid information. For reference purposes, the streamwise and spanwise grid resolution is + x/y = 18.75 in Leonardi & Castro (2010), + x/y/z = 7.8 to 15.6 in Coceal et al. (2006) and + z ≈ 8 to 13 in the recent work by MacDonald et al. (2018). The resolution of a DNS 920 A37-5 Mesh' is the total number of grids in million cells (M) and 'N' is the number of wall-mounted cubes in n cube,x × n cube,y , where n cube,x and n cube,y are the numbers of cubes in x and y directions, respectively. x/h should be proportional to the Kolmogorov length scale: where 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η 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. 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 2004). This code has been extensively used for wall-bounded flow calculations (Ma, Yang & Ihme 2018;Yang et al. 2019). Further details of the code can be found in Khalighi et al. (2011) and Bermejo-Moreno et al. (2014). 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 − u w + , the dispersive stress − ū w + and the viscous stress dū + /dz + to 1 − z/L z , where u, v and w are the velocity in the three Cartesian directions, denotes temporal fluctuation, · denotes horizontal averaging,· 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).

Basic flow phenomenology
In this section, we report the most essential flow quantities including the mean velocity profiles ū , 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 λ = 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 ū , R ij , D ij , z o and d are monotonic functions of λ 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) + ≈ 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 (Cionco 1965;).

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 τ S /τ R in our DNSs, where τ S is the friction force on the bottom wall and τ 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 τ S /τ 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 λ = 0.25. For λ > 0.6, the top surface friction dominates, and the flow transits from being fully rough to transitionally rough.
Next, we measure the effective roughness height z 0 and the displacement height d. We use a von Kármán constant of κ = 0.384 following Marusic et al. (2013). The friction velocity is u τ,C following Coceal et al. (2006). Jackson (1981) argued that the virtual ground must be located at the height where the resultant drag force acts (3.1) 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 (1981) 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.  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. (2007) of aligned cubes at a surface coverage density of 0.25, which gives a rough idea of the uncertainty in the data. 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 λ increases below z = h. Above the 920 A37-10 cube-occupied layer, the peak magnitude in R 11 increases as λ increases. The R 11 profiles do not depend significantly on the surface coverage density above z/h ≈ 1.6, suggesting a thin RSL. The dispersive component D 11 peaks at z ≈ h, and decreases with λ in the roughness-occupied layer. Above z = h, the dispersive stress decays rapidly to zero, again, suggesting a thin RSL.

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. 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.ū + =ū + − ū + , 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 2015). 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ū 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ū . Figure 13(a) shows a sketch of the principal vortical structures that arise in a cube-roughened boundary layer flow (Martinuzzi & Tropea 1993;Hussein & Martinuzzi 1996), and figure 13(b) show the visualization of these principal vortical structures in DNS via the Q-criterion at the surface coverage density of λ = 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 = ±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.

Secondary flow phenomenology
These secondary motions lead to a redistribution of the momentum/momentum flux in the RSL. Figure 15(a-d) shows contours of u w at a constant y location through the centre location between two rows of cubes, and figure 15(e-h) show 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).

Budget
Next, we examine the physical processes responsible for the generation and the destruction of the secondary motions by analysing the budget of ū iū i and u i u i . In this work, turbulent secondary motions refers to LMPs and HMPs, which manifests asū , so studying the generation and destruction ofū 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) 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: 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, Π t , is the turbulent transport and the fourth term, Π 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 ). In addition, the budget of the DKE reads (Raupach & Shaw 1982) 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 λ (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 λ and vanish as λ approaches one (at 100 % surface coverage density the dispersive stresses are zero).
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. (2015), where large-eddy simulation (LES) results agree well with the experiments although the small scales are modelled in LES.

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 τ,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 λ. We show that the method of Jackson (1981) 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.   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 1997). The reader is directed to Jimenez (1983) and Kitoh & Umeki (2008)

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 = δ, the streamwise velocity correlation u (x + r, y, z)u (x, y, z) is non-zero in our R cases as well as the L80 case and would remain non-zero at the distance r = 20δ (Sillero, Jiménez & Moser 2014), where δ is the half channel height. This puts a strict requirement on the domain size.
Lozano-Durán & Jiménez (2014) 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δ, L y > 3δ. This domain size, i.e. L x = 6δ, L z = 3δ, 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. (2006), 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 × L y = 4h × 4h. Our domain size of the R25-90 is L x × L y > 8h × 8h (in terms of half channel height, L x × L y > 1.67δ × 1.67δ, δ = 4h − h = 3h is the half channel height), which is more than twice the suggested domain size in Coceal et al. (2006), 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δ, 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. (2006). It is worth noting that the fact that we are using a larger computational domain than Coceal et al. (2006) does not mean that the minimum (rough wall) channel suggested in Coceal et al. (2006) is not adequate. The discussion here only shows that our domain has no effect on the resolved part of the flow. 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).

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. 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.