Effects of Wall Topology on Statistics of Cube-Roughened Wall Turbulence

In this work, we carried out roughness-resolved direct numerical simulations of cube-roughened turbulent channel flows to investigate the effects of aligned and staggered arrangements, element spacings, and element orientations on the statistics of rough-wall turbulence. The results show that the equivalent sandgrain roughness, Reynolds stresses and dispersive stresses are affected by element spacings and arrangements in different ways depending on the cube orientations. Placing the roughness elements in a staggered way in general increases the equivalent sandgrain roughness, unless when the element–element interaction is insignificant for the cube orientation with wakes of short length. As for the Reynolds normal stresses and the streamwise component of the dispersive stresses (both are normalized by the total friction velocity), placing the cube elements in a staggered way decreases their maximal values when compared with the aligned arrangements. As for the effects of element spacing on the equivalent sandgrain roughness and the maximums of the Reynolds normal stresses, similar trends are observed for different element orientations and arrangements when increasing the element spacing l/r (where r is the cube height) from 2.0 to 2.8. When further increasing the element spacing l/r from 2.8 to 3.5, however, different trends are observed for different element orientations. As for the dispersive stresses, greater maximal values of their streamwise components are observed for larger element spacing for all the considered cube-roughened surfaces. For the turbulence statistics in the wake of the cube element, certain similarities are observed for different element spacings.


