Neighbourhood-Scale Flow Regimes and Pollution Transport in Cities

Abstract

Cities intimately intermingle people and air pollution. It is very difficult to monitor or model neighbourhood-scale pollutant transport explicitly. One computationally efficient way is to treat neighbourhoods as patches of porous media to which the flow adjusts. Here we use conceptual arguments and large-eddy simulation to formulate two flow regimes based on the size of patches of different frontal-area density within neighbourhoods. One of these flow regimes distributes pollutants in counter-intuitive ways, such as producing pollution ‘hot spots’ in patches of lower frontal-area density. The regimes provide the first quantitative definition of the ‘urban background’, which can be used for more precisely targeted pollution monitoring. They also provide a conceptual basis for further research into neighbourhood-scale air-pollution problems, such as parametrizations in mesoscale models, and the transport of fluid constituents in other porous media.

1 Introduction

Air pollution is the greatest environmental threat to human health (World Health Organization 2016), particularly in urban areas (Landrigan et al. 2018), where the majority of the world’s population now lives (United Nations 2018). However, it is difficult to assess the efficacy of air-pollution policy because in situ observations are too sparse to monitor transport processes within a city. Permanent networks usually comprise a handful of sites within an area of hundreds of square kilometres (World Health Organization 2016; DEFRA 2019; US EPA 2019) and typically avoid observations near large obstacles such as buildings and trees (Muller et al. 2013). Because low-cost sensors require careful calibration (Kelly et al. 2017) and have high post-processing and maintenance costs (Kumar et al. 2015; Pope et al. 2018), they cannot simply fill the gaps in permanent networks.

Numerical simulation is a useful supplement to in situ measurements, provided the model captures the relevant physics of the problem. Numerical models provide flow information at every point in the simulated domain—sufficient to investigate the transport equations term by term, which is not usually possible with observations. Advances in affordable computing power and accessibility over the last few decades have meant numerical simulations have proliferated (Kumar et al. 2011; Tominaga and Stathopoulos 2013; Blocken 2014; Toparlar et al. 2017). However, because of the computational demands of simulating the turbulent boundary layer, numerical flow models are forced to represent urban areas differently depending on their scale of interest (Martilli et al. 2002; Hood et al. 2014; Zhong et al. 2015). The behaviour of pollutants, and of other scalar fluid constituents, around the scale of a neighbourhood (1–2 km) remains particularly difficult to interpret and model (Belcher 2005; Nikolova et al. 2018). This gap in our understanding is unfortunate because neighbourhood-scale processes disperse pollutants from peaks beside busy roads to levels treated as the ‘urban background’, and may link urban-pollution models with weather forecasts (Xie 2011). Indeed, numerical weather prediction models are beginning to resolve neighbourhood scales for certain processes, including the UK Met Office’s London Model, which has been used to forecast fog at a horizontal resolution of 333 m (Boutle et al. 2016).

As it is computationally unfeasible to resolve turbulent flow around individual obstacles at the neighbourhood scale, for most applications, researchers need a simplified approach. One route is to resolve the main features of the urban form, in the sense that boundary conditions are imposed at the surfaces of the largest obstacles (usually only the buildings are resolved). A common approach when resolving buildings at the neighbourhood scale is to apply the Reynolds-averaged Navier–Stokes (RANS) equations, which require the turbulent flow to be parametrized. Such an approach has been widely used to investigate how the urban form affects flow and dispersion at the smaller end of neighbourhood scale, such as through small networks of streets (Wang and McNamara 2006; Letzel et al. 2008; Carpentieri and Robins 2015). Many of the earlier RANS models simulated flow and scalar transport around generic urban-like structures, such as networks of cuboids or idealized street canyons (Toparlar et al. 2017 and references therein). More recently, resolved-building RANS models have been used to investigate flow in real neighbourhoods (Toparlar et al. 2015; Antoniou et al. 2017; Juan et al. 2017; Gao et al. 2018). For more information on RANS models of flow in urban areas, see Kumar et al. (2011), Tominaga and Stathopoulos (2013), Blocken (2014), and Toparlar et al. (2017).

2 Porous Model of Neighbourhood-Scale Flow

2.1 Use of a Porous Model in Urban Areas

Another way of approximating neighbourhood-scale flow is to treat the urban canopy layer as a porous medium. Instead of resolving the buildings and other roughness elements in the urban area, all aerial parts of the urban area are represented indirectly through a distributed momentum sink. This is achieved by averaging the equations of motion over a volume that contains multiple roughness elements, which amalgamates the drag force from each roughness element into a continuous drag force throughout the urban canopy layer (UCL). This averaging process introduces additional terms into the momentum equation which do not appear in the equations for the freestream flow (Sect. 3 below, Coceal and Belcher 2004; Lien et al. 2005).

The porous-canopy method has been used extensively to study flow and scalar exchange in vegetation, particularly in homogeneous stands, from which a theoretical framework has emerged to explain certain statistically consistent features (Raupach and Thom 1981; Raupach et al. 1996; Finnigan 2000). A characteristic feature of flow through vegetation is that momentum is transferred over the depth of the canopy. The momentum transfer creates a strong inflection in the mean streamwise wind-speed profile, which is approximately exponential within the canopy and logarithmic above, and a layer of high shear around the top of the canopy. This shear generates Kelvin–Helmholtz type instabilities that dominate turbulence and ‘coherent’ motions, analogous to the dominant processes in a plane mixing layer (Raupach et al. 1996). Using the analogy of vorticity thickness in a mixing layer, Raupach et al. (1996) reduced canopy turbulence to a single length scale, $L_{\text{s}} = U_{\text{H}} /U'$, where $U_{\text{H}}$ is the mean streamwise velocity component U at the height of the canopy H, and $U^{\prime} = {\text{d}}U / {\text{d}}z$ at z = H. The shear length scale $L_{\text{s}}$ provides a rough estimate of the diameter of the turbulent eddies, equating to around 0.5H for medium-density vegetation. Eddies of diameter $L_{\text{s}}$ dominate the exchange of momentum and scalar quantities to and from the canopy. The inverse of the length scale $L_{\text{s}}$ also represents the wavenumber of the fastest-growing instability at the top of the canopy. The coherent eddies in the shear layer induce characteristic patterns in the spatially averaged vertical profiles of higher-order moments, for example, of the variances or the skewness (Finnigan 2000).

Several studies have used the porous-canopy approach in mesoscale models to parametrize the flow through and above urban areas (Martilli et al. 2002; Santiago and Martilli 2010). Coceal and Belcher (2004, 2005) extended the work of Belcher et al. (2003) to calculate how the boundary layer adjusts upon entering an urban area (Sect. 2.3 below). Their model simulated mean velocity profiles that compared well, before and after the flow had adjusted, with wind-tunnel and field measurements reported by Davidson et al. (1995, 1996) and Macdonald (2000). Lien et al. (2005) derived a modified two-equation ($k$ –ϵ) turbulence model, where $k$ is the turbulence kinetic energy and ϵ is the turbulence dissipation rate, to simulate the flow through and over an array of buildings, in which clusters of buildings in the array are averaged and treated as a porous medium. Compared to high-resolution RANS simulations of flow over resolved cubes, this porous-canopy model generated very similar vertical profiles of the mean wind speed, but less satisfactory profiles for turbulent variances, particularly in the regions around the tops of the buildings and very close to the ground (Lien and Yee 2005). Di Sabatino et al. (2008) adapted the porous model in Macdonald (2000) to simulate spatially averaged wind-speed profiles over London, Toulouse, Berlin, and Salt Lake City by varying the vertical profiles of the frontal-area density, $\lambda_{\text{f}}$ (i.e., the total frontal area per unit ground area). Hang and Li (2010) derived a modified $k$ –ϵ porous model to simulate the flow over arrays of aligned cubes. Compared to wind-tunnel observations, their porous model simulated macroscopic properties of the flow well, but under-predicted turbulence at the leading edges of the arrays. More recently, porous models have been used to investigate wider features of urban microclimates, such as heat-island effects (Hu et al. 2012; Wang and Li 2016) and the influence of urban trees on flow and pollutant dispersion (Krayenhoff et al. 2015; Salim et al. 2015; Wang et al. 2018).

2.2 Applicability of the Porous Model to the Urban Boundary Layer

