Fast Models for Predicting Pollutant Dispersion inside Urban Canopies

: A fast pollutant dispersion model for urban canopies is developed by coupling mean wind proﬁles to a parameterisation of turbulent diffusion and solving the time-dependent advection– diffusion equation. The performance of a simpliﬁed, coarse-grained representation of the velocity ﬁeld is investigated. Spatially averaged mean wind proﬁles within local averaging regions or repeating units are predicted by solving the three-dimensional Poisson equation for a set of discrete vortex sheets. For each averaging region, the turbulent diffusion is parameterised in terms of the mean wind proﬁle using empirical constants derived from large-eddy simulation (LES). Nearly identical results are obtained whether the turbulent ﬂuctuations are speciﬁed explicitly or an effective diffusivity is used in their place: either version of the fast dispersion model shows much better agreement with LES than does the Gaussian plume model (e.g., the normalized mean square error inside the canopy is several times smaller). Passive scalar statistics for a regular cubic building array show improved agreement with LES when wind proﬁles vary in the horizontal. The current implementation is around 50 times faster than LES. With its combination of computational efﬁciency and moderate accuracy, the fast model may be suitable for time-critical applications such as


Introduction
Fast and accurate prediction of pollutant dispersion is important for practical applications such as operational air-quality modelling and emergency dispersion modelling. Computational fluid dynamics (CFD), which has been extensively applied to pollutant dispersion in idealised [1] and realistic [2,3] urban areas, features relatively high accuracy; however, it is too computationally demanding for time-critical applications. Hence, there is a need for a fast dispersion model that can provide useful predictions without significant computational resources, extensive training data or many adjustable parameters. The simplest dispersion model is the Gaussian plume model (GPM) [4], which is still widely used. In the original formulation, the pollutant concentration downwind of a point source is estimated by assuming flat terrain and an eddy diffusivity that increases linearly with downwind distance [5]. The GPM is therefore not well-suited to the urban canopy layer, where pollutants are emitted and many people live and work. The intrinsic limitations of the GPM have helped spur the development of fast dispersion models that are more appropriate for urban areas. The Quick Urban & Industrial Complex (QUIC) model is based on a variational solution for a steady velocity field [6] and a Lagrangian stochastic model for the time-dependent dispersion [7]. Although QUIC has been successfully applied in a variety of contexts (e.g., [8]), accuracy is potentially limited by the initial guess for the velocity and the nonlocal mixing parameterisation. The urban dispersion models SIRANE [9] and MUNICH [10] assume constant concentration on the building scale (i.e., within a street canyon) and parameterise the mean velocity profiles and turbulent transfer. They show good agreement with in situ observations [11,12]; however, information about individual buildings is lost as only the canyon dimensions and aspect ratio are considered.
The complexity of the urban environment is such that any fast dispersion model must rely on assumptions, approximations and empirical relationships. Hence, models such as QUIC, SIRANE and MUNICH may not be well-suited to all applications or geometries. In this paper, we introduce a new fast dispersion model based on coarse graining of the velocity field, i.e., spatial averaging over subregions or repeating units. This is a potentially attractive strategy because accurate and computationally efficient models of spatial averages are more easily obtained than ones for the fine-scale turbulent velocity field. We estimate the spatially averaged mean wind profiles within the urban canopy by approximating the vorticity field [13] and couple them to different turbulence parameterisations. Although the time-dependent three-dimensional (3D) advection-diffusion equation is solved, mean wind profiles within subregions are used in place of a continuously varying 3D velocity field. For one of the parameterisations, turbulent fluctuations within each subregion are prescribed; for the other one, the associated effective diffusivity is used instead.
This study introduces a new multiscale fast dispersion model based on estimates of the spatially averaged wind profiles and applies it to the canonical case of a cubic building array. The fast dispersion model (FM) is presented in Section 2: the estimation of the mean wind profiles within the array units is reviewed and the two turbulence parameterisations are described. The methodology section details the implementation of two versions of the FM, which are distinguished by their turbulence parameterisation, as well as the LES model that is used for calibration and assessment (Section 3). The performance of the fast model is assessed in Section 4 through comparisons with LES and the GPM. The sensitivity to mean wind speeds and the specification of the diffusivity tensor is described in Sections 5 and 6, respectively. Limitations of the fast dispersion model are discussed in Section 7. Practical implications for operational dispersion modelling are given in Section 8.

Derivation of the Fast Dispersion Model
The fast dispersion model requires mean wind profiles, a turbulence parameterisation and a numerical model for the solution of the advection-diffusion equation. It includes the following steps: 1.
Estimation of mean wind profiles using the vortex method [13].
• FM1 (direct approach): velocity perturbations are specified using the mean wind profiles. • FM2 (indirect approach): the effective diffusivity is specified using the temporal velocity variance.

3.
Solution of advection-diffusion equation using the mean wind profiles and turbulence parameterisation.
Mean wind profiles within urban canopies can be estimated by assuming that the velocity is controlled by intense layers of vorticity, i.e., by solving the three-dimensional Poisson equation for a set of discrete vortex sheets. The detailed procedure is summarised in Section 2.1. Briefly the method yields steady, spatially averaged wind profiles, where i = {x, y, z}, · R l denotes a volume average over a region R l (e.g., a canyon) and the overbar is a time average. Henceforth the overbar is dropped for convenience. Since urban dispersion is strongly affected by the turbulent component of the flow, a turbulence parameterisation is required. Two versions are examined. First, turbulent fluctuations are introduced via stochastic velocity perturbations: this is referred to as FM1.
Second, the effect of turbulent fluctuations is mimicked by defining an effective diffusivity: this is referred to as FM2. Note that the effective diffusivity [14] differs from standard prescriptions of the turbulent diffusivity.
The numerical solution of the advection-diffusion equation is described in Section 3.2.

