Development and evaluation of a physics-based windblown dust emission scheme implemented in the CMAQ modeling system.

A new windblown dust emission treatment was incorporated in the Community Multiscale Air Quality (CMAQ) modeling system. This new model treatment has been built upon previously developed physics-based parameterization schemes from the literature. A distinct and novel feature of this scheme, however, is the incorporation of a newly developed dynamic relation for the surface roughness length relevant to small-scale dust generation processes. Through this implementation, the effect of nonerodible elements on the local flow acceleration, drag partitioning, and surface coverage protection is modeled in a physically based and consistent manner. Careful attention is paid in integrating the new windblown dust treatment in the CMAQ model to ensure that the required input parameters are correctly configured. To test the performance of the new dust module in CMAQ, the entire year 2011 is simulated for the continental United States, with particular emphasis on the southwestern United States (SWUS) where windblown dust concentrations are relatively large. Overall, the model shows good performance with the daily mean bias of soil concentrations fluctuating in the range of ±1 μg m-3 for the entire year. Springtime soil concentrations are in quite good agreement (normalized mean bias of 8.3%) with observations, while moderate to high underestimation of soil concentration is seen in the summertime. The latter is attributed to the issue of representing the convective dust storms in summertime. Evaluations against observations for seven elevated dust events in the SWUS indicate that the new windblown dust treatment is capable of capturing spatial and temporal characteristics of dust outbreaks.

Introduction

Windblown dust emitted from the earth’s surface to the atmosphere has significant impacts on atmospheric phenomena and consequently on air quality. It contributes to the alteration of solar radiative budgets through absorption and reflection [Sokolik et al., 1998], modification of cloud properties and precipitation [Huang et al., 2006; Rosenfeld et al., 2001], reduction of the atmospheric visibility [Engelstaedter et al., 2003], and long-range transport of organic chemicals, airborne bacterial species, and trace metals [Ren et al., 2014; Li et al., 2013; Maki et al., 2014]. Air quality and human health can be directly affected by the changes in dust emission. Respiratory and cardiovascular disorders, meningococcal meningitis, conjunctivitis, and skin irritations are among the health problems that have been associated with exposure to dust [Goudie, 2014; Crooks et al., 2016]. The fraction of the cardiopulmonary deaths caused by dust is estimated to be approximately 1.8% globally, however this value increases to 15–50% for countries located in the so-called “dust belt” that stretches from North Africa across the Middle East and up to central Asia [Giannadaki et al., 2014].

Considering the strong impact that atmospheric dust can have on weather and climate, air quality, and human health, understanding and quantifying the factors that control windblown dust emissions is important. Following this need, numerous theoretical studies, as well as field and laboratory measurements, have focused on identifying the physical mechanisms behind dust generation through Aeolian (wind) erosion. For instance, Gillette et al. [1980, 1982] and Gillette [1988] conducted a series of measurements over arid, semiarid, and agricultural soils to generate a database for the threshold friction velocity (the surface friction velocity above which soil erosion begins) associated with dust emission. Gillette and Walker [1977] carried out a field measurement study to establish a relationship between vertical dust flux and surface friction velocity. Using wind tunnel measurements allowing control of aerodynamic conditions, Shao et al. [1993] showed that saltation bombardment is mainly responsible for dust emission. Furthermore, they developed a theory based on wind tunnel observations suggesting that the vertical dust flux is proportional to the third power of the friction velocity.

Based on these studies, several mathematical models parameterizing various physical steps in the windblown dust generation have been developed [Marticorena and Bergametti, 1995; Shao et al., 1996; Alfaro and Gomes, 2001; Ginoux et al., 2001; Zender et al., 2003; Kok et al., 2014]. These mathematical parameterizations have been implemented in several global or regional atmospheric/atmospheric-chemistry modeling systems [Pérez et al., 2011; Wang et al., 2012; Astitha et al., 2012; Xi and Sokolik, 2015], and satisfactory predictions have been obtained for many cases. Nevertheless, several studies still demonstrate considerable differences in the dust generation between model predictions and observations, pointing to a need to improve dust emission modeling. Furthermore, intercomparison studies showed that large diversity exists among dust emission models themselves. Uno et al. [2006] compared eight dust emission/transport models over Asia. Although the models captured the onset of dust events, the predicted peak concentration from the models were up to 4 times different. Todd et al. [2008] obtained significant differences (at least an order of magnitude) in estimates of dust concentrations between regional models that use the same dust emission parameterization. A similar degree of uncertainty was reported by Kang et al. [2011] who implemented three different vertical dust flux schemes (with other parameters being kept the same) in the Weather Research and Forecasting with Chemistry (WRF-Chem) modeling system. An intercomparison of 15 global aerosol models as part of the Aerosol Comparisons between Observations and Models (AeroCom) project [Huneeus et al., 2011] showed that the yearly global dust emission flux differed by an order of magnitude between the models.

The high level of uncertainty in modeling dust emission is attributed mainly to two issues. First, the dust emission involves several complex and nonlinear processes that are governed by the meteorology, land surface, and soil texture. It is well-known that the dust flux is closely related to wind speed (friction velocity), surface roughness, soil moisture, vegetation fraction, soil type and texture, and air density. Therefore, small errors in modeling these parameters can result in large dust flux uncertainties. Darmenova et al. [2009] presented an extensive assessment of the sensitivities of important parameters to two dust emission schemes developed by Marticorena and Bergametti [1995] and Shao et al. [1996] (with some more recent improvements) over Central and East Asia. Their results indicated that in both schemes the dust flux is most sensitive to the wind friction velocity, while surface parameters such as roughness and vegetation fraction become increasingly important at low to moderate wind speeds. Another important finding in Darmenova et al.’s [2009] work is that the two schemes produce similar horizontal fluxes over dry smooth surfaces, whereas their differences are large over rough surfaces. The second source of uncertainty is related to the implementation of the dust emission schemes into the atmospheric-chemistry modeling systems. In general, the physically based dust emission schemes have been developed based on detailed analyses of microphysical phenomena in erodible environments to predict saltation flux and the resulting dust emission. Therefore, attention needs to be paid in implementing these fully microphysical schemes into mesoscale models. In particular, the input parameters and physics packages used by mesoscale models should be consistent with the input data required by the dust emission schemes. This necessitates the understanding of the scientific foundation of each step in the parameterization of dust emission and the careful implementation of the scheme in the modeling system.

In this study, we present the development of a new treatment for windblown dust emission in the Community Multiscale Air Quality (CMAQ) [Byun and Schere, 2006] modeling system. The details of the model formulation and integration within CMAQ, as well as evaluation against observations are included. The Version 5.1 of the CMAQ model with the new windblown dust parameterization is used to simulate the entire year 2011 over the continental United States (CONUS), with particular emphasis on the southwestern U.S. over the spring season when several dust storms occurred. Finally, we discuss the uncertainties and limitations of the model and provide several recommendations for future model developments which are of interest to the community.

