Numerical simulation of particle transport in the urban boundary layer with implications for SARS-CoV-2 virion distribution

This paper presents a description, verification, and application of a model simulating the transport of aerosol particles with different properties in the urban boundary layer by using a Lagrangian approach. The model takes input fields of air flow characteristics from arbitrary external models or analytical solutions and allows estimating the movement, sedimentation, and decay of particles. In this paper, the accuracy of the model is successfully estimated on the basis of exact analytical solutions. Simulations are made for a series of urban canyons under different conditions of stratification and wind speed to assess the effects of these meteorological parameters on particle transport in urban areas. Under similar conditions, the transport of particles simulating SARS-CoV-2 coronavirus particles is calculated.


Introduction
This work is devoted to the development, verification, and application in practical problems of the Lagrangian model of transport of aerosol particles in the urban atmospheric boundary layer with a high spatial resolution.
The functioning of modern cities, especially of the most populous ones, is difficult to imagine without vehicles, industry, energy -sources of pollutants in the atmosphere. Such substances include various aerosols -particles suspended in the air. High concentrations of anthropogenic aerosols can lead to effects ranging from light haze to smog, from allergies to cancer. In conditions of high population density and developed international communication in cities, there is a high probability of transmission between people of various infections, many of which are transmitted by airborne droplets -biologically active atmospheric aerosol.
In urban environments, it becomes necessary to monitor particle concentrations with a high spatial resolution, but contact measurement networks have low station density and problems with representativeness, and remote sensing methods do not provide accurate measurements within built-up areas. The study of the topic and the development of the model are motivated by the need to improve the forecast of air quality for health and economy, as well as by lack of research tools, among which physical and mathematical modeling plays a key role. It is important to realize in such a model the Within the framework of this work, the goal is to study the transport of particles in the boundary layer of the atmosphere over urban built-up areas depending on the stratification of the atmosphere and the properties of particles by means of numerical modeling. The main tasks are the development of a model, its verification on analytical solutions under idealized atmospheric conditions, and further application of the model under conditions of a realistic turbulent flow over urban development.

Model description
In the developed model, the interaction of spherical aerosols with an air flow and solid surfaces is reproduced inside a limited three-dimensional area with a given geometry of boundaries and characteristics of the medium. There is a possibility of tracking particle trajectories, calculating concentrations in the air and accumulating on surfaces. The Lagrangian approach is applied to particles, while the characteristics of the air environment are specified by the Eulerian method, i.e. as three-dimensional fields on a discrete grid. These fields and their change over time are specified by the input data. The influence of the transported particles on the thermohydrodynamic state of the air when interacting with it is assumed to be negligible in comparison with external factors, initial and boundary conditions.

Mathematical model of Lagrangian particle transport
To develop the model, the Lagrangian approach to modeling particle transport was chosen. With this approach, each investigated particle is tracked along its trajectory, the velocity and acceleration are calculated for it at each moment of time; this allows one to study both the behavior of an individual particle and the propagation of a large number of particles in terms of concentrations. Thus, the Lagrangian approach is very informative, although it requires a large number of calculations to obtain statistically significant estimates of concentrations.

