Numerical approximation to the effects of the atmospheric stability conditions on the dispersion of pollutants over flat areas

Using the Computational Fluid Dynamics technique (CFD), we explored the effects of the atmospheric stability conditions on the dispersion of solid and gas-phase pollutants emitted from an area source located on a flat region. As an application, the dispersion of pollutants emitted from roads located on flat terrains was considered. Toward that end, we set up a model that describes the dispersion of air pollutants in a small region (< 1 km long) near the ground surface (< 250 m high). It consists of a neutrally stratified model modified to account for the atmospheric stability effects by imposing the near-ground stratification through the Monin–Obukhov similarity theory and the k–ε turbulence model adjusted for each atmospheric stability condition. Using this model, we simulated the dispersion of pollutants emitted from the road and plotted the resulting downwind concentrations in terms of dimensionless numbers. Results from our CFD-based model were highly correlated (R2 > 0.95) with the SF6 concentrations measured downwind a line source of this trace gas by the U.S. National Oceanic Atmospheric Administration in 2008 under different conditions of atmospheric stability. Numerical and experimental results showed that, under any of the stability conditions explored, the near-road pollutant concentrations are highly correlated (R2 > 0.87) to the concentrations observed under neutral conditions. When the atmosphere is extremely stable, those concentrations were up to 12 times higher than those observed under neutral conditions. We report the constant of proportionality obtained for every stability condition.

. Parameters used to describe atmospheric stability conditions. Richardson number Ri = . Where , ρ, and C p are the air thermal conductivity, density, and specific heat, respectively. T o and H o are the surface temperature and heat flux, respectively. z is height, and g is the gravity constant. u and u* are the wind and friction velocity, respectively. According to the Pasquill-uifford stability classification classes D applies to heavily overcast skies, at any wind speed day or night. www.nature.com/scientificreports/ of 2-5 from the highly unstable to the stable atmospheric conditions 12 . All these experimental works agree that the pollutants concentrations are higher under stable atmospheric conditions than under unstable conditions. However, the magnitude of this effect varies among experimenters. Furthermore, these variations could be due to the sole effect of variations in wind speed rather than variations in the atmospheric stability conditions. The effects of the atmospheric stability conditions on pollutant dispersion can also be studied analytically. The physics phenomena occurring under the different atmospheric stability conditions are described by the mass, momentum, and energy conservation equations 5 . Nevertheless, the task of solving those equations is challenging because they are nonlinear partial differential equations.