Physical Parameterization of Windblown Dust Emission

In the present study, the main physical mechanism behind dust emission is considered to be saltation (sandblasting). As shown in Figure 1, wind blowing over a surface results in an aerodynamic force uplifting the sand particles. However, the gravity force on the particles ultimately exceeds the uplifting force, carrying the particle downward back to the surface. This hopping movement of the particles can be seen as a horizontal flux in an average sense. The striking of these returning particles to the surface, known as saltation bombardment, can overcome the interparticle cohesive forces and result in vertical dust flux [Kok et al., 2012]. It has been shown [Loosmore and Hunt, 2000] that the dust flux due to saltation can be orders of magnitude larger than the direct aerodynamic uplifting, i.e., upward transportation of fine particles by turbulent eddies. Therefore, the latter is neglected in the present dust emission scheme.

Figure 1.

Schematic of the mechanism of dust emission by saltation. The forces acting on a stationary sand particle include drag (FD), lift (FL), gravitational (Fg), and cohesive (Fc) forces. The particle starts to move when the aerodynamic drag and lift exceed gravitational and interparticle cohesive forces. The strike of the particle back to the surface results in the vertical emission of dust.

The physics of dust emission due to saltation can be divided into three processes: (1) the onset of saltation by movement of sand particles due to wind, (2) the horizontal flux of sand particles, and (3) the vertical flux of dust due to the impact of sand particles and the surface. All of these processes are clearly affected by the strength of the wind, as well as the surface conditions, including soil texture, soil moisture, and surface roughness. Figure 2 depicts the schematic diagram of the dust emission processes considered in the present model. The details of the model parameterization and formulation are presented in following sections.

Figure 2.

Schematic diagram of the new windblown dust emission treatment in CMAQ.

Threshold Friction Velocity

It has long been known [Bagnold, 1941] that the aerodynamic lifting of surface particles can occur only above a certain level of the wind speed. Therefore, the initiation of saltation is associated with some threshold stresses at the surface. The threshold friction velocity is the minimum surface shear stress (divided by the flow density and square rooted) due to wind, below which the saltation cannot occur. It is a key parameter in dust emission models, and depends not only on the size of the particles, but also on the presence of moisture and/or nonerodible roughness elements such as vegetation and pebbles. For modeling purposes, we consider the threshold friction velocity u*, as follows [Shao et al., 2011]: where u*, is the ideal threshold friction velocity on a dry, smooth, and bare surface, and f and f are the correction factors for soil moisture and surface roughness, respectively. The increase in soil moisture and surface roughness results in protecting the surface from erosion, therefore both f and f should increase the threshold friction velocity, hence they are always greater than or equal to one. It should be noted that the threshold friction velocity can be influenced by other factors, including the salt concentration and surface crustiness. However, due to the lack of adequate information regarding these properties, the correction factors for these parameters are considered to be one [Shao et al., 2011].

The ideal threshold friction velocity u*, can be derived by considering the balance of entraining forces including the aerodynamic drag and lift against the stabilizing forces including the gravity and cohesive inter-particle forces acting upon a particle. Iversen and White [1982] used a combination of theory and wind tunnel measurements, and proposed a semi-empirical implicit formulation for u*, as a function of particle size. Marticorena and Bergametti [1995] presented an explicit form for the formulation of Iversen and White [1982] making it easier to implement in a numerical modeling setup. Shao and Lu [2000], on the other hand, focused on an expression that explicitly treats the cohesive interparticle forces, and proposed the following simple explicit formulation for the ideal threshold friction velocity based on the particle size: In equation (2), D is the sand particle diameter, and ρ and ρ are the particle and air densities, respectively. In the present model, we use equation (2) with the two constants A and Γ (accounting for the magnitude of the interparticle cohesive forces) being set to 0.0123 and 1.65 × 10−4 kg/s2, respectively [Kang et al., 2011]. Equation (2) yields minimum threshold velocities for particle sizes of 60–200 μm (see supporting information Figure S1). At these particle sizes, saltation is most likely to occur. For smaller particles, cohesive interparticle forces prevent saltation, while larger particles are harder to move due to the gravity force. Equation (2) (see supporting information Figure S1) also shows the reason why saltation is the dominant mechanism of dust emission compared to aerodynamic uplifting. Saltation can be initiated at friction velocities well below that required to directly lift the fine particles of a typical size of 1–10 μm. It should be noted, however, that the threshold friction velocity is a stochastic parameter and it is possible for very fine particles to have u*, values as low as medium-sized sand particles [Shao and Klose, 2016].

Soil moisture increases the threshold friction velocity through increasing the cohesive forces suppressing the mobilization of particles. The effect of the increased soil moisture is considered through a correction factor derived by Fécan et al. [1999]: with w being the gravimetric soil moisture. The limit value of the soil moisture w′ depends on the clay content of the soil [Fécan et al., 1999]: The presence of nonerodible elements including vegetation, rock, pebble, and gravel is known to protect the surface from wind erosion through the drag partitioning process. In this process, part of the air flow momentum is extracted by nonerodible elements resulting in less wind shear stress acting on the surface. This effect is taken into account through a surface roughness factor f correcting the threshold friction velocity. Marticorena and Bergametti [1995] developed a drag partitioning parameterization through relating the roughness density (a dimensionless variable defined as λ = nbh/S, where n is the number of nonerodible elements with width b and height h per unit area S) to the aerodynamic roughness length of the surface. Their formulation was later corrected by King et al. [2005]. Nevertheless, it is limited to small solid obstacles only. In this study, we adapt the surface roughness factor based on the double drag partitioning [Raupach et al., 1993; Darmenova et al., 2009], to include the effect of both solid elements (rock, pebble, etc.) and vegetation coverage: In equation (5), λ and λ are roughness densities based on solid (nonvegetation) and vegetation elements, respectively. λ can vary between 0.002 and 0.2 according to measurements of Marticorena et al. [2006]. Here we specify the value of λ for each land type following Darmenova et al. [2009] and Xi and Sokolik [2015] (see Table 1). The roughness density of vegetation elements λ is obtained by the parameterization of Shao et al. [1996]: where A is the vegetation cover fraction. The values of the coefficients in equation (5) are chosen based on the recommendation of Darmenova et al. [2009] and Xi and Sokolik [2015], i.e., σ = 1.45, m=0.16, β=202, σ=1.0, m=0.5, and β=90.

Table 1.

Monthly Vegetation Geometric Height (h), and Geometric Height (h) and Density (λ) of Nonvegetation (Solid) Elements Based on the Land Type

			hV (m)