The porous model allows urban areas to be parametrized with a relatively small number of morphological variables. Detailed information around individual obstacles is lost, but cities can be represented simply enough to make numerical models of the neighbourhood scale, and upwards, computationally feasible. However, the model and its supporting theoretical assumptions must be used with caution when approximating urban areas. Direct numerical simulation (DNS) of idealized building arrays has shown that the eddy structure in the urban canopy shear layer is more spatially complex than its counterpart in homogeneous vegetation (Coceal et al. 2007), for example, quasi-coherent eddies may form in the wakes of the bluff elements (Böhm et al. 2013). The spatially averaged mean velocity profile in urban-like canopies—i.e., between squat, sharp-edged obstacles—is also rarely exponential (Castro 2017), having instead a sharp shear layer between the top of the canopy and the flow aloft.

However, the studies revealing the more striking differences between urban and vegetative flows have tended to use exaggerated urban forms, such as the use of arrays of cuboids (Coceal et al. 2007; Leonardi and Castro 2010; Castro 2017 and references therein). Further, because of the constraints in computational resolution and wind-tunnel capacity, even models of real urban areas typically only resolve the simplified forms of buildings and the streets (Xie and Castro 2009; Carpentieri et al. 2012; Toparlar et al. 2015; Antoniou et al. 2017; Juan et al. 2017). In reality, however, most urban areas are not just networks of buildings with small-scale irregular forms, but also have trees, street furniture, and other obstacles, all of which absorb momentum. Models of urban flow omitting those features also omit their effects on the flow and the transport of scalar quantities. For example, trees strongly influence urban flow and dispersion, particularly if they are of comparable height to buildings (Wang et al. 2018). Trees act as a direct momentum sink for the mean flow and reduce downwards turbulent transport of high velocity air from above the canopy (Giometto et al. 2017) as well as upwards turbulent transport of high-concentration air pollutants emitted along streets (Jeanjean et al. 2015).

Many of the scaled vertical profiles of flow statistics are similar across vegetation and urban areas (Finnigan 2000; Christen 2005, Fig. 5.1). Several high-resolution studies of urban flow have also identified Kelvin–Helmholtz instabilities in the shear layer around the tops of the buildings (Letzel et al. 2008; Salizzoni et al. 2011; Li et al. 2015), which mix air intermittently with the air aloft (Louka et al. 2000; Cui et al. 2004; Cai 2012; Li et al. 2015) in a similar manner to the processes in vegetation. These features, together with the comparisons with field and wind-tunnel observations in the studies cited above, support the use of the porous model for investigations of flow in urban areas as well as in vegetation. However, the results need to be interpreted carefully. While the porous model simulates spatially averaged statistics well, it offers no information about the flow or dispersion around the individual obstacles, which may differ significantly from spatial averages (Coceal et al. 2007).

2.3 Spatial Inhomogeneity and Adjustment

For flow into a city, such as from a rural area, there is an adjustment over a streamwise distance $x_{A}$ in which the mean flow balances the aerodynamic drag of the urban roughness (Coceal and Belcher 2004). Once the flow has adjusted to the presence of the canopy, the net mean vertical velocity component almost totally disappears but the large vertical wind shear generates Kelvin–Helmholtz like instabilities around the top of the canopy. It is worth noting that the ‘coherent structures’ generated by these instabilities refer to motion transporting a range of smaller vortices, rather than necessarily meaning the well-defined rolls of fluid that appear in conceptual diagrams (Bailey and Stoll 2016). The spatial complexity of urban areas may prevent coherent structures of a single length scale $L_{s}$ from dominating turbulent exchange (Coceal et al. 2007). As a rough guide, Okaze et al. (2015) calculated a turbulent length scale $\ell \approx 0.1$ H using high-resolution large-eddy simulation (LES) of the flow over an array of cubes.

Most real urban areas are heterogeneous (i.e., patchy) in structure and density. For the flow within a city with varying roughness, flow adjustment occurs over a distance $x_{A}$ determined by the local characteristics of the urban roughness (Coceal and Belcher 2004, 2005). We have some idea of how urban morphology affects pollutant transport at the neighbourhood scale. For example, the mean age of air (Buccolieri et al. 2010) and average pollutant concentration (Yuan et al. 2014) increase with the plan-area density $\lambda_{p}$ (i.e., the total plan area per unit ground area). Urban areas are also better ventilated with increasing variation in the height of obstacles (Hang et al. 2012).

However, we do not have a clear conceptual picture of how changes in neighbourhood-scale urban form affect the transport of pollutants and other scalar quantities. Here we present a simple model of inhomogeneous flow within an urban area, in which air moves between areas of high and low frontal-area density. We use the LES approach to investigate (a) the dynamical adjustment of the flow as it moves over a neighbourhood of heterogeneous roughness; and (b) the effect of the adjustment on pollutant transport, particularly features that may lead to improved parametrizations in mesoscale models.

3 Method

3.1 Transport Equations

We use right-handed Cartesian tensor notation, with the Einstein summation convention, and indices $\left( {i,j,k} \right)$ take values (1, 2, 3) respectively. For example, $u_{i}$ is the velocity in the $x_{i}$ direction, with i = 1, 2, 3 representing the streamwise ($x$), spanwise ($y$) and vertical ($z$) directions. We denote $\varvec{x} = \left( {x,y,z} \right)$, $(u_{1} ,u_{2} ,u_{3} ) = \left( {u,v,w} \right)$, and time as $t$.

We adapt LES models previously used for vegetation canopies by modelling the urban area as a porous body, with the obstacles represented as a momentum sink. The representation of the urban form is informed by porous models of urban areas, particularly studies by Coceal and Belcher (2004, 2005), which used mixing-length closure, and the RANS study of Lien et al. (2005). We used the ‘superficial’ or ‘extrinsic’ averaging procedure; i.e., the spatial averaging operation was performed over the total volume $V_{t}$ of the UCL including both solid obstacles $V_{s}$ and fluid parts $V_{f}$, where $V_{t} = V_{f} + V_{s}$. For a flow quantity $\phi$ such as velocity or stress,

$$ {\phi = \frac{1}{{V_{t} }}\mathop \int \limits_{{V_{t} }} \phi dV.} $$

(1)

Other LES studies of urban areas have used different procedures—for example, the ‘intrinsic’ average where the quantities are averaged only over the fluid volume $V_{f}$ (e.g., Giometto et al. 2016). The superficial average treats pressure and velocity gradients as continuous in space, making it more straightforward to focus on the general dynamical effects of density changes, rather than discontinuities at the top of urban canopy layer. The superficial average also generates the simplest form of the averaged transport equations (Kono et al. 2010; Xie and Fuka 2018), making it easier to implement in a three-dimensional LES model than the intrinsic average. There are other situations where the intrinsic average operation may be more appropriate than the superficial, for example, when comparing simulated results to observations, because the intrinsic average produces results representative of local conditions within the fluid (Schmid et al. 2019). See Lien et al. (2005), Böhm et al. (2013), Xie and Fuka (2018), and Schmid et al. (2019) for further discussion of different averaging operations for urban areas.

Here we envisage an urban area comprising roughness elements of various shapes and sizes (buildings, trees and plants, signage, infrastructure, parked vehicles, etc.). The averaging is over an area that: (a) covers multiple roughness elements, but (b) is small compared with the distance over which the mean characteristics of roughness (e.g., frontal-area density) varies. The vertical volume averaging is very thin in order to properly resolve the flow gradients.

We used the LES mode of version 3.6.1 of the Weather Research and Forecasting model (WRF) (Skamarock et al. 2008) to solve the transport equations. The WRF model solves discretized forms of the spatially averaged momentum equations¹ using the Runge–Kutta time-integration scheme (Wicker and Skamarock 2002),

$$ \begin{array}{*{20}c} {\frac{{\partial \left\langle {u_{i} } \right\rangle }}{{\partial x_{i} }} = 0, } \\ \end{array} $$

(2a)

$${\frac{{\partial \left\langle {u_{i}} \right\rangle }}{\partial t} + \frac{{\partial \left\langle {u_{i}} \right\rangle \left\langle {u_{j}} \right\rangle }}{{\partial x_{j}}} = - \frac{1}{\rho}\frac{\partial \left\langle {p} \right\rangle}{{\partial x_{i}}} + \frac{{\partial \tau_{ij}}}{{\partial x_{j}}} + B_{i} + f_{c} \epsilon_{ij3} \left({{\left\langle {u_{j}} \right\rangle } - U_{g,j}} \right) + f_{i},} $$