Atmospheric Stability Class
Authors have faced this challenge by solving these governing equations with a high time and space resolution to capture the turbulence phenomena (high-frequency fluctuations in speed) at the mesoscale range. In this way, they grasp the short-and long-range turbulent eddies of the different atmospheric stability conditions. This approach is known as direct numerical simulation (DNS). The CFD engineering and geophysical modeling communities have converged in that this is the best approach to model the atmosphere dynamics 13 . However, this approach is computationally expensive and has not been used to study the dispersion of pollutants systematically.
There are several challenges when attempting to use this DNS alternative for the study of the dispersion of pollutants. The main ones are related to the differences in the spatial and temporal scale of the atmospheric dynamics and of the dispersion phenomena: • The area of interest in the study of pollutant dispersion is near the emission source, where concentrations are the highest. i.e., the interest is at the micro-scale level (< 1 km). Then, the challenge is to reproduce the atmosphere dynamics, which is a regional phenomenon (~ 100 km), in a reduced computational domain. • The time scale of interest in the study of pollutant dispersion is of hours or days. Time average concentrations are of interest rather than instantaneous concentrations. Then, the challenge is to obtain long-term concentrations at a reasonable computational time.
Some authors have proposed the coupling between meso-and micro-scale models where results from the mesoscale model are used as boundary conditions to the micro-scale model 14 . However, the issues described above remain still unresolved, at least for practical applications 15 .
As an alternative to the DNS approach, Reynolds (1895), working at the micro-scale range, described the instant velocity as a short-time-average velocity plus a fluctuating velocity. Then, he transformed the instant governing equations in terms of these two velocity components, which produced additional unknown variables termed as Reynolds stresses. Those stresses represent turbulent transport of momentum and energy. Usually, it is assumed that those stresses can be expressed as a linear combination of the velocity gradients. However, additional equations are still needed to close the description of the turbulence problem. Turbulence models based on this approach are known as Reynolds Average Navier Stokes (RANS) models. Among them, the k-ε turbulence model has been the most widely used. The k-ε model adds an equation for the turbulent kinetic energy (k) and one equation for the turbulent energy dissipation rate (ε). Adding these turbulence models to the momentum and energy equations, CFD engineers have been able to reproduce the experimental observations related to the interaction of the flow field with solid bodies in terms of pressure gradients (aerodynamics) and in terms of heat and mass transfer.
All these turbulence models do not reproduce the physics of the turbulence phenomena, but they simulate the average turbulence effects on the flow field. The main drawback of this approach is that there is not a single model for every application. Furthermore, every turbulence model must be adjusted to the specific application via experimental calibration.
Following this approach, the k-ε turbulence model has been used to reproduce the effects of the atmospheric dynamics on the average flow field on the SBL 16 . Most of the work developed simulating near-road air pollution at the microscale level (< 1 km) adopts this alternative. However, these studies assume uniform temperature and isotropic turbulence, which are features of a neutral atmosphere [17][18][19] . Using this alternative, Huertas and Prato 17 used meteorological data as input to their Near-Road CFD (NR-CFD) model and simulated the hourly dispersion of particles near two unpaved roads located on a flat region. They obtained daily and monthly average values of particle concentration that were highly correlated to measurements obtained at several points downwind the roads.
The CFD based models that adopt the neutral atmosphere assumption have been criticized for not including the atmospheric dynamics and its different conditions of stability. It has also been counter-argued that in the very near ground surface, the flow field is highly influenced by the mechanical turbulence generated by the presence of physical obstacles (e.g., buildings), which disturb the free wind flow, while the atmospheric stability conditions are mesoscale or regional phenomena 20 , and therefore its influence on the dispersion of pollutants is minor at the very near ground surface.
CFD-based models can be used to evaluate the effects of the different atmospheric stability conditions on the dispersion of pollutants systematically. However, few works have been conducted with this objective in mind. Using CFD, Pieterse et al. 21 evaluated the effects of three atmospheric stability conditions on the wind flow over flat terrain. They found the parameters for the k-ε turbulence model that best describe those atmospheric conditions. However, they did not use their model to study the effects of atmospheric stability on the dispersion of pollutants. Steffens et al. 18 replicated in CFD the experimental results of the NRTS08 campaign described above. They modeled the atmospheric stability conditions by implementing the Monin Obukhov Similarity (MOS) theory and evaluated the performance of two turbulence models: Reynolds Average Navier Stokes (RANS) and Large Eddie Simulation (LES). They found that both models produced similar results. Nevertheless, they did not quantify the effect of the atmospheric stability conditions on the dispersion of the SF 6  www.nature.com/scientificreports/ of the MOS theory with the k-ε turbulence model produces a rough description of the atmospheric dynamics in the SBL, which could be useful for applications likes the study of the dispersion of pollutants and for the wind engineering community. However, we highlight that this approximation is not intended for the study of the atmospheric dynamics. We propose to advance Pieterse's work, and use said description of the atmosphere dynamics to study the effects of the different stability conditions on the dispersion of pollutants. Specifically, the objective of this work is to quantify the effects of atmospheric stability conditions on the dispersion of pollutants downwind an emission source, such as roads, located in flat regions.
In the process of pursuing this objective, the following contributions to new knowledge were developed: • The proposed approach was implemented in the NR-CFD model. Results from said model reproduces the experimental measurements obtained in the NRTS08 campaign studying the dispersion of a line source of SF 6 under different atmospheric stability conditions (Fig. 4b). In this manuscript, the NR-CFD model is described in a way that it could be reproduced by others. • An approximation to the quantification of the effects of the atmospheric stability conditions on the ground pollutant concentration downwind an area emission source was obtained by systematically using the NR-CFD model and experimental data (Fig. 5f). • This manuscript reports that the ground pollutant concentration downwind an emission source under any atmospheric stability condition is highly correlated to the concentration observed under neutral conditions ( Fig. 5e-g) when it is expressed in terms of dimensionless numbers (C* vs. x*). Furthermore, that said C* vs.
x* profile is independent of variations in the wind speed, mass emission rate, and of the pollutant's nature.

Methodology
Aiming to study the effects of the atmospheric stability conditions on the dispersion of the pollutants emitted from an area source located over a flat surface, (i) we selected the simplest possible representative case. i.e., a road on a horizontal area without any obstacle to the wind flow. Then, (ii) an approximation to the physics occurring in the atmosphere very near the ground surface (< 250 m high), in a small domain (< 1 km long), was implemented, by solving via CFD the mass, momentum, and energy equation, coupled to the MOS theory and the appropriate k-ε turbulence model for each atmospheric stability condition. (iii) Using this model (NR-CFD model), the dispersion of pollutants emitted from a road was simulated and the resulting downwind concentrations were expressed in terms of dimensionless numbers. (iv) Results from the NR-CFD model were compared with the experimental measurements obtained in the NRTS08 campaign. (v) Finally, the experimental and simulated results were used to compare the concentrations downwind the road under different atmospheric stability conditions with the ones observed under neutral conditions. Next, we will describe each of these steps.