Estimating Mean Wind Profiles
The mean wind profiles are estimated by the vortex method [13], in which a set of vortex sheets (i.e., thin layers of intense vorticity) are used to determine a set of basis functions. This method has been successfully applied to a variety of urban canopies, including realistic ones, by [15].
A brief summary of the vortex method now follows. For further details, see [13,15].

1.
Definition of the vortex sheets. The vorticity field is strongly localised near solid walls for urban canopy flows. Therefore, the continuous vorticity field is approximated by a set of discrete vortex sheets at the walls, namely, at the top of the tallest building ( t ), street level ( b ), or along side walls ( s ). Approximating the vorticity field in this way is a standard technique in vortex dynamics [16,17] that has been applied to different areas of fluid dynamics [18,19]. A given vortex sheet location may not include all three vorticity components. A vorticity sheet is defined by the vortex sheet location j and vorticity component ω k . The set of vorticity sheets, {( j , ω k )}, for a region R l is denoted by Ω l .

2.
Solution of Poisson equation. The flow induced by each vortex sheet (a Green's function or, in two dimensions, streamfunction) is obtained by solving a 3D Poisson equation. More precisely, velocity basis functions, V j,ω k i R l , are obtained for the entire set of vorticity sheets. The basis functions represent the effect of a specific vorticity sheet for the specified building geometry.

3.
Synthesis. The basis functions are used to estimate mean wind profiles. Summing over vorticity sheets, yields the mean velocity profiles for region R l . The α j,ω k are weights that determine the contribution of each vorticity sheet to the estimated velocity profile.

4.
Calibration. The weights are calculated through a calibration procedure. Using reference CFD data, the α j,ω k are defined to be proportional to the strength of each vorticity sheet, i.e., where γ j,ω k is the circulation of the vorticity sheet. C j,ω k is a geometric constant, which depends on the geometry only. The C j,ω k are obtained by minimising the residual between the reference CFD profile and the estimated profile (i.e., Equation (1)). By construction, the C j,ω k are calculated for a reference wind direction; however, they can be successfully applied to other wind directions for various urban canopies, justifying their interpretation as geometric constants [15]. The circulation for the reference wind direction is calculated from the reference mean wind profiles; γ j,ω k for other wind directions is obtained by rotating the mean velocity.

5.
Matching. Since inviscid vortex dynamics is assumed, the preceding steps only yield an estimate of the mean velocity profiles away from boundaries or above the ground. The vertical profiles near the ground, i.e., below a height z b , are parameterised with a log-law profile and the friction velocity.

Turbulence Parameterisations
The velocity components can be decomposed with respect to the space-time average [20,21], i.e., where i ∈ {x, y, z}, denotes the velocity component and R l label different subdomains or averaging regions. Above the canopy, only a single region is required. The tilde denotes fluctuations with respect to the spatial mean. In turbulence parameterisations, spatial fluctuations are rarely considered. Strictly speaking, this is justified only for homogeneous turbulence; however, it is reasonable to assume that ū i has no preferred spatial structure, in which case u i effectively represents temporal and spatial fluctuations.
Turbulent fluctuations are parameterized in terms of the mean flow by discretising in space and time. At times t ∈ [t k , t k+1 ], where t k ≡ kT update,i , k ≥ 0 is an integer, U R l i (z) ≡ ū i R l are the mean profiles obtained from the vortex method, β R l i is a turbulence coefficient, and N(0, 1) is a normal distribution with mean 0 and standard deviation 1. The spatially random fluctuations are essentially frozen between updates but temporally uncorrelated from one update to another.
Although specifying turbulent fluctuations explicitly is unconventional in CFD, Lagrangian stochastic models for turbulent diffusion work similarly [22]. We assume that the turbulent fluctuations in each direction are linearly related to the corresponding mean velocity components. As with the mean wind profiles, the β R l i are calculated for a reference wind direction and applied to other directions.
The turbulence coefficients are determined from reference LES velocity data, v i , rather than from a theoretical model because this study is concerned with the extent to which simple turbulence parameterisations may be applied to fast dispersion modelling. Vertical profiles are calculated as follows: where the temporal root-mean square (rms) is taken in the numerator and denominator. For simplicity, values inside and above the canopy, β R l i,canopy and β i,above , respectively, are considered. They are defined as where L z denotes the height of the computational domain. Turbulent diffusion is introduced by Equation (5). Its magnitude depends on the timescale at which the fluctuations are updated, T update,i , which may be interpreted as a correlation timescale for the turbulence. Since pure diffusion is associated with white noise (see Section 2.2.2), it is natural to use the shortest timescale resolved by the model, namely, the model timestep. Thus T update,i ≡ ∆t in the first instance.

Fast Model 2 (FM2)
Parameterising the effect of the turbulent fluctuations on the scalar field is more convenient. This can be achieved by appealing to standard results for stochastic differential equations (SDEs) to derive an effective diffusivity [22]. Consider the SDE where x i is the displacement vector, dW j is an incremental Wiener process and κ ij is the noise amplitude. This SDE is equivalent to the advection-diffusion equation if the concentration C is interpreted as the probability distribution function of x i [14,23]. Physically, κ ij is an effective diffusivity tensor that represents the effects of temporal fluctuations absent from the velocity, U i . The preceding theory is idealised and cannot be applied exactly to urban flows. First, it is assumed that the fluctuations are Gaussian and completely independent of the mean flow. Second, the boundary conditions at the wall are neglected. Nevertheless, this theory represents a natural starting point for estimating the effective diffusivity. Assuming that turbulent fluctuations can be modelled by a stochastic process (i.e., by κ ij dW j ), the vertical profile of the effective diffusivity within each repeating unit is estimated from reference LES velocity data: where the temporal standard deviation of the Reynolds stress, and T R l ij is the corresponding correlation time. In the first instance, by analogy with molecular diffusion. As with FM1, T R l ij ≡ ∆t and averages inside and above the canopy, κ R l ij,canopy and κ ij,above , are taken: The advection-diffusion equation is written in dimensional form so as to facilitate comparison with the SDE, Equation (8). Nevertheless, the effective diffusivity does depend on the choice of length and velocity scales. Nondimensionalising Equation (9) yields where the Peclet number Pe = UL/κ ij , with U and L being characteristic velocity and length scales, respectively. Hence, the effective diffusivity can be scaled by fixing Pe. Alternatively, the nondimensional diffusivity, is fixed.

Methodology
There are many models for predicting pollutant dispersion [24]. The most widely used theoretical model is the Gaussian plume model (GPM) [5], which incurs essentially no computational cost, neglects the presence of buildings but works well in the far field [1]. The GPM is summarised in Appendix A. The most accurate approach is building-resolving CFD, but it is technically complicated and computationally expensive, especially for LES [25]. LES and GPM results will be compared to the fast dispersion models derived in Section 2.

Computational Domain
A regular cubic array, which is more representative of urban areas than the 2D street canyon [26], constitutes the topography for this study. For a sufficiently large array, the flow inside the canopy is approximately homogeneous. Hence, there are three types of repeating units: channels (R1), where the buildings are aligned in the y direction; canyons (R3), where the buildings are aligned in the x direction; and intersections (R2), which lie between them. See Figure 1 for an illustration.
Properties of the mean flow and turbulence differ qualitatively among the repeating units [21]. For example, for a flow parallel to the x axis (i.e., at 0°), there is approximately rectilinear motion through the channel and intersection units, but a pair of counterrotating vortices inside the canyon units [27]. The flow patterns change with the wind direction, θ, but the inhomogeneity persists [28]. Hence, the repeating units are used to define the regions over which the mean velocity profiles are calculated (see Equation (2)).
The topography extends throughout the computational domain. A plan view is shown in Figure 1. It has dimensions L x × L y × L z or 24H × 24H × 6H, where H = 20 m is the length of each side of the cubic buildings. The same computational domain is used for the LES and FM.

LES
The fast models were calibrated and assessed using building-resolving, implicitly filtered, LES and OpenFOAM, a finite-volume library for the Navier-Stokes equations [29]. For all simulations described herein, the model configuration is as follows. The advection terms are discretised using a Gauss scheme and second-order Gauss linear upwinding is used for the velocity divergence. The Crank-Nicolson scheme is used for time-stepping. The subgrid scale (SGS) stress tensor, τ ij , is parameterized by the wall-adapting local eddy-viscosity (WALE) model [30], ) is the resolved-scale strain rate tensor. The Spalding [31] wall function (implemented in OpenFOAM as nutUSpaldingWallFunction) is applied to solid surfaces: it specifies the turbulent viscosity as a function of the wall-normal distance by correcting the shear stress at the wall. An isotropic grid spacing ∆ = 1 m is used away from the walls. Mesh refinement is applied in the immediate vicinity of the walls where the finest horizontal resolution is ∆ = 0.5 m. The flow is forced by a constant external forcing. A wind speed of 3 ms −1 is obtained at the upper lid. Periodic boundary conditions are applied in the horizontal. A free-slip boundary condition is imposed at the top of the computational domain and there is no slip on all solid surfaces.
The passive scalar is released at ground level (actually the centre of the first cell) from a steady point source (P3) located inside the canyon unit (R3). To avoid the accumulation of scalar, a sponge layer was added to the lateral boundaries: the concentration is manually reset to zero after each time step. The turbulent Schmidt number is set to Sc t = 0.7 [1,32]. Since the turbulent diffusivity is several orders of magnitude greater than the molecular diffusivity, the contribution from the latter is neglected. The calibration relies on velocity statistics only (Section 3.3.3), but passive scalar statistics are needed to evaluate the performance of the fast models (Section 4).
The model was spun up for approximately 4000 s, which is sufficient for the mean flow and turbulence to reach approximate statistical stationarity; after the spin up, data were collected for another 1000 s. The variable time step ∆t LES ∼ 0.2 s and the average CFL number is less than unity at every time step. Henceforth, only time-averaged data are considered. For convenience, the overbar is omitted for the mean wind profiles, U R l i . The validation of the LES model against wind-tunnel data for velocity and scalar statistics are presented in Appendix B. Validation metrics for the mean streamwise velocity and pollutant concentration indicate acceptable agreement. We conclude that the LES model is capable of providing a useful baseline truth for the fast model.

