Urban boundary layers over dense and tall canopies

Wind tunnel experiments were carried out on four urban morphologies: two tall canopies with uniform-height and two super-tall canopies with a large variation in element heights (where the maximum element height is more than double the average canopy height, $h_{max}$=2.5 $h_{avg}$). {The average canopy height and packing density were fixed across the surfaces to $h_{avg} = 80$ mm, and $\lambda_{p} = 0.44$, respectively.} A combination of laser doppler anemometry and direct drag measurements were used to calculate and scale the mean velocity profiles {within the boundary layer depth, $\delta$}. In the uniform-height experiment, the high packing density resulted in a `skimming flow' regime with very little flow penetration into the canopy. This led to a surprisingly shallow roughness sublayer ($z\approx1.15h_{avg}$), and a well-defined inertial sublayer above it. {In the heterogeneous-height canopies, despite the same packing density and average height, the flow features were significantly different.} {The height heterogeneity enhanced mixing thus encouraging deep flow penetration into the canopy. A deeper roughness sublayer was found to exist and extend up to just above the tallest element height (corresponding to $z/h_{avg} = 2.85$)}, which was found to be the dominant lengthscale controlling the flow behaviour. {Results points toward the existence of an inertial sublayer for all surfaces considered herein despite the severity of the surface roughness ($\delta/h_{avg} = 3 - 6.25$)}. This contrasts with previous literature.


Introduction
Understanding flow around urban environments is becoming of increasing importance as cities and their populations grow in size (Department of Economic and Social Affairs, 2019). Although surface energy balance models have recently improved, accurate models for aerodynamic parameters are still poor, particularly for nonconventional roughness geometries -e.g. tall canopies with heterogeneous height, with urban flows that remain poorly described by both theoretical and empirical models (Kanda et al., 2013). The surface layer within the atmospheric boundary layer (ABL) is customarily split into two regions. Firstly, the roughness sublayer (RSL) which acts from the wall surface up to a certain point above the roughness elements. Customarily, the RSL is the region where the flow still encounters the effect of the individual roughness elements (Kent et al., 2017). Secondly, the inertial sublayer (ISL), which covers a region above the RSL (Raupach et al., 1991). The ISL begins at the top of the RSL and is characterised as a region of constant momentum flux (Kent et al., 2017). Considerable debate has taken place over defining the boundaries of both the RSL and ISL. The height of the RSL has been commonly quoted between 2 -5 times the average height (h) of the roughness surface (Raupach et al., 1991), and more recently as low as 1.1 -1.2 h (Florens et al., 2013). The ISL, and its existence, has been subject to even more debate. Traditional theory (Stull, 1988) suggests that in a sufficiently developed boundary layer an ISL should form. In contrast Jiménez (2004) postulated that for δ/h < 80 (where δ is the boundary-layer height), the RSL will increase in depth and effectively replace the ISL. Jiménez (2004) idea has been supported by several studies (Rotach, 1999;Cheng and Castro, 2002;Cheng et al., 2007;Hagishima et al., 2009). However, Leonardi and Castro (2010) and Cheng et al. (2007) argue that the relative boundary layer depth quoted by Jiménez (2004) may not be accurate.
Typically, in the standard characterisation of the effect of wall roughness, a surface (at least close to the wall) is characterised by the aerodynamic parameters and the friction velocity. It is of increasing importance to be able to determine the zero-plane displacement (d) and the roughness length (z 0 ) to effectively predict and calculate wind flow in and above the growing urban environments (as discussed by Stull (1988) in Sec. 9.7). In fully-rough conditions, these parameters are usually obtained, following Cheng and Castro (2002), by (i) calculating the friction velocity (u * ) and then (ii) applying logarithmic law fitting procedure in the constant flux layer (Eq. 1).
The friction velocity can be determined in two ways. Directly, by calculating the drag generated by the rough wall e.g. using fully instrumented elements with static pressure ports or by floating element force balance (Cheng et al., 2007;Hagishima et al., 2009;Zaki et al., 2011). These methods are both acceptable in the fully rough regime, where the viscous drag is almost negligible in comparison with form drag (Leonardi and Castro, 2010). The second approach is indirect, where the friction velocity is evaluated from the Reynolds shear stress in the ISL (Cheng and Castro, 2002;Cheng et al., 2007), so that u * = τ 0 ρ −1 ≈ −u ′ w ′ , where τ 0 is the wall shear stress and ρ is the density of the medium. The roughness length z 0 is the height where wind speed reaches zero, and is used as a scale for the roughness of the surface (Stull, 1988). The zero-plane displacement, d, is understood as a correction factor to the logarithmic profile (Cheng and Castro, 2002;Kanda et al., 2013), while physical reasoning tends to describe it as central height were drag occurs on a rough wall (Jackson, 1981). The aerodynamic parameters vary due to surface properties and, therefore, it should be possible to examine their impact on flow by systematically varying the surface characteristics. Many experiments have taken place over uniform-height cube arrays with different h, λ p , and λ f (see Grimmond and Oke (1999) for λ p and λ f definitions) to examine how roughness affects the aerodynamic parameters (Cheng and Castro, 2002;Cheng et al., 2007;Jackson, 1981;Grimmond and Oke, 1999;Florens et al., 2013;Macdonald et al., 1998;Sharma and García-Mayoral, 2019). However, the discrepancy in these studies with additional evidence from some recent numerical and experimental work, have demonstrated the inaccuracy of using just these variables to examine the sensitivity of the aerodynamic parameters (Hagishima et al., 2009;Zaki et al., 2011;Kanda et al., 2013;Kent et al., 2017;Nakayama et al., 2011;Xie et al., 2008). Previous work has also highlighted how λ p or λ f -in isolation -are insufficient to characterise non-cubical roughness (Carpentieri and Robins, 2015), and advocated for the need to decouple the two solidity ratios Ganapathisubramani, 2015, 2017), however, this is outside the scope of this work. Kanda et al. (2013) pointed out the importance of two additional parameters when describing a canonical regular surface, the standard deviation (σ h ) and maximum element height (h max ). Others have attempted to determine the aerodynamic parameters for various environments by deriving semi-empirical relationships based on roughness geometry (Kanda et al., 2013;Macdonald et al., 1998;Kent et al., 2017). These have the advantage of being able to quickly determine the parameters without the need for field observations, wind tunnel tests, or computational experiments, allowing for fast prediction of flow in urban environments.
This article further explores how the aerodynamic parameters behave in extremely rough surfaces with heterogeneous heights in comparison with uniformheight cases at matching average height. The block arrays in this study are based on simplified attributes of super-tall grid cities, which feature a large standard deviation of element heights with large aspect ratios. A combination of variables, that to the authors' knowledge, have not yet been explored. The structure of this paper is as follows: the experimental set-up is discussed in Sec. 2. Sections 3.1.1 -3.1.4 and 3.2.1 -3.2.4 examine the depths of the boundary layer and surface layers. The aerodynamic parameters are then presented in Sections 3.1.5 and 3.2.5. Finally, conclusions are drawn in Sec. 4.
2 Experimental facility and details 2.1 Experimental facility Experiments were conducted in the 'Aero' tunnel within the EnFlo laboratory at the University of Surrey. This is a closed-circuit wind tunnel with a maximum speed of 40 m s −1 . The free-stream velocity, measured by a Pitot tube located upstream of the model, was set to 10 m s −1 for all cases presented here. The tunnel's test section is 9 m long, 1.27 m wide, and 1.06 m tall. The streamwise, spanwise, and vertical directions are identified with the x, y, and z axis, respectively. The z-axis is set from the top of the baseboard of the model (i.e. the actual wall), while the y = 0 is set in the centre of the test section.
The position x = 0 is considered as the beginning of the tunnel test section. Time-and spanwise-averaged mean, and fluctuating velocities are denoted as (U, V, W ), and (u ′ , v ′ , w ′ ), respectively. In the space between the beginning of the test section and the model a 1 m-long ramp rises from the floor of the tunnel to the average canopy height (h avg = 80 mm, Fig. 4). The ramp creates a smooth transition between the wind tunnel test section inlet and the roughness surface, thus minimising the flow disruption at the beginning of the roughness fetch (Cheng and Castro, 2002) and allowing an equilibrium boundary layer to form. The most upstream measurement station is at x = 3600 mm, where the flow is already fully developed.

Rough-wall models
Four surface roughnesses, all representing idealised tall and super-tall urban environments, were used in this study. Two of the surfaces have elements of uniform-height, while two have elements with varied-height. The individual roughness elements are sharp-edged cuboids of average height 80 mm. Urban buildings are classified by their aspect ratio, AR = h/w (where h and w are the height and width of the building, respectively). Buildings with 3 < AR < 8 fall in the tall regime (generally 100 m < h < 300 m), whilst buildings with AR > 8 (and h > 300 m) are referred to as super-tall (Council on Tall Buildings and Urban Habitat, 2019). Based on this criterion, the surface morphologies examined here are classified as tall or super-tall. Zero-pressure-gradient conditions were used in this work, as the acceleration parameter, calculated for both uniform-and varied-element height cases, was 4.85×10 −8 and 9.48×10 −8 , respectively. The surfaces in examination are further described in the following sections.

Homogeneous-height model
To create an idealised urban environment in the case of uniform-height canopy, two urban features were studied: the packing density, λ p , and the element aspect ratio. Large urban areas with high-density buildings of 14 dense cities were used to calculate a characteristic urban packing density, as reported in Table 1). Google Maps was used to measure the area plots and building base sizes. A range between λ p = 0.33 to 0.57 was determined, which is in line with values cited for real cities in Kent et al. (2017). The packing density of λ p = 0.44 was within this range, and it was selected to describe densely packed cities. The Council on Tall Buildings and Urban Habitat (2015) calculates a mean height of all buildings in Manhattan (New York) above 100 m to be 145.7 m. Using this average height and the bases of buildings from Google (n.d.) gives an aspect ratio of approximately 3.4. Guided by this criterion, the effects of walls and fetch, a scaled model with elements of size 80 mm × 20 mm × 20 mm was designed. The roughness was mounted on five base plates. Each base plate when rotated allowed the roughness pattern to vary from aligned to 50 % staggered, as in Fig. 2. These surfaces are referred to as uniform-height aligned (UHA), and uniform-height staggered (UHS). In the homogeneous height roughness cases, a total of 5,775 elements were used.

Heterogeneous-height model
A varied-height model was also designed that differed from the uniform-height model in only one aspect, the standard deviation of element height (σ h ). The two different roughness configurations are shown in Fig   (2012), when the information was lacking in the database. The real σ h of cities is 17.3 m and was scaled down to σ h = 0.049 m. As in the homogeneous model, the average height is 0.08 m, hereafter denoted as h avg . Other than the σ h = 0.049 m, element heights were selected by matching the distribution of building heights in the data set, resulting in an increased number of short elements with only a few elements taller than the h avg ; h max being 0.2 m. The rough surfaces were constructed by assembling elements into modules. Herein a module is: 3 × 3 elements consisting of five different heights, randomly placed (Fig. 2). In the heterogeneous model, a repeating module is needed to achieve statistically relevant statistics, as described by Cheng and Castro (2002). Each module contains elements with heights: 30, 30, 60, 60, 60, 80, 80, 110, and 200 mm, distributed as in Fig. 3. The purpose of this randomisation is to avoid creating preferential corridors for the flow between the super-tall elements facing the wind direction and to create a more realistic city layout. The tallest elements in the varied-height surface are within the super-tall regime, with AR = 10 (Jianlong et al., 2014). The varied-height surfaces have a total of 5,445 elements. Similarly to the homogeneous canopy, each base plate when rotated allowed for the roughness pattern to be modified from aligned to 50 % staggered. Due to additional geometrical constrains, the total number of elements is slightly different for the heterogeneous and homogeneous canopies, however, this is of little consequence as the downstream flow is spanwise self-similar over the different repeating units and fully developed in a central section of the wind tunnel (i.e. at the measurement stations). The surfaces are referred to as varied-height aligned (VHA), and varied-height staggered (VHS). A summary of the surfaces can be found in Table 2.

Instrumentation
A two-component Laser Doppler Anemometer (LDA), FiberFlow from Dantec, was employed to measure two velocity components of the flow simultaneously. The laser beams were converged by a 300 mm focal length lens, which facilitated measurements in between elements. In the varied-height model, two elements of 200 mm height are located adjacent to each other in some configurations. In these situations, no measurements can be acquired lower than z = 40 mm due to laser obstruction by the elements. The model was sprayed with black matt paint to minimise light reflections. The LDA probe was mounted onto a traverse system which could move the measurement volume in threedimensions within the tunnel with sub-millimetre accuracy. An elliptical mirror was attached to the LDA probe to rotate one of the laser-beam couples by 90 • , effectively changing the component of velocity measured by the LDA. Additionally, a fully instrumented pressure tapped element was used in the case of the uniform-height surfaces to measure the differential force across an element, hence its drag. To allow for these measurements one element was removed and replaced with an identical 3D-printed plastic element, fitted with a total of 25 static pressure ports on one of its faces. The element could then be rotated to allow for the differential pressure to be measured. The drag force was determined by integrating pressure distribution over the front and back of the pressure tapped element, allowing to measure the shear stress due to pressure (u * (p)) as in Cheng and Castro (2002). The u * (p) was only determined for the uniform-height surface given the impossibility of replicating all random permutation of the different height elements for the heterogeneous roughness case.

Experimental details
LDA was used to measure velocity and Reynolds stress profiles within the ISL and RSL above the elements, but also within the canopy. For each surface, different types of LDA profiles were acquired as shown in Fig. 4. For both homogeneous and heterogeneous models x-profiles were collected to study the development of the boundary layer. In the homogeneous model 9 vertical profiles were taken over a repeating element (see dots in Fig. 2) to study the depths of the ISL and RSL. Over the varied-height surfaces 81 vertical profiles were taken over one module to study the statistical convergence of the Reynolds shear stress profiles within the RSL (i.e 3 × 3-module in Fig. 4). Finally, over the heterogeneous model 18 vertical profiles across the width of the tunnel were taken to study the ISL depth (further details on y-planes are contained in Fig. 2).

Results and discussion
This section presents and discusses results for the cases of uniform-height in §3.1 followed by results for the varied-height cases in §3.2.

Depth of the boundary layer
The height of the boundary layer, δ, is commonly estimated as the height where is measured at wind-tunnel inlet, ahead of the rough wall models. Since z 0 is a strong function of roughness fetch in a developing flow (Cheng et al., 2007), the boundary layer must be in equilibrium with the surface below it and fully rough before representative measurements can be collected. It was determined that this occurs over a fetch of x = 4000 mm (≈ 15 δ) for both varied-and uniformheight models. The results past this location are shown in Fig. 5, where δ is shown to be up to the height of 3.25 h and 3 h in the UHS and UHA cases, respectively. The increase in boundary layer thickness in the UHS case is likely due to increased frontal blockage of the staggered array. When the elements are staggered, 'skimming' flow regime (Grimmond and Oke, 1999), is less likely to occur, as the streamwise distances between the element increase, allowing for the development of the 'wake interference' flow regime, with associated enhanced turbulent structures and an increased boundary layer thickness. The current data (h avg = 80 mm and λ p = 0.44) is compared with the previous literature in Table 3. The taller elements in this experiment occupy a much larger fraction of the boundary layer (δ/h max = 3.25) than those of Cheng and Castro (2002) and Cheng et al. (2007), where δ/h max = 12 and δ/h max = 7, respectively; i.e. the surfaces described here have a higher relative roughness height. At the same freestream velocity, the elements of h = 20 mm of Cheng and Castro (2002) and Cheng et al. (2007) show a much larger δ/h max than the one found here, indicating that -not surprisingly -the height of the elements can influence the boundary layer depth. For the sake of brevity, the reader is referred to Andrieux (2017) and Thorpe (2018) for a more indepth discussion of the boundary layer variation across stations in both the streamwise and spanwise direction due to the surface heterogeneity.  Cheng and Castro (2002) and C20A-25%, C20S-25%, C20A-6.25% and C20S-6.25% are from Cheng et al. (2007).

Mean and fluctuating velocity profiles
Mean streamwise velocity profiles for flow over uniform-height cases are shown in Fig. 6a. These profiles have been spatially averaged across one repeating unit to offer a fair representation of the surfaces under investigation and are taken in the region of fully developed flow. The units are non-dimensionalised with the skin friction velocity u * and the roughness length z 0 . For both cases in Fig. 6a, a linear region within the profiles is evident; this offers strong support for the existence of a well defined ISL. Within this region, the collapse of the rough cases onto the smooth theoretical dashed line is an indication of the validity of the assumptions underpinning the methodology used to estimate the roughness parameters. This is further discussed in §3.1.5. The logarithmic region seems to extend well into what would be expected to be the roughness sublayer, as suggested in Cheng and Castro (2002). It must be stressed that data from different streamwise stations were considered for this analysis. The characteristics of the linear region visible in Fig. 6a, and the roughness parameters describing it, were found to be fully converged and independent from the streamwise location (i.e. from the boundary-layer depth) once the flow developed over approximately 15 boundary layer height. This offers a further proof that 15 δ is the minimum required roughness fetch for the flow to be in equilibrium with the underlying surfaces considered herein. The parameters reported in Table 3 are referring to these conditions, where the flow is fully developed. Next the streamwise turbulent fluctuations are discussed in Fig. 6b in the form of the diagnostic plot (Alfredsson andÖrlü, 2010), which removes the need to define the roughness parameters and the friction scaling.
The U e used to normalise the mean velocity in Figures 6 and 10 is the local velocity at the edge of the boundary layer. The benefit of the diagnostic plot approach in this work is twofold, firstly it allows one to assess the universality and the self-similar character of the turbulent statistics, secondly it provides an independent assessment on whether the flow is in fully-rough conditions. For these reasons, the smooth (Alfredsson et al., 2011) and the fully-rough (Castro et al., 2013) asymptotes are reported in Fig. 6b (and Fig. 10b). The data for uniform-heights (both for the aligned and staggered arrays) shows satisfying collapse onto the fully-rough line (black solid line). This is a strong indication that the turbulence statistics are self-similar across cases and that the boundary layer developing above the urban canopies in examination respects the classical scaling laws.

The roughness sublayer
The top of the RSL is considered to be the point where all the effects of the individual roughness elements on the flow cease (Cheng et al., 2007;Kent et al., 2017). It follows that the flow inside the RSL is not spatially homogeneous.
The nine vertical profiles taken over the surface are presented in Fig. 7a,b for the UHA and UHS cases, respectively. These are shown to converge, unexpectedly, at about 1.15 h avg for both cases. The uniform-height surface results for this tall-element canopy differ from the small-cube cases analysed by Cheng and Castro (2002), where the RSL was found to be much deeper (≈ 2 h avg ). Despite the comparable Reynolds number, the RSL depth is much shallower here, closer to those found by Florens et al. (2013). Given the similar conditions between the current data and those of Cheng and Castro (2002), we attribute this difference to the discrepancy in the packing density and the much deeper canopy layer. The current results also vary significantly from the commonly cited 2 -5 h avg (Kent et al., 2017;Flack et al., 2007;Roth, 2000) due to the combination of high aspect ratio elements and dense canopy.
Furthermore, across Fig. 7a and Fig. 7b a clear difference can be seen between cases. Figure 7b show similarity and collapse between all profiles. The flow inside the canopy behaves as predicted in literature, similar to that of a densely forested canopy, with a region of severe velocity deficit up to the elements' average height; this is accordance with Fig. 9.7 in §9.7.3 from Stull (1988). These results are qualitatively in agreement with those in Nepf (2012) for vegetating canopies for which an inflection point in the velocity profiles appears just above the roughness height for dense canopy (i.e. λ f > 0.23). The tightly packed roughness generates a strong shear layer at the top of the elements. Cheng et al. (2007) argues that this region is the main contributor to z 0 , as demonstrated in the 20 mm cube arrays. Figure 7a, however, highlights the effect of aligned street canyons. Profiles 1, 4, and 7 exhibit flow penetration into the canopy where the street canyons are aligned. The rest of the profiles follow a trend very similar to that in Fig. 7b. This demonstrates that, for λ p = 0.44, the flow cannot penetrate deeply into the canopy, unless the streets are aligned. Even so, velocity profiles 1, 4, and 7 only begin to have a significant velocity increase at around 0.5 h avg , indicating that at this λ p , even in aligned street canyons, the penetration into the canopy is limited beyond the depth of 40 mm. A further study could be conducted on varied canopy depth models to verify how deep the flow can penetrate in aligned surfaces with varying h at different packing densities. In summary, uniform tall arrays, concerning the RSL, behave similarly to forested canopies, showing an inflection point in correspondence with the top of the canopy. Lastly, for large λ p and h/δ, even when the elements are aligned, the flow cannot penetrate significantly within the canopy.

The inertial sublayer
In the ISL, the wall-normal variation of shear stress may be neglected, hence the constant-flux region denomination (Kent et al., 2017). Here, we defined the ISL as the region where the vertical variation of the spatially averaged profiles of Reynolds shear stresses is below ± 10 %, as in Cheng et al. (2007). The profiles are reported in Fig. 8a,c. The base of the ISL is assumed to be the top of the RSL found in Section 3.1.3. Unlike the work by Cheng and Castro (2002) and Cheng et al. (2007), where the determination of an ISL was found to be challenging, here an ISL of significant depth can be observed in both cases. Jiménez (2004) reported δ/h ≈ 80 as the lower limit for a canonical ISL to exist. For the uniform-height model, δ/h ≈ 3, well below the limit where an ISL should form. Arguably, in Fig. 8 an ISL forms for these cases. Cheng et al. (2007) also argued that the ISL depth is not constant which is also demonstrated here. As discussed in Section 3.1.1, it is likely that a deep ISL forms due to the large λ p . It is intuitive to imagine that, in the limit of λ p → ∞, a new raised smooth wall surface (at z = h) will form, and recover the canonical turbulent boundary layer structure. Additionally, it was discussed in §3.1.2 how both the mean velocity profiles and the turbulent statistics were found to conform to the canonical scaling for turbulent boundary layers, and how the roughness parameters were invariant with streamwise location once the flow was fully-developed. These findings strengthen the argument that, despite the relative roughness height considered in this work, an ISL does indeed develop over these urban surfaces.

The aerodynamic parameters
The method used in this study to calculate d and z 0 uses a logarithmic law fitting. A common problem faced with the uncertainty in the fitting of the log-law is that there are three free parameters, u * , d and z 0 (Castro, 2007;Segalini et al., 2013). More generally, other unknowns are given by the wake parameter, Π, and by not considering the Von Kármán 'constant' a constant (Castro et al., 2017), however, this is not relevant for the current discussion.
To reduce the uncertainty in the fitting procedure, a common method involves first fixing u * . This is done here with two methods: (i) by using the Reynolds shear stress value in the ISL to determine u * , and (ii) by measuring the drag directly by use of pressure-tapped elements. In the first method, the horizontally-averaged vertical velocity profile is used in either of two ways. Firstly, one can simply compute the average value from all the points within the ISL; secondly, a best-fit extrapolation onto the height z = d or z = h can be used (Cheng and Castro, 2002;Cheng et al., 2007;Florens et al., 2013). An independent method to compute u * is with the use of a pressure instrumented element, as described in Cheng and Castro (2002). Results are summarised in Table 3. Once the friction velocity has been determined, the aerodynamic parameters are computed by best fitting the velocity profiles considering the Von Kármán's constant κ = 0.4, which is within the limits suggested by Marusic et al. (2013) for high Reynolds number flows. The discrepancy between the friction velocities calculated via averaged and extrapolated ISL is below 4 % (Fig. 8b,d). A further comparison extrapolating to z = h (Florens et al., 2013) gives results more closely resembling the pressure tapped results (i.e. ≈ 1.5 % and ≈ 15 % for d and z 0 respectively), suggesting that the extrapolation technique does yield to more robust results. . Black dotted lines show the ISL boundaries determined by the ± 10% from the averaged ISL value. Blue line indicates u ′ w ′ determined by average ISL value, whilst the red line uses best-fit extrapolation to the havg to determine u ′ w ′ (as in Florens et al. (2013)). On the right: Eq. 1 is rearranged to the form of z = z 0 e kU/u * + d, where u * is calculated using u ′ w ′ . Using this rearranged form, z 0 represents the gradient and d is the axis intercept. A least-mean-square fit is used to extrapolate onto the axis within the ISL boundaries. Case UHA for plots (a) and (b), and case UHS for plots (c) and (d).

Depth of the boundary layer
In the varied-height models, similarly to the previous case, an initial investigation was carried out to determine the fetch at which the boundary layer is fully developed and this was found to be at x = 4470 mm (≈ 9 δ), as in Fig.  9. Although the h avg is the same in both models, the boundary layer depth is approximately double in the varied-height model (δ UHA /h avg = 3 versus δ V HA /h avg = 6.25). It is likely that the increased σ h also increases the drag generated by the surface, hence creating deeper boundary layers.

Mean and fluctuating velocity profiles
As for the case of uniform-height previously discussed, mean and fluctuating velocity profiles are considered next. Mean streamwise velocity profiles for flow over varied-height cases are shown in Fig. 10a. As for the previous case, a region of near-linear behaviour is noticeable, however, both the degree of collapse onto the smooth wall dashed line and the extent of the linear regime are reduced from the previous cases in Fig. 6a. Both these aspects are discussed in the following sections, however, these seem to indicate that the methodology for the evaluation of the roughness parameters (see §3.2.4) has a higher uncertainty for surfaces with height heterogeneity and that the roughness sublayer has grown at the expense of the inertial sublayer (see §3.2.3), respectively. It is still important to highlight that, even for this cases, data from different streamwise stations yielded converged roughness parameters; a strong evidence of the fully-developed nature of the flow at the measuring stations. The diagnostic plot for the streamwise fluctuations is shown in Fig. 10b. Both cases for varied-heights are found to be self-similar, however interestingly, these sit above the fully-rough trend line. This can be interpreted as an indication of the enhanced turbulent mixing due to the heterogeneous surface morphology; an aspect further discussed in §3.2.3. This discrepancy with previous data on fully-rough walls reported in Castro et al. (2013) can also be attributed to the possibility that this fully-rough asymptote within the diagnostic plot is effected by several other roughness parameters, as highlighted  (Castro et al., 2013), while the grey line is the smooth-wall limit (Alfredsson et al., 2011). originally in Castro et al. (2013). Two parameters of particular relevance to this work are the much higher relative roughness height, δ/h, and the standard deviation in elements height, σ h , which distinguish this work from that data in Castro et al. (2013).

The roughness sublayer
In the varied-height experiment, 81 profiles were taken per repeating 3 × 3module (Fig. 2). These are shown in Fig. 11a, b for the two configurations. There is a clear collapse of mean vertical profiles, indicating the existence of a limited RSL depth (Fig. 11a,b). The RSL height is the same in both VHA and VHS configurations, which is of interest; it seems that the dominant feature in determining the height of the RSL is the tallest element in the 3 × 3module, as a collapse is found for z/h avg ≈ 2.5 (z/h max ≈ 1). The evident collapse of vertical profiles in the RSL contradicts Cheng and Castro (2002), Cheng et al. (2007) and Jiménez (2004), who predicted an RSL region that expands and 'squeezes' the ISL for rough walls with significantly tall and heterogeneous elements. For the VHA surface (Fig. 11a) not all profiles collapse at the same height. Several vertical profiles converge near the average height of the elements (i.e. z/h avg = 1). The remaining profiles converge just above the maximum-height element (z/h avg = 2.5 or z/h max ≈ 1). This trend possibly occurs for the same reasons discussed for the uniform-height cases. Where street canyons are aligned in the wind direction, the flow penetrates deeper into the canopy, allowing for higher velocities and enhanced mixing. Likewise, in the VHS configuration (Fig. 11b), the flow does not seem to fully converge until above the tallest element (z/h avg = 2.85). In the varied-height situation skimming no longer occurs since there is a large spread of velocities below the h avg and h max . As seen in both the VHA and VHS cases, the large spread of velocities within the canopy indicates significant flow penetration likely due to the increased σ h . Considering this possibility, cities designed with large σ h , could introduce much higher rates of mixing, with positive outcomes on urban air quality, and natural ventilation. As argued by Zaki et al. (2011), λ p is no longer a good indicator of the flow regime when a large σ h is introduced, which is counter to the studies conducted by Sharma and García-Mayoral (2019). The argument made here suggests that the surface parameters used to describe flow over uniform-height models are no longer sufficient to characterise canopies with significant height variations. Figures 11c and d show horizontally-averaged profiles within the RSL for the different experimental runs in the varied-height model. The spatiallyaveraged profiles taken in different areas of the array (across the wind tunnel span) collapse reasonably well onto each other, but there are some discrepancies around the region 1.5 < z/h < 2.5. The profiles measured across the y-planes collapse better than the profiles taken over the 3 × 3-module. Possibly, it is more efficient and easier to get a representative averaged profile from taking measurements across the width of the model rather than over a single repeating module. Another reason for this better collapse could be that the average height of the elements in-front and behind the y-planes did not always match h avg -see further discussion on this topic in Cheng and Castro (2002). However, the clear collapse in y-plane profiles indicates that h avg may not be the dominant lengthscale in models with large σ h . Comparing uniformto varied-height canopies across Fig. 7 and Fig. 11, a prominent difference appears. For the varied-height models, the wind velocities vary greatly within the canopy, in contrast, they remain close to zero until the very top of the canopy for the uniform-height experiments. This is likely due to the physics of the 'skimming flow' regime (Grimmond and Oke, 1999). Even below h avg , there are much higher velocities in the varied-height canopy compared to the uniform-height. The spread of velocity profiles in the varied-height model suggests that reasonable mixing can occur in the near-wall region for tall and dense canopy (large λ p and h avg ), providing that σ h is significant.

The inertial sublayer
Several thresholds have been used in the literature for defining this region, where ± 5, ± 10, and ± 20 % variation are all examined by Cheng et al. (2007). Kanda et al. (2013) used a region for logarithmic law fitting which is a function of the tallest and average roughness elements' height. Sharma and García-Mayoral (2019) noted how the ISL depth in their work is a function of the spacing between the roughness elements. These ISL fitting methods, however controversial, can prove very useful in accurately estimating aerodynamic parameters. However, a convincing relationship between the roughness geometry and the upper limits of the ISL is still elusive, particularly for tall canopies.
The depth of the ISL region, herein based on Kent et al. (2017), is shown in Fig. 12. Previous studies predicted that an ISL region would vanish as the RSL rises through the boundary layer due to increasingly rough surfaces (Cheng et al., 2007;Cheng and Castro, 2002;Rotach, 1999;Jiménez, 2004;Hagishima et al., 2009). The clear collapse of the RSL in Fig. 11 contradicts this theory. Furthermore, the ± 10 % variation definition does allow for an ISL with significant depth to be singled out, as seen in Fig. 12. This may be an indication that the ISL does, indeed, exist. It is important to point out, however, that the Reynolds shear stresses never reach a full plateau, but show a small -yet visible -gradient across the ISL. This is possibly related to the fact that the boundary layer in the cases with varied-height occupies a significant height of the wind tunnel. The acceleration parameter is fairly small (9.48×10 −8 ) indicating a zero-pressure-gradient, however, it must be acknowledged that it is nearly double the same quantity in the uniform-height case. An ISL within the ± 10 % variation does form even though δ/h = 6.25 is still well below the limits established by Jiménez (2004). Florens et al. (2013) argued that this ratio should be lowered, and additionally Cheng et al. (2007) demonstrates that the ratio δ/h should not be the only criterion for the existence of the ISL. Furthermore, Leonardi and Castro (2010) argued that the Reynolds number previously thought to be necessary for the ISL development may be much lower than expected. The finding from the current work, therefore, offers support to the idea that Jiménez's (2004) ratio is perhaps too stringent. The ± 10 % horizontal variation may not, however, be the best suit for determining the ISL depth. The top limit of the ISL can be hard to determine, and it is here tempting to soften the criterion to include the region where the gradient of the slope is more noticeable (i.e. 2.5 < z/h < 5 in Fig. 12c). Indeed, in the uniform-height experiment discussed earlier, the ± 10 % horizontal variation is quite liberal, and a more constrained ± 5 % could be used. The depth of the ISL could perhaps be better determined by the change in slope in the . Black dotted lines show the ISL boundaries determined by the ± 10% from the averaged ISL value. Blue line indicates u ′ w ′ determined by average ISL value, whilst the red line uses best-fit extrapolation to the havg to determine u ′ w ′ (as in Florens et al. (2013)). On the right: Eq. 1 is rearranged to the form of z = z 0 e kU/u * + d, where u * is calculated using u ′ w ′ . Using this rearranged form, z 0 represents the gradient and d is the axis intercept. A least-mean-square fit is used to extrapolate onto the axis within the ISL boundaries. Case VHA for plots (a) and (b), and case VHS for plots (c) and (d).
spatially-averaged Reynolds shear stresses, in a procedure similar to that in Kanda et al. (2013). The relevance of this discussion stems from the fact that, if sufficiently accurate d and z 0 estimates can be derived using the constant flux approximation, then wind-tunnel experiments could be greatly simplified, omitting pressure tapped elements or floating force balances, without compromising on the accuracy of the inner scaling. A summary of the depths of ISL and RSL for both cases of heterogeneous heights can be seen in Table 3. The heights of RSL, ISL, and boundary layer seem to be strongly dependent on the height of the tallest element within the surfaces. Even with appropriate normalisation (e.g. by h avg or h max ), there is no apparent universal scaling between results from varied-and uniform-height experiments. The significant change in the boundary layer, RSL, and ISL depths across the two cases seem to be uncorrelated. Thus, VHA and VHS surfaces could only be compared to models of similarly sized standard deviation and average height.

The aerodynamic parameters
The aerodynamic parameters and friction scaling were determined as in Section 3.1.5, but direct measurement of the friction velocity by an instrumented element were not available for these cases for reasons discussed in Section 2.3.
Here we follow the same procedure explained in Section 3.1.5, i.e. extrapolating to z = h. The results are reported in Fig. 12. The aerodynamic parameters calculated herein are compared to the predictions of Kanda et al.'s (2013) and Macdonald's (2000) morphometric methods which considered elements with varied-height, to assess their performance. There is a large discrepancy between the varied-and uniform-height results, where d is found to increase 1.5 times from uniform-to varied-height, whilst z 0 increases nearly 3.5 times. For the uniform canopy, Kanda et al.'s (2013) method predicted reasonable d while Macdonald's (2000) performed worse (≈ 5 % and 30 %, respectively). Neither produced a value for normalised z 0 that was close to that obtained herein (with both methods predicting a roughness length almost an order of magnitude larger). For the varied-height experiment, the extrapolated results from the shear stress plot produced a value of d above h max , whilst the average shear stress method resulted in values below h max (see Fig. 12). Both Kanda et al. (2013) and Millward-Hopkins et al. (2012) methods predict a value of d between h and h max , which differ ≈ 36 % and ≈ 43 %, respectively from those extrapolated here from the average shear. This discrepancy may be due to the fact that Tokyo (heavily used in Kanda et al. (2013)) does not have many super-tall structures such as Hong Kong (largely due to the seismic restrictions in the former). The lack of similarity of results continues when comparing z 0 . The z 0 values calculated with the two methods were roughly twice as large those evaluated here. The findings of this work strengthen the argument by Kanda et al. (2013) that new parameters are necessary to accurately characterise an urban rough wall, particularly when the elements are tall and heterogeneous in height if one wants to accurately estimate the aerodynamic parameters. Jiang et al. (2008) and Zaki et al. (2011) suggest that d increases with σ h , which is in line with findings presented in this work. Zaki et al. (2011) also demonstrated that the values of d and z 0 almost doubles for cases with the same average height but significant σ h variation. Xie et al. (2008) also acknowledges that the tallest elements have a disproportionate contribution to the drag across a surface. Kanda et al. (2013) further suggested that both λ p and σ h have an effect on d and z 0 . A lack of similarity between results of uniform-and varying-height strengthens the above arguments; particularly the fact that a large σ h produces significant changes in the flow field. Kanda et al. (2013) considered h max as the upper limit of d, and also advocate the use of standard deviation of element height (σ h ) as an important measure when calculating the aerodynamic parameters, which is corroborated by the findings of this work. This seems to support the argument that the sole use of σ h and h avg when designing an urban model are not representative of the flow physics in a real city as suggested by Kanda et al. (2013). Since there are usually a few very tall elements amongst many shorter elements, other important parameters to consider would be the distribution and ratios of tall to short elements.

Conclusions
Wind tunnel experiments were conducted at the University of Surrey on four dense (λ p = 0.44) and tall (δ/h avg ≈ 3) urban arrays; two canopies were of uniform-height and two of varied-height, whilst the average height was kept fixed at h avg = 80 mm in both cases. All canopies aimed to represent idealised modern cities. The experiments examined the differences in the flow features across canopies with homogeneous and heterogeneous heights by means of laser doppler anemometry and direct drag measurements. All cases examined were fully-developed, fully-rough surfaces in zero-pressure-gradients.
In the uniform-height experiments, the surfaces with a large λ p inhibited deep penetration of the wind into the canopy, thus hindering mixing at street level. For λ p = 0.44, 'skimming flow' regime seemingly occurred, and the rough surface started to recover smooth-like properties above the elements, producing an inertial sublayer region; this is in contrast to Cheng and Castro (2002) and Cheng et al. (2007) who questioned the existence of this region as the roughness influence grows. The RSL depth was found to extend approximately up to z = 1.15 h, which is much shallower than the typically expected 2 -5 h (Flack et al., 2007). Conventional turbulent scalings laws were found to be applicable to both mean velocity profiles and turbulent quantities despite the severity of the roughness.
For heterogeneous-height canopies, the usefulness of λ p in describing the wall properties became more questionable. Despite the same h avg , the boundary layer grew to almost double the thickness of the cases with uniform-height, highlighting the significance of σ h and h max in dictating the flow features. Surprisingly, even with such a heterogeneous canopy, a clear collapse of vertical profiles of Reynolds shear stress was observed, forming a coherent RSL extending to just above the height of the tallest element (z/h avg = 2.85 and z/h max ≈ 1), which is much shallower than anticipated. Moreover, more significant wind penetration was observed within the canopy when compared with the uniform-height array. This is responsible for an enhanced turbulent mixing resulting in velocity fluctuations which are higher than previously reported in surfaces with homogeneous heights throughout the boundary layer depth (Castro et al., 2013). These findings strengthen the need to include informa-tion regarding σ h and h max when describing flow over tall and heterogeneous canopies, supporting the conclusions highlighted in Kanda et al. (2013). A comparison of the aerodynamic parameters for the cases with heterogeneous height considered herein with existing morphometric methods has highlighted the inaccuracies of these in the case of tall canopies with significant variation in height between elements. Further work is therefore needed on this topic.