Numerical Modelling of Air Pollutant Dispersion in Complex Urban Areas: Investigation of City Parts from Downtowns Hanover and Frankfurt

: Hazardous gas dispersion within a complex urban environment in 1:1 scaled geometry of German cities, Hanover and Frankfurt, is predicted using an advanced turbulence model. The investigation involves a large group of real buildings with a high level of details. For this purpose, Computer Aided Design (CAD) of two conﬁgurations are cleaned, then ﬁne grids meshed in. Weather conditions are introduced using power law velocity proﬁles at inlets boundary. The investigation focused on the e ﬀ ects of release locations and material properties of the contaminants (e.g., densities) on the convection / di ﬀ usion of pollutants within complex urban area. Two geometries demonstrating di ﬀ erent topologies and boundaries conditions are investigated. Pollutants are introduced into the computational domain through chimney and / or pipe leakages in various locations. Simulations are carried out using Large Eddy Simulation (LES) turbulence model and species transport for the pollutants. The weather conditions are accounted for using a logarithmic velocity proﬁle at inlets. CH 4 and CO 2 distributions, as well as turbulence quantities and velocity proﬁles, show important inﬂuences on the dispersion behavior of the hazardous gas. display an instantaneous concentration of CH 4 and CO 2 at the ground, respectively. The CO 2 concentration is one order of magnitude larger than that of CH 4 . This concentration di ﬀ erence demonstrates that the deposition of pollutants is mainly controlled by the density ratio and not by the air-velocity.


Introduction
A growing concern about air pollution in urban environments results in stringent legislation for devices certification and energy utilization [1]. Several studies have been interested in the pollutant emissions reduction and the optimization of combustion systems by the use of new burners [2,3] or by capturing target species [4]. To assess the problem in detail, accurate prediction of the transport and dispersion of airborne contaminants in urban environments is needed. Several methods have been used to investigate air pollution [5,6]. Wind tunnels only provide an understanding of the dynamics and scaled data. Such approaches are able to perform parameters sampling points and to deliver well-documented data sets for the evaluation of numerical models [7]. Yet, due to scale limitations, it is a challenge for such models to fully be tested in large scale atmospheric turbulent flow conditions [8]. Numerical approaches have numerous advantages, e.g., high controllability, lower cost, three-dimensional (3D) data output, etc. [9].
Many preceding studies showed that Computational Fluid Dynamics (CFD) simulations, based on the Reynolds averaged Navier-Stokes (RANS) equations, delivered unaccrued results in reproducing the wind-flow pollutant dispersion around buildings [10][11][12][13][14]. Britter and Hanna [11] pointed out that the RANS turbulence model can be used in simple geometry (isolated building) to produce reasonable qualitative results for mean flows. When compared with laboratory or field experiments, this model showed a better performance than that of simple operational models. Blocken et al. [13] tested three cases of tracer gas dispersion in a neutrally stable atmospheric boundary layer and compared their results against wind tunnel full-scale measurements. For all three cases, they concluded that RANS simulations significantly underestimate the cross-wind (lateral) dispersion of the pollutants. Tominaga and Stathopoulos [14] numerically studied the dispersion around an isolated cube. They found that RANS models underpredict horizontal concentrations at the roof of simulated buildings.
To reach more accurate results, several researchers have turned to the Large-Eddy Simulation (LES) approach, which demonstrated its capability and potential in microscale atmospheric dispersion modelling [15][16][17][18][19][20][21][22][23][24][25]. Liu and Barth [15] studied the flow in 3D streets using scalar transport. Their results showed the limitations of eddies developments with sizes larger than 0.5H (H height of the building) in the so-called free surface layer and produce less turbulent kinetic energy (TKE) above the rooftop. Sada and Sato [18] simulated the stack-gas diffusion around a cubical building in order to predict the instantaneous concentration fluctuation of a plume. Their results were compared with data from wind tunnel experiments and showed a good agreement. Liu et al. [25] used the one-equation subgrid-scale (SGS) parameterization and the LES model to study the flow and pollutant transport in 2D urban street canyons. Their results showed that the combined model is an appropriate method for predicting wind field and pollutant dispersion in the crowded urban area.
Recent works [26][27][28][29][30] focused on the study of more complex geometries. Gausseau et al. [26] used LES to investigate the effect of two wind directions on near-field pollutant dispersion from stacks located at a low-rise building in downtown Montreal Canada. Their results showed some similarities with the case of dispersion around an isolated building. Van Hoof and Blocken [28] studied the CO 2 dispersion in a semi-enclosed stadium and validated their results using onsite measurements. The validated model was used to evaluate the carbon dioxide concentration in the stadium. Amorim et al. [29] chose two European urban areas (Lisbon and Aveiro, Portugal) to evaluate the effect of trees over the dispersion of carbon monoxide (CO) emitted by road traffic. They found that trees are responsible for an increase of 12% CO concentration in the first configuration (nonaligned wind 45 • ). In the second configuration (aligned wind), they observed a decrease of 16% of CO concentration.
In light of presented studies, there are many qualitative, as well as quantitative parameters, which need to be studied when investigating pollutants dispersion in an urban area, especially for complex geometries. Hence, the present study aims at providing more insight into the dispersion process by analyzing the effect of buildings arrangement and the material properties (such as density) on the dispersion of hazardous gases. For this purpose, two different configurations, which exhibit complex geometries of Hanover and Frankfurt cities, Germany, are investigated. The most challenging point is the high density of the surrounding buildings, which makes the carrier phase flow and pollutants concentration field complex, anisotropic, and difficult for prediction. Different scenarios of pollutant injection into the computational domain are tested. Dense buildings are included up to a distance of 600 m in diameter and a high-resolution refined grid with over 40 million cells was used. LES simulations are performed, then results are compared against published data from Leitl, B. and Schatzmann. [30] in order to validate the numerical models [31].

