Large-Eddy Simulations of Oil Droplet Aerosol Transport in the Marine Atmospheric Boundary Layer

In this study, a hybrid large-eddy simulation (LES) model is developed and applied to simulate the transport of oil droplet aerosols in wind over progressive water waves. The LES model employs a hybrid spectral and finite difference method for simulating the wind turbulence and a bounded finite-volume method for modeling the oil aerosol transport. Using a wave-following coordinate system and computational grid, the LES model captures the turbulent flow and oil aerosol fields in the region adjacent to the unsteady wave surface. A flat-surface case with prescribed roughness (representing a pure wind-sea) and a wavy-surface case with regular plane progressive 100 m long waves (representing long-crest long-wavelength ocean swells) are considered to illustrate the capability of the LES model and study the effects of long progressive waves on the transport of oil droplet aerosols with four different droplet diameters. The simulation results and statistical analysis reveal enhanced suspension of oil droplets in wind turbulence due to strong disturbance from the long progressive waves. The spatial distribution of the aerosol concentration also exhibits considerable streamwise variations that correlate with the phase of the long progressive waves.


Introduction
In the aftermath of an offshore oil spill accident, the spilled crude oil can stay on the ocean surface for a long period of time before being restored. During this period, frequently-occurring events at the sea surface, such as wave breaking and bubble burst, can generate a considerable amount of small aerosolized oil droplets [1], which can then be further transported upward and downwind by turbulent wind in the marine atmospheric boundary layer (MABL). This process can induce serious threats to public health if the aerosols remain suspended and carried by wind over coastal ocean near high population urban cities. In the past several decades, considerable progress has been made in understanding and modeling the marine aerosol generation and transport in the MABL flows [1][2][3][4][5]. Despite these progresses, the detailed turbulent flow physics in the wind adjacent to the ocean surface is not fully understood due to the challenges associated with modeling the aerosol transport close to unsteady wavy surfaces. It is worth mentioning that several recent studies using a high-fidelity two-phase flow direct numerical simulation method were able to model the detailed small-scale physical process for the generation of sea spray droplets from breaking waves, providing valuable insights for understanding and quantifying the generation of the sea spray aerosol [6,7].
In recent years, large-eddy simulation (LES) has become a valuable tool for modeling turbulent flow and material transport in the ocean and atmosphere. A number of studies have developed numerical algorithms based on a wave surface-fitted coordinate and computational grid system and applied to model the wave-phase-correlated flow phenomena in the near-surface wind field of the MABL [8][9][10][11][12][13]. For the ocean side, LES models combining turbulent flow and particle transport solvers have been developed and successfully applied to study the oceanic particle transport under an idealized flat surface [14][15][16][17][18]. LES has also been applied to simulate oil plume dispersion related to subsea oil spills [19][20][21][22].
In this study, many of these recent developments in the LES model technique are adopted to create a numerical modeling framework for oil droplet aerosol transport in the MABL in a wave-phase-resolved context. In particular, the current LES model uses a wave surface-fitted computational grid system that follows the instantaneous sea surface geometry to capture the turbulent flow transport phenomena near the sea surface. The LES model utilizes different spatial discretization methods for different aspects of the flow problem, i.e., a hybrid pseudo-spectral and finite-difference method for simulating the turbulent wind and a finite-volume method with a combination of an upwind scheme and central difference scheme for simulating the oil aerosol transport. Using this LES model, a set of simulations are performed to study the effects of surface waves and the droplet settling process on the transport of oil droplet aerosol in the MABL. In particular, two different sea surface conditions are considered: one case with a flat surface combined with a prescribed sea surface roughness to model the aerosol transport over pure wind-sea in which the wind-generated short waves are not directly resolved by the LES; and the other case with long wavelength regular plane-progressive waves together with the prescribed roughness to model the effects of ocean swell waves generated by distant storm fronts [8]. For each flow condition, the transport of oil droplets with four different diameters (thus different settling velocities) are simulated. Both instantaneous and statistically-averaged aerosol concentration fields are analyzed to help quantify and understand the effects of surface waves and droplet settling on the vertical and downwind transports of oil aerosols. This paper is organized as follows. Section 2 discusses the model equations and numerical schemes used in LES. Section 3 shows the simulation and statistical analysis results. Finally, conclusions are given in Section 4.

