Turbulence Characteristics Across a Range of Idealized Urban Canopy Geometries

Good representation of turbulence in urban canopy models is necessary for accurate prediction of momentum and scalar distribution in and above urban canopies. To develop and improve turbulence closure schemes for one-dimensional multi-layer urban canopy models, turbulence characteristics are investigated here by analyzing existing large-eddy simulation and direct numerical simulation data. A range of geometries and flow regimes are analyzed that span packing densities of 0.0625 to 0.44, different building array configurations (cubes and cuboids, aligned and staggered arrays, and variable building height), and different incident wind directions (0∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0^\circ $$\end{document} and 45∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$45^\circ $$\end{document} with regards to the building face). Momentum mixing-length profiles share similar characteristics across the range of geometries, making a first-order momentum mixing-length turbulence closure a promising approach. In vegetation canopies turbulence is dominated by mixing-layer eddies of a scale determined by the canopy-top shear length scale. No relationship was found between the depth-averaged momentum mixing length within the canopy and the canopy-top shear length scale in the present study. By careful specification of the intrinsic averaging operator in the canopy, an often-overlooked term that accounts for changes in plan area density with height is included in a first-order momentum mixing-length turbulence closure model. For an array of variable-height buildings, its omission leads to velocity overestimation of up to 17%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$17\%$$\end{document}. Additionally, we observe that the von Kármán coefficient varies between 0.20 and 0.51 across simulations, which is the first time such a range of values has been documented. When driving flow is oblique to the building faces, the ratio of dispersive to turbulent momentum flux is larger than unity in the lower half of the canopy, and wake production becomes significant compared to shear production of turbulent momentum flux. It is probable that dispersive momentum fluxes are more significant than previously thought in real urban settings, where the wind direction is almost always oblique.

Much of the inspiration for urban canopy turbulence and drag parametrization development has come from the vegetation canopy literature. The vegetation canopy-top shear layer is viewed as a plane mixing layer (Raupach et al. 1996), where the drag of the canopy elements reduces the flow velocity below the canopy-top and causes an inflection in the double-averaged velocity profile, which in turn generates Kelvin-Helmholtz waves. The latter themselves undergo secondary instabilities that lead to the formation of large threedimensional eddies. Vegetation canopy elements are usually either very close together (for example in wheat crops) or form a porous mesh (for example the foliage in a rainforest crown), so the mixing-layer eddies generated at canopy top are much larger than the separation of canopy elements (Raupach et al. 1996). These eddies are not local to the canopy elements and are the dominant turbulent motions in the canopy, a point that supports the use of a constant turbulent length scale in vegetation canopy models (Finnigan 2000;Harman and Finnigan 2007;Finnigan et al. 2015).
There are distinct differences between urban and vegetation canopies. The shear generated at the top of urban canopies is more local to the canopy elements than in vegetation canopies (Coceal et al. 2007a). Also, flow separation occurs at the edges of buildings causing vortex shedding and sometimes also turbulent flapping motions at canopy top (Coceal et al. 2007b;Perret and Savory 2013). In current urban canopy models that solve the double-averaged equations, a constant turbulent length scale is often used throughout the entire depth of the canopy (Santiago and Martilli 2010;Schoetter et al. 2020), and the mixing-layer analogy is sometimes used as justification (Hamdi and Masson 2008). To evaluate this assumption, the extent to which canopy top mixing-layer-type turbulence contributes to momentum transport within different urban canopy geometries requires investigation.
Additionally, the solid fraction in urban canopies is typically much larger than in vegetation canopies. The flow has to deviate more around the obstacles, and time-averaged flow patterns exhibit appreciable dispersive transport of momentum, which is distinct from turbulent transport. A recent study of a realistic urban canopy geometry conducted by Giometto et al. (2016) found that the dispersive momentum flux (DMF) can be as significant as the turbulent momentum flux (TMF), which is greater than previously thought. Also, Schmid et al. (2019) argue that, due to the non-negligible solid fraction of the urban canopy, more attention needs to be paid to the formal definition of the spatial average operator. When the intrinsic spatial average (an average only within the fluid volume) is applied to vertical derivatives of flow quantities, a generally overlooked term appears that accounts for changing solid fraction with height (see Sect. 3.2). The importance of including this term in urban canopy turbulence closures for velocity prediction needs to be investigated.
Despite its well-known limitations, a first-order momentum mixing-length turbulence closure approach has been successfully applied in vegetation canopies (Poggi et al. 2004;Harman and Finnigan 2007), and is simple and flexible, requiring no extra prognostic equations. We take inspiration from this pragmatic approach and in Sect. 3.4 formulate a momentum mixinglength closure suitable for urban canopies, which, following Schmid et al. (2019), accounts for the non-negligible urban solid fraction.
There have been various investigations of turbulent and dispersive momentum transport in urban canopies (e.g., Roth 2000;Belcher 2005;Coceal et al. 2007a;Takimoto et al. 2011;Nazarian et al. 2020), but this study is unprecedented in the range of large-eddy simulation (LES) and direct numerical simulation (DNS) datasets available. We take advantage of this rich source of data to: (1) characterize the momentum mixing-length behaviour in the urban surface layer and propose a general momentum mixing-length profile suitable for a wide range of urban flows (Sects. 4.1-4.5), (2) determine the impact on predicted velocity when the term accounting for solid fraction height variation is included within the turbulence closure Note that only a subsection of the domains is shown here (Sect. 4.6), (3) test the mixing-layer analogy for turbulence within urban canopies (Sect. 4.7), and (4) explain the differences in DMF characteristics across geometries to motivate future DMF parametrizations (Sect. 4.8).

Large-Eddy Simulations
The analysis draws primarily on LES datasets described in Nazarian et al. (2020). A suite of simulations were carried out for idealized urban-like staggered and aligned arrays of cubes of uniform height. The flow was oriented perpendicular to the array, and the stratification was neutral. The cases we study have packing densities (planar area fraction occupied by obstacles) λ p = 0.0625, 0.11, 0.16, 0.25, 0.35, 0.44 for staggered cubes, and λ p = 0.0625, 0.44 for aligned cubes. A summary of the LES datasets is given in Table 1. Schematics in Fig. 1a, b are illustrative of the aligned and staggered geometries respectively. The simulations were performed using the Parallelized Large-Eddy Simulation Model (PALM) which solves the non-hydrostatic incompressible Boussinesq equations (Raasch and Schröter 2001;Letzel et al. 2008;Maronga et al. 2015). The subgrid turbulence parametrization is based on a 1.5-order closure (Deardorff 1980) that involves solving the subgrid scale (SGS) turbulence kinetic energy (TKE) equation.
The flow is driven by a constant pressure gradient force F p = ρu 2 τ /H , where ρ is the density of dry air (taken to be 1 kg m −3 ), u τ ≈ 0.21 m s −1 is the wall friction velocity, H = 7.4h is the domain height, and h is the mean average height of the obstacles. The grid length in the x, y, and z directions, xyz , is equal to h/32, except above z = 4h where z slowly increases.
The Reynolds number of the flow is Re = U h/ν ≈ 10 6 , which is in the fully rough regime. Here U is the velocity at the top of the domain and ν is the molecular viscosity of air. It has been demonstrated using LES that flow through cube arrays only has a weak dependency on Re between 5 × 10 3 and 5 × 10 6 (Xie and Castro 2006). For momentum the boundary conditions are no-slip at the ground and building surfaces, free-slip at the domain top, and periodic in the horizontal. Spin-up time is 125h/u τ , and data are output every 20 timesteps, which corresponds to a time interval of 40 s or 0.46h/u τ . The time-averaging interval is 250h/u τ . The time-averaged data are also ensemble averaged over repeating units to effectively increase the averaging time.  Nazarian et al. (2020), Branford et al. (2011), andCastro et al. (2017) respectively. Under "Case" numbers correspond to λ p , S and A denote staggered and aligned respectively, and cb and cd denote cubes and cuboids respectively

