Dispersion Modeling of Particulate Matter in Different Size Ranges Releasing from a Biosolids Applied Agricultural Field Using Computational Fluid Dynamics

This paper proposes a methodology using computational fluid dynamics (CFD)-FLUENT to simulate the dispersion of particulate matter releasing from a biosolid applied agricultural field and predict the particulate concentrations for different ranges of particle sizes. The discrete phase model (Lagrangian-Eulerian approach) was used in combination with each of the four turbulence models: Standard kε (kε), Realizable kε (Rkε), Standard kω (kω), and Shear-stress transport k-ω (SST) to predict particulate matter size concentrations for distances downwind of the agricultural field. In this modeling approach, particulates were simulated as discrete phase and air as continuous phase. The predicted particulate matter concentrations were compared statistically with their corresponding field study observations to evaluate the performance of turbulence models. The statistical analysis concluded that among four turbulence models, the discrete phase model when used with Rkε performed the best in predicting particulate matter concentrations for low (u < 2 m/s) and medium (2 < u < 5 m/s) wind speeds. For high (u > 5 m/s) wind speeds, Rkε, kω, and SST showed similar performances. The discrete phase model using Rkε performed very well and modeled the best concentrations for the particle sizes (μm): 0.23, 0.3, 0.4, 0.5, 0.65, 0.8, 1, 1.6, 2, 3, 4, and 5. For particle sizes: 7.5 and 10, the performances of Rkε, kε, kω, and SST were similar.

science and engineering. For air pollution studies, researchers have applied CFD to simulate the dispersion of contaminants in complex urban environments [1]- [9]. A few CFD studies also involved modeling of the dispersion of particulate matter emissions. Alizadehdakhel et al. [10] studied the dispersion of pollutants from a stack of a cement plant and used the Large Eddy Simulation (LES) method for modeling the particulate matter concentrations for different downwind distances. Pospisil and Jicha [11] used kε RNG turbulence model for simulating the dispersion of PM 10 emissions from moving vehicles along roadways in an urban area.
A literature review showed that the researchers have applied CFD for modeling particulate matter emitting from area-source configurations. Diego et al. [12] applied k-ε model to estimate the downwind concentration of dust due to emissions from an open storage pile and studied the effect of orientation of piles with respect to wind speed on the total dust eroded. Silvester et al. [13] applied realizable k-ε model using discrete phase method to particles releasing during rock blasting operations and estimated their mass percentages depositing within and outside of open pit quarry boundary for different wind directions. Seo et al. [14] estimated the distance of the damage area over an agricultural field due to fugitive dust emissions from a reclaimed land by modeling the dust emissions using k-ε model. Hill et al. [15] studied near-field dispersion of pollutants from an industrial area using CFD and compared the modeled concentrations with those ones measured from wind tunnel experiments. In the study by Bhat et al. [16], dispersion of the particulate matter emitting from a biosolid applied agricultural land was simulated over an open field using k-ε model and the modeled concentrations were statistically compared with the observed ones. None of the above studies did not apply the discrete phase method for modeling particulate matter concentrations.
The primary focus of previous studies involving particulate matter emissions from biosolid applied locations included developing analytical dispersion models and applying them for predicting particulate matter concentrations at downwind distances [17] [18] [19]. However, none of them involved numerical solutions for modeling the particulate concentrations of different size ranges. To the authors' knowledge, no study in the literature involved numerical modeling applying discrete phase in combination with each of the turbulence models for modeling the dispersion of particulate matter emissions from a biosolid applied agricultural field. In addition, the authors could not find studies showing the statistical comparison of the modeling performance between turbulence models for different wind speeds and particle size ranges. Taking into account these research gaps, the present study applies the discrete phase model (Lagrangian-Eulerian approach) using CFD software-FLUENT for predicting the concentrations of particulate matter for different size ranges by simulating particulate emissions from a biosolid applied agricultural field (modeled as an area-source).
Four different turbulence schemes: Standard k-ε, Realizable k-ε, Standard k-ω, and Shear-stress transport (SST) k-ω were used for modeling atmospheric turbulence. The discrete phase model was combined with each of the four turbulence models to predict different size particulate matter concentrations for near-field downwind distances. A methodology for simulation of particulate matter dispersion using the discrete phase method with each of the four turbulence schemes was developed. The discrete phase model accounted for the effects of atmospheric turbulence and drag force which is dependent on particle physical characteristics (diameter, density, and velocity), gravitational velocity, and air viscosity for predicting the trajectories of particles. The modeled particulate matter concentrations were compared statistically with field study observations to evaluate the modeling performance of the turbulence schemes.

Turbulence Schemes and Dispersion Model
Four different turbulence schemes were used in this study. A more detailed explanation and additional corresponding equations can be found in FLUENT manual [20]. A brief description of each model is given below:

Standard k-ε
Based on Reynolds-Averaged Navier-Stokes (RANS) equations, Launder and Spalding [21] proposed one of the simplest turbulence models commonly known as Standard k-ε (hereafter known as kε). The governing equations of continuity and momentum are given by Equations (1) and (2) respectively: It is a semi-empirical model consisting of separate transport equations for calculating turbulence kinetic energy (k) and dissipation rate (ε) as shown in Equations (3) and (4) respectively:

Realizable k-ε
One of the variants of the k-ε model is the Realizable k-ε (hereafter known as Rkε) which has a new transport equation for the dissipation rate and a new formulation for the turbulent viscosity. The transport equations for calculation of the turbulence kinetic energy (k) and the dissipation rate (ε) involved in the Rk-ε model are given by Equations (5) and (6) The model constants are C 1ε = 1.44, C 2 = 1.9, σ k = 1.0, σ ε = 1.2.

Standard k-ω
The standard k-ω model (hereafter known as kω) is an empirical model based on transport equations for turbulence kinetic energy (k) and specific dissipation rate (ω). These models incorporate effects of low Reynolds number, compressibility, and shear flow spreading and are thus applicable to wall-bounded flows and free shear flows. The transport equations for k and ω are given by Equations (7) and (8):

Shear-Stress Transport k-ω
The Shear-stress transport k-ω model (hereafter known as SST) incorporates robustness of both kω and kε models by multiplying them with a function and adding them together. This gives effectiveness by activating only kω model in the near-wall region and only kε models in the far field regions. Also, SST model accounts for the transport of the turbulent shear stress. k and ω are given by Equations (9) and (10):

Dispersion Model: Discrete Phase Model
In this study, particulate matter concentrations were modeled using the discrete phase model built in the FLUENT software. This method involves defining a discrete phase and a continuous phase. Therefore, particulate matter emitting from the biosolids applied agricultural field was defined as the discrete phase and air as the continuous phase. phase on the particles. The trajectory of a discrete phase particle is predicted by integrating the force balance on the particle. This force balance equates the particle inertia with the forces acting on the particle and is given by Equation (11): where, u is the continuous phase velocity, u p is the particle velocity, µ is the molecular viscosity of the fluid, ρ is the fluid density, ρ p is the density of the particle, d p is the particle diameter, and g is the acceleration due to gravity. F D is drag force per unit mass and is given by Equation (12) for sub-micron particles:

Design
The GAMBIT program was used to design the geometry of the modeling domain. The first step in the design of geometry involves the creation of a ground surface and an agricultural field. A rectangular face having dimensions 300 m × 150 m (Length × Breadth) along X and Y axes respectively was designed to represent the ground surface ( Figure 1). Another face with four edges having dimensions 100 m × 100 m representing the agricultural field was created on the ground surface which was designed in step one ( Figure 1). In the second step, the agricultural field was separated from the ground surface. This separation process created the source inlet for particulate matter to enter the geometry domain. In GAMBIT, the "subtract" function was used for the separation process and geometry of the agricultural field was subtracted from the ground surface ( Figure 2). This separation process also allowed assigning different types of Figure 1. Geometry of ground surface (bigger rectangle) and agricultural field (smaller rectangle).
boundary conditions for the ground surface and inlet (see boundary condition discussion). From this step onwards, the source inlet represents the ground-level agricultural field. The third step involved the process of designing the volume of the modeling domain. Using the geometry in step two, the volume of the modeling domain was designed with a height of 50 m. The completed geometry of the modeling domain had following dimensions: Length (L) = 300 m; Breadth (B) = 150 m; and Height (H) = 50 m ( Figure 3).

Boundary Conditions
In the fourth step, types of zones for the modeling domain were specified using the "zones" function in GAMBIT. These zone-type specifications define the physical and operational characteristics of the model at its boundaries and within specific regions of its domain. There are two classes of zone-type specifications: (a) Boundary types (b) Continuum types. Boundary-type specifications define the characteristics of the model at its external or internal boundaries. Continuum-type specifications define the characteristics of the model within specified regions of its domain. The rear-face of the modeling domain near the agricultural field was chosen as velocity inlet through which wind magnitude entered the domain in a direction perpendicular to the plane of the respective face i.e. along X-axis ( Figure 4). The variation of wind speed with height above the ground was incorporated using power law profile given by Equation (13): The exponent "m" depends on atmospheric stability condition. "u 1 " is the measured wind speed at a height "z 1 " and "u z " is the calculated wind speed at a

Data
The data used in this study was taken from a field work by Akbar-Khanzadeh et al. [22] and Bhat et al. [18]. In the field work, the particulate matter concentrations Advances in Chemical Engineering and Science were measured at 10 m and 20 m distances downwind of a biosolid applied agricultural field. The GRIMM Series 1.108 Aerosol Spectrometer instrument was used for measuring particulate matter concentrations. This instrument has an accuracy of ± 5% for measured concentrations and ≤2% for particle size. Measurements were taken on five different days and they consisted of one-houraverage particulate matter concentrations for different particle sizes (in μm): 0.23, 0.3, 0.4, 0.5, 0.65, 0.8, 1.0, 1.6, 2.0, 3.0, 4.0, 5.0, 7.5, and 10. There are 28 one-hour-average particulate matter concentration observations for each particle size. The atmospheric stability condition for all the observations was unstable (either B or C), except for one where it was neutral (D). The wind speed (u) measured at the study site was in different categories: low (u < 2 m/s), moderate (2 m/s < u < 5 m/s), and high (u > 5 m/s). The maximum and minimum wind speed inputs were 0.13 m/s and 8.93 m/s respectively. The emission rates were calculated for different particle sizes using their concentrations measured inside the field and volume of the mixing box over the field. A more detailed explanation on the calculation of emission rates is given in a study by Bhat and Kumar [23].

Methodology
Particulates entered the modeling domain through the inlet (i.e. agricultural field) in a direction perpendicular to the plane of the inlet along Y-axis. At the same time, wind entered the domain through the velocity inlet in a direction perpendicular to the plane of the inlet along X-axis ( Figure 5). According to the FLUENT manual, a turbulent intensity of 1% or less is generally considered low and greater than 10% is considered high. Therefore, a turbulent intensity of 2% and 7% were assumed for the source inlet and velocity inlet respectively. Green-Gauss node solver was used to improve the computational accuracy of results [24]. The methodology developed for setting up discrete phase model, turbulence models, boundary conditions, and other inputs in FLUENT was described in Figure 6.

Statistical Evaluation
The statistical measures were used for evaluating the performance of turbulence P. Nimmatoori, A. Kumar Advances in Chemical Engineering and Science

Grid Sensitivity Analysis
The grid size used for meshing the modeling domain was determined after performing the grid sensitivity analysis. The size of grid was selected in such a way that there is no quantitative difference between modeled concentrations between the two grid sizes and the numerical solution is stable in nature. Therefore, trail runs were conducted using different grid sizes for all turbulence models. Figure   7 shows the modeled particulate matter concentrations by Rkε and kω for different grid sizes and they were unstable for the grid sizes less than 4 m and stable Advances in Chemical Engineering and Science for sizes greater than 4 m. The statistical analysis of modeled concentrations for the grid sizes (m): 4, 4.5, and 5 were performed using t-test with a null hypothesis that there is no significant difference between modeled concentrations if t df < t critical = 1.70 (p < 0.05) and the calculated values of t df did not reject the null hypothesis (Table 1). Therefore, a value of 4.5 m (between 4 and 5) was chosen as the size of grid. The modeling domain was meshed with "submap" using the hexahedral type of meshes which are computationally more efficient [27]. The meshing process created a total of 19,520 hexahedral cells within domain volume.

Statistical Evaluation
The discrete phase modeling method was combined with each of the four turbu-    The performance of each of four turbulence models: kε, Rkε, kω, and SST with the discrete phase model for predicting the concentrations of particulate matter emissions from a ground-level area source were compared by conducting both qualitative and quantitative analysis. The qualitative analysis involved plotting the observed concentrations and their corresponding modeled ones in different ways and visually analyzing the performance of each model. The qualitative analysis results were supported by quantitative analysis which involved determining statistical performance measure criteria discussed in the statistical evaluation section.
The scatter residual (defined as the ratio of model predicted concentration to observed concentration) plots were also used for the qualitative analysis. The performance of turbulence models was compared by conducting scatter residual analysis as a function of important model independent variables such as wind speed, emission rate, and particle size. The turbulence models over-estimated the particulate matter concentrations for low wind speeds. Since the transport of particulate matter in the atmosphere is highly affected by the driving force of wind, a low wind speed condition would lead to small turbulence scales and consequently the spreading rate or dispersion of particles in atmosphere will be less. In addition, low wind speeds tend to decrease the friction velocity and decelerate the transport of particulate matter leading to high concentration levels in the downwind direction near the source location [28]. Figure 8 shows the phenomenon of particulate matter transportation in atmosphere and their settling on the ground for low and high wind speed conditions for approximately equal emission rates. At a low wind speed condition, the cloud representing the high concentration levels of particulate matter dispersed for a short distance in the atmosphere and settled on the ground surface near the agricultural field (Figure 8(a) and Figure 8(b)). From the contour plots, the levels of particulate matter concentrations were high for the downwind distances up to50 m from the edge of the agricultural field. On the other hand, the distribution of particulate matter concentration in atmosphere as well as on the ground surface was seen for longer downwind distances at a high wind speed condition (Figure 8(c) and Figure 8(d)). Figure 9 shows scatter residual plots for different wind speeds (u) which included low (u < 2 m/s), moderate (2 m/s < u < 5 m/s), and high (u > 5 m/s). The dotted lines represent factor of two (FA2) criteria. From Figure 9, it can be noticed that the turbulence models over-estimated the particulate matter concentrations for u < 2 m/s with most of the scatter residual points falling outside FA2. The results of the statistical analysis also support these findings. The statistical measures were calculated with respect to low, medium, and high wind speed conditions after taking in account the modeled particulate concentrations of the turbulence models irrespective of the particle sizes. Table 2 shows the results of statistical analysis for different wind speeds. Table 3 presents the interpretation of statistical criteria which were satisfied (denoted by "") by the turbulence models for respective wind speeds.    (Table 2 and Table 3). Hence, the statistical results clearly indicate that Rkε was the best performing turbulence model for u < 2 m/s. Medium wind speeds (2 < u < 5 m/s): It was determined that more than 65% of modeled particulate matter concentrations were within the factor of two and hence the turbulence models satisfied the criteria of FA2. Again, it was found that Rkε modeled the best particulate concentrations (Table 2). However, the satisfied statistical criteria indicated similar performances of Rkε and kε when the statistical measure of |FB| was not considered (Table 3). Therefore, 95% confidence limits (CLs) were determined on all statistical measures in order to further examine whether Rkε was performing significantly different from kε. The CLs on all statistical measures did not overlap zero (both limits are positive or negative) for Rkε-kε which indicated that Rkε was performing significantly different from kε and hence Rkε was the best performing turbulence model for 2 < u < 5 m/s ( Table 3). The CLs were calculated for Rkε-kω and Rkε-SST in order to show that Rkε was performing significantly differently from kω and SST and it was the best performing among all the turbulence models.
High wind speeds (u >5 m/s): For u > 5 m/s, Rkε satisfied the criteria of statistical measures FA2, NMSE, VG, and MG ( Table 2 and Table 3). When the statistical measure MG was not considered, then kε, kω, and SST showed performance similar to Rkε. Hence, the CLs were calculated to further examine differences between the turbulence model performances. The results of CLs indicated Rkε performing significantly differently from kε and therefore Rkε was the best among them. However, the CLs showed no differences between the performances of Rkε, kω, and SST (Table 3). Therefore, the statistical analyses concluded that the performances of Rkε, kω, and SST were similar for high wind speeds and these results indicate the need for more data points in future research studies. From Table 3, it can be noted that the values of statistical measures of Rkε, kω, and SST were relatively better than those for kε and they also support the results of CLs showing similar performances of those models. Figure 10 shows the turbulence kinetic energy (TKE) and dissipation rate (DR) predicted by each turbulence model for a particle at low and high wind speed conditions. The size of turbulent eddies in the atmosphere decreases with dissipation of their kinetic energy. The turbulent eddies exist in the atmosphere until their kinetic energy is completely dissipated. Since the dispersion of particulate matter in the atmosphere takes place through eddies, they remain in the air as long as eddies possess the kinetic energy. Due to the dissipation of eddy energy, the particulate matter in the atmosphere loses energy and begins to settle on the ground surface. As DR increases, the rate at which the particulate matter settles on the ground surface also increases leading to high concentration levels.
From Figure 10, it can be noticed that kε, kω, and SST had high DR when compared to that of Rkε for low wind speeds. Consequently, kε, kω, and SST modeled high particulate matter concentrations than Rkε for low wind speeds and they did not satisfy the criteria of statistical measures. On the other hand, there were no big differences between the DR of the turbulence models for high wind Advances in Chemical Engineering and Science speeds and hence they showed similar performance which was also indicated by the statistical analysis ( Table 2).
The scatter residual plots of modeled concentrations as a function of particle size are shown in Figure 11. The turbulence models predicted relatively better concentrations for particle sizes larger than 1 μm. In the case of particle sizes 3, 4, 5, 7.5, and 10, their modeled concentrations were much higher for few observations which were represented as unfilled points in Figure 11. Further analysis on the unfilled points was conducted by comparing their corresponding measured concentrations with respect to the size of particles. The measured concentrations showed a drastic decrease in the range of 50% to 90% between sizes 2 µm and 3 µm. Again, they showed a decrease in the range of 50% to 70% between sizes 5μm and 7.5 μm. However, their corresponding calculated emission rates did not show a decrease in their magnitudes accordingly. Furthermore, the wind speed conditions of the measured concentrations corresponding to the unfilled points were low (u < 2 m/s). Consequently, the particle size concentrations predicted by the turbulence models were significantly higher for these unfilled residuals. The unfilled residual points had some effect on the values of statistical measures for particle sizes 3, 4, 5, 7.5, and 10. The results of the statistical measures and their criteria satisfaction (denoted by "") by the turbulence models with respect to particle size showed Rkε performing consistently better than kε, kω, and SST for every particle size ( Table 4). The statistical criterion of FA2 was satisfied by Rkε for all particle sizes (except for 0.4 μm) and other turbulence models satisfied it for the following sizes: kε, 0.65, 0.8, 1.0, 1.6, 2, 4, 5, 7.5, and Advances in Chemical Engineering and Science 10; kω 1.6, 2, 3, 7.5, and 10; SST 0.65, 1.6, 3, 4, 5, 7.5, and 10 (Table 4). Among the turbulence models, Rkε had satisfied the criterion of |FB| for maximum number of particle sizes: 0.8, 1, 1.6, 2, 3, 4, 5, 7.5, and 10 indicating that it is less biased in predicting particle size concentrations.
The results of statistical criteria satisfaction of particle sizes, except for 7.5 and 10, indicated that Rkε and kω were the best and least performing turbulence models. In case of particle sizes 7.5 and 10, the modeled concentrations satisfied all the statistical measures criteria (except MG) and hence the turbulence models showed similar performances. Therefore, further analysis was conducted by determining the confidence limits on statistical measures between the model pairs to examine if one of the turbulence models was performing significantly different from another for the respective particle sizes. The results of confidence limits showed significant differences in the performances of both kε and kω from SST (kε-SST and kω-SST in Table 5). Though the confidence limits on few statistical measures did show a significant difference in the performance of Rkε from kε, kω, and SST (Rkε-kε, Rkε-kω, and Rkε-SST in Table 5) as well as kε from kω (kε-kω in Table 5), it cannot be concluded that one specific turbulence model was the best performing for 7.5 and 10. Based on the results of statistical analysis and confidence limits, it can be concluded that all turbulence models  had similar performance for particle sizes 7.5 and 10 which indicate the need for more data points for the statistical analyses. The performances of the turbulence models were compared with the particulate matter deposition (PMD) model [19] which is an analytical area-source dispersion model. The PMD incorporates dry deposition while predicting concentrations for different particle sizes. The statistical performance measures for the PMD model are also shown in Table 4.
Among four turbulence models, the performance of Rkε was comparable to the PMD model.
The modeled concentrations of turbulence models could not satisfy statistical measures criteria given Kadiyala and Kumar [29] and this indicates the need for larger database consisting of a greater number of observations which can be ob-Advances in Chemical Engineering and Science tained through the future field data collection studies.

Conclusions
The following are the conclusions of this study: 1) The discrete phase modeling method was used successfully with each of the four turbulence models: kε, Rkε, kω, and SST for predicting particle trajectories and their concentrations for distances downwind of a biosolid applied agricultural field. Therefore, it can be concluded that the modeling domain and the dispersion methodology developed in the present study can be applied to future research studies involving the modeling of particulate matter emissions from similar area-source type configurations.
2) Each of the turbulence schemes: kε, Rkε, kω, and SST was used in combination with the discrete phase model for predicting respective particulate matter concentrations. The performance of turbulence schemes was determined from the statistical analysis by comparing their respective modeled concentrations with the observed ones as a function of different wind speed categories: low (u < 2 m/s), medium (2 < u < 5 m/s), and high (u > 5 m/s) irrespective of the particle sizes. For u < 2 m/s, the statistical analyses concluded that Rkε was the best performing among the turbulence models. Also, the results of 95% confidence limits (CLs) on statistical measures showed Rkε performing significantly different from kε, kω, and SST for u < 2 m/s. At low wind speeds, the eddy energy dissipation rates of kε, kω, and SST were much higher than that of Rkε and hence they over-estimated particulate concentrations. For 2 < u < 5 m/s, the results of both statistical analyses and confidence limits concluded that Rkε was the best performing turbulence model. For high wind speeds (u > 5 m/s), the statistical results and confidence limits showed similar performances of Rkε, kω, and SST indicating the need for more particulate matter data sampling studies for further statistical evaluation to determine the best performing turbulence model. Also, there was no significant difference in the eddy energy dissipation rates of the turbulence models at u > 5 m/s. When other important model inputs such as emission rates were taken into account, their effect on modeled particulate concentrations was less than low wind speed condition (u < 2 m/s).
It is recommended that further field studies should emphasize the sampling of particulate matter with respect to particle size for various downwind distances as well as different meteorological conditions. Modeling work should consider Brownian motions because they are significant for particles' sizes less than 1 micron. Also, attempts should be made to consider the interaction of particles and the medium during modeling.