## 1. Introduction

Understanding turbulent flows over rough walls is paramount to a wide range of practical problems. In engineering applications, they manifest in internal flow systems such as in pipes and ducts, marine transportation or turbomachinery (Shockling, Allen & Smits Reference Shockling, Allen and Smits2006; Bons Reference Bons2010; Monty *et al.* Reference Monty, Dogan, Hanson, Scardino, Ganapathisubramani and Hutchins2016). They are also ubiquitous in natural environments such as in the atmospheric boundary-layer (ABL), urban areas or in fluvial processes (Nezu, Tominaga & Nakagawa Reference Nezu, Tominaga and Nakagawa1993; Garratt Reference Garratt1994; Cheng & Castro Reference Cheng and Castro2002; Nikora & Roy Reference Nikora and Roy2012). In all the aforementioned examples, roughness plays a fundamental role, influencing the overall drag, turbulent mixing and transport properties (Jiménez Reference Jiménez2004).

One of the long-standing questions remains the characterisation of the drag by means of the surface properties; however, there is a plethora of arrangements and shapes under which roughness can be present. Following Grinvald & Nikora (Reference Grinvald and Nikora1988) and Stewart *et al.* (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019) definitions, roughness can be classified into two categories. It can be considered discrete (regular), when presented as an amalgamation of discrete elements defined by a combination of linear length scales such as height, width and pitch (e.g. cube roughness). Alternatively, it can be interpreted as continuous (irregular), when the roughness is randomly distributed among a wider range of scales and defined by its statistical moments such as standard deviation, skewness and flatness (e.g. naturally corroded pipes). The seminal work of Nikuradse (Reference Nikuradse1933) followed by the comprehensive reviews by Schlichting (Reference Schlichting1937) and Colebrook *et al.* (Reference Colebrook, Blench, Chatley, Essex, Finniecome, Lacey, Williamson and Macdonald1939) served as the basis for the early drag predictions performed by Moody (Reference Moody1944) for pipe flows, using the equivalent sandgrain roughness height $h_{s}$. The latter has since become the common currency for predicting drag, with the aim of correlating surface properties to this aerodynamic roughness length scale $h_{s}$. Several studies have attempted to provide relations that would predict drag at high Reynolds numbers based on either laboratory measurements or numerical simulations. However, the existence of a universal correlation that would generically represent any rough surface is unlikely, due to its simultaneous dependence on various parameters such as front and plan solidities (Van Rij, Belnap & Ligrani Reference Van Rij, Belnap and Ligrani2002), effective slope (Napoli, Armenio & De Marchis Reference Napoli, Armenio and De Marchis2008), roughness skewness (Flack & Schultz Reference Flack and Schultz2010), roughness directionality (Nugroho, Hutchins & Monty Reference Nugroho, Hutchins and Monty2013), spanwise heterogeneity (Anderson *et al.* Reference Anderson, Barros, Christensen and Awasthi2015) among other parameters.

There exists an extensive body of research dedicated to the study of flows over idealised two- and three-dimensional topographies such as cubes, bars, pyramids or spheres (see e.g. Cheng & Castro Reference Cheng and Castro2002; Schultz & Flack Reference Schultz and Flack2009; Volino, Schultz & Flack Reference Volino, Schultz and Flack2009). However, these surfaces have typically a limited range of scales through which they can interact with the turbulent flows. Fractal-like geometries on the other hand, have an inherent capacity to interact with turbulent flows through a much broader range of length scales. These flows have received special attention in the previous two decades from both the numerical and experimental communities, such as in flows encountering fractal grids or trees (see e.g. Meneveau & Katz Reference Meneveau and Katz2000; Mazzi & Vassilicos Reference Mazzi and Vassilicos2004; Chester, Meneveau & Parlange Reference Chester, Meneveau and Parlange2007; Hurst & Vassilicos Reference Hurst and Vassilicos2007; Gomes-Fernandes, Ganapathisubramani & Vassilicos Reference Gomes-Fernandes, Ganapathisubramani and Vassilicos2012; Graham & Meneveau Reference Graham and Meneveau2012). Particularly, realistic surfaces such as urban areas and natural landscapes are generally multiscale, exhibiting a fractal-like behaviour with a height spectral slope ranging between $-$1 and $-$3 (Huang & Turcotte Reference Huang and Turcotte1989; Passalacqua *et al.* Reference Passalacqua, Porté-Agel, Foufoula-Georgiou and Paola2006; Wan & Porté-Agel Reference Wan and Porté-Agel2011). These have only until the last decade started earning a wider attention from researchers, especially the modelling community. The readers are encouraged to refer to the recent study by Zhu & Anderson (Reference Zhu and Anderson2018) which offers an excellent and comprehensive review of the different contributions in the literature. One of the challenges when examining these flows is the determination of the smallest spatial resolution needed in order to accurately predict the bulk-flow quantities over surfaces featuring a multitude of scales. In fact, unless direct numerical simulations (DNS) or well-controlled laboratory measurements are accessible, computational fluid dynamics (CFD) methods such as large-eddy simulations (LES) or Reynolds-averaged Naviers–Stokes (RANS) require modelling of the unresolved roughness features.

Both ‘discrete’ and ‘continuous’ types of rough surfaces have recently been documented within the multiscale roughness framework. Using LES data obtained from an ABL flow over various irregular topographies, Wan & Porté-Agel (Reference Wan and Porté-Agel2011) examined the effect of the subgrid-scale roughness on the effective aerodynamic roughness length scale, and found the latter quantity to linearly increase with the root-mean-square (r.m.s.) of the unresolved roughness topography. In an LES study of a channel flow over multiscale fractal-like surfaces, Anderson & Meneveau (Reference Anderson and Meneveau2011) developed a dynamic roughness model whereby the unresolved small-scale topography fluctuations are modelled using a local equilibrium wall model. The effective roughness length scale is parameterised through the r.m.s. height and a roughness parameter determined *a posteriori* and validated with an invariance resolution test. They showed that, by including more and more subgrid-scale roughness, the wall shear stress increased, indicating the importance of correctly modelling the unresolved roughness fluctuations.

In an effort to examine the effect of the intermediate and small topographical length scales, Mejia-Alvarez & Christensen (Reference Mejia-Alvarez and Christensen2010, Reference Mejia-Alvarez and Christensen2013) experimentally studied a turbulent boundary-layer flow over a highly irregular surface roughness replicated from a turbine-blade damaged by deposition of foreign materials, from which they created filtered models. By comparing the original surface to its lower order representations, they observed the mean and turbulent statistics in the outer region to be self-similar as they are predominately governed by the large topographical scales. However, they reported that the filtered models failed to appropriately reproduce the contributions of the most intense Reynolds shear stress events (sweeps and ejections) within the roughness sublayer (RSL). In a similar approach to Anderson & Meneveau (Reference Anderson and Meneveau2011), Barros, Schultz & Flack (Reference Barros, Schultz and Flack2018) investigated a channel flow over systematically generated multiscale rough surfaces with varying spectral slopes, with predefined statistical moments. They showed that the surfaces with shallower spectral slopes led to higher drag, suggesting that high wavenumbers (small scales) can noticeably contribute to the overall drag. Stewart *et al.* (Reference Stewart, Cameron, Nikora, Zampiron and Marusic2019) have examined three different self-affine bed roughness in two different open-channel flows, and reported similar findings to Barros *et al.* (Reference Barros, Schultz and Flack2018), such that decreasing the roughness spectral exponent leads to higher drag. They explained this behaviour by analytically demonstrating the existence of a relation between the effective slope with the roughness spectral exponent (also validated with experimental data), which increases as the spectral exponent decreases. They suggested that this drag increase is due to contributions from an increase in the flow separation around steeper, scaling-range roughness features.

For the discrete roughness type, Yang & Meneveau (Reference Yang and Meneveau2017) used LES to examine turbulent boundary layers over randomly distributed prism-like roughness elements with a fractal-like size distribution, and investigated the flow response to the addition of smaller roughness generations. To account for the drag that rises from the unresolved scales, they developed a dynamic-roughness model which assumes the flow to become scale-invariant when it interacts with the scale-invariant roughness, which was then validated by the independence of the mean-flow profiles on grid resolution. They also proposed an analytical model that explicitly accounts for the sheltering effect within the canopy flow; however, its predictions seem to depend on the number of generations included in the roughness. In a similar approach to Yang & Meneveau (Reference Yang and Meneveau2017), Zhu & Anderson (Reference Zhu and Anderson2018) performed a series of LES simulations of a turbulent channel flow over cube-like fractal topographies, with varying fractal dimensions (spectral slope). They assessed the drag penalty associated with changing the spectral slope and the number of smaller generations. They showed that the drag associated with the unresolved generations can be parametrised with the first few generations, and proposed a logarithmic law-based model to model their contributions. They further added that the turbulence statistics are predominantly controlled by the first generation of elements.