Introduction
Rough-wall turbulence ubiquitously exists in nature and engineering applications (Jiménez 2004), such as atmospheric boundary layer flows over urban and vegetation sites (Giometto et al. 2016), river flows over gravel beds (Mignot et al. 2009), ocean currents over the hulls covered by barnacles and algae (Mejia-Alvarez and Christensen 2010), and the flow over a propeller eroded by cavitation (Sezen et al. 2021). One important issue, which remains unresolved, is how the drag on the surface is related with the roughness topography (Flack and Schultz 2010;Jouybari et al. 2021). This work is devoted to investigate the effects of wall topology on drag and turbulence statistics for cube-roughened surfaces using roughnessresolved direct numerical simulations.
Research on rough-wall turbulence dates back to the pioneering work by Colebrook (Colebrook 1939), Nikuradse (Nikuradse 1950) and Moody (Moody 1944). Early researchers carried out such studies mainly based on laboratory experiments. With the exponential growth of the computing power of super computers, simulations are playing a more and more important role in the study of rough-wall turbulence. Rough walls considered in the literature can be classified into two types, rough walls formed by organized roughness elements, e.g., sinusoidal roughness Ma et al. 2020), cubical roughness (Coceal et al. 2007;Leonardi and Castro 2010), riblets (Modesti et al. 2021), hemispherical roughness (Chatzikyriakou et al. 2015), and rough walls formed by roughness elements with certain degree of randomness, e.g., uniform ellipsoidal roughness randomly rotated (Yuan and Piomelli 2014a, b;Yuan and Jouybari 2018), ellipsoidal roughness of different size, aspect ratio and inclination (Jouybari et al. 2021), Fourier transform of white noise Jouybari et al. 2021), and real-life rough surface (Yuan and Jouybari 2018).
Cubes are one of the most widely employed roughness elements in the literature (Castro et al. 2006;Lee et al. 2011;Djenidi et al. 1994). Studies have been carried out on how the topology of cube arrays affects the roughness length and turbulence statistics, e.g., effects of element spacing (Cheng et al. 2007;Ahn et al. 2013;Yang et al. 2019;Xu et al. 2021), effects of element height, aspect ratio and arrangement (Yang 2016;Hamed et al. 2017;Sadique et al. 2017;Choi et al. 2020). Effects of the frontal and plane solidities on the drag of surfaces formed by LEGO bricks were investigated experimentally by Placidi and Ganapathisubramani (2015). For a fixed frontal solidity λ F = 0.15, they found that the drag monotonically decreases when the plane solidity is changed from 0.11 to 0.44. Using direct numerical simulations of staggered cube arrays, Lee et al. (2012) observed that the roughness length achieves its maximum value at streamwise element spacing of l/r = 4.
Dynamics of the coherent structures in the logarithmic layer and their interactions with the flow structures at the roughness element scale were investigated in the literature (Kanda et al. 2004;Kanda 2006;Volino et al. 2011;Basley et al. 2018). Inagaki and Kanda (2010) investigated the coherent structures in the logarithmic layer over an aligned cube array using outdoor scale model experiments, in which the height of the cube is 1.5 m. Perret and Kerhervé (2019) observed large-scale elongated coherent structures of low or high momentum with length of several boundary layer thickness from the analysis of wind tunnel experimental data using the spectral proper orthogonal decomposition method. Coceal et al. (2007) found that the characteristic scale of the coherent structures increases linearly with the distance from the wall in the logarithmic region. Basley et al. (2019) showed that the coherent structures in the inertial layer are independent of the configuration of the rough wall for a boundary layer over a staggered cube array for three different values of plane solidity (which is defined as the ratio of the area occupied by the roughness to the total area), i.e., 6.25%, 25%, 44.4%.
From the wind tunnel experiments of a turbulent boundary layer over a staggered cube array, Ferreira and Ganapathisubramani (2021) found that the uncorrelated, intermediate and smallscale pressure events are important to the overall drag fluctuations, which are modulated by large-scale structures in the outer layer. The interaction between the large-scale momentum regions and small-scale structures due to roughness elements was studied by Blackman and Perret (2016) using stochastic estimation for velocity decomposition, in which the data were from the experiments of cube-roughened boundary layers. Using large-eddy simulation, Anderson (2016) showed that the outer layer dynamics can modulate the streamwise velocity fluctuations in the roughness sublayer. The modulation effect of the large-scale motion on the small-scale turbulence in the roughness sublayer was also investigated in wind tunnel experiments of a boundary layer over staggered cube arrays by Basley et al. (2018).
The existence of outer layer similarity is one major focus in the study of rough-wall turbulence. Using direct numerical simulation (DNS) of turbulent channel flows with two-and three-dimensional roughness, Orlandi and Leonardi (2006) observed a rather good collapse of the roughness function plotted against the rms (root-mean-square) of the normal velocity on the plane of the crests. Measurements at field-scale of an aligned cube array were carried out by Roth et al. (2015), which showed that the spatial inhomogeneity almost disappears at 1.5r above the ground (where r is the cube height). Wind tunnel measurements by Amir and Castro (2011) indicated that for the ratio of roughness height to boundary layer thickness less than 0.15, the Reynolds stresses collapse well with each other, which is 0.2 for the mean velocity profile. The inner-layer scaling similarity was examined for an aligned cube array at outdoor scale by Inagaki and Kanda (2008), and was found being robust for the wallnormal fluctuations but not for the horizontal velocity fluctuations. From the measurements in a hydraulic flume of aligned cube arrays, Macdonald et al. (2002) found that the rms of velocity fluctuations can be scaled by the local values of the shear stress. In the study by Placidi and Ganapathisubramani (2018), the authors observed the lack of similarity for the streamwise and wall-normal Reynolds stresses and Reynolds shear stresses, and the criterion based on the "effective shelter area" (Raupach and Shaw 1982) can capture the departure from the outer layer similarity.
Engineering models for predicting the aerodynamic or hydrodynamic properties of cuberoughened surfaces have been developed in the literature. Lettau's relation (Lettau 1969) is one of the earliest models proposed in the literature. An improved method based upon Lettau's relation for estimating the roughness length was proposed by Macdonald et al. (1998). Different geometry parameters of the surface roughness, such as the mean roughness height, first-order moment of height fluctuations, skewness and kurtosis of the roughness height fluctuations, are employed to relate the surface topology with the drag of the surface . Models on predicting the drag of rough surfaces (often using the equivalent sandgrain roughness, a hydraulic scale defined by the drag) based on the roughness geometrical features were developed in the literature Flack and Schultz 2010). Theoretical model for predicting the roughness length for cube-roughened walls was developed by Jia et al. (1998). The drag law for different roughness geometries was derived by Wooding et al. (1973). An analytical model based on the von Kármán-Pohlhausen integral method was proposed by Yang et al. (2016) for predicting mean velocity and drag for flow over rectangular-prism roughness elements. In a recent work, Jouybari et al. (2021) developed data-driven models using deep neural network and Gaussian process regression, based on the data from high-fidelity simulation and experiments. In total, 45 different rough surfaces were employed for training the model in their study, with eight different geometry parameters of the surface roughness as the input features. Their results showed that the maximum error is less than 30%.
Although various engineering models for the sandgrain roughness k s exist in the literature, the prediction uncertainty is still high . Employing the big data and the machine learning method is promising for developing the models for k s . In order to reduce the amount of the required dataset and to improve the applicability of the data-driven models to a wide range of rough surfaces, deep understanding on the underlying flow physics is necessary. In nature and engineering applications, the arrangements (aligned or staggered) and the element orientation will be changed when the flow direction changes. To the best of the authors' knowledge, such effects have not been systematically investigated in the literature. In this work, we focus on the cube-roughened surfaces. The objective is to examine how element orientations and arrangements affect flow structures in the roughness sublayer and turbulence statistics for different element spacings. Specifically, we carry out roughnessresolved direct numerical simulations of turbulent channel flows with cube-roughened walls. For both the aligned and staggered arrangements, three different element orientations and element spacings are considered. As this work is focused on the drag of the rough surfaces and the turbulence statistics in the near surface region, we employ a minimal-span channel in the simulations (Jiménez and Moin 1991;Chung et al. 2015;MacDonald et al. 2017) to reduce the computational cost, which has been shown being able to accurately capture the mean streamwise velocity profile and the flow structures in the logarithmic region, while will be less accurate if the focus is on the large-scale coherent flow structures in the outer layer.
The rest of the paper is organized as follows. The employed numerical method is first described in Sect. 2, followed by the setup of the simulated cases in Sect. 3. The results from the simulations are then presented in Sect. 4. At last, conclusions are given in Sect. 5.