Governing equations.
The equation of motion is used to calculate the trajectory of a particle suspended in the air. It is based on Newton's second law, which in the model has the following form: (1) (2) where u p is the particle velocity, x p is its coordinate, t is time, g is gravitational acceleration, ρ p is the particle density, ρ is the air (medium) density, u is the ambient flow velocity, and F D is the drag coefficient.
Particle collisions are not taken into account here, since they are rare at not very high concentrations. The first term in (1) corresponds to the buoyancy force acting on the particle, and the second term characterizes the drag force of the medium. When describing the motion of particles in laminar flows, the Stokes law is used to calculate the drag force of the medium; however, this solution is not suitable for the boundary layer of the atmosphere, where the flow is turbulent. For a turbulent flow, F D in (1) can be specified, for example, by the formula [1]: where μ is the dynamic air viscosity, C D is the drag coefficient, Re is the Reynolds number for the particle, and d p is the particle diameter. C D is calculated using empirical formulas for spherical particles from the work [2]: (4) where the values of the coefficients a 1 , a 2 , a 3 depend on the range of values of the Reynolds number.  (3), the effect of stochastic vortices is taken into account through the Reynolds number for a particle calculated as follows: The averaged component is read from the input data on the average flow velocity or is set analytically, while the second one must be calculated separately using turbulent parameterizationsdescriptions of subgrid eddies that are not explicitly resolved in the flow velocity field. Taking into account the effect of turbulence is important when describing the motion of particles in this model, since aerosols have a size much smaller than the grid step of the input data and can stay inside one cell for a long time, falling under the influence of subgrid vortices.
In this model, three parameterizations are implemented: a simple Gaussian one, a discrete random walk model, and a random displacement model. The model allows using any of the presented turbulence parameterizations. In the first two cases, the fluctuation component is a Gaussian random variable ξ with a standard deviation determined by the turbulent kinetic energy (TKE): ( 7) where σξ is the standard deviation, and k is the turbulent kinetic energy.
In the simple Gaussian parameterization, a random variable ξ is generated at each time step of the Lagrangian model, which corresponds to the effect of various subgrid vortices at each iteration of the finite-difference scheme. In the discrete random walk model, the interaction time of a particle with a turbulent vortex t int is introduced, during which the fluctuation component of the velocity acting on the particle remains constant, characterizing the effect of a particular vortex [3]: (8) where t e is the eddy lifetime, and t R is the time for the particle to cross the eddy (transit time scale).
The eddy lifetime t e is estimated as follows: (9) where l e is the characteristic size of a randomly sampled vortex, which is assumed as a dissipation length scale given in [3].
The transit time scale t R is estimated from the following expression: (10) where τ is the droplet "relaxation" time defined in [3]. The random displacement model is similar to the simple Gaussian model in that a random normally distributed variable ξ is generated at each time step of the model. But the standard deviation and the pulsation component have a different form [4]: where ξ i is the random variable for the i-th axis, Δt is the time step, and K s is the turbulent diffusion coefficient.
Finally, turbulent fluctuations also take place in the air density value, , which is included in the buoyancy acceleration; however, they make a negligible contribution to the difference for "heavy" particles (with a significant density), so as .

Particle lifetime implementation.
The model includes the possibility of calculating the transfer of bioaerosols, which, unlike most inorganic particles, have a lifetime, which is different for different types of biologically active particles. Before the introduction of such possibility into the model, particles could be removed from the atmosphere in only two ways: either they flew out of the boundaries of the investigated region, or they settled on solid surfaces. For bioaerosols, the possibility of decay of each particle at each iteration was included.
In the implemented algorithm, the particle lifetime is presented as a half-life which, if necessary, can be calculated from other indicators. The general formula for the number of "alive" particles depending on time is as follows: (12) where N is the number of particles at time t, N 0 is the initial number of particles, and τ is the half-life for these particles.
From formula (12), it is possible to calculate the probability of particle decay P decay at each iteration with a time step Δt: At each time step, for each "alive" particle, a random variable r uniformly distributed over the interval (0, 1) is generated. The particle decays if r satisfies the condition .

Numerical implementation and operation scheme
Various explicit and implicit finite-difference schemes are available for approximate solution of differential equations. In this model, to solve the particle equations of motion (1) and (2), an explicit 1order Euler method is used. The equations are as follows: where the subscript n denotes the value of the variable in the position of the particle at time t n , the subscript n + 1, in the position of the particle at time t n+1 , a is the total acceleration of all forces except for the drag force of the medium, and τ D is the period of the drag force of the medium. The main algorithm of the model implements a time sequence of iterations with a given time step limited by the simulation period. Input data include thermohydrodynamic characteristics of the air (medium) and geometric parameters of a limited three-dimensional area at certain moments in time. The used thermohydrodynamic characteristics of the atmosphere at the moment include: wind speed components along three coordinate axes, turbulent kinetic energy, dissipation rate of turbulent kinetic energy, and turbulent diffusion coefficient. The model uses the computational grid "A" according to the classification of Arakawa: the values of meteorological parameters refer to the centers of the cells.
The boundary conditions for particles in the model in the absence of a solid surface at the boundary of the computational domain can have the form of a transparent boundary when the particles are considered to have escaped from the region (they are not calculated further), or a periodic boundary when a particle escaping from the area is transferred and flies into the area from the opposite side. The periodicity condition means that the simulation area is considered to be infinite horizontally, and all thermohydrodynamic fields and particle trajectories are periodic along the corresponding horizontal coordinate with a period equal to the size of the computational area. When a particle comes into contact with a solid surface, either its adhesion or elastic rebound occurs: the choice depends on the settings for this surface.
The model is implemented in the C++ programming language.

