Drag Distribution in Idealized Heterogeneous Urban Environments

Large-eddy simulations of nine idealized heterogeneous urban morphologies with identical building density and frontal area index are used to explore the impact of heterogeneity on urban airflow. The fractal-like urban morphologies were generated with a new open-source Urban Landscape Generator tool (doi:10.5281/zenodo.3747475). The vertical structure of mean flow and the dispersive vertical momentum transport within the roughness sublayer are shown to be strongly influenced by the building morphologies. The friction velocity and displacement height show high correlations with the maximum building height rather than the average height. Well-known roughness parametrizations of the logarithmic layer cannot adequately capture the large spread observed in the large-eddy simulation data. A generalized frontal area index Λf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Lambda }_f$$\end{document} is introduced that characterizes the vertical distribution of the frontal area in the urban canopy. The vertically distributed stress profiles, which differ significantly per simulation, are shown to roughly collapse upon plotting them against Λf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\Lambda }_f$$\end{document}. The stress distribution representing urban drag can be fitted with a third degree polynomial. The results can be used for more detailed and robust representations of building effects in the development of urban canopy models.

(building plan area per unit plan area) and frontal area index λ f (frontal area per unit plan area), and the vertical extent of buildings is often simply represented by an average building height z H , which may refer to the mean building height, or a (frontal-area-) weighted average height (Grimmond and Oke 1999).
These geometric parameters are used in NWP to describe the aerodynamic effect of buildings, commonly in terms of a logarithmic wind profile with a displacement height z d and roughness length z 0 (e.g., Porson et al. 2010). However, there are large uncertainties associated with estimating the functional relation between the surface geometry and the aerodynamic roughness parameters (Grimmond and Oke 1999;Hagishima et al. 2009;Kanda et al. 2013;Kent et al. 2017). Moreover, real urban canopies have a vertical structure, where building density and frontal area vary with height. Several studies have highlighted that an average building height is an insufficient representation of the real urban environment and that the maximum building height, even in case of an isolated tall building, has a disproportionately large impact on the velocity and building drag of urban areas (Xie et al. 2008;Millward-Hopkins et al. 2011;Kanda et al. 2013;Hertwig et al. 2019). Alternative parametrizations for the logarithmic wind profile incorporate extended morphology statistics such as the maximum height z max and the height standard deviation σ H (Millward-Hopkins et al. 2012;Kanda et al. 2013), but developing robust parametrizations of the aerodynamic roughness parameters remains challenging (Kent et al. 2017).
To provide insight into the aerodynamic effects of buildings, two approaches are typically adopted: (1) a bottom-up approach, which starts at idealized urban morphologies and then gradually adds complexity; (2) a top-down approach, in which real or realistic urban morphologies are considered. Studies pursuing the bottom-up approach have a particular strength in forming an understanding of the fundamental processes governing these flows, and have brought insight on the effects of building cover, building alignment, and wind angle (Coceal et al. 2006;Kanda 2006;Neophytou et al. 2014;Cheng and Porté-Agel 2015;Castro et al. 2017). Case studies pursuing the top-down approach create specific insights for that particular morphology. They show the complexity and rich variety of flow features at a specific location, and highlight the spatial variability of the flow field over urban sites (Carpentieri et al. 2009;Giometto et al. 2016;Hertwig et al. 2019). However, the two are not necessarily compatible, as concepts based on idealized geometries may not match realistic morphologies well (cf. Barlow 2014). And although much progress has been made on the representation of impacts of idealized buildings, the effects of heterogeneity of urban morphologiesdifferently sized and shaped buildings -are not well understood.
There is a need to improve our fundamental understanding of the effect of heterogeneity on flow in urban morphologies of intermediate complexity. That is, realistic urban heterogeneity needs to be introduced without adding too many new degrees of freedom to the parameter space, such that fundamental studies remain tractable. Our first aim is to introduce a tool, the Urban Landscape Generator, that can generate idealized heterogeneous urban environments with identical λ p and λ f . The second aim is to investigate how the drag of these idealized heterogeneous surfaces is distributed in the vertical direction.
In the present study, we investigate the impact of heterogeneous morphologies on urban airflow. Building-resolving large-eddy simulations (LES) are performed for different heterogeneous urban morphologies. The simulations have a similar building density and frontal area index, but vary in complexity with different building heights, plan areas, and street geometries. We explore effects on mean flow structure and turbulence statistics and analyze how they can be quantified with an extended range of building statistics. Aerodynamic roughness parameters for the logarithmic wind profile are estimated, compared with values from the standard parametrizations of Macdonald et al. (1998) and Kanda et al. (2013), and correlated to building statistics. The variation in vertical momentum transport across the simulations is explored, including the separate contributions of turbulent and dispersive momentum fluxes. A generalized frontal area function is proposed that characterizes the vertical structure of urban morphologies. Using this representation of building morphology and the LES data, we derive a parametrization model that describes the vertical distribution of drag inside the urban canopy. The urban drag represents the momentum loss due to buildings and can be incorporated as an additional stress term in the momentum equations.
The paper is organized as follows: Sect. 2 reviews the governing equations and how building effects can be quantified using the spatially-averaged momentum-budget equation. Section 3 presents the numerical set-up and generation of heterogeneous morphologies using the Urban Landscape Generator. Section 4 analyzes the simulation data. First, instantaneous and mean-flow structures are evaluated and logarithmic wind profiles are analyzed; thereafter, vertical momentum transport is explored. Section 5 derives a parametrization for vertical drag distribution based on building morphology. The main conclusions are presented in Sect. 6.