Numerical Method
The Virtual Flow Simulator (VFS-Wind) code (Yang et al. 2015b;Calderer et al. 2015;Yang et al. 2015a), which has been successfully applied to simulate flows over complex geometries Zhou et al. 2021), is employed in this work. A further development of the VFS-Wind code for simulating particle-laden turbulence in the presence of complex boundaries can be found in our recent work (Qin et al. 2022).
The governing equations are the three-dimensional, incompressible Navier-Stokes equations in non-orthogonal, generalized, curvilinear coordinates, which read in compact tensor notation (repeated indices imply summation) as follows (i, j, l = 1, 2 and 3 representing x, y and z directions, respectively): in which J is the Jacobian of the geometric transformation, U i = (ξ i m /J )u m is the contravariant volume flux, t is the time, x i and ξ i are the Cartesian and curvilinear coordinates, respectively, ξ i l = ∂ξ i /∂ x l are the transformation metrics, u i is the i th component of the velocity vector in Cartesian coordinates, g jk = ξ j l ξ k l are the components of the contravariant metric tensor, μ is the dynamic viscosity, ρ is the density, and p is the pressure.
The curvilinear immersed boundary (CURVIB) method (Ge and Sotiropoulos 2007) is employed for representing the roughness elements. In the CURVIB method, the flow is simulated based on a curvilinear (or Cartesian) grid, while the boundaries are discretized using triangular meshes independent of the background grid for flow simulations. To apply boundary conditions for the flow simulations, the background grid nodes are classified into fluid nodes and solid nodes. The fluid nodes, which have at least one neighbor in the solid, are further identified as the immersed boundary (IB) nodes. The boundary conditions for the flow simulations are then applied at the IB nodes by interpolating in the wall normal direction using the boundary conditions at the wall and the values of flow quantities at the fluid nodes. The governing equations are discretized in space using a second-order accurate central differencing scheme, and integrated in time using the fractional step method (Ge and Sotiropoulos 2007). An algebraic multigrid acceleration along with a GRMES solver is used to solve the pressure Poisson equation. A matrix-free Newton-Krylov method is used for solving the discretized momentum equations.

Case Setup
In this section, we present the computational setup of the simulated cases. Cubical roughness elements are employed to generate the rough surface topology. In total, 18 cube-roughened surfaces, with different element spacings l (the element spacing is the same in the streamwise (x) and spanwise (z) directions) and different ways of rotations and arrangement (aligned and staggered), are employed in the simulations as listed in Table 1. The height (width) of the cube is r = 0.0714h for all the simulated rough surfaces, where h is the height of the channel.  The basic geometry for generating the rough wall and the naming rules are described. In Table 1, the first letter of the case name stands for the shape of the roughness element, i.e., letter "C" for cube. The number after the first letter denotes the element spacing l, for instance, 20 in C20D0A means l = 2.0r . The letter "D" denotes the direction of the element relative to its original direction, with 0 for the original orientation of the roughness element. For instance, the letter-number combination "s45" stands for rotation along the spanwise axis by 45 • as shown in Fig. 1b, and "n45" stands for rotation along the vertical axis by 45 • as shown in Fig. 1c. The last letter "A" stands for aligned arrangement of the cubical roughness elements as shown in Fig. 1d, and letter "S" for staggered arrangement as shown in Fig. 1e. The so-generated rough surfaces for the simulated cases are shown in Fig. 2. The geometrical parameters of the simulated rough surfaces are listed in Table 2, which include the roughness crest k c = max(k), the mean roughness height k avg = 1 A t x,z k d A, the root mean square (rms) of the height fluctuations k rms = 1 A t x,z k − k avg 2 d A, the skewness of the roughness height fluctuations S k = 1 A t x,z k − k avg 3 d A/k 3 rms , the frontal solidity λ f = A f /A t , the plan solidity λ p = A p /A t . In the above expressions, k(x, z) is the height of rough surface at different positions, A t is the total planar area, and A f is the frontal area of all roughness elements, A p is the plan area of roughness elements.
The Reynolds number Re τ based on the friction velocity u τ = √ τ w /ρ (where τ w is the wall shear stress) and the height of channel (h) is approximately 1000. Slight differences are observed for different cases, with the actual Reynolds number shown in Table 1. The size of the computational domain is 3h ×h ×h in the streamwise x, vertical y and spanwise z directions, respectively. The first off-wall grid node is located at y + = 0.7 from the wall. In the other two directions, the grid spacings in wall units are approximately 6. The total number of grid On the bottom wall, the no-slip boundary condition is employed. Free-slip boundary condition is applied at the top boundary. Periodic boundary condition is applied in the horizontal directions. The flow is driven by a mean pressure gradient to maintain a constant mass flux. For all the cases, we first run the simulation to achieve a fully developed state and then average the flow for approximately 60 h/U b to compute the turbulence statistics.

Results
In this section, we present the results from the simulated cases, which include the instantaneous and time-averaged flow field, the temporally and horizontally averaged turbulence statistics, and the turbulence statistics in the wake of the cube element.