Case of study.
In this work, we studied the dispersion of the pollutants emitted from a road located on a flat region without the presence of any obstacle to the wind free flow (Fig. 1). As stated before, this case is an illustrative example of a more general case, which is the dispersion of the pollutants emitted from an area source located over a flat surface. We are interested in the region near the source where the highest concentrations occur. This means that the interest is located very near the ground surface (z < 250 m) and close to the emission source (length < 1 km). For this particular case, we defined the computational domain shown in Fig. 1, which is a box of dimensions 5 m, 330.5 m, and 50 m in the crosswind, along-wind, and vertical direction, respectively. These dimensions were selected according to Franke et al. 22 , who recommended the minimum dimensions of the computational domain in the way that the walls do not interfere with the dispersion of the pollutants. www.nature.com/scientificreports/ Simulation of the atmosphere dynamics. Next, we describe our approach to simulate, at microscale level (length < 1 km), the physics occurring in the near-surface atmospheric boundary layer (height < 250 m), under different stability conditions, using state of the art CFD solvers. The first step in the CFD simulation of pollutant dispersion near roads is the simulation of the atmospheric dynamics. The physics phenomena occurring in the atmosphere are described by the well-known mass, momentum, and energy governing equations 5 . These equations can be solved numerically using the computational fluid dynamics (CFD) technique. In this technique, the space domain is divided into small volumes, and these equations are solved over each volume. Besides the discretization of the computational domain into a large number of finite volumes, the simulation of the physics phenomena occurring in the atmosphere requires the specification of the conditions under which those equations need to be solved.
In the study of dispersion of pollutants, as in several other applications, researchers are interested in short (1 h) and long-term (~ 1 year) average results. Implicitly, this interest involves a pseudo-steady-state approximation where the daily evolution of the atmosphere is studied as a set of successive states where at each time interval, the atmosphere is assumed to be under a steady-state condition. This assumption is acceptable if the time step of the simulation is smaller than the response time of the SBL to changes in radiative forcing. The SBL's response to changes in surface radiative forcing is a large-scale phenomenon with long response times (~ 1 h) 6,23,24 . In practice, 1-h intervals are suitable, as meteorological data is reported in this way.
When modeling the atmospheric dynamics over a flat surface without any perturbation to the wind free flow, this pseudo steady state assumption implies: • The condition of horizontal homogeneity. The steady-state assumption of the atmosphere over a flat surface within a small computational domain implies that the value of any property remains the same at any position with the same height. It requires that, in the absence of any perturbation, the vertical wind and temperature profiles at the inlet remain the same at the outlet. This condition of horizontal homogeneity agrees with the observation that at any given time, for example, the ground temperature is the same everywhere within the microscale domain. • The inlet boundary condition is known and they satisfy the governing equations. The steady-state assumption requires that at every vertical position at the inlet, the value of each variable is a known input data. The use of constant values or of any arbitrary function for the physical properties at the inlet will lead to a violation of the horizontal homogeneity condition. Some researchers 18 have measured the vertical wind and temperature profiles and used them as the inlet boundary condition. However, in practice, researchers only have available meteorological data, where wind speed and temperature are measured at a single height at each time interval. Later in this section, we will describe how the inlet input data can be obtained from the meteorological data. • Agreement between the inlet and the ground boundary conditions. The condition of horizontal homogeneity requires that the ground boundary conditions for momentum and heat transfer agree with the inlet boundary condition. Otherwise, wind speed and air temperature will change along with the computational domain.
Inlet boundary condition. The objective is to specify the inlet boundary condition based on the measurements of temperature and wind speed at 10 m above ground, in a way that the result of the simulation of the atmospheric dynamics, under the steady-state assumption, satisfies the homogeneity condition. That is, the vertical inlet profile for wind speed and temperature at the inlet and outlet should be the same. The first alternative is to start with any arbitrary vertical profile for wind speed and temperature and simulate the wind flow over an extremely long flat surface subject to a desired constant ground heat transfer condition. The resulting profiles will satisfy the conditions specified above. This process is cumbersome, especially since it should be carried out for each time-step of meteorological data.
The second alternative is to use the knowledge gained on the atmospheric dynamics to specify the inlet boundary condition. In the 1920s, studying the fluid flow on a flat surface, Von Karman observed that the fluid flow on that surface exhibited a parabolic profile and that when the said profile is expressed in terms of dimensionless numbers, it does not depend on the downwind position. Monin-Obukhov, in the 1950s, extended Von Karman's concepts of fluid mechanics to the modeling of the atmosphere dynamics. They specified the shape of wind speed profiles for several stability conditions (Eqs. [1][2][3][4][5]. Today that approach is known as the Monin-Obukhov similarity (MOS) theory. It states that any mean flow quantity in the SBL (e.g., momentum or energy), when normalized by an appropriate scaling parameter, is a unique function of ζ = z/L where ζ is termed as the buoyancy parameter 5 , and L is the Monin-Obukhov length, which describes the atmospheric stability condition (Table 1).
Based on the observation that (i) in the SBL the wind speed (u) varies logarithmically with height and that (ii) the surface roughness forces the mean wind speed to be zero at the ground surface, the MOS theory states that the vertical profile for wind speed in the SBL follows Eq. (1) 23 .
where u* is the friction velocity, which is calculated from u * = √ τ w /ρ ; k is the Von Karman constant which has a value of 0.41 5 ; and z o is the surface roughness that for the case of short grass, the recommended value is 0.03 m 5 . ψ m is the similarity function for the flux of momentum (Eqs. 3 and 4). For 0 < z < z o , the wind speed is zero. In Eqs. (4) and (5), α = (1 − 15z/L) 1/4 . We selected a set of values for L and used the resulting vertical profiles of wind speed as the inlet boundary condition. www.nature.com/scientificreports/ The MOS theory also specifies the vertical profiles for the ambient temperature (Eq. 2, 25 ). In this case θ is the potential temperature which is calculated from θ = T(P/P 0 ) R/c p ; θ o is the ground potential temperature; θ* is the potential temperature which most of the time is equal to T 0 , and ψ h is the similarity function for heat flux described by Eqs. (3) and (5). We used these profiles as the inlet boundary condition of the computation domain.
Under stable atmospheric conditions, both profiles tend to become linear for large values of ζ . Under unstable conditions ψ m and ψ h are positive, and therefore, the velocity and temperature profiles in the surface layer are more curvilinear as the condition of instability increases 5 .
Use of the wall function for the ground boundary condition. The flow field over the ground surface should satisfy the non-slip condition. The first alternative to implement this condition is to specify a wind speed equal to zero at the surface and let the model to solve for the vertical wind speed. This alternative requires finite volumes much smaller than the boundary layer thickness, which is expensive computationally speaking.
Instead, best practices in CFD recommends the use of the near-wall treatment. This second alternative uses again the knowledge gained on fluid dynamics to specify the vertical wind speed profile as the ground boundary condition. It reduces substantially the number of elements required near the surface. Best practices in CFD suggests that the size of the elements in the ground surface should be at most half of the estimated boundary layer. CFD modelers developed the metric y + to define the size of the elements near the surface to assure congruency of the boundary condition with the solution above the ground surface. Following this practice, the simulation of wind flow over a grass surface required finite volumes of ~ 30 cm near the ground surface, which is equivalent to a y + value of 214.2. This value satisfies the log-law of modeling that recommends 30 < y + < 1000 26 . To fulfill the horizontal homogeneity condition, we used as ground boundary conditions the same velocity profile used at the inlet boundary condition. Similarly, we used the vertical profile of temperature as the ground boundary condition for the energy equation.
Other boundary conditions. As described before, the height of the computational domain should be selected as the minimum where vertical interaction (other than turbulences) is negligible. Then, following the best CFD practices, we selected as the upper boundary condition the symmetry condition (zero gradients). When simulations were carried out in 3D, we proceeded in a similar way and used the periodic boundary condition, which is used when the solution in the crosswind direction remains the same with periods equal to the width of the computational domain. There was no need to specify any surface roughness at the road because it was modeled as a source of pollutants. That is, we specified a uniform velocity inlet at the road. Physically it means that the incoming wind flow faces a perpendicular ascending wind flow with a high pollutant concentration.
Turbulence model. As stated in the introduction section, the governing equations for momentum and energy need to be solved with a high time and space resolution to capture the turbulence phenomena (high-frequency fluctuations in speed). This approach is known as direct numerical simulation (DNS), which is computationally expensive. As an alternative, a model that describes the turbulence phenomena is added to the momentum and energy equations. Turbulence models for the SBL should account for both shear and buoyancy produced turbulence 21 . The k-ε turbulence model meets this requirement and has been widely used in SBL-related studies 16 . The k-ε model adds an equation for the turbulent kinetic energy (k, Eq. 6) and an equation for the turbulent energy dissipation (ε, Eq. 7). Both ensure the closure of the system of equations that describe the physics phenomena occurring in the SBL.
In these equations, G k (Eq. 8) represents the generation of turbulent kinetic energy due to the mean velocity gradient, and G b (Eq. 9) is the generation of turbulent kinetic energy due to buoyancy. Y M represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. µ t is the turbulent www.nature.com/scientificreports/ viscosity, and C 1ε , C 2ε and C 3ε , are constants. σ k and σ ε are the turbulent Prandtl numbers for k and ε , respectively. S k and S ε are source terms of k and ε 27 .
τ ij is the shear stress in the direction i and perpendicular to the plane j; u' i is the fluctuating wind speed in the i direction; β is the fluid thermal expansion coefficient; T is temperature, and Pr t is the energy turbulent Prandtl number 27 . In Eq. (9), g i /c p is the adiabatic lapse rate. An unstable condition will occur when the actual ∂T/∂x i is greater than the g i /c p . When ∂T/∂x i is less than g i /c p , turbulence will be suppressed from the atmosphere, and the SBL will exhibit a condition of stability 16 . Pieterse and Harms 22 measured k and ε at different heights for different atmospheric stability conditions on flat terrain without any obstacle to the wind flow. Then, they simulated the wind flow via CFD. They reproduced those experimental results fixing average values for k and ε at the inlet boundary conditions ( Table 2) and specifying values for C 3ε within the computational domain that depends only on the vertical direction. Alinot and Masson 16 expressed the C 3ε obtained as a polynomial expression (Eq. 10) that depends on the height (z). In Eq. (10), L is the Monin-Obukhov length, n is an integer, and a n is the coefficient found by Alinot and Masson 16 . C 1ε and C 2ε remained constant. We adopted this methodology in our work to describe the variations of turbulence with atmospheric stability conditions, and we also considered the k and ε average values reported by them for the inlet boundary condition.
Simulation of the dispersion of pollutants. The transit of vehicles over roads generates pollutants that are emitted from the vehicle's tailpipe and from the road-wheels interaction. Due to the complexity of the nearroad air pollutions, researchers have divided its study into 3 phases 15 : • From the source (tailpipe or wheel-road interaction) to the ambient near the vehicle, where solid and gas phase pollutants maintain their characteristics with a dilution ratio of up to 1000:1. Our work concentrates on the last phase. We, as all studies focused on the third phase, assume uniform concentration on the road as a net result of the multiple vehicles crossing by the same position over a long time, regardless of the exact source of the pollutants (tailpipe or road-wheel interaction). For comparative purposes, we chose an arbitrary emission of E b = 1 g s −1 m −2 and verified that the results do not depend on the emission rate.
Chemical reactions were not considered in the model, which is an acceptable assumption for the case of inert pollutants or pollutants with a mean lifetime much longer than the residence time within the computational domain. For a wind speed of 1 m s −1 , the residence time of the pollutants within the computational domain ( Fig. 1) is ~ 5.5 min.
Based on the experimental work of Huertas et al. 28 , we modeled total suspended particles (TSP) that follow a Rossin-Ramler particle size distribution with a maximum diameter of 34 μm, a minimum diameter of 1 μm, an average diameter of 8.5 μm and a dispersion parameter of 3. We also used the Discrete Random Walk (DRW) model to observe the particle dispersion. This model is a stochastic zeroth-order tracking method that starts with the velocity flow field obtained by solving the Navier-Stokes equations. Then, the DRW integrates Eqs. (11) and (12) to predict particle trajectory using the local continuous phase conditions as the particle moves through the flow, for a sufficient number of representative particles 27 .  www.nature.com/scientificreports/ In Eqs. (11) and (12), u and u p are the fluid and particle velocity vectors, F d is the drag force acting on the particle, g z is the gravity vector acting in the vertical direction, ρ and ρ p are the fluid and particle density, F is any external force acting on the particle, d is the aerodynamic particle diameter, Re is the Reynolds number based on the local speed and particle diameter, µ is the fluid viscosity, and C d is the particle drag coefficient. Bold symbols are vectors. The DRW model does not include particle diffusion effects since said phenomena is non-significant for the size of the particle considered in this study (0.1 < d < 30 µm) 27 . When particles become smaller than the ones considered in this work (ultra-fine particles, d < 100 nm), Eqs. (11) and (12) are modified to consider molecular effects. This DRW model determines the effect of the particles on the velocity flow field of the continuous phase via the source terms of momentum, and iterates, coupling the impact of the discrete and the continuous phases, on the fluid flow. We used the discrete phase module of ANSYS -Fluent v17, to simulate the dispersion of the particles described above. This module has implemented the DRW described above 27 . We used 50.000 time-steps to track the trajectory of each particle within the computational domain. Finally, we established that particles that reach the ground surface get trapped in the short grass that covers the near road surface.
Convergence and grid independence analysis. In this study, the pressure-based solver of ANSYS-Fluent v17.1 was used, with the following options: (i) Double precision, as recommended for large geometries and significant pressure and speed variations 29 . (ii) The steady-state transport equations. (iii) The second-order interpolation scheme. (iv) The semi-implicit method for pressure-linked equations (SIMPLE), which uses a combination of continuity and momentum equations to derive an equation for pressure (or pressure correction).
Quad-type elements were used to discretize the domain shown in Fig. 1. Aa structured mesh was applied because it is easy to implement, requires less computing time, and facilitates particle tracking. We used 3.4 million of computational cells for the simulations with refinement in the regions of interest, such as the sections adjacent to and downwind the roadway. This mesh was selected after a grid independence analysis from 0.11 to 3.9 million elements. On a server with 16 parallel processors, solution times ranged from 30 min to 3 h for the finest mesh.