LES Model Formulation
In this study, we simulate the aerosol transport by the turbulent wind in a neutrally-stratified marine atmospheric boundary layer. The coordinate system is denoted as x i (i = 1, 2, 3) = (x, y, z), where x is for the streamwise direction, y is for the spanwise direction, and z is for the vertical direction with the origin z = 0 defined at the mean water surface level. The corresponding velocity components in the x-, y-, and z-directions are denoted as u i (i = 1, 2, 3) = (u, v, w), respectively. In the LES model, wind turbulence is governed by the filtered Navier-Stokes equations, which are written using the index notation as: Here, the tilde indicates the variables resolved at the LES grid scale ∆; ρ a is the density of air; τ ij = u i u j − u i u j is the subgrid-scale (SGS) momentum flux per unit mass, with τ d ij = τ ij − τ kk δ ij being its deviatoric part and δ ij being the Kronecker delta; p * = p + ρ a τ kk /3 − p ∞ is the filtered modified pressure, with p being the resolved dynamic pressure. In this study, we consider the direction of the mean wind being along the x-direction. The pressure gradient dp ∞ /dx imposed in the x-direction models the effect of geostrophic wind forcing that drives the wind field [23]. In Equation (2), the molecular viscous effect is neglected because in LES of the MABL flows, the effects of unresolved SGS terms dominate over the molecular viscous terms [23]. In this study we consider a relatively low oil droplet aerosol release rate, so that the feedback force from the suspended droplets to the wind is negligible.
The transport of oil droplet aerosols in the wind turbulence is simulated using a Eulerian approach. The mass concentration of monodispersed oil droplets of diameter d is described by the continuous Eulerian mass concentration function C d (x i , t). The LES model solves the filtered transport equation of C d [17,20,24], where w s is the droplet aerosol slip (settling) velocity in the air and φ c j = ( u j C d − u j C d ) is the SGS oil droplet aerosol mass concentration flux.
The LES model Equations (2) and (3) are mathematically closed by modeling the SGS momentum flux τ d ij and aerosol concentration flux φ c j with proper models. In particular, τ d ij is modeled using the Lilly-Smagorinsky eddy viscosity-type model [25,26], τ d ij = −2ν τ S ij , where the SGS eddy viscosity is modeled as ν τ = (c s ∆) 2 | S|. Here, S ij = (∂ u i /∂x j + ∂ u j /∂x i )/2 is the resolved strain-rate tensor, | S| = 2 S ij S ij is the strain-rate magnitude, ∆ is the LES grid (filter) scale, and c s is the Smagorinsky model coefficient. In the current LES model, c s is determined dynamically in real time during the simulation using the Lagrangian-averaged scale-dependent (LASD) dynamic model [27], which has been successfully applied in several prior LES studies of the atmospheric boundary layer flows [11,23,28]. The SGS aerosol concentration flux φ c j is parameterized using an eddy diffusivity closure, φ c j = −(ν τ /Sc τ )(∂ C d /∂x j ), where ν τ is provided by the LASD model and the SGS turbulent Schmidt number is set to be Sc τ = 0.4 [20,21,24].
At the particle scale, the aerosolized oil droplets may exhibit various non-spherical shapes and instantaneous deformations, which can be different from those for the natural sea spray droplets. The difference may become more significant if dispersants are applied for oil spill remediation, which can reduce the oil surface tension to enhance droplet breakup [29,30]. These effects are not included in the current LES considering that the coarse grid resolution used in a typical LES of the MABL flow cannot provide sufficient information for directly determining the detailed droplet shapes. For the sake of simplicity, in this study the settling velocities of the aerosolized oil droplets are modeled based on the parameterization for spherical particles that has been widely used for modeling sea spray droplets [1]. In particular, the settling effect of aerosolized oil droplets can be characterized based on the particle Reynolds number Re d = ρ a w s d/µ a , where µ a is the dynamic viscosity of air. At low particle Reynolds numbers Re d < 1, the droplet slip velocity can be accurately modeled using the Stokes law, while at larger particle Reynolds numbers, additional turbulence effects need to be taken into account [1]. Here, the widely-used empirical relationship is adopted to model the droplet slip velocity, which is valid for solid spherical particles with up to Re d ∼ O(10 5 ) [1,31]: where w s,0 is the slip velocity given by the standard Stokes law for low particle Reynolds number condition, where ρ d is the density of the oil droplets, g is the gravitational acceleration, and C f is an empirical correction coefficient as a function of the particle Reynolds number [31], Equation (3) accounts for the main effects acting on the oil droplets, including the Stokes drag, gravitational force, and buoyancy. Note that the droplet response time when transported by the turbulent wind can be estimated as [32]: where R = 3ρ a /(2ρ d + ρ a ) is the density ratio parameter. On the other hand, the smallest time scale resolved by the LES grid scale can be estimated as τ ∆ ∼ | S| −1 . For the MABL flow and oil droplet parameters considered in this study, the corresponding LES grid-scale Stokes number is ). This small Stokes number allows us to neglect other secondary effects (e.g., the SGS fluid stress force, the Saffman lift force, and the droplet inertia) in the current LES for simplification purpose [2,3,20,32].