Numerical Solution
The numerical solution of the advection-diffusion equation is essentially identical to that for the full LES model (Section 3.2). The numerical schemes, boundary conditions and configurations are identical. The (variable) timesteps are, ∆t FM ∼ ∆t LES , and the integration lengths are identical (i.e., 5000 s). The only difference is that the velocity field is prescribed (obviating the need for a spin up) and the turbulent diffusivity is set to zero. FM1 and FM2 use the same mean velocity profiles, U R l i (Section 3.3.2), but differ in their turbulence parameterisations.

Mean Wind Profiles
The first step in solving for the mean wind profiles, U R l i , is to define the vorticity sheets. From the vorticity field at 0° (Figure 2), strong vorticity appears at the roof, street levels and sidewalls, as with street canyons [13] and realistic canopies [15]. This suggests the following arrangement of vorticity sheets: where Ω i denotes the set of vorticity sheets used to predict velocity component u i . Note that each Ω i includes all vorticity components that influence velocity component u i : by definition, ω i has no effect on u i . The procedure can be simplified by taking a subset of the vorticity sheets [15], but this is not necessary to test the usefulness of the FM. The flow induced by each vorticity sheet is obtained by solving the associated Poisson equation with a geometric-algebraic multi-grid solver and a free-slip boundary condition at the walls. For all cases considered herein, the calibration was performed using data for 45°; hence, the geometric coefficients for 45°are applied to other directions (e.g., 0°).
At z/H = 0.25, the estimate from the vortex method is matched to a log profile that extends to the ground.
The vortex method performs well for idealised [13] and realistic [15] domains. Here we summarise its performance for the cubic building array. Figure 3 compares estimated and LES wind profiles for θ = 45°. There is generally good agreement for the velocity components, though relative errors are larger for the vertical velocity on account of the small magnitudes (typical dimensional errors are ∼0.01 ms −1 ). For the horizontal velocity components, relative errors inside the canyon are ∼10%. Agreement is maintained up to z/H = 2, but for clarity, the figure extends to z/H = 1 only.