(2b)

where the kinematic mean stress tensor, $\tau_{ij}$, represents the subgrid-scale (SGS) stresses; $\left\langle p \right\rangle$ is the spatially-averaged pressure; $\rho$ is the air density; $B_{i}$ is the buoyancy force: $B_{i} = - \delta_{i3} g\theta '/\bar{\theta }$, where $\delta_{ij}$ is the Kronecker delta, $\bar{\theta }$ is the potential temperature for hydrostatic balance, and $\theta '$ is the temperature variations with respect to $\bar{\theta }$; $f_{c}$ is the Coriolis parameter; ϵ_ij3 is the alternating unit tensor; and $U_{g,j}$ is the geostrophic velocity. The term $f_{i}$ is related to the drag force (see Sect. 3.2 below). The angled brackets $\left\langle \cdot \right\rangle$ denote the spatially averaged quantities, for example, $\left\langle {u_{i} } \right\rangle$ is the i component of the averaged velocity field. Equation 2b is closed by parametrizing the stress $\tau_{ij}$ as

$$ \begin{array}{*{20}c} {\tau_{ij} = - 2\nu_{SGS} S_{ij} ,} \\ \end{array} $$

(3)

$$ \begin{array}{*{20}c} {S_{ij} = \frac{1}{2}\left( {\frac{{\partial \left\langle {u_{i} } \right\rangle }}{{\partial x_{j} }} + \frac{{\partial \left\langle {u_{j} } \right\rangle }}{{\partial x_{i} }}} \right),} \\ \end{array} $$

(4)

$$ \begin{array}{*{20}c} {v_{SGS} = c_{k} \sqrt {\left\langle {e_{SGS} } \right\rangle } \left( {\Delta x\Delta y\Delta z} \right)^{{\frac{1}{3}}} ,} \\ \end{array} $$

(5)

where $\left\langle {e_{SGS} } \right\rangle$ is the SGS turbulence kinetic energy (TKE) and $c_{k} =$ 0.10 is a modelling constant. The prognostic equation for the evolution of the term $\left\langle {e_{SGS} } \right\rangle$ is,

$${\frac{{\partial \left\langle {e_{SGS}} \right\rangle}}{\partial t} + \frac{{\partial \left\langle {u_{j}} \right\rangle \left\langle {e_{SGS}} \right\rangle}}{{\partial x_{j}}} = \nu_{SGS}\frac{\partial}{{\partial x_{j}}}\left({\frac{{\partial \left\langle {e_{SGS}} \right\rangle}}{{\partial x_{j}}}} \right) + P + F - \frac{{C_{\epsilon} \left\langle {e_{SGS}} \right\rangle}}{{\left({\Delta x\Delta y\Delta z} \right)^{{\frac{1}{3}}}}}, } $$

(6)

where P represents the shear- and buoyancy-production terms (Skamarock et al. 2008), F is a cascade term to account for the SGS TKE lost due to the interaction of SGS eddies with the obstacles in the urban area (see Eq. 11 below), and C_ϵ is the dissipation coefficient (Moeng et al. 2007).

After model spin-up (see Sect. 3.4 for details), we introduce continuous sources of passive scalars to represent traffic fumes at the horizontal centres of each grid cell at z =0.15H. Pollutant fluxes from traffic fumes—such as ultrafine particles—vary spatially and temporally (Levy and Hanna 2011). However, at the neighbourhood scale, scalars at a particular location can be treated as the sum of releases of many sources (Belcher 2005). Here, we focus on ‘road-to-ambient’ scalar processes, up to the neighbourhood scale (Harrison et al. 2018). We specify a unity-emission source term, Q (in μg m⁻³ s⁻¹), in the filtered advection–diffusion equation that the LES model solves to model the transport of a passive scalar,

$$ \begin{array}{*{20}c} {\frac{\partial \left\langle c \right\rangle }{\partial t} + u_{i} \frac{\partial \left\langle c \right\rangle }{{\partial x_{i} }} = \frac{\partial }{{\partial x_{i} }}\left( {K_{c} \frac{\partial \left\langle c \right\rangle }{{\partial x_{i} }}} \right) + Q,} \\ \end{array} $$

(7)

where $\left\langle c \right\rangle$ is the filtered concentration for a passive scalar, and K_c is the SGS eddy diffusivity. We applied the ‘superficial’ averaging operation to all of the variables, meaning the pollutant is diluted in the same volume independently of changes to the frontal-area density (see Sect. 3.3 below).

3.2 Approximating the Drag Force

The parametrization of the drag force $f_{i}$ in Eq. 2b is performed by spatially averaging the localized drag from the individual elements of the urban area, assuming the drag force is proportional to the square of the flow speed. The viscous component of the drag is assumed negligible compared with the much larger inertial component (Hamlyn and Britter 2005). For the form-drag component, consider an array of roughness elements distributed over a total area $A_{t}$, where each element has mean height $H$, frontal area $A_{f}$ and drag coefficient $C_{d} \left( z \right)$. The drag force over a thin layer $dz$ of each element at height $z$ is

$$ \rho{f_i}{A_t}{dz} = {\rho U^{2} \left( z \right)C_{d} \left( z \right)A_{f} \frac{dz}{H}.} $$

(8)

The thin averaging volume at height $z$ is $A_{t} dz$. The total force per unit volume (including both solid and fluid elements) is, therefore,

$$ \begin{array}{*{20}c} {\rho f_{i} = \rho \frac{{C_{d} \left( z \right)\Sigma A_{f} }}{{HA_{t} }}\left( {\left\langle {u_{j} } \right\rangle \left\langle {u_{j} } \right\rangle } \right)^{{\frac{1}{2}}} u_{i} ,} \\ \end{array} $$

(9)

which amalgamates the spatially discontinuous drag force from each roughness element into a continuous resistive body force throughout the urban area (Belcher et al. 2003). In Eq. 9, $C_{d} \left( z \right)$ is the sectional drag coefficient, which is largely unknown for real urban areas, for which bulk transfer coefficients such as $u_{*} /U$ are easier to measure (Nordbo et al. 2013; Peng and Sun 2014). In arrays of cubes, the drag coefficient $C_{d} \left( z \right)$ decreases with height because the wind speed increases from the surface (Cheng and Castro 2002; Kono et al. 2010; Castro 2017), although the value of $C_{d} \left( z \right)$ can vary by almost an order of magnitude depending on whether the cubes are aligned or staggered (Coceal et al. 2006). As a simplification, we adopted the common approach of taking a height-averaged sectional drag coefficient $\bar{C}_{d}$ (Macdonald 2000; Martilli et al. 2002; Coceal and Belcher 2004; Yang et al. 2006). We assumed $\bar{C}_{d}$ = 1.3, which is derived from wind-tunnel observations (Coceal and Belcher (2004²), who analyze the results of Cheng and Castro (2002)), and is within the range of values measured in urban-like canopies (e.g., Castro 2017). Taking the frontal-area density $\lambda_{f} = \sum {\text{A}}_{\text{f}} /{\text{A}}_{\text{t}}$, $\left| U \right| = \left( {\left\langle {u_{j} } \right\rangle \left\langle {u_{j} } \right\rangle } \right)^{{\frac{1}{2}}}$, and cancelling the density $\rho$, Eq. 9 simplifies to

$$ \begin{array}{*{20}c} {f_{i} = - \frac{{\bar{C}_{d} \lambda_{f} }}{H} \left| U \right|u_{i} .} \\ \end{array} $$

(10)

We also add a cascade term,

$$ \begin{array}{*{20}c} {F = - 2\frac{{\bar{C}_{d} \lambda_{f} }}{H}\left| U \right|\left\langle {e_{SGS} } \right\rangle } \\ \end{array} $$

(11)

into the transport equation for the SGS TKE, $\left\langle {e_{SGS} } \right\rangle$, in Eq. 6 to account for additional dissipation of kinetic energy from air–obstacle interactions at scales below the spatial filter, bypassing the usual inertial cascade (see discussion in Lien et al. 2005 in the context of RANS modelling of urban areas; and Shaw and Schumann 1992, Shaw and Patton 2003, and Dupont and Brunet 2008 in the context of LES investigations of vegetation canopies). Without this SGS sink, the model overestimates the kinematic turbulent momentum flux within the UCL, as well as the TKE above the UCL (Lien and Yee 2005).

3.3 Simulated Cases

