# CFD simulation of pollutant dispersion using anisotropic models: Application to an urban like environment under neutral and stable atmospheric conditions

Boulos Alam, Rayan Nkenfack SOPPI, Amir-Ali Feiz, Pierre Ngae, Amer Chpoun, Pramod Kumar

<figure><img src="image_1.png" /></figure>

► To cite this version:

Boulos Alam, Rayan Nkenfack Soppi, Amir-Ali Feiz, Pierre Ngae, Amer Chpoun, et al.. CFD simulation of pollutant dispersion using anisotropic models: Application to an urban like environment under neutral and stable atmospheric conditions. Atmospheric Environment, 2024, 318, pp.120263. 10.1016/j.atmosenv.2023.120263. hal-04321151

HAL Id: hal-04321151

https://hal.science/hal-04321151v1

Submitted on 15 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
CFD simulation of pollutant dispersion using
anisotropic models: Application to an urban like
environment under neutral and stable atmospheric
conditions

Boulos Alam<sup>a,*</sup>, Rayan Nkenfack SOPPI<sup>a,b</sup>, Amir-Ali Feiz<sup>a</sup>, Pierre Ngae<sup>a</sup>, Amer
Chpoun<sup>a</sup>, Pramod Kumar<sup>c</sup>

Abstract

Modeling the atmospheric dispersion of pollutants emitted from different
types of sources under various atmospheric conditions is an essential prereq-
uisite for risk assessment studies and emergency preparedness. In this study,
we evaluate the performance of different turbulence models and velocity-scalar
correlation models implemented in the *Code_Saturne*. The evaluation is done
with observations of four trial cases of the Mock Urban Setting Test (MUST)
campaign in an urban type environment and in neutral and stable atmospheric
conditions. For all of the trials studied, the CFD model with a first-order clo-
sure model (k - ε) predicts 61.1% of the concentrations within a factor of two
of the observations, which is higher than the percentage of predicted points
(58.8%) when a second-order closure model (R<sub>ij</sub> - ε) is used. Overall, the
CFD model underestimates observed concentrations, regardless of the turbu-
lence model used. For the trial with slightly stable conditions, the results show
that the k - ε model combined with an algebraic SGDH model predicts 75%
concentrations within a factor of two of the observations. The performance of
the k - ε model is compared to that of the R<sub>ij</sub> - ε model when used with the alge-
braic SGDH, GGDH models and with the scalar flux transport equation (DFM
model). Under slightly stable atmospheric conditions, the DFM model predicts
69% of concentrations within a factor of two of observations, showing promise
for modeling under these conditions, despite its relatively high computational
cost.

Keywords: CFD modeling, MUST experiment, algebraic models, anisotropic,
stable conditions

<footnote><sup>a</sup>Université Paris-Saclay, Univ Evry, LMEE, 91020, Evry, France</footnote>

<footnote><sup>b</sup>ONERA-Université Paris-Saclay, Châtillon 92322, France</footnote>

<footnote><sup>c</sup>Université Paris-Saclay, CNRS, CEA, UVSQ, Laboratoire des sciences du climat et de l'environnement, 91191, Gif-sur-Yvette, France</footnote>

<footnote>*Corresponding author.</footnote>

<footnote>Email address: boulos.alam@univ-evry.fr (Boulos Alam)</footnote>
1

# 1. Introduction

2

3

4

5

6

7

8

9

10

11

12

The dispersion of gaseous pollutants and small particles in an urban area is characterized by complex turbulence phenomena, such as flow detachment and high shear zones (Yu and Thé, 2017). In addition to this dynamic turbulence, thermal air turbulence tends to dampen the turbulent fluctuations resulting from flow shear under conditions of stable or very stable stratification, thus reducing pollutant dispersion. The complexity of the physical phenomena involved in these conditions results in high anisotropy of the atmospheric boundary layer (ABL), which limits the ability of standard isotropic models to correctly predict the dispersion field. Computational fluid dynamics (CFD) models are one of the most promising tools to describe the dispersion of pollutants in the near field of obstacles.

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

The CFD approaches used to solve the turbulent flow equations in most studies mainly involve large eddy simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) methods. Large eddy simulations (LES) are considered a promising step in atmospheric flows, as they provide a better understanding of the various flow phenomena that occur around obstacles (e.g. (Tominaga and Stathopoulos, 2012, 2011, 2010; Salim et al., 2011)). However, this type of simulation requires higher computing power due to the time dependence of the simulations and the finer grids that have to be used (Lateb et al., 2016). Reynolds-Averaged Navier-Stokes (RANS) simulations use Reynolds decomposition to filter turbulent fluctuations from the mean flow field (Kavian Nezhad et al., 2022; Bahlali et al., 2019; Tominaga and Stathopoulos, 2018, 2017; Yu and Thé, 2017; Miliez and Carissimo, 2008, 2007; Kumar and Feiz, 2016; Kumar et al., 2015). The turbulence on the mean flow field is fully incorporated by a turbulence model and has an effect on both the dynamic and scalar fields. As a result, the accuracy and reliability of field modeling results are highly dependent on the quality of the turbulence model applied. Researchers are constantly investing in RANS simulations to model atmospheric fluxes, trying to improve the submodels included. In most simulations of atmospheric dispersion in urban areas, the standard $k-\epsilon$ model (Launder and Spalding, 1974) is the most widely used. The standard model is usually modified or calibrated to adapt to the problem of dispersion in different terrain topographies and atmospheric stabilities (see for example (Kumar et al., 2015; Miliez and Carissimo, 2007; Alinot and Masson, 2005), etc.). In addition to the standard $k-\epsilon$ model, the Reynolds Stress Model (RSM) has been used in numerous test cases for pollutant dispersion in cities and urban areas. This model, originally developed for high Reynolds number flows, has sometimes proved to outperform the standard $k-\epsilon$ isotropic model. For instance, Knaus et al. (2018) compared three turbulence models - standard $k-\epsilon$, $k-\epsilon$ renormalization group (RNG) and RSM model- to simulate wind flow over complex terrain. Only the RSM model proved capable of predicting turbulence intensity with a high degree of accuracy. Bahlali et al. (2019) showed that the use of the $R_{ij}-\epsilon$ turbulence model developed by Launder et al. (1975) gives a more physical representation of the obstacles' influence on the concentration profiles.
Regarding the turbulent flux in the scalar transport equation, it is often solved using the Simple Gradient Diffusion Hypothesis (SGDH). There are, however, other models, such as the Generalized Gradient Diffusion Hypothesis (GGDH) (Daly and Harlow, 1970), the Differential Flow Model (DFM), etc., more commonly used in an industrial context. For atmospheric dispersion problems, Rossi and Iaccarino (2009) showed that gradient-transport models are improved by the use of algebraic closures (GGDH and HOGGDH) and predict atmospheric dispersion in complex flows better than the standard eddy-diffusivity model (see also (Rossi et al., 2009)).

The objective of this work is to perform a comprehensive validation study of different anisotropic turbulence models for flow and scalar fields. Simulation results are evaluated using observations from the Mock Urban Setting Test (MUST) campaign (Yee and Biltoft, 2004; Biltoft, 2001), which is a quasi-urban experiment involving the release of a neutral gas from a point source inside a set of rectangular shipping containers. Four trials were chosen to simulate flow and dispersion patterns near the source and at greater distances. The stability conditions in these trials range from neutral to extremely stable. Numerical simulations are performed using the three-dimensional (3-D) code *Code_Saturne*. The $k-\epsilon$ (Launder and Spalding, 1974) and $R_{ij}-\epsilon$ LRR (Launder et al., 1975) models are used and compared to predict the flow field, while the SGDH and GGDH models are initially adopted to relate the scalar flux to the mean gradients of the concentration field. An additional analysis of closure for turbulent scalar flux based on the differential flow model (DFM) transport equation is also presented, and the results are compared with the predictions provided by the SGDH and GGDH algebraic models. The capabilities and limitations of the tested models for both wind flow and concentration fields are rigorously evaluated using multiple validation metrics.

The paper is organized as follows: section 2 describes the experimental data sets and tests selected for the numerical simulations. Section 3 presents the CFD code, including the solved equations and turbulence models used in this study. Section 4 presents the numerical setup for the simulations in the code. Section 5 presents the numerical results of the simulations and the qualitative and quantitative analysis of the predictions against the experimental observations. Section 6 summarizes the conclusions of this study.

## 2. Description of the experimental data sets

### 2.1. *The MUST experiment*

The experimental campaign MUST (Mock Urban Setting Test) (Yee and Biltoft, 2004; Biltoft, 2001) was conducted at the US Army Dugway Proving Ground (DPG), in the Utah desert, in the western United States (40°12.606′N–113°10.635′W) from September 6 to 27, 2001, by the Defence Threat Reduction Agency. The main objective of this campaign was to obtain an experimental database on the dispersion of a pollutant at the scale of a pseudo-city in order
88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

• Line 1: receptors 1 through 12,

111

• Line 2: receptors 13 to 21,

112

• Line 3: receptors 22 to 30,

113

• Line 4: receptors 31 to 40.

114

115

116

to develop new models or validate existing numerical models. Figure 1 shows a photograph of the experiment site.

To highlight the presence of obstacles in the urban landscape (buildings, trees, etc.), the buildings were rectangular containers of dimensions 12.2m × 2.42m × 2.54m. The urban canopy is thus formed by 120 containers uniformly distributed in 10 rows of 12 containers. The terrain of the experimental site was mostly flat, with an average roughness length of 0.045±0.0005m (Yee and Biltoft, 2004).

Receptors are placed between rows of containers on four masts, a central tower, and four sampling lines to collect information on the dynamic field and concentration field of the pollutant, as well as outside the canopy to measure the velocity field and temperature. Figure 2 shows the location of the containers, pollutant source, and measurement devices.

The experiment consisted of 63 continuous releases of a non-reactive gas into the air under different meteorological and atmospheric stratification conditions (stable, unstable and neutral). The gas released is propylene ($C_3H_6$) and each release lasts for 15min, from which quasi-stationary periods of 200 seconds were extracted and used to evaluate the simulated dispersion results in each trial of the experiment (Yee and Biltoft, 2004). The concentration field is mainly measured using 40 PID (Photoionization Detectors) placed at 1.6m from the ground and arranged horizontally along 4 rows. According to Figure 2, the sampling lines for the concentration measurements are defined as follows:

Concentrations were also measured at 8 levels (1, 2, 4, 6, 8, 10, 12, and 16m) vertically on the 32m tower above the ground surface, and at 6 levels on masts A, B, C, and D.

117

## 2.2. Selected trials

118

119

120

121

122

123

124

125

For the present work, 4 trials are selected for the validation of the dispersion model. The characteristics of each trial, used as input data for the model simulations, are summarized in Table 1. Average wind speed and direction are taken at a height of 4m above the ground surface from the south tower upstream of the urban canopy (Yee and Biltoft, 2004). Observations of surface friction velocity ($u*$), Obukhov length ($L$), and turbulent kinetic energy ($k$) are taken at a height of 4m from a 32m tower (tower T) in the center of the domain. The selected trials are the following:
<figure><img src="image_3.png" /><figcaption>Figure 1: Photograph of the MUST experimental campaign site (Biltoft, 2001).</figcaption></figure>

* Trial No.2681829: this trial was performed under neutral atmospheric conditions, i.e. with zero vertical potential temperature gradient for a very high value of the Obukhov length (L = 28000m). The ABL in this trial has a stability class D according to the Pasquill classification (Pasquill, 1972). The average wind direction is -41° and the wind is oriented from southwest to northeast (Figure 2). Of the twenty trials in the experimental campaign, which are often used to validate models, this trial recorded the strongest wind.

* Trial No.2692157: this trial was conducted under weakly stable atmospheric conditions, with a moderate wind velocity (2.1 < U < 5m/s). It is of stability class F according to the classification of (Pasquill, 1972) (100 < L < 200m). In this trial, the wind blows from southeast to northwest with a direction angle of 43° (Figure 2).

<table><caption>Table 1: Physical characteristics of the selected trials, obtained at a reference height of 4m: U is the mean wind speed, α is the incident wind direction measured from the x-axis, k is the turbulent kinetic energy, Q is the pollutant flow rate, z_s is the height of the release source, u_* is the frictional velocity, L is the Obukhov length, and ε is the dissipation rate.</caption><thead><tr><th>Trial No.</th><th>U (m s<sup>-1</sup>)</th><th>α (°)</th><th>k (m<sup>2</sup> s<sup>-2</sup>)</th><th>Q (L min<sup>-1</sup>)</th><th>z<sub>s</sub> (m)</th><th>u* (m s<sup>-1</sup>)</th><th>L (m)</th><th>ε (m<sup>2</sup> s<sup>-3</sup>)</th></tr></thead><tbody><tr><td>2671934</td><td>1.63</td><td>-48</td><td>0.148</td><td>200</td><td>1.8</td><td>0.08</td><td>5.8</td><td>0.0013</td></tr><tr><td>2672033</td><td>2.69</td><td>-26</td><td>0.251</td><td>200</td><td>1.8</td><td>0.17</td><td>4.8</td><td>0.012</td></tr><tr><td>2681829</td><td>7.93</td><td>-41</td><td>1.46</td><td>225</td><td>1.8</td><td>1.1</td><td>28000</td><td>0.8</td></tr><tr><td>2692157</td><td>2.98</td><td>43</td><td>0.505</td><td>225</td><td>2.7</td><td>0.39</td><td>130</td><td>0.042</td></tr></tbody></table>

* Trial No.2672033: this trial was conducted under highly stable atmospheric conditions, with moderate winds (2.1 < U < 5m/s). It is of stability class F according to the classification of (Pasquill, 1972) (0 < L < 10m). Among all the trials of the experimental campaign, the lowest value of Obukhov length was obtained for this trial. The wind blows from southwest to northeast with a direction angle of 26° (Figure 2).
<figure><img src="image_4.png" /><figcaption>Figure 2: Schematic representation of the MUST experimental site: the triangles represent the masts A, B, C and D, the blue square represents the central tower of 32m height. The position of the sources for the 4 selected trials is represented by yellow stars.</figcaption></figure>

• Trial No.2671934: this trial was carried out under very stable atmospheric conditions, with low winds (U < 2.1m/s). It is of stability class F according to the classification of (Pasquill, 1972) (0 < L < 10m). The lowest velocity value was observed for this trial. The wind blows from southwest to northeast with a direction angle of 48° (Figure 2).

### 3. Description of Code_Saturne model for flow and dispersion

#### 3.1. Governing equations