Quantification of Momentum Loss in Urban Environments
Urban settlements modify the flow by forcing it around and over buildings, thus generating complicated flow structures in the wake of the buildings and producing turbulent motion around them (Belcher 2005;Oke et al. 2017). Below we review how buildings affect the mixing and vertical transport processes throughout and above the urban canopy. A momentum-budget equation may be derived by integrating the governing flow equations over the horizontal (x and y) directions, which describes the balance of horizontally-averaged forces acting on the fluid, such as a large-scale pressure gradient accelerating the flow, and a drag force by the building obstacles opposing it.
The horizontal spatial average of a flow quantity φ(t, x, y, z) over the fluid area in the x-y plane, denoted by , is defined by where A a (z) = dA is the fluid area at height z. For a statistically steady-state process, the Reynolds average can be obtained by taking a time average of φ(t, x, y, z) over a time span T longer than the typical flow-fluctuation time scale; the average is defined as A framework for studying flow over spatially heterogeneous surfaces has been established for vegetated canopies (e.g., Raupach and Shaw 1982;Finnigan 2000;Nepf 2012), where flow variables are decomposed to isolate the effects of spatial inhomogeneity by means of a triple decomposition where · is a space-time mean, · represents spatial variations of the time mean, and · represents the turbulent fluctuations. This decomposition is a further refinement of classical Reynolds averaging, since φ = φ + φ , where φ is the Reynolds average.

Momentum-Budget Equation over Spatially Heterogeneous Surfaces
Ignoring buoyancy effects by assuming constant density, the airflow in and above cities is described by the incompressible Navier-Stokes equations where the Einstein summation notation has been used. In the equation above, u i represents the velocity in the i-direction; p =p/ρ 0 + F x x + F y y + gz is the kinematic pressure deviation, wherep is pressure; ρ 0 is a constant density; F i denotes large-scale kinematic pressure forcings in the x and y directions; and g is the standard acceleration due to gravity. We will assume a constant pressure-gradient forcing F i = −ρ −1 0 ∂ P/∂ x i = constant(> 0). Time-averaging, spatial integration over the x-y planes , and division of (4) by the total area of the x-y plane, A T , yield d dz where For simplicity, we have assumed that the flow is in a statistically steady state, which implies that the z-direction variable is the only independent one. The line integral is defined along the building contours ∂ at height z with the normal vector n i pointing into the fluid region.
The volumetric aerodynamic drag f D,i (z) describes the air resistance of the obstacles. The pressure contributions are referred to as form drag [first integral term in (7)], the frictional force related to molecular diffusion when air is moving parallel to the surface of the obstacles is the skin (or viscous) drag [second term in (7)]. Triple decomposition of the velocity variables using (3) and introducing the spatial-average notation (1), gives where Since the spatial average is defined as an intrinsic quantity of the flow variable, a correction term of the area ratio of fluid to total plane area A a (z)/A T arises in Eq. 9. Changes in fluid area with height are substantial in urban areas, and therefore this term cannot be neglected. In the following we denote A a (z) A T · (z) = · c (z), which is also known as superficial (Nikora et al. 2007) or comprehensive (Xie and Fuka 2018) spatial average. This simplifies the notation to The term τ i (z) represents the averaged total kinematic shear stress (turbulent and viscous), which describes the vertical momentum transfer. It consists of the dispersive flux w u i c (z), which represents the vertical transport due to the spatial inhomogeneity in the mean flow, the turbulent momentum flux w u i c (z), and the viscous momentum flux −ν ∂ u i c ∂z (z).