Instantaneous and Time-Averaged Flow Field
To have an intuitive picture on the flow over different rough surfaces, the contours of the instantaneous streamwise velocity field are examined in Fig. 3. As seen, the flow in the roughness sublayer is featured by patches of high-speed and low-speed regions. The locations of these patches are affected by the arrangement and orientation of the roughness element. For instance, the instantaneous flow fields are featured by long streaks of high-speed velocity located in-between rows of roughness elements through the whole channel for the cases with aligned roughness elements (e.g., C28D0A, C35D0A), which, on the other hand, are shorter and located in an oblique way between the roughness elements for staggered cases when l/r = 2.8 and l/r = 3.5. For the cases with larger spacings between roughness elements, the high-speed patches are larger with higher velocity magnitude (e.g., C35D0A compared with C28D0A). On the vertical slice, the interaction between the wake of the element and the flow above is observed, being stronger for larger element spacings for the simulated cases. Particularly, the wake flows upward when the roughness elements are rotated along the spanwise axis by 45 • . Vertically rotating the element, on the other hand, does not significantly affect the fluid motion in the vertical direction. The flow structures identified using the λ 2 criterion are then examined in Fig. 4. It is seen that the complexity of the flow structures is strongly related with the element orientations, arrangements and spacings. For larger element spacings, the flow structures become more Fig. 3 Contours of the instantaneous streamwise velocity on a vertical slice and a horizontal slice passing through the roughness elements for flows over different rough surfaces. The domain in the box is zoomed in to show the flow structure around roughness elements clearly complex especially for staggered cases, influencing a higher and wider region around the roughness element. For most cases, the flow structures around neighboring roughness elements are relatively independent from each other. For some cases, on the other hand, strong interactions between the flow structures of different roughness elements are observed, for instance, the staggered cases with the roughness element rotated along the vertical axis by 45 • . One important feature is the upward motion of the flow as indicated by the inclination of the vortex structure shown in Fig. 4, which involves in the momentum exchange between the wake region and the above high momentum region. The pattern of this upward motion depends on the element orientation as well as the element spacing. One observation is that such upward motion persists when the roughness elements are rotated along the spanwise axis by 45 • , which is consistent with that shown by the instantaneous streamwise velocity in Fig. 3, although the rotation essentially reduces the projected area facing the upstream. The upward motion is not strongly affected by the way of arrangements when the roughness elements are rotated along the spanwise axis by 45 • . For other element orientations, on the other hand, differences are observed between the aligned and staggered layouts for large element spacings. Overall, differences are observed among the flow patterns around the roughness element, although certain similarities are observed. This highlights the difficulty in modeling the effect of individual roughness elements on the upward flow motion in analytical models for rough-wall turbulence.
To demonstrate the three-dimensional geometrical feature of the roughness element wake, the iso-surfaces of time-averaged streamwise velocity with u = 0 are shown in Fig. 5 for different rough surfaces. It is seen that the element spacing is the most important factor affecting the overall shape of the element wakes. The type of element arrangement plays a key role on the interaction of wakes from neighboring rows, that the wakes are connected when the elements are placed in a staggered way, while are isolated for the aligned arrangements for l/r = 2.0, 2.8. For the element spacing l/r = 3.5, on the other hand, the wakes from neighboring rows are isolated for both aligned and staggered arrangements. In general, it is observed that the frontal area of the element facing the forward flowing fluid increases with the increase of element spacing, which is also affected by the element orientation as expected.