Turbulence Model
In LES, the flow variables are decomposed into the resolved-scale components and the (SGS) components using a filtering operation [32]. The filtering process effectively filters out eddies whose scales are smaller than the filter width or grid spacing used in the computations. Resulting in filtered, continuity and momentum equations, the incompressible Navier-Stokes equations are given as follows: where bars and over-bars express mean and filtered quantities. The variables u, ρ, p, g, and δ ij denote the velocity, the density, the hydrostatic pressure, the gravitational force, and the kronecker delta, respectively. The quantity ν represents the molecular viscosity and S is the source term. The Smagorinsky SGS model is applied to close the system of equations and determine the SGS stresses τ SGS ij using the Equation (3)

Dispersion Model
Numerous models are accessible to predict the dispersion of airborne material. In this study, pollution dispersion was modelled using the species transport approach by calculating convection-diffusion as described by Equation (4). This equation is solved for N-1 species, where N is the total number of species in the fluid flow. The sum of the mass fractions of all species must be equal to unity, the Nth mass fraction is given from the sum of the N-1 solved mass fractions.
The conservation equation takes the following general formula: where ρ is the density, R i is the net rate of production of species by chemical reaction, S i is the rate of creation by addition from the dispersed phase, D i is the diffusion flux of species, Y i is the mass fraction of specie I, and u j is the flow velocity.

Numerical Procedure
The commercial CFD software ANSYS FLUENT 17.0 is employed to simulate the airflow and pollutants dispersion. For all configurations, the dispersion of the pollutant is simulated using LES and dynamic Smagorinsky, SGS, along with the species transport model. The bounded central-differencing scheme is used to discretize the momentum equation. A second order upwind scheme is applied for the pollutant concentration and energy equation. Time integration is bounded to second order implicit and pressure interpolation is second order. For both cases, the time step advancement is set to ∆t = 10 −3 s, which produces a Courant number "C" less than one within the computational domain. The presented results are averaged over a period of 220 s, which corresponds to five flow-through times.