We used this model to simulate five different cases comprising:

• 1. a homogeneous control case—uniform frontal-area density across the entire neighbourhood; and

• 2. five cases of alternating values of $\varvec{\lambda}_{\varvec{f}}$ (4H, 6H, 8H, 12H, and 16H cases)—alternating patches of high and low frontal-area density spanning the entire domain east–west and, respectively, extending 4H (80 m), 6H (120 m), 8H (160 m), 12H (240 m), and 16H (320 m) south–north.

We simulated these patches of high and low frontal-area density by alternating the value of $\lambda_{f}$ in Eqs. 10 and 11. For each simulated case, we calculate a mean frontal-area density,$\lambda_{m} = \sum {\text{A}}_{\text{f}} /{\text{A}}_{\text{t}} = 0.25$ across the entire domain, which we apply uniformly in the homogeneous case. For the alternating cases, we simulated alternating patches of high and low frontal-area density, which we respectively refer to as ‘dense’ ($\lambda_{f} = \lambda_{d} = 1/3$) and ‘sparse’ ($\lambda_{f} = \lambda_{s} = 1/6$) (Fig. 1). Each of the alternating-density cases contained equal numbers of sparse and dense patches. The horizontal averages of the ratio $\bar{C}_{d} \lambda_{f} /H$ at a given height z in Eqs. 10 and 11 are, therefore, equal across all the simulated cases. Here we represent the effective heterogeneous frontal-area density $\lambda_{f}$ in neighbourhoods where the plan-area density of the buildings $\lambda_{p}$ remains approximately constant. This is not possible in the case of flow perpendicular to aligned cubes, the subject of many previous studies of urban-like areas, for which $\lambda_{f} = \lambda_{p}$. However, by considering the effect on the flow of all of the obstacles in an urban area, not just the buildings, the value of $\lambda_{f}$ varies over small areas even within a neighbourhood (e.g., neighbouring streets having very different amounts of vegetation, or signage and infrastructure). In these cases, the presence of non-building obstacles may contribute significantly to the total drag exerted on the flow while occupying relatively little plane-surface area and volume in comparison with the buildings (i.e., the quantities $\lambda_{p}$ and $V_{s}$ are approximately constant between patches). Alternating variations of the frontal-area density $\lambda_{f}$ may also arise due to variations in the size, form, or arrangement of obstacles with respect to the wind direction. Formulating the model in terms of alternating variations of the frontal-area density $\lambda_{f}$ raises the question of whether the drag coefficient $\bar{C}_{d}$ also varies between patches. Wind-tunnel (Hagishima et al. 2009; Zaki et al. 2011) and LES (Nakayama et al. 2011) studies across a variety of building geometries indicate the value of $\bar{C}_{d}$ does appear to vary slightly with the frontal-area density $\lambda_{f}$, with a peak in the range $\lambda_{f} \approx$ 0.15–0.2. However, there is considerable scatter in the results depending on the geometry of the modelled urban areas and the value of $\bar{C}_{d}$ varies relatively little for the range of $\lambda_{f}$ values considered here. We, therefore, assume $\bar{C}_{d}$ is constant across the simulated domain, while acknowledging this slightly overestimates the difference in drag between the sparse and dense patches.

Fig. 1

Schematic of the simulated domain, showing alternating density patches within an urban neighbourhood, with $\lambda_{d}$ and $\lambda_{s}$ referring to the frontal-area densities of the dense and sparse patches, respectively. The 12H case is shown from the post-processing frame-of-reference with the x-direction aligned south–north, approximately with the mean wind direction

3.4 Numerical Details and Post-Processing

The simulated domain is 191 $\times$ 191 $\times$ 79 grid cells in the streamwise, spanwise, and vertical directions, respectively (Fig. 1). The horizontal grid resolution is 10 m in each direction, and increased vertically from around 0.67 m within the lower half of the UCL to around 60 m at the top of the domain. The mean height of the roughness elements (and therefore the UCL), H, is set to 20 m. We simulated the flow under neutral conditions. We included a dampening layer of 300 m at the top of the domain to minimize wave reflection because the LES model does not always maintain a deep neutral boundary layer over long simulations (Nottrott et al. 2014). The potential temperature $\theta$ was held constant at 283.15 K within the neutral layer at the bottom of the domain (up to z = 475 m), increasing to $\theta =$ 291.7 K at the top of the simulated domain. The geostrophic velocity components are set to $U_{g}$ = 10 m s⁻¹ and $V_{g}$ = 5 m s⁻¹, giving a mean wind speed of 1.7 m s⁻¹ from a flow direction of 158^o at H = 20 m—i.e., at the top of the canopy—and at an angle of 22° to the patches (see Fig. 1). The difference in mean wind direction with respect to the geostrophic wind direction is due to the Coriolis force, with the steepest change occurring across the inversion layer at the top of the boundary layer. The spin-up time was 9 h, with cyclic boundary conditions for all dynamical variables (u, v, w, $e_{SGS}$, T, etc.) in both the x- and y-directions. Boundary conditions for the scalars were non-cyclic; at the inlet boundaries, all scalars were set to zero and the outlet boundaries were open. In the model set-up, we used the meteorological convention where the x-direction is aligned west–east and the y-direction south–north. For post-processing, we rotate our frame-of-reference so that the $x$-direction is aligned south–north, approximately with the mean wind direction, as shown in Fig. 1.

After the spin-up, we introduced the scalars and ran the simulations for a further 120 min, taking samples at intervals of 3 s. We time averaged over the latter 100 min (denoted by T) of these samples (i.e., $t_{0} = 20$ min to $t_{0} + T = 120$ min). Sensitivity testing (not shown) indicated scalar concentrations reached local pseudo-equilibrium 15–20 min after being introduced, where pseudo-equilibrium was defined as the time at which mean concentrations at z = H deviated within $\pm 5\%$ from those at t = 120 min. The concentrations before the time $t_{0}$ were therefore excluded from the analysis. This process generated a three-dimensional time series of 2000 resolved samples in the form $\phi \left( {x,y,z,t} \right)$. We derived the resolved turbulent statistics using (a) time averages over the sampling period; and (b) spatial averages along the y-direction ($L_{y}$ = 1910 m), over which the turbulent statistics are homogeneous. For each resolved variable, $\phi$ this generated a two-dimensional function,

$$ \begin{array}{*{20}c} {\left[ {\bar{\phi }} \right]\left( {x,z} \right) = \frac{1}{{TL_{y} }}\mathop \int \limits_{0}^{{L_{y} }} \mathop \int \limits_{{t_{0} }}^{{t_{0} + T}} \phi \left( {x,y,z,t} \right) dtdy, } \\ \end{array} $$

(12)

where the resolved fluctuating component of $\phi$ around $\left[ {\bar{\phi }} \right]$ is defined as $\phi '\left( {x,y,z,t} \right) = \phi \left( {x,y,z,t} \right) - \left[ {\bar{\phi }} \right]\left( {x,z} \right)$. The covariance of two resolved variables $\phi$ and $\varphi$ is in turn defined as

$$ \begin{array}{*{20}c} {\left[ {\overline{{\phi^{\prime}\varphi^{\prime}}} } \right] = \frac{1}{{TL_{y} }}\mathop \int \limits_{0}^{{L_{y} }} \mathop \int \limits_{{t_{0} }}^{{t_{0} + T}} \phi '\left( {x,y,z,t} \right)\varphi '\left( {x,y,z,t} \right)dtdy,} \\ \end{array} $$

(13)

with the autovariance $\left[ {\overline{{\phi^{\prime}\phi^{\prime}}} } \right]$ similarly defined. The time-averaged resolved TKE is defined as $\left[ {\bar{e}} \right] = \frac{1}{2}\overline{{u_{i} 'u_{i} '}}$. We use square brackets [$\cdot$] to denote the spatial averaging along $y \in \left[ {0,L_{y} } \right]$, and we omit the angled brackets denoting the volume averaging below.

4 Results—Dynamics

4.1 Comparison with Published Data