Wave Surface Boundary Conditions
In this study, we model the effects of long-crest long-wavelength ocean swell waves on the aerosol transport using idealized regular plane-progressive waves [8]. The wave surface elevation and orbital velocities are prescribed based on the Airy wave solution [8,13], Here, η is the wave surface elevation, (u w , v w , w w ) are the wave surface orbital velocities, a is the wave amplitude, k = 2π/λ is the wavenumber, λ is the wavelength, and c is the wave phase speed. For a finite water depth h, the wave dispersion relationship gives c = (g/k) tanh(kh). Equations (8)- (11) provide the wave surface conditions for the LES of the wind turbulence.
In the LES of the marine atmospheric boundary layer flow, the very high Reynolds number prevents the direct resolution of the viscous boundary layer in the air flow adjacent to the water surface. As a result, it is impractical to prescribe the wave surface orbital velocities (u w , v w , w w ) as the Dirichlet boundary conditions to the momentum equation. Instead, an equilibrium surface-layer model is applied to model the SGS stresses of the wind turbulence at the water surface, which provides a Neumann-type boundary condition given by [8,27,33]: where: In Equations (12) and (13), κ = 0.4 is the von Kármán constant, z 0 is the SGS sea surface roughness, ( ...) denotes variables filtered at the test-filter scale 2∆, and ( u r , v r , w r ) are the 2∆-filtered wind velocities relative to the water surface at the first off-surface grid point (where z 1 denotes its instantaneous vertical distance to the local wave surface), Following previous LES studies of the atmospheric boundary layer flows [23,27,28], in Equations (12) and (13), the 2∆-filtered velocity components are used in order to ensure numerical stability when applying the equilibrium logarithmic similarity law-of-the-wall as the instantaneous local boundary condition in LES.
Note that the bottom boundary conditions described above are idealized to simplify the configuration of the LES while capturing the dominant effects of the surface waves on the wind turbulence and aerosol transport. Other effects such as the wave growth/decay, the temporal and spatial variations of the SGS sea surface roughness, and the thin surface layer of drift currents generated by the wind shear are not considered in the current study. As pointed out in previous studies [8,10,13,34], these effects, excluded in the current study, are also of great interest and demand further attention in future research.