Land Type	hS (m)	λ S	January	February	March	April	May	June	July	August	September	October	November	December
Shrubland	0.02	0.03	0.05	0.05	0.15	0.15	0.12	0.12	0.10	0.10	0.10	0.05	0.05	0.05
Shrub/grass	0.02	0.04	0.05	0.05	0.05	0.10	0.20	0.15	0.12	0.12	0.10	0.05	0.05	0.05
Barren or sparsely vegetated	0.02	0.002	0.05	0.05	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.10	0.05	0.05
cropland	0.02	0.15	0.05	0.05	0.05	0.05	0.10	0.30	0.50	0.50	0.30	0.10	0.05	0.05

Friction Velocity

The friction velocity, representing the surface shear stress exerted by the wind, is a critical parameter in all dust emission schemes. In fact, it has been shown [Darmenova et al., 2009] that dust models are most sensitive to the value of the friction velocity. It not only dictates the onset of saltation but also affects the total dust flux. Considering the logarithmic wind profile, valid in the surface layer of the planetary boundary layer, the friction velocity can be calculated as: where u is the wind speed at reference height z, κ is the von Kármán constant, and z0 is the surface roughness length. It should be noted that in using equation (7) the planetary boundary layer is assumed to be neutrally stable. The neutral stability assumption is reasonable when turbulent mixing overpowers atmospheric instability which is often the case during dust events. When using dust parameterization schemes in mesoscale atmospheric models (e.g., WRF), the value of u is obtainable from the model calculations. Here we use the value of modeled 10 m wind speed, hence z = 10 m. Therefore, the only remaining parameter to obtain the friction velocity is the roughness length z0. In the next section, we discuss the surface roughness length in details.

It is important to note that while the friction velocity u* is calculated and available in mesoscale atmospheric models, we argue, similar to Darmenova et al. [2009] and Menut et al. [2013], that this modeled friction velocity should not be used in dust emission calculations. Instead, a u* value applicable to the dust emission scheme should be recalculated using equation (7) with a surface roughness length z0 relevant to the physics of the dust generation. The friction velocity in the meteorological models is calculated using a predefined aerodynamic roughness length (often in the form of a table for various land types), which is significantly higher than the “microscopic” length relevant to dust generation. In fact, Menut et al. [2013] showed that there is an order of magnitude difference between these two values in a domain including North Africa, the Middle East, and East Asia. Therefore, using this value in a dust emission model would result in unrealistic predictions.

Surface Roughness Length

It was previously discussed that nonerodible elements protect the surface from erosion by damping the wind. However, they also produce an opposing effect which increases the total surface shear stress acting on both roughness elements as well as the ground. The physics behind this effect can be explained by considering the wake region behind each element resulting in a local acceleration of flow and a higher shear stress. Integrating over the whole surface, the net effect is to increase the friction velocity as represented mathematically by equation (7). The surface roughness length z0 is the parameter that reflects the net effect of nonerodible elements on increasing the friction velocity and depends on the height and density of the roughness elements. Previous models used different approaches to address the roughness length. Some studies [Zender et al., 2003] assumed a constant value for z0, while others used values retrieved from observations [e.g., Laurent et al., 2006, who used POLDER data] or empirical formulations [e.g., Zhao et al., 2006, who used the Marticorena et al.’s, 1997 relationship].

In the present study, we collected the surface roughness length data from several previous field and laboratory measurement studies including those in King et al. [2005], Marticorena et al. [2006], and Hébrard et al. [2012], in order to generate a comprehensive database of surface roughness lengths. These data are provided in Figure 3, which shows variations of the roughness length, normalized by the roughness height, along with the roughness density. From the figure, it is evident that the effective roughness length increases up to a point with increasing roughness density (number of roughness elements). However, at excessive roughness densities this parameter decreases. The physics involved can be described using Figure 4. For a smooth surface (Figure 4a), drag is minimal. Once a moderate number of roughness elements are introduced at the surface (Figure 4b), drag increases due to the wake region that forms behind each element. The number and size of these wakes determine the increase in drag acting on the surface. Nevertheless, introducing too many elements (Figure 4c) may result in decreasing the number and size of the wake regions since they are being filled by the elements. The net effect of this is a decrease in the effective roughness length, as suggested by Figure 3. A correct representation of these phenomena may be important in predicting dust emission over vegetated surfaces such as shrubland and grassland which are dominant in western U.S. drylands, where dust storms are common.

Figure 3.

The roughness length (normalized by the mean roughness elements height) as a function of the roughness density collected from several field and laboratory data summarized by King et al. [2005], Marticorena et al. [2006], and Hébrard et al. [2012]. The fitted line represents equation (8).

Figure 4.

Coherent structures in turbulent flow over smooth, moderately rough, and very rough surfaces.

From a linear fitting of the data in Figure 3, the surface roughness length relevant to the Aeolian processes is estimated as: Unlike Marticorena et al.’s [2006] formula that keeps the value of z0/h constant above a critical λ value, equation (8) accounts for the decrease in roughness length in highly covered surfaces. Additionally, the critical λ value (0.045 in Marticorena et al. [2006]) is found to be around 0.2, which is in agreement with the finding of Shao and Yang [2005]. Xi and Sokolik [2015] stated that the value of 0.045 causes too high z0 in their study, and used the formula of Marticorena et al. [2006] with the critical λ value of 0.1.

In equation (8), λ is the total roughness density, i.e., λ=λ + λ, and h is the effective height of roughness elements which is calculated by weighting the height of solid (h) and vegetation (h) elements based on their densities as: Here we consider a temporally variable surface roughness. In this approach, the density and height of solid elements λ and h are considered to be time-invariants, and a predefined value for each land type is used [see Pierre et al., 2012, for another approach]. The vegetation roughness, on the other hand, is considered to be dynamic, and vegetation height h is changed monthly to take the effect of vegetation growth/decay into account. The density of vegetation elements λV is also temporally variable and inherits its dynamic nature from temporal changes in the vegetation cover fraction according to equation (6) (discussed further in section 3). Table 1 summarizes the values of h, λ, and h used in the present model, following the studies of Marticorena et al. [2006], Darmenova et al. [2009], and Xi and Sokolik [2015]. It should be noted that all of these parameters are treated the same when calculating the surface roughness factor f to make sure that there exists a consistency between the roughness length in calculating the friction velocity and the surface roughness correction factor presented by equation (5).

Dust Flux

Once the friction velocity exceeds a threshold value, sand particles may be lifted from the surface and moved in ballistic trajectories [Kok and Reno, 2009]. The mass flux of these particles, known as the horizontal or saltation flux, can be computed using White’s [1979] relation. For each sand particle size D: where ρa is the air density and the constant c is set to 1.0 [Darmenova et al., 2009].

The total horizontal flux for all particles is calculated as: where S(D) is the relative surface area covered with particles with diameter D. To calculate S(D), we assume sand particles are spherical and each soil type consists of coarse sand, fine-medium sand, silt, and clay. Each of these four soil populations is associated with a mean mass diameter D and a standard deviation σ following the study of Menut et al. [2013]. Table 2 summarizes the composition of each soil type and size-distribution characteristics of each soil population.