We examined a selection of flow statistics: the mean streamwise velocity component $U = \left[ {\bar{u}} \right]$, the square root of the magnitude of the shear stress $F_{u} = \left( { - \left[ {\overline{{w^{\prime}u^{\prime}}} } \right]} \right)^{{\frac{1}{2}}}$, and the streamwise turbulence intensity $\sigma_{u} = \left[ {\overline{{u^{\prime}u^{\prime}}} } \right]^{{\frac{1}{2}}}$. Figure 2 presents (a) the vertical profiles of U for the homogeneous control case and the sparse patches of the 16H case; (b) $F_{u}$ for the homogeneous case; and (c) $\sigma_{u}$ for the homogeneous case. The vertical axis is normalized by z = H = 20 m, and the velocity moments with either (a) $u_{H}$ = |U| at $z/H$ = 1; or (b) friction velocity $u_{ *} = \left( {\left[ {\overline{{w^{\prime}u^{\prime}}} } \right]^{2} + \left[ {\overline{{w^{\prime}v^{\prime}}} } \right]^{2} } \right)^{{\frac{1}{4}}}$ at $z/H$ = 1. Also plotted are data from Castro et al. 2006 (hereafter CCR06), a wind-tunnel study; Coceal et al. 2007 (hereafter CDT07) and Leonardi and Castro 2010 (hereafter LC10), both DNS studies; and Yang et al. 2016 (hereafter YSMM16), an LES study. Each of these studies models the flow over arrays of cubes. In both CCR06 and CDT07, measurements were taken on top of a cube (which the authors denote position 0), behind (position 1), and in front of a cube (position 2), and in the gap between two adjacent cubes (position 3)—see Fig. 1 of CCR06 and Fig. 1b of CDT07. The measurement position 2 and the mean quantities (average of positions 0–3) provide the most meaningful comparison to our simulations because the others are strongly influenced by small-scale features that the porous model cannot capture. The spatially averaged vertical profiles in Fig. 2a, using data from LC10 and YSMM16, were generated using the superficial averaging operation.

Fig. 2

Canopy spatially averaged mean vertical profiles of the a mean streamwise velocity component $U/u_{H}$, b shear stress $F_{u} /u_{*}$, and c streamwise turbulence intensity $\sigma_{u} /u_{*}$. Values of the frontal-area density $\lambda_{f}$ are given in the legends. The sparse profile in a is the spatial average of the sparse patches of the 16H case. The wind-tunnel (WT), DNS and LES cases in the legends denote the type of model used in each study. The m and p2 symbols in the legends for CCR06 and CDT07 respectively denote the mean (average of positions 0–3) and position 2 quantities

Our simulated U component above the canopy, where the profile is logarithmic, agrees well with both the observations at position 2 and the mean profile from CCR06.³ Within the canopy, the vertical profile of U for our homogeneous case is approximately exponential, and closely matches that found by CCR06 for position 2 in the top two-thirds of the canopy. However, it deviates from the CCR06 mean and from LC10 for a frontal-area density $\lambda_{f} =$ 0.25, which both reflect the non-exponential behaviour observed in flow over cubic arrays (Castro 2017). Our sparse profile ($\lambda_{f} = 1/6$) shows qualitative agreement with LC10 for a frontal-area density $\lambda_{f} =$ 0.16 in that the ratio $U/u_{H}$ increases at each height z compared to the case where $\lambda_{f} =$ 0.25, but again our model does not simulate the non-exponential behaviour seen in the cubic arrays. We attribute these differences to the different model configurations; i.e., we are simulating flow through a porous medium rather than over cubes, and therefore expect a less pronounced shear region at the top of the canopy. We also suspect the filter resolution accounts for some of the difference. Our LES mesh is quite fine in the vertical (H/20 within the canopy), but coarse in the horizontal (H/2) compared with the DNS mesh of H/32 in LC10, for example. We note that our homogeneous profile is very similar to that in YSMM16, who simulate the flow over cubes but using a relatively coarse LES mesh (H/8).

Our model simulates greater $F_{u}$ and $\sigma_{u}$ values above the canopy than were observed in CCR06 and CDT07 (Fig. 2b, c), which we attribute to the very low boundary-layer heights ($H_{i}$) in CCR06 ($H_{i}$= 7.4H) and CDT07 ($H_{i}$= 8.0H), whilst ours is much higher at $H_{i} \approx$ 25.0H. If the parameters $F_{u}$ and $\sigma_{u}$ are scaled vertically by ($H_{i} - H$) rather than H (not shown), the shapes of the vertical profiles above the UCL differ less than in Fig. 2b, c. We also attribute the smaller $\sigma_{u}$ values above the UCL in CCR06 and CDT07 to the reduced amount of resolved TKE, which decreases with decreasing ratios of $H_{i} /H$ (Grylls et al. 2020). Our simulated vertical profiles of these quantities do eventually decrease with height, as would be expected, for $z/H >$ 5 (not shown). The most notable difference in the profiles between our simulations and CCR06/CDT07 occurs for the turbulence intensity $\sigma_{u}$ in the lower part of the canopy (Fig. 2c) where our model predicts less turbulence than was observed in CCR06 and CDT07. The LES model does not expressly resolve the smallest eddies (i.e., those $< 2dx = H$ in our case), which likely account for a high proportion of the total turbulence in the lower part of the canopy. Figure 2c includes only the resolved turbulence intensity $\sigma_{u}$, and therefore some of the ‘missing’ turbulence is captured by the SGS scheme. Because our mean wind direction is not completely perpendicular to the canopy, there is a small spanwise ($y$) component for each variable that we do not present in these results. However, the spanwise components exist along the homogeneous y-direction with a cyclic boundary condition, and the impact on the processes in the streamwise direction is, therefore, negligible. We also bear in mind that Fig. 2b, c are not entirely direct comparisons. Our results are superficially averaged values throughout a porous medium and the published values are vertical profiles from a single position (position 2) in an array of cubes. Overall, our simulated profiles show a reasonable agreement with the published data, given that we are not simulating flow over cubes, but rather a porous medium of the same mean frontal-area density. Our choice of the $\bar{C}_{d}$ value in the drag parametrization in Eq. 10 appears reasonable in that our simulations do not strongly under- or overestimate any of the parameters U, $F_{u}$, or $\sigma_{u}$ at each height in the domain. Our simulated profiles differ most from the flow over cubes at the bottom of the canopy, where reversed-flow regions can occur around bluff objects, such as cubes, but only occur at quite high densities in porous media (Cassiani et al. 2008). It is also in the lower part of the canopy that our assumption of a height-independent drag coefficient is probably least appropriate, through analogy with flow through arrays of cubes (Cheng and Castro 2002; Kono et al. 2010; Castro 2017), the subject of the studies to which we compared our simulations. Of course, our results and those in the studies cited here are simplified models of urban areas, and therefore can only guide the interpretation of flow through real cities.

4.2 Adjustment of the Flow

As the flow adjusts between the patches of different frontal-area density, we define the adjustment distance $x_{A}$ as proposed by Belcher et al. (2003) and Coceal and Belcher (2004, 2005), as the distance in the streamwise direction at which mean vertical velocity component $W = \left[ {\bar{w}} \right] \approx$ 0 at z = H. Figure 3 shows the perturbation of mean vertical velocity component $W$ by the alternating density patches, with the values of W normalized by $u_{ref}$ = |U| at $z/H$ = 7.5. The first four patches are pictured (Fig. 3a–e), with the streamwise distance normalized by the patch length $L_{p}$. Here $x_{A}$ ≈ 8–9.5H, which can be seen most clearly in the position of the zero contour lines in the 12H and 16H cases (Fig. 3d, e) that occur at x/ $L_{p}$ $\approx$ 0.8 and 0.6, respectively, downstream of each patch edge. The adjustment distance $x_{A}$ is labelled in the first patch of Fig. 3e. In the 8H case (Fig. 3c), the magnitude of $x_{A}$ is visible as the zero contour lines that occur around the end of each patch.

Fig. 3

Effect of patchiness on the mean vertical velocity component, W. a–e Two-dimensional patterns of W in each of the alternating-density cases. Only the first four patches are pictured, with the streamwise distance normalized by patch length $L_{p}$ and the values of W normalized by $u_{ref}$. The dense patches and the sparse patches are bounded by the solid and dashed lines, respectively; f vertical profiles of W at x =3H downstream of each patch. The solid and dashed lines represent profiles at x =3H inside the dense patches and x  = $L_{p}$ +3H inside the sparse patches, respectively; the magnitude of $x_{A}$ is marked on the first patch of e

Another measure of the minimum distance over which the flow adjusts uses the urban canopy length scale $L_{c}$, which can be interpreted as the half-distance the flow takes to accelerate or decelerate (Belcher et al. 2003; Coceal and Belcher 2004). The canopy element drag in Eq. 10 can be expressed as