Force Balance for Streamwise Velocity and Constant Pressure Gradient
For an illustration of the force balance, the airflow is assumed only in the x direction with the streamwise velocity component . Equation 10 then simply reads A pressure gradient results in acceleration of airflow from higher pressure towards lower pressure. The averaged forcing F(z) is therefore an acceleration term and a source of momentum. Aerodynamic drag describes air resistance due to the obstacles, which is a force acting opposite to the fluid motion, decreasing the velocity, and therefore f D (z) is a deceleration term and a momentum sink. The change of momentum flux dτ dz (z) describes the transport of momentum from higher altitudes, where the wind speed is greater, driven by large-scale flows and undisturbed from the influence of surface elements, down towards the surface, where surface friction and building drag generate resistance and lower velocities. This is a shear force acting throughout the atmosphere with equal acceleration and deceleration, and therefore a momentum-transport term. The implications of an assumption of steady flow are that the time-averaged forces balance and are unchanging over time. The force balance also provides a quantitative description of the different sublayers in the urban boundary layer (cf. Belcher 2005) ISL: RSL: Note that above the buildings · (z) = · c (z), and that we have ignored viscous terms by assuming high Reynolds number flow. Far above the surface in the inertial sublayer (ISL), the only force acting on the fluid is the constant pressure-gradient force, which is balanced by the constant rate of downward momentum transport, initiated by the surface.
Coming closer towards the buildings in the roughness sublayer (RSL), the mean flow becomes inhomogeneous due to the presence of obstacles below, and spatial patterns emerge in meanflow quantities. Vertical transport initiated by the obstacles on the surface corresponds to a dispersive flux in the force balance that changes with height. Within the bottom urban canopy layer (UCL), all downwards-transported momentum from the upper layers is absorbed by the building drag. An illustration of the balancing forces in the urban sublayers is shown in Fig. 1a.

Cumulative Stresses and Volume Forces
Integrating the streamwise momentum budget (11) from some height z within the urban boundary layer to its top at height h we obtain where the non-negative cumulative stress functions for drag and pressure are given by Figure 1b illustrates the kinematic stresses in Eq. 15. The profile of the constant pressure force τ F (z) does not extend as a straight line due to the volume occupied by buildings within the UCL. The function τ D (z) accumulates the drag force exerted by the buildings from the highest building point throughout the urban canopy layer, reaching the total canopy drag at the ground surface. There, total net momentum sinks and sources are equal, which means that τ D (0) = τ F (0) ≡ τ 0 , where the kinematic surface stress τ 0 represents the effective drag of the urban canopy. The total drag force F D due to the urban canopy is thus given by where F(z) = F x A a (z)/A T and V air is the volume of air inside the domain.