Direct Numerical Simulations
Several DNS datasets are included in the analysis as summarized in Table 1. Unlike in the LES approach, there is no subgrid parametrization in DNS and turbulent scales are resolved through the majority of the dissipation spectrum. As seen in Table 1, a compromise is made by reducing Re due to computational demands. Branford et al. (2011) conducted simulations of flow through an aligned array of λ p = 0.25 cubes (see Fig. 1a for a schematic). Castro et al. (2017) documented simulations for an aligned array of λ p = 0.33 cuboids, with cuboid length-to-width ratio of 1/2 and height-to-length ratio of 1, where the long faces were orientated perpendicular to the flow (see Fig. 1c for a schematic). Both studies had neutral stratification, and forcing flow directions of 0 • and 45 • to the normal of the building faces. Throughout the entirety of each domain xyz = h/32. All time-averaged data are also ensemble-averaged over repeating units.

Double-Averaging Theory
In this section the double-averaged momentum equation and first-order momentum mixinglength closure parametrization of the TMF term are given. The definitions of the spatial average operator and spatial averaging theorem are first presented, since they are necessary in deriving budget equations and parametrizations. The intrinsic first-order momentum mixinglength closure is used in later sections to investigate its use as an urban canopy turbulence closure model and to explore the turbulence behaviour across different canopy geometries. Equations are given using the intrinsic average and comprehensive average (where the spatial average is done over the fluid and solid volume). Following Schmid et al. (2019), the impor-tance of representing λ p height variation in the turbulence closure for predicted values of velocity is investigated in the context of both horizontal averaging approaches (see Sect. 4.6).

Intrinsic and Comprehensive Averaging
In NWP models flow through canopies is often modelled at grid lengths much larger than the typical turbulent eddies and obstacle wakes. An increasingly common approach in such circumstances is to apply to the instantaneous equations temporal averaging over a longer interval than the time scale of the slowest eddies and spatial averaging over a length scale greater than the largest spatial deviations in the flow (Raupach and Shaw 1982;Nikora et al. 2007;Martilli and Santiago 2007). A one-dimensional representation is taken, assuming quasi-horizontal homogeneity of the flow and canopy geometry within vertically thin horizontal-averaging slabs. Terms arise representing the effects of turbulence, form drag, and persistent flow structures on vertical transport of momentum and scalars, and can subsequently be parametrized.
In vegetation canopies the obstacles are considered to occupy a negligible fraction of the canopy, and the comprehensive average (also known as the superficial or extrinsic average) and intrinsic average are equivalent (Finnigan 2000). In porous media hydrodynamics, as in urban canopies, the volume fraction occupied by obstacles can be more significant. At any given height it is proportional to λ p and it is well known that vertical changes in λ p should be treated when applying the spatial average operator (Nikora et al. 2007).
The use of the operators φ , φ,φ, and φ , on a variable φ defined only within the fluid region, denote its horizontal spatial average, time average, dispersive fluctuations (φ = φ − φ ), and turbulent fluctuations (φ = φ − φ ), respectively. The intrinsic and comprehensive spatial operators are defined by Whitaker (1999) as and respectively. Here x, y, and z are spatial coordinates with z in the vertical, and t is time: is the fluid fraction: the ratio of fluid volume within a thin horizontalaveraging slab, V f , to the total volume of the slab including the solid volume, V . The relationship between the comprehensive (subscript C) and intrinsic (subscript I) average of a quantity φ is given by The application of the comprehensive operator twice results in φ C C = φ C , whereas for the intrinsic operator φ I I = φ I . This means factors of appear in parametrizations of comprehensive spatially-averaged quantities when they are functions of spatially-averaged parameters.
The comprehensive averaging approach has a disadvantage when used to compute dispersive fluctuations (φ) of a variable φ, since interpretation of quantities involving dispersive fluctuations can become more difficult (Schmid et al. 2019). When φ C is non-zero,φ is non-zero in the building region. The dispersive flux due to the correlation of two such variables can then also be non-zero in the building region, which is physically unintuitive. In the context of parametrizing vertical dispersive momentum (and scalar) fluxes the issue is not important. They involvew, and w C = 0 due to mass continuity, so they are zero in the building region. However, some terms in higher-order turbulence parametrization schemes (e.g., wake production of TKE) involve dispersive fluctuations, and in the comprehensive averaging framework would be non-zero in the building region, thus complicating their parametrization.
There has been some debate in the recent literature on whether the intrinsic or comprehensive spatial average should be used in urban canopy modelling (Xie and Fuka 2018;Schmid et al. 2019;Sützl et al. 2021). The double-averaged momentum equation (and budget equations in general) derived using either method are equally valid if obtained correctly. Equivalently to solving the intrinsic double-averaged momentum equation, one could solve the comprehensive double-averaged momentum equation using parameters that have also been derived using the comprehensive average, and then convert the predicted comprehensive velocity to the intrinsic velocity via division by (z). Intrinsic average values are used here since they are more directly representative of the fluid and hence easier to interpret physically. However, neither spatial averaging approach is being advocated in a parametrization context.

Spatial Averaging Theorem
The relation between averages of spatial derivatives and spatial derivatives of averages is given by the spatial averaging theorem (Whitaker 1999). For intrinsic averages and for comprehensive averages where x i (x 1 = x, x 2 = y, x 3 = z) is spatial location, S int is the interface between the fluid and solid regions, and n i is the unit vector normal to the interface directed from fluid to solid. The first term on the right-hand side (r.h.s.) of Eq. 3 is the spatial derivative of the intrinsic average φ. The second and third terms account for variations in (i.e., changes in averaging volume) and discontinuities in flow properties over the surface of obstacles respectively. The second term does not appear in Eq. 4 since the comprehensive averaging region does not change with height. The first and second terms in the comprehensive spatial averaging theorem have equivalent interpretations to the first and third terms in the intrinsic spatial averaging theorem.
For horizontal derivatives (i = 1, 2), the second term on the r.h.s. of Eq. 3 is zero, since the canopy is assumed horizontally homogeneous within the averaging region, so that is only a function of z. The third term also vanishes for horizontal derivatives if φ is constant at the fluid-solid interface, as is the case for velocity due to the no-slip boundary condition. However, for example, pressure is variable at the interface, tending to be larger at the windward than leeward obstacle faces. It is the discontinuity in the pressure field over obstacles that gives rise to form drag.
For vertical derivatives (i = 3), the second term on the r.h.s. of Eq. 3 is non-zero everywhere that (z) varies and φ I = 0. The third term on the r.h.s. is only zero for vertical derivatives of quantities like velocity, where φ = 0 at the vertical facing interfaces. Unlike when calculating horizontal derivatives for the x−z and y−z plane interfaces, when calculating vertical derivatives for the x−y plane interfaces there are no adjacent fluid-solid interfaces. Therefore, the third term on the r.h.s. cannot have cancellation of adjacent interfaces, and is non-zero so long as φ > 0.