$${f_{i} = \frac{{\bar{C}_{d} \lambda _{f} }}{H}\left| U \right|\left\langle {u_{i} } \right\rangle = \frac{{\left| U \right|u_{i} }}{{L_{c} }},}$$

(14)

so that $L_{c} = H/\bar{C}_{d} \lambda_{f}$. At the transition of the frontal-area density, the mean velocity component adjusts over a total streamwise distance $x_{A} \approx 3L_{c} { \ln }\left( K \right)$, where $K = \left( {U_{H} /u_{*} } \right)/\left( {H/L_{c} } \right)$ (Coceal and Belcher 2004). Using our results, this gives $x_{A} \approx$ 7.8H and 9.2H for the dense and sparse patches, respectively. This matches the range obtained using our alternative definition surprisingly well given that Coceal and Belcher (2004) derived this formula based on linear-scaling arguments by Belcher et al. (2003) and a very approximate linear regression of the variation of the normalized adjustment distance $x_{A} /L_{c}$ with ${ \ln }\left( K \right)$. Their formula also assumes the incident wind-speed profile is logarithmic, which is not the case here because the presence of the canopy causes the profile to inflect around $z/H$ =1.

4.3 Dynamical Patterns Induced by Adjustment to Density Changes

The alternating-density patches induce clear dynamical patterns in the flow, which affect the exchange between the UCL and the air aloft. Figure 4 shows the perturbation of the mean streamwise velocity component U by the alternating-density patches, relative to the homogeneous control case, with positive and negative values representing higher and lower U than the homogeneous case. At the top of the UCL, the value of U is higher in sparse patches than in the dense (Fig. 4a–e). However, the flow adjusts over a distance $x_{A}$ ≈ 8–9.5H inside the UCL as it moves between patches of different density, as indicated by the location of the maxima and minima in U values for the 12H and 16H cases, for which the mean flow fully adjusts in each patch. The adjustment distance $x_{A}$ is marked on the second patch in Fig. 4e (the 16H case). The patch length $L_{p}$ relative to $x_{A}$ determines the magnitude of the perturbations and the effect of density transitions on the roughness sublayer above the canopy. For the cases we investigated, the magnitudes of the perturbations increase with increasing values of $L_{p}$ / $x_{A}$, and vice versa.

Fig. 4

Effect of patchiness on the mean streamwise velocity component U. a–e Two-dimensional perturbations of U in each of the alternating-density cases. Only the first four patches are pictured, with the streamwise distance normalized by the patch length $L_{p}$ and the values of U normalized by $u_{ref}$. The dense patches and the sparse patches are bounded by the solid and dashed lines, respectively; f mean vertical profiles of U in the dense and sparse patches. The solid and dashed lines represent profiles at x =3H inside the dense patches and x =$L_{p}$+3H inside the sparse patches, respectively; the magnitude of $x_{A}$ is marked on the second patch of e

For the 4H and 6H cases, where $L_{p} < x_{A}$, the flow cannot fully adjust between patches, which means that each patch reflects the characteristics of the previous patch, particularly in the lower half of the canopy (i.e., the value of U is high in the dense patches and low in the sparse). However, for the 8H, 12H and 16H cases, where $L_{p} \ge x_{A}$, the patches are sufficiently extensive for the mean flow to adjust within each patch. This is reflected in the vertical profiles of U (Fig. 4f) where, for the 12H and 16H cases, the value of U is high in the sparse patches and low in the dense at nearly all heights in the canopy. The profiles for the 4H and 6H cases, however, invert at around z =0.5H, while the profile for the 8H case shows a transition between the two extremes, with the magnitudes of the perturbations in U values increasing with the patch size (Fig. 4f).

The adjustment of the streamwise flow affects vertical exchange. Figure 3f shows the peak magnitudes of the vertical velocity component W at the upstream edge of each patch increase with patch size for the 4H and 6H cases, where $L_{p} < x_{A}$. However, the profiles are very similar for the 8H, 12H, and 16H cases, where $L_{p} \ge x_{A}$, because the patches are sufficiently extensive for the mean flow to adjust within each patch, generating strong mean vertical exchange with the UCL. Generally, the value of W is negative (downwards) in and above the sparse patches and positive (upwards) in and above the dense patches, aside from the small reversal caused by the sudden pressure changes at the upstream edge of each patch (Fig. 3a–e). Patch density affects the value of W above the UCL up to around z =5H in the 4H and 6H cases, rising to around z =15H in the 12H and 16H cases.

The alternating density of the urban form also affects the turbulent components of the flow. Figure 5 presents normalized vertical profiles of resolved TKE, [$\bar{e}$], and $F_{u} = \left( { - \left[ {\overline{{w^{\prime}u^{\prime}}} } \right]} \right)^{{\frac{1}{2}}}$ as the perturbation of each variable by the alternating-density cases relative to the homogeneous control case. For a resolved variable $\phi$, profiles of $\phi_{alt} - \phi_{hom}$ are shown, where subscripts alt and hom refer to the alternating-density and homogeneous cases. Online Resources 1 and 2 present two-dimensional plots of perturbations of the parameters $\left[ {\bar{e}} \right]$ and $F_{u}$ analogous to those for the velocity components W and U in Figs. 3 and 4.

Fig. 5

Vertical perturbation profiles of a $\left( {\left[ {\bar{e}} \right]_{alt} - \left[ {\bar{e}} \right]_{hom} } \right)/u_{*}^{2}$ and b $(F_{{u_{{alt}} }} - F_{{u_{{hom}} }} )/u_{*}$. The solid and dashed lines represent profiles at x =3H inside the dense patches and x =$L_{p}$+3H inside the sparse patches, respectively

Figure 5a shows that the resolved TKE $\left[ {\bar{e}} \right]$ is high in the sparse patches and low in the dense, with peak positive and negative perturbations occurring just below z = H. The profiles for the sparse and dense patches are not symmetrical. In all cases, positive perturbations (high $\left[ {\bar{e}} \right]$) in the sparse patches extend through the depth of the UCL, while negative perturbations (low $\left[ {\bar{e}} \right]$) are pronounced only at the top of the UCL. The profiles show slightly higher peak magnitude perturbations for the 6H and 8H cases, because, for those cases, the spatial average includes the strongest perturbations in turbulence that occur at the upstream edge of each patch (Online Resource 1) but only very little adjusted flow, where the perturbations are weaker (as in the 12H and 16H cases towards the downstream edge of each patch). The value of $F_{u}$ is high in the sparse patches and low in the dense patches at most heights within the UCL (Fig. 5b). The height at which the greatest mean perturbations of $F_{u}$ occur decreases with patch size from around z =0.8H for the 4H case to z =0.5H for the 16H case. As previously observed for the TKE $\left[ {\bar{e}} \right]$, the profiles show slightly higher peak magnitude perturbations for the 6H and 8H cases, again due to the proportion of adjusting versus adjusted flow that is captured in the spatial average.

5 Results—Behaviour of Pollutants

5.1 Two-Dimensional Patterns in Pollutant Concentration

Figure 6 presents two-dimensional plots of normalized pollutant concentration $\left[ {\bar{c}} \right]/\left( {QH/u_{ *} } \right)$ in each of the alternating-density cases. We used non-cyclic boundary conditions for the passive scalar. The scalar has a fixed inlet of zero at the upstream edge of the domain and is open at the downstream edge, meaning the scalar concentration increases across the domain. Figure 6 presents the final four patches of each domain where the scalar concentration is highest and by which point the patterns in concentration and fluxes have had the greatest distance to establish themselves downstream of the fixed inlet.

Fig. 6

Effect of patchiness on pollutant concentration. a–e Two-dimensional plots of normalized pollutant concentration $\left[ {\bar{c}} \right]/\left( {QH/u_{ *} } \right)$ in each of the alternating-density cases. The final four patches of each domain are pictured, with the streamwise distance normalized by the patch length $L_{p}$. The dense patches and the sparse patches are bounded by the solid and dashed lines, respectively; f Normalized mean vertical profiles of pollutant concentration in the sparse patches subtracted from that in the dense patches $(\left[ {\bar{c}} \right]_{d} - \left[ {\bar{c}} \right]_{s} )/\left( {QH/u_{ *} } \right)$. Profiles include the entirety of each patch rather than a single x location