Numerical Schemes
Due to the presence of surface waves at the bottom boundary of the wind field, the simulation domain of the LES has an irregular geometry that varies in time. To capture the effects of the surface waves on the wind turbulence and aerosol transport, the current LES model uses a boundary-fitted grid system that follows the instantaneous wave surface motion, as illustrated in Figure 1a. Discretizing and solving the LES governing equations on this moving grid system is challenging. In the current LES model, an algebraic mapping method is used to transform the irregular wave surface-bounded domain in the physical space (x, y, z, t) to a right rectangular prism in the computational space (ξ, ψ, ζ, τ), as illustrated in Figure 1b. In particular, the coordinate transformation is given by: Here, H(x, y, t) is the local instantaneous height of the physical domain, which is decomposed into an average height H 0 and a wave-induced variation −η(x, y, t). The spatial coordinates in the computational space can also be denoted using the index notation as ξ i (i = 1, 2, 3) = (ξ, ψ, ζ). For the algebraic mapping method defined in Equation (20), its Jacobian matrix of the spatial coordinate transformation is: Because the wave surface oscillates in time, a transformation is also required for the temporal derivative terms, i.e., ∂ ∂t where: By applying the chain rule together with Equations (21) and (22), the governing Equations (1) and (2) for the velocity field are transformed into the following form in the computational space [10,34]: where J ij is the index notation of the matrix J. The aerosol concentration transport Equation (3) is transformed into a strong conservative form [8,9] in the computational space to ensure the mass concentration for the oil droplet aerosol, which is written in index notation form as: where |J| is the determinant of the Jacobian matrix.
In the current LES model, the transformed Navier-Stokes Equations (24) and (25) are discretized by a pseudospectral method on a collocated grid in the ξ-and ψ-directions and a second-order central difference method on a staggered grid in the ζ direction. The equations are advanced in time using a fractional-step method. First, Equation (25) is integrated in time by the second-order Adams-Bashforth scheme to obtain a prediction of the new velocity vector field. Then, the divergence-free condition in Equation (24) is achieved by constructing and solving a Poisson equation for a pseudo-pressure field, which is then used to project the predicted velocity onto a divergence-free space. More details of this fractional-step projection method used in the current LES model can be found in [34,35]. Differently, the transformed oil aerosol transport equation is discretized by a finite-volume algorithm with a bounded third-order upwind scheme for the advection term [24,36,37]. A similar LES flow and particle transport modeling strategy has been successfully applied in several recent simulations of oil transport in the ocean boundary layer with simple flat boundaries [19][20][21][22]. In this study, we adopt this finite-volume methodology, but implement it in a more complicated LES model framework with an unsteady deforming boundary and coordinate transformation. This allows us to simulate the oil aerosol transport in the lower portion of the marine atmospheric boundary layer over progressive surface waves. Figure 2 shows a representative simulation result obtained using the current LES model. More details of the simulation parameters are discussed in the next section.

Simulation Setup
In this study, we used a computational domain of size (L x , L y , H 0 ) = (600, 400, 100) m discretized using 128 evenly-distributed grid points in each direction. Similar to many previous LES studies of the MABL flows, we considered the open-sea condition and imposed periodic velocity boundary conditions in the xand y-directions. The top boundary of the simulation domain was kept flat, with free-slip conditions for u and v and impermeable condition for w. The bottom of the simulation domain was bounded by the water surface. Two bottom conditions were considered in the current study. In the pure wind-sea case, the bottom surface was kept flat, and the effect of unresolved sea surface wind-waves was parameterized using the SGS sea surface roughness height z 0 . We used a typical open-sea value of z 0 = 2 × 10 −4 m that has also been used in several previous LES studies [8,9,11,12]. In the plane-progressive (swell) wave case, a two-dimensional downwind wave train with wavelength λ = 100 m and steepness 2πa/λ = 0.1 was imposed on the water surface, with the same SGS roughness z 0 = 2 × 10 −4 m as the pure wind-sea case for the unresolved small-scale waves. For a water depth of h = 30 m, the corresponding wave phase speed was c = 12.2 m/s, and the wave period was T = 8.2 s. Due to the relatively high computational cost associated with the current boundary-fitted LES model, the current study only considered one SGS surface roughness value and one wave condition chosen based on previous studies. Changing the values of z 0 , λ, and 2πa/λ can also affect the aerosol transport due to their effects on the wind turbulence field. These effects are not considered in this paper and may be investigated in future study.
In the simulations, the wind flow was driven by an imposed streamwise pressure gradient dp ∞ /dx as included in Equation (2), which is related to a wind friction velocity as u * = −(dp ∞ /dx)H 0 /ρ a . In this study, we considered a wind friction velocity of u * = 0.37 m/s, which corresponds to a characteristic wind speed of U 10 = (u * /κ) ln(10/z 0 ) = 10 m/s at a 10 m height for the pure wind-sea case. In the plane-progressive wave case, the effective value of U 10 is smaller due to the additional pressure form of the drag effect induced by surface waves. Note that u * was set to be the same in both the flat-surface and plane-progressive wave cases so that the turbulence statistics can be compared. The wind and wave conditions in the progressive wave case correspond to a wave age of c/u * = 33, which falls into the mature sea state [38]. At this wave age, the mean momentum transfer from the wind to the waves was expected to be small, and the distortion of the fast waves to the wind turbulence field was expected to be strong [8,34,38,39]. It is worth mentioning that Sullivan et al. [8] considered a lower wind speed (i.e., a geostrophic wind of 5 m/s over similar 100 m plane-progressive regular waves) and found the existence of the low-level jet (i.e., the increased wind speed in the lower portion of the wind field) as the wind field was driven by the fast waves.
In this study, the air density was set to be ρ a = 1.1845 kg/m 3 , and the air viscosity was µ a = 18.444 × 10 −6 kg/(m s). The oil droplet density was set to be ρ d = 895.5 kg/m 3 [19,20]. Four different droplet diameters were considered, i.e., d = 2.5 µm, 40 µm, 60 µm, and 100 µm. The corresponding slip velocities based on Equations (4)-(6) were w s = 1.6 × 10 −2 cm/s, 3.9 cm/s, 9.1 cm/s, and 21.4 cm/s, respectively. For each oil droplet size, the oil aerosols were released from 29 surface computational cells surrounding the point (x 0 , y 0 ) = (42, 0) m, with a total surface release area of 424.8 m 2 and a mass release rate of 1 g/s. The suspended oil droplets exited the simulation domain by deposition back to the water surface or through the downstream end of the domain via outflux condition.

Instantaneous and Averaged Flow and Aerosol Fields
Figures 3 and 4 illustrate the instantaneous flow fields on the (x, z)-plane across the center of the aerosol releasing region for the flat-surface and plane-progressive wave cases, respectively. The LES results showed two main effects induced by the downwind long waves on the wind field. First, the wave-induced form of drag on the wind field caused the wind velocity to reduce (Figures 3a  and 4a). Second, the fast surface orbital velocity of the swell waves induced strong disturbance to the wind field, resulting in some strong upwelling events (e.g., at x = 380 m in Figure 4b) that may lift the oil aerosol to a considerable elevation.  Figures 5 and 6 show the instantaneous oil aerosol concentrations of the four droplet diameters for the flat-surface and plane-progressive long wave cases, respectively. For both surface conditions, the wind turbulence caused considerable vertical suspension and downwind transport of oil droplets. The droplets with d = 2.5 µm had negligible slip velocity and were transported like passive tracers. For larger oil droplets, the settling effect caused considerable deposition of oil aerosols back to the water, resulting in lower oil aerosol concentration towards downstream locations. For d = 40 µm and d = 60 µm, a large fraction of the oil droplets still remained suspended after being transported through the simulation domain and exited at the outflux boundary at x = 600 m. The concentration of the suspended oil droplets of d = 100 µm was much lower than those of other droplet sizes due to the significant settling effect. Note that the upwelling events (with positive vertical velocity fluctuation) in the turbulent boundary layer played a major role in suspending the oil droplets, while the downwelling events (with negative vertical velocity fluctuation) together with the settling effects helped to deposit the oil droplets back to the ocean. The intermittency in the upwelling/downwelling events together with the strong settling effect resulted in the obvious intermittency in the instantaneous distribution of the large oil droplets with d = 100 µm (Figures 5d and 6d).  It is worth mentioning that for the fast progressive wave case (Figure 6), strong vertical transport of the oil droplets may occur as a result of the upwelling event event caused by the strong wave disturbance (e.g., see the positive vertical velocity fluctuation in Figure 4b and the high oil concentration pocket in Figure 6 at around x = 380 m), which was capable of lifting even the large 100 µm droplets to a high elevation (e.g., 40 m in the snapshot shown in Figure 6d). A similar upwelling event can also be observed in the flat-surface case, but the magnitude appeared to be relatively weak (e.g., the positive vertical velocity fluctuation in Figure 3b and the high oil concentration in Figure 5 at around x = 340 m). The average effects of the plane-progressive wave and slip velocity on the transport of oil aerosol can be obtained by performing the ensemble average of the instantaneous flow field snapshots. In the current analysis, 200 snapshots of the three-dimensional instantaneous flow fields (containing both velocity and oil droplet concentration information) were sampled. These samples were separated by a constant time interval of one wave period (i.e., T = 8.2 s), which were then averaged to obtain the time-averaged quantity denoted by · . For an instantaneous physical quantity f , its turbulent fluctuation is defined as f = f − f . Note that for the plane-progressive wave case, because all the instantaneous snapshots were sampled at the same fixed wave phase (see, e.g., Figure 6), the time average f can also be referred to as the wave phase average [13,39]. Performing the average at the fixed wave phase allowed us to keep the surface wave elevation information in the averaged field, which helped to reveal the effect of waves on the wind turbulence and aerosol transport [12,34].
Using the time average analysis, the enhanced vertical wind motion due to the long plane-progressive waves can be seen clearly from the vertical profiles of the root-mean-square (rms) vertical velocity fluctuation w rms shown in Figure 7. At most of the heights in the simulated domain, the wind turbulence over swell waves exhibited stronger vertical velocity fluctuation, and thus was expected to be able to suspend the oil droplet aerosols to higher elevations under the same pressure gradient forcing condition (i.e., with the same friction velocity of u * = 0.37 m/s).  Figures 8 and 9 show the time-averaged oil aerosol concentrations of the four droplet diameters for the flat-surface and plane-progressive wave conditions, respectively. In all cases, the oil aerosol concentration was diluted towards the downstream direction due to the vertical and lateral turbulent dispersions that caused the aerosol plumes to expand (see Figure 2 for example). In addition, the settling motions further reduced the oil droplet concentration, with the effect stronger for larger droplets. For the flat-surface case, the concentrations of the three smaller droplet sizes (i.e., d = 2.5 µm, 40 µm and 60 µm) exhibited a gradual increase of height for the upper edge of the high concentration region caused by the turbulent entrainment at the upper edge of the oil aerosol plume. The upper edge of the average aerosol plume for d = 100 µm remained almost flat due to the large droplet slip velocity that overcame the upward turbulent entrainment. Differently, Figure 9 showed that the presence of plane-progressive long waves not only enhanced the mean vertical entrainment of oil aerosols (thus resulting in larger vertical extension of the aerosol plumes), but also caused streamwise variations of the averaged concentration field that correlated with the wave phase.
To quantify the downwind development of the oil aerosol plumes over a flat surface and plane-progressive long waves, the center-of-mass height of the plume z c was calculated based on: where C d is the average concentration. Figure 10 shows z c (x) for all the simulated cases. For the flat-surface case without the plan-progressive long waves, z c increased almost linearly with x, with the slope decreasing as the droplet diameter increased. On the other hand, in the plane-progressive wave case, z c exhibited clearly wave-correlated variation along the x-direction for all four droplet diameters. Similar to the flat-surface case, z c decreased monotonically as the droplet diameter d increased. Except for x < 90 m, the value of z c in the plane-progressive wave case was higher than the corresponding value in the flat-surface case for all four droplet diameters considered in this study.  Moreover, it should be pointed out that the instantaneous oil droplet aerosol concentrations shown in Figures 5 and 6 and the time-averaged oil concentrations shown in Figures 8 and 9 were obtained based on a prescribed continuous source release rate of 1 g/s (see Section 3.1 for the simulation setup). For realistic oil spill accidents, the surface oil slicks may release droplet aerosols in an intermittent manner with a release rate different from the one used in this study depending on the wind condition and sea state. Caution should be given when estimating the air quality based on the idealized simulation setup used in this study. Nevertheless, the statistical results reported here still provide some useful insights for understanding the dispersion of the oil droplet aerosols in the MABL. At the very least, one may be able to estimate the order of magnitude of the oil droplet concentration by rescaling the results reported in Figures 8 and 9 based on the actual aerosol release rate under similar wind and sea conditions.

Conclusions
In this study, a hybrid LES model was developed and applied to simulate the oil droplet aerosol transport in the MABL. The model featured a wave surface-fitted grid to help capture the detailed flow physics in the wind adjacent to the wave surface. Aerosol transport over pure wind-sea was modeled by simplifying the sea surface as a flat boundary, on which the effect of unresolved short waves was parameterized as the surface roughness. In another case, the effect of ocean swell waves was modeled by directly imposing the wave elevation and surface orbital velocities as the bottom boundary conditions of the LES. For oil droplets released from a source region on the water surface, the simulation results showed considerable suspension and downwind transport of aerosol by the turbulence in the MABL. When plane-progressive long waves (swells) were present, their highly-organized wave motions induced strong distortion to the lower portion of the MABL and enhanced the vertical turbulence fluctuation, resulting in enhanced aerosol suspension and wave-phase-correlated streamwise variations.
Author Contributions: M.L. performed the simulations, analyzed the results, and contributed to writing the paper. Z.Z. implemented the finite-volume particle transport module into the LES solver. Y.P. and G.V.I. contributed to interpretation of the results and writing of the paper. D.Y. formulated the problem, interpreted the results, and wrote the first draft of the paper.