Double-Averaged Momentum Equation
The approach of solving the double-averaged momentum equations in canopies was first developed by Wilson and Shaw (1977) and Raupach and Shaw (1982). The influence of changing solid fraction was treated as negligible in the context of vegetation canopies, where comprehensive and intrinsic averaged values are approximately equal.
The double-averaged momentum equation can be obtained by first applying the timeaveraging operator to the instantaneous momentum equation to obtain the time-averaged momentum equation, then applying the spatial averaging operator, taking particular care with use of the spatial averaging theorem. The intrinsic and comprehensive double-averaged momentum equations for a statistically stationary, horizontally homogeneous flow at height z in the canopy are given by (Schmid et al. 2019) and respectively. Here ρ is the density of dry air, p is pressure, F p is a volumetric body force driving the flow and u i are velocity components (u 1 = u, u 2 = v, u 3 = w). It has been assumed that atmospheric stability is neutral, that the Coriolis force and molecular transport are negligible, and that the viscous drag can be neglected since it is usually only a few percent of the total drag, as discussed in Leonardi and Castro (2010). The first and second terms on the left-hand sides (l.h.s.) of Eqs. 5 and 6 are the gradients of the TMF and DMF respectively. They are generally negative within the canopy, while above the canopy the TMF gradient is positive and the DMF gradient can be positive or negative but becomes negligible in the inertial sublayer.
The third and fourth terms on the left-hand side of Eq. 5, which are not usually included in urban canopy literature, appear due to the second term on the r.h.s. of Eq. 3. When (z) increases with height, vertical momentum transport occurs through a larger cross-sectional area, and consequently the intrinsic averaged TMF and DMF decrease. The extra terms act to increase the contribution of the intrinsic averaged vertical TMF and DMF gradients to the double-average momentum equation, when momentum transport occurs through increasing cross-sectional area with height.
The first term on the r.h.s. of Eqs. 5 and 6 is a body force applied uniformly throughout the fluid and drives the flow in the computational fluid dynamics (CFD) simulations. It is noted here that a mean streamwise pressure gradient is not actually present within the fluid in the CFD simulations, so F p does not arise from spatial averaging of the time-averaged pressure gradient term, −(1/ρ)∂ p/∂ x i , and should be included as a separate term in the instantaneous momentum equation when deriving Eq. 5. The second term on the r.h.s. of Eqs. 5 and 6 is the form drag exerted by the obstacles and is the sink of momentum in the canopy. It appears after spatial averaging of the time-averaged pressure gradient term.

Parametrization of Momentum Fluxes
K -theory is often used to parametrize the TMF term in Eqs. 5 and 6. The momentum eddydiffusivity K m,i can be represented using a first-order momentum mixing-length-closure approach, where l m,i is the momentum mixing length. It is well known that Prandtl-mixinglength-type turbulence closures break down when applied to the local flow in urban canopies, for example, due to counter-gradient momentum transport in building wakes. The momentum mixing-length closure is therefore formulated in terms of horizontally-averaged properties of the flow, so that for intrinsic averages and for comprehensive averages K m,i and l m,i are parameters representing the global properties of the flow at height z. The modulus operators ensure K m,i is positive and that turbulent diffusion is only downgradient. The terms in Eq. 7 occur due to the second term on the r.h.s. of the intrinsic spatial averaging theorem; −3 occurs in Eq. 8 to balance the number of comprehensive spatial averaging operator applications on the l.h.s. and r.h.s., accounting for the zeros introduced into the horizontal average by the building volume. Schmid et al. (2019) derived the same intrinsic and comprehensive mixing-length closures (although without the modulus operator). However, they started from application of the horizontal averaging operator to parameters representing the local flow. Following use of their Eqs. 10 and 11, the derivation involved neglecting correlations between spatial fluctuations of l m (l m ), and betweenl m and spatial fluctuations of ∂ u /∂z. This implies that our approach of starting from global properties of the flow assumes spatial fluctuations are small or terms involving them largely cancel.

Application to Urban Canopy Geometry
The i = 1 and i = 2 components of double-averaged flow properties (e.g., u and u w , and v and v w ) in general take different values, except in the special case when the canopy has lateral symmetry and the driving flow has equal components in the x and y axes. As an approximation one could use the same turbulent parameters in both axes, and derive them with the x axis aligned with the driving flow direction. However, given urban morphology information such as the frontal area variation with driving flow angle, it is conceivable that parametrizations might be developed that account for flow direction relative to the two axes. We therefore analyze separately the i = 1 and i = 2 TMF components.
In real-world urban canopies the double-averaged flow properties are never exactly the same along the two axes, but the 0 • driving flow and B-25-A-45 • simulations are special cases for which subscript i can be dropped from l m and K m . In the 0 • driving-flow simulations turbulence is not the same in the two axes but v w ≈ 0 so the i = 2 component of the turbulent momentum flux need not be represented. For B-25-A-45 • the canopy is symmetric about the driving-flow direction and the driving flow has equal components in the x and y axes, so that l m,1 = l m,2 and K m,1 = K m,2 . In this study the i = 1 and i = 2 components only need to be retained for C-33-A-45 • where the geometry is asymmetric about the streamwise axis.
A pragmatic approach to parametrizing the DMF is to include it with the TMF on the l.h.s. of Eq. 7 (Simón-Moral et al. 2017;Nazarian et al. 2020). Horizontally-averaged DMF can be negative or positive at different heights within the canopy, but when summed with the horizontally-averaged TMF the total tends to be negative. The DMF can thereby be parametrized using K -theory but avoiding the unphysical scenario where K m becomes negative, which would occur if the DMF was parametrized separately.
K -theory assumes that turbulent fluctuations are related to local gradients. Although it seems plausible that downward transport of momentum from higher in the canopy leads to increased velocity lower in the canopy, it has not been demonstrated that there is a strong relation between DMF and local velocity gradients in urban canopies. Time-averaged motions might largely be determined by non-local forcings. For example recirculations within building wakes often span the entire depth of the canopy, and are likely largely driven by strong shear at roof level behind the buildings, rather than local velocity gradients within the wake. An analogy can be drawn from two-dimensional street canyon literature where it is well known that recirculation within the canyon can be modelled as solid body rotation (Caton et al. 2003). Recirculation is driven by the flow above the canyon and occurs due to the requirements that tangential stresses be continuous at cavity top, and mass within the cavity be conserved. Flows in urban canopies consisting of cuboidal buildings are more three-dimensional in nature, and streamlines can diverge in the spanwise direction within building wakes, but it seems reasonable to assume that canopy recirculations share some properties with two-dimensional street canyons.
Since the aim is for the first-order turbulence parametrization presented here to be applicable to multi-layer urban canopy models in general and to analyze the turbulent flow properties using l m , it is chosen not to include the DMF in the turbulence closure; l m can be interpreted as a length scale that describes only the turbulence. The horizontal-average and three-dimensional characteristics of the DMF are investigated across geometries in Sect. 4.8, as a first step towards developing separate parametrizations of DMF.

Mixing Length Within the Canopy
It is investigated whether l m I (calculated using Eq. 7) shares similar characteristics within the canopy across different idealized urban geometries and driving-flow directions. Figures 2 and 3 correspond to the 0 • and 45 • driving flow datasets, respectively. In both figures, plots a-e correspond to intrinsic average profiles of velocity, TMF, DMF, the ratio of DMF to TMF and momentum mixing length, respectively. Figure 2f shows l m I normalized by the maximum mixing length l m max I within the canopy. The DMF profiles are discussed later in Sect. 4.8.

0 • Driving Flow Mixing Length Characteristics
There is some collapse in the 0 • simulation l m I profiles from occurs at z/h ≈ 0.6, since the minimum gradient of double-averaged velocity tends to occur near the middle of the canopy (Fig. 2a), but the magnitude of TMF increases with height within the canopy (Fig. 2b). This description of l m I is consistent with Coceal et al. (2006) for a λ p = 0.25 staggered cube dataset and is qualitatively similar to the turbulent length scales presented in Fig. 9a of Nazarian et al. (2020).
From Fig. 2e it can be seen that l m max I does not vary monotonically with λ p for either the aligned or staggered geometries. This has also been observed for staggered geometries when turbulent length scales are calculated with the sum of DMF and TMF, as presented by Nazarian et al. (2020); l m higher within the canopy than the 0 • staggered simulations, since the double-averaged velocity gradient (or shear) within the canopy does not increase significantly until very close to canopy top in the 0 • aligned simulations. For the same λ p , the aligned geometries have larger double-averaged velocity and TMF than staggered geometries. The aligned geometries have unobstructed channelling regions (see Fig. 1a) where flow is fast and there is large shear, which promotes shear production of TMF.
The magnitude of TMF generally decreases with increasing λ p throughout the canopy as seen in Fig. 2b, particularly for the staggered geometries. A possible explanation is that increasing λ p results in less mean kinetic energy penetrating into the canopy, particularly for the staggered arrays where there is more flow obstruction, and therefore less conversion into TKE. TMF and TKE have broadly similar dependence on λ p and aligned versus staggered canopy geometry as seen from Fig. 6 of Nazarian et al. (2020).
Between z/h = 0.34−0.44 the N-44-S simulation TMF just becomes positive so that flow is counter-gradient, and l m I cannot be computed. A first-order momentum mixing-lengthclosure approach is less suitable for very high-density idealized staggered arrays. Future investigations are required to explore whether this is the case in more realistic staggered type urban geometries, where there is more horizontal randomness in positioning of buildings.
The The contribution from such regions to the spatial average results in positive overall double-averaged velocity near the surface. The gradient of double-averaged velocity between the surface and the first grid point above the surface is negative in the N-25-S, N-35-S and N-44-S simulations, and along with TMF being positive above the surface (Fig. 2b), means that l m I at the first grid point above the surface cannot be computed, since transport is counter-gradient.