Temporally and Horizontally Averaged Turbulence Statistics
In the following, the turbulence statistics, including the mean streamwise velocity, the Reynolds stresses and the dispersive stresses, are examined quantitatively. The dispersive stresses describe the spatial heterogeneity of the time-averaged flow field, which is computed by decomposing the instantaneous velocity u i as follows: where · denotes the temporal averaging, and · s for spatial averaging in the horizontal directions. It is noted that the spatial averaging is carried out in both solid and fluid region in the roughness sublayer. The profiles of the mean streamwise velocity are shown in Fig. 6. A clear logarithmic region is observed for all the simulated cases. The effect of the rough wall topology on the roughness function, which is the difference between the simulated mean streamwise velocity and the logarithmic law for a smooth wall, is examined. The orientation of the roughness element affects the roughness function, but in a different way for different element spacings and element arrangements. Placing the roughness elements in a staggered way increases the Fig. 4 Near-wall time-averaged flow structures visualized using the λ 2 criterion (λ 2 = −0.5). The iso-surfaces are colored using the time-averaged streamwise velocity. Only a part of the domain is shown in the figure roughness function, which, however, is not clearly observed when the roughness element is rotated along the spanwise axis by 45 • (e.g., C28Ds45A vs. C28Ds45S). Increasing element spacing in general increases the roughness function for the staggered layouts, which is not clearly shown for the aligned layouts and the staggered layouts with elements rotated along the spanwise axis (i.e., C20Ds45S, C28Ds45S and C35Ds45S). After showing the mean streamwise velocity profiles, the variations of the sandgrain roughness k s and the zero-plane displacement d via the element spacing l are examined in Figs. 7 and 8, respectively. The zero-plane displacement d is meant to describe the position where the mean streamwise velocity is zero, which, however, might not be the case . In this work, the values of d and k s are computed by fitting the mean streamwise  Table 1 velocity using the logarithmic law, which is in the form of u + = 1 κ ln y−d k s + 8.5. When computing d and k s , the friction velocity u τ for normalization is computed using the mean pressure gradient in the streamwise direction, which is the driving force for the flow. In Fig. 7,  Fig. 7 Variations of sandgrain roughness k s via element spacing for different rough surfaces it is seen that the variations of k s for the cube-roughened surfaces are fairly complex. Two different kinds of trends are observed. For the surfaces with cubes vertically oriented and arranged in a staggered way (D0S) and the surfaces with the cubes rotated along the vertical axis by 45 • and arranged in an aligned way (Dn45A), k s gradually increases with l for l/r from 2.0 to 3.5, which increases for l/r from 2.0 to 2.8, and decreases for l/r from 2.8 to 3.5, respectively, for other cases. Among all the simulated cube-roughened surfaces, the values of k s from the surfaces with cubes rotated along the vertical axis by 45 • and arranged in a staggered way (Dn45S) are the highest for element spacing l/r = 2.8, 3.5. Arranging the roughness elements in a staggered way increases the values of k s when compared with the aligned arrangement for the D0A and D0S, and Dn45A and Dn45S cases. For the Ds45A and Ds45S cases, in which the roughness element is rotated along the spanwise axis by 45 • , the values of k s from the aligned and staggered arrangement, on the other hand, are close to each other for the considered element spacings. This can be explained by the relatively short wake for this cube orientation, that the cube wake barely affects the incoming velocity of the downstream cube for the considered element spacings, such that element spacing instead of element arrangement is the more dominant factor affecting k s for this cube orientation. The different trends of k s for l/r ranging from 2.0 to 2.8 and from 2.8 to 3.5 observed in most cases are associated with different flow regimes. For small element spacings (high density of roughness elements), the sheltering effect (Millward-Hopkins et al. 2011;Yang et al. 2016) dominates. The increase in k s with l/r is associated with the increase in the area of the surface facing the forward flowing fluid. For large element spacings, the area of the surface facing the forward flowing fluid barely changes when increasing l/r . The decrease in k s with l/r is caused by the decrease in the density of roughness elements. The flow regimes varying with the element spacing are classified as the skimming flow, the wake interference flow and the isolated roughness flow in the work by (Oke 1988).
The variations of zero-plane displacement d are shown in Fig. 8 for different rough surfaces. Complex variations of d are observed. For the surfaces with cubes non-rotated and arranged in an aligned way (D0A) and the surfaces with cubes rotated along the vertical axis by 45 • and arranged in a staggered way (Dn45S), the value of d decreases for l/r from 2.0 to 2.8, and increases for l/r from 2.8 to 3.5, respectively. For the surfaces with cubes rotated along the vertical axis by 45 • and arranged in an aligned way (Dn45A), on the other hand, the value of d increases for l/r from 2.0 to 2.8, and decreases for l/r from 2.8 to 3.5, respectively. For other surfaces, the value of d gradually decreases with l.
The Reynolds stresses from the cases with different rough surfaces are compared in Fig.  9. It is seen that the three components of Reynolds normal stresses from high to low ranked by their magnitudes are u 2 + s , w 2 + s and v 2 + s . The peaks of the maximal u 2 + s and w 2 + s are located at positions close to the roughness crest, while the peaks for the maximal v 2 + s are found at locations slightly further from the surface. The maximal values of v 2 + s and w 2 + s from different cases are fairly similar with each other for different rough surfaces. For the maximal values of u 2 + s , it is seen that they generally decrease when increasing the roughness element spacing (e.g., from C20D0A and C28D0A to C35D0A).
The maximal values of u 2 + s are higher for the aligned arrangement when compared with the staggered arrangement (e.g., C20D0A vs. C20D0S) for most cases, except for the one with the cubical element rotated along the vertical axis by 45 • for l/r = 2.0 (i.e., C20Dn45A vs. C20Dn45S). The maximal values of u 2 + s are decreased when rotating the roughness elements by 45 • along the spanwise axis for both aligned and staggered layouts for l/r = 2.0, which are barely changed for l/r = 2.8r , 3.5, except for the staggered layout with l/r = 3.5, in which the maximal value of u 2 + s from the surface with inclined roughness elements is higher (C35D0S vs. C35Ds45S). When rotating the cubical roughness elements by 45 • along the vertical axis, the maximal values of u 2 + s remain approximately the same (e.g., C28D0A vs. C28Dn45A), except for the aligned layout with l/r = 2.0 (i.e., C20D0A vs. C20Dn45A), in which the maximal value of u 2 + s from the case C20D0A is higher than that from C20Dn45A.
Below the roughness crest, one obvious difference between the aligned and staggered cases are the values of u 2 + s and w 2 + s . The values of u 2 + s are larger than that of w 2 + s for aligned cases, while the values of u 2 + s are close to w 2 + s for the staggered cases. The different patterns of the flow structures in the roughness sublayer are the major reason for the above differences. For the aligned cases, the high-speed streaks along with the streamwise direction located in the gap between two rows of roughness elements (as shown in Fig. 3) are Fig. 9 Vertical profiles of the normal Reynolds stresses computed from the cases with different rough surfaces. The parameters for each rough wall are shown in Table 1. The vertical grey dash lines show the location of the roughness crest responsible for the streamwise velocity fluctuations of high intensity. The high-speed streaks in the staggered cases, on the other hand, are short and placed in an oblique way, which are effective in inducing the velocity fluctuations in the spanwise and vertical directions.
The vertical profiles of the dispersive stresses ũ 2 + s are shown in Fig. 10. The maximal values of ũ 2 + s are approximately one order of magnitudes higher than the other two components. The dispersive stresses measure the spatial heterogeneity of the time-averaged velocity field. Significant differences are observed among cases with different rough surfaces.
It is seen the ũ 2 + s deceases in an abrupt way as moving upward from the top of the roughness elements. The magnitudes of ũ 2 + s are higher for the aligned layout when compared with the corresponding staggered layout (e.g., C20D0A vs. C20D0S). For the aligned layouts (C20D0A vs. C20Ds45A and C20Dn45A), rotating the roughness elements decreases the magnitudes of ũ 2 + s no matter vertically or along the spanwise axis. For the staggered layouts (C20D0S vs. C20Ds45S and C20Dn45S), on the other hand, rotating the cubical roughness element does not have significant effects on the magnitudes of ũ 2 + s . Increasing the element spacing increases the magnitudes of the ũ 2 + s for all the considered orientations and layouts of the cube-roughened surfaces. The other two components of the dispersive stresses, i.e., ṽ 2 + s and w 2 + s , are shown in Fig. 10. It is seen that the magnitudes of ṽ 2 + s and w 2 + s are very small when compared with ũ 2 + s for some rough surfaces, e.g., C20D0A. The magnitudes of the spanwise component w 2 + s are higher than the vertical component ṽ 2 + s . Magnitudes of the dispersive stresses are associated with the way how the roughness elements affect the flow in the roughness sublayer. The profiles of the dispersive stresses (Fig. 10) and the flow structures in the roughness sublayer (Fig. 3) show that the magnitude of the streamwise component of the dispersive stress is high when the high-speed and low-speed streaks are long in the streamwise direction, which is the case for the aligned arrangement and large streamwise element spacings. The magnitudes of the spanwise and vertical components of the dispersive stress, on the other hand, are high when the flows are disturbed in the corresponding directions, which happens for the staggered arrangement. The magnitudes of ṽ 2 + s and w 2 + s are significantly larger for the staggered layouts when compared with the aligned layouts, which is consistent with the high-speed streaks observed in Fig. 3.
How the maximal values of different components of Reynolds normal stresses vary with the element spacing l is quantitatively examined in Fig. 11. It is seen that all the three components of the Reynolds normal stresses in general decrease for the element spacing l/r from 2.0 to 2.8, except for v 2 + max from the cases Dn45A. When increasing l/r further from 2.8 to 3.5, different trends are observed for different orientations and arrangements. For the surfaces with cube rotated along the spanwise direction by 45 • and placed in an aligned way (i.e., Ds45A), the values of u 2 + max , v 2 + max and w 2 + max decrease at a very low rate, which is approximately the same for l/r from 2.0 to 2.8 and from 2.8 to 3.5 for u 2 + max and v 2 + max . For the rough surfaces with cubes placed vertically (D0A and D0S), the values of u 2 + max , v 2 + max and w 2 + max decrease at a lower rate for l/r from 2.8 to 3.5 when compared with that for l/r from 2.0 to 2.8. For other cases, i.e., Ds45S, Dn45A and Dn45S, the values of u 2 + max , v 2 + max and w 2 + max increase for l/r from 2.8 to 3.5, except for w 2 + max the cases Ds45S. Placing the elements in a staggered way in general decreases all the three components of the Reynolds normal stresses when compared with the corresponding aligned arrangement. One exception is the Dn45A and Dn45S cases, in which the cube is rotated along the spanwise Fig. 10 Vertical profiles of the dispersive stress ũ 2 + s , ṽ 2 + s and w 2 + s computed from the cases with different rough surfaces. The left vertical axis is for ũ 2 + s , and right vertical axis is for ṽ 2 + s and w 2 + s . The parameters for each rough wall are shown in Table 1. The vertical grey dash lines show the location of the roughness crest axis by 45 • . It should be noted that the sandgrain roughness lengths for the Dn45A and Dn45S cases are approximately the same with each other. However, placing in a staggered way does decrease the maximums of the streamwise and the vertical components of the Reynolds stresses. On the other hand, the maximums of the spanwise Reynolds normal stresses from the Dn45A and Dn45S cases are approximately the same for the considered element spacings.
The variations of the maximal values of the streamwise component of the dispersive stresses ũ 2 + max via the element spacing l are shown in Fig. 12. It is seen that the values of ũ 2 + max increase with l for all the considered rough surfaces. The increase rates are lower for the D0S and Dn45S cases when compared with other cases. The highest values of ũ 2 + max are observed in the D0A cases for all the considered element spacings. Similar to the Reynolds normal stresses, placing roughness elements in a staggered way in general decreases the maximums of the streamwise dispersive stresses.