Despite the burgeoning interest in this area, there is still a need for experimental studies that relate to flows over such complex surfaces. To our knowledge, the systematic investigation of the multiscale roughness hierarchy effect on turbulent flows in regular roughness has only been reported in a limited number of works (e.g. Yang & Meneveau Reference Yang and Meneveau2017; Zhu & Anderson Reference Zhu and Anderson2018; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2019), and remain experimentally largely unexplored. It is also worth mentioning that when simulating turbulent flows over these types of surfaces, the previous LES studies usually rely on channel flow conditions, in which case the conclusions may not be fully applicable, for instance, to external flow systems such as an ABL flow. Therefore, wind-tunnel experiments in addition to the modelling techniques can prove beneficial with the different lines of enquiries mentioned above. In fact, understanding the dynamics of turbulent flows over multiscale rough surfaces will result in appropriate modelling strategies and better predictions of flows such as in urban environments.

To this end, a comprehensive experimental study is designed aimed at examining the effect of roughness-scale hierarchy on turbulent boundary layers. The design used is similar to the one employed by Yang & Meneveau (Reference Yang and Meneveau2017) and Zhu & Anderson (Reference Zhu and Anderson2018) which first builds large-scale cuboids, subsequently adding smaller iterations where their size decreases with a power-law as the number increases (see e.g. the study of Zhu & Anderson (Reference Zhu and Anderson2018) for a detailed description of an iterated function system). The main difference between the current investigation and the former studies lies in the layout of the roughness in the horizontal plane, currently chosen to be uniform instead of randomly distributed. To uncover the effects of the multiscale roughness on the turbulent boundary layer, direct wall drag as well as flow field measurements are conducted using an in-house floating-element drag balance (Ferreira, Rodriguez-Lopez & Ganapathisubramani Reference Ferreira, Rodriguez-Lopez and Ganapathisubramani2018) and particle image velocimetry (PIV). The remainder of the manuscript is presented in three sections. The experimental methodology is described in § 2 depicting the multiscale roughness topography, drag and flow-field measurement techniques. The results and discussion are presented in § 3, analysing the near-wall and outer-flow regions. Finally, a summary and concluding remarks are provided in § 4.

## 2. Experimental methodology

### 2.1. Facility

Measurements were carried out in an open-circuit suction-type wind tunnel at the University of Southampton. The working section follows a 7:1 contraction and extends 4.5 m in the streamwise direction, with a cross-section of $0.6\ \textrm {m} \times 0.9\ \textrm {m} \times 4.5\ \textrm {m}$ in wall-normal and spanwise directions, respectively. The test section was designed with a weak diverging cross-section to allow a constant free stream along the streamwise direction and a growth of a turbulent boundary layer with a nominally zero pressure gradient. The acceleration parameter $K = (({\nu }/{U_{\infty }}) ({\textrm {d}U_{\infty }}/{\textrm {d}x}))$, where $\nu$ and $U_{\infty }$, the air kinematic viscosity and free-stream velocity, respectively, have been observed in the range $1\text {--}5\times 10^{-8}$ from various recent studies performed in this facility (see e.g. Ferreira *et al.* Reference Ferreira, Rodriguez-Lopez and Ganapathisubramani2018; Medjnoun, Vanderwel & Ganapathisubramani Reference Medjnoun, Vanderwel and Ganapathisubramani2018). The turbulent boundary layer grows over a flat surface composed of five equally sized wooden boards onto which the roughness was mounted. In order to assess the skin friction, the boundary-layer plate at the measurement location was cut out to allow the floating-element drag balance to be inserted in. The boundary-layer plates were preceded by a ramp of 0.2 m long inclined by four degrees to the horizontal, ensuring a smooth transition of the flow from the bottom floor of the test section. The free-stream velocity can reach up to 30 ms$^{-1}$, with a turbulence level less than 0.5 $\%$, and was monitored and acquired using the Micromanometer FCO510 by Furness Controls. To account for air-density variations, temperature was also acquired, and its standard deviation for an average run was less than ${\pm }0.5\,^{\circ }\textrm {C}$. The facility schematic as well as the experimental procedures employed in this study are shown in figure 1.

### 2.2. Multiscale roughness

The surface roughness considered in the current investigation is based on a multiscale distribution pattern as illustrated in figure 2. The design is inspired by a Sierpinski carpet model, which essentially builds the roughness by superimposing size-decreasing self-similar cuboid elements. The model offers an easy control of the geometrical characteristics of the roughness, allowing numerical simulations to be replicated in wind tunnels and *vice versa*. To generate the fractal-like roughness, a similar method employed by Zhu & Anderson (Reference Zhu and Anderson2018) is used, which considers an iterated function system to construct the successive generations with $h^{(i+1)} = rh^{(i)}$, with $h^{(i)}$ and $r$ being the cuboid height at the $i$th iteration and the scale-reduction factor, respectively. The cuboid-side width is defined as $w^{(i)} = A_{r} h^{(i)}$, with $A_{r}$ being the width-to-height aspect ratio. These roughness elements are uniformly distributed in the horizontal $(x,z)$-plane with each individual cuboid of a given iteration being equidistant from its neighbouring ones, that is $S_x^{(i)} = S_z^{(i)}=S^{(i)}$, with $S^{(i)}$ being the streamwise/spanwise spacing between two successive elements at the *i*th iteration. For the first iteration, $S^{(1)} = 2 w^{(1)}$, whereas for the following iterations, $S^{(i)} = w^{(i)}$. It should be noted that for scaling purposes, $S^{(1)}$ will exclusively be referred to as $S$ in the remainder of the manuscript. The subsequent self-similar cuboids following the first iteration are overlaid both at the surrounding as well as on top of the previous parent cuboid.

Due to both roughness manufacturing constraints and flow considerations, the height of the first iteration roughness element was set to $h^{(1)} = 8$ mm, so that $\delta /h^{(1)} \ge 10$. This constrained the number of iterations to three due to the lowest cuboid height that can be manufactured. The scale-reduction factor $r$ and the width-to-height aspect ratio $A_{r}$ were fixed to 1/4 and 4.2, respectively, whereas the number of elements was set to 1, 36 and 576 for the three successive iterations, namely, large, intermediate and small cuboids, respectively. The roughness elements were confined in a $100\ \textrm {mm}\times 100\ \textrm {mm}$ repeating unit tile, which was 3-D printed then replicated using a moulding-casting manufacturing process to cover the floor of the wind-tunnel test section.

In order to investigate the effect of the multiscale roughness on the turbulent boundary-layer flow and assess the contributions of different length scales, four combinations are considered. (i) Iter$_{1}$: which considers the presence of a single iteration with a cuboid height $h^{(1)}$. (ii) Iter$_{12}$: which looks at the presence of two iterations with cuboid heights $h^{(1)}$ and $h^{(2)}$. (iii) Iter$_{13}$: which examines the presence of two iterations with cuboid heights $h^{(1)}$ and $h^{(3)}$. Finally, (iv) Iter$_{123}$: which investigates the effect of the presence of all three iterations with cuboid heights $h^{(1)}$, $h^{(2)}$ and $h^{(3)}$. Schematics for the different cases are illustrated in figure 2 while a summary of their different statistical attributes are presented in table 1.

### 2.3. Floating element

The wall shear stress was directly measured by means of a floating-element drag balance, developed and validated by Ferreira *et al.* (Reference Ferreira, Rodriguez-Lopez and Ganapathisubramani2018) in both smooth- and rough-wall flow conditions. This method relies on the determination of the streamwise net force $F_{d}$ acting on a finite and structurally independent area $A$, representative of the surface being investigated. It makes use of the parallel-shift linkage principle, through pairs of bending beam transducers which allow the isolation of both the streamwise load as well as the induced pitching moment. The wall shear stress is then deduced via $\tau _{w} = F_{d}/A$. The floating element used consists of a $200 \times 200$ mm$^{2}$ surface area, on top of which four roughness tiles are placed. The balance was flush-mounted with the wind-tunnel floor and positioned around 3.1 m downstream of the leading edge. Additional care was taken when setting up the surrounding roughness owing to the tight tolerance of the air gap which is only 0.5 mm wide.

To assess whether these surfaces can reach the fully rough regime, the floating element was exposed to a series of nine free-stream speeds ranging from 9 up to 27 ms$^{-1}$. Each acquisition lasted 30 s with a sampling rate of 150 Hz (corresponding to 2500 eddy turnover times at the lowest speed), with a total of five repetitions per speed. Pre- and post-calibrations were performed for each configuration without notable discrepancies. A more detailed discussion on the complete design, acquisition procedure and uncertainty of the measurement technique can be found in Ferreira *et al.* (Reference Ferreira, Rodriguez-Lopez and Ganapathisubramani2018).

### 2.4. Particle image velocimetry

The turbulent boundary-layer flow was diagnosed in both the streamwise-wall-normal plane $(x,y)$ as well as in the cross-plane $(\kern-0.05pt y,z)$, using planar two-dimensional two-component (2D2C) and stereoscopic two-dimensional three-component (2D3C) PIV measurements, respectively. Both types of PIV measurements were performed at approximately 3.2 m downstream of the contraction. As illustrated in figure 1($c$), two 2D2C-$(x,y)$ measurements were conducted at the spanwise symmetry planes $z/S = 0,0.5$ (at the roughness valley and ridge), and two additional 2D3C-$(\kern-0.05pt y,z)$ measurements were also acquired at two successive streamwise symmetry planes, at ${x= 3.2}$ and 3.25 m. The flow was seeded with vaporised glycerol-water particles generated by a Magnum 1200 fog machine, then illuminated with a laser light sheet sourced by a two-pulse Litron Nd:YAG laser operating at 200 mJ. A LaVision optical system for the beam focus/expansion of the light sheet was used, which comprised convex and concave lenses in order to focus the beam, and a cylindrical lens in order to expand the sheet with relatively constant thickness in the measurement plane ($\approx$1 mm thickness).