Validation of the Numerical Model
In order to validate the numerical model, the (CEDVAL B1-1) configuration [30] (Figure 1) is investigated using the same boundaries conditions and physical models. This study was conducted in the BLASIUS wind tunnel at the Meteorological Institute of the University of Hamburg [30]. A modified Standing-Spires and uniform LEGO-roughness was used to generate a neutral atmospheric boundary layer ABL at a scale of 1:200 with Re = 37,252. First, the boundary layer flow was validated based on detailed flow measurements and comparison with full-scale data. Then, the building models were placed in the test section. Inside the wind tunnel, a 3 × 7 array of buildings composed of rectangular blocks with the same dimensions H = 0.125 m are used to model the dispersion of carbon dioxide (CO 2 ). The CO 2 was emitted by four sources, located at the leeward of building source, with exhaust velocity we = 0.033 m/s. Validation is achieved by comparing LES results against the published measurement [30]. The horizontal plane at z = 1.5 m (full-scale) height was used as a measurement plane in which pollutants data were acquired. The same experiment is also simulated by Longo et al. [31]. Identical computational domain, as well as boundary conditions, are used for the present work. The inlet boundary is set to 1 m upstream of the buildings, in which the ABL profiles are measured within the wind tunnel. The outlet boundary is located at 4 m downstream of the last array of the buildings. The width, depth and height of the domain are 1 m, 5 m, and 1 m, respectively, which correspond to the wind tunnel size. In the simulation, a power law velocity profile is used at the inlet, a zero-pressure gradient is used at the outlet and top, no-slip velocity is used at the ground and walls, and symmetry is used for lateral boundary conditions. All the boundaries and initial conditions are summarized in Table 1. Although the temperature is isothermal (300 k), due to thermal stratification, the buoyancy is not present. Yet, we do have stratification due to density difference between carbon dioxide and air. An unsteady, 3D, double precision, pressure-based solver with second order schemes for momentum, turbulence quantities, and concentrations are used. The coupled algorithm is selected for pressure-velocity coupling. A transport equation approach is used for the dispersion study. A grid-independent performance of the numerical simulation results is verified by comparing the results of test cases having different grid numbers. Three meshes with different levels of resolution are investigated. A fine resolved mesh of 10 million elements is used for all computations. Simulations are performed using parallels computation of 40 CPUs. The total flow period is set to 12.0 s, which allows the flow to sweep the domain five times. Figure 2 shows a comparison between the dimensionless turbulent kinetic energy, dimensionless velocity components, and the dimensionless concentration of CO 2 obtained by numerical simulation against the wind tunnel experimental data.
This figure demonstrates a good agreement between LES numerical results and measurements. The values of turbulent kinetic energy, velocity components, and CO 2 concentration are well-predicted and follow the experimental measurements, demonstrating the ability of LES Smagorinsky, SGS, and species transport models in predicting the pollutants dispersion in atmospheric conditions.