Numerical Model and Set-up
Large-eddy simulations were performed using the Dutch Atmospheric LES model (DALES; Heus et al. 2010), which has recently been extended for the urban environment (uDALES, Tomas et al. 2016;Suter 2018;Grylls et al. 2019Grylls et al. , 2020, providing the capability to model buildings using an immersed boundary method. The eddy-viscosity subgrid scheme follows Vreman (2004) and subgrid-scale (SGS) dynamics close to the building walls and street surface are modelled by logarithmic wall functions (Uno et al. 1995;Suter 2018). The numerical schemes are second-order central differences on a staggered Arakawa C-grid for spatial discretization and an explicit third-order Runge-Kutta time-integration scheme. Nine different urban morphologies S1-S9 are considered, all with a building density of λ p = 0.45-0.46 and frontal area index of λ f = 0.22, but with different numbers, shapes and heights of buildings. The method by which these are generated is discussed in Sect. 3.2. Layouts S1-S3 have uniform building heights. Note that S1 and S2 have approximately the same structure but at different scales. However, the external length scales bounding the flow (in this case the model domain, or the inversion height in the environment) produce different flows, as will be evident later. Layouts S4-S9 have heterogeneous building heights and differ in the street networks. The layouts are illustrated in Fig. 2 and geometry statistics are presented in Table 1. We distinguish three different definitions of mean height: a (plain) average building height z H , an average height where each building height is weighted by its frontal area z H , f , and an average height where buildings are weighted by their plan area z H , p . Weighted standard deviations σ H , f and σ H , p are defined with the respective weighted means. Immediately, these morphologies illustrate the wide variation in urban layout which may be encompassed within near-identical values of λ p and λ f .
The simulation domain size is 480 m × 240 m × 234 m with a horizontal resolution of 2 m. The vertical resolution is constant at 1 m for the lower half of the grid cells and stretched thereafter with a stretching ratio of 1.01. Each layout describes a 120 m × 120 m region and is repeated eight times within the domain (see Fig. 3). This ensures a domain, in which the integral length scales of the turbulence can develop without interference, and has the Table 1 Building statistics for morphologies S1-S9: number of buildings n b , average building height z H , maximum building height z max , average building height weighted by the frontal area z H , f , average building height weighted by the plan area z H , p , building height standard deviation σ H , building height standard deviation weighted by frontal area σ H , f and by plan area σ H , p

Fig. 2
Urban morphologies for simulations S1-S9, generated with the Urban Landscape Generator added benefit of improving convergence of the averages within the RSL. Periodic boundary conditions are applied in the lateral directions, which implies that the simulations capture the atmospheric flow in a larger urban area. The domain-top boundary condition for momentum is free-slip. The numerical simulations are performed under neutral atmospheric stability, and for each of the simulations the domain-average wind speed is kept constant at U = 2 m s −1 . Since the drag will be different for each of the urban morphologies, this implies that the associated pressure-gradient forcing will be different for each simulation.
The simulations are spun up for 10,000 s. An averaging time of 20,000 s is sufficient to produce converged statistics for the momentum budget. For the correct attribution of dispersive stresses, a much longer averaging time is required (Coceal et al. 2006). We found that an averaging time of 80,000 s is sufficiently long to average out slowly evolving mean streamwise circulations. Simulation S3 experiences more persistent slow-transient flow structures, and we therefore increased the averaging time to 100,000 s for this simulation to obtain converged statistics. Later, in Sect. 4.5 some examples of dispersive flux profiles are plotted in context. On the scale of those plots, the averages at 20,000 s would be indistinguishable from those at 80,000 s. Convergence times vary, and other simulations were similar at 60,000 s.

Urban Landscape Generator
To generate a range of idealized heterogeneous urban morphologies we developed the Urban Landscape Generator (ULG) tool, which is a procedural algorithm generating randomized urban landscapes with the morphological parameters λ p and λ f as input. The tool is available as open-source code (Sützl and van Reeuwijk 2020). The ULG tool was designed to gradually introduce more realistic features to idealized urban environments, such as variable building heights and shapes, and more complex street networks. The ULG creates fractal-like street networks by dividing building blocks into smaller blocks until the target building density is reached. Building heights are then chosen to reach the target frontal area density. The resulting buildings are cuboid-shaped and their faces are aligned in the streamwise and cross-stream directions. The ULG has parameters to control the level of randomness in the urban morphology. Layouts with little randomness are similar to the set-ups of generic block simulations (i.e., similar shaped and aligned blocks) as previously used in studies of idealized urban geometries (e.g., Coceal et al. 2006). A more organic city layout can be obtained by increasing the randomization. The resulting layouts are of self-similar structure, this means that they have similar statistical properties at many scales and therefore they are also useful to describe urban landscapes on different scales, for instance they can describe both the layout

Instantaneous Flow Fields
Instantaneous flow fields reveal significant features of the fully-developed turbulent flow around buildings. Figure 4a shows a vertical section of the instantaneous horizontal wind speed u(t, x, y, z) 2 + v(t, x, y, z) 2 for the heterogeneous simulation S4. The influence of the buildings extends well above the maximum height, and particularly the wakes of tall buildings show wind sheltering effects (Xie et al. 2008). The flow displays great spatial variability despite the repeated layouts. Figure 4b shows a horizontal section of the same wind field with velocity components u(x, y), v(x, y). Within the urban canyon, street corridors are important for redistributing the displaced air, with much higher wind speeds than more sheltered areas of the domain. The complex layout of buildings leads to spatially diverse patterns of wind speed and direction. Figure 5 shows profiles of the average wind speed u (z), with height and velocity normalized by the weighted mean height z H , f and mean wind velocity at this height, respectively. The mean wind profiles of simulations with uniform building heights (S1-S3) are evidently distinct from simulations with non-uniform building heights (S4-S9). Wind profiles of the simulations with uniform building heights show a clear separation in within-canopy flow and above-canopy flow with an inflection point around the average building height. Within the canopy, airflow is obstructed by buildings and the velocity profile has a concave shape. Above the canopy, a typical boundary-layer profile develops quickly and the velocity profile changes to a convex shape. Simulations S4-S9 with heterogeneous building heights lack a clear separation in above-and within-canopy flow, instead the velocity profiles show a gradual change from approximately linear (within) to strictly convex (above the canopy); see in particular the difference between the wind profiles of uniform S2 (b) and non-uniform S4 (c) in Fig. 5.

Estimation of Aerodynamic Roughness Parameters
It is commonly assumed that the wind profile in the inertial sublayer can be approximated by a logarithmic profile (Tennekes 1973) where u * is the friction velocity, κ = 0.4 is the von Kármán constant, z d is the displacement height, and z 0 is the aerodynamic roughness length. The parameters z d and z 0 are of interest for weather or climate models, as they are often used to determine wind profile or exchange fluxes for near-surface flows over urban areas. An advantage of LES is the availability of complete vertical profiles of mean wind speed u (z), and friction velocity by u * = √ τ 0 (recall τ 0 is a kinematic stress), which is directly  (17) The friction-velocity values for simulations S1-S9 are shown in Table 2.
Since u * is known a priori, z 0 and z d can be obtained by rewriting (18) as Indeed, given z and u (z), linear regression can be used to estimate parameters z d and z 0 as the slope and intercept of a least-squares fit. The log law is assumed to be valid in the ISL, therefore only wind speeds from within the ISL are appropriate for the linear regression. We tested several fitting regions that are likely to be within the ISL: [2z max , 2z max + 2z H ], [3z max , 3z max + 2z H ], and [4z max , 4z max + 2z H ] for each simulation. However, these ranges may exceed an upper limit to the log law with respect to the boundary layer height, as argued in Jimenez (2004) for example. The resulting logarithmic wind profiles generally match well the wind profile data above the height of 2z max . For simulation layouts with uniform building heights, there is little difference between the different fitting ranges, and the logarithmic wind profiles yield a good approximation to the LES wind profiles even within the RSL. For layouts with non-uniform building heights there can be significant differences in the logarithmic wind profiles, depending on the fitting range. The estimated logarithmic wind profiles vary strongly below the height of 2z max , particularly for simulations with higher friction velocity and larger building-height variation (larger σ H ), showing that the log-law relation breaks down near the surface for strongly heterogeneous building heights. These strong variations in the lower parts of the logarithmic profiles suggest that the region where a log profile may be suitable is at least a length scale of height variation above the buildings. The difference between uniform and non-uniform layouts is illustrated in Fig. 5b and c, where the ranges of logarithmic wind profiles are indicated as grey-shaded areas.
A higher fitting range generally results in greater displacement heights and smaller roughness lengths, and a better agreement between the estimated logarithmic and LES wind profiles higher up. Conversely, a lower range results in lesser displacement heights and higher roughness lengths, with better agreement of the wind profiles closer to the canopy top. A similar trend is reported by Kanda et al. (2013) for their tests of various fitting regions. Ranges of the estimated roughness parameters are presented in Table 2.
An alternative calculation for the displacement height (e.g., used in Kanda et al. 2013;Castro et al. 2017;Giometto et al. 2017) follows a hypothesis by Jackson (1981), that the displacement height corresponds to the centre of moment of the forces acting on the buildings. In other words the displacement height is assumed to be the level of the mean momentum sink. The centre of mass of the volumetric drag function is Rewriting the log law to yields an alternative form for linear regression, where the parameters κ/u * and ln(z 0 ) are slope and intercept respectively of a least squares fit. Using (20) as an estimate for z d and consequently also u * = √ τ 0 , we must assume that κ is the variable parameter. Linear regression with fitting region [3z max , 3z max + 2z H ] yields new values of the von Kármán constant ranging from 0.29-0.37 and consistently lower displacement heights compared to the previous method. A reason for this may be that using the centroid of drag does not consider any displacement effects from dispersive fluxes in the RSL, opposed to the purely statistical fit of z d . The smaller displacement heights are partly compensated for by higher estimates for roughness length, compared to the previous log-law fits. However, a welldefined displacement height as the centre of drag comes at the cost of losing a fixed value for the von Kármán constant, as the variability in κ shows. Values for estimated z d , z 0 , and κ using (21) are listed in Table 2. Figure 6 shows the displacement heights z d and roughness lengths z 0 normalized by the average height z H for simulations S1-S9 estimated with (19) (Macdonald et al. 1998 and other parametrizations, see citations in Kent et al. 2017); however, also z max and σ H are suggested to be important indicators for roughness parameters in the logarithmic law (Millward-Hopkins et al. 2012;Kanda et al. 2013;Kent et al. 2017). In the following, we investigate potential scaling parameters and, more generally, the effects of building morphology on aerodynamic roughness. This is achieved by correlating Fig. 6 Comparison of estimated roughness parameters. Normalized displacement height (a) and roughness length (b) for simulations S1-S9 with fitting range [3z max , 3z max + 2z H ] (markers), and estimated using Macdonald et al. (1998) parametrization (grey dot). Parameter curves for Macdonald et al. (1998) (dashed line) and Kanda et al. (2013) (grey shaded) the urban morphology statistics from Table 1 to z d , z 0 and the friction factor c f = 2τ 0 /U 2 as a non-dimensional and velocity-independent representation of friction. The values of z d and z 0 are taken from the mid-range fitting region [3z max , 3z max + 2z H ]. It should be noted that our aim is to elicit which parameters correlate strongly with c f , z d , and z 0 ; not to provide improved parametrizations of these quantities since these would require a much larger dataset. Figure 7 shows the correlations of a linear regression between the roughness parameters c f , z d , and z 0 , and average height z H , maximum height z max , and z H , z max combined with the standard deviation of heights σ H . Perhaps the most striking aspect of this figure is that average building height z H correlates very weakly with c f , z d , and z 0 (Fig. 7a-c). The weighted average heights z H , f and z H , p also show no significant correlation, with the exception of z H , f correlating to z d (not shown). This finding is consistent with Hagishima et al. (2009), Millward-Hopkins et al. (2011), Zaki et al. (2011, and Kanda et al. (2013). Furthermore, the maximum height z max correlates strongly with both the friction factor and the displacement height, as the only parameter from Table 1. Higher correlations are obtained for combining parameters to a multiple linear regression. Both average height z H combined with the building height standard deviation σ H (g-i), and the maximum height z max combined with height deviation σ H (j-l) show strong correlations for c f , z d , and z 0 .

Comparison and Scaling of Roughness Parameters
In agreement with Kanda et al. (2013), we conclude that z max is more likely to be the relevant scaling length for z d . The maximum height shows a better correlation than the average height parameters z H and z H , p (only the frontal-area weighted average height z H , f has a comparable correlation to the displacement height), and also correlates for c f and z 0 , which is not the case for the average heights. The height structure of buildings clearly affects roughness parameters, which suggests that z max and σ H are important indicators for the length scale of surface variation. The large spread in roughness parameters for similar λ p and λ f (cf. Fig. 6) and the uncertainties of characterizing friction with morphology statistics, demonstrate the challenges of using the logarithmic wind profile to represent momentum transport over heterogeneous urban surfaces. We suggest that it might be more productive to consider the momentum balance directly, from which a velocity profile could be inferred.  a, d, g, j), displacement height z d (b, e, h, k) and roughness length z 0 (c, f, i, l). On the horizontal axes are the urban morphology parameters with the respective best fit coefficients c i . Note that the coefficients differ for each of the subplots. Uniform height simulations S1-S3 are shown with circles, non-uniform S4-S9 with triangles Figure 8a shows the vertical momentum fluxes τ (z), which have all different slopes due to the difference in the pressure-gradient forcings (recall that the volume flux was prescribed). The integral effect of drag is given by the different kinematic surface stresses τ 0 , which vary greatly between the simulations. From simulation S3 with the lowest canopy drag to simulation S4 with the highest canopy drag, the surface stress roughly doubles. This once more highlights that the interaction of urban morphology and the atmospheric boundary layer is not determined solely by the building density and frontal area index.

Vertical Transport
Several physical processes in urban environments initiate vertical momentum transport. Mixing from turbulent eddies is the main contribution to the vertical momentum transport, although dispersive momentum fluxes, which represent the effects of spatial inhomogeneity of the mean flow, also play a role in the roughness sublayer and canopy layer. In simulations S1-S9, both turbulent and total momentum fluxes reach their maximum at the top of the urban canopies. The overall contribution of dispersive momentum fluxes is about 5% of the total momentum fluxes. Subgrid momentum fluxes, which represent unresolved momentum fluxes in regions of high shear, contribute a maximum of 1% of the total momentum fluxes and are only present at heights with large horizontal surfaces (i.e., the ground and building tops). Figure 8b and c show a detailed view of the momentum fluxes of two simulations with uniform (S2) and non-uniform (S4) building heights in the lower levels of the urban boundary layer.
Locally, dispersive fluxes are significant and can make up to 50% of the total fluxes at some heights within the canopy. Except for S1, all dispersive fluxes peak within the canopy layer. Above the canopy, the dispersive fluxes quickly reduce and turbulent fluxes become dominant. Non-zero dispersive fluxes range up to 2-6 times the maximum building height, suggesting a roughness sublayer up to these heights. For uniform building heights (S1-S3) the RSL does not extend further than three building heights. The shape of the dispersive flux varies with the building layout. In simulations with uniform building heights (S1-S3), the dispersive downward momentum transport increases continuously with height up until the canopy top. Heterogeneous building heights (S4-S9) show more variation in the dispersive transport, due to several height levels where the flow is displaced above the building tops and decelerates. Dispersive fluxes within the urban canopies are larger for heterogeneous simulations with larger total canopy drag (S4, S7, S8). The downward dispersive momentum transport peaks typically below the canopy top, and often around the frontal-area weighted average height. Some simulations have re-circulation zones with strong flows in the opposite direction to the driving flow near the bottom surface, so that the dispersive fluxes change sign.

Drag Parametrization
The cumulative drag function τ D (z) describes the accumulation of building drag in the urban canopy layer towards the ground. Figure 9a shows the cumulative drag profiles of all simulations S1-S9. The profiles vary strongly for each of the simulations with no apparent trend. The drag increases quickly at the top of the buildings, which is evident from the profiles of simulations S1-S3. We may further assume that the buildings' surface area higher up contributes more drag, due to the generally higher wind speeds at higher altitudes. In order to measure the effects of the wind-facing building surfaces relative to a reference level, a descriptor of the vertical structure of the building layouts is needed, which will be a generalized version of the frontal area index λ f .
Let L(z, θ) be the total projected building width at height z for wind direction θ . With this definition, the frontal area index is A generalized frontal area index f can be defined as which is essentially the normalized projected frontal area above height z. Clearly, f (0, θ) = λ f and f (z > z max , θ) = 0. In the simulations S1-S9, the wind orientation and the buildings are aligned with the x axis. This means that θ = 0 and the total building width is denotes the width of building B i at height z. Note that for the current set-up B w,i is constant with height. The scaled frontal area shown in Fig. 9b can be used to serve as an alternative height representation of the buildings where ζ = 0 represents the top of the buildings and ζ = 1 represents ground level. Because the height coordinate is scaled as z/z max , the ζ profiles of layouts S1-S3 with uniform buildings heights overlap to a linear function with constant slope. Layouts with non-uniform building heights have piecewise linear functions ζ(z).
The LES data reveal a strong relationship between the cumulative drag and the vertical building structure. Figure 9c shows that plotting τ D (z)/τ 0 as a function of 1 − ζ(z) results in a practically full collapse of all the data. The collapsed data can be fitted reasonably well using a single third-order polynomial of the form with A = 1.88 and B = −3.89. Testing the polynomial regression with different orders shows that the third-order polynomial captures sufficiently well the vertical shape of τ D (z); increasing the polynomial order does not result in significantly more accurate approximations. Note that (25) enables direct evaluation of the cumulative canopy-drag profile provided that the net total canopy drag τ 0 is known. Figure 9 shows the estimated profiles in (d), which closely resemble the data in (a). It shows that a distribution of the urban canopy drag can be reconstructed using only the building layout and an estimation for the total canopy drag. By differentiating (26), the volumetric drag term can be determined as The term L(z)/A F expresses the ratio of the frontal area above height z to the total frontal area, and it shows that the wind-facing area of the buildings play a key role in the production of drag. That the production of drag is not necessarily constant across the urban canopy is represented by s (ζ (z)), a nonlinear function of the building layout ζ(z). Equation 27 is a generalization of the canopy-drag model for uniform buildings of Coceal and Belcher (2004) and Belcher (2005). In Coceal and Belcher's model with N uniform buildings, L(z) = N B w , from which it follows that ζ(z) = 1 − z/z H for z ≤ z H . Belcher (2005) assumes constant canopy drag, in which case s(ζ ) = ζ and therefore s (ζ ) = 1. The volumetric drag is then given by f D (z) = τ 0 /z H with an appropriate estimation for the net total canopy drag. Figure 10 shows two generic examples, one of a single uniform building and one of several heterogeneous buildings, which illustrate the drag distribution throughout the canopy given by Eq. 27. Because the cumulative drag has been fitted with a cubic polynomial function, the volumetric drag profile is quadratic in shape. For a single building, the volumetric drag initially decreases above the ground, then increases nonlinearly towards the top of the building, where most drag is produced. The volumetric drag for multiple buildings with heterogeneous shapes and heights combines the drag of individual buildings in a nonlinear way, with spikes in volumetric drag clearly indicating each building top.

Conclusions
Urban areas are intrinsically heterogeneous. Parametrizations of urban areas often rely on simple characterizations based on building density and frontal area index, implicitly assuming homogeneous surface conditions. In this study, the differences in mean-flow structure and turbulence statistics resulting from idealized heterogeneous morphologies were explored. The urban morphologies were generated using a new tool, the Urban Landscape Generator, which is capable of generating idealized heterogeneous urban morphologies with identical λ p and λ f .
The heterogeneity of an urban site was shown to strongly influence the vertical structure of mean flow and dispersive vertical momentum transport. At equal average wind velocity U , the total surface drag varies as much as almost a factor two between a uniform and strongly heterogeneous morphology. Regression analysis showed that the variations are correlated to the height structure of buildings, particularly to the maximum building height z max and height variation σ H , rather than average height z H .
Developing parametrizations for heterogeneous building morphologies based on the logarithmic wind profile is challenging. First, even with a detailed velocity height-profile and known friction velocity, the estimated displacement height z d and roughness length z 0 for heterogeneous morphologies cannot be uniquely identified, because they depend on the selected fitting range of the data. Whilst this may partially be attributed to the limited domain height used in the simulations, it demonstrates intrinsic uncertainty in estimates for z 0 and z d . Second, our results indicate that current well-known parametrizations of the log law such as Macdonald et al. (1998) and Kanda et al. (2013) cannot sufficiently capture the large spread of roughness-parameter values, even with additional building statistics such as the maximum height z max and height variation σ H .
An exploration of the vertical structure of the surface drag revealed that, at first sight, there is very little that the different simulations have in common. However, the generalized frontal area index f (z) (23), which encapsulates the height distribution of the surface area, was shown to be able to collapse all distributed stress profiles on a single curve. This paved the way to a parametrization of the drag distribution via a third-order polynomial (25). Importantly, this relationship is not linear and the top of the urban canopy, i.e. the highest building in the morphology, produces significantly more drag than the buildings below.
The new drag-distribution model can be used to guide the development of urban canopy models. Numerical weather prediction models may benefit from a distributed-drag approach, because it allows for a more robust and detailed representation of building effects. This is especially important for urban areas with high-rise buildings, i.e. with a large subgrid heterogeneity. However, the distributed-drag model developed here still relies on the friction velocity as an input parameter. The current study investigated the correlations between friction velocity and urban morphology which showed particularly high correlations with the maximum height z max . Extensive simulation/experimental campaigns are necessary to provide improved parametrizations of the friction velocity that are statistically significant.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: Urban Landscape Generator
This appendix describes the technical details of the Urban Landscape Generator (ULG) tool, which is open-source and can be downloaded from https://github.com/bss116/citygenerator (Sützl and van Reeuwijk 2020). The Urban Landscape Generator is a procedural algorithm that builds fractal-like street networks based on urban morphology parameters. The algorithm divides building blocks into smaller blocks until a target plan-area density is reached. Building blocks are then given heights to match the target frontal-area density. A chosen level of randomness influences the street widths, intersection points and building heights. The ULG can generate arbitrarily many urban landscapes with a fixed set-up of λ f and λ p , and gradually increases the complexity of the layouts based on the degree of randomness in the procedural generation.
At the start the Urban Landscape Generator initializes the full domain as first block. The algorithm then chooses a street intersection within the block to subdivide the block into four new blocks. For as long as the sum of the newly created block areas is larger than the target plan-area, the algorithm continues subdividing blocks. On the street level, heterogeneity is introduced by specifying a fractal type, which determines how a block is chosen for subdivision. The 'hierarchical' type chooses the blocks in a set order, such that all blocks on the same level of hierarchy are divided first, before going into a lower hierarchy level. This type produces fairly idealized structures, the simplest of which is a regular grid. The 'random' type selects a random block out of all available blocks. If the selected block is too small for further subdivision, the algorithm selects the next suitable block. This type is designed to better resemble the organic growth of cities.
Two distinct randomness parameters control the street widths and intersections (layout randomness r l ), and the height of blocks (height randomness r h ). The randomness parameter Fig. 11 Fractal types and randomness of the Urban Landscape Generator. a Hierarchical layout, no randomness (r l = 0). b Hierarchical layout with randomness (r l = 0.6). c Random layout with randomness (r l = 0.7) r i , i ∈ {l, h}, defines a probability distribution across a given interval or set of points, from which the ULG draws sample values. No randomness r i = 0 returns a pre-set value from the interval, e.g. the mid-points of the block sides for a new intersection. With randomness r i = 1 the algorithms chooses freely from the given interval. For 0 < r i < 1, the algorithm chooses points such that the probability of obtaining the pre-set value is r i , the remaining probability is uniform across the rest of the interval. For new intersections, the probability of obtaining the mid-points of the block sides is therefore r l , the remaining probability is uniform across the rest of the block sides, obeying certain margins at the edges. A similar distribution is defined over the range of potential street widths. Note that the probability density function is written as a separate module to the code and can therefore easily be modified or replaced e.g. with a p.d.f. based on morphology data.
Realistic values for street widths of various types (e.g. high street, residential street; following, Department for Transport 2007) are defined and chosen according to the level of hierarchy, i.e. how often the algorithm has already subdivided the domain. After the final layout is set, the algorithm adds heights to each of the blocks, until the prescribed frontal area is reached. The height randomness parameter r h defines how a height is selected from a range of possible heights, where a uniform height is added to each of the blocks if r h = 0. Figure  11 illustrates three typical layouts with different configurations. Algorithms 1-3 describe the layout-generation code.  Select the next block on same level of hierarchy. If there is none, pick the first block on the next lower level. return b j function generate new blocks(b j , r l ) # Divides block b j into 4 new blocks b j 1 , b j 2 , b j 3 , b j 4 . Define possible street widths W and standard width W S . Select widths with randomizer(W , W S , r l ). For both block sides L j get mid-points M j . Select a street intersection point i j with randomizer(L j , M j , r l ). Divide block b j at intersection i j into new blocks b j 1 , b j 2 , b j 3 , b j 4 . return b j 1 , b j 2 , b j 3 , b j 4 Algorithm 3 Probability density function function randomizer(I , S, r ) # Samples R from interval I and standard output S with randomness 0 ≤ r ≤ 1. Define uniform distribution U (I ) over interval I . Define probability density function p(x ∈ I ) = (1 − r )S + rU (I ). Draw a random number R from p(x). return R