A Monte Carlo Model for Soil Particle Resuspension Including Saltation and Turbulent Fluctuations

Abstract

This article describes a stochastic model for resuspension that combines both aerodynamic entrainment and momentum transfer from saltating particles. Two case studies are performed based on soil conditions for the topmost layer of soil in Los Angeles County, CA, and Allegheny County, PA. Wind friction velocity, u∗, and soil size distribution were found to be the most important factors in predicting mass and number flux. Under a broad range of wind conditions mass and number fluxes agree to within an order of magnitude with the empirical models of Marticorena and Bergametti (1995) and Ginoux et al. (2001) at u∗ ≤ 0.4 m/s. For u∗ ≤ 0.60 m s⁻¹ and u∗ ≥ 0.85 m s⁻¹ aerodynamic forces and splash were the dominant resuspension mechanisms, respectively. Flux was sensitive to wind speed but was not proportional to u∗³. The mass and number distributions with height peaked at heights corresponding to the maximum concentration of saltating particles and the maximum concentration of suspended particles, respectively. Particles that are most likely to resuspend in the absence of saltation are < 10 μm or > 100 μm in diameter. The average particle diameter increases with height but is consistently less than the average particle diameter of the parent soil. Simulations reached steady state in approximately 0.01 seconds, and an alternative method of predicting the reduction in near-surface wind speed as a result of saltation is presented as a component of the model.

INTRODUCTION

The first significant look at soil resuspension was taken by Bagnold (1941) and later Chepil (1945a, b). Bagnold (1941) identified three types of particle movement: surface creep, suspension, and saltation. Surface creep is defined by the rolling of the largest particles along the surface. Saltation is the movement of particles characterized by smooth, approximately parabolic trajectories. Particles in suspension reach significant heights above the surface and are free from surface effects. In this study we focus on suspension and saltation.

The mechanisms responsible for initiating particle motion can be either aerodynamic effects from wind or momentum transfer from particles in saltation. It is known that significantly more particles resuspend in the presence of saltating particles than in the absence (e.g., Shao et al. 1993), but the relative importance of the two resuspension mechanisms has yet to be established. Several authors have hypothesized that particle motion downstream from an actively saltating bed is due exclusively to saltation (Shao et al. 1993; Anderson and Haff 1988; Chepil 1945b; Bagnold 1941). In the case that no large particles are available for saltation, direct wind pick up has been observed (Kjelgaard et al. 2004; Braaten et al. 1990; Garland 1982). Most existing resuspension models assume a single launch mechanism—either saltation or direct pick up by wind. Additionally, in the saltation layer the wind field is a source of considerable uncertainty; the role of coherent structures in the viscous sublayer during saltation in launching particles has yet to be explored.

Most models characterizing the mass flux of resuspending grains are empirical and hold only for the range of field and meteorological conditions under which experiments were run. Frequently models are based on monodisperse particle sizes; there is inadequate knowledge about the movement of polydisperse beds. In addition to the mass flux, little is known about the number flux of crustal grains, which is a significant source of uncertainty in global climate models.

The first goal of this work is to build a physics-based, mathematical model of saltation within 0.4 m of the surface and within the first three seconds after a gust of wind has blown to gain an accurate representation of the initiation of resuspension on a minute scale. The potential to scale up to heights in excess of 0.4 m is possible making this model a resource for climate models and global and regional transport models. To scale up, either results at a height of 0.4 m can be used as an initial condition for other models or greater accuracy in the trajectories of suspended particles can be incorporated into the current model theoretically allowing for flux predictions to 10 m or more (Harris 2007). The simplified suspended particle trajectories that limit model validity to heights of 0.4 m or less were necessary to limit computing time.

From this model we seek to determine (1) the relative importance of saltation and aerodynamic forces as resuspension initiation mechanisms, (2) the dependence of mass and number flux on wind speed as compared with experimental results and empirical models, (3) the mass and number concentrations of crustal material as functions of height above the surface including the average particle diameter, (4) the probability of single particle resuspension in the absence of saltation, and (5) the impact of saltation on wind shear.

METHODS

We present a sequential Monte Carlo simulation for a bed of polydisperse soil particles distributed randomly. The simulation runs as follows. At time t _o all particles are at rest. At time t _o + Δ t, wind entrains the particles most susceptible to resuspension following our earlier model of aerodynamic entrainment of a single particle (Harris and Davidson 2008). Based on kinetic energy analyses each particle is classified as suspended or saltating, and energy thresholds determine if saltating particles carry enough momentum to resuspend other bed particles. An order of magnitude analysis leads to the conclusion that particles in surface creep lack sufficient momentum to launch other particles. We assume neutral atmospheric conditions, spherical particles, and dry, unvegetated soil surfaces throughout. A schematic of the model is shown in Figure 1.

FIG. 1 Schematic of Monte Carlo simulation.

Modeling Domain