For the 12H and 16H cases, where $L_{p} > x_{A}$, concentration maxima occur in the dense patches and minima in the sparse (Fig. 6d, e and vertical profiles in 6f). However, for the 4H and 6H cases, where $L_{p} < x_{A}$, pollutant-concentration maxima occur in sparse patches and minima in the dense patches (Fig. 6a, b and vertical profiles in 6f). This pattern is consistent whether a sparse or a dense patch occurs at the beginning of the modelled domain. In the 8H case, where $L_{p} \approx x_{A}$, maxima and minima occur on the boundaries at the trailing edge of the dense and sparse patches, respectively (Fig. 6c). At the very top of the canopy, the pollutant concentration in the sparse patches is slightly higher than that in the dense patches for all cases (Fig. 6f). We attribute this to the higher turbulent exchange in the sparse patches, meaning the scalar transported from the release height at the bottom of the canopy is more apparent against the low background concentrations at this height.

5.2 Two-Dimensional Patterns in Kinematic Turbulent Pollutant Flux

Figure 7 presents two-dimensional patterns of the kinematic turbulent concentration flux, $F_{c} = \left[ {\overline{{c^{\prime}w^{\prime}}} } \right]$ across the final four patches of each domain. The values of $F_{c}$ are normalized by $u_{*} c_{*}$ where the concentration scale $c_{*} = \sqrt {\left[ {\overline{{c^{\prime}c^{\prime}}} } \right]_{H} }$, corresponding to the square root of the scalar variance at z = H. Figure 7a–e shows the difference between the flux $F_{c}$ in the alternating cases versus the flux $F_{c}$ in the homogeneous case; i.e., positive and negative values indicate regions where the magnitude of $F_{c}$ is, respectively, higher and lower in the alternating case than in the homogeneous case. In each of the alternating-density cases, there is a clear peak in the value of $F_{c}$ at the upstream edge of the sparse patches, at around x =2.5H downstream of the dense-to-sparse transition. Regions of low $F_{c}$ values occur at the upstream edge of the dense patches. The magnitude of the positive perturbations, the regions of high $F_{c}$ values centred at x =$L_{p}$ +2.5H inside the sparse patches, are around three times those of the negative perturbations, the regions of low $F_{c}$ values centred at x =2.5H inside the dense patches (Fig. 7f). The magnitudes of the perturbations are similar for each of the cases, other than the 4H case, for which they are slightly lower. The regions of high and low values of $F_{c}$ extend from z $\approx$ 0.4H within the UCL to z $\approx$ 2H for low values of $F_{c}$ to z $>$ 4H for high values of $F_{c}$.

Fig. 7

Effect of patchiness on the kinematic turbulent concentration flux $F_{c}$. a–e Two-dimensional plots of perturbations of $F_{c}$ in each of the alternating-density cases. The final four patches of each domain are pictured, with the streamwise distance normalized by the patch length $L_{p}$. The dense patches and the sparse patches are bounded by the solid and dashed lines, respectively; f normalized mean vertical profiles of perturbations of the flux $F_{c}$. The solid and dashed lines represent profiles at x =2.5H inside the dense and x =$L_{p}$+2.5H inside the sparse patches, respectively

5.3 Adjustment of Scalar Concentrations and Fluxes

Figure 8 presents streamwise profiles of (a) pollutant concentration and (b) turbulent fluxes of momentum (solid line) and pollutants (dashed line) from the 16H case at z =0.5H. Again, we consider the variables as adjusted once they are in equilibrium; i.e., their values lie within $\pm 5\%$ of those at the downstream edge of each patch, indicated by the two shaded areas in Fig. 8a. Scalar concentration maxima and minima (marked with the vertical dotted lines in Fig. 8a) occur at around x =8H =0.5 $L_{p}$ in the dense patch and x =$L_{p}$ +0.69H =1.43 $L_{p}$ in the sparse patch. Here $L_{p}$ = 16H, so the maxima and minima occur at $x \approx x_{A}$. However, equilibrium is reached slightly later (marked with the vertical dot-dash lines in Fig. 8a), occurring at $x \approx 12H$ = 0.73 $L_{p}$ and $x \approx L_{p}$ +12H = 1.77 $L_{p}$, respectively, for the dense and sparse patches. As turbulent scalar fluxes adjust to flow variations more slowly than momentum fluxes, the turbulence tends to transfer momentum more efficiently than scalar quantities (Kanani-Sühring and Raasch 2015; Li and Bou-Zeid 2019; Ma et al. 2020). This can be seen in Fig. 8b, where the momentum flux has adjusted by $x \approx 10 \;{\text{to}}\; 12$ H = 0.62 to 0.75 $L_{p}$ (marked with the vertical dotted lines in Fig. 8b), but the turbulent scalar flux only fully adjusts around $x \approx 16$ H, although it takes very small absolute values before this point. To summarize Sect. 5.3, maxima and minima of scalar concentration occur at around $x \approx x_{A}$ downstream of the density transitions, and scalar concentration almost fully adjusts by $x \approx 1.2x_{A}$.

Fig. 8

Adjustment of a the normalized scalar concentration $\left[ {\bar{c}} \right]/\left( {QH/u_{ *} } \right)$ at z = 0.5H; and b normalized turbulent fluxes of momentum $\left( { - \left[ {\overline{{u^{\prime}w^{\prime}}} } \right]} \right)/u_{*}^{2}$ (solid line) and scalars $F_{c} /u_{*} c_{*}$ (dashed line) at z = 0.5H. Both fluxes in b are presented as perturbations relative to the homogeneous control case. Only two patches are shown, with streamwise distance normalized by the patch length $L_{p}$. Blue shaded areas show values $\pm$ 5% of those at the downstream edge of each patch

6 Discussion and Conclusions

6.1 Neighbourhood-Scale Flow Regimes

Using these results, we propose two conceptual regimes of neighbourhood-scale flow that emerge from the interaction of the flow dynamics with the patch size. The first is the neighbourhood-ventilation regime, which occurs in neighbourhoods whose frontal-area density varies in patches for $L_{p} > x_{A}$. In neighbourhoods of this type, the mean streamwise velocity component fully adjusts to the change in density for flow between patches, inducing a strong vertical velocity component for $x < x_{A}$ (Fig. 3f) and weak vertical velocity component for $x \gg x_{A}$ (see downstream edge of the patches in Fig. 3d, e). The second is the neighbourhood-percolation regime, which occurs in neighbourhoods whose frontal-area density varies in patches such that $L_{p} < x_{A}$. In these neighbourhoods, the patches are too small for the mean flow to adjust fully within each patch, which constrains the magnitude of the vertical velocity component around the top of the UCL (Fig. 3f). Consequently, the transport of emitted pollutants within one patch is dominated by advection to the next patch. In our case, we take the 4H case as an example of a neighbourhood for which the percolation regime applies, and the 16H case as one for which the ventilation regime applies. Figure 9 presents a schematic of the percolation and ventilation regimes and the competing dynamical processes, illustrating enhanced turbulent vertical exchange in the sparse patches in both regimes. The dominant eddies are always small relative to the patches because turbulent exchange in urban areas is dominated by eddies much smaller than the patches, even in the 4H case (Okaze et al. 2015).

Fig. 9

Inflected wind-speed profile between dense and sparse patches of an urban neighbourhood. Top shows the ventilation regime and the bottom the percolation regime. Turbulent exchange is high in the sparse patches, regardless of patch length

6.2 Pollutant Concentrations

The emergence of the percolation and ventilation regimes profoundly affects the distribution of pollutants. In neighbourhoods subject to the percolation regime, maximum and minimum pollutant concentrations occur in the sparse and dense patches, respectively, which can be seen by comparing the 4H and 16H cases in Fig. 6a, e and f. This is contrary to the common interpretation that ventilation is lower, and the pollution concentration higher, in dense urban patches (Buccolieri et al. 2010). The counter-intuitive behaviour arises in neighbourhoods with the percolation regime (i.e., density changes in small patches) because turbulent exchange is suppressed in the dense patches. Therefore, pollutant accumulates over the length of the dense patch and is advected into the neighbouring sparse patches. Conversely, in the sparse patches, turbulent exchange is high and the patches better ventilated (Fig. 7f), meaning less pollutant is advected downstream into the neighbouring dense patches.

