Lattice Boltzmann Method-Based Simulations of Pollutant Dispersion and Urban Physics

Mesocale atmospheric flows that develop in the boundary layer or microscale flows that develop in urban areas are challenging to predict, especially due to multiscale interactions, multiphysical couplings, land and urban surface thermal and geometrical properties and turbulence. However, these different flows can indirectly and directly affect the exposure of people to deteriorated air quality or thermal environment, as well as the structural and energy loads of buildings. Therefore, the ability to accurately predict the different interacting physical processes determining these flows is of primary importance. To this end, alternative approaches based on the lattice Boltzmann method (LBM) wall model large eddy simulations (WMLESs) appear particularly interesting as they provide a suitable framework to develop efficient numerical methods for the prediction of complex large or smaller scale atmospheric flows. In particular, this article summarizes recent developments and studies performed using the hybrid recursive regularized collision model for the simulation of complex or/and coupled turbulent flows. Different applications to the prediction of meteorological humid flows, urban pollutant dispersion, pedestrian wind comfort and pressure distribution on urban buildings including uncertainty quantification are especially reviewed. For these different applications, the accuracy of the developed approach was assessed by comparison with experimental and/or numerical reference data, showing a state of the art performance. Ongoing developments focus now on the validation and prediction of indoor environmental conditions including thermal mixing and pollutant dispersion in different types of rooms equipped with heat, ventilation and air conditioning systems.