The spanwise width of our modeling domain is that of a single streamwise streak or 100 ν/u∗ where ν is the kinematic viscosity and u∗ is the friction velocity. This has been identified as the minimum width that will sustain turbulence (Jiménez and Moin 1991). We place a limit on the height of 0.4 m because simplified particle trajectories used in this model are not consistent with expected trajectories at heights above 0.4 m. We model saltation over the relatively short length of 5 m using the method of Anderson and Haff (1988) and McEwan and Willetts (1991) who use a modeling domain with periodic boundaries in which any particle that bounces out of the box is assumed to bounce in on the other side. In this manner, we can drastically reduce computation time by avoiding the issue of an increasing number of particles entrained in a flow with downwind distance also called “fetch.”

Characterizing the Soil

As case studies we look at silt loam as a characteristic soil from Los Angeles County, California and clay loam as a characteristic soil of Allegheny County, Pennsylvania (USDA-NRCS, 2003). These two distributions have mean diameters of 102 μ m and 53.5 μ m, respectively. Performing chi-squared tests on exponential, gamma, lognormal, normal, and Weibull distributions shows that the last best represents both distributions. Parameters for the distributions as determined by the method of moments are α = 0.396 and β = 42.5 for Los Angeles County and α = 0.622 and β = 23.9 for Allegheny County (Harris 2007). For modeling purposes we break the size distribution into 100 discrete size bins; each contains 1% of the mass.

Characterizing the Wind

The mean wind profile in the absence of saltation is given by:

where u is the mean horizontal wind profile for a rough surface with a displacement height of zero, k is von karman's constant, z is height, and z _o is the roughness height. The roughness height is taken as 0.27 D _p50 below the top of the roughness elements, ε (Blinco and Partheniades 1971), where D _p is particle diameter, the number in the subscript denotes the percentile of the mass distribution, and ε is taken as 2 D _p50 (Wu and Chou 2003; Cheng and Chiew 1998; van Rijn 1985a). The resulting value of z _o is of the same order as the height of the viscous sublayer sometimes defined as δ = 5ν/u∗ (Davidson and Wu 1990).

We generated a two-dimensional wind field using a bivariate normal distribution modified by a Gram-Charlier expansion. This probability distribution function (pdf) has been shown to reproduce accurately the coherent structures in the near-wall region (Wu and Jiang 2007; Nakagawa and Nezu 1977). According to this pdf, turbulent bursts occur on average 33% of the time at any given location. For the spatial distribution of coherent structures we look to Jiménez et al. (2005) who observe that bursts are distributed at distances approximately 1200 ν/u∗ downstream from one another. To accommodate both the spatial and temporal distribution of bursts, we break the modeling domain into grid cells that are 400 ν/u∗ in the streamwise direction by 100 ν/u∗ in the spanwise direction. The actual area over which a turbulent burst acts is approximated as π (10ν/u∗)² following the work of Cleaver and Yates (1973) and Yung et al. (1989). A burst is defined here as the situation in which turbulent wind fluctuations in the viscous sublayer exhibit a characteristic pattern of u′ < 0 and w′ > 0 where u′ and w′ are the turbulent components of the instantaneous wind speed in the horizontal and vertical directions, respectively. Bursts are also called “ejections” in some of the literature.

Particle Trajectories

From the force balance model described in Harris and Davidson (2008), particle velocities after liftoff are characterized by a Gaussian distribution in both the horizontal (x) and vertical (z) directions. Here we truncate the distribution so that particles cannot be launched from the surface in the negative z-direction. The horizontal component (standard deviation = 0.68 m s⁻¹, R ² = 0.96) and vertical component (standard deviation = 0.04 m s⁻¹, R ² = 0.90 and 0.99, respectively) of the initial particle velocity are (Harris and Davidson 2008):

where u _p and w _p are the mean particle velocities in the x- and z-directions, respectively in m s⁻¹, and D _p is in microns. Thus the model predicts that particles of 100 μ m will achieve the greatest vertical speeds. Aluko and Noll (2006) observed this peak liftoff diameter at approximately 80 to 100 μ m in experiments and approximately 30 to 50 μ m in theoretical calculations.

Full suspension and saltation are distinguished by a force balance; particles large enough for the gravitational force to exceed the drag force remain in saltation while particles experiencing a stronger drag than gravitational force are considered in full suspension (Greeley and Iversen 1985). Following Nishimura and Hunt (2000) and others (Lee and Hsu 1994; Nalpanis et al. 1993) we use a ballistic trajectory for saltating particles governed by their initial velocities. This slightly overestimates the height of a particle, but is consistently within 10% of observed trajectories (Anderson and Hallet 1986). In contrast, particles in suspension are assumed to follow a straight line trajectory at constant velocities equal to their initial velocities. Initial suspension velocities are determined from the force balance and the instantaneous two-dimensional wind field described in Harris and Davidson (2008). The straight line trajectory is a reasonable approximation at heights relatively close to the surface, but such a simplification limits the validity of the model to heights below 0.4 m.