Use of dimensionless numbers to study the effect of stability conditions on pollutants dispersion.
Wind speed has a strong influence on pollutants dispersion. Aiming to isolate the effects of wind speed from the effects of the other variables that influence the atm(ospheric stability condition, results are reported in terms of the dimensionless numbers developed by Huertas and Prato 17 . They demonstrated that, for neutral atmospheres, when pollutant concentrations downwind the road are expressed in terms of the dimensionless numbers for concentration (C*, Eq. 13) and distance to the road (x*, Eq. 14), the resulting profiles are independent of wind speed (U), emission rates per unit area (E) and the nature of the pollutant through the Schmidt number (Sc). Therefore, variations on the obtained universal profile are attributed to variations on the atmospheric stability conditions rather than just variations on wind speeds.
In Eqs. (13) and (14), w is the road width, and S c is the ratio of momentum diffusivity (kinematic viscosity, υ) and mass diffusivity (D). It is used to characterize fluid flows in which there are simultaneous momentum and mass diffusion-convection processes.
Comparison with experimental results. Aiming to validate the implemented CFD-based model for the dispersion of pollutants near roads (NR-CFD model), results obtained from this model were compared with the experimental data obtained in the Near Road Tracer Study (NRTS08). This study was carried out in 2008 by the NOAA, INL, ARL, and EPA in Idaho Falls, USA 30 and reported by Finn et al. 11 and Steffens et al. 18 . Figure 2 illustrates the experimental setup implemented. They released a 1 g s −1 constant flow of a tracer gas (SF 6 ) distributed uniformly along a 54 m line and used an array of bag samplers to measure the 15 min average SF 6 downwind concentrations at 3,4,6,8,11,15,20, and 30 times the road width. They adopted an equivalent road width of 6 m in their experiments. This gas was chosen due to its negligible background concentration in the atmosphere. They placed anemometers at several vertical positions to measure wind speed, wind direction, turbulence characteristics, friction velocity, heat flux, and atmospheric stability. Tests were repeated under several atmospheric conditions when the wind flowed nearly perpendicular to the line emission source.
Even though their interest was the study of the effects of a side road barrier on the dispersion of pollutants, they carried out 31 tests without any barrier to the wind free flow under different atmospheric stability conditions. We used this subset of data obtained without any barrier for our comparative analysis. A Schmidt number for SF 6 of 0.207 was used to express their experimental results in terms of C* (Eq. 13). www.nature.com/scientificreports/

Results and discussion
As a first step, we ran the NR-CFD model without any emission from the road for the stability conditions considered in this study, observed the evolution of the wind speed and vertical temperature profiles along with the computational domain, and confirmed that they remained essentially unaltered. This fact confirms that the condition of horizontal homogeneity for wind speed and temperature were met. We highlight that this condition of horizontal homogeneity is met when the inlet and ground boundary conditions are congruent with each other and when those boundary conditions satisfy the governing equations of momentum and energy. Then, we added the emissions of solid and gas-phase pollutants emitted from the road and observed their dispersion under different cases of atmospheric stability conditions. For the case of neutral conditions, Fig. 3a,b show the variations on the wind speed and turbulent kinetic energy, respectively, resulting from the mixing of the pollutants emitted from the road and the incoming wind flow. They show that the pollutants emitted at the road disturb the nearby wind flow and that those perturbations tend to disappear as the pollutants move away from the road.
Similarly, Fig. 3c illustrates the downwind TSP concentration resulting from the action of the wind dispersing the pollutants emitted from the road. This figure shows that the pollutant concentrations at the edge of the road reach their maximum values and then reduce toward a stable value far away from the road. Figure 3d shows the downwind ground TSP concentrations as a function of wind speed, keeping the emission rate and the condition of atmospheric stability constant. It shows that the pollutant concentration (in this case, TSP) at the road edge reduces from ~ 15 to ~ 2.5 g m −3 when the wind speed increases from 0.75 to 5 m s −1 . These results demonstrate the strong influence of wind speed on the resulting downwind pollutant concentration. Figure 3e shows the same results but now expressed in terms of C* and x*. It shows that all the profiles of downwind ground concentration shown previously in Fig. 3d now fall into a single universal profile. Previous results indicate that that C*, at any distance downwind the road, is independent of wind speed. It means that according to Eq. (13), the pollutant concentration (C) is inversely proportional to wind speed, which agrees with experimental observation. This result was first reported by Huertas and Prato 17 . Furthermore, they showed that C* is independent of the emission rate and the nature of the pollutant under consideration.
Comparison with experimental results. Figure 4.a shows the results obtained in the NRTS08 campaign when there was no barrier to the wind free flow for the cases when the atmospheric conditions were neutral (− 500 < L < − 100) and slightly stable (20 < L < 100). We used Sc = 0.21 for SF 6 . This figure confirms that, for neutral atmospheric conditions, when the downwind pollutant concentration is expressed in terms of dimensionless numbers, the observed profiles of C* versus x* are the same regardless of the wind speed or mass emission rate.
We then used the NR-CFD model to simulate the dispersion of SF 6 under the same conditions that the experimental tests were carried out, using as input only the data points of speed and temperature measured by the meteorological station at 10 m high. Then, we compared these experimental results of C* with the results of C* obtained by the NR-CFD model for identical positions downwind the road and under similar atmospheric stability conditions (similar L´s). Figure 4.b compares the SF 6 experimental results with the obtained by the NR-CFD model simulating the dispersion of SF 6 , NOx and CO 2 for the case of neutral atmospheric conditions. Vertical bars indicate the range of variation of experimental measurements. We used Sc = 0.946 for NOx and Sc = 0.86 for CO 2 . In agreement with a previous work 17 , these results indicate that the C* versus x* is independent of the gas nature.
We repeated the comparison for results obtained under different atmospheric stability conditions and observed that the C* simulated and C* experimental are highly correlated (R 2 > 0.95) and with slopes close to one www.nature.com/scientificreports/ (Fig. 4c). These results demonstrate the capacity of the NR-CFD model of reproducing short-term experimental measurements of pollutant dispersion under different atmospheric stability conditions.
Effects of the atmospheric stability conditions on pollutant dispersion. Finally, we used the experimental results and the NR-CFD model to study the effects of the atmospheric conditions on the dispersion near roads of solid and gas phase pollutants.
Previously it was shown that when the downwind concentrations are expressed as C* vs x*, the resulting profile is independent of the wind speed, mass emission rate and the nature of the pollutant. That is, the use of the dimensionless numbers for concentration C* isolates the effect of wind speed, on pollutant dispersion. Then, from now on we will express results only in terms of these dimensionless variables. Figures 5a, b show the downwind ground concentration profiles (C* vs. x*) for solid (TSP) and gas-phase pollutants, respectively, obtained under different atmospheric stability conditions. Figures 5a,b show differences in www.nature.com/scientificreports/ the horizontal concentration profiles due to variations of the atmospheric stability conditions. Again, we highlight that these differences are due to the heat transfer processes that occur in the atmosphere under different stability conditions rather than the predominant effect of wind speed since the effects of wind speed are isolated from the analysis when the pollutant concentration is expressed as the dimensionless number for concentration (C*). Aiming to quantify these differences, we selected the results obtained when the atmosphere was under neutral conditions and specifically when L = − 500 as a base case of comparison. Figures 5c shows that the downwind TSP concentration near roads obtained under different conditions of atmospheric stability is highly correlated with the one obtained under the neutral condition of atmospheric stability (R 2 > 0.95). This figure shows that TSP dispersion under extremely stable atmospheric conditions (L = 10) leads to concentrations up to 2 times higher when compared to neutral conditions (slope ~ 2.05). Meanwhile, when the atmosphere is in a slightly stable condition (L = 310), TSP concentrations are ~ 48% greater than when the atmosphere is in a neutral condition. Finally, we observed that under any neutral condition (L < − 100), these slopes remained approximately equal to 1. We named these slopes, obtained from each correlation analysis, as the atmospheric stability dispersion factors (f s ).
Similar results were obtained for the case of gas-phase pollutants. Figure 5d shows the results of the correlation analysis carried out among the gas-phase downwind concentrations obtained under different stability conditions compared to the ones obtained under a neutral condition of reference (L = − 500). For all cases, we observed a high correlation between the C* and the C* considered as reference (R 2 > 0.99). Figure 5f shows the slopes (fs) obtained from the correlation analysis as functions of L. It shows that fs remains constant (fs ~ 1) for all simulations performed under neutral conditions (L < − 100). It increases to fs = 1.6 for the simulations carried out under slightly unstable conditions (− 100 < L < − 10) but decreases to fs ~ 0.66 for the simulations conducted under extremely unstable conditions (− 10 < L < 0). When the atmosphere is under extremely stable conditions Aiming to double-check that these results were consistent with experimental observations, we revisited the subset of experimental data obtained in the NTRS08 campaign for the tests that were carried out without any barrier. We selected the test carried out when the atmosphere was under the most neutral conditions (L = − 500) and used the observed profile C* versus x* of this test as the base case of comparison. Then, for each test, we performed a correlation analysis between the C* measured during each test and the C* selected as the base case (Fig. 5g). As in the NR-CFD results, results showed a high correlation (R 2 > 0.82) between the C* and the C* considered as the most neutral for all cases (Fig. 5e). Then, the fs obtained from the correlation analysis were plotted as function of L (Fig. 5f). This figure shows, qualitatively, that the experimentally obtained fs exhibit the  www.nature.com/scientificreports/ same behavior obtained using the NR-CFD model. Finally, aiming to quantify the level of agreement between the values of fs obtained experimentally and numerically, we carried out a final correlation analysis among them. We found (no shown) a high agreement between them (slope = 1.11) with high consistency (R 2 > 0.95).
Use of the results obtained in this work for practical applications. As described before, for practical applications in environmental science, the results reported by the NR-CFD model are of great interest since they allow the prediction of short (~ 1 h) and long term (1 day-1 year) average concentration of a given pollutant at a given distance (x) from the road, under the varying conditions of wind speed (U), ambient temperature (T), and atmospheric stability conditions measured with a standard meteorological station. Nonetheless, the major drawback in the use of the NR-CFD model is the need of well-trained personnel to set up the NR-CFD model and the need of expensive software and hardware to run it for each time step. We propose, as an approximation, the use of the universal solution reported in Fig. 3e (C* vs x*) and of the atmospheric stability dispersion factors, fs (Fig. 5f), as follow.
• For a given distance to the road (x), calculate the normalized distance to the road (x*) using Eq. (14). • For each time step of data reported by the meteorological station: • Read from Fig. 3e the C* for the x* calculated in the previous step.
• Read from Fig. 5f the value for fs according to the atmospheric stability condition.
• Obtain the desired concentration C (x) using Eq. (15) • Average the values obtained at each time step for C (x) .
When the wind blows in a direction different than the direction of the line perpendicular to the road, we suggest to use the strategy reported in Huertas et al. 17 , which is based on the consideration of the effective distance to the road x e = x/cos θ, where θ is the angle between wind direction and the line perpendicular to the road.

Conclusions
This work aimed to evaluate the influence of the atmospheric stability on the dispersion of pollutants emitted from an area source over a flat surface without any obstacle to the free wind flow.
Initially, a near road CFD based model to simulate the dispersion of pollutants in a reduced computational domain was implemented (NR-CFD model). It: (i) solves via CFD the governing equation that describe the dispersion of pollutants, (ii) uses the Monin-Obukhov Similarity (MOS) theory to specify the vertical profiles for wind speed and temperature as the inlet and ground boundary conditions, and (iii) uses the k-ε standard turbulence model calibrated experimentally for each atmospheric stability condition by Pieterse and Harms. This approach provides a rough approximation to the atmospheric stability conditions. However, it is not appropriate for the study of the atmosphere dynamics but can be used to obtain an approximation to the effects of the atmospheric dynamics on the dispersion of pollutants.
Then, aiming to isolate the effects of wind speed on pollutant dispersion, results were reported in terms of the dimensionless numbers of concentration (C*, Eq. 13) and distance to the road (x*, Eq. 14).
Results from the NR-CFD model agree with the experimental results obtained by the NOAA in 2008 observing the SF 6 concentrations at several distances downwind a source line of SF 6 under different atmospheric stability conditions. Experimental and analytical results show that when pollutant concentrations downwind the road are expressed as C* versus x*, the resulting profile is independent of wind speed, emission rates and the nature of the gas-phase pollutant. Therefore, variations observed on that profile could be attributed to variations on the atmospheric stability conditions rather than the sole effect of wind speed.
Experimental and analytical C* obtained under different atmospheric stability conditions are highly correlated (R 2 > 0.82) to the C* obtained under neutral conditions. We named the slope obtained from the correlation analysis as the atmospheric stability factor (fs). Results show that fs remains constant (fs ~ 1) for neutral conditions (− 500 < L ≤ − 100). It increases to fs = 1.6 for slightly unstable conditions (− 100 < L ≤ − 10) and slightly stable conditions (L > 20) but decreases to fs ~ 0.66 for extremely unstable conditions (− 10 < L < 0). When the atmosphere is under extremely stable conditions (0 < L ≤ 20), fs varies with L and reached values of up to fs = 12.
Finally, for obtaining long-term average concentrations downwind the emission source, we proposed to assume a pseudo-steady-state approximation where at each time step (~ 1 h) the atmosphere is assumed to be under steady-state condition. For each time step, we suggest to use the results obtained under neutral conditions and affect the C* by f s (Fig. 5f) to include the effects of the atmospheric stability conditions on pollutants concentration, where C* continues being the dimensionless number for pollutant concentration near roads obtained under neutral atmospheric conditions. Even though these results indicate that they are valid for any non-reactive solid or gas-phase pollutant, we highlight that they are limited to the case of emission sources located on flat terrain without any obstacle of the free wind flow. www.nature.com/scientificreports/ Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.