The particle images were recorded by high-resolution LaVision Imager LX 16 MP CCD cameras fitted with 200 and 300 mm AF Micro Nikon lenses for the 2D2C- and 2D3C-PIV measurements, respectively. One camera was used for the planar-PIV set-up, and positioned at nearly 1 m away from the object plane, whereas two cameras mounted on Scheimpflug adapters to account for the oblique view angle ($\pm 42^{\circ }$) were needed for the stereo-PIV measurements, and placed at nearly 1.3 m from the object plane at either sides of the test section, also depicted in figure 1($a$). A single and a double-sided dual plane calibration target aligned with the laser light sheet were used to determine the mapping function for each set-up, using a third-order polynomial fit. This resulted in a field of view of approximately $0.9\delta \times 1.4\delta$ in the $(x,y)$-plane and $1.3\delta \times \delta$ in the $(\kern-0.05pt y,z)$-plane for the planar and stereoscopic PIV, respectively. In terms of the largest roughness length scale $S$, this was equivalent to $1.1S \times 1.4S$ in the $(x,y)$-plane, whereas in the $(\kern-0.05pt y,z)$-plane it spanned $1.5S \times 1.1S$, for the planar and stereoscopic PIV, respectively.

A total of 3000 statistically independent realisations of image pairs were acquired for each case at 0.6 Hz, with a time delay between pulses of 50 and 20 $\mathrm {\mu }$s for the 2D2C and 2D3C PIV, respectively. In order to allow a comparison between the different surfaces, the PIV measurements were performed at different speeds to obtain a relatively matched roughness function $\Delta U^{+}$, and ensure they reached a fully rough regime (although it will be seen later that Iter$_{1}$ never reached the fully rough condition). The free-stream speeds were adjusted from case to case, with $U_{\infty } = 18.7$, 10.3, 18.5 and 10.2 ms$^{-1}$ for the cases Iter$_{1}$, Iter$_{12}$,Iter$_{13}$ and Iter$_{123}$, respectively. The velocity vector fields were then obtained by interrogating particle images using a decreasing multipass scheme starting from $48\ \textrm {pixels} \times 48\ \textrm {pixels}$ down to a final pass of $16\ \textrm {pixels} \times 16\ \textrm {pixels}$ for the 2D2C and $24\ \textrm {pixels} \times 24\ \textrm {pixels}$ for the 2D3C. Using a 50 $\%$ interrogation window overlap, the resulting effective vector spacing ranged between 0.4 and 0.5 mm for the different cases both in the planar and stereoscopic-PIV measurements.

## 3. Results and discussion

This section focuses on the analysis of results reported from the floating-element drag balance with the data from both planar- and stereoscopic-PIV measurements. Section 3.1 discusses the results of the direct-drag estimates along with the flow topology in the canopy region, with the assessment of the aerodynamic parameters associated with themean-flow profiles. Section 3.2 reveals the multiscale-roughness effects on the flow topology in cross-plane, focusing on the turbulence properties in the outer region as well as the assessment of the outer-layer similarity hypothesis. The results presented in this study are openly accessible on the roughness database http://roughnessdatabase.org/ and the University of Southampton repository at: https://doi.org/10.5258/SOTON/D1765.

### 3.1. Inner region

#### 3.1.1. Surface drag

Results from the floating-element drag balance are presented in figure 3. They describe the response of the wall shear stress to the four multiscale rough surfaces investigated. Figure 3($a$) examines how the skin-friction coefficient $C_{f}$ varies with respect to $Re_{x}$ (with varying free-stream speed) in a linear-logarithmic scale. At the lowest tested Reynolds numbers ($Re_{x} <3 \times 10^{6}$), the roughest surface Iter$_{123}$ is distinctively shown to experience the largest frictional drag, and is noticeably higher than Iter$_{12}$. On the other hand, the surfaces Iter$_{1}$ and Iter$_{13}$ exhibit relatively similar values and are the weakest among the four cases. Beyond a certain Reynolds number ($Re_{x} >3 \times 10^{6}$), the skin-friction coefficients are seen to weakly vary for the cases Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$ with an average variation about the mean less than $\pm$2 $\%$, which is within the uncertainty of the measurements. Conversely, a substantial decay of $C_{f}$ is observed for the Iter$_{1}$ case, perhaps with a steeper slope in comparison with the smooth-wall curve represented by Schlichting's power-law $C_{f}=b\,Re_{x}^{a}$ (Schlichting Reference Schlichting1979), highlighted with the black solid line.

It is also observed that at the highest Reynolds numbers of the facility, Iter$_{12}$ experiences a reduction in $C_{f}$ of approximately 4 $\%$ from its mean value across Reynolds numbers. These two different behaviours demonstrate that Iter$_{1}$ and maybe Iter$_{12}$ still remain transitionally rough, whereas the cases Iter$_{13}$ and Iter$_{123}$ have reached the fully rough regime, i.e. Reynolds number invariance of $C_{f}$ (Jiménez Reference Jiménez2004; Flack & Schultz Reference Flack and Schultz2014). As demonstrated by other studies (see e.g. Napoli *et al.* Reference Napoli, Armenio and De Marchis2008; Yuan & Piomelli Reference Yuan and Piomelli2014), the fully rough regime stems from pressure-drag contributions prevailing compared with the viscous drag. This suggests that the viscous-drag contributions still represent a considerable fraction in the Iter$_{1}$ case (at least within the range of $Re_{x}$ tested herein), as opposed to the other cases, whose pressure-drag contributions dominate owing to the presence of additional roughness scales. Although this result seems counterintuitive since the large cuboid is essentially a sharp-edged bluff body, it should be recalled that the frontal and plan solidities of Iter$_{1}$ are noticeably smaller in comparison with the other cases. This is further investigated in § 3.1.2.

The magnitude of the skin-friction coefficient is also considerably affected by the roughness content of the surface, such that the addition of subsequent smaller scales enhances the shear stress acting at the wall. This is clear when examining the skin-friction coefficient averaged over the different Reynolds numbers $C_{F}$ as shown in figure 3($b$). It indicates that the magnitude of the skin-friction coefficient $C_{F}$ increases monotonically with increasing $\alpha$ (with $\alpha = ({\bar {h}}/{h_{rms}}) \lambda _{f}$ being a non-dimensional roughness parameter proportional to the frontal solidity). There are many alternatives and possible ways that can be used to correlate the drag with roughness statistics. In this instance, the use of this surface parameter is justified by the fact that its variation approaches a linear trend with the drag, while other parameters such as $\lambda _{p}$ would fail to capture a clear trend.

To quantify the difference between these surfaces, a relative drag-increase coefficient $\beta$ (where $\beta = 100\,\% \times C_{F}^{(i)}/C_{F}^{(4)}$ with $i=1\text {--}4$ for Iter$_{1}$, Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$) is shown in figure 3($c$) for the different cases. We observe that the largest cuboid alone is responsible for over 80 $\%$ of the drag of the full surface. It is also shown that the drag increase is higher for Iter$_{12}$ than Iter$_{13}$ (93 $\%$ against 87 $\%$ of drag of the full surface), indicating that the drag increment that stems from the intermediate iteration is more prominent than the one caused by the smaller iteration despite Iter$_{12}$ and Iter$_{13}$ having a similar plan solidity $\lambda _{p}$. Therefore, it can be implied that the contributions to the full multiscale surface from the smallest cuboid roughness alone amount to $\beta ^{(3)} = \beta ^{(123)}-\beta ^{(12)} \approx 7\,\%$. Similarly, the contributions from the intermediate roughness scale alone amount to $\beta ^{(2)} = \beta ^{(123)}-\beta ^{(13)} \approx 12\,\%$. These results are in agreement with the observations reported in the previous LES studies, highlighting the importance of small roughness features and emphasising the attention needed when modelling the unresolved subgrid-scale roughness in numerical simulations (Yang & Meneveau Reference Yang and Meneveau2017; Zhu & Anderson Reference Zhu and Anderson2018).

#### 3.1.2. Flow topology in the canopy

The effect of the multiscale rough surfaces on the flow topology is explored by examining the normalised mean streamwise velocity maps from the planar-PIV measurements, at the peak symmetry plane ($z/S = 0.5$). Two cases, Iter$_{1}$ and Iter$_{123}$, are illustrated for comparison in figure 4. Results show that, away from the wall, the flow is unaffected by the surface condition, while in the roughness-affected layer it undergoes strong changes in the streamwise direction (as shown by the streamline contours). Specifically within the surface canopy, the maps highlight the presence of a separation bubble past the large cuboid, whose size and intensity appear to depend on the surface condition. This is illustrated in figure 4($b$) which shows that by adding the intermediate and smaller scales, Iter$_{123}$ produces a weaker streamwise velocity deficit within the canopy.

The impact of the multiscale roughness on the recirculation region is further assessed by examining the separation length as shown in brown in figures 4($a$) and 4($b$), and highlighted in an appropriate scaling in figure 4($c$). The results show the extent of the zero contour-level curves to be relatively conserved for the different cases, irrespective of the surface condition, and is approximately $3h^{(1)}$. It is, however, reported that the wall-normal extent subtly weakens with addition of iterations. This behaviour is believed to be caused by the increased turbulent mixing and wall drag at the top of the large cuboid, owing to the interaction of a broader range of roughness scales with the separating shear layer. To explore further the effect of the multiscale roughness on turbulence, the normalised Reynolds shear stress $-\overline {uv}^{+}$ and its associated turbulence production $P_{xy} \delta /U_{\tau }^{3}$ (with $P_{xy} = -\overline {uv} \,\textrm {d}U/{\textrm {d} y}$) fields are examined, as shown in panels (*a*–*d*) and panels (*e*–*h*) of figure 5, respectively. For Iter$_{1}$, the shear layer formed at the top surface of the large cuboid separates at its trailing edge, from which strong vortices are shed towards the canopy. By comparing the area within the encapsulated contour level in panels (*a*–*d*) of figure 5, we observe that the separated shear layer carries a weaker vortical activity as more roughness length scales are imposed. However, a plausible explanation for this behaviour can be related to the range over which turbulent and roughness length scales interact with each other. In fact, when a turbulent flow encounters a broader range of scale, such as in the case of Iter$_{123}$, the $h^{(1)}$-scale vortical structures as seen for Iter$_{1}$ break down to smaller eddies. This means that while the overall drag increases when increasing the roughness content, the shear stress activity is redistributed among scales whereby the intermediate and smaller scales are responsible for a larger portion of drag. Consequently, this results in a relatively shorter extent and a weaker separation bubble and more drag emanating from these additional scales. These observations are further substantiated by the turbulence production maps which show stronger shear layers to be associated with substantial turbulence production. In contrast, weaker shear layers are accompanied with weaker turbulence production.

Using the streamwise-wall-normal PIV measurements, the mean pressure distributions are additionally estimated at the peak symmetry plane ($z/S$ = 0.5), by means of the two-dimensional Reynolds-averaged Naviers–Stokes equations. For more details on the methodology as well as the numerical integration schemes employed, the reader should refer to the study by Ferreira & Ganapathisubramani (Reference Ferreira and Ganapathisubramani2020). The mean pressure is expressed by its non-dimensional form as $C_{p} = (P-P_{\infty })/q_{\infty }$, with $P$ and $P_{\infty }$ being the mean static and free-stream pressure, while $q_{\infty }$ is the free-stream dynamic pressure.

An example of a pressure field is shown in figure 6($a$) for the Iter$_{1}$ case. The coefficient of pressure field is shown to be dominated by an alternating high- and low-pressure region in the canopy, whereas it quickly recovers to the free-stream pressure in the outer region ($y \ge 0.2\delta$). In fact, the high pressure is seen to be associated with accelerating flow regions while low pressures are accompanied with decelerating flow regions. More specifically, the highest pressure regions are recorded at the windward side of the cuboid, while the lowest pressure magnitude is shown to occur right past the leading edge of the upper surface of the cuboid, reflecting the presence of a strong shear layer shedding off from the blunt leading edge. This negative pressure is shown to trail downstream, forming a weak patch of low pressure within the recirculation region until the reattachment point. The pressure transitions to a positive magnitude as the flow reaccelerates again past the reattachment point.

The profile of the streamwise pressure difference across the largest cuboid is additionally examined and shown in figure 6($b$). The pressure difference is expressed as $\Delta P (\kern-0.05pt y) = P_{w}(\kern-0.05pt y)-P_{l}(\kern-0.05pt y)$ (with the $w$ and $l$ subscripts referring to the windward and leeward sides of the cuboid, respectively), normalised by the wall-normal-averaged pressure $\langle \Delta P \rangle$, and plotted against the wall-normal distance normalised by the cuboid height. Although the profile is not fully resolved down to the wall (due to the spurious region below 1 mm height), the drag profile shows similarities with the results of urban-like roughness (Ferreira & Ganapathisubramani Reference Ferreira and Ganapathisubramani2020), with a maximum at nearly two-thirds of the cuboid height, and decreasing when getting close to the wall. The profile is subsequently used to get an estimate of the pressure drag produced by the large cuboid. This is done by integrating the streamwise pressure-difference profile over the cross-sectional area $h^{(1)} \times W^{(1)}$ (assuming the pressure distribution around the cuboid is uniform in the spanwise direction). The drag force and the friction velocity over the large cuboid are expressed as

*a*,

*b*)\begin{equation} F^{(1)} = W^{(1)} \int_{0}^{h^{(1)}} \Delta P \,\textrm{d} y, \quad U_{\tau}^{(1)} = \left( \frac{F^{(1)}}{\rho h^{(1)} W^{(1)}} \right)^{1/2}. \end{equation}

with $\rho$ being the air density. The results reveal that the ratio between the friction velocity estimated from the largest cuboid to the total friction velocity estimated with the drag balance is approximately $60\,\%$ for Iter$_{1}$. Assuming this ratio to remain constant throughout the different cases, the pressure-drag contributions from the addition of intermediate and small roughness scales ultimately leads to the fully rough regime. In fact, for the other three surfaces, the intermediate and small roughness length scales proportionally add form-drag contributions to the overall drag, resulting in the Reynolds number invariance of $C_{f}$ observed in figure 3. In contrast, the pressure drag contribution for Iter$_{1}$ that stems from the largest cuboid alone is insufficient to result in a flow that is fully rough. Due to the viscous drag which amounts to nearly $40\,\%$ of the total drag, the skin-friction coefficient $C_{f}$ remains dependant of $Re_{x}$ as seen in figure 3. It should be noted that despite the inherent degree of uncertainty that rises from the PIV-based pressure and the assumption of the spanwise uniformity of the profiles, this method gives a reasonable indication of the expected form drag.

#### 3.1.3. Aerodynamic parameters

In order to estimate the aerodynamic quantities that characterise the rough-wall flow, both the peak and valley symmetry planes ($z/S$ = 0, 0.5) PIV data are used. Figure 7($a$) compares the wall-normal distribution of the mean streamwise velocity at the symmetry planes for Iter$_{1}$. The angled brackets $\langle \,\rangle$ in figure 7($b$) refer to the horizontally averaged velocity between both planes, comparing profiles of the different cases with the smooth-wall. Substantial differences between the peak and valley profiles can be observed from figure 7($a$). Below the canopy layer, the velocity is shown to be higher at the valley than at the peak as expected due to the blockage of the large obstacles. At $y/\delta \approx 0.15$, it is observed that $U_{Peak}=U_{Valley}$, with the subscripts ‘peak’ and ‘valley’ referring to the velocity measured at $z/S= \pm 0.5$ and $z/S = 0$, respectively. However, beyond this point until almost two thirds of the boundary-layer thickness, the velocity above the peak becomes higher than that at the valley. This indicates that the outer flow presents a degree of spanwise heterogeneity, probably caused by the surface condition. It is also observed that at the individual planes, the velocity profiles do not exhibit a clear wall-normal logarithmic distribution. However, when both planes are averaged, a velocity range appears to vary logarithmically, and occurs approximately between $0.2$ and $0.3\delta$. It is further shown from the profiles of the different surfaces plotted in figure 7($b$), that all the cases deviate from the smooth-wall behaviour. More specifically, the richer the roughness content is, the higher the profile deficit becomes.

The mean streamwise velocity profile over the rough surfaces can also be expressed using the modified law-of-the-wall,

where $\kappa$ and $B$ represent the slope and intercept of the logarithmic region, respectively, similar to a smooth-wall. In contrast with a flat smooth wall, rough walls will give rise to the quantities $d$ and $\Delta U^+$ in the form of wall-normal and velocity shifts in the logarithmic region, which are termed the zero-plane displacement and the roughness function, respectively. The former is interpreted as the ‘virtual’ origin representative of the height at which the mean drag acts (Jackson Reference Jackson1981). On the other hand, the latter provides a quantification of the momentum loss (if positive) or gain (if negative) due to surface roughness, and depends on the roughness Reynolds number $h_{s}^{+}$. The wall-normal distribution of the inner-normalised mean streamwise velocity profiles are presented in figure 8, compared with the DNS turbulent boundary-layer profile of Sillero, Jiménez & Moser (Reference Sillero, Jiménez and Moser2013). As expected, the profiles are shown to be affected by the surface condition as they exhibit both wall-normal as well as velocity shifts in comparison with the smooth-wall. The rough-wall profiles appear to have relatively similar vertical shifts. However, it is worth recalling that the free-stream speeds $U_{\infty }$ were not kept constant among cases.

The zero-plane displacement $d$ is estimated by making use of the modified law-of-the-wall. This is achieved by taking the derivative of (3.2) with respect to $y$, to obtain what is called the indicator function, expressed as:

This equation is classically used to determine the extent of the inertial sublayer as well as the logarithmic slope $\kappa$ for smooth-wall flows, once $U_{\tau }$ is known *a priori* (Österlund *et al.* Reference Österlund, Johansson, Nagib and Hites2000). This consequently means that if $\kappa$ is assumed to be universal between smooth and rough walls, $d$ should be only a function of the velocity gradient and the friction velocity. Therefore, to avoid solving a two-point fit equation in the present study, the value of $\kappa$ is assumed constant.