Model verification on analytical solutions
Practical use of the model is impossible without understanding its accuracy in problems with a known solution. Such analytical solutions in terms of concentrations can only be obtained by using the Euler approach. At the same time, it was proved that exact solutions of the equations of the Lagrangian and Eulerian methods are equivalent in terms of the concentration field when using the random displacement model and a time step approaching zero [5]. Thus, at a sufficiently small time step it is possible to compare the concentrations from the Lagrangian model with the analytical solution of the Euler approach.
Analytical solutions for two idealized cases were obtained from the advection and diffusion equation for the Euler approach [4]: where is the Reynolds-averaged concentration, u i is the i-th component of the flow velocity, x i is the i-th coordinate, and Q s is the sum of sources and sinks. It was assumed that the flow is stationary and horizontally homogeneous in all thermohydrodynamic quantities, that the average vertical velocity is zero. At the lower boundary the particle flux from the surface was specified, and at the upper boundary -the condition of zero particle concentration.
Two cases were considered as exact solutions: a boundary layer and a Couette flow over a flat surface. When calculating the transfer of particles in the model, the corresponding vertical profiles of the turbulent diffusion coefficient were set. In the first idealized case -the case of a logarithmic profile in a neutrally stratified boundary layer -the diffusion coefficient was specified as follows: ( 17) where is the friction velocity (set to a constant m/s), κ is the Von Kármán constant, , an z is the vertical coordinate. In the second idealized case, the Couette flow, the diffusion coefficient is given by a more complex formula: where h is the upper boundary height.
In the absence of an average vertical velocity, the diffusion coefficient and its vertical derivative make the main contribution to the nature of particle advection; therefore, the diffusion coefficient profile in each case will determine the shape of the concentration profile.
The calculation results for the case of the boundary layer are shown in Figure 1a. The model (without decay) and analytical profiles overlap each other and have a high correlation coefficient (0.96), small discrepancies are observed only at the boundaries of the computational domain. In general, this indicates an almost complete agreement between the Euler and Lagrangian approaches in the described case. To study the behavior of the profile in the case of bioaerosols, three additional calculations were performed using the Lagrangian model for particles with half-lives of 100, 10, and 1 second, respectively. As can be seen from Figure 1a, differences in half-lives by one order of magnitude lead to the formation of completely different profiles; at similar concentrations near the emission height, a 2-fold difference is observed already at a quarter of the height of the domain, increasing with height.
The results of calculations for the case of a Couette flow are shown in Figure 1b. The model (without decay) and analytical profiles again have a single form and a high correlation coefficient (0.99), small discrepancies are observed only at the lower boundary of the computational domain, where the emission of particles occurs. Similarly to the first experiment, we can talk about almost complete correspondence between the two approaches to calculating the particle concentration in this case. For the Couette flow, three additional calculations were also performed with similar half-lives (100, 10, and 1 second). As can be seen from Figure 1b, the difference between the cases is much ENVIROMIS 2020 IOP Conf. Series: Earth and Environmental Science 611 (2020) 012017 IOP Publishing doi:10.1088/1755-1315/611/1/012017 6 greater, especially without decay and with a 100-second half-life. The reason for this difference is that it takes longer for the particles to propagate vertically to the same heights as compared to the boundary layer case.