Turbulence Statistics in the Roughness Element's Wake
In this section, the turbulence statistics in the wake of the roughness element are examined. Instead of using the friction velocity for normalization, the bulk velocity U b (which is a representative of the outer flow scale) is employed for normalizing the mean streamwise velocity, the turbulence kinetic energy, and the Reynolds shear stresses. The time-averaged vertical profiles are further averaged over all the elements and compared at four different locations x r downstream of the cube center as shown in Fig. 13. The comparison is focused on the effect of element spacing on flow statistics that the results from the cases with different element orientations are compared separately.
The vertical profiles of the time-averaged streamwise velocity from the aligned cases are compared in Fig. 14. The region with reverse flow is clearly shown, with its height decreasing as traveling downstream. As a result, increasing the element spacing increases the frontal area and the velocity magnitude of the downstream element facing the forward flowing flow, and thus increases the form drag of the element. At x r = 1.0r , the height of this recirculation region is somewhat smaller for smaller element spacing. Nevertheless, the overall distributions are similar to each other for different element spacings, except when approaching the downstream element as it happens at different locations for different element spacings. For the comparison between different element orientations, rotating the element along the spanwise direction by 45 • significantly change the vertical extent of the wake region because of the decrease in the vertical height.
The comparisons of the TKE (turbulence kinetic energy) profiles are shown in Fig. 15 for the aligned cases. It is seen that the vertical distributions of TKE from cases with different element spacing are close with each other at x r = 1.0r for the D0A and Dn45A cases. For the Ds45A cases, on the other hand, the maximal value of TKE from the l/r = 2.0 case is smaller when compared with those from the cases with l/r = 2.8 and l/r = 3.5, which are close with each other. For the cases with element spacing l/r = 2.8 and l/r = 3.5, the magnitudes of TKE increase as traveling downstream from x r = 1.0r to x r = 2.0r . It is interesting to see that the distributions of TKE from the cases with l/r = 2.8 and l/r = 3.5 are fairly close to each other at x r = 1.0r and x r = 2.0r . Comparison between cases with different element orientations shows that the maximum of the TKE from the Ds45A case is slightly higher than the D0A and Dn45A cases.
The vertical profiles of the Reynolds shear stresses are compared in Fig. 16 for the aligned cases. It is seen that the maximum magnitude of u v occurs just above the element crest as expected. The region with high magnitude u v expands as traveling in the element's   downstream. Similar to the TKE and the mean streamwise velocity, the vertical distributions are approximately the same with each other at x r = 1.0r , 1.5r , 2.0r for the cases with l/r = 2.0, 2.8. That the vertical profiles from the case with l/r = 2.0r are different from the other two cases at x r = 1.5 is simply because this location is close to the edge of the downstream roughness element. At locations where the influence of the downstream element is less significant, the vertical profiles from the case with l/r = 2.0r are close to the other two cases as well (not shown here).
After showing the results from the aligned cases, in the following we examine the turbulent statistics in the wake of the roughness element for the staggered cases. The vertical profiles from x r = r to x r = 5r are examined as the streamwise spacing between an element and the one in its direct downstream is two times of the spacing between two neighbor rows (Fig. 13). In Fig. 17, we compare the vertical profiles of the mean streamwise velocity. The same as the aligned cases that the incoming velocity for the downstream rough element increases with the increase in the roughness element spacing, as a result of the wake recovery. The recirculation region is clearly shown at x r = r for all the staggered cases, then almost disappears at x r = 2r , and reappears at x r = 3D and x r = 5D for the D0S cases with l/r = 2.0 and l/r = 2/8, respectively. The disappearance at x r = 2r for the D0S cases with l/r = 2.0 (i.e., the C20D0S case) is caused by the flow acceleration because of the two elements on In Fig. 18, the vertical profiles of TKE in the wake of the roughness element are compared for the staggered cases. It is seen that the magnitudes of TKE below the roughness element crest are relatively higher when compared with those of aligned cases. The other local peak is observed below the element crest in addition to the one above for the D0S and Dn45S cases at x r = 2r , 3r , which is not clearly observed in the aligned cases. For the cases with roughness elements non-rotated (D0S) and rotated along vertical axis by 45 • (Dn45S), the TKE magnitudes from the l/r = 2.0 cases are significantly lower than those from the l/r = 2.8 and l/r = 3.5 cases, which are close for the later two cases at x r = 2r , 3r . For the Ds45S cases, on the other hand, the TKE distributions are fairly close to each other for The Reynolds shear stresses are compared in Fig. 19. In general, the magnitude of the Reynolds shear stresses increases with the increase in element spacing for Dn45S cases. For the D0S cases, the vertical distributions of the Reynolds shear stresses are close to each other for the cases with l/r = 2.8, 3.5. For the Ds45S cases, on the other hand, the vertical profiles of u v are close with each other for all the considered element spacings at all the considered downstream locations. The magnitude of the Reynolds shear stresses is maximal at x r = 2r of the considered five locations for the l/r = 2.8 and l/r = 3.5 cases. Overall, the profiles of turbulent statistics from the staggered cases deviate more from each other for different element spacing when compared with the aligned cases.