Geometry and Boundaries Conditions
In order to represent real case scenarios, two computational domains are generated, one for each case (Figure 3a,b). In case 1 (Hanover city), the focus is placed on the effect of buildings arrangement, whereas in case 2 (Frankfurt city), the behaviors of two hazardous gases with different densities are studied. At the inlet, flow is perpendicular to the injection plane. The streamwise, spanwise, and vertical coordinates are represented by x, y, and z, respectively. The buildings are modeled based on the existing full-scale data. To limit the number of cells, vegetation is not considered.  In case 1, shown in Figure 3a, results highlight the effect of buildings arrangement. For this purpose, two different scenarios dealing with different release locations are investigated. The first one is an accidental leakage of methane. The hazardous gas (methane) is introduced through an opening located at a confined geometry (x = −27 m, y = 26 m, z = 1 m). The choice of this location is based on its closeness to the buildings and/or distance to the exit. CH 4 is injected toward the wall of a building from an inlet having 1 m 2 , a velocity of 8 m/s, and a temperature of 310 K. The second leakage scenario focuses on the dispersion of methane within an open surface, i.e., along the street and far from buildings (x = 110 m, y = −29 m, z = 1 m). Identical boundary conditions for methane are applied at the location 2. For both leakages, CH 4 is injected with 100% mass-fraction. Figure 3b shows the effect of gas density for case 2. Two gases having different densities, i.e., methane and carbon dioxide, are used for the investigation. The choice of these gases is based on the density ratio to the ambient air. CO 2 is heavier than ambient air (ρ (CO 2 ) = 1.84 kg/m 3 ), while the methane is lighter than ambient air (ρ (CH 4 ) = 0.71 kg/m 3 ). The CO 2 and CH 4 gases are injected from a chimney sited at (x = 0 m, y = 0 m, z = 20 m) by applying an inlet velocity of 2 m/s and a temperature equal to 310 K.
For the carrier phase boundary conditions, a power law profile is used to describe the variation of the inlet wind speed as a function of the height given in Equation (5) [33] where U ref = 9 m/s is the velocity reference at z ref = 10 m and a = 0.21.
The turbulence was generated using the inlet boundary condition. For LES simulation, the vortex method was used with a number of vortices Nv = 190 in order to establish a time-dependent inlet profile for the inlet velocity. In this method, a number of vortices Nv are generated and convected randomly at each time step.
With this approach, a perturbation is added on a specified mean velocity profile via a fluctuating vorticity field (i.e., two-dimensional in the plane normal to the streamwise direction). The vortex method is based on the Lagrangian form of the 2D evolution equation of the vorticity and the Biot-Savart law. A particle discretization is used to solve this equation. These particles, or "vortex points',' are convected randomly and carry information about the vorticity field. If N is the number of vortex points and A is the area of the inlet section, the amount of vorticity carried by a given particle i is represented by the circulation Γ i and an assumed spatial distribution η: where k is the turbulence kinetic energy. The parameter σ provides control over the size of a vortex particle. The resulting discretization for the velocity field is given by where → z is the unit vector in the streamwise direction. Originally, the size of the vortex was fixed by an adhoc value of σ. To make the vortex method generally applicable, a local vortex size is specified through a turbulent mixing length hypothesis.
The turbulent kinetic energy, k, and the turbulence dissipation rate, ε, used for RANS simulations are determined using the Equations (9) and (10) [26]: I u is the streamwise turbulence intensity (I u = σ u /U, with σ u the standard deviation of u) The variable u * is the friction velocity, κ = 0.41 represents the von Karman constant and z 0 = 1 m (full scale) is the aerodynamic roughness length. All the boundaries and initial conditions are summarized in Table 2.  Figure 4a exhibits the mesh used for Hanover city simulations. A mesh having 40 million elements is generated by the patch independent algorithm. This algorithm is based on the top-down approach, in which volumes are meshed first and cells are projected to faces and edges. It is mainly used for complex geometries that demonstrate various features having different tolerances. In order to well-capture the boundary layer and well=predict the continuous phase physics, an appropriate refinement near the source of CH 4 and close to the ground was performed. The grid quality is improved by increasing the cell number. Figure 4b shows the grid used for the Frankfurt city case. In the building zones, the high-resolution computational grids are composed of prism and tetra elements only. Faraway from the building areas, a hexahedral mesh is applied. For both meshes, the resolution is determined using a grid-sensitivity analysis. It is worth mentioning that grid quality varies between 0.05 and 0.95. The metrics used to quantify the quality is the orthogonality, which relates the angles between adjacent elements. The optimum angle is 90 • for quadrilateral faced cells and 60 • for triangular faces elements. The mesh quality plays a significant impact on the accuracy of the numerical prediction, as well as the stability of the simulation.   (Figure 5b,c). The maximum velocity magnitude, 7 m/s, is at the central location section, which is the narrowest passage along the street. This is explained by the conservation of mass, i.e., if a fluid flow having a constant flow rate passes through smaller section, it will be accelerated. As the flow is incompressible, the transport and dispersion of pollutant species is governed by the wind velocity. Behind buildings, zones of recirculation are noticed. Figure 5b,c depict the resulting vertical profiles of the mean streamwise velocity <u> and the turbulence intensity (I), obtained at (y = −250 m) using LES model. The low power profile is well-reproduced and clearly shows the frictional effect of the ground. For the turbulence intensity, one notices a low variation from 4% (z = 4 m) to 6% (z = 200 m), which are smaller values compared to the literature [23] that predicted 12%. Consequently, one concludes that inlet turbulence has a relatively low influence on the dispersion of pollutants. In fact, the turbulence at the release location is mainly governed by the presence of the surrounding buildings, which generate velocity gradients, in turn enhancing the turbulence.  Figure 6 shows qualitative results of methane dispersion. Two different injection locations are investigated. The release location significantly affects the dispersion of pollutants. When methane is injected for a time period of 220 s, the burnable methane-air mixture, represented by an iso-surface equal 4% of CH 4 -mass fraction, reaches the wall of the nearby building (Figure 6b). The burnable mixture shows a well-developed vertical plume. Parts of this mixture are in touch with the ground. In this zone, the buildings interfere with the wind velocity, delivering a flow that is mainly characterized by the urban canopy and loosely dependent on atmospheric wind flow. As a result, the methane concentration increases in the wakes of the building where wind speed is considerably reduced and the recirculation zone becomes remarkably larger. Figure 6a depicts the concentration on buildings surfaces and demonstrates a critical value of 4% close to the wall of the nearby building. From a practical point of view, if a leakage occurs at this location, the ventilation intakes of these buildings should preferably not be located on this face. Indeed, pollutants can contaminate the indoor air of the nearby building if the windows of the leeward faces are open.  The turbulence intensity along the y-direction is more important than that in all other directions, making the species iso-surfaces show a narrowly opened cone (Figure 6c). Figure 6d shows the streamlines starting from the CH 4 source (location 2), in which the velocity of the carrier phase is the highest and exceeds 7 m/s. CH 4 mass fraction shows different pattern compared to location 1. Pollutants are transported horizontally more than that in the vertical direction. The important carrier phase momentum demonstrates a dominant effect. The horizontal time scale of the mixture momentum is clearly smaller than that in the vertical direction. Despite the importance of buoyancy forces, it appears that the dynamic/momentum of the carrier phase is the parameter that controls the motion and dispersion of CH 4 . Figure 7a,b depicts the quantitative result of the CH 4 dispersion obtained in location 1-simulation. The hazardous zone of an explosion and/or toxicity (mass fraction between 0.04 and 0.14) is bounded between (29 m < x < −28 m), (26 m < y < 27 m) and (2 m < z < 10 m). At the level z = 5 m, the curve bandwidth is less important in "y-direction" than that in "x-direction". The flow diffusion (dominated by the turbulence) is clearly not isentropic. This can be explained by the vortex effect behind buildings, which transports methane in the positive x-direction. Faraway from the source (z = 10 m), the concentration of methane decreases significantly to reach 0.1. Yet, it is strong enough to cause an explosion. Figure 7d,c shows the quantitative distribution of CH 4 dispersion produced in the simulation of location 2-case. The hazardous zone is bounded between (109 m < x < 111 m), (−29 m < y < 24 m) and (0 m < z < 2 m). Moving toward the outlet in the positive y-direction, the concentration decreases considerably and falls below 10% after 5 m. This demonstrates that the wind flow transports rapidly CH 4 and prevents a pollutants accumulation.  Figure 8a shows an instantaneous velocity contour of air at level z = 10 m. The velocity magnitude varies between zero and 9 m/s. Close to the buildings, obvious stagnation zones are registered. Velocity gradients are also observed downstream the buildings. Figure 8b depicts the streamlines obtained at level z = 10 m. The streamlines demonstrate a separation in front of the first raw of buildings and an acceleration of the carrier phase. The form of streamlines changes as we move toward the center of the geometry, and recirculation zones are observed behind buildings. Figure 8c,d display an instantaneous concentration of CH 4 and CO 2 at the ground, respectively. The CO 2 concentration is one order of magnitude larger than that of CH 4 . This concentration difference demonstrates that the deposition of pollutants is mainly controlled by the density ratio and not by the air-velocity.