Code_Saturne (https://www.code-saturne.org/cms/web/documentation/v6) is a 3-D CFD parallelized code which solves the Navier-Stokes equations by a stationary or non-stationnary algorithm for incompressible or weakly dilatable fluids, isothermal or not, with consideration of passive scalars (Archambeau et al., 2004). The spatial discretization is obtained using a finite volume approach for a collocated arrangement of all variables. In this study, the atmospheric module of the code is used along with an Eulerian RANS approach to solve the average form of the equations for both turbulent flow and scalars. The equations of conservation of mass (continuity), momentum, potential temperature and species concentration are as follows:

• Continuity equation (with the anelastic approximation):

$$ \frac{\partial \bar{\rho} \overline{U_i}}{\partial x_i} = 0 \qquad (1) $$

• Conservation of momentum:

$$ \frac{\partial \overline{U_i}}{\partial t} + \overline{U_j} \frac{\partial \overline{U_i}}{\partial x_j} = -\frac{1}{\bar{\rho}} \frac{\partial \bar{P}}{\partial x_i} + \frac{\partial}{\partial x_j} \nu \left( \frac{\partial \overline{U_i}}{\partial x_j} + \frac{\partial \overline{U_j}}{\partial x_i} \right) - \frac{(\bar{\rho} - \rho_{ref})}{\bar{\rho}} g \delta_{i3} - \frac{\partial}{\partial x_j} \overline{U_i' U_j'} \qquad (2) $$
• Conservation of potential temperature:

$$ \frac{\partial \bar{\Theta}}{\partial t} + \overline{U_j} \frac{\partial \bar{\Theta}}{\partial x_j} = \frac{1}{\bar{\rho} C_p} \frac{\partial}{\partial x_j} \left( \lambda \frac{\partial \bar{\Theta}}{\partial x_j} \right) - \frac{\partial}{\partial x_j} \overline{\Theta' U'_j} + \overline{S_{\Theta}} \quad (3) $$

• Conservation of concentration (Advection-diffusion equation for a scalar):

$$ \frac{\partial \overline{C}}{\partial t} + \overline{U_j} \frac{\partial \overline{C}}{\partial x_j} = \frac{\partial}{\partial x_j} \left( D \frac{\partial \overline{C}}{\partial x_j} \right) - \frac{\partial}{\partial x_j} \overline{C' U'_j} + \overline{S_C} + \overline{R_C} \quad (4) $$

where, $\rho$ is the air density, $U$ is the velocity vector, $p$ is the pressure, $\Theta$ is the potential temperature, $\rho_{ref}$ the air density at hydrostatic reference state. $g$ is the gravitational acceleration vector, $C_p$ is the specific heat at constant pressure, $\lambda$ is the effective thermal conductivity, $D$ is the molecular diffusion coefficient for the species, $S_{\Theta}$ is a thermal source term, $S_C$ is the species source term, $R_C$ is the term of production/destruction of the species by chemical reaction. The Reynolds tensor $\overline{U'_i U'_j}$, thermal turbulent flux $\overline{\Theta' U'_j}$ and concentration turbulent flux $\overline{C' U'_j}$ are modeled using turbulence closure models.

## 3.2. Turbulence models

### 3.2.1. The $k-\epsilon$ model

This model is based on the Boussinesq hypothesis which assumes a alignment between the average strain rate tensor and the anisotropic part of the Reynolds tensor (Launder and Spalding, 1974). Its formulation is given by the following relation:

$$ \overline{U'_i U'_j} = \nu_t \left( \frac{\partial \overline{U_i}}{\partial x_j} + \frac{\partial \overline{U_j}}{\partial x_i} \right) - \frac{2}{3} k \delta_{ij} \quad (5) $$

The transport equations for the kinetic energy ($k$) and the dissipation rate ($\epsilon$) are written as:

$$ \frac{\partial k}{\partial t} + U_i \frac{\partial k}{\partial x_i} = \frac{1}{\bar{\rho}} \frac{\partial}{\partial x_j} \left[ \left( \nu + \frac{\nu_t}{\sigma_k} \right) \frac{\partial k}{\partial x_j} \right] + P_t + G - \epsilon \quad (6) $$

$$ \frac{\partial \epsilon}{\partial t} + U_i \frac{\partial \epsilon}{\partial x_i} = \frac{1}{\bar{\rho}} \frac{\partial}{\partial x_j} \left[ \left( \nu + \frac{\nu_t}{\sigma_\epsilon} \right) \frac{\partial \epsilon}{\partial x_j} \right] + C_{\epsilon1} \frac{\epsilon}{k} (P_t + C_{\epsilon3} G) - C_{\epsilon2} \frac{\epsilon^2}{k} \quad (7) $$

With:

$$ P_t = -\bar{\rho} \overline{U'_i U'_j} \frac{\partial \overline{U_j}}{\partial x_i}, \quad G = -\frac{1}{\bar{\rho}} \frac{\beta g}{\text{Pr}_t} \nu_t \frac{\partial \bar{\theta}}{\partial z}, \quad \nu_t = \bar{\rho} C_\mu \frac{k^2}{\epsilon} \quad (8) $$

where $P_t$ and $G$ are respectively the mechanical and buoyancy production rate of $k$, $\nu_t$ is the turbulent eddy viscosity, $\beta \sim \frac{1}{\Theta}$ is the thermal expansion coefficient, $\mu_l$ is the kinematic viscosity, $\sigma_k$, $\sigma_\epsilon$ and $\text{Pr}_t$ are the turbulent Prandtl numbers. $C_{\epsilon1}$, $C_{\epsilon2}$, $C_{\epsilon3}$ and $C_\mu$ are empirical constants. The constant $C_{\epsilon3}$ is set to zero for unstable stratification and to 1 for stable stratification.
It is known that the constants of the standard model may not be suitable to reproduce the turbulence of the atmospheric flow as they were originally calibrated for industrial flows (Sadek et al., 2013). Therefore, a new model with a set of calibrated constants, proposed by Duynkerke (1988), is used in our calculations and compared to the model with the standard parameterization. The model constants are given in Table 2. We denote by [LA83] the standard $k - \epsilon$ model of Launder and Spalding (1974) while by [DU88] the calibrated model of Duynkerke (1988).

<table><caption>Table 2: $k - \epsilon$ model constants.</caption><thead><tr><th></th><th>C<sub>ε1</sub></th><th>C<sub>ε2</sub></th><th>σ<sub>k</sub></th><th>σ<sub>ε</sub></th><th>C<sub>μ</sub></th></tr></thead><tbody><tr><td>$k - \epsilon$ [LA83]</td><td>1.44</td><td>1.92</td><td>1.0</td><td>1.3</td><td>0.09</td></tr><tr><td>$k - \epsilon$ [DU88]</td><td>1.44</td><td>1.92</td><td>1.0</td><td>2.38</td><td>0.033</td></tr></tbody></table>

In fact, the vertical gradient of potential temperature $\partial\bar{\theta}/\partial z$ is negative ($G \ge 0$) under unstable atmospheric conditions, therefore, the thermal stratification has an effect of production of the turbulence. Under stable conditions, however, $\partial\bar{\theta}/\partial z$ is positive ($G \le 0$), therefore, the thermal stratification has an effect of destruction of the turbulence.

### 3.2.2. The $R_{ij} - \epsilon$ model

This model solves six transport equations for the Reynolds tensors $R_{ij} = \overline{U_i'U_j'}$ and one equation for the dissipation rate $\epsilon$. The latter is very similar to that of the $k - \epsilon$ model, only the turbulent viscosity and gravity terms change (Launder et al., 1975). For each $R_{ij}$ component, the following transport equation is solved:

$$ \frac{\partial R_{ij}}{\partial t} + U_k \frac{\partial R_{ij}}{\partial x_k} = P_{ij} + D_{T,ij} + D_{M,ij} + G_{ij} + \phi_{ij} - \epsilon_{ij} \quad (9) $$

where $P_{ij}$ is the average shear production term, $D_{T,ij}$ is the turbulent diffusion term, $D_{M,ij}$ is the molecular diffusion term, $G_{ij}$ is the production/destruction term by thermal effects, $\phi_{ij}$ is the pressure-strain term and $\epsilon_{ij}$ is the dissipation term.

<table><caption>Table 3: $R_{ij} - \epsilon$ LRR model constants.</caption><thead><tr><th></th><th>σ<sub>k</sub></th><th>σ<sub>ε</sub></th><th>C<sub>ε1</sub></th><th>C<sub>ε2</sub></th><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>C'<sub>1</sub></th><th>C'<sub>2</sub></th><th>C<sub>L</sub></th><th>C<sub>μ</sub></th></tr></thead><tbody><tr><td>$R_{ij} - \epsilon$ [LA75]</td><td>0.82</td><td>1.0</td><td>1.44</td><td>1.92</td><td>1.8</td><td>0.6</td><td>0.5</td><td>0.3</td><td>0.39</td><td>0.09</td></tr><tr><td>$R_{ij} - \epsilon$ [SA13]</td><td>1.0</td><td>2.01</td><td>1.44</td><td>1.92</td><td>1.8</td><td>0.6</td><td>0.94</td><td>0.03</td><td>0.19</td><td>0.033</td></tr></tbody></table>

We note that the terms $G_{ij}$, $\epsilon_{ij}$, $D_{T,ij}$ and $\phi_{ij}$ are not closed and therefore must be modeled. The modeling of the $G_{ij}$ and $\epsilon_{ij}$ terms is identical to that of the $k - \epsilon$ model (Launder et al., 1975).
Similar to the $k-\epsilon$ model, most of the standard empirical constants in the $R_{ij}-\epsilon$ model have been calibrated to achieve the atmospheric turbulence levels of the experiment. The methodology for setting the values of the constants is presented by Sadek et al. (2013), who also draws on the work of Panofsky and Dutton (1984); Duynkerke (1988); Rodi (1993). The new calibrated $R_{ij}-\epsilon$ model is therefore used in the present work along with the standard form. The coefficients of the $R_{ij}-\epsilon$ model for the standard [LA75] and calibrated [SA13] versions are presented in Table 3.

### 3.2.3. Models for scalar turbulent fluxes

By averaging the transport equations for temperature and concentration, we obtain turbulent heat flux and turbulent concentration flux. These terms must now be modeled either by a transport equation or by algebraic models (GGDH and SGDH). The latter two are obtained from the initial transport equation assuming that the turbulence is in equilibrium, such that: $P_k+G_k = \epsilon$ and $P_{ss} = \epsilon_{ss}$, where $P_k = \frac{1}{2}P_{jj}$ corresponds to half the trace of the production tensor, while $G_k = \frac{1}{2}G_{jj}$ refers to half the trace of the production tensor resulting from volumetric forces. In addition, $P_{ss}$ represents the production term of the variance $S'^2$ of scalar $S$, and $\epsilon_{ss}$ denotes dissipation of the variance of the same scalar quantity. For simplicity, the expression for each model is written only for the turbulent concentration flux responsible for the dispersion of the passive tracer, knowing that those for the turbulent thermal flux can be deduced in a similar manner.

* Differential Flux Model (DFM)

This approach is the most suitable to represent the physical mechanisms leading to turbulent fluxes (Lecocq, 2008; Lecocq et al., 2008; Dehoux, 2012). The transport equation for the scalar flux is written as:

$$ \frac{\overline{DC'U'_j}}{Dt} = -\overline{C'U'_i}\frac{\partial\overline{U_j}}{\partial x_i} - \overline{U'_iU'_j}\frac{\partial\overline{C}}{\partial x_i} - \beta g_i \overline{C'^2} - \frac{P'}{\rho}\frac{\partial\overline{C'}}{\partial x_j} - (D + \nu)\frac{\partial\overline{C'}}{\partial x_i}\frac{\partial\overline{U'_j}}{\partial x_i} \\ + \frac{\partial}{\partial x_i}\left(-\overline{C'U'_iU'_j}\right) + \frac{\partial}{\partial x_i}\left(D\frac{\partial\overline{C'}}{\partial x_i}U'_j + \nu\overline{C'}\frac{\partial\overline{U'_j}}{\partial x_i}\right) + \frac{\partial}{\partial x_i}\left(\frac{\overline{C'P'}}{\rho}\delta_{ij}\right) \quad (10) $$

* Generalized Gradient Diffusion Hypothesis (GGDH)

The GGDH (Daly and Harlow, 1970) model is obtained from the algebraic flux model by removing the production terms due to velocity and buoyancy gradients. Thus, the equation of turbulent scalar flux (neglecting the molecular diffusion) is written as:

$$ \overline{C'U'_j} = -C_{\theta} \frac{k}{\epsilon} \left( \overline{U'_kU'_j} \frac{\partial \overline{C}}{\partial x_k} \right) \quad (11) $$

with $C_{\theta} = 0.3$ (Daly and Harlow, 1970; Dehoux, 2012).
• Simple Gradient Diffusion Hypothesis (SGDH)

The SGDH model considers that the turbulent flux of scalar is proportional to the gradient of scalar, according to Fick's law (the molecular diffusion is neglected). It is given by:

$$ \overline{C'U_j'} = -\frac{\nu_t}{Sc_t} \frac{\partial \overline{C}}{\partial x_j} \qquad (12) $$

where $Sc_t$ is the turbulent Schmidt number ($Sc_t = 0.7$).

# 4. Numerical setup for simulations

## 4.1. Computational domain

To model the experimental site, we first created the geometry using Salome software (www.salome-platform.org), with the goal of reproducing the experimental field as closely as possible as well as the containers present. Figure 3 shows an illustration of the computational domain. As the container grid occupies a footprint of 193m x 171m, we used a domain with a surface area of 240m x 240m in order to comply with the best practice rules proposed in the COST Action 732 guidebook (Franke et al., 2007; Schatzmann et al., 2010) for urban air dispersion modeling. In particular, the report of Franke et al. (2007) mentions that a minimum distance of 5Hmax to 15Hmax (where Hmax is the height of the tallest building in the domain, here equal to 2.54m) between the domain boundaries and the area of interest should be considered, depending on the type of boundary face (inlet, outlet, top or side face). Regarding the domain height, other studies, such as the one conducted by Reiminger (2020), have proposed a minimum height of 96m to be used in the absence of obstacles, so this minimum distance is also needed for full-scale modeling in the presence of obstacles. The domain height also depends on the thermal stratification of the atmosphere and the position of the pollutant source. In this work, we will only deal with dispersion cases in neutral, stable and very stable atmosphere with a source located at a height lower than two meters from the ground. Given the position of the source (quite close to the ground), and since in these atmospheric conditions the vertical dispersion of the pollutant is not as strong as in unstable atmosphere, a domain height of 50m is sufficient to reproduce completely the pollutant plume. The size of the computational domain is therefore 240m x 240m x 50m.

## 4.2. Mesh characteristics

In order to enable the numerical solution of the Navier-Stokes and concentration transport equations on the created domain, it is necessary to discretize the computational domain. To do this, the NETGEN-2D algorithm of the salome software is used to generate an unstructured mesh in the ground plane, with increasingly finer cells near the containers and increasingly coarser cells as one moves away from the container grid. Then, this two-dimensional (2-D) planar
<figure><img src="image_5.png" /><figcaption>Figure 3: Representation of the computational domain with boundary conditions.</figcaption></figure>

mesh is extruded in the vertical direction with a growth factor so that a mesh height of 0.25m near the ground and 5.5m at the top of the domain is achieved. The total number of cells in the domain is 1,210,582. Figure 4 shows a 3-D representation of the mesh of the computational domain. The mesh is chosen after a sensitivity study performed on three mesh sizes for trial No.2692157. Table 4 summarizes the details of each of the three grid sizes. The results showing the sensitivity of the dynamic variables to the mesh size, obtained for three different meshes are presented in figures 5 and 6. Overall, the results show no difference, except for the first few meters near the ground. Therefore, the reference mesh (Mesh2) is chosen for the study.

<figure><img src="image_6.png" /><figcaption>Figure 4: Mesh of the computational domain: (a) External view (b) Internal view.</figcaption></figure>
<table><caption>Table 4: Details of the mesh grids used for the sensitivity study. *h* is the height of the first mesh adjacent to the ground wall.</caption><thead><tr><td>Mesh</td><td>Total number of cells</td><td><i>h</i> (m)</td></tr></thead><tbody><tr><td>Mesh1</td><td>1 095 052</td><td>0.15</td></tr><tr><td>Mesh2</td><td>1 210 582</td><td>0.25</td></tr><tr><td>Mesh3</td><td>1 095 052</td><td>0.35</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 5: Sensitivity of dynamic variable profiles to mesh size at mast A.</figcaption></figure>

### 4.3. Boundary conditions and pollutant release condition

In principle, the ABL extends over the entire surface of the earth. However, for practical reasons of computational cost and considering the local character of the experimental field, we have reduced the extent of the computational domain to an area of a few hundred square meters. This computational domain is thus an integral part of the ABL and consequently, the conditions applied to its borders must ensure continuity with the surrounding ABL. To do this, the types of boundary conditions imposed on the faces of the domain are as follows:
<figure><img src="image_8.png" /><figcaption>Figure 6: Sensitivity of dynamic variable profiles to mesh size at mast B.</figcaption></figure>

• At the inlet of the domain: a Dirichlet-type condition is applied, based on velocity, temperature and turbulence profiles whose expressions are developed in the next section.

• At the outlet of the domain : a Neumann type condition is applied, except for the pressure (free outflow).

• On the ground and the containers : a logarithmic rough wall function with aerodynamic roughness of 0.04m is applied.

• On the side faces: automatic inlet or outlet condition is applied depending on the wind direction of the flow.

• At the top of the domain: an inlet condition (Dirichlet) is imposed, allowing a maximum wind speed to be maintained at the center of the highest cell in the domain. A condition of symmetry is inadequate in this case since the fluxes of heat and non-zero momentum created at ground level will tend towards a zero value at the top. Consequently, this behavior
is in disagreement with the physics of the surface boundary layer (SBL) which requires a constancy of the flow with the altitude. Figure 3 shows the type of boundary conditions on each face of the domain.

### 4.3.1. Input profiles and turbulence models

For this work, a series of simulations is carried out on selected trials using the atmospheric module of Code_Saturne. The profiles of wind speed, temperature and turbulence variables are specified at the input boundary. Table 5 represents the dynamic input profiles as well as the models used for each of the selected simulations.

The inlet wind velocity is well-described by a logarithmic vertical profile, derived from the Monin-Obukhov similarity theory (Stull, 1988). It is used to determine the mean velocity in the surface layer. It is given by the following equation:

$$U(z) = \frac{u*}{\kappa} \left[ \ln\left(\frac{z}{z_0}\right) - \psi_m(\zeta) + \psi_m(\zeta_0) \right] \quad (13)$$

Where $\kappa$ (0.41) is the von-Kármán constant (Dyer, 1974); $\zeta$ (z/L) is the Monin-Obukhov stability parameter; $\psi_m$ is the momentum stability correction function. In neutral stability regime, $\psi_m = 0$; in stable regime, $\psi_m(\zeta) = -5\zeta$ (Hicks, 1976; Businger et al., 1971; Dyer and Hicks, 1970).

In all simulations, the Monin-Obukhov vertical temperature profile is set as the inflow thermal boundary condition, such as:

$$T(z) = \Theta(z) + z\lambda_{adia} \quad (14)$$

With:

$$\Theta(z) = \Theta(z_r) + \frac{\Theta^*}{\kappa} \left[ \sigma_n \ln\left(\frac{z}{z_r}\right) - \psi_h(\zeta) + \psi_h(\zeta_0) \right]$$

Where $\lambda_{adia}$ (-0.009766°C/m) is the adiabatic lapse rate; $z_r$ is the reference (anemometer) level at which $\Theta$ is determined; $\Theta^*$ is the temperature scale; $\psi_h$ is the stability correction function for potential temperature. In neutral stability regime, $\Theta^* = 0$; in stable regime, $\sigma_n = 0.74$ and $\psi_h = -5\zeta$ (Hicks, 1976; Businger et al., 1971; Dyer and Hicks, 1970).

For the trials conducted under stable conditions, the following wind velocity profiles are also used:

* Profile of Gryning et al. (2007) :

$$\frac{U(z)}{u_*} = \frac{1}{\kappa} \left( \ln\left(\frac{z}{z_0}\right) + \frac{5z}{L} \left(1 - \frac{z}{2z_i}\right) + \frac{z}{L_{MBL}} - \frac{z}{z_i} \left(\frac{z}{2L_{MBL}}\right) \right) \quad (15)$$

With:
$$ \frac{u_*}{fL_{MBL}} = \left( -2 \ln \left( \frac{u_*}{fz_0} \right) + 55 \right) \exp \left( -\frac{(u_*/fL)^2}{400} \right) $$

• Profile of Peña et al. (2009) (see also (Konow, 2014)):

$$ U(z) = \frac{u_*}{\kappa} \left[ \ln \left( \frac{z}{z_0} \right) + b \frac{z}{L} \left( 1 - \frac{z}{2z_i} \right) + \frac{1}{d} \left( \frac{\kappa z}{\eta} \right)^d - \left( \frac{1}{1+d} \right) \frac{z}{z_i} \left( \frac{\kappa z}{\eta} \right)^d - \frac{z}{z_i} \right] \quad (16) $$

With

$$ \eta = \frac{\kappa z_i}{[d(1+d)^{1/d}]} \left[ \left( \left[ \ln \left( \frac{u_*}{fz_0} \right) - A \right]^2 + B^2 \right)^{1/2} + \ln \left( \frac{z_i}{z_0} \right) \right]^{-1/d} $$

$$ A = 1.5, B = 5.2 \text{ et } d = 1.25 $$

• Profile of Beljaars and Holtslag (1991):

$$ U(z) = \frac{u_*}{\kappa} \left[ \ln \left( \frac{z}{z_0} \right) - \Psi_m \left( \frac{z}{L} \right) \right] \quad (17) $$

With

$$ \Psi_m = -a \left( \frac{z - z_0}{L} \right) - b \left( \frac{z}{L} - \frac{c}{d} \right) e^{-d \frac{z}{L}} + b \left( \frac{z_0}{L} - \frac{c}{d} \right) e^{-d \frac{z_0}{L}} $$

$$ a = 1, b = 0.667, c = 5 \text{ et } d = 0.35 $$

Where $L_{MBL}$ is the length scale in the middle layer of the ABL, $f_c = 10^{-4}$ s$^{-1}$ is the Coriolis parameter, $z_i \approx 0.1u_*/f$ is the height of the ABL (Zilitinkevich and Baklanov, 2002).

For the turbulent kinetic energy and its dissipation rate, the following profiles are used:

• Profiles of Gorlé et al. (2009):

The evolution of $k(z)$ with altitude $z$ is modeled by a square root mathematical relationship whose parameters $C_1$ and $C_2$ are determined by a least squares method (Gorlé et al., 2009). The experimental measurements to which the $k(z)$ profile is fitted are those taken by three 2-D horizontal sonic anemometers at levels 4, 8 and 16m from a 16m telescopic pneumatic
mast located 30.5m upstream of the obstacle array. The input profile of turbulent dissipation is then deduced from the smoothing of $k(z)$. The values of $C_1$ and $C_2$, obtained for each trial, are given in Table ?? . The expressions of $k$ and $\varepsilon$ are as follows:

$$k(z) = \sqrt{C_1 \ln(z + z_0) + C_2} \qquad (18)$$

$$\varepsilon(z) = \frac{\sqrt{C_{\mu}} u_{*}}{\kappa (z + z_{0})} \sqrt{C_{1} \ln(z + z_{0}) + C_{2}} \qquad (19)$$

• Profiles of Kerschgens et al. (2000)

These profiles are used to model the turbulence in the ABL under stable atmospheric conditions. Under these conditions, their analytical expressions are given by the following relations:

$$k(z) = \frac{1}{2} (\overline{U'^2} + \overline{V'^2} + \overline{W'^2}) = \frac{1}{2} (\sigma_U^2 + \sigma_V^2 + \sigma_W^2) \qquad (20)$$

$$\varepsilon(z) = \frac{u_{*}^{3}}{\kappa z} \left(1 + 4 \frac{z}{L}\right) \qquad (21)$$

With :

$$\sigma_U = 2.4u_* \exp \left( -0.3 \frac{z}{z_i} \right),$$

$$\sigma_V = 2.0u_* \exp \left( -0.3 \frac{z}{z_i} \right),$$

$$\sigma_W = 1.3u_* \exp \left( -0.3 \frac{z}{z_i} \right)$$

As for the turbulence models used for the simulations, the $k - \varepsilon$ model is applied with the SGDH approach. In addition, the $R_{ij} - \varepsilon$ model is applied with the SGDH and GGDH approaches as well as with the DFM model (Table 5).

In order to evaluate the performance of different algebraic models (SGDH and GGDH), as well as that of the transport equation (DFM model), we select trial No.2692157 to perform the simulations. The results obtained with this trial, for modeling scalar concentration flux, are developed and analyzed in section 5.3 of this article.

For each turbulence model used, numerous simulations are performed by changing the parameterization of the dynamic profiles (wind speed and turbulent
<table><caption>Table 5: Input profiles and turbulence models used for all simulations on selected trials of the MUST experiment.</caption><thead><tr><th>Trial No.</th><th>Stability condition</th><th>Velocity profiles</th><th>Temperature profiles</th><th>Turbulence profiles</th><th>Turbulence models</th><th>Scalar fluxes turbulence models</th></tr></thead><tbody><tr><td rowspan="2">2671934</td><td rowspan="2">Very stable</td><td>Eqn (15)</td><td rowspan="2">Eqn (14)</td><td>Eqn (20 &amp; 21)</td><td>k-ε [DU88]</td><td rowspan="2">SGDH</td></tr><tr><td>Eqn (17)</td><td>Eqn (18 &amp; 19)</td><td>R<sub>ij</sub>-ε [SA13]</td></tr><tr><td rowspan="2">2672033</td><td rowspan="2">Highly stable</td><td>Eqn (15)</td><td rowspan="2">Eqn (14)</td><td>Eqn (20 &amp; 21)</td><td>k-ε [DU88]</td><td rowspan="2">SGDH</td></tr><tr><td>Eqn (17)</td><td>Eqn (18 &amp; 19)</td><td>R<sub>ij</sub>-ε [SA13]</td></tr><tr><td rowspan="2">2681829</td><td rowspan="2">Neutral</td><td rowspan="2">Eqn (13)</td><td rowspan="2">Eqn (14)</td><td>Eqn (18 &amp; 19)</td><td>k-ε [DU88]</td><td rowspan="2">SGDH</td></tr><tr><td></td><td>R<sub>ij</sub>-ε [LA75]</td></tr><tr><td rowspan="2">2692157</td><td rowspan="2">Weakly stable</td><td rowspan="2">Eqn (13)</td><td rowspan="2">Eqn (14)</td><td>Eqn (18 &amp; 19)</td><td>k-ε [DU88]</td><td rowspan="2">SGDH</td></tr><tr><td></td><td>R<sub>ij</sub>-ε [LA75]</td></tr><tr><td rowspan="2">2692157</td><td rowspan="2">Weakly stable</td><td rowspan="2">Eqn (16)</td><td rowspan="2">Eqn (14)</td><td>Eqn (20 &amp; 21)</td><td>R<sub>ij</sub>-ε [SA13]</td><td rowspan="2">GGDH</td></tr><tr><td></td><td>R<sub>ij</sub>-ε [SA13]</td></tr></tbody></table>

variables) at the domain entrance each time. For instance, the profile of Businger et al. (1971) is applied with the $k-ε$ and $R_{ij}-ε$ models (both with SGDH). The profiles of Na et al. (2009) and Gryning et al. (2007) are applied with the model $R_{ij} - ε$ (with GGDH and DFM, respectively). As for the turbulent variables, they are computed from the profiles of Gorlé et al. (2009) for the $k-ε$ and $R_{ij} - ε$ models (both with SGDH), and from the profiles of Kerschgens et al. (2000) for the $R_{ij} - ε$ (with GGDH and DFM). In what follows, the analysis of the simulation results will be based only on the turbulence models used in each simulation.

### Notes:

* Regardless of the type of turbulence model (isotropic or anisotropic) used in the simulations, the meteorological turbulence input profiles are those that describe the evolution of kinetic energy, isotropically, with altitude. In other words, in the case of the $R_{ij}-ε$ model, only the diagonal terms of the $R_{ij}$ tensors are considered to construct the TKE, and the other terms are considered to be zero. We assume that this is a reasonable assumption for operational applications.

* For trials involving stable atmospheric conditions, thermal turbulent flux is always modeled using the SGDH hypothesis, while concentration turbulent flux is modeled using different algebraic models (SGDH or GGDH) or
405

406

407

408

the concentration scalar transport equation (DFM model). Consequently, only the effect of the model from the point of view of dispersion is studied, focusing on the passive scalar, and the study excludes the influence of thermal turbulent flux.

409

### 4.3.2. Pollutant release condition

410

411

412

Regarding the release condition, the pollutant is injected through a scalar source term. In the injection cell, the pollutant flow rate in L/min (Table 1) is imposed for the studied trials.

413

### 4.4. CFD solver options

414

415

416

417

418

419

420

421

422

423

424

The simulations are performed in stationary mode in order to solve the wind and dispersion fields in the domain. Thus, two calculations have been implemented. The first calculation is used to determine the dynamic field, i.e. only the wind field, and stops when the steady state is reached. The second calculation is used to compute the dispersion of the passive scalar, in this pre-calculated wind field with frozen velocity, pressure and temperature. To ensure convergence of the results, the accuracy of the solver is set to $10^{-4}$ for all the computed variables. The spatial numerical scheme for each variable (convection scheme) is set automatically by the solver. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm is used for the pressure-velocity coupling with a reduction of the relaxation coefficient.

425

## 5. Results and discussion

426

427

### 5.1. Comparison between $k-\epsilon$ and $R_{ij}-\epsilon$ models: dynamic field in neutral condition

428

429

430

431

The neutral case is chosen to study the effect of the anisotropic turbulence model $R_{ij}-\epsilon$ both on the atmospheric flow and on the dispersion of the passive tracer. For this purpose, the scalar turbulent flux is modeled only by the SGDH model.

432

433

434

435

436

437

438

439

Figure 7 illustrates the velocity and TKE fields obtained using both the $k-\epsilon$ and $R_{ij}-\epsilon$ models at $z = 4m$. Figure 8 shows the lateral velocity and turbulent kinetic energy fields, obtained with the $k-\epsilon$ and $R_{ij}-\epsilon$ models at $x = 130m$. We find that the presence of the containers significantly influences the flow, even at heights greater than the height of a container ($H = 2.54$ m). The container grid acts as a rough surface, causing a deceleration of the flow downstream and amplifying turbulence generation through the shearing of the mean flow and the wake interference produced behind each container.

440

441

442

By comparing the two turbulence models, we find that the $R_{ij}-\epsilon$ model generates relatively low levels of turbulence in the near field of the container grid as well as in the far field (Figure 8b). Analysis of the numerical results
for the dynamic field revealed the recirculation phenomena in the flow, which both models ($k-\epsilon$ and $R_{ij}-\epsilon$) effectively capture at the rear of the containers. However, the length of the recirculation as well as the size of the recirculation bulbs are larger with the $k-\epsilon$ model (result not shown here).

<figure><img src="image_9.png" /><figcaption>Figure 7: Horizontal wind speed (left) and TKE (right) fields at $z = 4m$ for the neutral atmospheric stability trial No. 2681829 using the $k-\epsilon$ model (a, b) and the $R_{ij}-\epsilon$ model (c, d).</figcaption></figure>

Figures 9 and 10 show the velocity and TKE profiles, respectively, on masts A and D, and tower T. At masts A and D, the velocity and TKE profiles are essentially identical between the $k-\epsilon$ and $R_{ij}-\epsilon$ models and in agreement with the available experimental measurement points. The peak production of the TKE, a few meters above the canopy, is consistent with the TKE contours observed in Figure 8b, for both turbulence models. However, the small number of measurement points (two measures) does not guarantee a good correlation between the numerical and experimental profiles over the whole height of the
<figure><img src="image_10.png" /><figcaption>Figure 8: Lateral dynamic fields at $x = 130m$ (a) wind speed and (b) TKE for the neutral atmospheric stability trial No.2681829 with $k - \varepsilon$ (left) and $R_{ij} - \varepsilon$ (right) models.</figcaption></figure>

computational domain. On the central tower for which we have four measurement points, we observe an overestimation of the velocity by the numerical models and an underestimation of the TKE compared to the experimental data.

## 5.2. Comparison between $k - \varepsilon$ and $R_{ij} - \varepsilon$ models: dispersion field in neutral condition

The pollutant release in this trial is from a 1.8m high source with a flow rate of 225L/min. Figure 11 represents the pollutant concentration field at ground level for the $k - \varepsilon$ and $R_{ij} - \varepsilon$ models.

From a qualitative standpoint, we have noted a minor impact of the turbulence model on pollutant dispersion. The plume described by the $k - \varepsilon$ model appears to be slightly larger in comparison to the one predicted by the $R_{ij} - \varepsilon$ model. This distinction is primarily due to the slight overestimation of TKE by the $k - \varepsilon$ model, as depicted in Figure 10. Since the turbulent fluctuations are larger, the turbulent concentration flux is also larger. Therefore, we observe a more extensive plume with the $k - \varepsilon$ model. We can also note that the plume guideline does not coincide with the wind direction (-41° with respect to the x-axis). This is due to the presence of a channeling effect imposed by the containers on the plume, making it deviate from the wind direction. This plume deflection phenomenon has also been observed by Milliez and Carissimo (2007); Kumar et al. (2015); Bahlali et al. (2019) etc., who explained it by the channeling effect that the containers impose on the plume, causing it to deflect from the wind direction. In other words, the aspect ratio of the containers also influences plume dispersion.

Figure 12 shows the concentration profile over all 48 receptors, with the
<figure><img src="image_11.png" /><figcaption>Figure 9: Comparison between the observed velocity and the one calculated with the $k-\epsilon$ and $R_{ij}-\epsilon$ models for the neutral case No.2681829 at the position of (a) Mast A (b) Mast B (c) Mast D and (d) Tower T.</figcaption></figure>

last eight (41-48) placed vertically on the central tower T. For both models, we observe a main concentration peak measured at receptor No.8, followed by four secondary peaks recorded by the receptors placed in the direction of the plume. As the plume moves away from the source, it is diffused and its concentration decreases. This is the reason why the concentration at the secondary peaks evolves in a decreasing way. This phenomenon is also observed in the experiment, but with a main peak located at receptor No.9, followed by four secondary peaks. The correspondence between the numerical and experimental results is globally satisfactory for the position of the peaks. It can be noted that the underestimation of the concentration increases with distance from the source.
<figure><img src="image_12.png" /><figcaption>Figure 10: Comparison between the observed TKE and the one calculated with the $k-\epsilon$ and $R_{ij}-\epsilon$ models for the neutral case No.2681829 at the position of (a) Mast A (b) Mast B (c) Mast D and (d) Tower T.</figcaption></figure>

## 5.3. Evaluation of anisotropic models for the stable case

### 5.3.1. Dynamic field

Figure 13 shows the velocity profiles computed at different masts using the CFD code and compared to experimental observations. We note that the velocity profiles at the masts vary slightly when the velocity and turbulence input profiles are varied. However, the modification of the turbulence models leads to more noticeable variations of the velocity profiles, in particular in the first five meters above the ground, where the $R_{ij}-\epsilon$ model leads to a variation of the sign of the vertical velocity gradient. Since this model is anisotropic, it is supposed to better capture the physical phenomena present in the near wall (ground and containers). However, due to the small number of experimental measurements between 0 and 6 meters from the ground, it is quite difficult to identify if this change in sign of the gradient captured by the $R_{ij}-\epsilon$ model reflects the reality of the experiment.
<figure><img src="image_13.png" /><figcaption>Figure 11: Average pollutant concentration field (kg/kg) at ground level (z = 1.6m), with the models (a) $k - \varepsilon$ (b) $R_{ij} - \varepsilon$ (With SGDH) for Trial No.2681829.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: Comparison between the observed concentration profile and the one calculated using the $k - \varepsilon$ and $R_{ij} - \varepsilon$ models across all PID sensors for Trial No.2681829.</figcaption></figure>

Figure 14 shows the TKE profiles calculated at different masts using the CFD code, in comparison with the experimental measurements. A first observation can be made on the shape of the profiles, which is more or less the same, from a qualitative point of view. The profiles show the presence of three large zones in the flow:

* From 0m to about 5m above the ground (almost 2 times the height of the containers): the vertical magnitude of the TKE is positive. Turbulence
is generated by the flow's shearing action against the ground and the containers' walls, making it a zone of turbulence production.

* From 5m to 20m: the vertical gradient of the TKE is negative, the turbulent kinetic energy decreases gradually under the action of the viscosity of the air. It is thus a zone of dissipation of the TKE.

* From 20m to 50m: the vertical gradient of TKE is very weak (almost null in certain cases). The production and dissipation of turbulence compensate each other, thus the turbulence is in equilibrium.

Moreover, we notice that the $k-\epsilon$ model tends to predict an overproduction of the turbulence. For this model, the turbulence peak on all masts is the highest as well as the vertical gradient of TKE in the production area (steeper slope than that with the $R_{ij}-\epsilon$ model). This result highlights the difficulty of the $k-\epsilon$ model in predicting the anisotropy of the turbulence. The $R_{ij}-\epsilon$ model, being anisotropic, takes better into account the presence of the ground and the containers and allows to obtain a profile of the TKE which fits better the experimental measurements.

### 5.3.2. Dispersion field

Figure 15 shows the predicted mean concentration field at ground level with the $k-\epsilon$ and $R_{ij}-\epsilon$ models. Note that the plume size is significantly larger when employing the $k-\epsilon$ model, except in the case of the $R_{ij}-\epsilon$ model combined with DFM. This can be explained by the fact that the $k-\epsilon$ model predicts an overproduction of turbulence (Figure 14) close to the ground which increases the turbulent transport flux and thus leads to a larger lateral diffusion of the pollutant over the extent of the grid. As in the case of neutral stratification, we observe that the directional line of the plume does not coincide with the wind direction (43° from x-axis).

Figure 16 shows the comparison between observed and predicted concentrations at the receptors using different turbulence models. We observe a main peak followed by four secondary peaks. We also find that the concentration intensity at the main peak (7ppmv) is about twice as high as the main peak observed in trial No.2681829 with a neutral stability condition (3.5ppmv). This increase in concentration is primarily due to the atmospheric conditions under which the trial is conducted. Since trial No.2692157 is conducted under moderately stable atmospheric conditions, turbulence generation is attenuated by thermal stratification of the air, which reduces turbulent transport of the pollutant; therefore, concentrations are locally higher. We also observe, in this case, a good match between the numerical results and the measurements concerning the position of the peaks.

As for the magnitude of the highest peak, all the used models overestimate it, except the $R_{ij}-\epsilon$ model (with GGDH) which is the only one to predict the
<figure><img src="image_15.png" /><figcaption>Figure 13: Comparison between the observed velocity and the one calculated with the $k-\varepsilon$ and $R_{ij}-\varepsilon$ models for the stable case No.2692157 at the position (a) Mast A (b) Mast B (c) Mast C and (d) Tower T.</figcaption></figure>

maximum concentration measured by the third receptor. On the other hand, moving away from the source, this same model seems to underestimate the concentration the most among the other models. In general, the predictions obtained by all the models are close to the experimental measurements when moving away from the source. However, no model is able to correctly reproduce the profile corresponding to the concentrations observed by receptors 41 to 48 on the central tower.

The statistical assessment involves the measured concentrations from 40 receptors distributed along four sampling lines, as well as eight additional receptors located on central tower T, along with the corresponding average concentrations predicted by the simulations. The performance measures calculated in this study are the normalized mean square errors (NMSE), the fractional bias (FB), the correlation coefficient (COR), mean geometric bias (MG), the geometric variance (VG), and the fractional accuracy factor 2 (FAC2), which is the
<figure><img src="image_16.png" /><figcaption>Figure 14: Comparison between the observed TKE and the one calculated with the $k-\varepsilon$ and $R_{ij}-\varepsilon$ models for the stable case No.2692157 at the position of (a) Mast A (b) Mast B (c) Mast C and (d) Tower T.</figcaption></figure>

fraction of predictions within a factor of two of field measurements. The results of the statistical study, to assess the accuracy of each model, are presented in Table 6. A model is considered perfect if: NMSE = FB = 0 and COR = FAC2 = MG = VG = 1 (Chang and Hanna, 2004). A model is considered sufficiently efficient if it satisfies the following criteria: -0.3 < FB < 0.3, 0.7 < MG < 1.3, NMSE < 4, VG < 1.6 and FAC2 ≥ 0.5 (Chang and Hanna, 2004). We recall that these criteria are applicable to dispersion within urban areas.

All FB values are positive, indicating that all models tend on average to underestimate the concentration, with the $k-\varepsilon$ model being the least underestimating. Examination of the MG values in Table 6 shows that no model succeeds in satisfying the acceptance criterion, with the $k-\varepsilon$ model being the closest to the tolerable range. The same is true for the VG values. As for the NMSE values, they are all within the accepted range. Therefore, all models perform well for the statistical measures (FB and NMSE), but appear to have
a systematic bias at low concentrations based on the values for the logarithmic measures (MG and VG). All models succeed in having more than 50% of the predicted concentrations within a factor of two (FAC2) of the observations. The $k-\varepsilon$ model is the best performing with 75% of the predictions within a factor of two of the observations. As for the $R_{ij} - \varepsilon$ model, it performs significantly worse with about 63%, 60% and 69% with the SGDH, GGDH and DFM models, respectively.

Finally, good one-to-one correlations can be observed for all models through the high values of the correlation coefficient, the highest of which ($\text{COR} = 0.92$) is for the $R_{ij} - \varepsilon$ model (with GGDH and DFM). Comparing the $k-\varepsilon$ and $R_{ij} - \varepsilon$ models, we find that the latter, known to be anisotropic, does not necessarily improve the predictions over the former, regardless of the scalar turbulent flux model with which it is associated. However, a study by Bahlali et al. (2019) shows that using the $R_{ij} - \varepsilon$ model provides a more physical representation of the influence of obstacles on concentration profiles. The three velocity-scalar correlation models (SGDH, GGDH, and DFM), when employed alongside the $R_{ij} - \varepsilon$ model, achieve a score of three (indicating that three statistical measures meet their criteria) in the statistical evaluation, and therefore, all three models can be considered for retention. The GGDH hypothesis, which uses the Reynolds stress tensor does not improve the predictions compared to the isotropic SGDH model (lower FAC2 and more under-prediction). In fact, the GGDH model is obtained by considering assumptions of weak equilibrium ($\frac{DU_j'\theta'}{Dt} = 0$) and dynamic and thermal equilibrium of turbulence, compared to the DFM model. In addition, the production terms due to the velocity gradient $\overline{U_j'\theta'}\frac{\partial \overline{U}_i}{\partial x}$ and buoyancy $g_i\Theta'^2$ (Refer to Equation 10) are suppressed in the formulation of Daly and Harlow (1970). This model is not suitable for flows with natural convection and buoyancy, and does not give satisfactory results in the case of very strong stratification (Dehoux, 2012). By analogy, the same applies to atmospheric flows under stable stratification.

The DFM model, on the other hand, is more suitable for stratified atmospheric flows and gives better results than the GGDH model (less under-prediction and higher FAC2). This is probably due to the presence of a gravity term with scalar variance in the scalar turbulence transport equation.

<table><caption>Table 6: Statistical performance measures evaluated for trial No.2692157.</caption><thead><tr><th></th><th>FB</th><th>MG</th><th>NMSE</th><th>VG</th><th>FAC2(%)</th><th>COR</th></tr></thead><tbody><tr><td>$k-\varepsilon + \text{SGDH}$</td><td>0.13</td><td>1.39</td><td>0.37</td><td>2.23</td><td>75</td><td>0.91</td></tr><tr><td>$R_{ij} - \varepsilon + \text{SGDH}$</td><td>0.27</td><td>1.66</td><td>0.45</td><td>3.3</td><td>62.5</td><td>0.90</td></tr><tr><td>$R_{ij} - \varepsilon + \text{GGDH}$</td><td>0.5</td><td>1.81</td><td>0.59</td><td>2.98</td><td>60.4</td><td>0.92</td></tr><tr><td>$R_{ij} - \varepsilon + \text{DFM}$</td><td>0.37</td><td>1.86</td><td>0.45</td><td>4.47</td><td>68.8</td><td>0.92</td></tr><tr><td>Fluidyn-Panache</td><td>0.2</td><td>1.28</td><td>0.34</td><td>1.5</td><td>87.5</td><td>0.91</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>Figure 15: Average pollutant concentration field (kg/kg) at z = 1.6m, with models (a) $k - \varepsilon$ (with SGDH) (b) $R_{ij} - \varepsilon$ (with SGDH) (c) $R_{ij} - \varepsilon$ (with GGDH) and (d) $R_{ij} - \varepsilon$ (with DFM) for trial No.2692157.</figcaption></figure>

5.4. Overall model performance

5.4.1. Quantitative evaluation of Code_Saturne models

In this section, we present a statistical evaluation study of the model for the selected set of simulated trials. Figure 17 shows the scatter plots of observed and predicted concentrations using the $k - \varepsilon$ and $R_{ij} - \varepsilon$ models. Table 7 presents the values of the statistical indicators evaluated for all trials with the $k - \varepsilon$ and $R_{ij} - \varepsilon$ models.

In general, the CFD model underestimates the observations, this underestimation being slightly more pronounced when the anisotropic model is used.
<figure><img src="image_18.png" /><figcaption>Figure 16: Comparison between the observed concentration profile and the one calculated with the $k-\epsilon$ and $R_{ij}-\epsilon$ models on all the PID receptors for trial No.2692157.</figcaption></figure>

This underestimation is illustrated by the values of FB in Table 7 and is clearly highlighted in Figure 17 which shows that the scatter plot tends to spread out below the acceptable zone of FAC2.

Regarding the values of the statistical indicators presented in Table 7, it can be seen that the numerical calculations predict with a good precision (correlation higher than 0.8) the physical phenomena present during these experimental trials. As for the average values of MG and VG, they exceed the acceptable margins whatever the turbulence model used. However, their values remain close to the threshold when a $k-\epsilon$ is used. The $R_{ij}-\epsilon$ model, although a second-order model, does not differ significantly from the $k-\epsilon$ (first-order) model in terms of modeling dispersion for this trial. Both models obtain a score of 3, which means that the criterion for only three of the six statistical measures is satisfied.

The analysis of the results of the stable trials, in particular the two trials No.2671934 (with an Obukhov length L=5.8m) and No.2672033 (L=4.8m), allows to make several observations:

* For trial No.2671934, no turbulence model is able to perfectly predict the position and the amplitude of the maximum concentrations, even if the shape of the predicted concentration profile is similar to the one observed in the trial (profile not shown here). Consequently, none of the models used achieves acceptable performance (score ≥ 3). Overall, a shift in peak position and an underestimation of the concentration measured by the receptors are observed across this profile, as shown in Figure 17. This shift in concentration values is not detected with the trial No.2672033, even
though it has a more stable atmospheric stability than trial No.2671934. However, this last one took place in the presence of wind at a speed (1.6m/s) lower than that of the first one (2.69m/s) at the entrance of the domain. This finding confirms the conclusions of Degrazia et al. (2011) that the usual Eulerian RANS dispersion models are not suitable for very stable atmospheric conditions for two main reasons: (i) these models have limitations in terms of modeling when the advection velocity in the transport equation approaches very low values (U<2m/s) due to the meandering phenomenon; (ii) these models take into account pollutant transport by convection, molecular and turbulent diffusion, but not the large-scale phenomena that are dominant during low wind sequences. This may also explain the large underestimation of concentration values observed by receptors placed far from the source.

* It appears from the predictions of the trial No.2672033 that the use of the constants suggested by Sadek et al. (2013) for the $R_{ij} - \varepsilon$ model (retained as a result) significantly modifies the concentration profile obtained with the standard model and allows to obtain results closer to the experimental observations.

### 5.4.2. Comparison with numerical results of the Fluidyn-Panache code

In another study on the MUST experiment, Kumar et al. (2015) carried out an evaluation of CFD simulation results using the Fluidyn-Panache code (Fluidyn-Panache, 2010), which consists of a three-dimensional flow model and an Eulerian dispersion model. Many combinations of wind inflow boundary conditions for the wind and turbulence profiles were considered for the simulations. The turbulence model employed in the simulations is the isotropic $k-\varepsilon$ model, which includes a modified value of the $C_{\varepsilon3}$ coefficient to account for buoyancy effects under stable atmospheric conditions. In order to compare the evaluation results of the present study, statistical measures of performance are computed for trial No.2692157 (Table 6) as well as for the 4 trials on average (Table 7). Only the best performing combinations, i.e. the simulations with the highest score, are taken into account in this statistical study.

For trial No.2692157, we find that the Fluidyn-Panache model provides equally good predictions as the Code_Saturne model and that the criteria for all six statistical indicators are met. However, an even higher percentage of predicted concentrations (87.5%) are within a factor of two of the observations (Table 6).

For all trials, the Fluidyn-Panache model predicts 77.1% of the points in a FAC2 of the observations, while the present model simulates 61.6%, 58.8% of the points in a FAC2 of the observations with the $k-\varepsilon$ and $R_{ij}-\varepsilon$ models, respectively. As for the FB values, only the Fluidyn-Panache model satisfies the criterion (0.138) while the present model underestimates the concentration
<figure><img src="image_19.png" /><figcaption>Figure 17: Scatter plots (on the left) and q-q plots (on the right) illustrating observed and predicted concentrations using the k-ε model (a,b) and the R<sub>ij</sub> - ε model (c,d) for the four trials.</figcaption></figure>
values with a value of 0.473 with the $k - \varepsilon$ model and 0.528 with $R_{ij} - \varepsilon$. The values of NMSE for the simulated concentrations at all receptors for the four cases are respectively 1.513 and 1.493 with the $k - \varepsilon$ and $R_{ij} - \varepsilon$ models. These values are almost 3.75 times higher than the value (0.41) reported by Kumar et al. (2015) with the Fluidyn-Panache model.

<table><caption>Table 7: Statistical performance measures assessed for all trials studied.</caption><thead><tr><th></th><th>FB</th><th>MG</th><th>NMSE</th><th>VG</th><th>FAC2(%)</th><th>COR</th></tr></thead><tbody><tr><td>Code_Saturne, k - ε</td><td>0.473</td><td>1.815</td><td>1.513</td><td>30.92</td><td>61.6</td><td>0.833</td></tr><tr><td>Code_Saturne, R<sub>ij</sub> - ε</td><td>0.528</td><td>2.215</td><td>1.493</td><td>106</td><td>58.8</td><td>0.835</td></tr><tr><td>Fluidyn-Panache, k - ε</td><td>0.138</td><td>1.118</td><td>0.41</td><td>1.895</td><td>77.1</td><td>0.933</td></tr></tbody></table>

# 6. Conclusions

This study presents 3-D simulations of the dispersion of a pollutant under neutral and stable atmospheric conditions. The Code_Saturne model is evaluated with concentration measurements obtained from four trials of the MUST experimental campaign in a complex urban field. One of these tests is performed under neutral conditions, while the other tests are performed under relatively stable to very stable atmospheric conditions. Simulations are done using two turbulence models, $k - \varepsilon$ and $R_{ij} - \varepsilon$ LRR, and scalar turbulent fluxes are evaluated using SGDH and GGDH algebraic models, as well as the differential transport equation (DFM). Qualitative and quantitative analyses show that the performance of the models is reasonably acceptable in terms of capturing multiple concentration peaks and the distribution of the hourly average concentration.

Numerical results from simulations on the moderately stable test show that the turbulence models have an overall underestimation of the mean concentration. The comparison between the $k - \varepsilon$ and $R_{ij} - \varepsilon$ models shows that the former tends to generate an overproduction of turbulence near the ground, resulting in a significantly wider plume spread. The isotropic model ($k - \varepsilon$) predicts approximately 75% of the concentration values within a factor of two, whereas the $R_{ij} - \varepsilon$ model predicts 62.5% and 60.4% when used with the SGDH and GGDH models, respectively. The use of the differential transport equation (DFM) to solve scalar turbulent fluxes improves predictions with nearly 69% of the concentration values within a factor of two (FAC2) of the observations. It also provides a better correlation with experimental measurements, compared to algebraic models such as SGDH and GGDH. Regardless of its computational cost, this model proves to be very adequate to study atmospheric dispersion in complex terrain with moderate to low stable stability conditions.

In all simulation cases, the CFD model consistently underestimates the observations. Overall, when using the isotropic first-order $k - \varepsilon$ model, the CFD model predicts a slightly higher percentage of concentrations within a factor of two (FAC2) of the observations. Additionally, it is noteworthy that the model's
722
723

performance tends to decrease with lower wind speeds, regardless of atmospheric stability.

724

725

726

727

728

Further studies will be carried out on the different algebraic models and on the scalar transport equation, in order to better evaluate the efficiency of these models, especially under extremely stable atmospheric conditions. It might also be beneficial to fine-tune the various constants employed in the formulation of these models to enhance their compatibility with dispersion-related challenges.

729

730

731

732

733

734

735

As the Reynolds stress tensor has been used in its isotropic form, relying solely on the tensor trace or kinetic energy, this may explain the discrepancies obtained with the $R_{ij} - \varepsilon$ model. Therefore, future investigations will include a re-evaluation of calculations using this model by essentially incorporating the other components of the Reynolds tensor into the boundary conditions, particularly the momentum flux $\overline{U'}\overline{W'}$, to better account for anisotropy and enhance model performance.

736

737

738

739

As the approach used to model thermal turbulent flux is the simple gradient model (SGDH), it will be necessary to test other algebraic models or even the transport equation to verify the effect of thermal turbulence anisotropy on the turbulent flux of the concentration scalar.

740

## Acknowledgement

741

742

743

744

745

The authors would like to thank the Defense Threat Reduction Agency (DTRA) for providing access to the MUST field experiment dataset. Authors gratefully acknowledge the DGA (Direction Générale de l'Armement), Ministère des Armées, France, for its financial support. Our sincere thanks to Stéphane Burkhart (DGA) for many helpful discussions.

746

## References

747

748

749

Alinot, C., Masson, C., 2005. $k-\varepsilon$ model for the atmospheric boundary layer under various thermal stratifications. Journal of Solar Energy Engineering-transactions of The Asme - J. Sol. Energy Eng. 127. doi:10.1115/1.2035704.

750

751

752

753

Archambeau, F., Méchitoua, N., Sakiz, M., 2004. Code Saturne: A Finite Volume Code for the computation of turbulent incompressible flows - Industrial Applications. International Journal on Finite Volumes 1. URL: https://hal.science/hal-01115371.

754

755

756

757

Bahlali, M., Dupont, E., Carissimo, B., 2019. Atmospheric dispersion using a lagrangian stochastic approach: Application to an idealized urban area under neutral and stable meteorological conditions. Journal of Wind Engineering and Industrial Aerodynamics 193. doi:10.1016/j.jweia.2019.103976.

758

759

760

Beljaars, A., Holtslag, B., 1991. Flux parameterization over land surfaces for atmospheric models. J. Appl. Meteor. 30, 327-341. doi:10.1175/1520-0450(1991)030<0327:FPOLSF>2.0.CO;2.
761 Biltoft, C.A., 2001. Customer report for mock urban setting test. DPG Document, 160-000.

763 Businger, J.A., Wyngaard, J.C., Izumi, Y., Bradley, E.F., 1971. Flux-profile relationships in the atmospheric surface layer. Journal of the Atmospheric Sciences 28, 181-189. URL: https://doi.org/10.1175/1520-0469(1971)028<0181:FPRTA>2.0.CO;2.

767 Chang, J., Hanna, S., 2004. Air quality model performance evaluation. Meteorology and Atmospheric Physics 87, 167-196. doi:10.1007/s00703-003-0070-7.

770 Daly, B.J., Harlow, F.H., 1970. Transport Equations in Turbulence. The Physics of Fluids 13, 2634-2649. URL: https://doi.org/10.1063/1.1692845, doi:10.1063/1.1692845.

773 Degrazia, G., Timm, A., Moreira, V., Roberti, D., 2011. Meandering Dispersion Model Applied to Air Pollution. doi:10.5772/16892.

775 Dehoux, F., 2012. Modélisation statistique des écoulements turbulents en convection forcée, mixte et naturelle. Ph.D. thesis. Université de Poitiers. URL: http://www.theses.fr/2012POIT2276.

778 Duynkerke, P.G., 1988. Application of the $e-\epsilon$ turbulence closure model to the neutral and stable atmospheric boundary layer. Journal of Atmospheric Sciences 45, 865-880. doi:https://doi.org/10.1175/1520-0469(1988)045<0865:AOTTCM>2.0.CO;2.

782 Dyer, A.J., 1974. A review of flux-profile relationships. Boundary-Layer Meteorology 7, 363-372. doi:10.1007/BF00240838.

784 Dyer, A.J., Hicks, B.B., 1970. Flux-gradient relationships in the constant flux layer. Quarterly Journal of the Royal Meteorological Society 96, 715-721. doi:https://doi.org/10.1002/qj.49709641012.

787 Fluidyn-Panache, 2010. User Manual. 7 ed. Fluidyn France/TRANSOFT International. Version 4.0.

789 Franke, J., Hellsten, A., Schlünzen, H., Carissimo, B.t., 2007. Cost action 732 quality assurance and improvement of microscale meteorological models, in: COST Action.

792 Gorlé, C., Beeck, J., Rambaud, P., van Tendeloo, G., 2009. CFD modelling of small particle dispersion: The influence of the turbulence kinetic energy in the atmospheric boundary layer. Atmospheric Environment 43, 673-681. doi:10.1016/j.atmosenv.2008.09.060.

796 Gryning, S.E., Batchvarova, E., Brümmer, B., Jørgensen, H., Larsen, S., 2007. On the extension of the wind profile over homogeneous terrain beyond the surface layer. Boundary-Layer Meteorology 124, 251-268. doi:10.1007/s10546-007-9166-9.
Hicks, B.B., 1976. Wind profile relationships from the 'wangara' experi-
ment. Quarterly Journal of the Royal Meteorological Society 102, 535-551.
doi:https://doi.org/10.1002/qj.49710243304.

Kavian Nezhad, M.R., Lange, C., Fleck, B., 2022. Performance evaluation of the RANS models in predicting the pollutant concentration field within a compact urban setting: Effects of the source location and turbulent schmidt number. Atmosphere 13, 1013. doi:10.3390/atmos13071013.

Kerschgens, M., Nölle, C., Martens, R., 2000. Comments on turbulence pa-
rameters for the calculation of dispersion in the atmospheric boundary layer.
Meteorologische Zeitschrift 9, 155-163. doi:10.1127/metz/9/2000/155.

Knaus, H., Hofsäß, M., Rautenberg, A., Bange, J., 2018. Application of different turbulence models simulating wind flow in complex terrain: A case study for the windfors test site. Computation 6. doi:10.3390/computation6030043.

Konow, H., 2014. Tall Wind Profiles in Heterogeneous Terrain. Ph.D. thesis.

Kumar, P., Feiz, A.A., 2016. Performance analysis of an air quality CFD model
in complex environments: Numerical simulation and experimental validation
with emu observations. Building and Environment 108, 30–46. doi:https:
//doi.org/10.1016/j.buildenv.2016.08.013.

Kumar, P., Feiz, A.A., Ngae, P., Singh, S.K., Issartel, J.P., 2015. CFD
simulation of short-range plume dispersion from a point release in an ur-
ban like environment. Atmospheric Environment 122, 645–656. doi:https:
//doi.org/10.1016/j.atmosenv.2015.10.027.

Lateb, M., Meroney, R., Yataghene, M., Fellouah, H., Saleh, F., Boufadel, M., 2016. On the use of numerical modelling for near-field pollutant dispersion in urban environments - a review. Environmental Pollution 208, 271-283. doi:https://doi.org/10.1016/j.envpol.2015.07.039. special Issue: Urban Health and Wellbeing.

Launder, B., Reece, G., Rodi, W., 1975. Progress in the development of a reynolds stress turbulence closure. Journal of Fluid Mechanics 68, 537 – 566. doi:10.1017/S0022112075001814.

Launder, B., Spalding, D., 1974. The numerical computation of turbulent flow
computer methods. Computer Methods in Applied Mechanics and Engineer-
ing 3, 269-289. doi:10.1016/0045-7825(74)90029-2.

Lecocq, Y., 2008. Contribution à l'analyse et à la modélisation des écoulements
turbulents en régime de convection mixte - Application à l'entreposage des
déchets radioactifs. Ph.d. thesis. Université de Poitiers. URL: https://
theses.hal.science/tel-00375441.
Lecocq, Y., Manceau, R., Bournaud, S., Brizzi, L.E., 2008. Modelling of the turbulent heat fluxes in natural, forced and mixed convection regimes, in: 7th ERCOFTAC Int. Symp. on Eng. Turb. Modelling and Measurements, Limassol, Cyprus, 2008., Limassol, Cyprus. pp. –. URL: https://hal.science/hal-00406130.

Milliez, M., Carissimo, B., 2007. Numerical simulations of pollutant dispersion in an idealized urban area, for different meteorological conditions. Boundary-Layer Meteorology 122, 321–342. doi:10.1007/s10546-006-9110-4.

Milliez, M., Carissimo, B., 2008. Computational fluid dynamical modelling of concentration fluctuations in an idealized urban area. Boundary-Layer Meteorology 127, 241–259. doi:10.1007/s10546-008-9266-1.

Na, A., Gryning, S.E., Hasager, C., Courtney, M., 2009. Extending the wind profile much higher than the surface layer. European Wind Energy Conference and Exhibition 2009, EWEC 2009 2.

Panofsky, H.A., Dutton, J.A., 1984. Atmospheric turbulence. models and methods for engineering applications. New York: Wiley.

Pasquill, F., 1972. Some aspects of boundary layer description. Quarterly Journal of the Royal Meteorological Society 98, 469–494. doi:https://doi.org/10.1002/qj.49709841702.

Peña, A., Gryning, S.E., Hasager, C., 2009. Comparing mixing-length models of the diabatic wind profile over homogeneous terrain. Theoretical and Applied Climatology 100, 325–335. doi:10.1007/s00704-009-0196-8.

Reiminger, N., 2020. 3D Modelling of Air Pollution at Urban Scale. Ph.D. thesis. doi:10.13140/RG.2.2.10037.91366.

Rodi, W., 1993. Turbulence models and their application in hydraulics. CRC Press.

Rossi, R., Iaccarino, G., 2009. Numerical simulation of scalar dispersion downstream of a square obstacle using gradient-transport type models. Atmospheric Environment 43, 2518–2531. doi:https://doi.org/10.1016/j.atmosenv.2009.02.044.

Rossi, R., Philips, D.A., Iaccarino, G., 2009. Numerical simulation of scalar dispersion in separated flows using algebraic flux models, in: Turbulence Heat and Mass Transfer 6. Proceedings of the Sixth International Symposium On Turbulence Heat and Mass Transfer, Begel House Inc.

Sadek, R., Soulhac, L., Brocheton, F., Buisson, E., 2013. Adaptation of the reynolds stress turbulence model for atmospheric simulations, in: HARMO 15 conference.
Salim, S.M., Buccolieri, R., Chan, A., Di Sabatino, S., 2011. Numerical simulation of atmospheric pollutant dispersion in an urban street canyon: Comparison between rans and les. Journal of Wind Engineering and Industrial Aerodynamics 99, 103–113. doi:https://doi.org/10.1016/j.jweia.2010.12.002.

Schatzmann, M., Olesen, H., Franke, J., 2010. COST 732 model evaluation case studies: approach and results.

Stull, R.B., 1988. An Introduction to Boundary Layer Meteorology. Springer Netherlands.

Tominaga, Y., Stathopoulos, T., 2010. Numerical simulation of dispersion around an isolated cubic building: Model evaluation of RANS and LES. Building and Environment 45, 2231–2239. doi:https://doi.org/10.1016/j.buildenv.2010.04.004.

Tominaga, Y., Stathopoulos, T., 2011. CFD modeling of pollution dispersion in a street canyon: Comparison between les and rans. Journal of Wind Engineering and Industrial Aerodynamics 99, 340–348. doi:https://doi.org/10.1016/j.jweia.2010.12.005. the Fifth International Symposium on Computational Wind Engineering.

Tominaga, Y., Stathopoulos, T., 2012. CFD modeling of pollution dispersion in building array: Evaluation of turbulent scalar flux modeling in rans model using les results. Journal of Wind Engineering and Industrial Aerodynamics 104–106, 484–491. doi:https://doi.org/10.1016/j.jweia.2012.02.004. 13th International Conference on Wind Engineering.

Tominaga, Y., Stathopoulos, T., 2017. Steady and unsteady rans simulations of pollutant dispersion around isolated cubical buildings: Effect of large-scale fluctuations on the concentration field. Journal of Wind Engineering and Industrial Aerodynamics 165, 23–33. doi:https://doi.org/10.1016/j.jweia.2017.02.001.

Tominaga, Y., Stathopoulos, T., 2018. CFD simulations of near-field pollutant dispersion with different plume buoyancies. Building and Environment 131, 128–139. doi:https://doi.org/10.1016/j.buildenv.2018.01.008.

Yee, E., Biltoft, C., 2004. Concentration fluctuation measurements in a plume dispersing through a regular array of obstacles. Boundary-Layer Meteorology 111, 363–415. doi:10.1023/B:BOUN.0000016496.83909.ee.

Yu, H., Thé, J., 2017. Simulation of gaseous pollutant dispersion around an isolated building using the k-ω SST (shear stress transport) turbulence model. Journal of the Air & Waste Management Association 67, 517–536. doi:10.1080/10962247.2016.1232667.

Zilitinkevich, S., Baklanov, A., 2002. Calculation of the height of the stable boundary layer in practical applications. Boundary-Layer Meteorology 105, 389–409. doi:10.1023/A:1020376832738.