Introduction
The capability to accurately predict urban physics and environmental quality for citizens via numerical simulation is nowadays a critical challenge, since it is a key tool for designing future optimized and sustainable urban areas. Such predictive models should account for a very broad range of physical mechanisms, ranging from large-scale meteorological effects to very small scale unsteady fluctuations of physical quantities such as temperature, air velocity, pressure, humidity, etc.
Among the different interacting physical phenomena occurring in the atmosphere at the meso-and microscales are orographic effects, land sea breeze, urban heat islands, deep convection, thunderstorms, convection, thermals, building wakes and turbulence (see Figure 1 in Schlünzen et al. [1]). Mesoscale atmospheric flows that develop in the boundary layer or microscale flows that develop in urban areas are thus very complex, especially due to multiscale interactions, multiphysical couplings, land and urban surface thermal and geometrical properties and turbulence.
Accurately predicting these different atmospheric phenomena and their underlying physical mechanisms is challenging. It is even more the case in the urban roughness sublayer due to the intricate patterns of cities, generally composed of heterogeneous and dense layouts of numerous buildings and trees, as well as various types of surfaces, heat and moisture sources. However, accurately predicting these flows in cities is of the utmost importance, especially as urban air flows are determining pollutant dispersion, pedestrian wind comfort, and structural and thermal loads on buildings.
Therefore, a key feature of predictive numerical models is the capability to handle realistic full scale configurations (including geometrical details and physical mechanisms at play) via high-fidelity unsteady simulation techniques well suited for turbulent flows. Large eddy simulations (LESs [2]) have appeared as a promising numerical approach for that purpose, whose main limitation is related to the simulation complexity and numerical cost [3,4]. To alleviate this problem, Lattice Boltzmann Methods (LBMs [5,6]) have recently been identified as one of the most efficient approaches for "revolutionnary computational fluid dynamics (CFD)" [7] since they allow for a drastic reduction in (i) the computational time compared to classical CFD approaches based on the Navier-Stokes equations and (ii) the preprocessing step including the volumic mesh generation thanks to the coupled use of embedded Cartesian grids and immersed boundary conditions techniques.
The efficiency of LBM-based solvers to handle full scale urban flow simulations has been demonstrated by several authors. As an example, they have been used recently to simulate the flow in a 19.2 × 4.8 × 1 km 3 area of Tokyo to evaluate the gust index at pedestrian level [8] and the turbulence statistics [9]. These simulations were carried out in a complex geometry thanks to the use of Cartesian grids and the immersed boundary conditions strategy that permit drastically reducing the preprocessing cost. Thanks to its high efficiency in terms of parallel computing, the LBM is an attractive method for implementation on Graphic Processing Units (GPUs), which reduce the computational cost. The use of the LBM for implementing GPUs allows real time simulations of the flow over a build area [10] and the dispersion of a pollutant inside Oklahoma City [11] to be achieved with quite good accuracy. The efficiency of the LBM compared to a Navier-Stokes (NS) approach has been analysed in the case of the cross flow through an open window in an isolated cubical building. The comparison of the velocity field inside the building shows a good accuracy of both NS-LES and LBM-LES simulations compared to experimental data and the LBM simulation on GPU was up to 700 times faster than the NS simulations [12,13].
The goal of the present paper is to illustrate and summarize an innovative CFD LES approach for the simulation of atmospheric and urban flows, developed in the framework of the ProLB software [14,15] and based on the LBM. Already used to study the aerodynamics of vehicles and airfoils [16][17][18], the approach was adapted to deal with atmospheric and urban physics problems. Thanks to its formulation and the treatment of boundary conditions, the approach is computationally effective. It is thus possible to perform efficient but detailed and high-fidelity simulations of complex built environments with multiscale interactions and multiphysical couplings, in order to study a large scope of atmospheric and urban physics problems.
To illustrate the ability of the developed LBM-LES approach to study mesoscale and urban microscale atmospheric flows, including transport and dispersion processes, the paper is organized as follows. First, Section 2 synthesizes the main numerical aspects of the developed method. Then, Sections 3-5 summarize recent developments and validated studies performed using the hybrid recursive regularized (HRR) collision model for the simulation of complex or/and coupled turbulent atmospheric flows, including a micrometeorological model, with different applications. In particular: • Section 3 reviews four applications of the developed method to the prediction of convective humid air flows with condensation with increasing meteorological likeness: a double diffusive Rayleigh Bénard convection (Section 3.1), rising two-and threedimensional moist thermal bubbles (Section 3.2) and a shallow cumulus convection (Section 3.3). • Section 4 reviews three use cases focusing on pollutant dispersion, considering different geometrical or physical complexity levels: neutral gas dispersion in a street canyon without or with tree planting effects (Section 4.1), neutral gas dispersion behind an isolated building without or with unstable thermal stratification (Section 4.2) and neutral or dense gas dispersion in a complex realistic urban environment (Section 4.3). • Lastly, Section 5 reviews the validation and study of velocities (Sections 5.1 and 5.3) and pressure distribution on building facades (Sections 5.2 and 5.4) in a complex urban environment including high-rise buildings, with uncertainty quantification towards a more relevant assessment of pedestrian wind comfort and wind loads on urban buildings.
To conclude, Section 6 discusses the ongoing developments, especially regarding indoor dispersion problems with thermal and moving people effects, and Section 7 closes the present paper by summarizing the main findings of the different studies performed, highlighting the benefit of the current approach to support the development of fast response models and decision making.

The HRR-LBM-WMLES Numerical Model for Atmospheric Flows
This section summarizes the main elements of the developed LB-based numerical model. More details about the general LBM can be found in Krüger et al.'s work [5].

Generalities about LBM
The LBM mainly aims at simulating the macroscopic behaviour of fluids. However, as compared to the macroscopic description of the flow underlying common Navier-Stokesbased approaches, the LBM is based on a mesoscopic description of the flow. The fluid dynamics is simulated through streaming and collision steps based on the lattice Boltzmann equation: where c i is a set of discrete velocities, usually a D3Q19 for three-dimensional problems (3 dimensions, 19 discrete velocities), f i (x, t) is the discrete distribution function, and Ω i (x, t) is the collision operator, i.e., the source term representing the redistribution of f i induced by collision. From this lattice Boltzmann equation, we can see that in practice only the first order neighbours are used in the algorithm, which permits increasing the parallel computation efficiency and allows the LBM to be faster than classical Navier-Stokes approaches. Multiscale expansions show that the three-dimensional weakly compressible Navier-Stokes equations can be recovered to the second order by the LBM.