Transport of aerosol particles in a series of urban canyons
After verifying the developed model on analytical solutions, one can proceed to the problem of studying the transport of particles in typical urban conditions, depending on the stratification of the atmospheric boundary layer and the characteristics of particles. In the following we show results of particle transport simulations in a series of two-dimensional urban canyons (a three-dimensional area with the same building geometry and flow characteristics along the canyons).
For calculations of thermohydrodynamic fields and particle transport, a three-dimensional region was set with buildings installed parallel to each other, imitating a series of urban canyons. The dimensions of the area were 288 x 60 x 120 meters along the X, Y, and Z axes, respectively. A computational grid was built with a resolution of 1.5 meters. To obtain the input data, the ENVI-met meteorological model is used, which allows calculating the fields of meteorological quantities taking into account the geometry of buildings [6]. When calculating the input fields in the ENVI-met model, periodic boundary conditions were established for temperature, humidity, and turbulence characteristics. In combination with the cyclical configuration of buildings, this makes it possible to obtain at the entrance to the region a flow that has already adapted to this geometry of urban canyons.
For the created area, numerical experiments were carried out to reproduce two typical meteorological scenarios: unstable stratification with a higher wind speed during daytime heating and stable stratification with a lower speed during nighttime cooling. In both cases, the initial conditions were the air temperature at a height of 2 meters equal to 20 °C, relative humidity of 70%, east wind (against the X axis). The wind speed at the right border was given a logarithmic profile with a value at the 10-meter level of 5 m/s for unstable stratification and 1 m/s for stable stratification. Since the ENVI-met model takes into account the change in the components of the radiation balance depending on the time of day, date and latitude, the coordinates of the center of Moscow and the dates 02:00 (for the night scenario) and 14:00 (for the daytime) on June 22 were set.
The emitter of particles (with a diameter of 1 μm) in these experiments is the entire lower surface of the second from the right canyon, which simulates the emission of particles from cars, the rise of dust, and the spread of bioaerosols from plants, people, and animals.

3.2.1.Inorganic particles transport.
The calculation results for unstable stratification are shown in Figure 2a, and for stable stratification in Figure 2b. The largest part of the ejected particles remains in the air of the emitting canyon; however, with unstable stratification the ejected particles spread over the entire height of the region and to the surface of the earth in the following canyons. In both cases, the capture of particles by vertical vortices inside the canyons is noticeable, which is more intense in the case of unstable stratification. Under the influence of such a vortex, large concentration gradients are formed within the canyon, especially between its leeward and windward sides. Based on the results of the described tests, it can be concluded that aerosol pollution of the streets (canyons) adjacent to the emitting canyon is possible with unstable stratification and is most pronounced at the level of the upper floors of buildings from the windward side.

3.2.2.
Bioaerosol particles transport: SARS-CoV-2 case-study. In the light of the global pandemic of SARS-CoV-2 coronavirus, the possibilities of transferring virus particles (virions) outside a canyon where emission occurred have been investigated. The lifetime of virions was set through a half-life of 90 seconds; such values can be observed at sufficiently high temperature and humidity and when the virions are exposed to ultraviolet radiation [7].
The calculation results for different stratifications are shown in Figures 3a and 3b. With unstable stratification, particles in large numbers reach the next canyons, and even their lower parts where people can be. The concentrations in the emitting canyon and the following canyons differ by 2 orders of magnitude or more. With stable stratification, fewer particles leave the initial canyon, and most of them decay above the nearest canyons due to low advection speed. The particles do not reach the lower parts of the following canyons, which indicates almost complete safety in these zones. The results confirm that in conditions of a dense urban area and rather unstable stratification the virus can freely spread to the territory of neighboring streets and courtyards, albeit in small concentrations. This

Conclusions
This paper presents the development and applications of a numerical, physical and mathematical model of particle transport in the boundary layer of the atmosphere over urbanized areas.
This study included a verification of the model on two exact analytical solutions, which showed good agreement. This allows using the thus created program for solving specific research and applied problems.
The model has the possibility of studying bioaerosols with a certain lifetime. This technology has been tested on the example of the SARS-CoV-2 coronavirus in a typical urban environment -a series of canyons. Also, under these conditions the influence of stratification in the atmospheric boundary layer on the transport of particles was studied. The results indicate a high potential for aerosol pollution far from a source located in a canyon and the development of atmospheric instability and spatial heterogeneity of such pollution.
The developed model can be used in the future to solve specific problems and perform full-fledged physical research. It is also possible to integrate the method into large eddy simulation (LES) models for a more accurate description of the processes.