It is assumed that no particle ≤0.1 μ m in diameter is capable of being resuspended since these particles are primarily clay and are bound very tightly to other soil grains (e.g., Reeks et al. 1988). Furthermore particles of this size with crustal origin are rarely observed in the atmosphere (Ginoux et al. 2001; Tegen and Lacis 1996). Soil particles are assumed to be well mixed by size.

Rebound and Splash

A particle will remain in saltation only if its kinetic energy is greater than zero after collision with the surface. This implies that its initial kinetic energy upon impact exceeds the sum of energy lost due to van der Waals force and deformation of roughness elements on the surface of the particle. Therefore a particle will rebound only if its velocity exceeds a critical impact velocity V _cr :

where m is particle mass, and E _cs and E _asp are the kinetic energies lost by the particle upon collision due to van der Waals force and deformation of surface asperities, respectively. Details for calculating E _cs and E _asp are given in Tsai et al. (1990). Our analysis shows that V _cr is nearly independent of both the particle's kinetic energy and material properties and depends almost exclusively on particle radius. Estimates of surface asperities were taken from the experiments of Leite et al. (2003) and Matsui et al. (1980) on quartz sand. Power law fits to Equation (3) over a wide range of observed kinetic energies yield:

for silicon dioxide and lead particles, respectively, where D _p is in microns and V _cr is in m s⁻¹. Both equations have R ² = 0.99.

A saltating particle striking the surface may launch other particles into motion, a process known as a splash or sandblast effect. In the modeling efforts of Grini and Zender (2004), excluding splash led to a significant underprediction of the quantity of resuspended aerosols. Their results were improved when splash was included. A single impact from a saltating particle typically resuspends one to six surface particles (Shao and Li 1999; Rice et al. 1995; Anderson and Haff 1988). The velocity of splashed particles is usually between 3 and 10% and is most often found between 7 and 9% of the impact velocity (Rice et al. 1995; Anderson and Haff 1988; Anderson 1987). Here we account for splash through energy conservation equations.

Saltation/Wind Feedback

The logarithmic profile for the wind does not hold in the saltation layer (Nishimura and Hunt 2000; Sekine and Kikkawa 1992; Owen 1964), although there is still significant disagreement about the shape of the mean wind profile. Particles in saltation exert a drag force on the wind, thereby reducing the wind speed in the saltation layer, which subsequently reduces the number of particles that become airborne. The net effect is a feedback cycle between fluid shear stress and the mass of saltating particles (Doorschot and Lehning 2002; Hopwood and Scott 1997; Anderson et al. 1991; Ungar and Haff 1987; Owen 1964). Thus the total number of saltating particles tends to fluctuate around an equilibrium value.

The mean wind profile within the saltation layer as a function of τ _p is (Raupach 1991):

where τ is the total shear stress, τ _p is the shear stress carried by particles, and ρ is the density of air. Equation (6) is integrated numerically to obtain u as a function of z. This equation correctly predicts a reduction of the wind speed in the saltation layer that is most pronounced at the top of the saltation layer and is virtually zero very close to the surface. The standard logarithmic profile, Equation (1), continues to hold above the saltation layer (Nishimura and Hunt 2000; Sekine and Kikkawa 1992; Owen 1964).

In the model, the boundary condition for Equation (6) is u(z = z _o ) = 0. The value of z _o is assumed fixed even though the roughness height increases during saltation since saltating particles emulate an extension of the surface. While the effective roughness height is an important consideration for characterizing the wind profile above the height of the saltation layer (Bintanja 2000; Déry and Taylor 1996; Sekine and Kikkawa 1992; Raupach 1991), to our knowledge it has not been applied to velocity profiles for z < h_s. where h_s is the average maximum height of saltating particles. Because we calculate trajectories based on initial launch velocities the implications for changes in the roughness height in this model are minimal. Similarly, friction velocities are reduced in a wind profile modified by saltation. The density of the fluid increases from an increased concentration of suspended particles while the total shear stress remains constant. For the classic relationship of τ = u∗²ρ to hold, we must observe a corresponding reduction in u∗ (Bintanja 2000; Déry and Taylor 1996). Throughout this article wind fields are characterized by their initial friction velocities.

In the absence of saltation all of the shear stress near a surface is the result of wind shear, τ _a . As saltation begins, momentum is transferred to airborne particles, and some of the shear stress is then associated with these particles, which we term τ _p . Thus, τ _p has some positive value within the saltation layer, defined by z _o < z < h _s and τ _p is 0 for z > h _s . The hypothesis first put forward by Owen (1964), that the total shear stress τ = τ _a + τ _p remains constant, has been used to good effect in several saltation models (Doorschodt and Lehning 2002; Shao and Li 1999; Hopwood and Scott 1997; Sekine and Kikkawa 1992). Particle shear stress is quantified by Raupach (1991) as:

where i denotes the particle size bin, n is the number flux of saltating particles (m⁻²s⁻¹), m _i is the median mass of all particles in the size bin, and u _↓i (z) and u _↑i (z) are the horizontal velocities of a particle on the descent and ascent, respectively, as functions of height z and the drag force exerted by the wind on the particle. Particles tend to land with a larger velocity than their initial velocities as they acquire momentum from the wind during flight. Equation (7) was chosen over empirically based equations (Doorschodt and Lehning 2002; Shao and Li 1999; Hopwood and Scott 1997) that require knowledge of the poorly defined threshold friction velocity for saltation.

To parameterize the length of time a typical saltating particle stays airborne, the ratio of particles resuspended by aerodynamic forces versus splash is taken as 3:11 as suggested by the data of Rice et al. (1995). This is used to predict the time dependence of the number of particles in saltation, i.e., n = n(t), in Equation (7).

There is considerable uncertainty in turbulent wind fluctuations within the saltation layer. We therefore adopt the bivariate Gaussian distribution modified by the Gram-Charlier expansion in Harris and Davidson (2008). A correlation coefficient of −0.45 for turbulent fluctuations in the horizontal, u′, and vertical, w′, directions which is used for the viscous sublayer (Nezu 2005; Pope 2000; Tennekes and Lumley 1972) is also adopted in the saltation layer. Raupach (1991) explored the idea of using forest canopies as proxies for the saltation layer, and preliminary evidence suggests that their mean wind profiles are not dissimilar (Raupach 1988, 1991). However the standard deviations of turbulent wind fluctuations within forest canopies are approximately 4 u and 2 u in the horizontal and vertical directions, respectively (Raupach 1988). Applying such large standard deviations to the saltation layer produces unreasonably large wind velocities. Within the saltation layer Nishimura and Hunt (2000) measured a standard deviation of wind speed (direction not specified) in a wind tunnel of 2u∗, which is the same order of magnitude as the standard deviations found in Wu and Jiang (2007) and adopted by Harris and Davidson (2008). Experimental measurements of turbulence in the saltation layer are needed.

Description of Model

Our model is executed in Mathematica 6.0 and runs as follows. Time is advanced with steps of 100 ν/u∗², and simulation times are on the order of 3.0 seconds to allow adequate time to reach equilibrium. The modified mean wind profile is calculated at the beginning of each time step by solving Equation (6) via Mathematica's numerical solver. The activity in each grid cell is calculated sequentially. For each a two-dimensional instantaneous velocity is generated randomly using an acceptance-rejection routine and the Gram-Charlier distribution described above. Second, to simulate additional wind pick up of surface particles, one particle size is chosen randomly from the Weibull distribution followed by another particle size chosen randomly, and continuing until the sum of particle diameters exceeds the diameter of a turbulent burst (20 ν/u∗). The instantaneous wind speed is subsequently assessed to determine if it is sufficient to overcome the adhesion force, calculated in the manner of Harris and Davidson (2008). If the wind speed in either the horizontal or vertical direction is not adequate to instigate liftoff, another particle size is generated. However, if both the x and z components of the wind are sufficient to launch a particle aerodynamically, then a force balance is used to classify the particle as either suspended or saltating and to calculate its launch velocity.

For a particle in suspension, its height at the end of the simulation is calculated based on its initial velocity and the current time step. The particle is added to the overall number flux matrix, binned by particle size. For a particle in saltation, it is determined whether its initial velocity exceeds V _cr calculated via Equation (4) and if the particle will strike the surface before the end of the simulation. If both conditions are not met, a new particle size is generated. However if these conditions are met, then the splash function is initiated.

To address splash it is first assumed that the original saltator strikes five bed particles, randomly distributed in size. For each of these five particles, velocities are assigned as 0.5, 0.6, 0.7, 0.8, and 0.9 times the velocity of the initial saltator. Each particle is classified as saltating, suspended, or incapable of leaving the surface. Those splashed particles that are either in saltation with velocities that exceed V _cr or in suspension are added to the number flux and their heights at the end of the simulation are returned. We calculate a probability of < 1% that particles launched by splash will have sufficient energy to rebound and therefore neglect any splash caused by particles launched by saltators.

The new velocity of the original saltator is calculated by subtracting kinetic energy lost during rebound and splash. If its velocity still exceeds V _cr and will strike again before the end of the simulation, the splash function is repeated. After several rebounds the original particle will either have a velocity that does not exceed V _cr or will not rebound before the end of the simulation. The splash function is then exited and the system returns to the aerodynamic analysis. All saltating particles are assumed to land within the model domain in accordance with the periodic model of McEwan and Willetts (1991) and Anderson and Haff (1988).