Since the friction velocity $U_{\tau }$ was directly measured through the drag balance, the zero-plane displacement is found as the value that minimises the difference between the left-hand side ($\varXi$) and the right-hand side ($1/\kappa$) of (3.3). The results are illustrated in figure 9($a$) and show the appropriate values that minimise the difference yield a good collapse of $\varXi$ for the different cases in the outer region ($y/\delta >0.2$). In the inner region, $\varXi$ is shown to depend on the surface condition with the occurrence of a peak in intensity coinciding with $y=h^{(1)}$. Unexpectedly, the current profiles are observed to result in negative values of $d$ for the cases Iter$_{1}$, Iter$_{12}$ and Iter$_{123}$, with the exception of Iter$_{13}$ which reported a positive value. These are shown to range between $-$0.2 and 0.3$h^{(1)}$. Although the values observed appear to conflict with the interpretation provided by Jackson (Reference Jackson1981), namely the height where the mean momentum sink occurs (thus $d$ cannot be negative), these values are merely a solution of the logarithmic distribution fitting. In fact, given that both the velocity and wall shear stress were directly measured, a potential reason for this behaviour can arise from the assumed value of the logarithmic slope $1/\kappa$, hypothesised herein to be universal between smooth and rough walls. However, the discussion on the universality of $\kappa$ and $B$ between smooth and rough walls is beyond the scope of the current study. It is further reported that the indicator function exhibits a plateau $\varXi \approx 1/\kappa$ in the range $0.2\le y/\delta \le 0.3$, suggesting the inertial sublayer to have shifted farther away from the wall, consistently for all the cases. This is in contrast with the smooth-wall boundary layers which are shown to not exceed an ISL upper limit of $0.15\delta$ (equivalent to the $0.15 Re_{\tau }$ with $Re_{\tau }=\delta U_{\tau }/\nu$ reported by Marusic *et al.* Reference Marusic, Monty, Hultmark and Smits2013). This reveals that the roughness topography may have severely altered the mean and turbulence-flow structure in the outer region.

Once the zero-plane displacement is determined, the difference between the logarithmic velocity distribution and the measured inner-normalised velocity profiles is examined to obtain

In the highlighted grey region of figure 9($b$), the profiles of $\varPsi$ reach a plateau at zero for the appropriate values of $\Delta U^{+}$. Similarly to the indicator function, the plateau of $\varPsi$ is clearly identified between $0.2$ and $0.3\delta$, indicating that the inertial sublayer of the multiscale rough surfaces occurs farther than that of a smooth wall. The values of $\varPsi$ are reported to be relatively similar between cases; however, it is recalled that the free-stream speeds were not kept constant. Their values are collated in table 2. It should be mentioned that despite the unexpected negative results of the zero-plane displacement, the values of $d$ have only a marginal influence on the roughness function. In fact, by forcing $d=0$, the values of $\Delta U^{+}$ change by less than 2 $\%$. Furthermore, the profiles exhibit a similar logarithmic slope as the smooth-wall, albeit with a noticeably less intense wake in the outer region. This observation suggests that the multiscale rough surfaces might have substantially affected the wake strength parameter $\varPi$, thereby, the outer-flow dynamics.

Once $\Delta U^{+}$ is determined, the equivalent sandgrain roughness height can be deduced using the fully rough asymptote relation of Nikuradse (Reference Nikuradse1933),

The constant $C$ can be empirically determined by exposing the flow to fully rough conditions, and is defined as the difference between the fully rough intercept $B_{r}$ and the smooth-wall intercept $B$ ($C=B_{r}-B$). In the current study, velocity measurements are only available at one free-stream speed (i.e. one Reynolds number). Hence, the value of $C$ is assumed to be 3.5 which means the flow should have met the fully rough conditions (Schlichting Reference Schlichting1979). Following the $C_{f}$ results reported in figure 3, the fully rough condition has been met for the three rough surfaces Iter$_{12}$, Iter$_{12}$ and Iter$_{123}$. At the same time, it is evident from the same plots that Iter$_{1}$ is yet to reach an asymptotic value, hence should not be ascribed a sandgrain roughness value. However, for the purpose of characterising the multiscale rough surface, an $h_{s}$ value will also be assumed for Iter$_{1}$ in order to allow a comparison across all cases.

Equation (3.5) is shown to yield $h_{s}^{+}$ values that range between 400 and 560, which are considerably beyond the $h_{s}^{+} =70$ suggested to be the lower limit for a flow to reach the fully rough regime (Jiménez Reference Jiménez2004). Their values are tabulated in table 2 as a fraction of $\delta$, and range between $5$ and $10\,\%$ of the boundary-layer thickness. Figure 10($a$) shows the mean velocity profiles $\langle U^{+} \rangle$ plotted against the normalised wall-normal distance $(\kern-0.05pt y-d)/h_{s}$, and highlights a good degree of collapse with the blue-dashed line which follows a logarithmic distribution. More specifically, the cases Iter$_{12}$ and Iter$_{123}$ are seen to adequately collapse within $2<(\kern-0.05pt y-d)/h_{s}<3$, whereas Iter$_{1}$ and Iter$_{13}$ collapse better within $3<(\kern-0.05pt y-d)/h_{s}<6$.

Two reasons are believed to cause this behaviour: (i) the friction Reynolds numbers for Iter$_{1}$ and Iter$_{13}$ are nearly a factor of two higher than those of Iter$_{12}$ and Iter$_{123}$, leading to a larger scale separation manifested in a slightly wider log region; (ii) the equivalent sandgrain roughness height for Iter$_{1}$ and Iter$_{13}$ is almost half that of Iter$_{12}$ and Iter$_{123}$, causing the profiles to shift away from the wall. The profiles are further shown to have a weaker departure from the logarithmic distribution in the outer region, confirming the influence of the multiscale rough surface on the outer region. Interestingly, all profiles are shown to have a fully rough intercept $B_{r} \approx 8$. It should be noted that the latter is classically reported to be 8.5 (Schlichting Reference Schlichting1979). However, this difference is believed to stem from the fact that the 8.5 value was originally derived from pipe-flow data, for which a logarithmic slope is usually taken as $\kappa \approx 0.41$. In contrast, boundary-layer flows have a relatively smaller value of $\kappa \approx 0.39$ which subsequently results in a smaller intercept (Nagib, Chauhan & Monkewitz Reference Nagib, Chauhan and Monkewitz2007).

The roughness function (hence the equivalent sandgrain roughness height) is known to depend on various parameters such as the surface texture, shape as well as the plan distribution. There have been numerous proposed correlations between the roughness function and the surface properties (see e.g. Flack & Schultz Reference Flack and Schultz2010, Reference Flack and Schultz2014). In the current study, we also attempt to find a relation between the aerodynamic quantity $h_{s}$ with an appropriate geometrical parameter of the roughness, such that a link between the different multiscale rough surfaces can be established. Similarly to the observations made from figure 3($c$), the variation of the relative roughness length scale increment $\rho _{h_{s}}$, with respect to the geometrical parameter $\alpha$ is presented in figure 10($b$). In line with the $C_{F}$ results, $\rho _{h_{s}}$ also shows a reasonable linear behaviour observed for the range of surfaces tested herein. This is expressed in a non-dimensional form as

where $h_{s}^{(i)}$ and $\alpha ^{(i)}$ being the equivalent sandgrain roughness height and the geometrical parameter for the $i$th surface (with $i = 1,4$ for Iter$_{1}$ to Iter$_{123}$), respectively. This means that for these regular multiscale rough surfaces, the aerodynamic roughness length scale of the successive multiscale surfaces is linearly related to that of the previous iterations, therefore to the first iteration. This suggests that knowledge of the geometrical properties along with the aerodynamic properties of the parent rough surface are sufficient to predict the drag of a given multiscale rough surface. However, it should be noted that in the general context of multiscale roughness, this linear relation may allow the prediction of $h_{s}$ only for a limited range of scales/iterations such as from Iter$_{1}$ to Iter$_{1234}$ (with a fourth generation of smaller cuboids). In fact, if we keep adding small-scale iterations, the expected drag (hence its associated roughness length scale) that can be generated by the surface will only weakly increase (as if it reaches a certain asymptotic value of $h_{s}$), and can perhaps be more appropriately described by a power-law. In the current study, the limitation in roughness manufacturing as well as the sensitivity of the drag balance to such small magnitudes prevented the assessment of the limits of the observed linear behaviour. It is important to stress that despite the apparent suitability of this relationship for the discrete-like roughness, it currently lacks generalisation over the broader range of rough surfaces (e.g. the highly irregular continuous-type roughness). Nonetheless, while more testing is required to assess its applicability across a wider range, this relationship remains a potential candidate for a subset of topographies such as the multiscale random/regular urban-like rough surfaces.

### 3.2. Outer region

#### 3.2.1. Secondary motions