Impact of Flow Incidence Angle on Mixing-Length Characteristics
The 45 • driving-flow simulations in Fig. 3e have larger l m I in the canopy than their corresponding 0 • driving-flow simulations in Fig. 2e. Figure 3e shows that l m max I is largest for B-25-A-45 • , and is equal to 2.2h at z/h = 0.25 (off the axis). This is unphysical considering that eddies cannot be so large due to the limiting effect of the ground, and the largest eddies are not expected to exceed h in any case. l m max I occurs at z/h = 0.31 and z/h = 0.75 for l m,1 I and l m,2 I respectively in the C-33-A-45 • simulation. Since l m,1 I and l m,2 I are parameters describing the turbulence in the x and y directions, it is evident that turbulence characteristics can be significantly different in the two axes, when the geometry of the urban canopy is not diagonally symmetric.
The 45 • driving flow l m I are large because the double-averaged velocity profiles (Fig. 3a) have small gradients (except near the ground and canopy top), but significant TMF (Fig. 3b). There are two possible explanations for significant TMF despite small shear (i) there is large transport of TMF down into the canopy (see Finnigan (2000), Eq. 2.11 for a definition) and (ii) there is large wake production of TMF. Either way, this indicates that a first-order momentum mixing-length-closure approach may be less appropriate for oblique flows, since the approach assumes TMF is related to (i) the local flow (implying TMF has not been transported) and (ii) vertical velocity gradients (but wake production involves vertical and horizontal gradients of spatial velocity fluctuations).
Shear production of TMF is given by The third term on the r.h.s. is a new term that does not appear in Eq. 2.11 of Finnigan (2000), and arises from accounting for changes in λ p with height, when applying the intrinsic spatial averaging operator to the time-averaged TMF budget equation. For uniform-height canopies it is only non-zero at canopy top. Wake production of TMF is due to correlations between spatial deviations in turbulent stresses and gradients of spatial deviations in velocity, typically caused by building wakes and deviation of the flow around buildings, and is given by Plotted in Fig. 4a, b are the shear and wake production of TMF respectively, and Fig. 4c shows their ratio (P w /P s ). The production of u w I is plotted for the 0 • driving flows, the production of u w I and v w I for C-33-A-45 • , and the average of the two components for B-25-A-45 • . All geometries have a large P s peak at canopy top, as expected due to the canopy-top shear layer. The staggered geometries and N-0625-A have a small P s peak near the surface because the double-averaged velocity gradient is large there, and is associated with flow reversal in the denser staggered cases. Between z/h = 0.2 − 0.8 the aligned geometries tend to have larger P s than the staggered geometries, particularly the denser ones. In general P w is positive near the ground (i.e., acts as a TMF sink), becomes negative and increases in magnitude with height within the canopy, then changes sign and has a positive peak at canopy top, before becoming negligible in the inertial sublayer. The largest outlier is N-44-S, which has positive P w in the middle of the canopy.
In vegetation canopies it is generally assumed that P s dominates P w (Raupach et al. 1986). However, from Fig. 4c it can be seen that this is not always the case for urban canopies, particularly for the 45 • driving flows, and away from the canopy top and ground where P s is smaller. The 45 • driving flows tend to have large DMF (as discussed in Sect. 4.8.2). It is unlikely that P w and DMF are both large by chance, since they are both associated with spatial fluctuations in the flow. Also, a key assumption in the arguments presented by Raupach and Shaw (1982) and Raupach et al. (1986) for P s dominating P w , is that terms involving spatial deviations in the time-mean flow are negligible. This is often not true for urban canopies (as demonstrated for the case of DMF in Sect. 4.8.2).
In addition to the 45 • driving-flow simulations, N-44-S and C-33-A-0 • have comparable magnitude P w and P s near the middle of the canopy. Like the 45 • driving-flow simulations (Fig. 3d), N-44-S has significant DMF compared to TMF (Fig. 2d) at heights where P w /P s is largest. C-33-A-0 • is the outlier in that it has small DMF (Fig. 2d) compared to TMF despite having large P w . This shows that large P w is not always related to large DMF. As demonstrated in Sect. 4.8.1, urban canopies often have large local dispersive velocity correlations, whose sign may partially cancel upon horizontal averaging.
It has been demonstrated that wake production can have a large contribution to TMF, particularly when the flow is oblique to the building faces. Investigation of TMF transport is out of scope here.

Mixing Length at Canopy Top
The behaviour of l m I at z/h = 1 can be most easily understood in terms of comprehensively averaged variables. Rearranging Eq. 8 for l m,i 2 C , and substituting into l m,i At z/h = 1 for a uniform-height canopy with perfectly flat roofs u i C and u i w C are continuous, and and ∂ u i C /∂z are finite but discontinuous. ∂ u i C /∂z increases at z/h = 1 due to the sudden contribution of shear along the roofs, as opposed to beneath the roofs where there is no shear since u i is not defined within the buildings. The discontinuities in and ∂ u i C /∂z do not cancel, so that l m,i I is also discontinuous at z/h = 1. Parametrizing l m (h) I for situations with flat-roofed buildings would be challenging. Rather than analyze l m I at h, we instead choose to focus on values just above and below. In practice roofs are not perfectly flat and there is a variety of building heights, so that l m (h) I might reasonably be approximated as a continuous function. Plotted in Fig. 5a and b is l m I on the grid levels just above and below canopy top respectively; l m I is small near canopy top since there are large velocity gradients. l m (1.03h) I and l m (0.97h) I are between 2 and 10 times smaller than l m max I . Also, l m (1.03h) I is approximately a factor of 2 smaller than l m (0.97h) I , because l m (1.03h) I has velocity gradient contributions from just above the roofs where there is large shear. Values of l m (1.03h) I and l m (0.97h) I tend to increase faster than linearly with decreasing λ p . It is expected that l m increases near canopy top with decreasing λ p since there are fewer buildings generating canopy-top shear. In the limit of an infinitely sparse canopy, flow should resemble that over a smooth wall, so that l m ≈ κh, where κ is the von Kármán coefficient. In the limit of a very dense canopy, flow should also resemble that over a smooth wall, except displaced by h so that l m ≈ κ(z − h), and l m (h) ≈ 0. Although neither limit is reached in the geometries examined here, it provides some explanation for the l m behaviour near canopy top.

Mixing Length Above Canopy Top
The dispersive momentum flux occurs due to spatial deviations in the mean flow, and where it falls to zero can be used as a definition of the roughness sublayer height (Coceal et al. 2007c). As seen in Figs. 2d and 3d, by z/h = 1.5 DMF is generally negligible (apart from N-44-A and C-33-A-0 • which is likely due to insufficient averaging time); l m I tends towards an approximately linear increase by z/h = 1.5, which is typical of inertial sublayer flow where l m I ≈ κ(z − d). By inspecting l m I , a linear portion was found to exist for z/h = 1.5−2.3 across simulations (not shown), and was identified as the inertial sublayer. Linear regression of l m I against z was performed in the inertial sublayer, where the gradient corresponds to κ. The κ results are plotted in Fig. 6. The N-44-A and C-33-A-0 • simulations are excluded since they were not linear by z/h = 1.5, due to large DMF. Parameter κ varies between 0.20 and 0.51, showing no clear relationship with canopy geometry, and has many values far from the classical value of 0.41. Alternative inertial sublayer height definitions were tested to investigate the sensitivity of κ. There are several studies which suggest that κ is not a universal constant in rough wall flows, as discussed by Leonardi and Castro (2010). For example, Frenzen and Vogel (1995) argue that roughness influences κ through the ratio of TKE production and dissipation in the inertial sublayer. For DNS of flow over packed beds (Breugem et al. 2006) and LES of a staggered array of cubes (Claus et al. 2012), values as low as 0.23 and ≈ 0.3 have been observed, respectively, by logarithmic fitting to simulated velocity in the inertial sublayer. However, to the authors' knowledge this is the first set of results to demonstrate such a range of κ values across a range of geometries in the urban canopy literature.