Key Features of the Present Lattice Boltzmann Method
Based on the general LBM framework, different developments were made in the ProLB solver to deal more efficiently with atmospheric and urban problems. In particular, the collision term was estimated using a regularized BGK model with a hybrid recursive procedure [19] as follows: is the nonequilibrium part of the distribution function computed from the projection of distribution function on Hermite polynomials and f neq,FD i (x, t) is the nonequilibrium part of distribution functions approximated by finite differences. This procedure allows for both improved robustness as compared to usual BGK or MRT collision models and accuracy for large scale simulations, including urban and atmospheric flows. Different forcing mechanisms (e.g., buoyancy, Earth rotation, etc.) can also be taken into account in the model following Guo et al. [20]. In particular, buoyancy forces (F g ), due to differences in gas composition or temperature, can be modeled using a usual Boussinesq approach, as follows: with: where C is the concentration, C 0 is the reference concentration that is equal to 0, β c is the expansion coefficient of concentration, T is the temperature, T 0 is the reference temperature and β T is the expansion coefficient of temperature. Note that thermal buoyant flows can also be modeled thanks to the perfect gas law. Several other external forcings such as mesoscale or Coriolis effects can also be taken into account. The domain boundaries are integrated using the cut cell method, which enables complex geometries to be modeled while substantially reducing the preprocessing costs (see Feng et al.'s work [21] for details), and a synthetic eddy method (SEM, [22]) was developed to generate inlet turbulence [23].
Since the LBM exhibits the same turbulent closure problems as the Navier-Stokesbased methods, classical subgrid models and wall models (WMs) for atmospheric flow simulations (including neutral, convective and stable cases) are used in the present solver, leading to the definition of the present HRR-LBM-WMLES tool for atmospheric flow simulation. Details are omitted here for the sake of brevity and can be found in Feng et al.'s work [24].
The developed approach also relies on a hybrid strategy to address more complex problems than purely aerodynamic problems, such as dispersion or thermal problems. This approach does not consider multidistributions, but the mass and momentum conservation equations are solved using the LBM, while the scalar transport equations are solved using a usual finite volume/finite difference method for the sake of efficiency. The corresponding temporal and spatial discrete coordinates are the same as the ones used for the LBM. Thus, it is possible to consider only one additional unknown per additional equation and optimize the numerical efficiency of the model.

Application of the Present HRR-LBM-WMLES to Convective Humid Boundary Layers with Cloud Formation
This section reviews different use cases studied to finally address large scale buoyant meteorological flows accounting for cloud dynamics thanks to a condensation scheme. More details can be found in Feng et al.'s work [25,26].
For these applications, the potential temperature (θ) is defined as: where R q is the air gas constant per mass unit, C p is the averaged mass heat capacity, p 0 is the reference pressure at ground and p 0 (z) is the height-dependent reference state pressure. The air was assumed as a mixture of liquid water (mass fraction q l ), water vapour (mass fraction q v ) and dry air (mass fraction q d ), the rate of phase change was assumed to be infinitely fast, and the two phases were assumed to be in thermo-chemical equilibrium. The vapour and liquid water were modeled following Equations (8) and (9).
whereQ is a source term (typically the condensation) and D q is the water fraction diffusivity coefficient.

Double Convective Rayleigh-Bénard with Humid Air
The first validation case aims at validating the capability of the method to capture buoyancy effects, including both temperature and humidity gradients. It addresses the flow in a 2D square domain with temperature and humidity differences between its bottom and top boundaries as shown in Figure 1. Two Rayleigh numbers of 10 4 and 10 5 as well as two grids with δx = 0.02 m or half and δt = 0.0115 s or half were considered. Condensation was neglected and details of the model settings can be found in Feng et al.'s work [25]. HRR-LBM simulation results were compared to the reference solution of Ouertatani et al. [27], which used a finite volume method discretized using the QUICK scheme in the momentum equation, and a second order central difference scheme in the energy equation, with a nonuniform grid of 256 2 points.
The comparison between the simulated and reference velocity profiles along the domain midlines ( Figure 2), as well as local Nusselt number through the hot wall, potential temperature and total water humidity contours when the steady state had been achieved (not shown here), highlights a good performance of the developed HRR-LBM.

2D and 3D Rising Moist Bubbles
The second and third validation cases focus on a 2D and a 3D moist bubbles rising in a 2D or a 3D rectangular domain as shown in Figure 3. In both cases, the domain height was 2400 m and the domain width was 3600 m.  In the 2D configuration, a uniform grid with δx = 5 m and δt = 0.034 s was used, while in the 3D configuration a nonuniform grid with finest δx = 6.25 m and δt = 0.18 s was considered. Other model settings can be found in Feng et al.'s work [25]. 2D simulations were performed for 10 min physical time and the results were compared to reference data at 3, 5 and 7 min. 3D simulations were run for 6 min physical time and the results were compared to reference data at 2, 4 and 6 min.
HRR-LBM simulation results were compared to benchmark solutions, which used a multidimensional positive definitive advection transport algorithm and an anelastic solver with a 2.5 m grid in 2D configuration and on a grid with a 6.25 m finest grid spacing in 3D configuration [28,29].
Regarding the 2D bubble, HRR-LBM and reference results were compared in terms of highest vertical position of the 20% maximum q l contours and vertical fluid velocity at the top central position of the interface at different times ( Table 1). As in the reference study, the HRR-LBM results show the rising and expansion of the moist bubble as the phase changes. HRR-LBM results are thus generally in good agreement with the reference data except for the highest vertical position of the 20% maximum q l contours and the vertical fluid velocity at the top central position of the interface, especially at 7 min, which could be explained by the Boussinesq and not the aneleastic approximation used in the current study. Regarding the 3D case, the good performance of the HRR-LBM model was highlighted by the comparison of the vertical velocity profiles along the midline of the domain (Figure 4) as well as liquid humidity profiles and the liquid and vapour contours (not shown).

Atmospheric Cloud Formation
Finally, the fourth validation case addresses a convective atmospheric boundary layer case with a shallow cumulus formation to assess the capability of the HRR-LBM-WMLES model to predict the moist thermodynamics and its interactions with meteorological flows. The configuration refers to the model intercomparison case of the Barbados Oceanographic and meteorological Experiment (BOMEX). The considered domain is a 5 km long and 3 km high square-based domain with temperature and humidity surface fluxes and underwent an altitude-dependent geostrophic wind. Table 2 gives the related initial conditions. Table 2. Initial conditions set for the shallow cumulus convection case. For this configuration, HRR-LBM-WMLES simulations considered a δx = 40 m grid with δt = 0.27 s. The Monin Obukhov similarity theory was used as the surface model of the horizontal momentum components, temperature and humidity, and additional sources terms were added to the model to represent the large scale effects that cannot be included in LES. Details of the model settings can be found in Feng et al.'s work [24,25].
Simulations were run for 6 h physical time and statistics were computed on 1 h. The comparison of HRR-LBM and reference results of Siebesma et al. [30] in terms of profiles of mean velocity ( Figure 5), potential temperature, vapour and liquid water (not shown) highlights a good performance of the developed HRR-LBM-WMLES model, as the different atmospheric layers (mixed, conditionally unstable and inversion, from the ground upwards) and the instantaneous formation of the cloud are well captured. Hence, the hybrid LBM-based atmospheric numerical model, including the HRR collision model with forcing terms, and a finite volume method for temperature and water transport equation appear well suited to predict atmospheric humid convection with cloud formation, even considering large scale sources.

HRR-LBM-WMLES Application to Urban Pollutant Dispersion
This section reviews different case studies performed to assess the performance and the benefits of the proposed approach regarding the prediction of outdoor urban pollutant dispersion. Geometric complexity (trees, urban patterns) and buoyancy effects (pollutant or thermals) were especially studied.

Dispersion in a Street Canyon Including Tree Planting
A first validation study deals with dispersion of neutral traffic such as pollutant emissions in an ideal street canyon planted or not along a centred row of more or less dense tree crowns (λ = 0, 80 and 200 m −1 , respectively) and undergoing a perpendicular wind ( Figure 6). The drag force (F por ) induced by porous media is modeled via a Forchheimmer formulation, as follows: where R = λ/2 is the drag force coefficient and Φ is the ratio of porous media immersed in the volume cell.
Simulations were performed at reduced scale (1:150, H = 0.12 m) as in the reference experiment [31], with a domain of 3 × 2 × 1 m 3 and the street canyon located 0.84 m from the inlet. Table 3 gives the different boundary conditions. Accounting for δx = 0.00125 m, δt = 1.44 × 10 −5 s and five refinement levels, the total number of grid points was 41 × 10 6 . Simulations were run on 240 cores for 25 s physical time, with the last 10 s being kept for postprocessing. Details can be found in Merlier et al.'s work [32]. Quantitative (Table 4) and qualitative (not shown) comparisons of experimental and numerical mean concentrations on the leeward and windward walls of the street canyon highlight a state of the art performance of the HRR-LB model as compared to other reference studies.   [33]. Nonetheless, the results show a better agreement between numerical and experimental concentrations on the leeward wall, where the concentration is higher than on the windward wall, and for less dense tree crowns (λ = 0 and 80 m −1 ). The analysis of concentration distributions on walls also highlights higher concentrations on the streamwise symmetry plane as in the reference experiment.
In addition, thanks to high-fidelity unsteady simulations, the analysis of concentration statistics at different locations in the street canyon highlighted that high peaks of concentration can occur notably in the presence of dense tree crowns, which could be particularly harmful for short time exposure problems.

Dispersion behind a Building under Unstable Stratification
A second validation study deals with the dispersion of a gas released from a ground source just downwind of an isolated building located in an unstable boundary layer ( Figure 7). Simulations were performed at a reduced scale (H = 0.16 m) as in the reference study [34], with a domain of 2 × 1.2 × 1 m 3 according to AIJ guidelines [35], the building being located 0.32, 0.6 and 0.84 m from the inlet, lateral and top domain boundaries, respectively. Table 5 gives the corresponding boundary conditions. Accounting for δx = 0.002 m, δt = 12.8 × 10 −5 s and five refinement levels for the coarse grid case and δx = 0.001 m, δt = 6.4 × 10 −5 s and six refinement levels for the medium grid case, the total number of grid points was 4.7 × 10 6 for coarse grid and 11.2 × 10 6 for the medium grid. Simulations were carried out on 28 or 56 cores for 16.25 s physical time, the last 6.4 s being kept to compute the statistics. The pollutant source with a diameter of 0.005 m is located 0.04 m downstream the building. A gas flux of q = 9.17 × 10 −6 m 3 s −1 with a temperature of 30.4°C was imposed. The inflow velocity, temperature and turbulent kinetic energy were interpolated from experimental data provided in the TPU database [36] and are plotted in Figure 8.  The velocity, temperature and concentration of tracer gas were compared to experimental data at four different locations downstream of the building. Velocity was normalized with the reference velocity at building height (U * = U/U H ), the temperature was normalized using floor temperature (T f and the difference between floor temperature and temperature at building height (T * = (T − T f )/(T f − T H )) and the concentration was normalized by release gas flux, building height and reference velocity (c * = (cU H H 2 )/q). Figures 9 and 10 present the comparison of the HRR-LBM results obtained with a coarse and medium grid with the experimental data of Yoshie et al. [37]. A good agreement was obtained for concentration field downstream of the building and the present results are more accurate than other literature LES simulations related to the same case [34,38]. A fairly good agreement was found for streamwise velocity profiles downstream of the building, although velocity is generally underestimated compared to experimental data.

Dispersion of Neutral and Dense Gas in an Urban Area
A third validation study deals with the dispersion of neutral and dense gas from a ground source in a realistic urban environment. Two different source locations and wind incidences (C1: neutral and dense gas emitted in a rather channeled flow in a large avenue, C2: neutral gas emitted in a surround built crossroad, see Figure 11) are considered. Simulations were performed at reduced scale (1:350, H moy = 0.078 m) as in the reference experiment [39], with a domain of 8.75 (C1) or 9.5 (C2) ×3.5 × 1.5 m 3 , with the model located 1.15 m from the inlet of the domain. Table 6 gives the different boundary conditions. Accounting for δx min = 1.75 × 10 −3 m and δt = 1.5 × 10 −4 s with six refinement levels, the total number of grid points was 175 × 10 6 (C1) or 220 × 10 6 (C2). Simulations were run on about 10 3 cores for 27.5 s physical time and statistics were computed over 5 s. More details can be found in the work of Merlier et al. [40].
The performance indicators given in Table 7 for concentration at street level highlight a good performance of the model, especially for configuration 2. Results match the different acceptance criteria for urban dispersion models suggested in the work of Hanna and Chang [41], although the dense gas configuration shows a worse, but still acceptable, accuracy. Indeed, the analysis of the spatial distribution of concentration (not shown) highlights that: • The neutral gas emitted in the channeled flow is carried mainly in the main street by the prevailing flow and progressively mixes upwards; • The dense gas emitted in the channeled flow follows a rather circular and flat dispersion all around the source; • The dispersion of the neutral gas emitted in the crossroad is notably vertical due to the presence of the obstacles built upstream and downstream.
The comparison of the HRR-LBM and reference vertical as well as horizontal concentration profiles above the canopy (not shown) are generally also satisfactory regarding both the mean and the standard deviation of concentration levels. These results suggest that the dynamics of dispersion processes are well reproduced by the model. Being capable of providing the temporal statistics of velocity and concentration, the developed approach appears well suited to support the design of fast response models.

HRR-LBM-WMLES-Based Urban Wind Prediction with Application to Pedestrian Wind Comfort and Building Wind Loads under Uncertainty
A next step toward an improved reliability of numerical simulations for realistic full scale urban applications is the capability to account for uncertainties that appear in the prescription of the atmospheric conditions, instantaneous wind conditions, surface roughness, pollutant source features, etc. The variability of the numerical solution with respect to uncertain parameters must be quantified in order to provide users with useful information, since a single fully deterministic solution is meaningless in such problems. To this end, efficient techniques for Uncertainty Quantification (UQ) must be used. A challenging issue is that most existing UQ techniques require the use of a significant number of samples, each sample being currently an unsteady high-fidelity LES of the case under consideration. Therefore, adequate UQ methods must be defined that minimize the number of required simulations while preserving the accuracy of the uncertainty propagation in the results. Such a method, the c-APK method, was developed by Margheri and Sagaut [42] with application to urban flow simulations in complex areas.
The c-APK method will not be detailed here for the sake of brevity, and the reader is referred to Margheri and Sagaut [42] for details. Nonetheless, to highlight two additional applications of the c-APK method, this section reviews two validation studies and two projected applications related to the prediction of urban air flows with UQ: pedestrian wind comfort and wind loads on buildings. The configuration studied for the projected applications is a new hypothetical tower located in a complex and dense urban area.

Prediction of Mean Flow Field in a Complex Urban Environment
The first validation study focuses on a full scale urban area of Tokyo, which includes an area of low-rise buildings upstream a cluster of towers. Simulations were performed at full scale, considering H max = 225 m, to evaluate the performance of the model at a realistic Reynolds number, as both full scale and reduced scale data are provided in the reference experimental data [43]. More precisely, the configuration studied corresponds to the case F of the open source Architectural Institute of Japan's database, which gathers full scale and wind tunnel measurements, wind tunnel tests corresponding to the 1977 full scale measurements campaign. Figure 12 details the different model dimensions and boundary conditions used for the simulations. Accounting for δx = 0.5 m and δt = 0.0075 s with seven refinement levels, the total number of grid points was 136 × 10 6 for the finest grid (respectively, δx = 1 m, δt = 0.015 s, six refinement levels and 54 × 10 6 grid points for medium grid). Simulations were run on 504 cores (240 for medium grid) for 2 h physical time, the last 1 h being kept for postprocessing. More details can be found in the work of Jacob and Sagaut [44].
The HRR-LBM and reference mean velocities reported in Figure 13 at 10 m high and for different locations in the area of interest evidence a good performance of the developed model and also for mean velocity.

Prediction of Surface Pressure on a High-Rise Building
To evaluate the accuracy of the developed approach also on building surface quantities before coupling the HRR-LBM-WMLES tool with a UQ technique, a second validation study was extensively carried out at a reduced scale isolated high-rise building (1:300, H = 0.49 m), for which detailed measurement data of velocity and pressure on the model surface are available [45]. Figure 14 details the different model dimensions and boundary conditions, which included an extension of the original incompressible SEM to reconstruct the inlet turbulence in the current LBM framework, given the importance of the turbulent inflow properties dynamically studying wind loading on isolated structures. With δx = 1.56 × 10 −3 m and δt = 1.5 × 10 −5 s and six refinement levels, the total number of grid points was 5.8 × 10 6 . More details can be found in the work of Buffa et al. [23].
The comparison of mean and standard deviation values of pressure coefficients at the building surface for the different faces of the building and at different heights ( Figure 15) shows a very satisfactory match for all the tested data, especially for the mean Cp value. Experimental and simulated data were also extensively compared in terms of spectral analysis and local instantaneous pressure maxima (not shown), showing the reliability of the developed approach for wind load prediction. Figure 15. Reference and simulated pressure coefficients for the isolated high-rise building case.

Wind Comfort Assessment with Uncertainty Quantification
Thanks to the possibilities offered by the developed dynamic model, pedestrian wind comfort at street level (H = 2 m) was studied considering Melbourne criteria [46] in an area of 400 × 400 m 2 inside Shinjuku area as shown in red in Figure 16. The domain used in this study is the same as in Section 5.1. Simulations were run on 1.5 h physical time and the statistics were computed over the last hour. As described by Jacob and Sagaut [44], two buildings were added in Shinjuku area: one in the middle of the area (in blue in Figure 16) for wind comfort assessment and one (in green in Figure 16) for the study of Section 5.4. The case presented by Jacob and Sagaut [44] with the wind coming from the north is here considered as the reference sample for the UQ analysis, and several other samples have been considered by changing the velocity magnitude by a factor α (0.6 ≤ α ≤ 1.4) and wind direction with an angle θ around the north direction (−30 ≤ θ ≤ 30°). A set of 57 simulations have been performed for this study permitting to estimate first the sensitivity of the mean velocity field to α and θ and to compute the mean value of mean velocity at pedestrian level using the c-APK method. Figure 17 shows the Sobol indices obtained for the mean velocity field at the pedestrian level, which permits quantification the influence of each parameter on the global variance of the system. The results highlight that the inflow velocity magnitude is generally the main influential parameter in most locations in the area of interest. In some locations, wind direction has a large contribution on the variance, especially in areas that can be located in wake for some inflow directions.  Figure 18 shows the classification obtained following the Melbourne criteria using the data of the reference sample (α = 1, θ = 0) and the mean data obtained from the 57 samples. The results highlight that by using the c-APK output, a large part of the area is located in Zones A, B and C where pedestrian comfort can be considered as good, whereas a part of the area is unclassified (here the pedestrian comfort is not so good) using only the reference sample.

Prediction of Structural Wind Loads with Uncertainty Quantification
The same configuration as in Section 5.3 was considered to study pressure distribution on the facades of the hypothetical building (green building in Figure 16) for further applications to the prediction of structural wind loads. The analysis is based on the same set of simulations changing the inflow velocity magnitude and direction.
The Sobol indices plotted on the different building faces in Figure 19 highlight that the inflow velocity magnitude is the most influential parameter on the pressure distribution on the building faces, except on a part of the north face, where the wind direction appears to influence pressure values more than the wind velocity magnitude. This is explained by the fact that a large part of the north face is located in a wake area when the flow is coming from the northeast direction, whereas it is directly exposed to inflow wind when it comes from the northwest direction. The second order term is very low compared to the others and does not significantly contribute to the global variance of the pressure on building faces. Figure 20 presents the pressure coefficients obtained in the case of the reference sample (α = 1, θ = 0) on the different building faces and the pressure coefficients computed from the average pressure given by the c-APK analysis. Few differences were observed on top, east, west and south faces, whereas the estimated pressure coefficient was lowered by the c-APK analysis on the north face compared to the reference sample. This result suggests that uncertainties on the inflow wind should be accounted for in urban simulations since it has an impact on the dynamics of the flow around high-rise buildings and on the related surface quantities.

Indoor Pollutant Dispersion with Thermal and Moving Body Effects
Another step toward the numerical simulation of realistic full scale configuration is the use of high-fidelity CFD for evacuation problems. In such cases, human agents present in the domain under consideration may have an effect on the pollutant dispersion because they trigger some physical mechanisms (mixing induced by the wakes of moving persons, natural convection due to the heat release by persons, breathing effects, etc.), leading to the definition of a two-way coupling problem. In such a problem, the behaviour of human agents has a direct influence on the flow, but it is governed by their responses to external parameters, which result from the conjunction of physical (e.g., heat) and psychological parameters. Therefore, a more complex level of modeling is required that couples classical physical CFD to psychological and behavioural models.
The present HRR-LBM-WMLES model has been coupled with a Social Force Model (SFM) [47][48][49][50] that allows us to evaluate individual human displacement, taking into account some psychological effects, and therefore to account for the effect of crowd evacuation on pollutant dispersion. The trajectories of people leaving a room are evaluated using SFM, then the drag force of each person is added into LBM equation using the actuator line method (ALM) [51] to account for the effect of human motion on the fluid, leading to the occurrence of wakes.
An application of this coupled method is shown here. It is related to an evacuation scenario inside a concert hall [52]. This hall (with dimensions 83.2 × 68.4 × 15.1 m 3 ) contains 12 exits, indicated in green in Figure 21, and is initially occupied by 6026 persons. The concert hall is equipped with ventilation systems located on the room ceiling that allow for the balancing of the heat release of the persons in the concert hall. We consider an instantaneous pollutant release from a cloud with a 20 m diameter in the centre of the concert hall, as shown in Figure 21. Evacuation starts at the same time as pollutant release and 299.2 s are necessary to allow all the people inside the concert hall to exit.
The pollutant field at human head level obtained at several times is displayed in Figure 22. We can see here that the pollutant is advected inside persons wake from the initial cloud to the two central upper exits. Pollutant dispersion towards the other exits is not significant since people leaving the concert all through these exits were not initially located inside the cloud.

Concluding Remarks
The versatility and efficiency of the HRR-LBM-WMLES approach for atmospheric flow simulations, urban physics and pollutant transport prediction have been illustrated by a broad range of applications including humidity effects with phase changes and evacuation prediction with coupling to a social force model.
The main advantages of the present approach are its efficiency in terms of computational time, which is due to the explicit nature of the lattice Boltzmann method, the compactness of the underlying stencil, and the preprocessing time which is drastically reduced thanks to the use of embedded uniform grids along with an immersed boundary approach to handle complex fully arbitrary geometries.
In all cases, a very satisfactory agreement with reference data (if they exist) is reported, demonstrating the accuracy of the simulation tool. The coupling with the c-APK method for uncertainty quantification was also illustrated, showing that the HRR-LBM-WMLES tool is fast enough to allow for the use of UQ tools in complex full scale configurations.
Hence, thanks to its reliability-highlighted through the different validation studiesrelevance-highlighted by the different physical analyses carried out-and its computational efficiency-induced by its numerical properties-the developed method appears very to be promising to support the design of fast response models and urban decision making. Acknowledgments: Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high-performance computing resources. This work was granted access to the HPC resources of TGCC/CINES under the allocation 2021-A0092A07679 made by GENCI.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: