Pollutant Dispersion Dynamics Under Horizontal Wind Shear Conditions: Insights from Bidimensional Traffic Flow Models

Abstract

Meteorological factors, specifically wind direction and magnitude, influence the dispersion of atmospheric pollutants due to road traffic by affecting their spatial and temporal distribution. In this study, we are interested in the effect of the evolution of horizontal wind components, i.e., in the plane perpendicular to the altitude axis. A two-dimensional numerical model for solving the coupled traffic flow/pollution problem, whose pollutants are generated by vehicles, is developed. The numerical solution of this model is computed via an algorithm combining the characteristics method for temporal discretization with the finite-element method for spatial discretization. The numerical model is validated through a sensitivity study on the diffusion coefficient of road traffic and its impact on traffic density. The distribution of pollutant concentration, computed based on a source generated by traffic density, is presented for a single direction and different magnitudes of the wind velocity (stationary, Gaussian, linearly increasing and decreasing, sudden change over time), taking into account the stretching and tilting of plumes and patterns. The temporal evolution of pollutant concentration at various relevant locations in the domain is studied for two wind velocities (stationary and sudden change). Three regimes were observed for transport pollution depending on time and velocity: nonlinear growth, saturation, and decrease.

1. Introduction

Road traffic is a significant global source of air pollution and greenhouse gas emissions, releasing particulate matter and nitrogen oxides that degrade air quality and pose health risks. These pollutants also damage vegetation and water bodies, contributing to acid rain and smog [1].

In the context of urbanization choices in certain cities like Paris, it is possible to find large trunk roads and main highways very close to houses. It is of the highest importance to characterize all the features of pollutant transport from different traffic situations connected to urban configurations. Several studies are dedicated to simulating different cases using CFD methods. In [2], the authors investigated the spatiotemporal distribution of traffic pollutants at a busy signalized intersection for different intersection geometries. In [3], the authors analyzed the effect of avenue trees on traffic pollutant dispersion in asymmetric street canyons.

However, it is known that meteorological factors, like wind direction, influence the dispersion of traffic pollutants, affecting their spatial and temporal distribution. Understanding these interactions is crucial for developing effective air quality management strategies to mitigate the health and environmental impacts of vehicular emissions [4].

Wind shear, which refers to changes in wind speed and/or direction over a relatively short distance or period of time, significantly affects pollutant dispersion by enhancing or hindering the mixing of air layers.

Previous studies have highlighted that pollutants released into the atmosphere are initially transported by the prevailing wind [5]. However, the presence of wind shear can alter this transport by stretching and tilting pollutant plumes. Moderate wind shear can enhance vertical mixing, diluting pollutants and reducing their concentrations near the source. This can lead to a more uniform distribution of pollutants in the vertical column of the atmosphere. The effects of wind shear on pollutant dispersion are complex and depend on the specific atmospheric conditions. For instance, strong wind shear enhances vertical mixing. Under certain conditions, it can create turbulence, leading to the formation of eddies and irregular mixing patterns. In [6], the authors showed that both the magnitude and direction of wind shear can significantly influence the spread of pollutants. Understanding the effects of wind shear is crucial for predicting air quality, especially in urban areas where pollution sources are concentrated. Accurate predictions can help manage and mitigate the adverse effects of pollution on public health and the environment. This underscores the importance of considering wind shear in atmospheric models to improve the accuracy of pollution dispersion forecasts [6]. The research contributes to a better understanding of the dynamics of pollutant transport and the factors that influence air quality.

There are many studies on horizontal shear, but only a few have studied the effect of horizontal shear on pollution. Actually, many have studied the effect of horizontal shear in the formation of plumes (see, for example, [7]).

Ref. [8] investigated the impact of horizontal transport and vertical mixing on nocturnal ozone pollution in the Pearl River Delta (PRD) region. The authors focused on understanding the mechanisms that influence ozone concentrations during the night. The study revealed that the horizontal transport of air masses and vertical mixing significantly contribute to elevated nocturnal ozone levels. Observational data with numerical modeling were used to analyze these processes. Understanding these interactions is essential for effective air quality management. The results indicated that horizontal transport plays a significant role in nocturnal ozone pollution in the PRD. Air masses originating from industrial and urban areas can elevate ozone levels in downwind regions during the night. The results also suggest that controlling emissions in upwind regions and considering the timing of vertical mixing processes are crucial for mitigating nocturnal ozone pollution.

Ref. [8] recommended further research to explore the impacts of different meteorological conditions and land-use changes on nocturnal ozone pollution. Additionally, improving the resolution of air quality models can enhance the accuracy of predictions and inform better regulatory policies.

Ref. [9] discussed horizontal dispersion parameters, which are crucial for modeling the long-range transport of pollutants in the atmosphere. The study focused on defining and evaluating parameters that influence the horizontal spread of pollutants over extended distances. The study addressed the need for accurate models to predict the dispersion of pollutants over long distances, which is essential for assessing environmental and health impacts. Horizontal dispersion parameters are critical for determining how pollutants spread laterally as they travel through the atmosphere. These parameters influence the width and shape of the pollutant plume. Ref. [9] utilized the Gaussian dispersion model as a basis for representing pollutant spread. This model assumes that pollutants disperse in a normal distribution pattern from the source and provides valuable horizontal dispersion parameters for long-range transport modeling, improving the accuracy of predictions related to the spread of pollutants. Their work provides a more reliable basis for assessing the impacts of atmospheric pollution, thereby aiding in improved environmental management and policy making.

Ref. [10] reviewed CFD modeling of wind fields and pollutant transport in street canyons, which has significantly enhanced our understanding and capability to manage urban air quality. The authors highlighted that CFD is crucial for understanding wind fields and pollutant dispersion in urban street canyons and that combining CFD studies with experiments is also very important. Their approach provided detailed insights into how traffic-related pollutants disperse over time and across different spatial locations, influenced by various environmental factors. Ref. [11] explored the dispersion patterns of traffic-related pollutants using CFD simulations, complemented by wind tunnel experiments and field measurements.

A recent study [4] examined how dynamic traffic and wind conditions impact the effectiveness of green infrastructure in improving air quality, providing valuable insights for more effective urban planning and pollution control strategies. The numerical simulations in this research were conducted using the ENVI-met model.

In numerical simulations of traffic-related pollutant dispersion, a primary challenge is accurately modeling the production of pollutant concentrations over space and time through a source term generated by traffic activity. In this study, pollution levels are simulated dynamically across all locations within the domain using a 2D traffic flow model. Specifically, our work examines the effects of horizontal wind shear on pollutant dispersion patterns from a bidimensional traffic flow. Various mathematical models exist for two-dimensional road traffic flow, and they differ based on the study’s objective, physical model choice, and solution method. Reviews in [12,13] cover several macroscopic models, including two-dimensional traffic models. Extensions of the LWR model are explored in studies like [14,15], while data-driven methods are increasingly common [16,17]. In [18], finite element methods are applied to a two-dimensional model incorporating disorder effects. Modeling traffic on a 2D plane provides valuable insight for regional analysis, especially in congested urban environments. In this study, we have chosen a flat urban topology for the traffic flow model, where the two-dimensional approximation is particularly relevant.

This work represents an initial step in a larger investigation of wind dynamics, with a specific focus here on horizontal wind shear. We concentrate on how wind velocity varies with time and space. Unlike most previous studies that use arbitrarily defined pollution sources, our approach directly links pollution dispersion to emissions generated by a 2D traffic flow model. This method provides a dynamic connection between pollution levels and vehicle movements. To maximize pollution production, we also introduce a traffic obstacle to create congestion. The objective is twofold: to analyze the spatial distribution of pollution patterns at different times and to assess how concentration levels evolve over time at varying distances from the traffic source.

The paper is structured as follows. In Section 2, we formulate a two-dimensional traffic-related air pollution problem coupled with a two-dimensional air pollution diffusion model. In Section 3, we concisely detail the numerical methods and approximations used to solve the coupled space-time state equations. The analysis of simulation results and the discussion for a portion of the French Périphériques of Paris are achieved in Section 4 before concluding.

2. Problem Formulation

In this section, we examine an urban area with a road network dense enough to be approximated as a continuum, meaning vehicles can move at any point x on the 2D network (such as on multi-lane streets). This road network is denoted by $\Omega$, with its boundary represented by $\Gamma$ [in km]. Here, ${\Gamma }_{\mathrm{in}}$ represents the inflow boundary where vehicles enter, ${\Gamma }_{\mathrm{out}}$ the outflow boundary where vehicles exit, and ${\Gamma }_0$ the boundary of inaccessible areas or obstacles. The characteristic length of the road network is denoted by L. Thus, $\Gamma ={\Gamma }_0\cup {\Gamma }_{\mathrm{in}}\cup {\Gamma }_{\mathrm{out}}$, as illustrated in Figure 1.

2.1. Two-Dimensional Traffic Flow Model

In urban areas, the dense traffic can be modelled with two-dimensional continuous and dynamic models. These models, based on a 2D conservation law, represent the traffic density ${\displaystyle \rho (\mathbf{x},t)}$ [in veh/km² as a variable in a 2D plane ${\displaystyle \mathbf{x}=(x,y)\in \Omega }$ depending on the time ${\displaystyle t}$ [in h]. This density can be deduced by solving the advection–diffusion dynamic equation [19,20]. Such an equation is written as follows:

$$\left\{\begin{array}{cc}{\displaystyle \frac{\partial \rho }{\partial t}}+\nabla .\mathbf{F}=0\hfill & {\displaystyle if\mathbf{x}\in \Omega ,t\in (0,T),}\hfill \\ \mathbf{F}=\rho \mathbf{v}-k\nabla \rho \hfill & {\displaystyle if\mathbf{x}\in \Omega ,t\in (0,T),}\hfill \\ \mathbf{v}=v_{max}\left(1-exp\left({\displaystyle \frac{c}{v_{max}}}\left(1-{\displaystyle \frac{{\rho }_{max}}{\rho }}\right)\right)\right){\mathbf{d}}_{\theta }\hfill & {\displaystyle if\mathbf{x}\in \Omega ,t\in (0,T),}\hfill \end{array}\right.$$

(1)

with the following conditions:

$$\left\{\begin{array}{cc}\rho (\mathbf{x},0)={\rho }^0\left(\mathbf{x}\right)\hfill & {\displaystyle if\mathbf{x}\in \Omega ,}\hfill \\ {\displaystyle k{\displaystyle \frac{\partial \rho }{\partial \mathbf{n}}}=f_{in}+\rho \mathbf{v}.\mathbf{n}}\hfill & {\displaystyle if\mathbf{x}\in {\Gamma }_{in},t\in (0,T),}\hfill \\ {\displaystyle k\frac{\partial \rho }{\partial \mathbf{n}}=0}\hfill & {\displaystyle if\mathbf{x}\in \left({\Gamma }_{out}\cup {\Gamma }_0\right),t\in (0,T),}\hfill \end{array}\right.$$

(2)

where ${\displaystyle \mathbf{F}}$ [in veh/km/h] is the traffic flow vector, ${\displaystyle k}$ [in km²/h] is a positive diffusion coefficient representing the permeability of the road network, ${\displaystyle v_{max}}$ is the maximum velocity, ${\displaystyle c}$ [in km/h] is a velocity that determines how rapidly the velocity magnitude decreases with increasing density, ${\displaystyle {\rho }_{max}}$ is the maximum density, ${\displaystyle {\mathbf{d}}_{\theta }(x,y)=(cos\left(\theta (x,y)\right),sin\left(\theta (x,y)\right))}$ is the direction vector with ${\displaystyle \theta }$ the angle between ${\displaystyle {\mathbf{d}}_{\theta }}$ and the x-axis, ${\displaystyle {\rho }^0}$ is the initial density, ${\displaystyle f_{in}(\mathbf{x},t)}$ [in veh/km/h] is the vehicle flow rate at the inflow boundary ${\displaystyle {\Gamma }_{in}}$, $\mathbf{n}$ is the outward unit normal vector to the boundaries of ${\displaystyle \Omega }$ and ${\displaystyle \mathbf{v}}$ [in km/h] is the advection velocity (of motion) of vehicle flow verifying ${\displaystyle \mathbf{v}.\mathbf{n}=0}$ on ${\displaystyle {\Gamma }_0}$, ${\displaystyle \mathbf{v}.\mathbf{n}<0}$ on ${\displaystyle {\Gamma }_{in}}$ and ${\displaystyle \mathbf{v}.\mathbf{n}>0}$ on ${\displaystyle {\Gamma }_{out}}$. The term ${\displaystyle \rho \mathbf{v}}$, in the second equation of (1), is an advective flow describing the displacement of vehicles at their velocity ${\displaystyle \mathbf{v}}$. This velocity is a priori via the last expression of (1) and it varies according to the density of road traffic. Moreover, ${\displaystyle -k\nabla \rho }$ represents a diffusive flow introduced to smooth out the sudden changes in density and velocity between the different regimes. It explains how vehicles adapt their speeds to the surrounding traffic conditions. Thus, the phenomenon of diffusion in urban traffic can be considered compatible with the behavior of users [20]. It is specified that the anisotropic property of the traffic flow is not adapted to the two-dimensional case in particular for multi-lane traffic [21], which is our case.

Notice the substituting expression, the first two equations of (1), and with certain regularity conditions, we obtain:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \rho }{\partial t}+\mathbf{v}.\nabla \rho +\rho \nabla .\mathbf{v}-k\Delta \rho =0}& {\displaystyle \mathrm{in}\Omega \times (0,T).}\end{array}$$

(3)

The constant ${\displaystyle k}$ defines the level of flow variation as a function of the density gradient.

This equation is then used for numerical approximation (see Section 3). The constant k determines the extent of flow variation as a function of the density gradient; if density increases, it tends to oppose the flow, reducing its magnitude. Several traffic factors, such as vehicle type, weather conditions, and driver behavior, affect the value of k, highlighting its importance in the two-dimensional road traffic model [20]. If k is significantly large compared to the flow speed, diffusion dominates Equation (3); otherwise, advection prevails [22]. The parameter c is crucial in setting the critical density ${\rho }_c$, the density threshold above which traffic jams begin to form. The direction $d\theta$ is fixed in specific areas (e.g., the inflow boundary ${\Gamma }_{\mathrm{in}}$) and depends on the network’s geometry. The last equation of (1) represents the concave Macroscopic Fundamental Diagram (MFD) proposed by Newell [23] and Francklin [24], while the second equation in (2) indicates a heterogeneous Neumann condition at ${\Gamma }_{\mathrm{in}}$. Finally, the last equation in (2) reflects vehicle exit at ${\Gamma }_{\mathrm{out}}$; if a vehicle reaches this boundary, the simulation stops (a homogeneous Neumann condition). Notice that when $k=0$, the model reverts to the well-known two-dimensional Lighthill–Whitham–Richards traffic flow model introduced by Mollier [25].

2.2. Macroscopic Pollution Model

Consider now a larger region, of characteristic length ${\displaystyle L}$ and width ${\displaystyle W}$, containing the road network ${\displaystyle \Omega }$ in order to visualize the distribution of the pollutant outside of ${\displaystyle \Omega }$. This region and its boundary are denoted by ${\displaystyle {\Omega }_p}$ and ${\displaystyle {\Gamma }_p}$ [in km], respectively (see Figure 2).

In this section, we combine the model of traffic flow with a macroscopic air pollution model whose pollution source is generated by vehicles. The source of pollution ${\displaystyle Q}$ [in kg/km²/h] on $\Omega$ depends mainly on the density ${\displaystyle \rho }$, which will be determined by solving the problem (1)–(2) as well as the traffic flow norm ${\displaystyle \parallel \mathbf{F}\parallel }$ computed via the second equation of (1). This source can be defined as the Radon measure (roughly speaking, the measure with respect to the two-dimensional Dirac delta) given by [26] (see Section 4.1.5):

$$\begin{array}{cccc}\hfill Q(.,t):& \mathcal{C}\left(\overline{\Omega }\right)\hfill & \hfill ⟶& \mathbb{R}\hfill \\ \hfill & v\hfill & \hfill ⟼& \langle Q(.,t),v\rangle ={\mathit{\int }}_{\Omega }\left(\eta \rho (\mathbf{x},t)+\nu \parallel \mathbf{F}\parallel \right)v\left(\mathbf{x}\right)d\Omega ,\hfill \end{array}$$

(4)

for each ${\displaystyle t\in (0,T)}$, where ${\displaystyle \mathcal{C}\left(\overline{\Omega }\right)}$ is the space of continuous functions equipped with the usual scalar product ${\displaystyle \langle ,\rangle }$, and ${\displaystyle \eta }$ [in kg/veh/h] and ${\displaystyle \nu }$ [in kg/veh/km] are weight parameters representing contamination rates.

To simulate the distribution pollutant, a mathematical model similar to the one proposed in Skiba [27] is used. The evolution of the pollutant concentration ${\displaystyle \varphi (\mathbf{x},t)}$ [in kg/km²], in the domain ${\displaystyle {\Omega }_p}$ and over the time interval ${\displaystyle (0,T)}$, can be estimated by solving a two-dimensional convection–diffusion pollution model. This model, with initial and boundary values, is written as follows:

$$\left\{\begin{array}{cc}{\displaystyle {\displaystyle \frac{\partial \varphi }{\partial t}}+\mathbf{u}.\nabla \varphi +\sigma \varphi -\mu \Delta \varphi =Q}\hfill & {\displaystyle if\mathbf{x}\in {\Omega }_p,t\in (0,T),}\hfill \\ {\displaystyle \varphi (\mathbf{x},0)={\varphi }^0(x,y)}\hfill & {\displaystyle if\mathbf{x}\in {\Omega }_p,}\hfill \\ {\displaystyle \mu {\displaystyle \frac{\partial \varphi }{\partial \mathbf{n}}}=\varphi \mathbf{u}.\mathbf{n}}\hfill & {\displaystyle if\mathbf{x}\in S^{-},}\hfill \\ {\displaystyle \mu {\displaystyle \frac{\partial \varphi }{\partial \mathbf{n}}}=0}\hfill & {\displaystyle if\mathbf{x}\in S^{+},}\hfill \end{array}\right.$$

(5)

where the field ${\displaystyle \mathbf{u}(\mathbf{x},t)}$ [in km/h] denotes the wind velocity (experimentally known and verifies ${\displaystyle \nabla .\mathbf{u}=0}$), ${\displaystyle \sigma >0}$ [in h⁻¹] is the extinction rate corresponding to the reaction term, ${\displaystyle \mu >0}$ [in km²/h] is the molecular diffusion coefficient [28], ${\displaystyle {\varphi }^0\left(\mathbf{x}\right)}$ is a function giving the initial concentration,

$${\displaystyle S^{-}=\{(\mathbf{x},t)\in {\Gamma }_p\times (0,T)suchthat\mathbf{u}.\mathbf{n}<0\}}$$

represents the inflow of ${\displaystyle {\Gamma }_p}$, that is, pollution flow is directed out of the region ${\displaystyle {\Omega }_p}$ and

$${\displaystyle S^{+}=\{(\mathbf{x},t)\in {\Gamma }_p\times (0,T)suchthat\mathbf{u}.\mathbf{n}\ge 0\}}$$

is the outflow of this boundary, that is, pollution flow is directed into ${\displaystyle {\Omega }_p}$, with ${\displaystyle {\Gamma }_p=S^{-}\cup S^{+}}$ (see Figure 2).

Equation (5) shows that the combined pollutant flux (diffusive and advective) is always zero at the inflow boundary ${\displaystyle S^{-}}$, as there is no pollution source outside ${\displaystyle \Omega }$. This condition means that at the boundary ${\displaystyle S^{+}}$, the diffusion pollutant flux ${\displaystyle \mu \frac{\partial \varphi }{\partial \mathbf{n}}}$ is negligible compared with the advection pollutant flux ${\displaystyle \varphi \mathbf{u}.\mathbf{n}}$ directed out of the domain ${\displaystyle {\Omega }_p}$ [29,30].

3. Numerical Method

This section is devoted to the numerical method used to solve the coupled space-time Equations (1), (2) and (5). An algorithm combining the characteristics method for the temporal discretization with the ${\displaystyle {\mathbb{P}}^1}$ Lagrange–Galerkin finite element method for the spatial discretization is used [31,32]. This algorithm presents good theoretical convergence results, guarantees the stability of the numerical solution and allows us to theoretically avoid the standard CFL (Courant–Friedrichs–Levy) condition on the time step [33,34].

In order to approach the time derivative, we subdivide the time interval ${\displaystyle (0,T)}$ into ${\displaystyle N\in \mathbb{N}}$ subintervals of length ${\displaystyle \Delta t=T/N}$ small enough. Subsequently, we define ${\displaystyle t^n=n\Delta t}$ and denote by ${\displaystyle {\rho }^n}$ the density of traffic at time ${\displaystyle t^n}$ and ${\displaystyle {\varphi }^n}$ the pollutant concentration, for ${\displaystyle n=0,\dots ,N}$. Applying the characteristics method detailed in [35] (see Section 4.2), the total derivatives of ${\displaystyle \rho }$ and ${\displaystyle \varphi }$ in Equation (3) and first equation of (5), respectively, at the instant ${\displaystyle t^{n+1}}$ can be approximated by:

$${\displaystyle \frac{D\rho }{Dt}=\frac{\partial \rho }{\partial t}+\mathbf{v}.\nabla \rho \approx \frac{{\rho }^{n+1}\left(\mathbf{x}\right)-{\rho }^n\circ X^n\left(\mathbf{x}\right)}{\Delta t},}$$

(6)

$${\displaystyle \frac{D\varphi }{Dt}=\frac{\partial \varphi }{\partial t}+\mathbf{u}.\nabla \varphi \approx \frac{{\varphi }^{n+1}\left(\mathbf{x}\right)-{\varphi }^n\circ X^n\left(\mathbf{x}\right)}{\Delta t},}$$

(7)

where ${\displaystyle X^n\left(\mathbf{x}\right)=X(\mathbf{x},t^{n+1};t^n)}$, i.e., the position of the particle (vehicle or pollutant) at instant ${\displaystyle t^n}$ that was in ${\displaystyle \mathbf{x}}$ at the instant ${\displaystyle t^{n+1}}$. These approximations allow us to write a semi-discretization of problems (1), (2) and (5). Given ${\displaystyle {\varphi }^0}$ and ${\displaystyle {\rho }^0}$, we must find functions ${\displaystyle {\varphi }^{n+1}}$ and ${\displaystyle {\rho }^{n+1}}$, for ${\displaystyle n=0,\dots ,N-1}$, respectively, verifying:

$$\left\{\begin{array}{ccc}& {\displaystyle \frac{{\rho }^{n+1}-{\rho }^n\circ X^n}{\Delta t}}+{\rho }^{n+1}(\nabla .{\mathbf{v}}^n)-k\Delta {\rho }^{n+1}=0\hfill & \hfill {\displaystyle \mathrm{in}\Omega ,}\\ & {\mathbf{v}}^n=v_{max}\left(1-exp\left({\displaystyle \frac{c}{v_{max}}}\left(1-{\displaystyle \frac{{\rho }_{max}}{{\rho }^n}}\right)\right)\right){\mathbf{d}}_{\theta }\hfill & \hfill {\displaystyle \mathrm{in}\Omega ,}\\ & k{\displaystyle \frac{\partial {\rho }^{n+1}}{\partial \mathbf{n}}}=f_{in}^{n+1}+{\rho }^{n+1}{\mathbf{v}}^n.\mathbf{n}\hfill & \hfill {\displaystyle \mathrm{on}{\Gamma }_{in},}\\ & k{\displaystyle \frac{\partial {\rho }^{n+1}}{\partial \mathbf{n}}}=0\hfill & \hfill {\displaystyle \mathrm{on}{\Gamma }_{out}.}\end{array}\right.$$

(8)

and

$$\left\{\begin{array}{ccc}& {\displaystyle \frac{{\varphi }^{n+1}-{\varphi }^n\circ X^n}{\Delta t}}+\sigma {\varphi }^{n+1}-\mu \Delta {\varphi }^{n+1}=Q^{n+1}\hfill & \hfill {\displaystyle \mathrm{in}{\Omega }_p,}\\ & \mu {\displaystyle \frac{\partial {\varphi }^{n+1}}{\partial \mathbf{n}}}={\varphi }^{n+1}{\mathbf{u}}^{n+1}.\mathbf{n}\hfill & \hfill {\displaystyle \mathrm{on}{S}^{-},}\\ & \mu {\displaystyle \frac{\partial {\varphi }^{n+1}}{\partial \mathbf{n}}}=0\hfill & \hfill {\displaystyle \mathrm{on}{S}^{+},}\end{array}\right.$$

(9)

where an exponent ${\displaystyle n\in 0,\dots ,N}$ is the corresponding function valued at ${\displaystyle t^n}$.

To complete the semi-discretization, we consider a triangular discretization of the computational domain ${\displaystyle {\Omega }_p}$ (containing ${\displaystyle \Omega }$) [36] and two finite element spaces: ${\displaystyle U_h}$ for density (on ${\displaystyle \Omega }$) and ${\displaystyle V_h}$ for concentration (on ${\displaystyle {\Omega }_p}$). Let ${\displaystyle {\rho }_h^0\in U_h}$ and ${\varphi }_h^0\in V_h$, the variational formulations of (8) and (9) written as follows: find ${\displaystyle {\rho }_h^{n+1}\in U_h}$ and ${\varphi }_h^{n+1}\in V_h$, for ${\displaystyle n=0,\dots ,N-1}$, such that

$$\begin{array}{cc}& {\int }_{\Omega }{\displaystyle \frac{{\rho }_h^{n+1}-{\rho }_h^n\circ X_h^n}{\Delta t}}{u}_{h}d\Omega +{\int }_{\Omega }{\rho }_h^{n+1}\left(\nabla .{\mathbf{v}}_h^n\right)u_{h}d\Omega +{\int }_{\Omega }k\nabla {\rho }_h^{n+1}.\nabla u_{h}d\Omega \hfill \\ & -{\int }_{{\Gamma }_{in}}\left(f_{in}^{n+1}+{\rho }_h^{n+1}{\mathbf{v}}_h^n.\mathbf{n}\right)u_{h}d{\Gamma }_{in}=0,\forall u_h\in U_h,\hfill \end{array}$$

(10)

and

$$\begin{array}{cc}& {\int }_{{\Omega }_P}{\displaystyle \frac{{\varphi }_h^{n+1}-{\varphi }_h^n\circ X_h^n}{\Delta t}}{v}_{h}d{\Omega }_P+{\int }_{{\Omega }_P}\sigma {\varphi }_h^{n+1}{v}_{h}d{\Omega }_P+{\int }_{{\Omega }_P}\mu \nabla {\varphi }_h^{n+1}.\nabla v_{h}d{\Omega }_P\hfill \\ & -{\int }_{S^{-}}{\varphi }_h^{n+1}{\mathbf{u}}^{n+1}.\mathbf{n}{v}_{h}d{S}^{-}={\int }_{{\Omega }_P}\left(\eta {\rho }_h^{n+1}+\nu \parallel \mathbf{F}\parallel \right)v_{h}d{\Omega }_P,\forall v_h\in V_h,\hfill \end{array}$$

(11)

where ${\displaystyle X_h^n}$ and ${\displaystyle {\mathbf{v}}_h^n}$ are approximations of ${\displaystyle X^n}$ and ${\displaystyle {\mathbf{v}}^n}$, respectively, and ${\displaystyle u_h}$ and ${\displaystyle v_h}$ are the test functions in the finite element spaces ${\displaystyle U_h}$ and ${\displaystyle V_h}$, respectively. For more details on the numerical method, see [37]. The numerical solutions ${\displaystyle {\rho }_h^{n+1}}$ and ${\displaystyle {\varphi }_h^{n+1}}$ are computed on each node of the dense mesh and at each time ${\displaystyle t^n}$ with open-source software FreeFem++ (Version 4.7) [38].

4. Numerical Results and Discussion

4.1. Parameterization

The chosen road network ${\displaystyle \Omega }$ (road + bifurcation) is a portion of the northern ring road of Paris (“Périphérique de Paris”) with a characteristic length ${\displaystyle L}$ = 1 km (see Figure 3).

The inflow boundary ${\displaystyle {\Gamma }_{in}}$ is constituted of four road lanes with a width of ${\displaystyle 14\times {10}^{-3}}$ km, i.e., ${\displaystyle 3.5\times {10}^{-3}}$ km per lane. ${\displaystyle {\Gamma }_{out}}$ measures ${\displaystyle 14\times {10}^{-3}}$ km for the exit at the top and ${\displaystyle 7\times {10}^{-3}}$ km for the exit at the bottom, i.e., two lanes. The rectangular obstacle is located in the middle of the network at 0.4 km of ${\displaystyle {\Gamma }_{in}}$ and of dimension ${\displaystyle (12.5\times {10}^{-3})\times (9\times {10}^{-3})}$ km² (see Figure 1). The mesh of ${\displaystyle \Omega }$ is shown in Figure 4.

The modelling period is 8 a.m. to 9 a.m., i.e., ${\displaystyle t\in (0,1)}$ with ${\displaystyle T=1}$ h. In order to have a good precision of the numerical solution, we choose a sufficiently small time step ${\displaystyle \Delta t=0.001}$ h. Initially, there is no traffic in ${\displaystyle \Omega }$, i.e., ${\displaystyle {\rho }^0=0}$ veh/km² and ${\displaystyle \parallel \mathbf{v}\left({\rho }^0\right)\parallel =0}$ km/h. The direction of road traffic in the last equation of (1) is fixed in ${\displaystyle {\Gamma }_{in}}$ per ${\displaystyle {\mathbf{d}}_0=(1,0)}$, i.e., for ${\displaystyle \theta =0}$, and is then given by the geometry of the network. In the second equation of (2), the vehicle flow rate function ${\displaystyle f_{in}}$ is given by ${\displaystyle f_{in}(x,t)=f_{max}\times f\left(t\right)}$, where ${\displaystyle f_{max}=9408}$ [in veh/km/h] is the maximum flow rate and $f\left(t\right)$ is a positive continuous function that controls the variation of the flow rate during the modeling period and is defined by:

$$\begin{array}{c}\hfill {\displaystyle f\left(t\right)=\left\{\begin{array}{cc}5t\hfill & \mathrm{if}t\in \left(0\mathrm{h},0.2\mathrm{h}\right)\hfill \\ 1\hfill & \mathrm{if}t\in \left(0.2\mathrm{h},0.4\mathrm{h}\right)\hfill \\ -3(t-0.6)+0.4\hfill & \mathrm{if}t\in \left(0.4\mathrm{h},0.6\mathrm{h}\right)\hfill \\ 0.4\hfill & \mathrm{if}t\in \left(0.6\mathrm{h},1\mathrm{h}\right).\hfill \end{array}\right.}\end{array}$$

(12)

The value of ${\displaystyle f_{max}}$ is chosen to generate a sufficiently dense traffic after ${\displaystyle 0.4}$ h. In the second equation of (1), the maximum velocity ${\displaystyle v_{max}=70}$ km/h, the parameter ${\displaystyle c=40}$ km/h and the maximum density ${\displaystyle {\rho }_{max}=800}$ veh/km². We solve the problem (1)–(2) for different values of ${\displaystyle k}$ (see Section 4.2). For more details on the parameter settings, we refer the reader to our article [37].

Concerning the atmospheric pollutant generated by road traffic, we are interested in the highly toxic gas of carbon monoxide (CO). It is the most abundant air pollutant in the lower atmosphere, except for carbon dioxide (CO₂) [39]. To visualize the distribution of the CO concentration outside of ${\displaystyle \Omega }$, we consider a larger region ${\displaystyle {\Omega }_p}$ of characteristic length ${\displaystyle L}$ = 1 km and width ${\displaystyle W=0.5}$ km (see Figure 2). Initially, we assume that there is no air pollutant in ${\displaystyle {\Omega }_p}$, i.e., ${\displaystyle {\varphi }^0(x,y)=0}$ kg/km², and simulate traffic-related emissions over a period of 1 h. In the first equation of (5), the CO extinction rate ${\displaystyle \sigma =0.6\times {10}^{-2}}$ h⁻¹ and the CO molecular diffusion coefficient ${\displaystyle \mu =3.5\times {10}^{-8}}$ km²/h. To compute the source ${\displaystyle Q}$ of pollution related to CO, the parameters of contamination rates in Equation (4) are fixed at ${\displaystyle \eta =3.16\times {10}^{-5}}$ kg/veh/h and ${\displaystyle \nu ={10}^{-6}}$ kg/veh/km [40]. The mesh of ${\displaystyle {\Omega }_p}$ is shown in Figure 5.

4.2. Traffic Flow Model Validation

This section aims to validate the proposed mathematical model by examining the effects of three different values of k on traffic density across distinct dynamic regimes. The selected road segment, denoted as $\Omega$, includes a rectangular obstacle to simulate congestion. Here, k represents the effective diffusion coefficient of the traffic flow, with values of $k=1$, $k=1.5$, and $k=3$ km²/h. Figure 6, Figure 7 and Figure 8 illustrate the traffic density $\rho$ over time t. After $0.4$ h, a high density is observed for $k=1$, especially before the obstacle, while a medium density occurs for $k=1.5$, and a low density for $k=3$. This highlights the obstacle’s effect on flow. By ${\displaystyle t=0.5}$ h, density begins to decrease across all regimes, stabilizing (indicating no congestion) at $t=0.6$ h. The impact of k on traffic density and flow is clear, with pollutant transport affected accordingly, as pollution source Q depends on these parameters. Observing density distribution alongside pollutant concentration patterns reveals that wind direction guides pattern evolution, while temporal characteristics are clearly defined. For the remainder of this study, k is set to 1.

4.3. Numerical Results for Different Space/Time Wind Velocity Function

The aim of this academic study is to emphasize the impact of time and spatial dependence of wind on pollution pattern evolution. For spatial dependence, we use a Gaussian wind model. This model provides a simplified, probabilistic approach to representing average wind conditions and minor fluctuations, making it valuable for preliminary analysis. It is particularly effective for short-term forecasting in stable, localized urban environments.

However, urban wind variability, driven by obstacles, turbulence, and asymmetries, often results in non-Gaussian behavior, limiting the accuracy of the Gaussian model, especially in extreme events and complex flow patterns. To achieve more precise results in such environments, alternative models like the Weibull or log-normal distributions are commonly used. For further details, see [41], where the authors discuss the challenges of modeling wind in urban areas using CFD, highlighting the limitations of simplified models like Gaussian distributions in turbulent, obstacle-rich environments and suggesting alternative approaches.

For the time dependence model, we use a linear increase (and decrease) time evolution to simulate a smooth and regular wind to study the pollutant dispersion.

4.3.1. Stationary Wind

Before showing the results for wind velocity as a function of time or space, we have to give a frame with a homogeneous, unidirectional, laminar and stationary wind (${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =2}$ km/h) to have in mind the effect of pollutant transport in this case. During all the simulations, we keep the same direction for the wind. We do not study the effect of a change of direction. Figure 9 shows the distribution of the CO concentration, computed via the source ${\displaystyle Q}$ generated by the traffic density presented in Figure 6, for ${\displaystyle -\pi /4}$ direction wind at different times ${\displaystyle t}$. The concentration is high at ${\displaystyle t=0.4}$ h and it decreases over time. This is due to a high density at this instant, which also decreases over time.

4.3.2. Gaussian Model

The chosen road network ${\displaystyle \Omega }$ (road + bifurcation) is a portion of the northern ring road of Paris (“Périphérique de Paris”). This system is situated between two rows of buildings. That is why we can assume that the maximum of velocity magnitude is in the center of the domain. We name Gaussian model the mathematical expression: ${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =2.}$$exp(-x^2-y^2)$ km/h. Figure 10 shows the time evolution of the concentration distribution for ${\displaystyle -\pi /4}$ direction Gaussian model wind with velocity magnitude at different times ${\displaystyle t}$.

If we compare Figure 9 and Figure 10, we can notice lots of differences for the time evolution patterns. The timescales are very different, with a short time evolution for the Gaussian model. The pollutant accumulation occupies various regions within the domain. This study highlights the significance of wind shear characteristics in the evolution of pollutant patterns. However, after one hour, the distribution pattern appears quite similar.

4.3.3. Time Dependence Wind Velocity Function

We study the time dependence of the wind velocity function magnitude while keeping the same direction. We study three linear cases of interest:

• Case 1: Linear increase in wind velocity magnitude (see Equation (13) for mathematical formulation).
• Case 2: Linear decrease in wind velocity magnitude (see Equation (14) for mathematical formulation).
• Case 3: “Sudden change” in wind velocity magnitude (see Equation (15) for mathematical formulation).

Specifically, for case 3, we want to the study the impact of the formation of pollutant patterns in response to a sudden change in wind speed over a short duration.

$$\begin{array}{c}\hfill {\displaystyle g\left(t\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}t\in \left(0\mathrm{h},0.2\mathrm{h}\right)\hfill \\ 2\hfill & \mathrm{if}t\in \left(0.2\mathrm{h},0.4\mathrm{h}\right)\hfill \\ 3\hfill & \mathrm{if}t\in \left(0.4\mathrm{h},0.6\mathrm{h}\right)\hfill \\ 4\hfill & \mathrm{if}t\in \left(0.6\mathrm{h},0.8\mathrm{h}\right)\hfill \\ 5\hfill & \mathrm{if}t\in \left(0.8\mathrm{h},1\mathrm{h}\right).\hfill \end{array}\right.}\end{array}$$

(13)

$$\begin{array}{c}\hfill {\displaystyle h\left(t\right)=\left\{\begin{array}{cc}5\hfill & \mathrm{if}t\in \left(0\mathrm{h},0.2\mathrm{h}\right)\hfill \\ 4\hfill & \mathrm{if}t\in \left(0.2\mathrm{h},0.4\mathrm{h}\right)\hfill \\ 3\hfill & \mathrm{if}t\in \left(0.4\mathrm{h},0.6\mathrm{h}\right)\hfill \\ 2\hfill & \mathrm{if}t\in \left(0.6\mathrm{h},0.8\mathrm{h}\right)\hfill \\ 1\hfill & \mathrm{if}t\in \left(0.8\mathrm{h},1\mathrm{h}\right).\hfill \end{array}\right.}\end{array}$$

(14)

$$\begin{array}{c}\hfill {\displaystyle i\left(t\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}t\in \left(0\mathrm{h},0.4\mathrm{h}\right)\hfill \\ 4\hfill & \mathrm{if}t\in \left(0.4\mathrm{h},0.6\mathrm{h}\right)\hfill \\ 1\hfill & \mathrm{if}t\in \left(0.6\mathrm{h},1\mathrm{h}\right).\hfill \end{array}\right.}\end{array}$$

(15)

Figure 11, Figure 12 and Figure 13 show the evolution of the pollutant concentration patterns for each case described above. Specifically, we analyze the concentration distribution for wind coming from the ${\displaystyle -\pi /4}$ direction, with increasing step functions of velocity magnitude ${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =g\left(t\right)}$, ${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =h\left(t\right)}$, and ${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =i\left(t\right)}$ at different times ${\displaystyle t}$.

We can directly compare Figure 11 and Figure 12, as the increase and decrease in the pollutant patterns evolve in a “parallel” manner. Notably, there are significant differences in both the spatial and temporal evolution of the patterns. The “accumulated pollutant” patterns (indicated in red in the figures) do not occupy the same areas of the traffic flow, reflecting a purely nonlinear evolution. The timescales for pattern development also differ. However, the same underlying process is observed, where pollutant accumulation grows from the obstacle in the traffic flow. While intuitive, it is important to highlight this point.

In Figure 13, we emphasize the formation of pollutant patterns in response to a sudden change in wind speed over a short period. This case is particularly dramatic in terms of pollutant spread across the domain over time. The higher pollutant concentrations grow from the obstacle (and extend around it) and accumulate rapidly throughout the southern part of the road.

In general, we observe that wind shear significantly influences the horizontal transport of pollutants by stretching and tilting pollutant plumes. This can alter both the direction and the distance over which pollutants are transported from their source. Under unstable atmospheric conditions, strong wind shear can lead to the formation of eddies, resulting in irregular mixing patterns and potentially higher pollutant concentrations in localized regions.

To further study the evolution of the pollution plume patterns, we examine the temporal evolution of pollutant concentration at several relevant locations within the domain.

4.3.4. Time Evolution Pollutant Concentration

To study the evolution of the CO concentration as a function of time, we chose different points at different distances below the road network ${\displaystyle \Omega }$ where there is a pollutant transported from the road by the wind. Figure 14 shows the position of the points ${\displaystyle A_i}$, ${\displaystyle B_i}$, ${\displaystyle C_i}$ and ${\displaystyle D_i}$ for ${\displaystyle i=1,2,3,4}$ in the region ${\displaystyle {\Omega }_p}$.

The letters ${\displaystyle A}$, ${\displaystyle B}$, ${\displaystyle C}$ and ${\displaystyle D}$ are located at a distance of ${\displaystyle 200\times {10}^{-3}}$, ${\displaystyle 400\times {10}^{-3}}$, ${\displaystyle 600\times {10}^{-3}}$ and ${\displaystyle 800\times {10}^{-3}}$ km, respectively, with respect to the left boundary of ${\displaystyle {\Omega }_p}$. The numbers ${\displaystyle 1}$, ${\displaystyle 2}$, ${\displaystyle 3}$ and ${\displaystyle 4}$ are located at a distance of ${\displaystyle 20\times {10}^{-3}}$, ${\displaystyle 30\times {10}^{-3}}$, ${\displaystyle 50\times {10}^{-3}}$ and ${\displaystyle 100\times {10}^{-3}}$ km, respectively, at the bottom of ${\displaystyle \Omega }$. For example, the point ${\displaystyle A_1}$ represents the intersection of the column ${\displaystyle A}$ (at ${\displaystyle 200\times {10}^{-3}}$ km from ${\displaystyle {\Omega }_p}$) and the line ${\displaystyle 1}$ (at ${\displaystyle 20\times {10}^{-3}}$ km from the road ${\displaystyle \Omega }$). The evolution of the concentration is presented at each point as a function of time ${\displaystyle t}$ in Figure 15 for a stationary wind ${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =2}$ km/h (a, b, c, d) and a time-dependent wind velocity magnitude ${\displaystyle \parallel \mathbf{u}(\mathbf{x},t)\parallel =i\left(t\right)}$ (b, d, f, h), respectively.

For each graphic, we can separate the curves into three different regions because of the feature similarity curves. For the stationary wind graphics (a, c, e, g), the regions are named and classified as:

• for $t\in \left(0\mathrm{h},0.2\mathrm{h}\right)$ nonlinear growth: Although the nonlinear evolutions are similar, the shape of the evolution depends on the distance from the Peripherique. The evolution remains steady due to the laminar and unidirectional regularity of the wind and its low speed.
• for $t\in \left(0.2\mathrm{h},0.6\mathrm{h}\right)$ saturation: In this region, a plateau is observed where the concentration no longer changes, indicating a stationary pollutant distribution dynamic. As one moves further from the road traffic, the characteristic time decreases.
• for $t\in \left(0.6\mathrm{h},1\mathrm{h}\right)$ decrease: The last part of the graph reflects the fact that the pollutant exits the area with a characteristic time that strongly depends on the distance from the obstacle, particularly downstream of the congestion.

For the abrupt wind change graphics (b, d, f, h), the regions are named and classified:

• for $t\in \left(0\mathrm{h},0.2\mathrm{h}\right)$ nonlinear growth: We observe the same trends as in Region I, which seems intuitive given that, during this time interval, the wind is steady and of low intensity.
• for $t\in \left(0.2\mathrm{h},0.6\mathrm{h}\right)$ rapid decrease: During the time interval associated with this region, there is a sudden change in wind speed amplitude. In terms of dynamic behavior, the response of the pollutant concentration is an almost instantaneous drop, indicating that the pollutant is quickly expelled from the area, as expected. However, after this immediate drop, three distinct time scales can be observed (except at the position farthest from the obstacle on the right), reflecting a ‘stepped’ decrease in pollutant concentration. This decline follows a nonlinear process.
• for $t\in \left(0.6\mathrm{h},1\mathrm{h}\right)$ increase: During this interval, the wind becomes laminar again, and the dynamics of the pollutant concentration increase up to a certain value. This characteristic growth time depends on the distance from the vertical relative to the Peripherique.

5. Conclusions

This study addresses the problem of horizontal wind shear and the time dependence of wind velocity magnitude. Few studies have focused on the effects of horizontal wind shear, as most research has been dedicated to vertical shear. In this work, we develop a coupled model (two-dimensional traffic flow coupled with pollutant concentration diffusion) where the wind velocity magnitude is an input that can vary. Our study specifically focuses on the time and scale dependence of wind velocity magnitude. Starting with a laminar, stationary, unidirectional wind, we can easily compare different scenarios. We confirm that wind shear significantly affects the horizontal transport of pollutants by stretching and tilting plumes and patterns, altering both the direction and distance pollutants travel from the source (in this case, the obstacle in the traffic flow). In unstable atmospheric conditions, strong wind shear can create eddies, leading to irregular mixing patterns and potentially higher pollutant concentrations in localized areas.

Furthermore, the pollutant concentration distribution, driven by traffic density as the source, is analyzed across a single wind direction with various wind magnitudes: stationary, Gaussian, linearly increasing or decreasing, and sudden shifts over time. The model takes into account the stretching and tilting of pollution plumes. Temporal changes in pollutant concentration are further examined at key locations within the domain under two wind scenarios: stationary wind and sudden wind shifts. For each case, three regimes are observed for pollutant transport, connected to different time scales: nonlinear growth, saturation, and decrease. The study of the time evolution of pollutant concentration confirms the nonlinear features behind pollutant diffusion. However, these time scales are not directly related to the characteristic time evolution of sudden wind changes. Furthermore, in the saturation regime, when studying the time evolution of saturation as a function of distance from the obstacle (i.e., comparing $A_1$ with $A_4$, $B_1$ with $B_4$, $C_1$ with $C_4$, and $D_1$ with $D_4$), we observe a similar reduction in time scale (approximately 20).

Future studies concern the change in the wind direction, the evolution of the patterns for a wind velocity magnitude viewed as a function of space and time (in the same time). We will take a coupled system with a wind given by turbulence models (RANS, LES). Then, we will be interested in data measuring wind, coming from Modeled Wind Resource Data or an experimental database. We conducted simulations that provide us with the concentration values in the different areas near the roads. As future perspectives for this work, we can include the comparison with data collected in real-life situations by sensors.

Future studies will focus on the impact of wind direction changes and the evolution of pollutant patterns with wind velocity magnitude as a function of both space and time. We plan to use a coupled system with wind data from turbulence models (RANS, LES). Additionally, we will explore data from measured wind, sourced from Modeled Wind Resource Data or experimental databases. Our simulations have provided concentration values for different areas near roadways. As a future direction, we aim to compare these results with real-world data collected from sensors.

Our results also reveal a pattern of pollutant accumulation at specific locations away from direct traffic flow, suggesting that green walls could play a critical role in these areas. As natural air filters, green walls absorb pollutants such as CO₂, volatile organic compounds, and particulate matter, thereby improving air quality in urban environments. Additionally, our framework indicates that a targeted road signage policy could help manage pollution dispersion more effectively during traffic disruptions, ultimately reducing the broader environmental impact.