Table 2.

Composition of Each Soil Type

		%
Number	Soil Type	Coarse Sand D = 690 μm σ = 1.6	Fine-Medium Sand D = 210 μm σ = 1.6	Silt D = 125 μm σ = 1.8	Clay D = 2 μm σ = 2.0
1	Sand	46	46	5	3
2	Loamy sand	41	41	18	0
3	Sandy loam	29	29	32	10
4	Silt loam	0	17	70	13
5	Silt	0	10	85	5
6	Loam	0	43	39	18
7	Sandy clay loam	29	29	15	27
8	Silty clay loam	0	10	56	34
9	Clay loam	0	32	34	34
10	Sandy clay	0	52	6	42
11	Silty clay	0	6	47	47
12	Clay	0	22	20	58

Upon hitting the surface, the saltating particles may splash and eject dust particles into the air, a process known as saltation bombardment or sandblasting. The mass flux of dust particles generated by this process, known as the vertical flux, is considered to be proportional to the horizontal flux [Shao et al., 1993; Marticorena and Bergametti, 1995]: where α is the so-called vertical-to-horizontal flux ratio or sandblasting efficiency. Marticorena and Bergametti [1995] derived an empirical formulation for α based on the data set of Gillette [1979]. However, this formulation is valid only for soils with clay content of up to 20%. Furthermore, Kang et al. [2011] showed that the vertical flux parameterization of Marticorena and Bergametti [1995] may result in considerable over-prediction of the dust emission.

In this study, we adopt a physics-based relation for the vertical-to-horizontal flux ratio following Lu and Shao [1999]. In this approach, the sandblasting efficiency is determined based on the small volume of the soil removed by impacting an individual sand grain, which can be obtained by solving the equations of particle motion. Through a detailed analytical solution together with a set of simplifying assumptions, Lu and Shao [1999] derived the following equation for the vertical-to-horizontal dust flux ratio: In equation (13), f is the fraction of fine particles contained in the soil volume, p is the surface deforming property or plastic pressure (large p corresponds to a hard surface), ρ and ρ are the bulk soil and soil particle densities, and Cα and Cβ are constants. All these parameters are functions of soil type. We implement this parameterization in our windblown dust emission model with all the parameters being set following the study of Kang et al. [2011].

According to equation (13), the sandblasting efficiency is not only a function of the fine-particle content in the soil, but also a function of the surface hardness and the impact velocity. For instance, a soil with high clay content has a large value of f and a potential to generate more dust, however, the surface plasticity p is also high suppressing the dust emission. Furthermore, a stronger wind (higher u*) results in more energetic particles impacting the surface and elevated dust generation. The scheme of Marticorena and Bergametti [1995] does not involve these physical processes and, therefore, presents limitations.

It should be noted that the inherent limitation of equations (12) and (13), as well as any other bulk dust emission schemes, e.g., Marticorena and Bergametti [1995], Alfaro and Gomes [2001], and Zender et al. [2003], is the lack of size information of emitted dust. Size-resolved schemes [Shao et al., 2011] have been developed to address this issue and satisfactory results have been obtained. The challenge, however, is to correctly represent the soil particle size data including the shape factor of the size distribution.

Model Implementation, Input, and Configuration