Scale Separation
Jiménez (2004) suggests that two conditions must be satisfied so that a region exists where the only relevant turbulent length scale is the distance from the wall and the velocity profile is truly logarithmic. First δ/h > 40-80, where δ is the boundary-layer depth, so that the flow is far enough from the surface that inner-layer turbulence associated with roughness elements does not interact with it. Second z/δ < 0.15, so that the flow is close enough to the wall to not interact with the outer-layer turbulence of scale δ. The first condition may not be met in urban areas under near-neutral atmospheric stability since the atmospheric boundary layer over deeper urban canopies may not grow to δ/h > 40-80. This might explain some of the large variability in relations for z 0 and d which are often found using experimental datasets (e.g., Kent et al. 2017).
In CFD half-channel flow investigations domain height H places a limit on δ. The first condition is not met in the majority of the simulations presented here where H /h = 8. The second condition becomes z/H < 0.15, which is also not met since the region defined earlier as the inertial sublayer (z/h = 1.5-2.3) corresponds to z/H = 0.19-0.29 for H /h = 8. However, N-35-S and C-33-A-45 • have H /h = 12, so that the region defined earlier as the inertial sublayer corresponds to z/H = 0.125-0.19. Roughly the lower half of the region is then close enough to the wall for the outer-scale turbulence to not interact with the flow according to Jiménez (2004). It is therefore interesting that N-35-S, the C-33-A-45 • x-component and the C-33-A-45 • y-component have κ values of 0.38, 0.39 and 0.5 respectively, which are generally larger than the other values, and, with the exception of the C-33-A-45 • y-component, κ ≈ 0.4. This provides tentative evidence that H used in urban canopy modelling can place limitations on δ, and lead to inadequate separation of the innerand outer-scales in the inertial sublayer resulting in modification of the turbulence.
That H might influence flow in the inertial sublayer does not affect the validity of results presented within the canopy. Turbulence statistics within the canopy have been demonstrated to converge with H in the DNS studies of Coceal et al. (2006) and Coceal et al. (2007a) for H /h values of up to 6 and 8 respectively. Also, large-eddy simulations with H /h ≈ 8 have been shown to agree well with experiments Castro 2006, 2009;Nazarian and Kleissl 2016).

Influence of Vertical Turbulent Momentum Flux Variation
Replacing l m I with κ(z − d) in Eq. 7 and assuming that momentum transport is downgradient, after rearranging one can obtain where τ i = − u i w I . In the κ analysis presented in the previous section, the exact value of τ i derived from the CFD simulations was used in Eq. 12, following Coceal et al. (2006). If instead τ i is taken to be constant then Eq. 12 has the form of the MOST flux-gradient relationship under neutral conditions. Its integral from z 0 + d to z results in the standard logarithmic law (Garratt 1994), where z 0 is the roughness length. When τ i is taken to be constant it is not obvious what value it should take. In the CFD simulations there is a body force applied uniformly throughout the fluid, which results in τ i decreasing linearly with height above the roughness sublayer. A linear fit would equal u 2 τ,i if extrapolated to the surface, where u τ,i = u τ cos θ and θ is the angle between the x axis and the driving-flow direction. Therefore, above the roughness sublayer τ 1/2 i = u τ,i √ 1 − z/H . The effective ground experienced by the flow in the inertial sublayer is d. Assuming that the TMF is constant in the surface layer (as is done in derivation of the logarithmic law), one might then take τ (Leonardi and Castro 2010). Taking this approach κ was estimated from the gradient of the linear regression between u τ,i √ 1 − d/H /(∂ u i I /∂z) and z. The fits were similarly good in the region z/h = 1.5-2.3 across simulations when compared to the linear regressions in the previous analysis. The κ values are plotted in Fig. 6b and vary widely between 0.27 − 0.50. When performing fitting analysis to obtain d and z 0 in the inertial sublayer, neither the approach of taking the exact value of τ i or constant τ i is advocated here as being correct. We instead note that in both approaches a large range of κ are found, and taking a fixed value in any fitting analysis would have a large influence on the values of d and z 0 obtained. If the domain height used in CFD simulations is increased this could have the twin benefits of improving scale separation and reducing τ i variation in the inertial sublayer.

General Characteristics of Mixing Length in Urban Canopies
The general features of the l m I profiles for uniform-height, flat-roofed geometries are summarized in a schematic in Fig. 7 and consist of: an increase from zero at the surface to a maximum around the middle of the canopy, a decrease until canopy top (where there is a discontinuity), and a rapid increase above z/h = 1 before transitioning to a linear increase by approximately z/h = 1.5. The vertical extent of the arrow reflects the fact that the aligned simulations have l m max I nearer canopy top with increasing λ p , and that l m max I can be below the middle of the canopy for 45 • driving-flow simulations; l m max I /h ≈ 0.3 is a reasonable approximation apart from the 45 • simulations which have much larger values. Xie et al. (2008) conducted an LES of flow through λ p = 0.25 staggered cuboids with equal length and width, and varying heights. The horizontal layout was the same as that in Fig. 1b. Within each repeating unit there were 16 cuboids (1 × 0.27h, 3× 0.63h, 7× 0.98h, 4× 1.33h and 1×1.68h). The dataset offers more realistic variability in (z) than uniform-height roofs. The influence of variable building height on the flow and the general applicability of the l m I schematic in Fig. 7 to urban canopies is investigated.