Figures 11($a$) and 11($b$) show the cross-stream fields of the normalised mean streamwise velocity for the cases Iter$_{1}$ and Iter$_{123}$, respectively. These maps have been obtained by averaging two cross-planes at the streamwise locations $x= 3.2$ and $x=3.25$ m downstream of the contraction. Unlike the streamwise-wall-normal maps shown in figure 4, the mean cross-plane flow exhibits significant spanwise heterogeneity above the canopy forming alternating high- and low-momentum pathways (HMP and LMP) between the peaks and valleys, respectively. This spanwise undulation in the mean flow is shown to be accompanied by two streamwise rotating cells flanking the sides of the large cuboids emphasised with the in-plane velocity vector plot superimposed on top of the maps. These time-averaged structures were identified using the vorticity-signed swirling strength $\lambda _{ci}$ as shown in figure 12 and indicate the presence of considerable streamwise circulation in nearly half the boundary-layer thickness irrespective of the multiscale rough surface.

Figure 12 highlights the rotational sense of the vortical structures, with upwash motions occurring at the valley, while a downwash motion is observed above the peaks. These cross-plane motions are known to be a manifestation of Prandtl's secondary flows of the second kind since they arise from turbulence anisotropy, as opposed to the first kind which stems from mean flow curvature (Prandtl Reference Prandtl1952). These features have recently come to scrutiny and are now well-documented for instance in boundary layers (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2013; Nugroho *et al.* Reference Nugroho, Hutchins and Monty2013; Willingham *et al.* Reference Willingham, Anderson, Christensen and Barros2014; Anderson *et al.* Reference Anderson, Barros, Christensen and Awasthi2015), channels and open-channel flows (Wang & Cheng Reference Wang and Cheng2006; Yang & Anderson Reference Yang and Anderson2017; Chung, Monty & Hutchins Reference Chung, Monty and Hutchins2018; Zampiron, Cameron & Nikora Reference Zampiron, Cameron and Nikora2020) as well as pipe flows (Chan *et al.* Reference Chan, MacDonald, Chung, Hutchins and Ooi2018).

A common feature in all these studies is that when the surface spanwise characteristic length scale is comparable to the dominant length scale of the flow (boundary-layer thickness, channel half-height or pipe radius), large-scale secondary motions manifest with generally upwash and downwash motions occurring above ‘low-’ and ‘high-roughness’ respectively. In the present investigation, the largest surface length scale corresponds to the largest phase in the multiscale roughness which is $S$ in both the streamwise and spanwise direction. This incidentally leads to $\delta /S \approx 1.1$ for all different multiscale rough surfaces, which appears to line up with the previous studies, hence, the observed large-scale secondary motions.

The location of the HMP and LMP has recently been questioned since there has been apparent disagreement between some studies that reported opposite trends, where upwash and downwash motions can occur above low- and high-roughness, respectively. In a recent study by Medjnoun, Vanderwel & Ganapathisubramani (Reference Medjnoun, Vanderwel and Ganapathisubramani2020), it has been highlighted that these conflicting trends are likely to arise from two different driving mechanisms representing non-equivalent types of surface heterogeneity conditions, which can result in opposite observations. It is argued that turbulent secondary flows of the second kind can be generated over two main types of heterogeneous surfaces, namely, topography-driven and skin-friction-driven heterogeneity.

For surfaces purely driven by skin friction heterogeneity (spanwise-uniform height distribution), Willingham *et al.* (Reference Willingham, Anderson, Christensen and Barros2014) and Chung *et al.* (Reference Chung, Monty and Hutchins2018) have demonstrated that HMP and LMP (hence downwash and upwash) systematically occur above regions of high and low skin friction, respectively. For topographically driven heterogeneity (alternating peaks and valleys), both scenarios have been observed. Several numerical and experimental studies showed upwash and downwash to occur above low- and high-roughness (Mejia-Alvarez & Christensen Reference Mejia-Alvarez and Christensen2013; Barros & Christensen Reference Barros and Christensen2014; Yang & Anderson Reference Yang and Anderson2017; Awasthi & Anderson Reference Awasthi and Anderson2018), which is in line with the current observations. Other studies on the other hand highlighted the opposite behaviour (Nezu & Nakagawa Reference Nezu and Nakagawa1984; Wang & Cheng Reference Wang and Cheng2006; Vanderwel & Ganapathisubramani Reference Vanderwel and Ganapathisubramani2015; Hwang & Lee Reference Hwang and Lee2018). However, one of the key differences between these two sets of studies is the presence/absence of streamwise heterogeneity in the topography. For surfaces where HMP and LMP have been reported above peaks and valleys, respectively, the roughness clearly exhibited wake producing protrusions which lead to form-drag contributions, similar to the present study. In contrast, studies that showed upwash and downwash occurring above the elevated and recessed regions, mainly used streamwise homogeneous surfaces. It is worth recalling that the location of the HMP (downwash) and LMP (upwash) is a consequence of the averaging procedure. The instantaneous structures on the other hand have been shown to be able to reverse this behaviour as both scenarios can occur in a non-periodic (chaotic) fashion, as recently highlighted by Anderson (Reference Anderson2019).

Hence, despite the current surfaces being topographically spanwise-heterogeneous, the presence of the streamwise heterogeneity leads the surfaces to act more like the skin-friction-type heterogeneous surfaces due to the additional pressure drag, unlike the other set of studies which used surfaces that are predominantly viscous drag producing surfaces. Therefore, care should be taken when designing numerical simulations or laboratory experiments, by accounting for both the roughness height spectral content as well as the phase information, which under certain conditions can dynamically interact with the $\delta$-scale flow structures, leading large modifications in the boundary-layer characteristics.

#### 3.2.2. Turbulent and dispersive stresses

The effect of the multiscale rough surfaces on the turbulence organisation is examined using the triple decomposition performed on the flow field in the cross-plane. This method is usually employed in flows where spatial heterogeneity is prevalent, in order to quantify the magnitude of the stresses originating from the mean flow inhomogeneity. In this scenario, the dispersive stresses can have substantial contributions to the total stresses and play a major role in the transport of momentum fluxes (Meyers, Ganapathisubramani & Cal Reference Meyers, Ganapathisubramani and Cal2019; Nikora *et al.* Reference Nikora2019). In this scheme, the streamwise velocity field can be written as

where $u_{i}$ is the instantaneous velocity field measured at a fixed streamwise location. Here, $\langle U(\kern-0.05pt y) \rangle _{S}$ is the time and horizontally averaged velocity profile over the spanwise wavelength $S$, $\tilde {u}(\kern-0.05pt y,z)$ is the time-invariant spatial deviation field and $u(\kern-0.05pt y,z,t)$ is the time- and space-dependant fluctuating part from the Reynolds double decomposition. Using this method, the different turbulent and dispersive stress tensor terms can be evaluated and compared. More specifically, the total shear stress term can then be computed as

which typically accounts for viscous, turbulent and dispersive shear stress contributions. The viscous contribution is only present very near the wall, and very negligible away from the wall. On the other hand, Nikora *et al.* (Reference Nikora2019) have shown that the dispersive stresses can be decomposed into roughness-induced and large-scale secondary motion-induced contributions. They have highlighted that the former can have a notable part in the total dispersive stress within the canopy region, whereas the secondary motion-induced part is predominant above the canopy layer. In the present study, we have noticed that these roughness-induced stresses are very negligible throughout the entirety of the resolved flow, as the heterogeneity is mainly driven by the secondary motion-induced stresses.

Contour maps of the turbulent, dispersive and total normal stress terms of Iter$_{1}$ are presented in columns of figure 13, while the shear stress terms are shown in figure 14, respectively. The different normal turbulent stress maps illustrated in figure 13 clearly exhibit spanwise heterogeneity due to the surface condition. In fact, two behaviours can be distinguished, in the near canopy region (at $y/\delta \approx 0.15$), the turbulence intensities are shown to be much higher above the roughness elements while weaker magnitudes are seen between the ridges. In the outer region (at $y/\delta \approx 0.4$), the opposite behaviour occurs with higher intensities occurring between the ridges while weaker magnitudes are located above the roughness elements. These heterogeneous cross-plane turbulence intensity distributions are seen to be related to the rotational sense of the secondary motion, which itself is imposed by the surface condition. The high-turbulence fluctuations observed near the ridges (at $y/\delta \approx 0.15$) are advected laterally towards the mid canopy, and upwards resulting in the LMP. This is further accompanied by the low-turbulence fluctuations being advected from the free stream towards the ridges, inducing the observed HMP above the ridges (clearly emphasised when examining the activity of $\overline {vv}$ distribution). This spanwise turbulence inhomogeneity has previously been shown to be associated with an imbalance between turbulence production and viscous dissipation, causing these amplified lateral motions and sustaining the streamwise vorticity observed in the form of secondary flows (Hinze Reference Hinze1973; Anderson *et al.* Reference Anderson, Barros, Christensen and Awasthi2015; Hwang & Lee Reference Hwang and Lee2018).

Furthermore, the normal dispersive stresses shown in panels (*d*–*f*) of figure 13 are shown to be relatively small with respect to the turbulent ones, however, not negligible. They are shown to extend for nearly half of the boundary-layer thickness although with a varying intensity along the spanwise direction. Their magnitude is shown to be relatively weak in comparison with the turbulent fluctuations (shown with a smaller colourmap scale), while their spatial distribution also exhibits spanwise variation different to that observed for the turbulent components. They are in fact shown to be more localised with upwash regions accompanied with intense activity whereas HMP are associated with smaller magnitudes.