When the sum of the diameters of particles launched aerodynamically within a single grid cell meets or exceeds the diameter of a turbulent burst, the model advances to the next grid cell. When the program terminates, the number of particles currently in saltation is summed. For u∗ = 1.0 m/s, this is approximately the number of saltating particles in the last 4 time steps. When u∗ = 0.75 m/s, 0.50 m/s, or 0.25 m/s, it is approximately the number of saltating particles in 2, 1, or 0 time steps, respectively. For small u∗, τ _p is set to 0 since our model indicates that saltation occurs very rarely at u∗ ≤ 0.25 m/s.

Results are returned as comma delimited text files and include: the total number and mass of particles entering suspension, the total number and mass of particles entering saltation, the sum of the heights achieved by particles at the end of the simulation, the number of particles still in saltation at the end of the simulation, and the average wind speed at z = 0.0015 m at each time step to capture the effect of saltation on the wind speed over time. With the exception of the wind speed, all results are sorted into 100 particle size bins.

RESULTS

Emission Rates

Results of the model indicate that at u∗ = 0.60 m s⁻¹ and below, virtually all of the resuspension is a result of aerodynamic forces. In fact the mass flux is slightly higher when splash is neglected at low friction velocities, presumably because the reduction in wind shear due to saltation is neglected when splash is neglected. For u∗ > 0.60 m s⁻¹, momentum transfer from saltating particles is responsible for some resuspension. At friction velocities ≥0.85 m s⁻¹, splash resuspends nearly all of the particles once equilibrium is reached. Emission rates by mass (kg m⁻²s⁻¹) from this study are compared with the literature in Table 1.

TABLE 1 Emission rates (kg m^–2s⁻¹) within the saltation layer except where noted. These data correspond to soils with different and, in some cases unknown, size distributions

	u∗(m/s)	1.0	0.85	0.75	0.60	0.50	0.25	0.10
	Los Angeles County	6.83 × 10^−3	1.22 × 10^−6	3.09 × 10^−9	8.42 × 10^−10	3.21 × 10^−10	6.40 × 10^−12	1.87 × 10^−14
	Allegheny County	2.53 × 10^−5	3.61 × 10^−6	1.16 × 10^−8	2.07 × 10^−9	1.08 × 10^−9	1.80 × 10 ^−11	2.14 × 10^−14
Huan et al. 2006	Wind tunnel		0.8	0.4	0.12
Shao and Mikami 2005	Takla Makan desert					∼0.001		0 to 0.01
Shao et al. 2003	Gobi desert (at rooftop height)	∼5 × 10^−6	∼3 × 10^−6	∼1 × 10^−6	∼1 × 10^−7	∼5 × 10^−8
Chen and Fryrear 2002	Big Spring, TX (at 0.25 m)	1.5 × 10^−2
Westphal et al. 1987	Model (PM 10 only)	∼1 × 10^−5	10^−7 to 8 × 10^−6	8 × 10^−7 to 10^−5	∼1 × 10^−6	3 × 10^−8 to 10^−6	2 × 10^−8 to 10^−7
Anderson and Hallet 1986	Model					0.05
Anderson and Hallet 1986	Bagnold model					0.07
Anderson and Hallet 1986	White model					0.11 to 0.13
Chepil and Woodruff 1957	Syracuse, KS (at 6 ft.)	1.03 × 10^−3	1.50 × 10^−4	5.51 × 10^−5		5.16 × 10^−6

These results agree qualitatively with the channel flow experiments of Niño and Garcia (1998) and the model of Doorschodt and Lehning (2002). In Niño and Garcia's experiments an estimated 60 to 92% of particles in water were resuspended as a result of hydrodynamic effects, while 8 to 40% were resuspended as a result of particle collisions. This was observed at relatively low friction velocities between 0.021 and 0.032 m s⁻¹. The model of Doorschodt and Lehning (2002) for air flow indicates that at friction velocities below about 0.4 m s⁻¹, aerodynamic forces are the dominant entrainment mechanism. Above ∼ 0.4 m s⁻¹ momentum transfer from saltating particles becomes important. These results suggest a middle ground between the “Owen hypothesis” (Owen 1964) in which saltating particles reduce the wind shear to just above the threshold velocity, and particles therefore continue to be entrained as a result of aerodynamic forces even during steady-state saltation and the other extreme suggested by the steady-state model of 10⁵ to 10⁹m⁻²s⁻¹ depending on the size and number of saltating Anderson and Haff (1988) in which aerodynamic forces played little role and particle-to-particle interactions dominated for 1000 μ m, monodisperse particles.

Many authors have estimated that the mass flux of resuspended particles is proportional to u ^a , where u is the wind speed and a is a constant on the order of 3 (Bagnold 1941; Garland 1979; Sehmel 1980; Nicholson 1988). Wind tunnel experiments performed by Fairchild and Tillery (1982) suggest that the resuspension rate is proportional to u∗^1.3 in the absence of saltating particles, but that the exponent increases to 1.8 and 2.1 in the presence of 100 μ m and 200 μ m saltating particles, respectively. The work of Sorenson (1985) shows a mass flux proportional to u∗^2.8. The saltation simulations of Ungar and Haff (1987) show a u∗² dependence for the mass flux. Measurements taken in the Takla Makan Desert scale somewhat to u∗³ but with considerable scatter (Shao and Mikami 2005). In contrast, Werner (1990) argues that mass flux does not have a simple dependence on u∗. For our model linear regressions of the friction velocity cubed (m³ s⁻³) against the predicted mass flux (kg m⁻² s⁻¹) are:

for Los Angeles and Allegheny Counties, respectively. This is shown graphically in Figure 2. R ² values improve if the exponent on u∗ is increased above 3.

FIG. 2 Mass flux rates within the saltation layer plotted against u∗³. Conditions were modeled to match Los Angeles and Allegheny County soil shown by diamonds and squares, respectively. For comparison purposes, flux estimates are shown from the models of Marticorena and Bergametti (1995) and Ginoux et al. (2001). The fraction of soil available for resuspension was assumed to be 0.2 in both cases. The Ginoux et al. lines are virtually identical for the two case studies illustrated.

The estimated emission rates for particles by number (m⁻²s⁻¹) are listed in Table 2. This compares well with the wind tunnel measurements of Fairchild and Tillery (1982) who, at u∗ = 0.43 m/s, found particle number fluxes in the range of 10⁵ to 10⁹m⁻²s⁻¹ depending on the size and number of saltating particles. Shao and Mikami (2005) observed peak number fluxes ranging from ∼ 300 m⁻²s⁻¹ for particles in the 380 to 401 μ m size bin to 1.3 × 10⁶ m⁻²s⁻¹ for particles in the 57 to 72 μ m size bin during field experiments in Qira, China.

TABLE 2 Simulated emission rates for particles within the saltation layer over a range of friction velocities. Units are m⁻²s⁻¹

u∗(ms^−1)	LA County	Allegh. County
0.10	5.22 × 10^5	3.40 × 10^5
0.25	3.61 × 10^7	4.71 × 10^7
0.50	6.20 × 10^8	1.08 × 10^9
0.60	1.29 × 10^9	1.89 × 10^9
0.75	3.19 × 10^9	5.98 × 10^9
0.85	4.90 × 10^9	7.36 × 10^9
1.0	9.21 × 10^9	1.71 × 10^10

Mass Concentrations

Results of the model show two peaks in the distribution of the mass with height. The first peak occurs at ∼ 3 mm or below, where mass concentration due to saltation is at a maximum. A secondary peak occurs at a higher height corresponding to the mean height that suspended particles reach, which is dependent on simulation time. Some example results for the mass concentration as a function of height are shown in Figure 3 for three second simulation times. Peak concentrations in the saltation layer are about 1.05 × 10⁶μ g m⁻³ for Allegheny soil and 335 × 10⁶μ g m⁻³ for Los Angeles soil—both at u∗ = 1.0 m s⁻¹ and z = 3 mm. Similarly, peak dust storm concentrations measured at a height of 1 m near Big Spring, Texas were 1.3 × 10⁶μ g m⁻³ for u∗ ∼ 1 m s⁻¹ (Chen and Fryrear 2002).

FIG. 3 Mass concentrations of suspended and saltating particles with friction velocities of 0.85 m s⁻¹ and 1.0 m s⁻¹. The simulation time was 3 seconds in all cases.

Graphing the concentration of suspended particles with height (saltating particles not shown) indicates mass concentrations that are essentially Gaussian according to chi-squared tests and this is consistent with the Gaussian distribution of initial particle velocities described above. Plots shown in Figures 4 and 5 show rough agreement with the saltation model of Anderson and Hallet (1986), which indicates a single peak concentration of suspended material around 0.02 m for u∗ = 0.53 m s⁻¹. Near surface concentrations shown in Figures 5 and 6 for u∗ in the range of 0.3 ≤ u∗ ≤ 0.6 m s⁻¹ are in the range of 0 and ∼ 130 μ g m⁻³, which agree to an order of magnitude with daily near-surface concentrations, for particles not distinguished between suspension and saltation in Puerto Rico (Grini and Zender 2004).

FIG. 4 The mass concentration of particles in full suspension for friction velocities of 0.50, 0.75, and 1.0 m s⁻¹. Soil is characteristic of Los Angeles County. Particles in saltation are not shown.

FIG. 5 The mass concentration of particles in full suspension for friction velocities of 0.50, 0.75, and 1.0 m s⁻¹. Soil is characteristic of Allegheny County. Particles in saltation are not shown.

FIG. 6 The number concentration of suspended and saltating particles entrained in flows with u∗ = 0.85 m s⁻¹ and 1.0 m s⁻¹. The simulation time was 3 seconds in all cases.

Number Concentrations

Number concentrations (m⁻³) are shown as a function of height in Figures 6–8. Figure 6 shows the number flux of both saltating and suspended particles at friction velocities of 0.85 and 1.0 m s⁻¹ over a simulation time of 3 s. Although no direct comparisons with experiment can be made at this time, concentrations in Figure 6 are comparable to measurements of saltating and suspended snow particles at a height of 1 m with a 10 m wind speed of 15 m s^{− 1}, which were found to be on the order of 10⁷ particles m⁻³ per micron (Déry and Taylor 1996). Figures 7 and 8 show the number concentration for suspended particles also for 3 s. A unimodal peak is observed in Figures 7 and 8 corresponding to the mass concentration peaks shown in Figures 4 and 5.

FIG. 7 The number concentration of particles in full suspension for u∗ = 0.50, 0.75, and 1.0 m s⁻¹. Soil is characteristic of Los Angeles County. Particles in saltation are not shown.

FIG. 8 The number concentration of particles in full suspension for u∗ = 0.50, 0.75, and 1.0 m s⁻¹. Soil is characteristic of Allegheny County. Particles in saltation are not shown.

Particle Size Distribution

Figure 9 shows particle size distribution results of simulations run at u∗ = 0.75 m s⁻¹ for Los Angeles and Allegheny County soils for suspended soil grains only. Mean diameters at all heights are significantly smaller than mean and median particle sizes of the parent soils. The mean particle diameter shows a marked increase with height up to ∼ 0.4 m. Above 0.4 m the mean diameter begins to decline with height. The initial increase is a direct result of the liftoff velocities scaling with particle diameter up to 100 μ m as described in the Methods section. At first glance this appears to contradict the logarithmic or power law decrease in mean diameter with height that is consistently observed for eroded material (e.g., Chen and Fryrear 2002; Tegen and Lacis 1996; Ahmed et al. 1987; Nickling 1983). However very few studies have measured the size distribution at heights less than 0.5 m from the surface. In one notable exception measurements taken at 2 cm height increments in Big Spring, Texas show an increase in mean diameter with height up to 0.35 m (Chen and Fryrear 1996) and indicate a trend similar to our simulation results. Note that the mean diameters below 0.35 m in the Texas study are significantly larger (∼ 200 μ m) than any of the mean diameters found in our case studies (6.9 and 8.4 μ m for Los Angeles County and Allegheny County, respectively, at 0.3 m and u∗ = 0.75 m s⁻¹). This is to be expected as only 14% of particles by mass had diameters less than 63 μ m in Big Spring soil. In Los Angeles County approximately 69% of the particle mass is smaller than 63 μ m and in Allegheny County about 84% of particles by mass fall below 63 μ m in the surface layer of soil. Measurements and simulations in Qira, China for particle sizes in the range of 40 μ m to 660 μ m did not show a clear trend below 0.3 m, but do seem to indicate increasing concentrations with height in the 40 to 100 μ m size range and decreasing concentrations with height in the > 100 μ m size range (Shao and Mikami 2005) consistent with our results.

FIG. 9 The mean diameter calculated by both particle mass and particle number for Los Angeles and Allegheny Counties. Friction velocity is 0.75 m s⁻¹.

Third order polynomial regressions of the mean diameter with height were found to agree reasonably well with model results. This equation is given as:

where D _p is in μ m, z is in m, and the coefficients a through d are given in Table 3.

TABLE 3 Regression coefficients and R ² values for Equation (9)

	a	b	c	d	R ^2
LA by mass	−568	340	−37.4	2.923	0.966
LA by number	−340	280	−43.7	1.698	0.890
Allegh. by mass	−554	367	−44.5	3.407	0.980
Allegh. by number	−87.3	179	−33.9	1.515	0.917

Probability of Resuspension

Figure 10 shows the results of simulations with aerodynamic entrainment forces only. The probability that a particle will be entrained into the flow increases with friction velocity as expected. The probability of resuspension for small particles (< 10 μ m) appears to be insensitive to wind speed in the range of u∗ = 0.25 m s⁻¹ to 1.0 m s⁻¹. However, for large particles (> 100 μ m) the probability varies from near zero to 0.19 over the same velocity range. Figure 10 exhibits the same basic shape as entrainment probabilities predicted for Reynolds numbers in the range of 0.1 to 1000 (Wu and Chou 2003) where the Reynolds number is proportional to diameter (Re = u∗ D _p /ν).

FIG. 10 Aerodynamic entrainment probabilities for a range of particle diameters and u∗ = 0.50, 0.75, and 1.0 m s⁻¹.

A threshold instantaneous wind speed is derived for a force balance between aerodynamic forces and adhesion forces. In its exact form the threshold wind speed u _t is a transcendental function with complicated dependencies on D _p and u∗. To produce the values shown in Figure 10, the instantaneous threshold wind speed was simplified somewhat:

where λ is the mean free path of atmospheric particles used in the slip correction, z _s is the separation distance between the bottommost point of the particle and the surface, and h is the Lifshitz-van der Waals constant. Equation (10) overestimates u _t (∼ 10%) at small diameters (< 20 μ m) and underestimates u _t (∼ 5%) at large diameters (∼ 1000 μ m).

Modified Wind Conditions

As described in the Methods section, the presence of saltating particles reduces the effective wind speed. One example of the modified horizontal mean wind profile with height is shown in Figure 11 for u∗ = 0.75 m s⁻¹. An ambient roughness height of z _o = 10⁻³ m was used in this simulation. At the top of the saltation layer (0.01 m) the mean velocity in the absence of saltation is about 4.21 m s⁻¹. In steady-state saltation this is reduced to about 1.51 m s⁻¹. For purposes of calculating the modified wind profile, it was assumed that only saltating particles reduced the effective wind velocity; particles in suspension were neglected.

FIG. 11 The mean horizontal wind profile with height during steady–state saltation is shown by the dotted line. A standard logarithmic profile is shown as a solid line for comparison. This is for u∗ = 0.75 m s⁻¹, z _o = 10⁻³ m, and a number flux of 3 × 10⁻⁹ m⁻² s⁻¹.

According to our model, the saltation flux reaches steady state ∼ 0.01 seconds into the simulation regardless of initial wind speed or soil type for the two soil types tested. The prediction is significantly smaller than previous estimates of 1–2 seconds (Anderson and Haff 1988; 1991) and 2.0 seconds (McEwan and Willetts 1991) for the time to steady state, which corresponds to several leaps of saltating particles. Our value of 0.01 seconds corresponds to about 5 time steps within the model where the time to steady state was presumed to be the time it takes for the horizontal wind speed to decline from its initial value to the lesser value around which the wind speed fluctuates for the remainder of the simulation. That we obtained such a small value for the time to steady state may be an artifact of our model because in each time iteration the splash function is carried through until all of the particles splashed in a given time step have come to rest. Therefore the results of a single time iteration are not necessarily representative of the number of particles in saltation at a given time.

CONCLUSIONS

A model of particle resuspension has been built that accounts for both aerodynamic entrainment and momentum transfer from saltating particles in a polydisperse bed. The predicted fluxes from this model have been analyzed for two case studies and are generally in rough agreement with experimental results.

The friction velocity and the particle size distribution of the parent soil are the two most important factors in predicting soil flux. Soils with many large particles show a much greater flux due to the presence of saltating particles. For friction velocities above 0.6 m s^{− 1}, accounting for particle splash increases the quantity of resuspended mass by several orders of magnitude. For friction velocities below this value, including the splash function leaves the flux virtually unchanged. This result meets the first goal of determining the relative importance of saltation and aerodynamic forces as entrainment mechanisms. Aerodynamic entrainment is almost exclusively responsible for entraining particles at low friction velocities but that entrainment due to momentum transfer from saltating particles is vitally important at greater friction velocities. The reduction in near-surface wind shear stress as a result of saltation is reproduced by the model.

In meeting the second goal of determining the dependence of mass and number flux on wind speed we observe the following. Mass emission rates scale with wind speed and range from virtually zero at u∗ = 0.1 m s⁻¹ to 6.83 × 10⁻³ kg m⁻² s⁻¹ at u∗ = 1.0 m s⁻¹. Emission rates by number range from 3 × 10⁵ to 9 × 10⁹ m⁻² s⁻¹ over the same range of wind speeds. Both mass and number flux rates agree qualitatively with experimental results of other investigators. However, the mass flux shows only weak evidence to suggest a proportionality to u∗³.

The conclusions with respect to the third goal of finding mass and number concentration with height and average particle diameter with height are shown above. The concentrations of both mass and number with height show peaks at two different heights. Mass and number concentrations are very high within the saltation layer, and reach a second peak at the mean height that suspended particles reach at the end of a simulation. The concentration of suspended particles alone has a single peak. The height of this peak varies with simulation time, but is approximately 0.09 m after 3 s. The mean particle diameter increases significantly with height to 0.4 m. Suspended particle size distributions are skewed towards small particle sizes compared to the parent soil.

Meeting goal four of finding the probability of particle resuspension in the absence of saltation, we found that particles most likely to resuspend particles are < 10 μ m and > 100 μ m in diameter. For particles < 10 μ m, the probability of resuspension is relatively insensitive to wind speed. The wind speed makes a dramatic difference for particles larger than 100 μ m in diameter.

Finally, the vertical wind profile was shown to change in the presence of saltating particles within the saltation layer. The mean wind speed was significantly reduced from a standard logarithmic profile and all heights less than h _s .

The authors acknowledge the helpful insights of Mitchell Small, Nadine Aubry, Peter Adams, Spyros Pandis, and Yunha Lee. The authors are also grateful for the invaluable discussions with others in the Center for Atmospheric Particle Studies. This work was funded by National Science Foundation Grant BES-9714162 and a National Science Foundation Graduate Student Fellowship.