Mixing-Length Characteristics of a Variable Height Building Array
Plotted in Fig. 8a, b are the (z) and l m I profiles, respectively. Below h the momentum mixing-length profile for the variable height array is similar to the schematic for the uniformbuilding-height arrays (Fig. 7); l m max I occurs at z/h = 0.49, which is almost identical to z/h = 0.50, the height at which l m max I occurs for N-25-S, a case with identical layout The variable height array l m (h) I /h was equal to 0.14, and can be compared with the N-25-S simulation at grid points just below and above z/h = 1, where l m I /h = 0.09 and 0.05, respectively. Smaller values occur for N-25-S because there is a larger proportion of roofs at h. Flow interacts with the roofs producing large double-averaged velocity gradients (i.e., shear) and small l m (h) I /h.
There are two less pronounced l m I minima above z/h = 1 compared to the minima at z/h = 1 in Fig. 8, and are a result of shear generated at the tops of the 1.33h and 1.68h buildings. The peaks are less pronounced because there is a smaller number of height 1.33h and 1.68h buildings than height h buildings. l m I minima do not occur at the tops of the 0.27h and 0.63h buildings, despite them occurring in similar proportion to the 1.33h and 1.68h buildings respectively. This is likely because the buildings shorter than h are sheltered from the flow above, due to shear layers shed from the taller buildings (Coceal et al. 2006). This would result in little time-mean flow penetrating to their heights, reducing the shear production of turbulence associated with the shorter building roofs. The turbulence would then become dominated by turbulence originating from the taller buildings. The taller building eddies are likely larger scale than the eddies generated by the roofs of the shorter buildings. This is because eddies generated at roof level are small in scale (as demonstrated in Sect. 4.7), and the shorter buildings are in the wakes of the taller buildings, and eddies in building wakes scale with building dimensions.
This explanation is consistent with Giometto et al. (2016) who conducted an LES of flow through a portion of Basel's urban canopy. They identified three modes in the building-height probability distribution function (PDF)-two above h and one below h (see their Fig. 2). The buildings with heights at the mode below h had little influence on the flow, and (upon inspection of their Fig. 4) a double-averaged velocity gradient maximum cannot obviously be identified at the mode below h. Yoshida and Takemi (2018) conducted LES of variable-height, flat-roofed buildings, with λ p = 0.25 between the ground and the height of the short buildings, and λ p = 0.0625 between the height of the short buildings and tall buildings. The ratio of the taller to the shorter building heights was 2.7 in one simulation and 7.25 in another. The number of l m local maxima and minima within the canopy corresponded to the number of building heights (i.e., 2). Therefore, unlike in the Xie et al. (2008) and Giometto et al. (2016) variablebuilding-height simulations, a l m minimum occurred below h at the height of the shortest buildings. One possible explanation is that the tall buildings in the Yoshida and Takemi (2018) simulations are sparse, so do not block flow from reaching the shorter buildings.
In contrast with the uniform-height, Xie et al. (2008) and Giometto et al. (2016) simulations, the Yoshida and Takemi (2018) simulations had two rather than one double-averaged velocity inflections. Again this is likely due to flow being able to penetrate to the heights of the shorter buildings in the Yoshida and Takemi (2018) simulations, and a clear separation of shear layers owing to the large difference in building heights.
As seen in Fig. 8, the transition to a linear l m I profile above z/h = 1 is much slower compared to the uniform-height buildings. The inertial sublayer was taken to be z/h = 3.0 − 4.0 since l m I was linear there, and a straight line was fitted in the region (see Fig. 8). Compared to the inertial sublayer fit there was generally a small decrease in l m I between z/h = 1-3. This is likely due to shear layers shed from the tops of the taller buildings.
Based on the results presented here it is proposed that the characteristic l m I shape in Fig. 7 can be applied generally to canopy geometries with uni-modal building height PDFs. A single l m I minimum can then be expected at the mode canopy height, where maximum shear is associated with the building roofs. The l m I minimum is smaller with increasing proportion of building roofs at the mode height. The transition to a linear l m I profile above the canopy is likely to occur for all canopy geometries, so long as the tallest buildings are significantly smaller than the boundary-layer depth, thus giving adequate scale separation of inner-and outer-layer turbulence for an inertial sublayer to develop. For urban canopies with multi-modal building height PDFs the schematic would require extension. Whether l m I minima exist at each mode depends on several factors such as whether there is a large proportion of buildings at the mode, if the mode occurs below h, whether the taller buildings are densely packed, and the height separation between modes.

Importance of Accounting for Solid Fraction Variation in Mixing-Length Turbulence Closures
The importance of accounting for (z) in the first-order momentum mixing-length closure and the consequent effect on predicted u was investigated. The variable-building-height LES of Xie et al. (2008) was used since of the datasets available, its (z) is most representative of that in a real urban neighbourhood. The detailed analysis is in the Appendix. It was found that if the solid fraction is not accounted for in the intrinsic mixing-length closure, then errors in the predicted double-averaged velocity of up to 17% occur. If the solid fraction is not accounted for in the comprehensive mixing-length closure, then errors in the predicted double-averaged velocity do not occur, so long as is not included in Eq. 8 when deriving l m,i C and solving the momentum equation.

Relationship Between Canopy-Top Shear Length Scale and Mixing Length
Based on similarity arguments Finnigan et al. (2015) argued that when most of the momentum is absorbed as drag by the canopy elements rather than by the ground, then there is only one relevant length scale in the canopy (S(z)C d (z)) −1 , which can be interpreted as a drag length scale L c (Belcher et al. 2003;Coceal and Belcher 2004). Here S(z) is the sectional obstacle area density (obstacle area facing the wind divided by canopy air volume) and C d (z) is the sectional drag coefficient. In dense vegetation canopies the approximation that S(z) and C d (z) are constant with height is also often made. This is consistent with the mixing-layer analogy, which suggests that there is one dominant turbulent length scale in the canopy, determined by the shear length scale L s (h) associated with the canopy-top mixing-layer.
With decreasing vegetation canopy density there is increasing flow penetration into the canopy, and there are length scales other than L s (h) that become important. For example those associated with eddies originating from shear near the surface (Watanabe and Kondo 1990), and penetration of inner layer eddies into the canopy and von Kármán street vortices (Poggi et al. 2004).
It is investigated here to what extent turbulence in urban canopies is dominated by mixinglayer eddies. If mixing-layer eddies are generated with length scale L s (h) due to the inflection in the double-average velocity profile at canopy top, and are dominant within urban canopies so that only they control turbulent mixing, it would be expected that L s (h) = α l m av across canopy geometries. α is a constant and l m av is the depth averaged value of l m within the canopy.
The shear length scale is defined here as This is the same as the definition given by Raupach et al. (1996) but with the addition of the second term in the brackets, which is included to account for the influence of solid fraction on the intrinsic spatial average of the velocity gradient. The component label i has been retained to allow the possibility of the shear being different in the x and y axes, while the shear of vertical velocity is zero since w I = 0. Figure 9 shows L s (h)/ l m av I plotted against λ p . For the 0 • aligned and staggered cube geometries it can be seen that with decreasing λ p the ratio L s (h)/ l m av I increases faster than linearly. α is not constant suggesting that mixing-layer eddies are not always the dominant turbulence within urban canopies. Similar L s / l m av I versus λ p relationships occurred when was 10% and 40% larger respectively than N-25-S, so increased with increasing λ p .
Inspection of individual velocity profiles reveals that inflections in the velocity profile occur at canopy top behind each roof, consistent with Letzel et al. (2008). The inflections are much weaker ∼ h downstream and extend little in the cross-stream. As suggested by Coceal et al. (2007a), the urban canopy top shear layer is highly heterogeneous and velocity profile inflections are local to individual buildings. As the canopy becomes sparser, fewer mixinglayer eddies are produced, and more flow penetrates into the canopy, making turbulence increasingly likely to be associated with other turbulence production mechanisms.
Ghisalberti (2009) argues that in 'obstructed shear flows', where part of a flow is obstructed via a permeable medium, the mixing-layer eddies generated at the interface between the permeable medium and unobstructed flow, cause environmental flows of this type to be dynamically similar. Ghisalberti (2009) observed that across canopy types (e.g., sediment beds, aquatic and terrestrial vegetation canopies, and coral reefs) the depth turbulence penetrates into the canopy is proportional to L c . However, urban canopies were found to be an outlier with smaller penetration depth. In moderately dense and dense vegetation canopies L s /h ≈ 0.1 − 0.5 (Raupach et al. 1996). For moderately dense and dense urban canopies (λ p ≥ 0.25) presented here, L s /h ranged between 0.06 and 0.16 (not shown). The mixinglayer eddies are generally smaller and more localized to the obstacles than in vegetation canopies, which explains why Ghisalberti (2009) found urban canopy turbulence to penetrate less deep into the canopy than in the other canopy types.
It can be seen from Fig. 9 that for the denser (λ p ≥ 0.25) 0 • staggered cube geometries that L s (h)/ l m av I varies less with λ p . This suggests that mixing-layer eddies likely make a large contribution to turbulent mixing in those geometries. As seen in Fig. 2a, b the velocity and TMF respectively are small within the canopy for the dense staggered geometries, since they have no unobstructed portions of the canopy. It is probable that most of the conversion of time-mean flow into turbulence occurs at canopy top, and is transported into the canopy.
The 45 • simulations have lower L s (h)/ l m av I values compared to the 0 • simulations with the same λ p (see Fig. 9). The 0 • and 45 • simulations with the same λ p have similar L s (h)/h (not shown), but the 45 • simulations have larger l m av I . Mixing-layer eddies generated by shear at the canopy top are therefore unlikely to be the cause of the differences in turbulence characteristics within the canopy between the 0 • and 45 • simulations.
In summary, the relation L s (h) = α l m av I does not hold across geometries, and turbulence is not dominated by mixing-layer eddies within the canopy, apart perhaps for the dense (λ p ≥ 0.25) 0 • staggered geometries.