However, in neighbourhoods with the ventilation regime (i.e., density changes in larger patches), maximum pollutant concentrations occur in the dense patches (Fig. 6e, f). The exact locations of the maxima and minima in the ventilation regime reflect the importance of both the mean flow and turbulent exchange. The turbulent pollutant flux is high at the upstream edge of the sparse patches and low at the upstream edge of the dense patches (Fig. 7). Therefore, pollutant-concentration minima at $z \approx 0.5$ H occur at shorter distances downstream of the edge of sparse patches (Δx $\approx$ 5.5H) compared with maxima downstream of the edge of the dense patches (Δx $\approx$ 7.5H) (Fig. 10). In neighbourhoods where the frontal-area density varies in patches where $L_{p} \approx x_{A}$ (the 8H case), we see a smooth transition between the two regimes. Pollutant-concentration maxima and minima occur around the leading edge of sparse and dense patches, respectively. The mean flow is perturbed more than in the percolation regime, but less than in the ventilation regime.

Fig. 10

Location of concentration maxima and minima at $z = 0.5$ H downstream of the leading edge of the dense (red) and sparse (grey) patches. Values of $x/L_{p}$ $>$ 1 indicate maxima/minima downstream of the transition occurring in the neighbouring patch (for example, scalar maxima in the sparse patches in the 4H case)

6.3 Limitations and Effect of Other Variables

Here we present a quasi-two-dimensional model of flow in an inhomogeneous urban area, by altering a single variable (frontal-area density), and the streamwise distance over which the density changes occur. The wind direction was held constant, and therefore our model does not capture the variations that the wind direction can induce in scalar fields even in idealized models of urban areas (Philips et al. 2013). For example, we would not expect our results to hold for flow parallel to the alternating-density patches in Fig. 1. In real urban areas—particularly those with deep street canyons—the wind direction and urban form interact, for example, regions of high ventilation may occur at different locations depending on the prevailing wind direction (Yim et al. 2014). We also assume a constant mean obstacle height. In real urban areas, the heights of the buildings and the other obstacles vary, which affects the distance over which the flow adjusts from one patch or neighbourhood to another (Coceal and Belcher 2005) and the rate at which the UCL is ventilated (Hang et al. 2012). We also approximate the drag coefficient $\bar{C}_{d}$ as constant (Sect.  3.2), which would not be a suitable simplification for a very tall canopy, for which a more empirical approach based on a real urban structure is likely to be more appropriate (Gutiérrez et al. 2015). We considered only neutral atmospheric conditions. Atmospheric stability can profoundly affect the flow, even at the neighbourhood scale; for example, in unstable conditions, scalars are transferred more efficiently with increasing instability (Wang et al. 2014). Further, while our model provides a helpful way to approximate the dynamics of air pollutants at the neighbourhood scale, careful consideration of the urban form and the behaviour of pedestrians, such as their mode of transport, is needed to assess the exposure of the population to those pollutants (Kingham et al. 2013; Tenailleau et al. 2015).

Our approximation of pollutants as a homogeneous source of passive scalar may also be inappropriate for certain applications, such as for investigating species with short lifetimes, which may require a coupled chemistry scheme (Zhong et al. 2015), or for considering species with very heterogeneous sources and sinks, such as water vapour and CO₂, which behave very differently even in the same neighbourhood under the same conditions (Nordbo et al. 2013; Ramamurthy and Pardyjak 2015). In real urban areas, emission rates of even passive pollutants also vary spatially, depending on factors such as the location of busy roads. This patchiness in pollutant emissions affects concentrations within the UCL. As an illustration, we simulated a further case with the same urban form as the homogeneous control case, varying the emission rate of a passive pollutant between high (Q =1.5 μg m⁻³ s⁻¹) and low (Q =0.5 μg m⁻³ s⁻¹) rates in patches with $L_{p} = 16H$. Figure 11 shows the normalized concentration $\left[ {\bar{c}} \right]/\left( {QH} \right)/u_{ *}$ at z = 0.5H. Because the urban area is homogeneous in form, there are no dynamical patch effects induced by density transitions as in the previous subsections. However, the concentration still adjusts as it moves between patches of different emission rates. The adjustment occurs over a distance of around $x_{A}$, the point at which the concentration comes within 5% of that at the downstream edge of each patch (marked with the vertical dotted lines in Fig. 11). Many cities worldwide have implemented measures to improve air quality, creating, for example, low emission zones in which the traffic of older vehicles is limited, with the inevitable leakage of pollutants from surrounding areas. The adjustment distance $x_{A}$ appears to provide an approximate distance over which this contamination occurs, even in the absence of dynamical effects induced by changes in density.

Fig. 11

Adjustment of pollutant concentration between patches of different emission rates into an area with a homogeneous urban form. Only two patches are shown, with the streamwise distance normalized by the patch length $L_{p}$. Solid line shows normalized concentration $\left[ {\bar{c}} \right]/\left( {QH} \right)/u_{ *}$ at z = 0.5H, blue shaded areas show concentrations $\pm$ 5% of those at the downstream edge of each patch, and grey shading and no shading, respectively, show the high- and low-emission patches

6.4 Applications and Outlook

While mean pollutant concentrations and fluxes are spatially highly variable between and within different parts of a city (Figs. 6, 7), in patchy neighbourhoods subject to the ventilation regime (i.e., where the frontal-area density varies in large patches), point observations at a distance $> x_{A}$ downstream of a density transition in the direction of the prevailing wind direction should represent pseudo-equilibria well. The adjustment distance $x_{A}$ therefore represents an approximate lower bound for fetch requirements for ‘urban background’ micrometeorological and flux observations. Further, our results suggest neighbourhood ‘hot spots’ of pollution, such as contamination downstream of a small dense patch with heavy traffic, may not be captured by relying on only roadside or ‘urban background’ observations.

A porous-medium approach to neighbourhood-scale processes provides a good balance between computational efficiency and flow resolution, provided the porous model is formulated to capture the characteristics of the urban area in question. The percolation and ventilation regimes (a) capture the mean and turbulent components of the transport terms for pollutants and other scalar quantities; and (b) can be described using a range of geometrical variables simple enough to be deployed in high-Reynolds-number flow models. Further, in the quasi-two-dimensional case we present here, the broad characteristics of the regimes are insensitive to the frontal-area density in the quantity $f_{i}$ and the SGS sink in Eq. 11. A greater density contrast between the sparse and dense patches would lead to more pronounced examples of the regimes but would not materially change the adjustment distance $x_{A}$ either side of which the regimes emerge. A lesser contrast would simply blur the transitions between sparse and dense patches, eventually approaching flow in a homogeneous urban area. The effect of stochastic porosity and patch-scale changes on the flow dynamics would depend on the patch size. In large patches, once the flow has adjusted, it has very little ‘memory’ of the preceding urban form. To approximate flow dynamics in neighbourhoods with small patches, we would need to know the size of the patch of interest and its density relative to the neighbouring patch upstream.

Neighbourhood-scale processes provide potential boundary conditions to couple micro-scale urban-transport models and mesoscale meteorological forecasts (Xie 2011). However, parametrizations of neighbourhood-scale processes must be flexible enough to capture the effect of inhomogeneity in the urban form on flow and dispersion. Here, the specifics of the urban form were captured entirely in the drag-force parametrization in Eq. 10. This approach adapts LES models previously used for vegetation canopies, with the representation of the urban form informed by porous models of cities, specifically a vertically homogeneous frontal-area density $\lambda_{f}$ and a higher drag coefficient than in vegetation models. Further in situ observations and experiments in wind tunnels (preferably both) would help to determine the extent to which this approach holds for real urban areas. For example, we suspect our approach would be more appropriate for irregular cities comprising vegetation and obstacles of various shapes and sizes, rather than regularly spaced buildings with few other features, for which a vertically varying value of $\lambda_{f}$ or the drag coefficient may yield better results. Further testing may allow the percolation and ventilation regimes to be adapted to link between neighbourhood and city/mesoscale models. Urban air quality motivated this study, and therefore we concentrated on scalar quantities representing pollutants emitted near the ground in urban areas. However, with modifications, this framework could be used to model processes in other environments, such as forests, or other scalar quantities, such as temperature. For example, in forests, after adding a scalar sink throughout the canopy, this framework could be used to investigate the transport of pollutants, such as ozone, into a patchy forest (Ma and Liu 2019; Ma et al. 2020). This type of investigation is extremely difficult using observations, even with an extensive network of sampling equipment (Aubinet et al. 2010).