Discussions on Engineering Models for k s
In this section, we apply three engineering models for k s to the simulated cases, and evaluate different geometrical parameters of the rough surface on scaling k s , with the objective to provide some insights on developing a model of k s applicable to a wide range of rough surfaces.
Different empirical formulae were proposed in the literature . The empirical formulae for k s employed in this paper are listed as follows: where E S = 2λ f is effective slope. 2. The empirical formula proposed by Chan et al. (2015), who studied 3D sinusoidal roughness and related the obtained sandgrain roughness with k avg and E S, in the following form: 3. The empirical formula proposed by Macdonald et al. (2002) for cube arrays with aligned and staggered arrangements, with d/k c = 1 + A −λ p λ p − 1 , where β = 0.55 and A = 3.36 for aligned arrangement, β = 1 and A = 4.43 for staggered arrangement.
The values of k s predicted by the empirical formulae are compared with the DNS data obtained in this work. As seen in Fig. 20, the overall variations of k s as a function of the element  Chan et al. (2015) model for almost all the cases, especially for l/r = 2.8 and l/r = 3.5. Such underprediction is probably due to the fact that Chan et al. (2015) model was formulated for the rough wall formed using the sinusoidal roughness elements, for which the wake effects may not be as strong as for cubes especially for large element spacings. For Macdonald et al. (2002) model, a very good agreement is observed for the D0A cases with l/r = 2.8 and l/r = 3.5, while discrepancies are observed for almost all other cases, that further development are needed to better take into account the effects of element arrangement and orientations.
A characteristic length defined using the topography parameters of the rough surface may do a better job for scaling k s from different cases, which were normalized using h in Fig. 7. In Fig. 21, the variations of k s via l, which are normalized using different length scales, are presented. The considered length scales include k avg , k c and k rms (listed in Table 2) and their product with S k and λ f , i.e., k avg S k , k c S k and k rms S k , and k avg λ f , k c λ f and k rms λ f , respectively. It is seen that none of the considered length scales can collapse all the data, with one exception for the Ds45A and the Ds45S cases that the variations of k s with l overlap with each other independent of the employed length scale. One interesting observation is that Fig. 19 Vertical profiles of Reynolds shear stresses at different locations downstream of the roughness element for the staggered cases. The first row shows the cases with non-rotated elements (D0S); the second row shows the cases with elements rotated along the vertical axis by 45 • (Dn45S); the third row shows the cases with elements rotated along spanwise axis by 45 • (Ds45S). The vertical grey dash lines indicate zero Reynolds shear stresses. The horizontal grey dash lines show the location of the roughness crest employing the length scales k avg λ f , k c λ f , and k rms λ f collapse the D0A and Dn45A data, and the D0S and Dn45S data, respectively, indicating the effect of rotating along the vertical axis by 45 • is properly accounted for by the three length scales.
Overall, we have seen that none of the considered empirical formulae and length scales work for all the cases. On the other hand, some results are encouraging, including the good agreement for some empirical formulae (e.g., Forooghi et al. (2017)'s model for the overall trends and Macdonald et al. (2002)'s model for the D0A cases with l/r = 2.8 and l/r = 3.5) and the possibility of employing k avg λ f , k c λ f or k rms λ f to account for the element orientation effect.