Dispersive Momentum-Flux Characteristics
As a preliminary step towards developing DMF parametrizations, the behaviour of the three-dimensional DMF fields and horizontally-averaged profiles are investigated across geometries.
From Figs. 2c and 3c, it is seen that generally the horizontally-averaged DMF magnitude increases with height in the canopy until z/h = 0.7-0.95. The maximum magnitude of DMF tends to occur lower in the canopy for the sparser geometries. For the 0 • driving flows, DMF is typically negative in the middle and just below the top of the canopy, and positive near the surface (Fig. 2c). For the 45 • driving flows, DMF is negative throughout the canopy, with increasing magnitude from the surface until approximately z/h = 0.8 (Fig. 3c). These trends are understood in terms of the three-dimensional DMF field in the following sections.

0 • Driving Flows
Horizontal cross-sections of DMF through canopy geometry repeating units at z/h = 0.75 are shown in Fig. 10. Dispersive velocities are calculated using intrinsic averages to simplify interpretation as explained in Sect. 3.1. Inspecting the 0 • simulations in Fig. 10a, c, d, e reveals that DMF is generally negative along the leeward face of the buildings (ũ < 0,w > 0) and positive along the windward face (ũ < 0,w < 0), in agreement with Coceal et al. (2007c) and Yoshida and Takemi (2018). For λ p < 0.44 the negative fluxes on the leeward side have significantly larger streamwise extent than the positive fluxes on the windward side (Fig. 10a,  d), and therefore dominate the spatial average flux. This flow pattern exists in the middle and top portions of the canopy for these cases, and it explains why the horizontal average DMF is generally negative there in Fig. 2c.
For λ p = 0.44 (Fig. 10c, e) the streamwise extent of the negative fluxes on the leeward side is much reduced, since recirculations in wakes are very limited by the downstream building separation. In Fig. 10c it can be seen that for N-44-A the extent of negative fluxes from the leeward wall is still slightly larger than the positive fluxes at the windward wall. However, the magnitude of the positive fluxes is larger, and positive fluxes dominate the spatial average (as seen in Fig. 2c). From Fig. 10e it can be seen that for N-44-S the positive and negative fluxes occupy similar fractions of the canopy, with comparable magnitude, so that the horizontal average DMF is small (as seen in Fig. 2c). Figure 2d shows the ratio of DMF to TMF for the 0 • simulations. The λ p = 0.0625-0.35 geometries have ratio between 0.2-0.6 for z/h = 0.4-0.9. These values are consistent with other studies (Kanda et al. 2004;Xie et al. 2008;Castro 2017). DMF for λ p ≤ 0.25 tends to be convex within the canopy but TMF across all λ p tends to be concave. The ratio is generally largest at z/h ≈ 0.5. Near the surface the ratio is negative because DMF is positive there. The regions of negative flux in Fig. 10a  That the ratio of DMF to TMF is far from constant across the range of geometries presented here is further evidence that they are not governed by the same processes. DMF is related to the correlation of u and w (since u w = ũw ), and so requires knowledge of the threedimensional time-averaged flow field.

45 • Driving Flows
Compared to the 0 • simulation DMF cross-sections in Fig. 10, B-25-A-45 • in Fig. 10b has distinctly different characteristics, with negative fluxes dominating the cross-section. There are strong negative fluxes in the wakes of the buildings as in the 0 • driving flows, but at the windward walls (x/h = 2, y/h = 0 − 1 and x/h = 0 − 1, y/h = 2) there are no strong positive fluxes. There are two separate streams either side of each building where flow is diverted around them as also noted by Claus et al. (2012). At the windward walls they have downward motion (w < 0), like the flow at the windward walls in the 0 • driving flows. The flow streams have large positive horizontal velocity at the windward walls (ũ,ṽ > 0) so correlation withw < 0 gives negative fluxes, unlike the 0 • driving flows whereũ < 0 at the windward wall so that there are positive fluxes. This flow pattern extends throughout the canopy and explains why the B-25-A-45 • horizontally-averaged DMF is positive within the canopy in Fig. 3c.
The 45 • datasets have particularly large ratios in the lower canopy with values larger than and C-33-A-45 • y-component respectively. The ratio of 129.1 is off the axis limit and the particularly large value is due to very small TMF. It has already been shown by Castro et al. (2017) using the C-33-A-45 • dataset that DMF is a significant and at times dominant fraction of the total momentum flux. It is demonstrated for the first time here in the B-25-A-45 • case, suggesting that this might be a general result for 45 • driving flows. Using LES, Giometto et al. (2016) found for a portion of Basel's urban canopy (λ p = 0.5 near the surface), with two prevailing wind directions, that TMF and DMF are approximately equal up to h. This is consistent with our finding that DMF can be much larger when flow is oblique to building faces. They found the DMF to increase approximately linearly with height up to h, above which it starts to decrease, whereas the TMF peaks above h. The TMF profiles were slightly convex and reasonably similar in shape to the DMF profiles.
It may be expected that the magnitude of DMF increases with height (if λ p (z) is approximately constant) since |u| tends to increase with height. It is less obvious that DMF might increase approximately linearly with height for realistic urban canopy geometries. It is possible at high λ p in realistic geometries where buildings have varying orientations to the flow, that there is a mixture of convex and concave shaped DMF in different regions of the canopy. It was found here that when λ p ≥ 0.25, the 0 • and 45 • driving flows produce concave and convex profiles, respectively. Upon spatial averaging this could give rise to approximately linearly increasing spatially-averaged DMF with height. More CFD and experimental evidence of realistic urban geometries is required to see if this is the case.

Summary
The influence of different surface morphologies on urban canopy turbulence characteristics was investigated for numerous urban canopy LES and DNS datasets covering a range of idealized geometries (λ p = 0.0625 to 0.44, cubes and cuboids, uniform and variable building heights, driving flows at 0 • and 45 • to the building faces, and aligned and staggered arrays). A first-order momentum mixing-length turbulence closure was formulated and momentum mixing-length profiles examined across geometries, with a view towards developing a new urban-canopy turbulence model. The dispersive momentum-flux profiles were calculated and interpreted in terms of the three-dimensional flow field.
The intrinsic average of vertical gradients in velocity gives rise to a term that accounts for changes in λ p with height. Without it there are large discontinuities in the gradient whenever there are step changes in λ p , which cause spikes in the momentum mixing length. For a variable-height building geometry it was shown that excluding the extra term in the momentum mixing-length closure results in overestimation of predicted velocity by up to 17% near the top of the canopy. For denser urban canopies (λ p > 0.25) the error is expected to be larger. The intrinsic averaging approach is used commonly in urban canopy models (Martilli and Santiago 2007), and when K -theory turbulence closures do not account for the extra term large errors can be expected, since they involve horizontal averaging of velocity gradients.
It was found that the shear length scale at the canopy top does not vary with depth-averaged momentum mixing length within the canopy in a systematic way. This suggests that the mixing-layer analogy (Raupach et al. 1996) does not apply to turbulence in urban canopies, unlike in vegetation canopies, since it is not dominated by one length scale associated with mixing-layer eddies generated at canopy top. Shear-generated eddies at the top of urban canopies are local to the buildings (Coceal et al. 2007a), and do not penetrate as far down into the canopy as in vegetation canopies. Flow velocity can be considerable near the surface in unobstructed regions of urban canopies so that significant shear can be generated there. The momentum mixing length might be better interpreted as a combination of two length scales-the distance from the ground and from canopy top, both of which act to block eddies (Coceal et al. 2006).
The momentum mixing length demonstrated common characteristics across urban canopy geometries. For uniform-height, flat-roofed obstacles it increases from zero at the surface to a maximum roughly in the middle of the canopy, decreases until canopy top and increases rapidly above z/h = 1, before transitioning to a linear increase by approximately z/h = 1.5. The 0 • driving flows typically have a maximum momentum mixing length in the canopy of ≈ 0.3h.
There were some outliers to the trend, specifically the λ p = 0.44 staggered geometry and the 45 • driving flows as discussed in Sect. 4.1. Within the canopy these tended to have similar order magnitude wake production compared to shear production of turbulent momentum flux (away from the ground and canopy top), and larger magnitude ratios of dipersive to turbulent momentum flux than the other simulations. Wake production and dispersive momentum flux are both associated with building-scale spatial fluctuations in the flow, and are not directly related to the vertical gradient of the double-averaged velocity, as assumed by the mixinglength parametrization.
Using a variable-building-height dataset, it was demonstrated that the momentum mixinglength profile shape and magnitude below mean building height is similar to that with uniform-height geometries. Above mean building height, mixing length increases approximately linearly with height, but is slightly reduced from what would be expected in a logarithmic region, due to multiple shear layers associated with the tops of the buildings. A minimum in the mixing-length profile occurs at the mean building height but is larger compared to uniform-height geometries.
Von Kármán's parameter (constant) κ was calculated from the gradient of linear regressions fitted to momentum mixing-length profiles in the inertial sublayer. It was found that κ varies between 0.20 and 0.51, and there is no clear relationship between canopy geometry and κ. The simulations with larger domain heights tended to have larger κ. Whether the ratio of roughness height to boundary-layer depth is insufficient for adequate scale separation in the logarithmic region, as suggested by Jiménez (2004), requires further investigation.
Dispersive momentum-flux profiles for the 0 • driving flow datasets were in line with those in the literature (Kanda et al. 2004;Xie et al. 2008;Castro 2017), with dispersive to turbulent momentum-flux ratios of 0.2-0.6 for z/h = 0.4-0.9 in the λ p = 0.0625-0.35 simulations. The 45 • datasets have larger dispersive momentum flux, and ratios to turbulent momentum flux that exceed 1 for 0.05 < z/h < 0.5. As ratios of approximately unity have also been found for realistic geometries (Giometto et al. 2016), dispersive momentum flux in oblique flows may be more important than previously thought.

Implications for Parametrization
Many urban canopy models assume turbulent length scales that are constant up to mean canopy height and linearly increasing above (Hamdi and Masson 2008;Santiago and Martilli 2010;Nazarian et al. 2020). In NWP the main role of the urban canopy model is to provide a boundary condition for the flow above. The use of a constant turbulent length scale is not necessarily a poor approximation, if its value is determined so that the correct fluxes are provided to the levels above the canopy. However, if vertical profiles of velocity and scalars are of interest within the canopy, then using a turbulent length scale that varies with height in the correct way is important. TMF is proportional to l m 2 in first-order momentum mixinglength closures and proportional to the turbulent length scale in commonly used k − l TKE closures (Martilli et al. 2002;Schoetter et al. 2020), so the turbulent length scale (of which momentum mixing length is an example) is an important parameter.
Models that use a constant turbulent length scale often do not enforce the no-slip condition at the surface, so velocity gradients at the surface are very small. Canopy dispersion parametrizations based on the double-averaging approach involve velocity gradients, so it is necessary to correctly predict velocity for use in determining the vertical scalar distribution. This is especially the case near the surface where there are strong velocity gradients collocated with various scalar sources (e.g., pollution emissions, and anthropogenic and sensible heat).
Based on sensitivity tests (not shown), changing the turbulent length scale alone in urban canopy models does not always improve the predicted velocity profile, if the drag term is not treated properly. It is well known that the C d (z) tends to be larger near the surface than close to the tops of buildings (Macdonald 2000;Cheng and Castro 2002;Coceal et al. 2006;Santiago et al. 2008;Leonardi and Castro 2010), so it follows that the drag length scale L c is variable in urban canopies. Future multi-layer canopy models will likely need to account for both turbulent length scales and drag length scales that vary with height in the canopy.
It was found that the 45 • driving flow simulations have large turbulent momentum flux despite having small double-average velocity gradient. Turbulent transport was proposed as a possible explanation. When turbulent transport of turbulent momentum flux is important then first-order closures are insufficient. A possible pragmatic approach to a higher-order closure is to add an extra term to the first-order momentum mixing-length closure, so that u w = − l m 2 | ∂u/∂z | ∂u/∂z + N l , where N l represents turbulent transport (i.e., the non-local turbulence). Such an extension to the first-order momentum mixing-length approach has already been formulated by Wang (2014). This approach could also offer a method of parametrizing buoyancy effects which drive non-local momentum transport, in a way analogous to some convective-boundary-layer schemes (e.g., Lock et al. 2000), whilst maintaining first-order momentum mixing-length closure of mechanical turbulence.
It was argued in Sect. 4.8 that dispersive and turbulent transport are not governed by the same processes. Rather than representing dispersive momentum flux as part of the turbulence closure it could be treated separately. For example, if it is found more generally that dispersive momentum flux is an approximately linearly increasing function across different realistic geometries, then it could be taken as a function that is scaled with velocity and building morphology parameters such as λ p . Recirculations in building wakes cause the majority of dispersive momentum flux since they are where large magnitudeũ andw occur. Recirculations are likely predominantly driven by the flow at canopy top as in street canyons, so that u(h) is the relevant velocity scale. 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/. Fig. 11 (a) Profiles of l m I calculated according to Eq. 7: including the term (black), excluding the term (dashed blue), and excluding the term followed by removing the spikes and interpolating (dotted red). (b) l m C calculated according to Eq. 8: including (solid) and excluding (dashed) the 3 factor. (c) Solution to Eq. 15 using the full l m I from (a) (solid black), using the interpolated l m I from (a) and excluding the term (red dotted), alongside the solution to Eq. 16 excluding the 3 factor and using the dashed grey l m C from (b) (dotted grey). The solution to Eq. 16 including the 3 factor and using the solid grey l m C from (b) is plotted in (c) and is identical to the dotted grey curve. (d) The ratio of the dotted red and solid black profiles from (c)

Impact of Accounting for Solid Fraction Variation on Predicted Velocity
Here we investigate the impact on predicted velocity by accounting for solid fraction variation in the momentum mixing-length closure formulation. To do so, TMF and l m are taken to be known from the LES, and the corresponding double-averaged velocity is solved by discretising the momentum mixing-length closure. A simple discretization of the intrinsic averaged closure (Eq. 7) is given by u j+1 = u j + z ±(u w ) j+1/2 l m, j+1/2 where time and spatial averaging symbols, and spatial indices i have been dropped for brevity. j labels the LES grid levels and +1/2 indicates their midpoints, with linear interpolation of values at the adjacent levels where necessary. The ± is to be taken as negative (positive) when the value within the modulus in Eq. 7 is positive (negative). Discretization of the comprehensive averaged closure (Eq. 8) is given by where ± is negative (positive) when the value within the modulus in Eq. 8 is positive (negative). The solution to Eq. 15 using l m I calculated including the term is plotted as solid black in Fig. 11c. As a simple matter of consistency, this solution reproduces the true u I , since this was used along with u w I to derive l m I . Plotted as dotted red in Fig. 11c is the solution using l m I calculated when the (z) term is excluded and spikes are interpolated. The second term on the r.h.s. of Eq. 15 was dropped in the calculation for consistency with that approach. Inspecting the ratio of u I values when the (z) term is excluded and included (Fig. 11d), reveals that excluding the (z) term results in velocity overestimation of up to 17% just above h. Each time there is an (z) discontinuity the percentage error increases, and between discontinuities the percentage error decreases, since the magnitude of the error remains constant with increasing height, but the magnitude of the two velocities increases.
So long as the (z) 3 factor is dropped in Eq. 16 as well as in the l m C calculation, then the factor cancels, and the solution is the same as when the (z) 3 factor is included. Multilayer canopy model turbulence closures based on comprehensive averaging therefore do not encounter errors in predicted velocity due to incorrect application of the comprehensive spatial-averaging theorem to velocity gradients. The two solutions are plotted as one dotted grey curve in Fig. 11c.