Frankfurt City Results
The maximum concentration of both gases is observed behind the source building. A significant pollutant concentration is observed in the vicinity of computational domain outlet. This is due to the horizontal convective transport of the pollutant toward the outlet due to high wind velocity and negligible diffusion.  Figure 9 depicts the CO 2 and CH 4 mass fraction variations along vertical lines defined by (x, y) coordinates. Six different locations are chosen to investigate the effect of the density on the dispersion of pollutants. They are located behind the building source sited at (x = 0 m, y = 0 m). All profiles of CO 2 mass fraction show larger concentration than CH 4 . This is caused by the buoyancy force. The heavier gas is dragged downward (to the ground) by the effect of gravitational force. Carbon dioxide CO 2 is heavier than ambient air. As a result, its higher concentration is close to the ground. Pollutant transport is also influenced by the vortex located behind the building source (see Figure 9c). This zone is created by the separation of fluid flow from the edges of the buildings. The appearance of this zone strongly affects the transport and diffusion of pollutants. Due to the recirculation of the carrier phase (air), methane is dragged downward, and CH 4 is accumulated at the ground.
A maximum concentration of CO 2 mass fraction is recorded at (x = 11 m, y = 0), which is the closest location to the wall boundary. In this region, the effect of the non-slip condition reduces the velocity magnitudes, resulting in the formation of stagnating zone in which pollutants are collected. Faraway from the recirculation zone (x = 13 m, y = −17 m) and (x = 13 m, y = 17 m), the concentration of pollutants, CO 2 and CH 4 , is very low and mostly negligible. Figure 10 shows the effect of the density ratio on the dispersion of CH 4 and CO 2 . The density of pollutant-air mixture increases (CH 4 )/decreases (CO 2 ) rapidly and reaches the density of ambient air at a height of z = 5 m, far from the injection source (Figure 10a). The density ratio controls the shape and the volume of the pollutant plume. Figure 10b,c show iso-surfaces of CH 4 and CO 2 concentration. In order to compare the dispersion of pollutant gases, an iso-value of 4% mass fraction is represented. The plume of CH 4 is first transported vertically by the inlet flows and then carried by the wind horizontally to cover the roof of nearby buildings. In the case of CO 2 , the dispersion of the plume is located immediately at the top of the building that include the pollutant source. Compared to CH 4 , the dispersion is remarkably reduced.