The streamwise and wall-normal dispersive stresses are shown to vanish at spanwise locations between HMP and LMP which coincide with weak upwash/downwash motions. On the other hand, the spanwise dispersive stress is shown to be negligible at the LMP and HMP except in between. This region is seen to coincide with enhanced spanwise fluctuations, albeit with a relatively weaker magnitude in comparison with the normal and streamwise dispersive stresses. The total normal stresses expressed as the sum of the turbulent and dispersive stress (with the viscous stress contributions above the canopy being negligible) is illustrated in panels (*g*–*i*) of figure 13. The results show that globally, the total stresses remain relatively similar to the turbulent ones except with an increased activity between the large roughness elements corresponding to the upwash regions, which are amplified because of the dispersive stress contributions.

For the shear stress terms presented in figure 14, their spatial distributions are also seen to undergo a spanwise modulation due to the wall condition. The spatial distribution of $\overline {uv}$ presents enhanced shear stress activity above the ridges (in the near canopy) while a weaker activity is shown in the outer region. Between the ridges, the opposite trend happens, with weak activity near the canopy and increased shear stress in the outer layer. At the same time, the dispersive shear stress distribution shows a similar behaviour to that reported for the normal dispersive stress term $\overline {\tilde {u}\tilde {u}}$, with localised patches above the ridges and in valleys, while its magnitude weakens in between. Interestingly, a positive dispersive shear stress activity is reported near the canopy, more specifically above the ridges and slightly less at the valley.

The reason for this behaviour comes from the sign of the product of both quantities $\tilde {u}$ and $\tilde {v}$. These quantities ($\tilde {u}$ and $\tilde {v}$) are shown to have opposite signs at the HMP and LMP, with $\tilde {u} > 0$ and $\tilde {v} < 0$ at the HMP while with $\tilde {u} < 0$ and $\tilde {v} > 0$ at the LMP. However, the latter is true except near the canopy where $\tilde {u}$ exhibits a sign reversal such that $\tilde {u} < 0$ right above the ridge (flow experiences increased local friction) and, at a similar height at the valley $\tilde {u} > 0$ (flow experiences weaker local friction), resulting in the observed $\overline {\tilde {u}\tilde {v}} >0$. This leads to a relative increase of the total shear stress $\overline {\tau _{xy}}$ in the outer region and a reduction near the canopy. The turbulent, dispersive and total shear stress distributions of the streamwise-spanwise and the wall-normal-spanwise shear components are shown in figures 14($b$,*e*,*h*) and 14($c$,*f*,*i*), respectively. Similarly to the previous stress terms, their spatial distributions are shown to experience spanwise heterogeneity with varying intensity along the span. The sign of the stress maps is shown to be switching at the symmetry planes ($z/S = 0$ and ${\pm }0.5$) which coincide with locations of zero spanwise fluctuations. The extent of the dispersive shear stress terms is seen to reach at least half of the boundary-layer thickness but with a weaker contribution to the total stresses.

#### 3.2.3. Outer-layer similarity

Following the observations made above, outer-layer similarity is also examined. The similarity between the rough- and smooth-wall structural paradigms stems from an established concept of flow universality known as Townsend's wall-similarity, which states that the influence of viscosity is limited to the near-wall region (Townsend Reference Townsend1976). In this analogy, outer-layer similarity extends this concept to include surface roughness with its validity relying on two main conditions: (i) a large scale separation ($Re_{\tau } \gg 1$) and (ii) a small relative roughness height to the depth of the flow ($h/\delta \ll 1$), as emphasised by Jiménez (Reference Jiménez2004). This is generally expressed in the form of the streamwise velocity deficit and the turbulence quantities which are hypothesised to have a universal form of

*a*,

*b*)\begin{equation} \frac{U_{\infty} - \langle U \rangle}{U_{\tau}}= f\left(\frac{y}{\delta}\right), \quad \frac{\langle\overline{u_iu_j} \rangle}{U^2_{\tau}}= g_{ij}\left(\frac{y}{\delta}\right). \end{equation}

This implies that apart from affecting the near-wall region, the main role of roughness is to set the wall drag ($U_{\tau }$) and adjust the boundary-layer thickness ($\delta$), while the mean and turbulence structure remain unaffected by the surface condition. To examine the impact of the multiscale rough surfaces on the outer-layer similarity hypothesis, the horizontally averaged mean streamwise velocity profiles are plotted in defect form and shown in figure 15($a$), while the streamwise turbulence intensity profiles are depicted in figure 15($b$). The mean velocity profiles indicate a clear lack of collapse in the outer region, showing a weaker deficit than that of the smooth-wall. However, it is interesting to point out that between the cases, a reasonable degree of self-similarity is maintained. The absence of similarity is also shown in the inner-normalised turbulence intensity profiles, both with respect to the smooth-wall profile as well as among cases. The profiles are shown to recover similarity only near the edge of the boundary-layer thickness.

In order to quantify the changes in the outer flow due to the multiscale rough surfaces, the wake strength parameter $\varPi$ is examined. To this end, it is common to assume the flow in the outer region to be fully described by the composite log-wake profile. This is expressed as

where $U_{log\text {-}law}^{+}$ represents (3.2) and $w$ an analytical expression known as Cole's wake function (Coles Reference Coles1956). There are numerous expressions that have been proposed in the literature, and each one yields a value specific to that function (see e.g. Jones, Marusic & Perry Reference Jones, Marusic and Perry2001; Junke, Julien & Meroney Reference Junke, Julien and Meroney2005; Nagib *et al.* Reference Nagib, Chauhan and Monkewitz2007). The wake strength parameter is then determined through fitting the measured velocity profiles to the assumed ‘universal’ expression. In order to circumvent the use of a particular expression for the outer region, the wake strength parameter is instead determined using (3.4), and varies proportionally with the maximum of $\varPsi$. This is expressed as

Their values are summarised in table 2 and are shown to be considerably weaker than that of a smooth-wall (0.57), with $\varPi$ = 0.18, 0.27, 0.20 and 0.23 for the cases Iter$_{1}$, Iter$_{12}$, Iter$_{13}$ and Iter$_{123}$, respectively. These values of $\varPi$ seem to agree with the defect velocity profiles which were shown to have a relatively similar deficit across cases, while being smaller than that of a smooth-wall. The results suggest that universal outer-layer similarity cannot hold for flows developing over this type of surface. In other words, flows that harbour large-scale secondary motions, yielding a highly three-dimensional flow are unlikely to satisfy similarity, as they could have substantially altered the wake intermittency and the uniform momentum zones which together result in the wake strength parameter (Krug, Philip & Marusic Reference Krug, Philip and Marusic2017).

Using the dispersive stress components, the inner-normalised profiles of the total stress terms $\langle \tau _{xx} \rangle ^{+}$ and $\langle \tau _{zz} \rangle ^{+}$ are examined in figure 16($a$), while $\langle \tau _{yy} \rangle ^{+}$ and $\langle \tau _{xy} \rangle ^{+}$ are shown in figure 16($b$). The spanwise-averaged shear stresses $\langle \tau _{xz} \rangle ^{+}$ and $\langle \tau _{yz} \rangle ^{+}$ are null because of the spanwise antisymmetry at the symmetry plane $z/S =0$, as illustrated in figures 14($b$) and 14($c$). Despite the contributions of the dispersive stresses, overall, the normal and the shear stress profiles show a smaller magnitude when compared with their equivalent smooth-wall. This indicates that the increase in the stresses due to the turbulent and dispersive fluctuations caused by the roughness is not necessarily accompanied with a proportional increase in skin friction, since the profiles do not collapse with the smooth-wall (i.e. lack of outer-layer similarity). A degree of collapse seems to be recovered only beyond $y/\delta >0.6$, which coincides with the wall-normal distance at which the dispersive stresses are mitigated as reported in figures 13 and 14.

Figure 16($c$) illustrates the inner-normalised dispersive shear stress profiles with respect to the wall-normal distance normalised by the roughness length scale $h_{s}$. It shows that the maximum magnitude of the dispersive shear stresses rises up to nearly 0.16 $\%$ of the friction velocity, whereas its location changes across cases. It is observed that the wall-normal location of this maximum reduces as the roughness content ($h_{s}$) increases, with $\langle \tilde {u} \tilde {v} \rangle ^{+}_{max}$ occurring at $y \approx$ 2$h_{s}$ for Iter$_{123}$ while its occurrence is seen at $y \approx$ 5$h_{s}$ for Iter$_{1}$. The reason for the latter behaviour rises from the increase in $h_{s}$ when increasing the roughness content. This leads to a lack of similarity between cases suggesting that the increment in the aerodynamic roughness length scale does not necessarily induce a proportional change in the secondary motions. However, it is worth noting that $\langle \tilde {u} \tilde {v} \rangle ^{+}_{max}$ remains nearly constant across cases. Beyond their maximum value, the dispersive shear stresses appear to have reasonably similar decay rates in the outer region, despite the offset in the wall-normal direction.