Turbulence Parameterisation
The numerical solution of the advection-diffusion equation is essentially identical to that for the full LES model (Section 3.2). The numerical schemes, boundary conditions and configurations are identical. The (variable) timesteps are comparable. The only difference is that the velocity field is prescribed (obviating the need for a spin up) and the turbulent diffusivity is set to zero. FM1 and FM2 use the same mean velocity profiles, U R l i (Section 3.3.2), but differ in their turbulence parameterisations.

Fast Model 1 (FM1)
The turbulence coefficients, β R l i,canopy and β i,above , were calculated from reference LES data using Equations (6) and (7). Analogously to the estimation of the mean wind profiles in Section 2.1, reference LES data for 45°were used. The β R l i,canopy and β i,above so obtained (Table A6) are then applied to other wind directions, e.g., 0°. The sensitivity to the choice of reference wind direction is considered in Section 5.

Fast Model 2 (FM2)
The effective diffusivities, κ R l ij,canopy and κ ij,above , were determined from reference LES data using Equations (10) and (13). As with the turbulence coefficients, reference wind data for 45°were used to calculate the effective diffusivities, which are applied to 45°and 0°. Normalising the effective diffusivities via Equation (15) to yield κ R l * ij,canopy and κ * ij,above (Table A7) facilitates the application to stronger or weaker mean flow. The reference velocity scale is U ≡ U ref = 3 ms −1 (i.e., free-stream speed at the upper lid, 5H) and the reference length scale L = 20 m (i.e., building height, H). Since the choice of 45°for the reference data is arbitrary, the sensitivity of FM2 to the wind direction is examined in Section 6.

LES and GPM
The performance of the fast models is now evaluated against LES passive scalar data. For simplicity, a vertical average is taken. Figure 4a,b show mean concentration fields, C H from LES. The spatial patterns are consistent with many previous numerical studies (e.g., [27,28,33]). At θ = 0°, pollutants preferentially disperse in the x direction, i.e., along the rows, and are easily trapped within the canyon units. The plume width does not increase noticeably downstream of the source, likely on account of spanwise channelling of pollutants. At θ = 45°, there is more crossstream dispersion. Predictions from the GPM are shown in Figure 4c,d. The GPM does a poor job of capturing the plume structure and generally underestimates concentrations as the cross-stream dispersion in the direction is overestimated. Since the GPM does not include a representation of buildings, trapping within canyons is not captured. Hence, the dependence on the wind direction is also missed.

FM1
Scalar fields generated by FM1 are plotted in Figure 5. The same turbulence coefficients and geometric constants (obtained at 45°) were used for both wind directions, i.e., 0°and 45°. Compared to the GPM ones, the FM1 fields inside the canopy show much better agreement with LES: in particular, the strong downwind and weak cross-stream dispersion are partially captured. However, the cross-stream dispersion is too weak in FM1 as the localisation of the plume along the 45°line is exaggerated. This suggests that turbulent diffusion is too weak in FM1. To analyse the accuracy of FM1 as a function of the downwind distance, the concentration is averaged over the canopy and along the spanwise direction: Lateral profiles, i.e., C H,y versus x ≡ x − x s , where x s is the streamwise location of the point source, are plotted in Figure 5c,d. The LES and FM1 profiles show reasonable agreement, though the step-like changes in the concentration between units are smoothed out, especially at 0°. The GPM underestimates the concentration everywhere. Far from the source, the nature of the dispersion should become independent of the precise arrangement of buildings since pollutants will have passed many buildings, leading to central-limit-type behaviour. This idea is captured in the notion of the radius of homogenisation (RAD) [34], which demarcates the near-(x < RAD) and far-field (x > RAD) regions. The RAD can be estimated from concentration autocorrelations [1].
Statistical performance measures (Appendix C) for different averaging regions are listed in Table 1. Compared to the GPM, FM1 shows significantly smaller errors; for example, the normalised mean-square error (N MSE) for the entire field decreases by ∼60-90% inside the canopy. However, the spatial dependence differs. Whereas the far-field error for the GPM is always smaller, this is not necessarily true for FM1 as the N MSE at 0°actually increases slightly in the far field. On the other hand, FM1 at 45°performs rather poorly in the near field as the correlation coefficient (R) is close to zero. Above the canopy, the N MSE is significantly smaller but the other metrics do not show obvious improvements. The preceding results may be explained as follows. For the GPM, the accuracy of a diffusion equation for modelling the pollutant dispersion improves in the far field. Inside the canopy, FM1 is similar to the GPM insofar as dispersion is driven by the mean flow and a turbulence parameterisation; however, the GPM assumes a constant streamwise velocity and a turbulent diffusivity that is a function of x (Appendix A). At 0°, the FM1 approximation of the mean flow by the averages for each repeating unit evidently works well as the near-and far-field N MSE values are similar, but at 45°, the flow within the canyon and intersection units is of greater importance (e.g., see the horizontal streamlines in [1]). The N MSE is much smaller above the canopy, likely on account of the greater homogeneity of the scalar field. Figure 6 shows scalar fields and lateral profiles inside the canopy for FM2. The plume structures (panels a,b) are very similar to those for FM1, although cross-stream dispersion at 45°increases slightly for FM2, in agreement with the LES. The lateral profiles (panels c,d) show that, like FM1, FM2 underestimates concentration differences between the units.

