Qualitative and Quantitative Investigation of Multiple Large Eddy Simulation Aspects for Pollutant Dispersion in Street Canyons Using OpenFOAM

Air pollution is probably the single largest environment risk to health and urban streets are the localized, relevant hotspots. Numerous studies reviewed the state-of-the-art models, proposed best-practice guidelines and explored, using various software, how different approaches (e.g., Reynolds-averaged Navier–Stokes (RANS), large eddy simulations (LES)) inter-compare. Open source tools are continuously attracting interest but lack of similar, extensive and comprehensive investigations. At the same time, their configuration varies significantly among the related studies leading to non-reproducible results. Therefore, the typical quasi-2D street canyon geometry was selected to employ the well-known open-source software OpenFOAM and to investigate and validate the main parameters affecting LES transient simulation of a pollutant dispersion. In brief, domain height slightly affected street level concentration but source height had a major impact. All sub-grid scale models predicted the velocity profiles adequately, but the k-equation SGS model best-resolved pollutant dispersion. Finally, an easily reproducible LES configuration is proposed that provided a satisfactory compromise between computational demands and accuracy.


Introduction
Air pollution is a well-known problem of modern cities and despite the great on pollutant control technology, it is still very difficult to ensure the elimination of the anticipated consequences in environment and human population.Therefore, the need for studying air pollution becomes apparent, since it is recognized as one of the main reasons of a large number of deaths [1].Air pollution can be studied systematically by adopting various experimental and numerical methods, such as small air quality models (AQMs) [2], wind tunnel experiments [3], field measurements [4], computational fluid dynamics (CFD) models [5][6][7], and semi-empirical models [8,9].
Turbulence modelling is a core issue to address in CFD and it poses unique challenges for flows in the urban infrastructure [10].The large size of buildings and high velocities of the atmospheric boundary layer result in highly turbulent flows [11], which are computationally expensive to simulate [12].Depending on the selected CFD technique, the whole spectrum of turbulent motions can either be resolved, modeled, or a combination of both.
Solutions of the Navier-Stokes equations for turbulent flows can be obtained using various numerical and analytical approaches, providing with different level of accuracy in each case.The most 'accurate' approach is the direct numerical simulation (DNS), which involves the numerical solution of the full spectrum of turbulence without the need of any modeling [13].Even though DNS provides a high-fidelity representation of the physical phenomena, it is unfeasible for urban-scale flows and it is usually executed for flows with low Reynolds numbers owe to its high-computational cost [14].On the other end, the Reynolds-averaged Navier-Stokes (RANS) simulation is a more feasible method that provides a time-averaged representation (or steady-state solution) of the studied phenomena with less demands in terms of computer speed and memory [15].RANS (and unsteady-RANS) approach employs the Reynolds decomposition, a mathematical technique that separates the flow properties into the average and fluctuating components.This decomposition introduces the Reynolds-stress tensor, creating the well-known 'closure problem', which is solved by the addition of a turbulence model.
Another modelling approach, which demands less computational resources than DNS and resolves turbulence with more accuracy than unsteady-RANS, is the large eddy simulation (LES) method.LES begins with the filtering of Navier-Stokes equations, which divides turbulence into the large-scale resolved part and the sub-grid scale (SGS) modelled part [13].The LES approach provides a time-dependent solution of the problem, with more details than RANS, as LES resolves a larger percentage of turbulence spectrum [16].RANS and LES are the most common techniques to simulate atmospheric flows between buildings [17,18].The general notion is that LES can capture the simulated phenomena with greater details, but it requires increasing computational power and time [19].Due to their transient solution, LES can provide better prediction for the maximum values of concentrations, which is important for simulations involving toxic materials [20].Moreover, LES is proven to provide results that are validated with experimental data from water channels [21,22], as well as wind tunnels [17,[23][24][25].Regarding the quasi-2D LES setup, there are many applied LES configurations.For example, Chung and Liu [26] employed the classic Smagorinsky SGS model [27], with the selected value of the Smagorinsky constant (C s ) being taken equal to 0.135.Moreover, Li studied the effect of the aspect ratio (AR) and heated walls to the flow and dispersion, for AR from 2 to 10 [21], using the one-equation SGS model by Moeng and Sullivan [28], which solves an additional transport equation for the SGS turbulence kinetic energy.At the solid boundaries, a 1/7 power-law wall model [29] was applied.The flow field was validated against water channel experimental data [30], while the concentration field was compared with wind tunnel measurements [31][32][33], presenting good agreement in all cases.Furthermore, Cai performed LES with the RAMS code [34] to investigate street canyon flow for quasi-2D [35] and 3D geometries [36].The standard Smagorinsky model was adopted for the SGS modelling, with C s equal to 0.1.The value of the constant was further tested by Cui et al. [37] to show that larger C s values can cause excessive dissipation, while smaller ones create numerical errors.The obtained numerical results were compared with available wind tunnel experimental data [3,38,39] for the quasi-2D geometry [40] and water channel experimental data [41,42] for the 3D geometry.
Additionally, LES has been coupled with simplified and computationally affordable custom models of photochemical reactions under neutral atmospheric conditions in order to validate simple chemical AQMs [43].For a deep street canyon (AR = 2), Zhong et al. [22] implemented the typical boundary conditions and compared the results with water channel data by Li et al. [30].The k-equation eddy-viscosity SGS model [44] along with the logarithmic law of the rough wall [45] were used for the LES method.The obtained datasets of the spatial and temporal distributions of the reacting pollutants were used for validating a simple chemical box model [46].
Recently, Kikumoto and Ooka [47] performed a high resolution LES study for the investigation of turbulent dispersion in a unity street canyon, with the aim to reproduce a physical experiment of passive tracer dispersion.They predicted the mean and root-mean-square (RMS) velocity and concentration levels, exhibiting fair agreement with the experimental data, except from slight underestimation for the concentration fluctuations near the source.The same work suggests that the LES method could eventually replace wind tunnel experiments, although with the burden of extremely high grid resolution.
In all these studies, the main objective was to simulate a street canyon flow using LES, while only in very few of them was there an attempt to describe and fully justify the appropriate LES configuration (e.g., SGS models, grid resolution, time-step and boundary conditions, among others) for the considered problem.To spend more time studying the employed phenomena, available CFD commercial codes are treated as black boxes and past works or guidelines are followed without adequate explanation or justification.
Therefore, in this study we selected a simple, yet important, street canyon geometry (described in Section 2) to conduct a comprehensive investigation, with verification and validation, of various parameters affecting LES.The open-source code OpenFOAM v5.0 [48] has been utilized for this investigation (details discussed in Section 3), for a typical quasi-2D geometry in 'real' scale, i.e., the height of buildings is set to 10 m.The investigated parameters are the domain height, the grid resolution, time step, Reynolds number, source height, four of the available SGS models, and three different wall functions.The obtained results for each parameter are qualitatively and quantitatively analyzed to understand its impact and identify the most suitable modelling configuration (presented in Section 4).Results are also compared extensively with available experimental and numerical data from the literature.Finally, the most important findings are summarized along with potential future studies (listed in Section 5).

Physical Problem
Many cities have complex layouts, created by a diverse variety of structures, ranging from the usual cubical buildings to skyscrapers and stadiums.Depending on the computational power and the scope of the work, the computational grid may contain whole cities, building blocks or street canyon(s), which is the basic structural unit of the city.Figure 1a presents the special (ideal case) of an 'infinite length' (L) street canyon, formed by two (rows of) rectangular parallelepiped buildings with height (H) and width (B), and a perpendicular wind direction (left to right).Because of the infinite length and the resulted symmetry, this is considered as a quasi-2D geometry.
Both the height and width of the buildings are set to 10 m, while the street width (W) can vary to depending on the studied AR (here in AR = 1).In the streamwise direction, the computational domain begins from the second half of the upwind and ends at the first half of the downwind building.In the spanwise direction, half (H/2) of the buildings' width is contained, effectively simulating a 'slice' of the street canyon.The upper boundary extends at five building heights (5H); more details in a following section.
Since the studied area is the smallest possible and the geometry is simple and ideal, the creation of a high resolution, fully structured, hexahedral mesh is easy and applicable.Structured hexahedral meshes discretize the computational domain more efficiently [49] and can have cells aligned with the flow direction, thus achieving smaller numerical errors [50] and better spatial convergence [51].The height and width of the buildings are discretized with 120 and 60 cells, respectively, and the depth of the street is set to 60 cells, satisfying the demand for 3D LES.The cells are smaller near the solid boundaries, to resolve the flow characteristics of the wall boundary layer, reaching a 0.0038 m width.The Reynolds number (Re = U ref •H/v) is around 3,200,000, calculated using the reference velocity U ref at the top (5 m•s −1 ), the kinematic viscosity of air, v, and the building height H.

Governing Equations
The current study is carried out by using the LES methodology [52].The unsteady, threedimensional incompressible filtered Navier-Stokes equations [16] can be written in Cartesian tensor notation as follows.
Filtered continuity equation where i u is the filtered mean velocity component in the i-direction, j u is the filtered mean velocity component in the j-direction, xi,j is the distance in the i-and j-directions, p is the filtered pressure, v is the kinematic molecular viscosity, ρ is the density, t is the time and the subscript indices i, j (= 1, 2, 3) denote the three space coordinates.τ ij is the SGS stress tensor and is analysed by the Leonard decomposition [53] as

Governing Equations
The current study is carried out by using the LES methodology [52].The unsteady, threedimensional incompressible filtered Navier-Stokes equations [16] can be written in Cartesian tensor notation as follows.
Filtered continuity equation Filtered momentum equations where u i is the filtered mean velocity component in the i-direction, u j is the filtered mean velocity component in the j-direction, x i,j is the distance in the iand jdirections, p is the filtered pressure, v is the kinematic molecular viscosity, ρ is the density, t is the time and the subscript indices i, j (= 1, 2, 3) denote the three space coordinates.τ ij is the SGS stress tensor and is analysed by the Leonard decomposition [53] as where L ij is the Leonard term, C ij is the Cross term and R ij is the Reynolds term.Each term represents physical interactions between the resolved and modelled scales of the flow that arise from the LES approach [13].
The coupled pressure and velocity are solved with the PISO (Pressure-Implicit with Splitting of Operators) algorithm by Issa [54].The pressure equation is treated implicitly by the PCG solver using the DIC preconditioner, while the momentum equation is resolved by a simple solver, namely smoothsolver in OpenFOAM, along with the GaussSeidel smoother.The terms of the equations are discretized with second-order numerical approximations, according to the standard practice [55].The time first derivative is discretized with the backward difference scheme (backward in OpenFOAM), while the UMIST scheme [56] is employed for the convective terms.The gradient and the Laplacian terms are discretized using the central difference schemes, Gauss linear and Gauss linear corrected [48], respectively.
The sub-grid part of turbulence is represented by the SGS stress tensor (τ ij ) and can be modelled with various methods.A common approach is the Boussinesq approximation, which assumes that the deviatoric part of the SGS tensor is proportional to the resolved strain rate tensor S ij , through the eddy viscosity, v sgs , where τ d ij is the deviatoric part, τ ij is the normal component and τ kk is the shear component of the SGS stress tensor.δ ij is the Kronecker delta and 1  3 τ kk δ ij is the hydrostatic tensor, which is embedded in the modified filtered pressure term.The resolved strain rate tensor S ij is written as Doing so, the problem of unresolved modelled turbulence has been transformed to the calculation of eddy viscosity v sgs .In the present study, the considered SGS models are presented in Section 3.2.
Finally, the selected incompressible solver pisoFoam [48] was modified to solve the convectiondiffusion equation for a passive scalar, at every time-step.
Convection-diffusion equation where C is the concentration, P are the characteristics of the pollutant source (i.e., geometry and emission rate) and D total is the total diffusion coefficient where D is the molecular diffusion coefficient, v sgs is the eddy viscosity, and Sc t is the turbulent Schmidt number.The value of Sc t is taken equal to 0.7, in accordance with previous atmospheric dispersion studies [22,57], while D was equal to 1.64•10 −5 m 2 •s −1 , which is the molecular diffusion coefficient for, a typical pollutant, CO 2 at 25 • C [58].

Standard Smagorinsky Model
The standard Smagorinsky model [27] is the most well-known SGS model and constitutes the base for the development of many other SGS models.The main equation of the model based on the eddy-viscosity hypothesis can be written as where C s is the Smagorinsky coefficient, |S| = 2S ij S ij is the magnitude of the resolved strain rate tensor and ∆ is the filter width, usually associated with the cell dimensions The behavior of the standard Smagorinsky depends heavily on its constant, C s , and suffers from many weaknesses (e.g., backscatter phenomenon, tuning of C s and requirement of a damping function).The model also fails to eliminate the eddy viscosity near to the wall boundaries and the requirement of the van Driest damping function [59] is necessary to reduce the C s value to zero and ensure the no-slip boundary condition.

WALE Model
The variations of the Smagorinsky model calculate the SGS turbulence, based on its connection with the symmetric part of the velocity gradient tensor, i.e., the resolved strain rate tensor, effectively capturing the vortex deformation.The local strain and rotation rates of the smallest resolved turbulence scales are taking into account in the WALE model by Ducros et al. [60].It is reported that WALE reproduces the transition from laminar to turbulent flow and the correct scaling of turbulence close to the walls without the requirement of a damping function [61].According to the WALE model the eddy viscosity is given by where ∆ is the usual grid length and C w = 0.325 is the WALE constant.S d ij is the deviatoric symmetric part of the square for the resolved velocity gradient tensor and is defined as where g ij = ∂u i /∂x j are the velocity gradients.

k-Equation Model
The above-mentioned algebraic models connect v sgs to the local instantaneous resolved velocity gradients.On the other hand, there are regions with low resolved gradients, but high SGS kinetic energy because of convection from regions with high sub-grid scale turbulence [62].A common approach to capture the historic effects of production, dissipation, and diffusion of the SGS kinetic energy (k sgs ) is the resolving of its transport equation.The k-equation model by Yoshizawa and Horiuti [63] resolves the transport equation of and the eddy viscosity is given by where v is the laminar viscosity, C k (=0.094) and C e (=1.048) are the model constants.The terms in the right hand of Equation (11) represent the diffusion, production and dissipation of k sgs , respectively.The constants C k and C e are calculated with the assumption of balance between SGS energy production and dissipation rates [64], similar to the standard Smagorinsky model, presenting the same limitations.Consequently, improved versions of the k-equation model have been proposed to calculate the constants dynamically, depending on the local turbulence.

Dynamic k-Equation Model
The localised dynamic k-equation model by Kim and Menon [65] employs a similar operation to that introduced by Lilly [66] for the dynamic Smagorinsky model.A second test filter ∆ = 2∆ and the Germano identity are used to derive equations that relate the values of C k and C e to quantities resolved at the test filter scale [67] C k, dynamic = 1 2 where and

Boundary Conditions
In the streamwise direction, periodic boundary conditions were imposed (cyclic in OpenFOAM) for all the flow variables, i.e., U, p, v sgs , and C. The main flow was recycled from the outlet boundary patch to the inlet, ensuring that the boundary conditions are time-dependent, a requisite for LES [68].To initialize the turbulence, a single velocity vector with the value of the reference velocity (5 m•s −1 ) was imposed to the whole computational grid.The simulation of the wind flow begun and gradually the turbulence inside the street canyon dissipated forming the expected recirculating vortex, while the turbulence in the upper atmosphere created a realistic wind profile.The simulation ran for a sufficient time span, in which the obtained average velocity was statistically correct, when compared with experimental data and the flow was judged as realistic (i.e., the instantaneous snapshots of velocity showed no indications of laminar flow).Only after that, the emission of the pollutants from the source begun.This method is well-established in similar works, e.g., [21,22,37].The recycling of the concentration was countered by overwriting C at the inlet patch, at every time-step with the utility scalarFixedValueConstraint, available after OpenFOAM 3.0.
Another option for the streamwise boundaries could be the mapped boundary conditions, in which the user can choose which parameters will be recycled.These conditions were used extensively by the authors in OpenFOAM v2.3.1, with the results of flow being statistically valid for AR = 1 and 0.5.The method was also used successfully for AR = 0.33, but for AR = 0.25 a recirculating flow was found to develop just after the inlet patch above the upwind building.As solution progressed, this flow was amplified, eventually causing the simulation to diverge.The reason for this was not investigated further neither this method was tested in newer versions.Therefore, cyclic conditions have been employed which proved more stable in either versions of OpenFOAM, v2.3.1 and v5.0.At the top of the domain, the symmetry boundary condition (symmetryPlane in OpenFOAM) was set, following the common practice for the quasi-2D setup, e.g., [22,40,69].This means that the normal velocity on this patch is assigned to zero and for the scalars there is a zero-normal gradient.
Periodic conditions are used also for the spanwise direction, again as the common practice for the quasi-2D configuration.The recycling of the random small velocity vectors from patch to patch is not affecting the main flow.Furthermore, tests with symmetry for these patches resulted in unnatural flows.The wall boundaries were set to no-slip condition for the velocity field (fixedValue in OpenFOAM) where U is the filtered velocity component and i, j and k are the three x, y and z-coordinates, respectively.The fixed gradient (zeroGradient in OpenFOAM) was used for the scalar parameters.
where ϕ is the scalar quantity, i.e., p, v sgs , and C.
For high Re flows or relatively coarse grids, depending on the non-dimensional distance from the wall to the first grid point (y + ), special treatment should be employed for the flow characteristics on the wall boundaries [13].OpenFOAM provides a variety of wall functions to handle flow characteristics on walls [70].Herein, two functions studied combined with the standard Smagorinsky model to calculate the value of turbulent viscosity v sgs on the walls.The equation that relates v sgs to y + is where κ = 0.41 is the von Karman's constant and E = 9.8 is a constant that describes the velocity profile near the wall.The well-known wall function by Launder and Spalding [71] is implemented in OpenFOAM v5.0 as nutUSpaldingWallFunction, for which y + is calculated as a function of the dimensionless velocity u + A simple logarithmic wall function (nutkWallFunction) was also tested, that estimates v sgs using Equation (20) for values of y + > 11, or else assumes v sgs = 0.

Computational Details
The solution of LES provides dynamic time-series, which are usually averaged over a feasible and representative time-span.Previous works have specified this time-span, as a function of U ref and a characteristic length.For example, Li et al. [69] used 400 H/U ref to reach a "pseudo-steady state" and another 300 H/U ref to collect the appropriate statistics.Employing a similar logic for turbulence generation, Boppana et al. [72] used the streamwise length of the domain L x to determine the respective time-spans as 120 and 200 L x /U ref .For this work, the fully developed flow field was obtained after 150 H/U ref , while the statistics were collected over 650 H/U ref (or 325 L x /U ref , or 1300 s).The computational cost for these calculations requires the employment of high-performance computing (HPC).Indicatively, the fine grid required a computational time of three months using 240 CPU cores, which for the medium grid dropped to seventeen days, this time employing 192 CPU cores.

Results and Discussion
This section presents the detailed results of the investigation along with the qualitative analysis and quantitative comparison with literature data.The results relay to average and instantaneous values of the velocity and to averaged concentrations, for the wind flow and pollutant dispersion, respectively.At first, the investigation is based on a standard model configuration employing the Standard Smagorinsky model and the van Driest wall function, over a coarse resolution grid (80-cells across the street canyon; 768,000 in total) and a time-step of 1 × 10 −3 s.As the investigation progresses, the model configuration is updated.
Both experimental and numerical literature data are used to compare with own results.For the velocity components the experimental data by Li et al. [30] and the numerical data (also LES) by Li [73] are employed.For the pollutant concentrations the experimental data by Pavageau and Schatzmann [31], Meroney et al. [32], Pavageau [33], and Hoydysh and Dabberdt [74] are employed.The quantitative comparison follows well-known metrics (i.e., FAC2, FB, and NMSE) proposed by collaborative projects such as the COST Action 732 [75], according to the works of Hanna [76] and Chang and Hanna [77], and recently for the validation of LES, e.g., [25,78,79].

Domain Height
Best practice guidelines suggest that the external boundary patches should be far away from the investigated area, typically around five times the average building height [55].In similar studies, this area has been taken equal to 1H [80], 2H [22], 3H [69], or 4H [81] building heights.This dimension controls the available space for the wind flow (velocity profile) to be developed.It is important to be large enough, so the profile can resemble the 'stratification' of the urban boundary layer.The first option for a larger upper atmosphere height is to increase the number of cells proportionally, which consequently increases the computational cost for the less studied part of the flow.The second option is to keep the same number of cells and adjust the cell-to-cell ratio, so that the cells at the roof level inside the street canyon and the upper atmosphere match the same size.Here, the second option was selected.Also, note that the inlet, outlet, front, and back patches did not follow this recommendation owe to the cyclic boundary conditions.
Therefore, a sensitivity analysis was conducted using 2H, 3H, 4H, and 5H to determine the appropriate height.The grid resolution was kept at 80 cells per building height to compromise between accuracy and computational cost for the 2H-5H, however, a medium (120 cells) and a fine (180 cells) grid were also examined for 5H.More details on the grid resolution are presented in the next section.Figure 2 presents the average concentrations (normalised by C max ) on the two walls of the unity street canyon compared with experimental data by Pavageau and Schatzmann [31], Meroney et al. [32], Pavageau [33], and Hoydysh and Dabberdt [74].The LES results by Li [73] are also exhibited for comparison.The concentration results for the upwind wall are in fair agreement with the experimental data and the obtained values for the four upper atmosphere heights are nearly identical, while the largest variation is observed at the upper part of the canyon.For the downwind wall, the results for 2H, 3H, and 4H are different among them, while for 5H are marginally closer to the experimental data.The The concentration results for the upwind wall are in fair agreement with the experimental data and the obtained values for the four upper atmosphere heights are nearly identical, while the largest variation is observed at the upper part of the canyon.For the downwind wall, the results for 2H, 3H, and 4H are different among them, while for 5H are marginally closer to the experimental data.The results for the downwind wall are affected more than those for the upwind side, as this side is characterized by the entrainment of clean air from the upper atmosphere, which is in turn controlled by the wind velocity profile.For example, the transition of the flow from the top of the domain to the shear layer for the 2H is more abrupt and unnatural than that of 5H.On the other hand, the upwind side is influenced by the phenomena developed inside the canyon, mainly the recirculating vortex that dominates the flow.It is concluded that the experimental data are replicated better in the 5H case due to the larger available area which allows for a less-constrained development of the wind flow.

Grid Resolution
A grid sensitivity analysis was conducted to quantify the error induced by discretization, utilizing the grid convergence index (GCI) [82] and the LES index of quality (or LES_IQ) [83].
The GCI estimates an error in a 'loose statistical sense' for measuring the effect of a change on the grid resolution, i.e., the refinement/transition from one grid to another and can be calculated as where F s is the safety factor, r is the refinement ratio, and p is the order accuracy of the used method, while f 1 and f 2 are the solutions from each used grid.Three grids with increasing density, i.e., a constant refinement ratio of 1.5 with F s = 1.25, were tested following the common practice [12,20,84,85]: a coarse grid of 768,000 cells (∆x = ∆z = 0.125 m, ∆y = 0.063 m), a medium grid of 2,592,000 cells (∆x = ∆z = 0.083 m, ∆y = 0.042 m) and a fine grid of 8,748,000 cells (∆x = ∆z = 0.056 m, ∆y = 0.028 m).Note that the adopted numerical schemes for the discretization are second order accurate (p = 2).The grid-induced error was inspected at three areas, the relatively low-resolution area in the middle of the canyon and the denser upwind and downwind areas, where the polluted air removal and clean air entrainment occur, respectively.Thus, the GCI was calculated in the middle vertical plane of the grid (y/L = 0.5) for three profiles across the canyon-width at x/W = 0.25, 0.50, and 0.75, respectively.The GCI fine values for the medium to fine grid refinement are presented in Figure 3.
The GCI fine values are proportionally larger for U z /U ref than U x /U ref , although the absolute values are smaller for the first.The larger differences for U z /U ref are observed in the middle height of the canyon, where the grid is relatively coarse, as shown in Figure 3d,f.On the other hand, for U x /U ref the lower GCI fine values are calculated inside the canyon area, while abnormally high values are identified in the area above the buildings (Figure 3a-c), although the velocity at the top of the domain is virtually the same (~5.1 m/s) for all the grids.This could be linked with the observation that the velocity profile for the fine grid adapts more abruptly towards the street canyon.
The GCI medium (not presented here; for the coarse to medium grid refinement) was found to be in all cases larger than that of the fine grid.For example, the average GCI medium for the U z /U ref inside the canyon was calculated as 1.8, 0.5, and 1.1% for the three vertical profiles, while the values for GCI fine were found to be 0.6, 0.3, and 1%, respectively.This trend is also valid for the U x /U ref , with the values for GCI medium were found to be equal to 1.1, 1.9, and 1.3% for the three vertical lines, while the values for GCI fine were recorded as 0.7, 1, and 0.7%, respectively.As expected, the transition from the medium to fine grid gives smaller differences than the transition from the coarse to the medium.Nevertheless, the GCI fine is less than 1%, implying that the medium grid is sufficient for the wind flow solution.
The grid-induced error was inspected at three areas, the relatively low-resolution area in the middle of the canyon and the denser upwind and downwind areas, where the polluted air removal and clean air entrainment occur, respectively.Thus, the GCI was calculated in the middle vertical plane of the grid (y/L = 0.5) for three profiles across the canyon-width at x/W = 0.25, 0.50, and 0.75, respectively.The GCIfine values for the medium to fine grid refinement are presented in Figure 3.The GCIfine values are proportionally larger for Uz/Uref than Ux/Uref, although the absolute values are smaller for the first.The larger differences for Uz/Uref are observed in the middle height of the canyon, where the grid is relatively coarse, as shown in Figure 3d,f.On the other hand, for Ux/Uref the lower GCIfine values are calculated inside the canyon area, while abnormally high values are identified in the area above the buildings (Figure 3a-c), although the velocity at the top of the domain is virtually the same (~5.1 m/s) for all the grids.This could be linked with the observation that the velocity profile for the fine grid adapts more abruptly towards the street canyon.
The GCImedium (not presented here; for the coarse to medium grid refinement) was found to be in all cases larger than that of the fine grid.For example, the average GCImedium for the Uz/Uref inside the canyon was calculated as 1.8, 0.5, and 1.1% for the three vertical profiles, while the values for GCIfine were found to be 0.6, 0.3, and 1%, respectively.This trend is also valid for the Ux/Uref, with the values for GCImedium were found to be equal to 1.1, 1.9, and 1.3% for the three vertical lines, while the values The average concentration profiles of the numerical and experimental data for the two canyon walls are compared in Figure 2. As the grid resolution increases, the concentration of the downwind wall approaches the experimental and numerical values.The best results obtained from the fine grid, although the results from the medium grid are also in good agreement with the experiments.For the upwind wall and from the ground up to z/H = 0.6, the three solutions are identical and larger than the experimental results, while they differ for the distance until the roof height.For the part from z/H = 0.7 up to the roof, the results of the coarse grid are closer to the experimental, while the medium grid exhibits reasonable agreement.The results for GCI for C/C max are presented in Figure 4.
wall approaches the experimental and numerical values.The best results obtained from the fine grid, although the results from the medium grid are also in good agreement with the experiments.For the upwind wall and from the ground up to z/H = 0.6, the three solutions are identical and larger than the experimental results, while they differ for the distance until the roof height.For the part from z/H = 0.7 up to the roof, the results of the coarse grid are closer to the experimental, while the medium grid exhibits reasonable agreement.The results for GCI for C/Cmax are presented in Figure 4.For C/Cmax, the average GCImedium was calculated as 0.8% and 3.1% for the upwind and the downwind walls, respectively.The GCIfine was reduced even more for the downwind wall at 2.2%, but was increased for the upwind wall at 1.2%.Nevertheless, this discrepancy reveals the drawbacks of GCI, as the transition from medium to fine just moved the results closer to the experimental, as shown in Figure 2.
While, the use of GCI determined that the uncertainty induced by grid is low enough, other metrics are required to find out the percentage of the resolved turbulence.According to Pope [86], LES is properly employed, when 80% of the turbulence kinetic energy is resolved, a value that will be examined using the LES index of quality (or LES_IQ) by Celik et.al. [83], with the equation involving the laminar ν and turbulent viscosity νsgs For C/C max , the average GCI medium was calculated as 0.8% and 3.1% for the upwind and the downwind walls, respectively.The GCI fine was reduced even more for the downwind wall at 2.2%, but was increased for the upwind wall at 1.2%.Nevertheless, this discrepancy reveals the drawbacks of GCI, as the transition from medium to fine just moved the results closer to the experimental, as shown in Figure 2.
While, the use of GCI determined that the uncertainty induced by grid is low enough, other metrics are required to find out the percentage of the resolved turbulence.According to Pope [86], LES is properly employed, when 80% of the turbulence kinetic energy is resolved, a value that will be examined using the LES index of quality (or LES_IQ) by Celik et al. [83], with the equation involving the laminar ν and turbulent viscosity ν sgs In Figure 5, it is observed that a region of relatively lower LES_IQ v exists in front of the downwind building for the selected three grids.This is the area where occurs the entrainment of clean air and is characterized by intensive turbulence.Therefore, very coarse grids may have difficulties to resolve the turbulence kinetic energy at the acceptable 80% for this area, which does not occur for the selected grids.Another area of intense turbulence is the shear layer between the roofs, where LES_IQ v exceeds 85% for all resolutions, due to the dense grid.Outside the street canyon, LES_IQ v has significantly lower values than inside the canyon, reaching 70% for the coarse grid in Figure 5a.Similarly, in the same area, the medium and the fine grid in Figure 5b,c present their lower values at 75% and 80%, respectively.On the other hand, this should not affect the quality when using the medium grid, because the main area of interest is inside street canyon.In Figure 5, it is observed that a region of relatively lower LES_IQv exists in front of the downwind building for the selected three grids.This is the area where occurs the entrainment of clean air and is characterized by intensive turbulence.Therefore, very coarse grids may have difficulties to resolve the turbulence kinetic energy at the acceptable 80% for this area, which does not occur for the selected grids.Another area of intense turbulence is the shear layer between the roofs, where LES_IQv exceeds 85% for all resolutions, due to the dense grid.Outside the street canyon, LES_IQv has significantly lower values than inside the canyon, reaching 70% for the coarse grid in Figure 5a.Similarly, in the same area, the medium and the fine grid in Figure 5b,c present their lower values at 75% and 80%, respectively.On the other hand, this should not affect the quality when using the medium grid, because the main area of interest is inside street canyon.
For the coarse grid in Figure 5a, the LES_IQv is higher than 85% inside the canyon, while the lower values around 80% are in the outside area.As for the medium grid in Figure 5b, the turbulence is resolved at more than 85% inside the street canyon, while the percentage is over 90% for the walls and the shear layer.The value of LES_IQv is lower than the 80% (lower limit) at the outside area and approximately at one building height above the roofs.Again, this should be acceptable because turbulence is resolved at the areas that matter the most for this study.Finally, inside the street canyon the fine grid demonstrates values higher than 88%, while the values close to the wall reach to 95%, as shown in Figure 5c.
The use of LES_IQv provides valuable insights about the resolving of turbulence in the relatively high resolutions, usually utilized for the quasi-2D setup.The resolution of the coarse grid proved not sufficient for an LES study at the Re of 3,200,000.The medium grid performed satisfactory for the street canyon area, while showed some deficiencies for the upper region.The fine grid provided excellent compliance with the 80% limit for both areas.Nevertheless, the medium grid is judged to be a fair compromise between the required computational time and level of resolving turbulence.

Time-Step
Along with the spatial discretization, the selection of the appropriate time-step (Δt) must also ensure that the evolution of the small vortex structures will be described [68].Employing a structured grid generally leads to smaller time-step requirements.The small cells on the building walls combined with the hexahedral structured grid result in unnecessarily very small cells outside the street canyon, where velocity reaches its maximum.In turn, very small time-steps are required to maintain reasonably low Courant-Friedrichs-Lewy (CFL) number throughout the domain.A typical example of this situation can be found at Figure S1 of the supplementary material.
In this study, for the medium grid and Uref = 5 m‧s −1 , three different Δt values were examined, namely 6 × 10 −4 s, 8 × 10 −4 s, and 9 × 10 −4 s.Any value larger than 9 × 10 −4 s caused numerical instabilities, leading to diverging solutions.For the coarse grid in Figure 5a, the LES_IQ v is higher than 85% inside the canyon, while the lower values around 80% are in the outside area.As for the medium grid in Figure 5b, the turbulence is resolved at more than 85% inside the street canyon, while the percentage is over 90% for the walls and the shear layer.The value of LES_IQ v is lower than the 80% (lower limit) at the outside area and approximately at one building height above the roofs.Again, this should be acceptable because turbulence is resolved at the areas that matter the most for this study.Finally, inside the street canyon the fine grid demonstrates values higher than 88%, while the values close to the wall reach to 95%, as shown in Figure 5c.
The use of LES_IQ v provides valuable insights about the resolving of turbulence in the relatively high resolutions, usually utilized for the quasi-2D setup.The resolution of the coarse grid proved not sufficient for an LES study at the Re of 3,200,000.The medium grid performed satisfactory for the street canyon area, while showed some deficiencies for the upper region.The fine grid provided excellent compliance with the 80% limit for both areas.Nevertheless, the medium grid is judged to be a fair compromise between the required computational time and level of resolving turbulence.

Time-Step
Along with the spatial discretization, the selection of the appropriate time-step (∆t) must also ensure that the evolution of the small vortex structures will be described [68].Employing a structured grid generally leads to smaller time-step requirements.The small cells on the building walls combined with the hexahedral structured grid result in unnecessarily very small cells outside the street canyon, where velocity reaches its maximum.In turn, very small time-steps are required to maintain reasonably low Courant-Friedrichs-Lewy (CFL) number throughout the domain.A typical example of this situation can be found at Figure S1 of the supplementary material.
Following the works of Kornhass et al. [87] and Ai and Mak [88], the obtained velocity time-series for the three time-steps were transformed to the frequency domain using the fast Fourier transformation (FFT) technique.Figure 6 presents the spectra of turbulence energy generated for the U x and U z velocity components, for a point in the upwind side with x/W = 0.1 and z/H = 0.9.
Following the works of Kornhass et.al. [87] and Ai and Mak [88], the obtained velocity timeseries for the three time-steps were transformed to the frequency domain using the fast Fourier transformation (FFT) technique.Figure 6 presents the spectra of turbulence energy generated for the Ux and Uz velocity components, for a point in the upwind side with x/W = 0.1 and z/H = 0.9.For both velocity components in Figure 6a,b, the frequencies present only small differences in their minimum and maximum values for the three Δt values.This is expected, as the values of the Δt are close to each other.Six frequency scales from 10 −3 to 10 3 Hz are covered, reaching the significant high frequencies that represent the small-scale turbulence.Furthermore, the inertial sub-range of turbulence is covered between 1 and 10 Hz for both Ux and Uz, while turbulence is shown to dissipate faster for larger frequencies.
In Figure 6a, the power amplitudes for Ux exhibit small differences among the three selected Δt values.The coverage of the high frequencies occurs slightly more for the smaller Δt (2000 Hz), while the difference for the medium and large Δt is very small (1600 Hz instead of 1400 Hz).The same maximum frequencies are observed for Uz in Figure 6b.On the other hand, the calculated power amplitude is slightly larger for the smaller Δt = 6 × 10 −4 s, while the result for the medium and large Δt is the same.
For all cases of Δt and velocity components, the turbulence spectra present a coverage of the inertial sub-range of turbulence, up to higher frequencies in the dissipation range.As a fair compromise between the coverage of the small-scale turbulence and computational cost, the medium Δt = 6 × 10 −4 s is selected for the rest of the study.

Reynolds Number
The high values of wind velocity and building height result into high Reynolds numbers (Re) of atmospheric flows.Wind tunnel studies by Snyder [89] set the critical Re number at 11,000, based on the reference velocity of the free stream and building height.Later experiments either followed the aforementioned critical number [90], or directed to larger Re numbers, e.g., 19,000 [17], 37,000 [91], 44,000 [92], and 56,000 [3].Recently Cui et al. [93] reviewed the latest Re values reported in the literature and available CFD simulations with regard to the critical Re number, proposing a critical value of 32,000.At the same time, following such critical values do not result to validation issues since most available validation datasets have been produced in wind tunnels where 'low' Re is a constraint.To re-examine the effect of Re, we selected three reference wind velocities, namely 0.15 For both velocity components in Figure 6a,b, the frequencies present only small differences in their minimum and maximum values for the three ∆t values.This is expected, as the values of the ∆t are close to each other.Six frequency scales from 10 −3 to 10 3 Hz are covered, reaching the significant high frequencies that represent the small-scale turbulence.Furthermore, the inertial sub-range of turbulence is covered between 1 and 10 Hz for both U x and U z , while turbulence is shown to dissipate faster for larger frequencies.
In Figure 6a, the power amplitudes for U x exhibit small differences among the three selected ∆t values.The coverage of the high frequencies occurs slightly more for the smaller ∆t (2000 Hz), while the difference for the medium and large ∆t is very small (1600 Hz instead of 1400 Hz).The same maximum frequencies are observed for U z in Figure 6b.On the other hand, the calculated power amplitude is slightly larger for the smaller ∆t = 6 × 10 −4 s, while the result for the medium and large ∆t is the same.
For all cases of ∆t and velocity components, the turbulence spectra present a coverage of the inertial sub-range of turbulence, up to higher frequencies in the dissipation range.As a fair compromise between the coverage of the small-scale turbulence and computational cost, the medium ∆t = 6 × 10 −4 s is selected for the rest of the study.

Reynolds Number
The high values of wind velocity and building height result into high Reynolds numbers (Re) of atmospheric flows.Wind tunnel studies by Snyder [89] set the critical Re number at 11,000, based on the reference velocity of the free stream and building height.Later experiments either followed the aforementioned critical number [90], or directed to larger Re numbers, e.g., 19,000 [17], 37,000 [91], 44,000 [92], and 56,000 [3].Recently Cui et al. [93] reviewed the latest Re values reported in the literature and available CFD simulations with regard to the critical Re number, proposing a critical value of 32,000.At the same time, following such critical values do not result to validation issues since most available validation datasets have been produced in wind tunnels where 'low' Re is a constraint.To re-examine the effect of Re, we selected three reference wind velocities, namely 0.15 m•s −1 (Re = 96,000), 1.5 m•s −1 (Re = 960,000) and 5 m•s −1 (Re = 3,200,000).The normalized average magnitude of velocity (U mean /U ref ) and dimensionless concentration (C/C max ) are averaged over a time-span of 300 H/U ref , as shown in Figure 7, to employ equivalent mixing scales.
m⋅s −1 (Re = 96,000), 1.5 m⋅s −1 (Re = 960,000) and 5 m⋅s −1 (Re = 3,200,000).The normalized average magnitude of velocity (Umean/Uref) and dimensionless concentration (C/Cmax) are averaged over a timespan of 300 H/Uref, as shown in Figure 7, to employ equivalent mixing scales.The similarity criterion is generally satisfied for the normalised velocity, as shown in Figure 7ac, however, some differences also appear.In Figure 7b,c, the magnitude of the wind speed for the larger Re numbers at 960,000 and 3,200,000 is nearly identical, while for Re = 96,000 the height in which U/Uref = 0.5 is lower than the former two.As Re increases, the shear layer between the roofs is slight widening, while large vortex in the center of the canyon is reinforced, thus reducing the size of the three smaller ones.Moreover, the magnitude of velocity on the downwind wall is increasing along with velocity, meaning that the entrainment of clean air gets proportionally stronger with the increase of Re.
On the other hand, the impact on the pollutant dispersion is profound.The normalized concentration C/Cmax changes for the three Re numbers, as presented in Figure 7d-f.The patterns of pollution are similar, but the frequency of the small concentrations changes drastically from Re = 960,000 and 3,200,000 to Re = 96,000 (Figure S2 of the Supplementary Material).The larger values of C/Cmax at Re = 96,000 are observed near the ground, at the height of the source.At Re = 960,000 and 3,200,000 pollution gets more mixed inside the canyon, with higher concentrations observed at the upwind side of the canyon and the upwind wall, respectively.Finally, it is noted that for Re = 96,000, C/Cmax values close to zero are simulated for the downwind side of the canyon.The similarity criterion is generally satisfied for the normalised velocity, as shown in Figure 7a-c, however, some differences also appear.In Figure 7b,c, the magnitude of the wind speed for the larger Re numbers at 960,000 and 3,200,000 is nearly identical, while for Re = 96,000 the height in which U/U ref = 0.5 is lower than the former two.As Re increases, the shear layer between the roofs is slight widening, while large vortex in the center of the canyon is reinforced, thus reducing the size of the three smaller ones.Moreover, the magnitude of velocity on the downwind wall is increasing along with velocity, meaning that the entrainment of clean air gets proportionally stronger with the increase of Re.
On the other hand, the impact on the pollutant dispersion is profound.The normalized concentration C/C max changes for the three Re numbers, as presented in Figure 7d-f.The patterns of pollution are similar, but the frequency of the small concentrations changes drastically from Re = 960,000 and 3,200,000 to Re = 96,000 (Figure S2 of the Supplementary Material).The larger values of C/C max at Re = 96,000 are observed near the ground, at the height of the source.At Re = 960,000 and 3,200,000 pollution gets more mixed inside the canyon, with higher concentrations observed at the upwind side of the canyon and the upwind wall, respectively.Finally, it is noted that for Re = 96,000, C/C max values close to zero are simulated for the downwind side of the canyon.
The average normalized velocity for the three compared Re numbers showed discrepancies for the flow field, even though the notion of Re independency is generally confirmed.The results for concentration showed their clear dependency on Re number, as expected.Even though 3,200,000 is a high value of Re to simulate, the values for a real city geometry can be even higher, because of the large average height of the buildings.To keep the study as realistic as possible, the reference velocity at 5 m•s −1 was chosen.

Source Height
A brief discussion regarding the pollutant source is presented.Having to simulate the effect of traffic pollution, we assume a line source, even though area sources have also been reported in the literature [94,95].In OpenFOAM, the emission characteristics of the source are specified in the file fvOptions, inside the utility scalarSemiImplicitSource, while the topology is created with the topoSet utility [48].
The value for the mass emission rate was set to 620 g km −1 •h −1 and it represents a realistic emission rate of NO x for continuous canyon traffic [22].Li et al. [69] claimed that the laminar layer formed on the walls creates a bottleneck to the dispersion of pollutants.This resistance creates non-physical results, when the volume of the source is located entirely inside the laminar layer of the wall.Thus, the height and width of the source should be investigated.
In the present study, the source is located at the middle of the canyon, so its width is set to two cells, which differs slightly depending on the grid (see Figure 1).Six potential source heights were examined, namely 0.005, 0.01, 0.1, 0.5, 0.75, and 1 m.The mass emission rate is kept steady for all source heights, by setting the selection volumeMode as absolute inside the utility scalarSemiImplicitSource.The average normalized velocity for the three compared Re numbers showed discrepancies for the flow field, even though the notion of Re independency is generally confirmed.The results for concentration showed their clear dependency on Re number, as expected.Even though 3,200,000 is a high value of Re to simulate, the values for a real city geometry can be even higher, because of the large average height of the buildings.To keep the study as realistic as possible, the reference velocity at 5 m⋅s −1 was chosen.

Source Height
A brief discussion regarding the pollutant source is presented.Having to simulate the effect of traffic pollution, we assume a line source, even though area sources have also been reported in the literature [94,95].In OpenFOAM, the emission characteristics of the source are specified in the file fvOptions, inside the utility scalarSemiImplicitSource, while the topology is created with the topoSet utility [48].
The value for the mass emission rate was set to 620 g km −1 ⋅h −1 and it represents a realistic emission rate of NOx for continuous canyon traffic [22].Li et al. [69] claimed that the laminar layer formed on the walls creates a bottleneck to the dispersion of pollutants.This resistance creates non-physical results, when the volume of the source is located entirely inside the laminar layer of the wall.Thus, the height and width of the source should be investigated.
In the present study, the source is located at the middle of the canyon, so its width is set to two cells, which differs slightly depending on the grid (see Figure 1).Six potential source heights were examined, namely 0.005, 0.01, 0.1, 0.5, 0.75, and 1 m.The mass emission rate is kept steady for all source heights, by setting the selection volumeMode as absolute inside the utility scalarSemiImplicitSource.Initially, the source height was set to 0.1 m above the ground that clearly resides in the wall boundary, as shown in Figure 8a.As a result, the pollutant is weakly dispersed, but gathered at the bottom of the street.Then, the source height was increased to 0.5 m (Figure 8b), resulting in a dispersion pattern similar to the experimental data, such as those by Pavageau and Schatzmann [31].Finally, Figure 8c depicts the results of C/Cmax for the source height equal to 1 m, presenting a broader dispersion pattern.The results for 1 m improved the agreement with the experimental values for the downwind side, but increased the difference for the upwind side disproportionally.Thus, the source height of 0.5 m was selected as the optimal solution to both overcome the effect of the wall boundary layer and resemble the experimental data, even though the 1 m source height may also be considered a valid option.Initially, the source height was set to 0.1 m above the ground that clearly resides in the wall boundary, as shown in Figure 8a.As a result, the pollutant is weakly dispersed, but gathered at the bottom of the street.Then, the source height was increased to 0.5 m (Figure 8b), resulting in a dispersion pattern similar to the experimental data, such as those by Pavageau and Schatzmann [31].Finally, Figure 8c depicts the results of C/C max for the source height equal to 1 m, presenting a broader dispersion pattern.The results for 1 m improved the agreement with the experimental values for the downwind side, but increased the difference for the upwind side disproportionally.Thus, the source height of 0.5 m was selected as the optimal solution to both overcome the effect of the wall boundary layer and resemble the experimental data, even though the 1 m source height may also be considered a valid option.

Turbulence Modelling
Four SGS models (standard Smagorinsky, WALE, k-equation and localised dynamic k-equation) and three of the available wall functions (van Driest damping function, Spalding's law, and the simple logarithmic law) were tested in OpenFOAM v5.0.The considered scenarios are presented in Table 1.

Turbulence Modelling
Four SGS models (standard Smagorinsky, WALE, k-equation and localised dynamic k-equation) and three of the available wall functions (van Driest damping function, Spalding's law, and the simple logarithmic law) were tested in OpenFOAM v5.0.The considered scenarios are presented in Table 1.In Figure 9, there are no significant differences for the average Ux, between the applied models and the experimental data.The largest discrepancies are observed for Uz in the middle of the canyon (x/W = 0.50), which applies to all the employed combinations.Mismatch from the experimental data is also observed for Uz at x/W = 0.75, for the region between z/H = 0.4 and 0.8.For this case, the best prediction is obtained for KE-N, while the worst case represents WALE.Generally, the calculation for U was proved to be more accurate than that of W. This trend is also appeared for the velocity fluctuations, as shown in Figure 10.For this case, a disparity between the predictions and experiments appears as we move from x/W = 0.25 to 0.75.The capture of Ux,fluc at x/W = 0.25 is in a good agreement with the available experimental data for all combinations, while for Uz,fluc exists an over-prediction for the area between z/H = 0.1 and 0.7.For the line in the middle, the Ux,fluc is over-predicted and under-predicted at the regions above and below z/H = 0.7, respectively, while Uz,fluc is under-predicted for all the height.Finally, Ux,fluc at x/W = 0.75 is generally under-estimated along the height of the canyon, but close to the experiments.On the other hand, the Uz,fluc presents the worst discrepancies among the considered comparisons for the region higher than z/H = 0.7.In this area occurs the entrainment of clean air where the presence of turbulence is intensive and coincides with the area of lower LES_IQv (Figure 5b).The averages and fluctuations for Ux and Uz are compared with [30] in two more locations inside the canyon, at the Figure S3 and Figure S4 of the supplementary material.For the predicted averaged Ux and Uz fluctuations, DKE model's predictions present the largest mismatch compared to the experimental measurements, In Figure 9, there are no significant differences for the average U x , between the applied models and the experimental data.The largest discrepancies are observed for U z in the middle of the canyon (x/W = 0.50), which applies to all the employed combinations.Mismatch from the experimental data is also observed for U z at x/W = 0.75, for the region between z/H = 0.4 and 0.8.For this case, the best prediction is obtained for KE-N, while the worst case represents WALE.Generally, the calculation for U was proved to be more accurate than that of W.
This trend is also appeared for the velocity fluctuations, as shown in Figure 10.For this case, a disparity between the predictions and experiments appears as we move from x/W = 0.25 to 0.75.The capture of U x , fluc at x/W = 0.25 is in a good agreement with the available experimental data for all combinations, while for U z , fluc exists an over-prediction for the area between z/H = 0.1 and 0.7.For the line in the middle, the U x , fluc is over-predicted and under-predicted at the regions above and below z/H = 0.7, respectively, while U z , fluc is under-predicted for all the height.Finally, U x , fluc at x/W = 0.75 is generally under-estimated along the height of the canyon, but close to the experiments.On the other hand, the U z , fluc presents the worst discrepancies among the considered comparisons for the region higher than z/H = 0.7.In this area occurs the entrainment of clean air where the presence of turbulence is intensive and coincides with the area of lower LES_IQ v (Figure 5b).The averages and fluctuations for U x and U z are compared with [30] in two more locations inside the canyon, at the Figures S3 and S4 of the supplementary material.For the predicted averaged U x and U z fluctuations, DKE model's predictions present the largest mismatch compared to the experimental measurements, while KE-N model results are again approaching closer to the experimental data for all considered cases.
The most interesting observations obtained from the results of C/C max , some of which are presented in Figure 2. The concentrations for the downwind wall are very close to each other and under-estimated for all the used combinations.For the upwind wall, the results differentiate a lot, with the largest over-prediction obtained from D-KE and WALE.Results obtained from S-V, S-SP, K-VD, and KE-N present almost identical performance and closer to the experimental but still resulting in over-prediction of the concentration field, while the last model exhibits slightly better performance.The C/C max results of all the employed SGS models and wall functions can be found in of the Figure S5 supplementary material.
K-SP and KE-N models exhibit the best predictions, as shown in Figure 2. Their results for the downwind side are in fair agreement with the experiments, while the K-SP model gives better predictions for the heights below z/H = 0.4.However, KE-N achieved the best improvement for the upwind side, out-performing the other combinations.Hence, the use of the k-equation model, together with the simple logarithmic nutkWallFunction (KE-N), is proposed herein.

Validation Metrics
The validation metrics FAC2, NMSE, and FB are employed to obtain a quantitative measure of the quality of the selected KE-N simulation, which was determined herein as the most suitable option for the prediction of the average concentration (in Section 4.6).The averages and the fluctuations of velocity components, at 200 points, were compared with the available experimental data by Li et al. [30], while the LES results by Li [73] were also used for comparison.The concentration results were validated against the experimental data by Pavageau and Schatzmann [31], Meroney et al. [32], Pavageau [33], and Hoydysh and Dabberdt [74]; for a total of 24 points.For the following definitions, O i are the observed values in experiments, and P i are the predicted values from the obtained LES, while "[ ]" denotes the averaged values.
1. Factor of two or where n i = 1, for 0.5 ≤ P i /O i ≤ 2, or else n i = 0. 2. Normalized mean square error (NMSE) 3. Fractional bias (FB) All the employed metrics are summarized in Table 2.
FAC2 estimates the percentage of the predictions, which are within a factor of two of the observations, based on the ratio of the predicted and observed value [75].It was employed for U x and U z presenting a value of 0.85 for the horizontal velocity, the result for the vertical component is lower at 0.68 and the total FAC2 is 0.77.The ideal value is the unity, but values above 0.5 [77] and 0.8 [78] are also acceptable.Total NSME for both velocity components was found to be 0.25, while individually for U x and U z was calculated as 0.19 and 0.30, respectively.For this metric, the ideal value is zero, but NMSE < 1.5 [77] and NMSE < 4 [75] are likewise reported as acceptable limits.By these measures, the NMSE gives the higher level of confidence to our results.It is important to mention that the NSME does not have sense for parameters with negative and positive values [75], nevertheless, we tested it because 184 of the 200 points of our comparison have the same sign.
Regarding the concentration calculation, two ways are generally used.The first is the normalized concentration (C*), employed to account for Re and strength of the source, described by Schatzmann and Leitl [96] where C is the simulated concentration in Kg•m −3 , U ref is the reference wind velocity, H is the characteristic building height and Q/L is the strength of the source in Kg The second approach is the dimensionless concentration (C/C max ), where C is divided by the maximum value in the domain, usually residing at the upwind wall.C* is used most of the times, even though the normalizing quantities vary, such as H 2 instead of HL and volumetric flow of the source instead of mass rate [93].In most of our cases, we applied the dimensionless concentration C/C max , which was observed to reach a steady state value at around 500 H/U ref .
The three employed metrics for the dimensionless concentrations C/C max show excellent results.FAC2 is 0.8 for the upwind wall, while it reaches the absolute unity for the downwind wall.The NMSE is very close to the optimal zero, presenting 0.05 and 0.04 for the upwind and downwind walls, respectively.Finally, FB gave the largest values at 0.15 for the downwind wall, while it was around 0.06 for the upwind wall that satisfies the recommended levels.

Turbulence Statistics
The energy spectra for the middle of the canyon were computed to examine the quality of our LES approach.Several works have utilized energy spectra to examine the ability of LES to reproduce realistic results.These are compared either against wind tunnel data [17,24], or among different LES parameterizations (e.g., methods for the production of turbulence [97], grid resolutions, and Re number [98]).
parameterizations (e.g., methods for the production of turbulence [97], grid resolutions, and Re number [98]). Figure 11 presents the energy spectra for three heights in the middle of the canyon (x/W = 0.5), at z/H = 0.1, 0.5, and 1.0, where cells have relatively larger size compared to close the walls.For all the cases the amplitudes of the spectra correspond to the magnitude of the local velocity.The spectrum for Ux at the top of the canyon is parallel to the 5/3 slope described by Kolmogorov, between 0.3 and 2 Hz, indicating the existence of the inertial sub-range in a narrow band.On the other hand, for the regions with the dense grid near the wall (z/H = 0.1) and the lower velocities in the middle height (z/H = 0.5), the spectra follow the −5/3 slope throughout the six orders of frequency, except from an area between 1 and 10 Hz.The spectra for Uz are presented in Figure 11b and they show a very narrow inertial sub-range between 0.8 and 4 Hz.For higher frequencies and smaller turbulence scales the slope is steeper, implying that turbulence decays faster than expected.This may explain the inability of all combinations to predict Uz in the middle of the canyon (see Figure 9) and can be attributed either on grid resolution or the numerical schemes.

Conclusions
This study focused on the LES of wind flow and pollutant dispersion in an ideal quasi-2D street canyon, using the open-source CFD code OpenFOAM v5.0.Initially, the computational domain height, the grid resolution and time-step were studied, to examine the quality of the LES approach.Moreover, the effect of Re and the height of the pollutant source were investigated, to find out their possible impact to the validation results.After establishing the details about the grid, the Re number and the source setup, the effect of turbulence modelling was analyzed.Four well-known SGS models that are implemented in OpenFOAM were tested, together with three standard wall-functions, to cover for the excess values of y + .The results were compared with available experimental and numerical data.The main conclusions are: 1.The height of the domain indicated a small effect for the concentration in the downwind side, while had no effect for the upwind wall.Both effects diminished for the case with 5 H, which was selected for the height of the upper atmosphere.
2. The grid resolution with 120 cells per building height was selected.The use of the LES_IQv index showed that turbulence kinetic energy was resolved at more than 80% Figure 11 presents the energy spectra for three heights in the middle of the canyon (x/W = 0.5), at z/H = 0.1, 0.5, and 1.0, where cells have relatively larger size compared to close the walls.For all the cases the amplitudes of the spectra correspond to the magnitude of the local velocity.The spectrum for U x at the top of the canyon is parallel to the 5/3 slope described by Kolmogorov, between 0.3 and 2 Hz, indicating the existence of the inertial sub-range in a narrow band.On the other hand, for the regions with the dense grid near the wall (z/H = 0.1) and the lower velocities in the middle height (z/H = 0.5), the spectra follow the −5/3 slope throughout the six orders of frequency, except from an area between 1 and 10 Hz.The spectra for U z are presented in Figure 11b and they show a very narrow inertial sub-range between 0.8 and 4 Hz.For higher frequencies and smaller turbulence scales the slope is steeper, implying that turbulence decays faster than expected.This may explain the inability of all combinations to predict U z in the middle of the canyon (see Figure 9) and can be attributed either on grid resolution or the numerical schemes.

Conclusions
This study focused on the LES of wind flow and pollutant dispersion in an ideal quasi-2D street canyon, using the open-source CFD code OpenFOAM v5.0.Initially, the computational domain height, the grid resolution and time-step were studied, to examine the quality of the LES approach.Moreover, the effect of Re and the height of the pollutant source were investigated, to find out their possible impact to the validation results.After establishing the details about the grid, the Re number and the source setup, the effect of turbulence modelling was analyzed.Four well-known SGS models that are implemented in OpenFOAM were tested, together with three standard wall-functions, to cover for the excess values of y + .The results were compared with available experimental and numerical data.The main conclusions are: 1.
The height of the domain indicated a small effect for the concentration in the downwind side, while had no effect for the upwind wall.Both effects diminished for the case with 5 H, which was selected for the height of the upper atmosphere.

2.
The grid resolution with 120 cells per building height was selected.The use of the LES_IQ v index showed that turbulence kinetic energy was resolved at more than 80% for the area of the canyon, while the lower values at the street level are in the downwind side of the canyon.

3.
The use of energy spectra showed small differences in the coverage of the turbulence scales for the examined time-steps, as expected, because the time-steps had very close values.Smaller values were not tested.A CFL value below 0.1 was kept for the most period of the numerical simulation and the largest part of the computational domain.

4.
All the SGS models reproduced the average velocity at an acceptable level.Horizontal velocity U x is generally resolved better than vertical velocity U z .Some discrepancies were observed for the fluctuations of velocity, again the larger differences were occurred for U z . 5.
On the other hand, pollutant dispersion was most affected by the SGS model, with best results obtained with the k-equation model and nutKwallfunction.

6.
The height of the source can be critical to the dispersion of the pollutants and subsequently the validation process.If the selected height is located inside the laminar layer of the wall boundary layer on the ground, then dispersion will be artificially hindered.The selected source height overcame this issue and resulted to a better agreement with the experimental data.However, this is a known challenge of numerical models and a recurring issue in the literature, that requires better consideration in future experimental studies.7.
Three Re numbers were studied to examine its effect on the results.As expected, the normalized average velocity was nearly identical for the three cases.On the other hand, the average dimensionless concentration was affected and not proportionally to the Re number increase.
The y + at the ground reached an average value of 2.7 and a maximum value of 6, while the respective values for the lateral building walls at 3.7 and 21.These values do not cover completely the proposed maximum of y + = 5.On the other hand, the use of the selected wall function and the results of the validation metrics imply that these y + values can be acceptable for Re = 3,200,000.8.
The discrepancies for U z in the middle of the canyon (Figure 9) can be attributed to the local grid resolution or the use of dissipative numerical schemes.In our case, second order schemes were applied for all terms, while the cells are relatively larger for this area.Furthermore, compared to the medium the fine grid showed a clear improvement, moving the slope of the energy spectrum for U z , closer to the −5/3 slope described by Kolmogorov (Figure S6).This implies that a higher grid resolution may indeed improve the results in this region.
LES is a valuable tool and upon the selection and validation of the most appropriate configurations (e.g., turbulence schemes, grid resolution) could supplement and extended experimental work.Similarly to other modelling efforts, LES simulations could be used to enhance understanding of physical phenomena, and contribute on the development of new or improvement of existing numerical models and software.Such software could be capable to provide predictions or extensive hindcasts in almost real-time, thus their improvement could expand their applicability for emergencies.

Figure 1 .
Figure 1.The quasi-2d street canyon geometry: (a) overview of the physical geometry and actual buildings (light shaded block), (b) detail of the modelled geometry and boundary conditions (discussed later); and (c) xz cross-section of a computational grid example.

Figure 1 .
Figure 1.The quasi-2d street canyon geometry: (a) overview of the physical geometry and actual buildings (light shaded block), (b) detail of the modelled geometry and boundary conditions (discussed later); and (c) xz cross-section of a computational grid example.

27 Figure 2 .
Figure 2. Vertical dimensionless concentration profiles for the upper atmosphere height and grid resolution: (a) the downwind and (b) the upwind wall for the unity street canyon.

Figure 2 .
Figure 2. Vertical dimensionless concentration profiles for the upper atmosphere height and grid resolution: (a) the downwind and (b) the upwind wall for the unity street canyon.

Figure 3 .
Figure 3.The calculated GCIfine for the medium to fine grid refinement with the corresponding dimensionless Ux and Uz at x/W = 0.25, x/W = 0.50, and x/W = 0.75.

Figure 3 .
Figure 3.The calculated GCI fine for the medium to fine grid refinement with the corresponding dimensionless U x and U z at x/W = 0.25, x/W = 0.50, and x/W = 0.75.

Figure 4 .
Figure 4. Calculated GCIfine for the medium to fine grid refinement and corresponding averaged and dimensionless concentration C/Cmax at (a) the downwind and (b) the upwind walls.

Figure 4 .
Figure 4. Calculated GCI fine for the medium to fine grid refinement and corresponding averaged and dimensionless concentration C/C max at (a) the downwind and (b) the upwind walls.

Figure 6 .
Figure 6.Comparison of the turbulence energy spectra in the frequency domain of velocity: (a) Ux and (b) Uz at x/W = 0.1 and z/H = 0.9 on the upwind side of the street canyon.The solid red time series represent the Δt = 6 × 10 −4 s, blue the Δt = 8 × 10 −4 s and green the Δt = 9 × 10 −4 s.

Figure 6 .
Figure 6.Comparison of the turbulence energy spectra in the frequency domain of velocity: (a) U x and (b) U z at x/W = 0.1 and z/H = 0.9 on the upwind side of the street canyon.The solid red time series represent the ∆t = 6 × 10 −4 s, blue the ∆t = 8 × 10 −4 s and green the ∆t = 9 × 10 −4 s.

Figure 8
presents the normalized C/C max values after 900 s for three indicative source heights.Atmosphere 2018, 9, x FOR PEER REVIEW 16 of 27

Figure 8
presents the normalized C/Cmax values after 900 s for three indicative source heights.

Table 1 .
Combinations of the employed SGS models and wall functions.Combination Abbreviation standard Smagorinsky model and van Driest dumping function S-VD standard Smagorinsky model and Spalding wall function S-SP standard Smagorinsky model and nutkWallFunction S-N WALE model WALE k-equation model and van Driest dumping function KE-VD k-equation model and Spalding wall function KE-SP k-equation model and nutkWallFunction KE-N localized dynamic k-equation model DKE The results of the simulations are shown in Figures 9 and 10, presenting the averages and fluctuations of Uand W-components, respectively.Three vertical lines at x/W = 0.25, 0.50, and 0.75 of the canyon width were chosen for the comparison of numerical data obtained with the available experimental data by Li et al. [30].Atmosphere 2018, 9, x FOR PEER REVIEW 17 of 27

Table 1 .
Combinations of the employed SGS models and wall functions Combination Abbreviation standard Smagorinsky model and van Driest dumping function S-VD standard Smagorinsky model and Spalding wall function S-SP standard Smagorinsky model and nutkWallFunction S-N WALE model WALE k-equation model and van Driest dumping function KE-VD k-equation model and Spalding wall function KE-SP k-equation model and nutkWallFunction KE-N localized dynamic k-equation model DKEThe results of the simulations are shown in Figure9and Figure10, presenting the averages and fluctuations of Uand W-components, respectively.Three vertical lines at x/W = 0.25, 0.50, and 0.75 of the canyon width were chosen for the comparison of numerical data obtained with the available experimental data by Li et al.[30].

Figure 9 .
Figure 9. Normalized averages for the velocity components Ux and Uz for a unity street canyon at: (a) and (d) x/W = 0.25; (b) and (e) x/W = 0.50; and (c) and (f) x/W = 0.75.

Figure 9 .
Figure 9. Normalized averages for the velocity components U x and U z for a unity street canyon at: (a) and (d) x/W = 0.25; (b) and (e) x/W = 0.50; and (c) and (f) x/W = 0.75.

Figure 10 .
Figure 10.Normalized fluctuations for the velocity components Ux and Uz for a unity street canyon at: (a) and (d) x/W = 0.25; (b) and (e) x/W = 0.50; and (c) and (f) x/W = 0.75.

Figure 10 .
Figure 10.Normalized fluctuations for the velocity components U x and U z for a unity street canyon at: (a) and (d) x/W = 0.25; (b) and (e) x/W = 0.50; and (c) and (f) x/W = 0.75.

Figure 11 .
Figure 11.Turbulence energy spectra of the (a) Ux and (b) Uz velocity components, for three heights in the middle of the domain x/W = 0.5, at z/H = 0.1, 0.5, and 1.0.

Figure 11 .
Figure 11.Turbulence energy spectra of the (a) U x and (b) U z velocity components, for three heights in the middle of the domain x/W = 0.5, at z/H = 0.1, 0.5, and 1.0.

Table 2 .
Calculated and ideal values of the employed validation metrics and experimental datasets (200 points for the velocity components and 24 for the concentration).