The physics-based windblown dust scheme discussed in section 2 is implemented in version 5.1 of the CMAQ modeling system (v5.1) [Appel et al., 2016]. CMAQ is developed and maintained by the U.S. Environmental Protection Agency’s (EPA) Office of Research and Development (ORD). CMAQ v5.1 was publically released in November 2015 (http://www.cmaq-model.org), with numerous fixes, enhancements, and scientific updates in meteorology and transport, aerosol treatment, in-line photolysis and cloud modeling, atmospheric chemistry, and air-surface exchange. The details regarding new developments in CMAQ v5.1, model setup, and a comprehensive evaluation can be found in the release notes available along with the model source code as well as in Appel et al. [2016].

To obtain the total dust emission in each computational grid cell in CMAQ, equation (13) is integrated over all erodible land types: where A is the relative surface area of the grid cell covered by the land type L. Erodible lands are prescribed based on the land use data set being utilized. Here we use the Biogenic Emissions Landcover Database version 3 (BELD3) that combines the U.S. Geological Survey (USGS) data with the detailed tree and crop species information available in county-level forest and agricultural data sets to provide 1 km resolution for 230 different land use types [Vukovich and Pierce, 2002]. In the framework of BELD3, we consider USGS_shrubland, USGS_shrubgrass, USGS_cropland, and USGS_sprsbarren land use categories as erodible lands. The coefficient 1–A in equation (14) is included to account for the portion of the surface covered (and thus protected from the erosion) by vegetation.

The total dust emission is then distributed to four size bins of 0.1–1.0, 1.0–2.5, 2.5–5.0, and 5.0–10.0 μm. The first two bins are considered to be the fine mode with the geometric mean diameter of 1.3914 μm (and the standard deviation of 2.0), and the second two bins are considered to be the coarse mode with the geometric mean diameter of 5.2590 μm (and the standard deviation of 2.0). The mass distribution between fine and coarse mode should be determined. In this study, we distinguish between the freshly emitted dust and the transported dust in terms of the mass distribution, with the former being relevant in our model. The observations [Sugimoto et al., 2016] revealed that the fine-to-coarse ratio for local dust (freshly emitted) is less than 0.1 while this ratio ranges from 0.7 to 0.8 for the transported dust. Kok [2011] developed a theoretical expression for the emitted dust size distribution supported by the size-resolved measurements, and suggested that the previous size distributions of emitted dust aerosols used in the global circulation models overestimate the fraction of fine dust aerosols. Nabat et al. [2012] used the dust size-distribution scheme based on Kok [2011] in a regional climate model and showed a clear improvement. In this study, we split the generated windblown dust mass to 7% and 93% for the fine and the coarse mode, respectively, following the study of Nabat et al. [2012].

The soil type data are taken from the US State Soil Geographic (STATSGO) soil database [Miller and White, 1998]. Following Tegen et al. [2002], each soil type is considered to be a combination of four soil textures, namely coarse sand, fine-medium sand, silt, and clay as summarized in Table 2.

Correctly representing the spatial and temporal variations in surface vegetation is important due to its various effects on dust generation including drag partitioning, local wind acceleration, near-source removal, and protective coverage. In the present model, the fraction of absorbed photosynthetically active radiation (FPAR) based on Moderate Resolution Imaging Spectroradiometer (MODIS) observations are used to represent more realistic time-varying vegetation coverage. The global MODIS FPAR data in the MOD15A2 product at 1 km resolution and every 8 days [Myneni et al., 2011] are processed and regridded onto the computational domain. The gridded cell-averaged MODIS FPAR is used as a surrogate for vegetation cover fraction A [Mu et al., 2011; Ran et al., 2015]. Ran et al. [2016] incorporated MODIS products including FPAR, leaf area index (LAI), and albedo into the Pleim-Xiu land surface model (PX LSM) [Pleim and Xiu, 1995], and showed reduction in 2 m temperature, mixing ratio, and 10 m wind speed biases. Supporting information Figure S2 shows the monthly variation of the vegetation fraction plotted at the 15th of each month of the year 2011.

The required meteorological inputs for the present dust scheme and the CMAQ modeling system are provided by WRF v3.7 simulations with the physics packages and computational setup described in details in Appel et al. [2013] and Ran et al. [2016]. Specifically, the air density (in equation (10)), 10 m wind speed (in equation (7)), and soil moisture (in equation (3)) are directly used from WRF simulations. Attention should be paid in using the meteorologically calculated soil moisture in regional and global models as an input required in the dust emission scheme. Since windblown dust generation is a surface process, only the moisture of the active “topmost layer” matters. In fact, the empirical equation of Fécan et al. [1999] (equation (3)) that quantifies the effect of the soil moisture on the threshold friction velocity has been developed based on the wind tunnel studies in which the moisture representative of the topmost ~1 cm of the soil was changed. The “surface layer” in the regional or global models, however, can vary widely based on the land surface model (LSM) being used. The integrated soil moisture in this (sometimes thick) “surface layer” can differ significantly from the “topmost layer” relevant to dust generation. For instance, Darmenova et al. [2009] indicated that the soil moisture calculated by the Noah LSM [Chen and Dudhia, 2001] (with the “surface layer” being 10 cm) is quite high resulting in complete suppression of the dust emission. They suggested a factor of 0.1 to be used to reduce this calculated soil moisture. In this study, the WRF simulation is configured with the PX LSM which has the top soil layer with 1 cm thickness. Therefore, the meteorologically calculated value of soil moisture should be appropriate to the dust generation processes. Another issue regarding the soil moisture is the actual drying time in drylands, which is believed [Myhre et al., 2003] to be shorter than what LSMs predict. Consequently, some studies [Myhre et al., 2003; Grini et al., 2005] used hourly or daily precipitation thresholds as dust suppressing parameters. However, considering the lack of detailed evaluation and justification of these methods, we use the calculated value of the surface layer soil moisture by the PX LSM.

The CMAQ v5.1 model with the new windblown dust emission treatment is used to simulate the entire year 2011 on a domain that covers the contiguous United States (CONUS), as well as portions of Mexico and Canada with 12 km horizontal resolution and 35 vertical layers extending up to 50 hPa. The CMAQ is configured with the updated Carbon Bond mechanism (CB05e51) for gas-phase chemistry, the AERO6 aerosol module, in-line photolysis calculation, the bi-directional ammonia flux (bi-di) option for estimating the air-surface exchange of ammonia, and NO emissions from lightning (see http://www.airqualitymodeling.org/cmaqwiki for more details). The base emissions are based on version 2 of the 2011 National Emissions Inventory (NEI; www.epa.gov/ttnchie1/net/2011inventory.html). The input meteorological fields for CMAQ are provided by running version 4.2 of the Meteorology-Chemistry Interface Processor (MCIP) [Otte and Pleim, 2010] for the WRF v3.7 outputs. Lateral boundary conditions for the present CMAQ simulations are obtained from an annual 2011 GEOS-Chem [Bey et al., 2001] simulation. The CMAQ simulations are conducted from 21 December 2010 (with clean initial conditions) to 31 December 2011 with the first 10 days being considered spin-up (not included in the analysis) to minimize the effects of initial conditions on simulation results.

Model Evaluation

Observational Data and Evaluation Methodology

There are several national or regional monitoring networks in the United States that provide routine observations of particulate matter (PM). In this study, we utilize data from the Interagency Monitoring of Protected Visual Environments (IMPROVE; http://vista.cira.colostate.edu/improve/) network to evaluate the newly developed dust model. The daily averaged values of total PM2.5 and PM10, along with the chemical composition of PM2.5 are available at IMPROVE sites every third day. Here we consider the fine “soil,” as defined by the IMPROVE network, as the relevant parameter in our evaluations. It is calculated based on the mass concentrations of five most important mineral elements in dust. Specifically, the soil concentration [Soil] is calculated as: where [Al], [Si], [Ca], [Fe], and [Ti] are the measured concentrations of fine (PM2.5) particulate Aluminum, Silicon, Calcium, Iron, and Titanium, respectively. The IMPROVE network sites are typically located in rural areas with several of the sites being close to the major dust sources in the United States, and therefore the values of soil from these sites are good indicators of dust outbreaks. To study episodic and individual dust outbreaks where a higher temporal resolution of PM concentration observations is needed, the EPA’s Air Quality System (AQS; https://www.epa.gov/aqs) network is considered where hourly and daily PM2.5 and PM10 concentrations are measured.

Output from the CMAQ simulations with the new dust model is paired in space and time with observed data using the Atmospheric Model Evaluation Tool (AMET) [Appel et al., 2011] and the corresponding statistics are calculated. Specifically, the mean bias (MB), normalized mean bias (NMB), and median bias (MdnB) are calculated as: where [M] and [O] are modeled and observed concentrations, respectively, and N is the total number of model-observation pairs.

Overall Performance

Figure 5 compares the time series of observed and modeled soil concentrations averaged over all IMPROVE sites in the CONUS. The highest observed soil concentrations are seen in the springtime (March–May), with a peak reaching 2.2 μg m−3. The CMAQ-modeled soil concentrations are in good agreement with the IMPROVE observations, with a mean bias (MB) in the range of ±1 μg m−3 throughout the year. Particularly, the new windblown dust model shows a relatively good performance in capturing the occurrence of dust outbreaks and the daily variation of the soil concentration, although there are some over/under-predictions. The model, however, is not able to predict the peak level of soil concentrations in early July, which is attributed to the inability of the model to simulate the occurrence of the severe convective dust storm, so-called haboob, in Phoenix, AZ on 5 July 2011. Further discussion of this issue is provided in section 5. The domain-wide and annual MB is 0.11 μg m−3 (with a MdnB of 0.12 μg m−3) over 14,573 observations from 134 IMPROVE sites in 2011.

Figure 5.

(top) Time series of fine particulate soil and (bottom) the model bias (MB) for all IMPROVE sites in the CONUS in 2011.

To further investigate the performance of the model, spatial plots of seasonal MB for soil concentration are presented in Figure 6. It is seen that for the majority of the sites the MB for soil is less than ±0.5 μg m−3 over the entire year. There is a systematic moderate overprediction of the soil concentrations in the Midwest United States (37°N-44°N, 84°W-103°W). The airborne soil in this region is primarily due to anthropogenic sources (based on NEI) such as unpaved roads, commercial construction, and agricultural tilting, with a smaller contribution from windblown dust [Appel et al., 2013]. This suggests that the anthropogenic dust may be overestimated in the CMAQ model. The Southwest United States, on the other hand, is greatly affected by windblown dust, and this is where the greatest variability in seasonal bias is seen. Specifically, modeled soil concentration in this region varies between being slightly positively biased in the spring (March-May), moderately to highly negatively biased in the summer (June-August), slightly negatively biased in the autumn (September-November), and slightly positively biased in the winter (December-February). Table 3 summarizes the seasonal statistics of the soil concentrations in the Southwest United States obtained from the present simulations with the new windblown dust model.

Figure 6.

Soil seasonal mean bias (μg m−3) for the IMPROVE network sites in the CONUS.

Table 3.

Seasonal Statistics for Soil Concentrations in the Southwestern United States (SWUS)

	Number of Observations	Mean Observed (μg m−3)	Mean Modeled (μg m−3)	MB (μg m−3)	NMB (%)
Spring	1029	1.66	1.80	0.14	8.3
Summer	1065	1.52	0.52	−1.01	−66.1
Autumn	1014	0.81	0.57	−0.24	−29.5
Winter	993	0.46	0.56	0.11	23.3
Total	4101	1.12	0.86	−0.26	−23.2

The Southwest United States in Spring

The domain-wide performance of the model over the entire year 2011 was examined in the previous section. The dust outbreaks, however, are episodic, i.e., they generally occur over a limited region and on relatively small temporal scales (hours to days). Therefore, the focus of further discussion in this study will be on a subset region of the CONUS covering the Southwest United States (SWUS). The area of interest (22°N-43°N, 101°W-116.5°W), shown in Figure 7, is chosen based on the findings of several previous studies [Malm et al., 2004; Wells et al., 2007; Tong et al., 2012] on the geographical pattern of active dust sources in the United States. This region is associated with strong winds over arid and semiarid lands where the natural windblown dust emission is a major factor affecting the atmosphere’s composition and air quality. As predicted by climate models [Seager et al., 2007], a transition to a more arid climate is under way in the SWUS, therefore, correctly simulating the windblown dust emission is of high importance in this region. The focus here will be on the springtime, when high wind speeds, low soil moisture, and a lack of vegetation coverage over erodible land surface result in peak windblown dust emissions. The springtime maximum of the windblown dust emission over the SWUS is reported by several multi-year studies, including Tong et al. [2012].

Figure 7.

The southwestern United States (SWUS) domain investigated in this study. The background is the percentage of the CMAQ grid cell covered by erodible land indicating the potential for the windblown dust generation.

The model performance for PM2.5 and PM10 at IMPROVE sites in the SWUS during spring 2011 is summarized in Figure 8 and Table 4. The data are from 39 IMPROVE sites in the SWUS providing 1029 and 1074 daily observations of PM2.5 and PM10, respectively. To isolate the influence of the new windblown dust model, results of a baseline simulation with no windblown dust (and all other configurations being the same) are also included in the analysis. The “no dust” simulation underestimates the concentration of PM2.5 and PM10 by −32.7% and −60.1% respectively, indicating the importance of including a windblown dust parameterization in the model. For the simulation that includes the new dust model, estimates of PM2.5 and PM10 improve dramatically, with NMB values of −10.5% and −19.2% for PM2.5 and PM10, respectively.

Figure 8.

CMAQ evaluation against PM2.5 and PM10 observations during spring 2011 at IMPROVE sites in SWUS for simulations with the new dust model and no dust.

Table 4.

Comparison of PM2.5 and PM10 Observations and Models (New Dust and No Dust) Over the SWUS in Spring 2011

		Number of Observations	Mean Observed (μg m−3)	Mean Modeled (μg m−3)	MB (μg m−3)	NMB (%)
PM2.5	No dust	1029	4.2	2.8	−1.37	−32.7
	New dust			3.7	−0.44	−10.5
PM10	No dust	1074	12.6	5.0	−7.58	−60.1
	New dust			10.2	−2.42	−19.2

The contribution of soil to the concentration of total PM2.5 is examined in Figure 9, where the observed and modeled stacked bar plots of PM2.5 species composition from 37 IMPROVE sites (total of 903 daily observations) for the SWUS are illustrated. Note that “no dust” simulations still contain soil from the boundaries as well as contributions of nonwindblown dust emissions of Aluminum, Silicon, Calcium, Iron, and Titanium (see equation (15)). The fine particulate soil makes up the largest component of PM2.5 aerosol in the SWUS in spring, composing nearly 40% of the total PM2.5 aerosol mass concentration at IMPROVE sites. After incorporating the new windblown dust treatment, the predicted level of soil concentration is in much better agreement with the observations, resulting in an overall improvement in PM2.5 predictions. Figure 9 also suggests that an appropriate windblown dust module can help improve predictions of other aerosol constituents such as nitrate () and salt (NaCl) by more accurate representation of their thermodynamic portioning.

Figure 9.

Stacked bar plots of PM2.5 compositions at IMPROVE sites in SWUS; comparison between observations, CMAQ simulations with the new windblown dust parameterization, and CMAQ simulations with no windblown dust.

To evaluate the capability of the new windblown dust treatment in CMAQ to capture individual dust outbreaks, we use dust storms as reported in the U.S. National Weather Service (NWS) storm database (https://www.ncdc.noaa.gov/stormevents/faq.jsp). Following the methodology of Crooks et al. [2016] we select seven major dust storms that occurred over the SWUS during spring 2011 for a more detailed analysis. Descriptions of the events are extracted from the NWS storm database. Windblown dust emission is shown [Darmenova et al., 2009] to be most sensitive to wind speed, therefore before investigating CMAQ results, we evaluate WRF predictions of the 10-m wind speed. Figure 10 compares 10 m wind (m/s) from our WRF simulations and observations from the Meteorological Assimilation Data Ingest System (MADIS; https://madis.noaa.gov/) data set. For each of the seven dust storms (Figures 10a–10g), the time series of observed and modeled 10 m wind speed averaged over MADIS observation sites in the region are plotted. It can be seen that WRF simulations generally agree well with the observations in all seven cases. For specific hours, however, the peak wind speed may be over-/underpredicted by more than 20% (see Figures 10b and 10d), which should be considered as a source of uncertainty in the present dust modeling results.

Figure 10.

Time series of 10 m wind speed (m s−1) obtained from WRF simulations (Red) in comparison with observations from MADIS sites (Black). The values are averaged over (a) KLBB, KINK, KODO, (b) KABQ, KELP, KCNM, KSAF, KRTN, (c) KABQ, KELP, KCNM, KSAF, KRTN, KVEL, KSLC, KTUS, (d) KABQ, KELP, KCNM, KSAF, KRTN, (e) KABQ, KELP, KCNM, KSAF, KRTN, (f) KINW, KFLG, (g) KABQ, KELP, KCNM, KSAF, KRTN, KTUS, KPHX, KSDL, KDVT sites. See http://www.cnrfc.noaa.gov/metar.php for information about each site.

27 February 2011

Strong midlevel winds at the border of Texas and New Mexico on 27 February 2011 resulted in widespread blowing dust across the region. The dust storm is reported to have begun at 18:00UTC, and continue for 6 h (Figure 10a). The visibility at Lubbock’s Preston Smith International Airport (33.66°N, 101.82°W) dropped to less than one-quarter mile, while Van Horn (31.04°N, 104.83°W) and Pecos (31.42°N, 103.49°W) in Texas experienced nearly zero visibilities in the storm. This resulted in an accident involving several vehicles. Figure 11 shows contours of hourly PM10 concentrations averaged from 18:00UTC 27 February to 00:00UTC 28 February, overlaid with the locations of Lubbock, Pecos, and Van Horn in Texas. The new windblown dust model implemented in CMAQ is able to qualitatively capture the time and location of the dust outbreak. Comparing the concentrations of modeled daily averaged PM2.5 with observation from an AQS site located in Lubbock (not shown here) shows the model underestimated PM2.5 by 43%.

Figure 11.

Average predicted hourly PM10 surface concentrations (μg m−3) during 18:00UTC 27 February 2011 and 00:00UTC 28 February 2011. Also, the locations of Lubbock (black), Pecos (blue), and Van Horn (green), where very low visibilities due to dust were reported by the NWS.

7 March 2011

A large low-pressure system produced strong winds across a large portion of northern Mexico, far west Texas, and southern New Mexico on 7 March 2011. During 22:00UTC 7 March to 03:00UTC 8 March, the peak wind gust at Apache Point (32.78°N, 105.82°W) and White Sands Missile Range (32.24°N, 106.35°W) in New Mexico and El Paso (31.76°N, 106.46°W) in Texas reached 63 mph, 71 mph, and 66 mph, respectively. As a result, a large area of blowing dust was generated, affecting the air quality in New Mexico and western Texas. The results of the new dust model are depicted in Figure 12a where a large region of high PM10 concentrations is predicted over New Mexico and western Texas. In order to quantitatively examine the model performance in capturing the temporal variation of the dust storm, hourly values of PM2.5 and PM10 concentrations averaged over the 12 sites shown in Figure 12a are plotted in Figures 12b and 12c, respectively. Predictions from the new dust model are compared with the observations from AQS sites along with similar estimates from the “no dust” simulation. The model agrees quite well with the averaged concentrations over all sites, both in terms of the time and level of the peak PM2.5 concentrations. The model, however, clearly underpredicts the peak PM10 averaged concentration. The concentrations of PM2.5 and PM10 remain unchanged in the “no dust” simulation, indicating that natural windblown dust is indeed the cause of the peak PM in the simulation with the new dust model.

Figure 12.

(a) Average predicted hourly PM10 surface concentrations (μg m−3) during 22:00UTC 7 March 2011 and 03:00UTC 8 March 2011, overlaid with the locations of 12 AQS sites (black dots) in NM, (b) hourly concentrations of PM2.5, and (c) PM10 averaged for these 12 sites.

21 March 2011

Very strong midlevel winds overspread across New Mexico, eastern Utah, and northern Arizona on 21 March 2011. Sustained winds over 40 mph were reported at Clines Corners (35.01°N, 105.67°W) in New Mexico and Winslow Airport (35.02°N, 110.71°W) in Arizona for four hours, resulting in blowing dust and reduced visibility below two miles. The NWS reported the event started around 19:00UTC 21 March and lasted for more than 5 h (Figure 10c). As shown in Figure 13, the model captures the geographical extent of the dust storm. The hourly variation in the spatially averaged PM concentrations agrees well with observations averaged over all the AQS sites located in the region (32°N–41°N, 105°W–112°W), while the peak values of PM2.5 and PM10 are overestimated and underestimated, respectively. This issue, in general, may be attributed to the calculation of the deposition velocity in CMAQ. Further investigations are currently underway.

Figure 13.

(a) Average predicted hourly PM10 surface concentrations (μg m−3) during 19:00UTC 21 March 2011 and 00:00UTC 22 March 2011, (b) hourly concentrations of PM2.5, and (c) PM10 averaged for the AQS sites (black dots) indicated.

3 and 9 April 2011

The NWS reported continuous strong winds from 16:00UTC 1 April to 00:00UTC 10 April over New Mexico. These strong winds, generated over the Great Basin Desert and approaching New Mexico from west, reached 60 mph in speed and resulted in blowing dust that reduced the visibility. Visibilities of one-quarter mile were reported at the Deming airport (32.27°N, 107.72°W), while several highways in the region were closed for parts of the day due to significant blowing dust. Two major dust storms in this period occurred on 3 and 9 April 2011 (Figures 10d and 10e). Figures 14 and 15 show PM10 concentrations predicted by the model averaged for the hours of 17:00UTC 3 April to 06:00UTC 4 April, and 17:00UTC 9 April to 09:00UTC 10 April, respectively. There are 12 active AQS sites in the region for this period as shown in Figures 14a and 15a. Comparing the time series of observed and modeled hourly PM2.5 and PM10 averaged over these sites (Figures 14b, 14c, 15b, and 15c) shows that the CMAQ model with the new windblown dust scheme is able to capture the occurrence of these dust events and the associated averaged PM concentrations quite well.

Figure 14.

(a) Average predicted hourly PM10 surface concentrations (μg m−3) during 17:00UTC 3 April 2011 and 06:00UTC 4 April 2011, overlaid with the locations of 12 AQS sites (black dots) in NM, (b) hourly concentrations of PM2.5, and (c) PM10 averaged over these 12 sites.

Figure 15.

(a) Average predicted hourly PM10 surface concentrations (μg m−3) during 17:00UTC 9 April 2011 and 09:00UTC 10 April 2011, overlaid with locations of 12 AQS sites (black dots) in NM, (b) hourly concentrations of PM2.5, and (c) PM10 averaged over these 12 sites.

18 May 2011

An unseasonably strong and cold low-pressure system brought strong winds to northeastern Arizona between 17:00UTC and 21:00UTC on 18 May 2011. The resulting blowing dust along interstate highway 40 (I-40) between Flagstaff (35.20°N, 111.65°W) and Winslow (35.02°N, 110.70°W) reduced visibility to the point where the Department of Public Safety began escorting vehicles along a 5 mile stretch of the freeway. Considering the small spatial and temporal scales of this dust storm, the new dust model shows relatively good performance in capturing the time and location of the event (Figure 16). No observational site is available in the region for a further quantitative assessment of the model performance.

Figure 16.

Average predicted hourly PM10 surface concentrations (μg m−3) during 17:00UTC and 21:00UTC on 18 May 2011, together with the locations of Winslow (blue) and Flagstaff (black) where very low visibility was reported by the NWS.

29 May 2011

Another low-pressure system produced high winds and blowing dust over northern Arizona and New Mexico on May 29, 2011. The visibility was reported to be reduced to 25 ft at Many Farms (36.35°N, 109.62°W) in Arizona. Figure 17 depicts the contours of hourly PM10 concentrations averaged from 20:00UTC 29 May to 03:00UTC 30 May overlaid with the locations of 14 AQS sites in the region (32°N–40°N, 105°W–115°W) where hourly observations of PM2.5 and PM10 concentrations were available. Overall the model simulation with the new dust scheme compares well with the observations (Figures 17b and 17c) in an average sense, but misses the peak observed PM concentrations. However, the performance is greatly improved over the simulation without any windblown dust treatment.

Figure 17.

(a) Average predicted hourly PM10 surface concentrations (μg m−3) during 20:00UTC 29 May 2011 and 03:00UTC 30 May 2011, (b) hourly concentrations of PM2.5, (c) PM10 averaged for the AQS sites (black dots) indicated.

Uncertainties, Remaining Issues, and the Path Forward

There are various sources of uncertainties in windblown dust emission parameterization schemes, and it is important to identify these sources to understand the limitations of the model, as well as to determine necessary future developments. The physics-based schemes, in general, require several inputs, and therefore, any uncertainties in the input parameters would propagate through the model predictions. The newly developed windblown dust emission scheme is closely tied to the meteorological fields. Specifically, the wind speed, soil moisture, and air density are directly used from the meteorological driver, i.e., WRF. Hence, dust emission calculations are sensitive to the accuracy of WRF simulations (see Figure 10 and supporting information S3 for evaluations of WRF results in this study). Sensitivity studies [Menut, 2008; Darmenova et al., 2009] have concluded that dust emission models are most sensitive to the wind friction velocity. Through definition of a normalized sensitivity coefficient, Darmenova et al. [2009] showed that for a 10% change in the friction velocity, the total horizontal flux can change by more than 100%. Therefore, continuous improvement of the meteorological models is crucial for better simulation of the windblown dust emission.

Another key uncertainty relates to the fact that the results from regional and global atmospheric models are resolution-dependent [Liu and Westphal, 2001]. Specifically, the processes that are nonlinearly dependent on meteorological fields are affected by the resolution of the simulation due to subgrid variability. For instance, the present windblown dust model uses the grid-averaged values of wind speed, soil moisture, soil type, and vegetation coverage fraction (all nonlinearly contributing to the calculations of dust emission flux) and, therefore, results would be sensitive to the size of the grids. Attempts have been made to develop resolution-independent dust emission schemes [Ridley et al., 2013], but these models require initial high-resolution simulations to derive the required coefficients.

Finally, it should be noted that windblown dust emission parameterizations implemented in large-scale weather and climate models mainly represent the dust emission due to synoptic scale winds. Nevertheless, it has been shown that there are other processes in the atmosphere associated with dust emission, including cold pools of convective systems [Knippertz et al., 2009] and microscale dust devils [Klose and Shao, 2016]. In particular, moist convection is shown to be an important contributor to the atmospheric windblown dust, especially in the summertime [Heinold et al., 2013; Pantillon et al., 2016]. The “large-scale” atmospheric models, as presented herein, incorporate parameterization schemes for subgrid convection rather than explicitly representing the cold pool and its propagation and, therefore, lack the wind speeds required to generate convective dust storms. The same conclusion has been reached by several other studies [Garcia-Carreras et al., 2013; Heinold et al., 2013; Sodemann et al., 2015]. The most prominent limitation of the model in this context is the inability to capture the severe convective dust storm, also known as haboob, that hit Phoenix, AZ in the early evening hours of 5 July 2011 (see Figure 5). This complex event, caused by a severe thunderstorm and associated downburst winds, can only be accurately represented through high-resolution convection-permitting simulations [Vukovic et al., 2014].

The above mentioned limitations highlight the challenges that should be addressed in future studies. In general, attempts should be made to move away from “tuning” toward more physical representation of small-scale meteorological and land surface processes in operational models. Specifically, the subgrid scale wind variability and the role of convective processes [Cakmur et al., 2004; Klose et al., 2014; Zhang et al., 2016] should be represented in a windblown dust parameterization. We are currently working on implementing and testing subgrid variability of surface wind through probability density functions in our model, with the ultimate goal being the development of a scale-aware windblown dust parameterization with reduced sensitivity to grid resolution.

Conclusion

In this work, we presented the incorporation of a new windblown dust emission treatment in the CMAQ modeling system. The new scheme, built upon the physics-based models from the literature, incorporates the effects of wind speed, soil texture, soil moisture, and surface roughness in a physically sound manner. Particularly, we developed a new formulation for the surface roughness length relevant to the small-scale Aeolian processes by collecting surface roughness length data from several previous field and laboratory measurements. We argued that the WRF-modeled friction velocity should not be used as an input to the dust emission parameterization since it represents the surface roughness at a different scale. Instead, the friction velocity should be recalculated with a “microscopic” roughness length relevant to the dust generation phenomena. We linked this newly developed formulation to the MODIS FPAR observations to represent the dynamic nature of the surface roughness length. In the new scheme, the effects of solid and vegetative nonerodible elements in the local wind acceleration, drag partitioning, and protective coverage is formulated in a consistent manner.

To assess the performance of the new windblown dust treatment in CMAQ, an annual model simulation for 2011 was performed for the CONUS domain and results of the model were compared against observations from IMPROVE and AQS networks. Evaluation statistics indicate overall good performance of CMAQ with the new windblown dust treatment, with a domain-wide, annual MB for fine soil of 0.11 μg m−3. Seasonal investigations, however, showed moderate to high negative MBs for soil in the summertime, which is attributed to the inability of the model to simulate convective dust storms which are important contributors to dust generation in summer.

Considering the episodic nature of the windblown dust generation, we focused on the Southwest U.S. in the springtime (March–May) where the model showed good agreement with PM2.5 and PM10 surface measurements. Seven individual dust storms during spring 2011 were highlighted and performance of the model in capturing the timing and geographical pattern of these dust events was assessed. Comparison between model predictions and hourly PM2.5 and PM10 observations from AQS network sites along with event descriptions from the NWS storm database suggest that the new windblown dust emission treatment in CMAQ is able to capture the occurrence and temporal/geographical span of the dust outbreaks. The level of fine aerosol concentrations is simulated with relatively reasonable accuracy, while PM10 concentrations are systematically underpredicted, which requires further investigations.

As a summary, the present study indicates an overall good agreement between the newly developed windblown dust treatment in CMAQ and observations (with large improvement compared to the default dust parametrization in CMAQ v5.1; see supporting information Figures S4 and S5), and therefore, the model is believed to be able to serve as a useful tool to study the impact of dust emission in the atmosphere. Nevertheless, further developments including the effect of subgrid surface wind variability are still needed and are currently underway.