Conclusions
In this work, we carried out roughness-resolved direct numerical simulations of turbulent channel flows with cube-roughened walls. The roughness elements are arranged in two different ways, i.e., aligned and staggered. For each arrangement, the cube is rotated either along the vertical axis or along the spanwise axis by 45 • in addition to the vertically placed The results show that the orientations and arrangements of the cube elements affect the flow structures in the roughness sublayer, e.g., high-speed and low-speed streaks, upward motion at the windward side of the cube, and the cube wakes. Changing the element spacing affects the strength and extent of these flow patterns. By placing the elements in a staggered way, the high-speed and low-speed streaks become short, with intense fluid motions in the spanwise direction.
Increasing the element spacing affects the turbulence statistics in different ways for different element arrangements and orientations. Placing the roughness elements in a staggered way in general increases the roughness function. Increasing the element spacing in general increases the roughness function for the staggered layouts except for the one with the roughness elements rotated forward by 45 • . The normalized Reynolds stresses and dispersive stresses were then examined for different rough walls. The maximal values of the normalized Reynolds normal stresses from the aligned layout are in general higher than their staggered counterparts. Increasing the element spacing decreases the normalized Reynolds stresses, while increases the maximal values of the streamwise component of the normalized dispersive stresses ũ 2 + s and the width of the region with high magnitude ũ 2 + s . The magnitudes of Fig. 21 Equivalent Sandgrain roughness k s normalized by different length scales for a k avg , b k avg S k , c k avg λ f , d k c , e k c S k , f k c λ f , g k rms , h k rms S k , and i k rms λ f , respectively the other two components of the normalized dispersive stresses are lower than the streamwise component, and are noticeable only for the staggered layouts. For the vertical profiles of the mean streamwise velocity, turbulence kinetic energy and Reynolds shear stresses in the wake of the roughness element, differences are observed for different element arrangements and orientations, with certain similarities observed between cases with different element spacings for the same element arrangement and orientation at locations less influenced by the downstream element. The variations of the turbulence statistics are different for different ranges of element spacings. As for the normalized Reynolds normal stresses, their maximal values u 2 + s max , v 2 + s max and w 2 + s max in general decrease for the element spacing l/r from 2.0 to 2.8. For l/r from 2.8 to 3.5, different trends are observed for cubes with different orientations and arrangements. The maximal values of the streamwise component of the dispersive stresses, on the other hand, monotonically increase with the all the considered element spacings. The equivalent sandgrain roughness increases for l/r from 2.0 to 2.8 and in general decreases for l/r from 2.8 to 3.5, respectively. The underlying mechanism for the variations of k s with l is analyzed via two competing factors, i.e., the increase of the form drag on each individual roughness element via l, which depends on the area of the frontal surface facing the fluid flowing forward and the incoming velocity magnitude, and the decrease in the number of roughness element via l. As for the zero-plane displacement d, it in general decreases with the element spacing for the considered rough surfaces with three exceptions. Tests of three different engineering models and different characteristic length scales for k s show that none of them work for all the cases, e.g., capturing the variations of k s with l, and accounting for the effects of element arrangements and orientations. One interesting observation is that including λ f for defining the length scale (e.g., k avg λ f , k c λ f and k rms λ f ) can account for the effect of rotating the cube element along the vertical axis by 45 • .
Although the considered roughness element, i.e., a cube, is rather ideal, it is closely related to the practical applications, such as, the atmospheric turbulent flow over urban canopies. The obtained results highlight the importance of accounting for the effects of element arrangement and orientation for different element spacings when developing engineering models for rough wall turbulence. Studies to include these effects in engineering models and carrying out related experiments could be done in the future work.