Conclusions
High-resolution CFD simulations of near-and far-field pollutant dispersion in a building group in downtown Hanover and Frankfurt are predicted using LES turbulence model along with the dynamic Smagorinsky SGS. The simulations focus on the concentration of pollutants in the vicinity of the release sources and the spatial distribution of the transported species following different building structures and boundary conditions. The following remarks can be drawn: (1) Smagorinsky LES shows reasonable results of the velocities and turbulence quantities compared to experimental data of generic building configuration. The results agreement is used as validation of the flow model to investigate a larger full-scale 1:1 city geometry. (2) Methane dispersion is controlled by the air-CH 4 density ratio in cases of pollutant injection close to buildings. The buoyancy forces are responsible for the stratification of gases. CH 4 species are first transported in the vertical direction then dragged by the carrier phase toward the outlet. The velocity direction of pollutant injection and magnitude, as well as temperature difference between CH 4 and ambient air, do not play important roles. (3) If CH 4 is injected far from buildings, i.e., in open area in which carrier phase velocity is important, then the transport of pollutant showed different behavior. CH 4 was dragged horizontally much faster than in the vertical direction. It followed the street and kept almost a constant altitude until reaching the outlet.
(4) The main parameter that controls the dispersion in Frankfurt city configuration is the density ratio of CO 2 and CH 4 to ambient air. Despite the injection of both gases in the vertical direction, only CO 2 shows a significant concentration at the ground.
Future works should focus on the existence of vegetation, as well as the implementation of new sources that account for additional species creation, e.g., CO, NO x , SO x , which describe the transport means in the cities.