FM2
Statistical performance measures (Table 1) confirm that FM2 is more accurate than the GPM inside and above the canopy. For example, the N MSE inside the canopy is ∼70-80% smaller for FM2. The N MSE decreases in the far field at 0°and 45°for FM2.

Comparison of FM1 and FM2
FM1 and FM2 yield similar results. Although there are slight differences in the plume structures at 45°, the statistical performance metrics do not indicate that one model always performs better than the other. According to the N MSE inside the canopy, FM1 is ∼40% better at 0°, but ∼30% worse at 45°. Since FM1 explicitly parameterises turbulent fluctuations in terms of the strength of the mean flow, it is plausible that it should perform worse when turbulent fluctuations are weaker, i.e., at 45°(see [1] for analysis of the dependence on the wind direction). Since FM2 is based on a diffusion equation, its performance should be largely insensitive to the strength of the turbulent fluctuations. The differences between FM1 and FM2 are also seen in the near-and far-field errors: whereas the near-field FM1 errors are smaller at 0°, the near-field FM2 errors are larger, as with the GPM.

Wind Speed
The results of Section 4 were obtained using LES data for a specific reference wind speed. We now assess whether the geometric constants and the turbulence coefficients for the baseline configuration, i.e., U ref = 3 ms −1 (at z/H = 1), are applicable to other reference wind speeds, namely U ref,low = 1 ms −1 and U ref,high = 9 ms −1 .
The performance of the vortex method at low and high wind speeds is compared using a spatial average within each repeating unit, i.e., where A R i is the area of R i . Figure 7a shows   The dependence of the turbulence coefficients on the reference speed is quantified by recalculating them (Equation (6) (Table A3) and lateral profiles ( Figure A3). Applying the geometric constants and turbulence coefficients for U ref = 1 ms −1 to different reference wind speeds is justified.
The performance of FM2 is comparable. For U ref,low and U ref,high , the effective diffusivity is scaled by fixing the Peclet number (Equation (15)). Statistical metrics (Table A3) and lateral profiles ( Figure A4) confirm good agreement with the LES results for U ref = 1 ms −1 .

Mean Flow Only
The reasons for the superior performance of the fast model relative to the GPM ( Figure 5) are not immediately obvious. The GPM assumes an idealised mean flow, i.e., a constant streamwise velocity, as well as a specific form for the eddy diffusivity (Appendix A). The contribution of the mean flow is assessed by comparing the performance of the fast model with different prescriptions of the mean flow but no turbulence parameterization The prescriptions include: (i) constant (height-independent) velocity, (U 0 , 0, 0); (ii) uniform mean velocity profile for the entire domain, U i (z); and (iii) uniform mean velocity for each repeating unit, U R l i (z). In the first case, U 0 = U ref , and the reference wind speed at z/H = 5, as with the GPM; in the latter two cases, the mean profiles were obtained using the vortex method. The abbreviation 'FM' is used in this section because FM1 and FM2 are equivalent when the pollutant dispersion is driven by the mean flow only. Figure 8 compares spatially and temporally averaged concentration fields at 0°for the different cases. In all cases, the plumes are too narrow compared to LES; this can be attributed to the neglect of turbulent diffusion, which causes the scalar to spread in directions other than that of the mean flow. For the constant velocity (panel a), downstream variations are relatively weak inside the canopy. For the uniform velocity profile (panel b), the lateral profile (panel d) shows a slight asymmetry between the leeward and windward walls and little downstream decay. Allowing the mean profiles to vary between repeating units (panel c) leads to higher concentrations near the windward wall, in agreement with the LES (Figure 4). Moreover, concentrations near the centre of the plume decrease more rapidly than in cases with a horizontally uniform velocity field. Nevertheless, the spanwiseaveraged lateral profiles still decay too slowly.
Retaining horizontal dependence in the mean velocity field is essential if satisfactory performance is to be obtained from the fast model. Statistical metrics (Table A4) confirm that inside the canopy the FM agrees best with LES when U R l i (z) is used: although the N MSE in this case is O(1), which is large but acceptable according to standard air quality acceptance criteria [35], it increases by two to three orders of magnitude larger when the horizontal dependence of the mean flow is neglected. These results are consistent with those of [36], who found that the pollutant dispersion from a unit-aspect-ratio street canyon is largely controlled by the mean flow.

Sensitivity to the Turbulence Parameterisations
The standard configurations of FM1 and FM2 may be not optimal. The sensitivity to the specification of the turbulent fluctuations and effective diffusion is examined below.

FM1
Two key choices were made in the definition of the default FM1 configuration. First, we chose to update the turbulent fluctuations in FM1 at every timestep, with ∆t FM ∼ 0.2 s being equal to the LES timestep, ∆t LES . The validity of this prescription is not guaranteed as the appropriate turbulence timescale may differ from the numerical timestep. From the (temporal) autocorrelations, the timescale of the fluctuations, T R l i , may be estimated from the e-folding timescale, (Values of T R l i are listed in Table A8). Second, we chose the spatially averaged turbulence coefficients, β R l i,canopy and β i,above , for the canopy and overlying region, Equation (7). The effects of these two choices are now assessed. First, we consider vertical profiles of the turbulence coefficients, β R l i (z) ( Figure A5a). Figure 9a shows that the spatial structure of the mean concentration field at 0°is essentially unchanged from that for the standard configuration with β R l i,canopy and β i,above (Figure 5a). Since the differences between β R l i (z) and β R l i,canopy , β i,above are not very large (e.g., the maximum relative error is ∼50%), FM1 should not be especially sensitive to the inclusion of the vertical dependence of the turbulence coefficients. Indeed, the N MSE is 0.25 using β R l i,canopy , β i,above (Table 1) and 0.51 with β R l i (z) ( Table 2). It should be noted, however, that performance does not improve when the vertical dependence is specified. The reasons for this are unclear, but the change in the averaging region (from the canopy or overlying atmosphere to individual vertical levels) may play a role. The assumption that the turbulent fluctuations are linearly proportional to the mean velocity components (Equation (6)), which underlies the definition of the β, may break down on smaller vertical scales due to overfitting. Furthermore, the β(z) values are more subject to sampling effects. Second, the turbulent fluctuations are updated at the e-folding timescale, i.e., T update,i = T R l i . From the mean concentration field at 0° (Figure 9b), downwind dispersion is greatly reduced and the prediction is clearly worse. This is confirmed by the statistical metrics (Table 2): for example, N MSE increases from 0.25 (Table 1) to 1.96 when the turbulent fluctuations are updated at a time interval of T R l i rather than ∆t. Evidently choosing T update,i to be as short as possible improves the prediction. This is plausible inasmuch as the SDE, Equation (8), is equivalent to a diffusion equation only in the limit of white noise, i.e., a vanishing correlation timescale [22]. Physically, defining the turbulence timescale from the e-folding timescale, like the related integral timescale [37], amounts to specifying a length scale; however, the SDE, or more generally Brownian motion, has no spatial dependence or spatial scales. The choice T update,i = T R l i effectively imposes an artificially large length scale, which may be inconsistent with turbulent diffusion at the smallest possible scales. The lateral profiles (Figure 9c) confirm the preceding findings. Prescribing the vertical profile of the turbulence coefficients, β i = β R l i (z) yields nearly identical results to the standard configuration, T update,i = ∆t and β i = β R l i,canopy , β i,above . Updating the turbulent fluctuations less frequently, i.e., T update,i = T R l i , causes the profile to decay too rapidly with downwind distance. For brevity, the results for T update,i = T R l i and β i = β R l i (z) are not shown because they are very similar to those for the preceding case, T update,i = T R l i and β i = β R l i,canopy , β i,above . Aside from the vertical structure and correlation timescale, assumptions were also made about the horizontal structure and applicability to other wind directions. The dependence of β R l i,canopy on the wind direction is illustrated in Figure 10a. Values are not independent of the angle. To assess the validity of the choice of 45°for the standard configuration, the turbulence coefficients for 0° (Table A9) are applied to an FM1 calculation for 0°. The scalar field ( Figure A6a) shows enhanced cross-stream dispersion compared with the standard configuration; the lateral profile ( Figure A6c) shows enhanced contrast between the windward and leeward walls. A slight degradation of the performance is seen in the statistical performance measures (Table 2), e.g., the NMSE increases from 0.25 to 0.37. Using the turbulence coefficients for the same wind direction does not necessarily help.
The sensitivity to the spatial variation of the β R l i,canopy is assessed by considering the same turbulence coefficients for all repeating units β i,canopy . Compared with β R l i,canopy (0°), there are lower concentrations at the leeward walls ( Figure A6b) and the concentration decays more rapidly in the downwind direction ( Figure A6d). Neglecting the spatial variation of the turbulence coefficients has a negative impact: the N MSE increases from 0.25 to 0.74 (Table 2). Nevertheless, the sensitivity to the spatial variation of the turbulence coefficients is much smaller than that to a spatially varying mean flow.   Figure 10. Dependence of β R l i,canopy and κ * R l ij,canopy on wind angles. (a) β R l i,canopy ; (b) κ * R l ij,canopy .

FM2
The default FM2 configuration also relies on spatial averages over the canopy and the overlying region to determine the effective diffusivities at 45°. Different specifications of the effective diffusivities are now applied to simulations for 0°.
First, prescribing the vertical profiles, κ R l ij (z) (Figure A5b), has a small effect on the scalar field ( Figure A7a) and lateral profile ( Figure A7d). The statistical metrics indicate a small degradation when the vertical profiles are used (e.g., NMSE increasing from 0.40 to 0.46; Table 2), as with FM1. Second, the effective diffusivities are defined using another wind direction. The assumption that κ R l ij,canopy is independent of the angle is not accurate (Figure 10b). Hence, κ R l ij,canopy (0°) and κ ij,above (0°) (Table A9) are tested in place of the effective diffusivities at 45°. The scalar field ( Figure A7b) and lateral profile ( Figure A7e) show higher downwind concentrations and the errors decrease slightly (Table 2). Third, horizontally averaged diffusivities, κ ij,canopy , are used for all repeating units inside the canopy. Neglecting the spatial dependence lowers the concentrations along the leeward walls in the near-field region ( Figure A7c) and weakens concentration gradients ( Figure A7f). However, the statistical performance measures do not show much change (e.g., the N MSE decreases from 0.40 to 0.39).
The sensitivity to κ R l ij in FM2 is weaker than the sensitivity to β R l i in FM1. Not only does FM1 show great sensitivity to the correlation time, neglecting the dependence on the repeating units degrades performance. A possible explanation is that κ R l ij only affects diffusion, while β R l i affects advection and induces turbulent diffusion, which could lead to less robust behaviour.
The standard configuration of FM2 is limited to diagonal components of the effective diffusivity tensor, κ ij . Referring to the "microscopic" definition of the diffusivity, Equation (11), diagonal components, i = j, refer to motions arising from aligned fluctuations in the same direction. For a homogeneous medium and isotropic fluctuations, the off-diagonal components vanish, κ ij = 0, i = j = 0, as with the molecular viscosity of a Newtonian fluid, but this is not necessarily true for an inhomogeneous medium or one with boundaries. Hence it is possible that a better parameterisation could be obtained if the effective diffusivity included off-diagonal contributions.
Off-diagonal components can be calculated in the same way as the diagonal components (i.e., using Equations (10) and (11) with T ij = ∆T). Values for 0°and 45°are listed in Table A10. The sensitivity to the inclusion of off-diagonal components is illustrated in Figure 11: at 0°, the scalar field is asymmetric along the centreline, while at 45°, there are minimal differences with respect to the standard configuration. The asymmetric appearance for 0°may be related to the K xy and K yx components, which are always positive for repeating units inside the canopy: these components introduce correlations between x and y fluctuations in the scalar field, i.e., anisotropy. The inclusion of the off-diagonal components clearly degrades performance at 0°; the effects are much smaller at 45°. This is confirmed by the statistical metrics (e.g., the NMSE inside the canopy increases from 0.40 to 0.74 in the former case and from 0.60 to 0.68 in the latter case; Table A5). The negative impact of the off-diagonal components is perhaps surprising inasmuch as a parameterisation that takes fuller account of the turbulence statistics and does not require isotropy should perform better. A likely explanation is that the off-diagonal elements are assumed to be spatially uniform over a repeating unit. However, the non-zero off-diagonal components should be closely related to the local building geometry, which induces an inhomogeneous response. At 0°, u v differs in sign from one lateral boundary of a canyon unit to another (e.g., [38]); hence, K xy and K yx should be defined accordingly. Another contributing factor is sampling: at 45°, the values for the canyon and channel units are not identical, and the differences are comparable to or greater than those for the off-diagonal elements at 0°.

Discussion
The relatively good performance of the fast model can be attributed to the inclusion of spatially varying mean wind profiles. When it is assumed that the mean velocity is the same within all repeating units, the NMSE within the canopy is two to three orders of magnitude higher. Since air quality acceptance criteria (e.g., FAC2 0.5, FB 0.3, N MSE 0.3; [35]) are satisfied when the mean flow differs between repeating units, we conclude that the representation of the coarse-grained spatial variation of the mean flow is of primary importance. One may expect similar behaviour for domains that show comparable levels of spatial inhomogeneity. A simpler representation could be sufficient for an idealised street canyon; however, the importance of intersections to flow and dispersion within real cities is widely recognised [39].
The performance of the fast model is mostly insensitive to the choice of turbulence parameterisation. Temporally and spatially averaged fields from FM1 and FM2 show relatively small differences, as do the statistical metrics, when the correlation timescale is sufficiently short. Despite their different formulations, the models are largely equivalent: FM1 parameterises turbulent fluctuations directly; FM2 parameterises the turbulent diffusion arising from the turbulent fluctuations. Both approaches assume characteristic scales (e.g., the β i of FM1 and κ * ij of FM2 are normalised) and neglect off-diagonal diffusivity or turbulent velocity fluctuations associated with a mean velocity component in another direction. The general insensitivity to their precise specification, as demonstrated by the sensitivity tests of Section 6, suggests that the assumptions underlying them are robust; however, the approaches are formally equivalent only in the limit of a vanishing timestep. FM1 and FM2 could deviate significantly for large timesteps, but the need to maintain a small Courant number helps ensure that the differences between them are small in practice. A small timestep mitigates differences associated with FM1 introducing fluctuations in the velocity field (which could introduce a drift over short time interval).
The fast model has two inherent limitations. First, it requires calibration to determine the geometric constants for the estimation of the mean velocity profiles and the coefficients for the turbulence parameterisations. This necessitates additional calculations, preferably from reference CFD simulations for a specific wind direction, which could limit the model's applicability to other wind directions. The sensitivity tests for 0°and 45°, however, suggest that the turbulence coefficients and effective diffusivity for one wind direction can be applied to another. In any case, some form of calibration is almost unavoidable for simplified models that do not solve the governing equations: adjustable parameters are present in the GPM, as well as in fast dispersion models SIRANE [9] and QUIC [6]. Second, the spatial dependence of the mean velocity field inside each repeating unit is neglected while the turbulent fluctuations are spatially uncorrelated. The performance of the fast model suggests, however, that, at least for the cubic building array, fine-grained spatial structure is not essential. Averaging over repeating units should be an acceptable approximation so long as the spatial variability of the flow within them is not too great.
Although this study is limited to idealised configurations, its key components should be robust. The vortex method for the mean flow has already been extended to realistic building geometries [15]: the flow inside deep canopies is typically controlled by a few vortex sheets, which also facilitates calibration. The turbulence parameterisations have yet to be tested for more realistic domains, but they should perform comparably if the relationship between turbulent fluctuations and the mean flow does not differ significantly (cf. Equation (5)).
The fast model is unlikely to work well for a more realistic configuration if the nature of the turbulence should change from that for the cubic building array. One such scenario is unstable stratification (e.g., ground heating), for which convective plumes induce qualitative changes in the mean circulation and turbulence statistics [40]. Another is non-cuboidal buildings or ones with complex facade features for which the local turbulence will not be determined by the spatial average over a repeating unit. Nevertheless, the fast model should apply to realistic urban areas if the buildings are approximately cuboidal and the flow is neutral or stable.

Conclusions
The main finding of this study is that the dispersion of a passive scalar inside a cubic building array can be effectively modelled by combining estimated mean velocity profiles from a vortex method with a parameterization for turbulent diffusion based on introducing velocity fluctuations or specifying an effective diffusivity. Compared to reference LES simulations at 0°and 45°, the spatial structure of the steady (time-averaged) scalar field is largely captured while air quality acceptance criteria are satisfied. Although there is a relatively modest sacrifice in accuracy, computational efficiency is greatly increased. Table 3 compares the computational performance of both versions of the fast model with LES. FM1 and FM2 are ∼40-50 times faster than LES; the speed-up would be even greater if the LES spin-up were included in the comparison. FM2 is more computationally efficient because the numerical solution does not require the introduction of random fluctuations everywhere in the domain; moreover, sampling issues potentially arising from a finite timestep are avoided. Given its greater computational efficiency and simpler implementation, as well as the close agreement with FM1, FM2 is the better choice for most applications. Table 3. Comparison of computational time required to simulate scalar dispersion at 0°for 5000 s. The total CPU time consumed is given by product of the number of computational cores and the time elapses. The computational speed-up is calculated relative to the total CPU time for LES.

LES FM1 FM2
Number The most important implication is that urban pollutants can be modelled using a representation of the coarse-grained mean flow and turbulence. On one hand, this represents a significant simplification compared to direct CFD; on the other hand, the flow is represented far more accurately than in the Gaussian plume model or network models such as SIRANE [9] or MUNICH [10]. Any simplified or intermediate dispersion model necessarily involves simplification. In the fast model, dispersion inside the canopy is driven by the mean flow and turbulence varying on the scale of the building geometry. This is a multiscale modelling strategy which assumes that the flow is approximately constant over a finite region and that turbulent diffusion can be parameterised in terms of the mean velocity components: the use of the vortex method for the mean flow or the FM1 and FM2 turbulence parameterisations are not strictly required. Moreover, its applicability should extend beyond the cubic building array. While real urban canopies do not have repeating units, subdomains over which the mean flow and turbulence are approximately homogeneous could be defined empirically, though at the cost of added technical complexity. Nevertheless, the fast model should still work well for approximately grid-like neighbourhoods.
A secondary implication is that the fast model can be applied to time-critical applications such as operational air quality prediction or emergency dispersion modelling. Although calibration using reference CFD simulations is time-consuming, this could be performed in advance for problems in which the region of interest (e.g., a densely populated neighbourhood) is known, allowing for the FM's computational efficiency to be realised in practice. If the building geometry were approximately uniform, as with the cubic building array, calibration could be limited to a small number of directions (or perhaps just a single one) as the turbulence coefficients and effective diffusivities should show minimal sensitivity to the wind direction.
There are several directions for future work. Most obviously, the method should be tested for more realistic geometries and flow configurations. As explained in the Discussion, there are good reasons for expecting the fast model's applicability to extend beyond idealised flow over a cubic building array. Nevertheless, the extent of applicability needs to be determined. Furthermore, the turbulence parameterisations could be improved. A clear weakness of the fast model is that cross-stream dispersion is underpredicted for FM1 and FM2. This shortcoming could be partially ameliorated in the context of FM2 by allowing for correlations between turbulent fluctuations in different directions, which would amount to a proper specification of the off-diagonal components. For FM1 a more realistic stochastic parameterisation could be adopted [22]. Finally, more computationally efficient methods for solving the advection-diffusion equation could be explored. A Lagrangian stochastic method would be a natural choice for FM1.

Conflicts of Interest:
The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Gaussian Plume Model (GPM)
The GPM yields a solution for the advection-diffusion equation given a pollutant emission rate (for a point source), atmospheric stability, wind speed and direction. The main assumptions are as follows: • The velocity is a constant, u = (u, 0, 0). • The eddy diffusivities, K i (x), are isotropic and functions of the downwind distance x only. • Diffusion is weak compared to advection in the streamwise (x) direction. • The pollutant does not penetrate the ground. • There are Dirichlet boundary conditions at the lateral edges, C(0, y, z) = 0, There is a Neumann boundary condition ∂C ∂z (x, y, 0) = 0 at the ground. • The source is located at a height H ≥ 0. The steady solution for the concentration may be written as The standard deviations, σ 2 y and σ 2 z , depend on the atmospheric stability. Assuming neutral stability, σ y = ax b , σ z = cx d + f , where a = 156, b = 0.89, c = 108.2, d = 1.098 and f = 2.0 [41].

Appendix B.1. Velocity Statistics
The mean streamwise velocity predicted by the LES model described in Section 3.2 is validated against wind-tunnel measurements of flow over a regular array of cubic buildings [42]. Vertical profiles of the mean streamwise velocity measured at the center of canyon units are compared in Figure A1. There is good agreement above the ground, i.e., z/H > 0.1.
Statistical performance measures (Appendix C) confirm that the LES model simulates the mean wind profile accurately (Table A1). The NMSE and FB are close to 0 and the correlation coefficient is close to 1. Table A1. Statistical performance measures corresponding to Figure A1. the temporally averaged streamwise velocity. The mean streamwise velocity and rms streamwise velocity fluctuation at the centre of a canyon unit are compared to the wind-tunnel data of [42].

Appendix B.2. Concentration Statistics
Concentration statistics generated by a point source located at the centre of a unitaspect-ratio street canyon are validated using experimental data [43]. The lateral profiles show good agreement ( Figure A2). This is confirmed by the statistical metrics (Table A2). Hence, we conclude that the LES model is also capable of simulating the dispersion from a point source.  Table A2. Statistical performance measures corresponding to Figure A2. The scalar statistics comprise spanwise profiles at x 1 = x s + 2H and x 2 = x s + 6H, where x s is the source location. In order to validate the simulation data (D s ) against measurements (D e ), we calculate the following statistical performance measures [35]: fractional bias (FB),