The quantification of flow heterogeneity for the different surfaces is also examined in outer-scaling as presented in figure 16($d$), illustrating the relative contribution of the dispersive to the total shear stress. The profiles are shown to have a weak magnitude near the canopy, exhibiting even negative values right above the large cuboid for Iter$_{1}$, in line with the observations reported in figure 14. Farther away, the dispersive shear stresses reach their maximum between $0.2$ and $0.25\delta$, corresponding to the wall-normal location at which the secondary motions are most significant. In this scaling, the dispersive stresses are shown to amount to a substantial fraction of the total shear stresses with roughly 22 $\%$ at their maximum, while they seem to have a relatively better collapse in the outer region as opposed to the scaling presented in figure 16($c$). Their contributions rise up to approximately $0.6\delta$, corresponding to the turbulent stresses collapsing with the total shear stress profiles. This means that the inertial sublayer is completely submerged within the roughness sublayer (indicated in figure 16($d$) as the RSL). This behaviour is in contrast with homogeneous rough-wall flows, which generally exhibit a RSL as a small portion of the flow depth, and does not exceed the upper edge of the ISL. The current results also indicate that the dispersive stresses are mostly a function of the largest length scale features of the surface ($S$), while they only weakly vary with addition of small-scale roughness. The multiscale rough surfaces have therefore impacted not only the wall-normal location of the inertial sublayer as seen in figure 9, but also the overall structure organisation in the outer region.

Castro, Segalini & Alfredsson (Reference Castro, Segalini and Alfredsson2013) surveyed a wide range of rough-wall data published in the literature to examine the outer-flow behaviour using the diagnostic-plot method. This method has first been introduced by Alfredsson, Segalini & Örlü (Reference Alfredsson, Segalini and Örlü2011) for smooth-wall flows, and consists of examining the dependence of the streamwise turbulence intensity profile with respect to the mean flow. Their assessment revealed the existence of a region spanning from the inertial sublayer up to almost the entire wake which varies linearly, whose extent increases with increasing scale separation. Accordingly, Castro *et al.* (Reference Castro, Segalini and Alfredsson2013) examined the applicability of this method over a comprehensive range of rough-wall data to reveal whether smooth and rough walls also exhibit outer-layer similarity in this form.

Their results showed that just like the smooth-wall flow, rough-wall flows also presented a linear variation except with a higher slope, shown to be independent of the roughness morphology. This rise in the turbulence intensity level was shown to be primarily due to the roughness function $\Delta U^{+}$, with marginal effect from the wake strength parameter. Using a similar approach, we try to assess whether a self-similar region still exists for the current multiscale rough surfaces using this scaling.

The total normalised streamwise r.m.s. profiles $\langle \sqrt {\tau _{xx}}/ U \rangle = \sqrt { \langle \overline {uu} + \tilde u \tilde u \rangle }/\langle U \rangle$ are plotted against the normalised streamwise velocity $\langle U \rangle /U_{\infty }$, as shown in figure 17($a$). The dispersive stress component has been included to account for the presence of spanwise mean flow heterogeneity, which would have been non-existent in the case of a two-dimensional boundary-layer flow. The profiles of the different cases have been plotted along with a smooth-wall profile and two asymptotes: the smooth-wall asymptote of Alfredsson *et al.* (Reference Alfredsson, Segalini and Örlü2011) and the fully rough asymptote of Castro *et al.* (Reference Castro, Segalini and Alfredsson2013). The results show a good collapse both between the different cases as well as with the fully rough line, indicating that all the cases exhibit similar attributes as a fully rough flow, including Iter$_{1}$. This suggests that despite the lack of collapse observed in the profiles in figure 15, the surface condition influences both the turbulence intensity and the mean-flow proportionally, regardless of the presence of smaller or intermediate roughness length scales. As demonstrated by Castro *et al.* (Reference Castro, Segalini and Alfredsson2013), the increase in the streamwise stress level (slope) for rough-wall flows stems from the roughness function $\Delta U^{+}$. As such, they showed that by rescaling the turbulence-intensity profiles with the modified quantities $U' = U + \Delta U^{+}$ and $U'_{\infty } = U_{\infty } + \Delta U^{+}$, all profiles collapsed on to the smooth-wall asymptote.

After applying the modified scaling approach as seen in figure 16($b$), we notice in fact that the slopes of the different profiles line up relatively well with the smooth-wall asymptote, as reported by Castro *et al.* (Reference Castro, Segalini and Alfredsson2013). It should be noted that if the dispersive stress component $(\langle \overline {\tilde {u}\tilde {u}}\rangle )$ is omitted, the profiles will result in shallower slopes. Figure 17($b$) also suggests that the roughness function embodies both stress contributions, namely, the turbulence which stems from the temporal fluctuations and dispersive stresses which arise from spatial heterogeneities. This means that despite the large ramifications they cause to the base flow, the secondary motions are capable of reorganising the flow such that the mean and turbulence structures maintain a degree of local similarity between smooth and rough walls. However, it is worth recalling that despite this apparent collapse in the diagnostic plot, Townsend's similarity hypothesis remains clearly violated as a consequence of the large-scale secondary motions. Hence, the validity and application of outer-layer similarity hypothesis should be carefully considered when examining flows developing over such surfaces.

## 4. Conclusions

In this experimental study, the characteristics of a turbulent boundary-layer flow developing over regular multiscale rough surfaces have been examined. Four surfaces representing a multiscale distribution of roughness elements, based on a Sierpinski fractal-like pattern, made of superimposing size-decreasing self-similar cuboids are considered. A floating-element drag balance along with PIV measurements allowed us to assess the impact of these multiscale rough surfaces on the skin friction as well as the mean and turbulent flow properties.

Drag-balance measurements revealed that the skin-friction coefficients reached their asymptote with Reynolds number (i.e. fully rough regime) for surfaces containing at least the intermediate or smaller scales (Iter$_{12}$ or Iter$_{13}$), in contrast with the surfaces containing only the largest roughness element (Iter$_{1}$) that was shown to be Reynolds number dependent (despite the high value of $\Delta U^+$). The results showed that the magnitude of the skin friction is dependant of the scale hierarchy of the surface roughness, and increases with frontal solidity of the added roughness ($\alpha$).

The impact of the multiscale roughness on the turbulent boundary layer is also examined using velocity field measurements. In the canopy region, a strong shear layer formed above the largest cuboid separates at its trailing edge leading to the formation of a separation bubble past the obstacle. The streamwise extent of this bubble (reattachment length) was seen to be insensitive to the roughness content while its intensity decreases with the presence of smaller roughness length scales. This behaviour is caused by the increased turbulent mixing and wall drag at the top of the largest cuboid, owing to the increased interactions of the shear layer with a broader range of roughness scales. Additionally, the PIV-based pressure fields showed the canopy to be dominated by alternating high- and low-pressure zones associated with the flow accelerating and decelerating upstream and downstream of the first roughness iteration, respectively. This allowed the assessment of the cross-sectional drag, and revealed that these large-scale roughness elements are responsible for nearly 60 $\%$ of the surface drag.

The aerodynamic parameters are examined using the horizontally averaged mean velocity profiles. The virtual origin as well as the roughness function have been determined by means of the indicator function. The latter highlighted an inertial sublayer occurring between $0.2$ and $0.3 \delta$, which is farther in comparison with the classical smooth and rough walls. The roughness function estimates were then used to determine the equivalent sandgrain roughness height $h_{s}$ whose values indicated that all the cases can be considered in the fully rough regime, including Iter$_{1}$ ($h_{s}^{+} > 400$). A ratio of the aerodynamic roughness length scale between successive iterations ($\rho _{h_{s}}$) was defined and shown to vary linearly with a suitable geometrical parameter of the rough surface ($\alpha$). This relationship together with the aerodynamic properties of the first iteration can be sufficient to describe the roughness length scale of the subsequent multiscale surfaces (and hence the surface drag).

The cross-plane velocity fields revealed a significant spanwise heterogeneity in the outer layer for all surfaces examined. This spanwise flow modulation is shown to be caused by the presence of Prandtl's secondary flows of the second kind, which created HMP and LMP alternating above the ridges and the valleys, respectively. As shown from previous studies, this feature occurs as a consequence of the dominant length scale of the flow interacting with $\delta$-scale roughness features. In this study, the largest characteristic length scale of the roughness is the spacing $S$ between the largest cuboid which results in $\delta /S \approx O(1)$.

Triple decomposition of the velocity fields allowed the qualitative and quantitative assessment of the spanwise heterogeneity in the form of dispersive stresses. The results showed that, at their maximum, these form-induced stresses amount to more than $20\,\%$ of the total stresses and can extend up to 0.6$\delta$. Additionally, the validity of Townsend's similarity hypothesis was assessed by examining the velocity deficit as well as the turbulence intensity profiles. Lack of outer-layer similarity in both the mean and turbulence statistics is reported despite the required conditions being satisfied ($4600 \leq Re_{\tau } \leq 9700$ and $0.05 \leq h_{s}/\delta \leq 0.1$), while a good degree of collapse between cases is observed in the mean flow.

## Supplementary material and movies

All data presented herein as well as the supporting data of this study are openly available on the roughness database http://roughnessdatabase.org/ and the University of Southampton repository at https://doi.org/10.5258/SOTON/D1765.

## Acknowledgements

The authors would like to also express their acknowledgement for the insightful comments and pertinent observations of the referees during the review process.

## Funding

The authors gratefully acknowledge the financial support of the Engineering and Physical Sciences Research Council of the United Kingdom (EPSRC grant Ref No: EP/P021476/1). J.M. acknowledges the support from the Research Foundation – Flanders (FWO grant no V4.255.17N) for his sabbatical research stay at the University of Southampton.

## Declaration of interests

The authors report no conflict of interest.