Presented here is a new dust modeling framework that uses a backward-Lagrangian particle dispersion model coupled with a dust emission model, both driven by meteorological data from the Weather Research and Forecasting (WRF) Model. This new modeling framework was tested for the spring of 2010 at multiple sites across northern Utah. Initial model results for March–April 2010 showed that the model was able to replicate the 27–28 April 2010 dust event; however, it was unable to reproduce a significant wind-blown dust event on 30 March 2010. During this event, the model significantly underestimated PM2.5 concentrations (4.7 vs 38.7 μg m−3) along the Wasatch Front. The backward-Lagrangian approach presented here allowed for the easy identification of dust source regions with misrepresented land cover and soil types, which required an update to WRF. In addition, changes were also applied to the dust emission model to better account for dust emitted from dry lake basins. These updates significantly improved dust model simulations, with the modeled PM2.5 comparing much more favorably to observations (average of 30.3 μg m−3). In addition, these updates also improved the timing of the frontal passage within WRF. The dust model was also applied in a forecasting setting, with the model able to replicate the magnitude of a large dust event, albeit with a 2-h lag. These results suggest that the dust modeling framework presented here has potential to replicate past dust events, identify source regions of dust, and be used for short-term forecasting applications.
Dust generated from wind erosion is the largest single contributor to particulate matter found in the atmosphere (Forster et al. 2007; Rind et al. 2009). Particulate matter from dust erosion can be transported over large distances from its original source, with significant ecological and societal impacts at downwind sites (Steenburgh et al. 2012; Goudie and Middleton 2001; Prospero and Lamb 2003). Examples of impacts on society include disruption to transportation due to reduced visibility as well as significant health impacts (Alharbi et al. 2013; Pauley et al. 1996; Ichinose et al. 2008; Kang et al. 2013).
Wind-blown dust events often lead to elevated levels of particulate matter (PM), which is regulated by the U.S. Environmental Protection Agency (EPA). Under the Clean Air Act, the EPA has established the National Ambient Air Quality Standards (NAAQS) for PM with aerodynamic diameters less than 2.5 and 10 μm (PM2.5 and PM10). The 24-h standards for PM2.5 and PM10 are 35 and 160 μg m−3, respectively. During wind-blown dust events, the NAAQS levels can be exceeded for both PM2.5 and PM10, with potentially significant impacts on human health (EPA 2011). Prolonged exposure to elevated PM2.5 levels can cause significant health issues in the form of cardiovascular and respiratory diseases as well as premature death, with the young and elderly being the most sensitive to high PM2.5 concentrations (Beard et al. 2012; Ichinose et al. 2008; Kang et al. 2013).
In addition to adverse health impacts, severe dust storms can also have a significant impact on visibility, resulting in disruptions to transportation and daily human activities (Liu et al. 2007). For example, a wind-blown dust event over the Arabian Peninsula on 10–11 March 2009 left hundreds of people hospitalized with respiratory-related issues, as well as causing multiple accidents as a result of reduced visibility (Alharbi et al. 2013). In 1991, a dust storm over the San Joaquin valley in California resulted in one of the worst multivehicle collisions in U.S. history, with 151 injuries and 17 fatalities (Pauley et al. 1996). More recently, a wind-blown dust event over Salt Lake City (SLC), Utah, during 14 April 2015 led to 25 injuries and 1 fatality, closure of major highways, and the cancellation of flights at the SLC International Airport.
Most dust emissions to the atmosphere occur over arid and semiarid environments such as northern China, Mongolia, Sahara, the Middle East, and the western United States (Liu et al. 2007; Zhang et al. 2003; Hahnenberger and Nicoll 2012; Prospero and Lamb 2003). Satellite observations and numerical simulations have shown that the western United States is a major source region for wind-blown dust (Ginoux et al. 2001; Woodward 2001; Zender et al. 2003). Across the western United States, exceedances of NAAQS levels can sometimes take place as a result of wind-blown dust events (Steenburgh et al. 2012; Hahnenberger and Nicoll 2012). Wind-blown dust is a particularly large problem for the SLC area, which anchors a large population center along the Wasatch Mountains with over 2 million people. This populated region known as the “Wasatch Front” is located adjacent to geographic features prone to soil erosion such as the Great Basin, the Great Salt Lake Desert (GSLD), Tule Dry Lake, Lake Sevier, Delta, and Milford Flats (Hahnenberger and Nicoll 2012; Steenburgh et al. 2012). It is for these reasons that we choose SLC as the area of interest for this study.
An example of dust being emitted across the eastern Great Basin can be seen in a visible satellite image for a dust event on 4 March 2009 (Fig. 1), which shows large dust plumes originating from the aforementioned dust source regions. Hahnenberger and Nicoll (2012) found that dust events across Utah occurred approximately 4.7 times per year, with one of these events exceeding the NAAQS for PM10. On average, most dust storms occurred during the late afternoon, with the highest frequency being reported during the month of April and a secondary maximum in September (Steenburgh et al. 2012). The majority of the dust events (68%) in Utah were accompanied by an approaching midlevel trough while the remaining 32% resulted from airmass convection (Hahnenberger and Nicoll 2012).
Exposure to adverse air quality is expected to continue into the foreseeable future for SLC as a result of rapid population growth (Utah Foundation 2014) combined with an increase in dust production associated with climate change (Munson et al. 2011) and intensification of anthropogenic activities, for example, agriculture, energy exploration/development, and recreation (Field et al. 2010; Neff et al. 2008). In addition, the Great Salt Lake (GSL), which is adjacent to SLC, has been steadily declining in size since the late 1800s due to increased water consumption along the Wasatch Front, exposing more erodible land surfaces that may contribute to dust emissions (Wurtsbaugh et al. 2016). Thus, a pressing need exists for a modeling framework that can 1) identify regions of significant dust emissions and their potential impacts on population centers, 2) quantify model uncertainty as a result of transport errors, and 3) determine the predictability of such events for areas susceptible to wind-blown dust events. Such a framework would be useful for both air quality and land managers for identifying sources of poor air quality as a result of wind-blown dust, determining where dust mitigation techniques can be deployed, and for forecasting air quality and visibility for major population centers.
Numerous studies have investigated wind-blown dust events from the Eulerian modeling approach, which simulates dust concentrations in grid cells (Fu et al. 2014; Albani et al. 2014; Haustein et al. 2012; Liu et al. 2007; Pérez et al. 2011; Han et al. 2004; L. T. Wang et al. 2012; Zender et al. 2003; Pantillon et al. 2015; Huang et al. 2015; Dong et al. 2016). Others have adopted a Lagrangian modeling approach, in which dust is transported in air parcels following the wind field (Draxler et al. 2001; Draxler et al. 2010; Ashrafi et al. 2014; Escudero et al. 2006). Some advantages of the Lagrangian approach over the Eulerian modeling perspective include 1) physical realism (Lin et al. 2013), 2) numerical stability (Wohltmann and Rex 2009), 3) lack of numerical diffusion (Shin and Reich 2009), 4) mass conservation (Lin et al. 2013), and 5) the ability to resolve subgrid variability (Lin et al. 2013). However, there are several disadvantages of Lagrangian models, including computational cost (Lin et al. 2013), grid irregularity (Stein et al. 2000), and inconsistencies with Eulerian driving meteorological fields as a result of interpolation (Alam and Lin 2008).
A subset of Lagrangian models, the Lagrangian particle dispersion models (LPDMs), simulate an ensemble of trajectories that account for turbulence via stochastic processes. Some examples of widely used LPDMs include the Flexible Particle Dispersion Model (FLEXPART; Stohl et al. 2005), HYSPLIT (Stein et al. 2015; Draxler and Hess 1998), and the Stochastic Time-Inverted Lagrangian Transport (STILT) model (Lin et al. 2003). In addition to accounting for turbulence, LPDMs can also run backward in time. Using this methodology, time-reversed LPDMs such as STILT, FLEXPART, and HYSPLIT are able to simultaneously determine upwind source regions along with their contributions toward concentrations simulated at the point of interest (i.e., the receptor) while accounting for turbulence often seen in the planetary boundary layer (Lin et al. 2013). As a result, backward LPDMs can be a computationally cheaper alternative when compared with forward-Lagrangian and -Eulerian approaches. However, this only holds true for cases when the number of potential source regions exceeds the number of receptors (Lin et al. 2003). Important advancements have also been made to the latest generation of backward LPDMs such as STILT, which now can quantitatively assess the impacts of modeled wind errors on the simulated concentration of trace gases and other species (Lin and Gerbig 2005). These advancements make the receptor-oriented modeling framework, based on backward LPDMs, a natural tool for identifying the source regions of tracers and quantifying their impacts on downwind locations.
Here, we present a dust modeling framework that leverages the benefits of the receptor-oriented approach using backward-LPDM simulations coupled to a dust emission scheme, driven with high-resolution wind fields from a mesoscale model. This framework will quantify the impact of wind-blown dust events on PM2.5 concentrations while simultaneously determining the source regions of wind-blown dust. In addition, wind error statistics will be incorporated into the dust simulations via stochastic parcel motions, which can be used to determine the uncertainty in modeled dust as a result of transport errors.
The PM2.5 measurements used for model evaluation were obtained from the Utah Department of Environmental Quality’s (UDEQ) air quality monitoring sites at Salt Lake City, Tooele, Provo, Ogden, Brigham City, and Logan (locations indicated in Fig. 1). The Tapered Element Oscillating Microbalance Filter Dynamic Measuring System was used to obtain hourly PM2.5 concentrations and is regularly maintained in accordance with federal regulations and instrument specifications (UDEQ 2012). Meteorological data used to validate the dust model framework discussed in the following section were obtained through the MesoWest cooperative networks (Horel et al. 2002). A list of all the measurement sites (both meteorological and air quality measurements) used in this study and their geographic locations can be found in appendix A.
b. Dust modeling framework
To simulate wind-blown dust events, we coupled a time-reversed LPDM (STILT) (Lin et al. 2003) to a dust emission scheme (FENGSHA) (Dong et al. 2016; Fu et al. 2014; Huang et al. 2015) driven with the Weather Research and Forecasting (WRF) Model, Advanced Research version 3.4.1 (Skamarock et al. 2008). WRF provides three-dimensional meteorological fields that drive STILT trajectories while also providing soil moisture, surface wind speeds, and land cover data to FENGSHA, which calculates dust emissions (Fig. 2). Integrating STILT with WRF enables air arriving at the receptor site to be traced to source regions, thereby quantifying contributions from upwind dust sources (Nehrkorn et al. 2010). This then allows for the calculation of dust concentrations at the receptor site.
1) Atmospheric transport
WRF is a Eulerian nonhydrostatic atmospheric mesoscale model widely used in the atmospheric science community to simulate various meteorological phenomena (Skamarock et al. 2008). Our WRF setup consisted of two domains gridded at 12- and 4-km spacing, with two-way nesting (Fig. 3). Boundary conditions were derived from the North American Regional Reanalysis (NARR) (Mesinger et al. 2006) with simulations being reinitialized every 7 days with 6 h of spinup time. We selected identical WRF parameterizations to those tested in a previous study that was also centered on Utah (Lin et al. 2017).
The STILT-based receptor-oriented modeling framework has hitherto been widely used to interpret CO2, CO, and other trace gases, while also being applied to chemically active species (Lin et al. 2017; Mallia et al. 2015; Wen et al. 2012). STILT was built upon the widely used HYSPLIT Lagrangian model developed at NOAA’s Air Resources Laboratory (Stein et al. 2015; Draxler and Hess 1998). STILT is an LPDM that simulates atmospheric transport by following the evolution of an ensemble of stochastic particles (or trajectories) moving backward in time, released for a receptor location (Fig. 2). These trajectories are driven by mean winds provided by a meteorological analysis such as WRF, with subgrid-scale turbulence being represented as a stochastic process. Each STILT particle represents a “parcel of air” that is small enough where it has similar properties to the surrounding air so that it is unaffected by gravitational settling and buoyancy but large enough where the size of the parcel is much larger than the average distance between molecules (Lin et al. 2013). These air parcels are also able to follow a turbulent flow without becoming deformed.
As part of the receptor-oriented framework, STILT trajectories link the receptor with upwind source regions via the footprint f (xr, tr|xi, yj, tm), where xr and tr are the receptor location and time linked to an upwind source at (xi, yj) and prior time tm (Lin et al. 2003). The footprint is simply the measure of the upwind surface influences for a receptor, as determined by the STILT backward trajectories of Lagrangian particles. The footprint can be defined as
where mair is the molar mass of air, h is the height of the volume in which the surface fluxes are diluted over (surface influence volume), is the average air density of trajectories within the surface influence volume, Ntot is the total number of backward trajectories, and Δtp,i,j,k is the amount of time a particle p spends within the surface influence volume at location (xi, yj) and time tm (Lin et al. 2003; Wen et al. 2012; Kim et al. 2013; Lin et al. 2013). By combining the footprint with surface emissions, STILT can quantitatively link surface emissions with concentration changes at the receptor.
For this study, we chose an ensemble size of 1000 backward trajectories for each receptor simulation hour, as suggested by trajectory number sensitivity tests carried out in Mallia et al. (2015). STILT backward trajectories were transported backward in time for 3 days or to the outer extent of WRF domain 1 (Fig. 3).
(i) Dry and wet deposition
To account for dry deposition of PM2.5, we applied a size-segregated and land-cover- dependent scheme based on the work of Zhang et al. (2001):
where Vd is the dry deposition velocity, Vg is the gravitational settling velocity, and Ra and Rs are the aerodynamic and surface vegetation resistances, respectively. The dry deposition calculation was carried out forward in time, along the STILT backward trajectories, every 2 min. At each time step, the dry deposition scheme was used to estimate a mass survival rate along trajectories where PM2.5 was present and within the planetary boundary layer. For further details behind this scheme, see appendix B. Wet deposition within STILT was calculated using the GEOS-Chem wet scavenging scheme for soluble tracers (Liu et al. 2001). Similar to dry deposition, wet deposition was also calculated as a survival rate along each backward trajectory, every 2 min.
(ii) Transport error
Recent developments have been made to STILT, which now allows the model to utilize errors in wind fields to quantify the uncertainty of modeled concentrations (Lin and Gerbig 2005). Errors in modeled winds are directly incorporated within STILT trajectories as an additional stochastic motion (ε):
where u is the total velocity vector, is the mean component of the wind, and u′ is the turbulent component. For the u- and υ-wind components, ε is assumed to closely resemble a Gaussian distribution with minimal biases (Lin and Gerbig 2005). The strength of ε is a function of the size of the wind errors and error length scales (l) in time and space, with l being calculated using an exponential variogram to fit wind errors:
where γ is the variance of the difference between two quantities separated by L. STILT trajectories incorporating these additional stochastic motions as a result of transport errors are further dispersed (Fig. 4), which reflects model uncertainty as a result of wind errors that are assumed to be random. The source contributions for each STILT trajectory are then calculated to determine the variance for simulations with and without transport errors. The difference in the variance of source contributions to the receptor’s dust concentrations before and after incorporating STILT transport errors is an estimate of transport errors on the simulated concentrations. For more details of this methodology, we refer to the reader to Lin and Gerbig (2005).
2) Dust emissions
To determine PM2.5 emissions during the wind-blown dust events, we used the FENGSHA dust emissions model (Fu et al. 2014; Huang et al. 2015; Dong et al. 2016). This dust emission model is currently used within the latest version of the Community Multiscale Air Quality model (CMAQ 5.0; www.cmaq-model.org) and has been deployed successfully in several dust studies (Fu et al. 2014; Huang et al. 2015; Dong et al. 2016). Using the FENGSHA model as implemented within CMAQ, the vertical flux of dust (F) (g m−2 s−1) can be calculated by
where i is the land cover type and j is the soil type. Here, K represents the ratio of the vertical flux to horizontal sediment (Marticorena and Bergametti 1995), A is the particle supply limitation set as 32 (unitless), ρ is the air density, and g is the gravitational constant (Fu et al. 2014; Dong et al. 2016). The SEP is the soil erodibility factor, which is dependent on the fractions of clay, silt, and sand in the soil (Fu et al. 2014; Dong et al. 2016). In addition, u* is the friction velocity and is the threshold friction velocity that determines the u* needed to transfer dust from soil surfaces to the atmosphere. The quantity u*t is dependent on the land cover (i) and soil (j) type, as well as the soil moisture. In addition, the FENGSHA model also accounts for snow cover and precipitation, which will flag the emissions as 0 if either of these conditions are present within the grid cell.
Explicit clay, silt, and sand fractions used to calculate SEP were obtained from the Soil Information for Environmental Modeling and Ecosystem Management (http://www.soilinfo.psu.edu/index.cgi?soil_data&index.html), instead of using fixed percentages based on the soil-type category. A PM2.5/PM10 ratio of 0.15 was chosen for dust emissions (K. Wang et al. 2012). The WRF fields, which contain information for land cover, ρ, u*, soil moisture, and snow cover, were used to drive the dust emission model.
a. Overview of March–April 2010 simulations
We tested the WRF-STILT dust modeling framework from 1 March to 30 April 2010 for multiple locations along the Wasatch Front (Fig. 1). Multiple frontal passages and dust events were observed in SLC during this period, which enabled testing of model performance under various wind and precipitation conditions. Error statistics for the entire simulation period can be seen in Table 1, which shows summary statistics of model versus observations comparisons at over 30 locations across Utah ( appendix A). It is worth noting that these locations were carefully selected so that observations were located either near STILT receptors, dust emission sources, or along potential transport pathways between dust emission sources and STILT receptors. Finally, wind error statistics calculated at the aforementioned locations were incorporated within STILT in an effort to determine the uncertainties of modeled PM2.5 as a result of modeled wind errors. Correlation length scales lx, lz, and lt used for the STILT transport error analysis were calculated as 50 km, 200 m, and 2.8 h, respectively, which are comparable to values computed by Lin and Gerbig (2005).
During this time, two major dust events associated with frontal passages with observed PM2.5 concentrations exceeding 200 μg m−3 on 30 March and on 27–28 April 2010 (Fig. 5). Both of these events resulted in violations in the NAAQS standard for PM2.5 of 35 μg m−3 (24-h average) at SLC with the 30 March 2010 and 27–28 April 2010 events averaging 39 and 40.9 μg m−3, respectively (Table 2).
To assess the performance of the WRF-STILT dust modeling framework, we compared modeled dust-derived PM2.5 concentrations with observed concentrations at SLC (Fig. 5), along with Provo, Tooele, Brigham City, Ogden, and Logan (Table 2). Note that observations indicated a general increase in PM2.5 levels between 15–18 March and 10–12 April 2010, which were not captured by the model. Low wind speeds and strong atmospheric stratification during these events suggest that these enhancements were most likely due to a strengthening surface inversion in the Salt Lake valley, which often causes an accumulation of particulate matter from anthropogenic sources near the valley floor, rather than advection of dust into the valley. Since the modeling framework presented here only considered dust emissions, the model did not capture these events.
The first modeled dust event was the 30 March 2010 event (Fig. 5). Initial results showed that the model drastically underestimated PM2.5 concentrations during the duration of the event at SLC (39 vs 3.1 μg m−3). Gross underestimations in modeled PM2.5 (−90%) were also observed at the five remaining study locations (Provo, Tooele, Brigham City, Ogden, and Logan), suggesting that the WRF-STILT modeling framework either underrepresented dust emissions and/or that significant transport errors existed within the model. Model performance was improved for the 27–28 April 2010 event with the model being able to replicate PM2.5 concentrations (Table 2). However, these results were somewhat inconsistent across the remaining sites (Table 2), with large underestimations in PM2.5 (−60%) being reported at Brigham City and Logan while a large overestimation was observed at Tooele (+100%) (Table 2).
Finally, the model indicated enhanced dust contributions on 5, 13, and 22 April, which were not seen in the observations. A quick analysis showed that the 22–23 April event in particular was associated with larger wind errors across much of Utah, likely resulting in the misplacement of upwind source regions. Errors for the other two events likely stemmed from overestimation in near-surface winds and an underestimation of precipitation prior to the frontal passage, which can limit dust production. For the next section of this study, we will attempt to identify the sources of errors seen during the 30 March (sections 3b and 3c) and 27–28 April 2010 (section 3e) events using the receptor-oriented framework of STILT. In addition, we will also test the modeling framework presented here in forecast mode for the 30 March event in section 3d while the forecast case for the 27–28 April 2010 case can be found in appendix C.
b. 30 March 2010
A significant dust event occurred on 30 March 2010 with a maximum measured hourly PM2.5 concentration of 250 μg m−3 and an average concentration of 39 μg m−3 throughout the duration of the event (Fig. 6a; Table 2). This event coincided with a strong frontal passage at 2100–2200 UTC associated with a trough moving in from the Pacific Northwest. Strong southerly winds gusting over 20 m s−1 were observed out ahead of the cold front with winds shifting to the northwest during the frontal passage at 2100–2200 UTC (Fig. 6a). Observed PM2.5 in SLC displayed two peaks, with elevated concentrations occurring before and after the frontal passage at 2100–2200 UTC.
Initial WRF-STILT model simulations for the 30 March 2010 event failed to replicate the spike seen in the observed PM2.5, with modeled concentrations remaining close to 0 μg m−3 throughout the duration of the event (Fig. 6a). Large underestimations in PM2.5 were also present at five other STILT receptor locations indicating that the underestimation of PM2.5 over SLC was not an isolated incident (Table 2). To identify the source of the model error, we first compared modeled wind directions with observed wind directions at the SLC airport to determine whether underestimations in PM2.5 were related to model transport errors. Overall, the modeled wind direction in Fig. 6a compared favorably to the observed wind directions between 1500 and 0000 UTC, with winds shifting from the south to the north-northwest by 0000 UTC as a result of the cold frontal passage. Granted, the modeled wind direction was a bit slower at veering toward the north when compared with observations by approximately 2 h. The modeled uncertainty in PM2.5 due to transports errors (Fig. 6a) was also minimal throughout the duration of the event, further indicating that the large underestimation in PM2.5 was likely not related to wind errors.
WRF-STILT footprints averaged between 2200 and 0000 UTC (Fig. 6b) were oriented to the west, which included major dust source regions such as the GSLD and western Nevada. Curiously, there were no dust contributions from the GSLD (Fig. 6b), despite the footprint being oriented over this region. This result contradicts findings in previous work that suggested that the GSLD is one of the major sources of dust emissions in Utah (Steenburgh et al. 2012; Hahnenberger and Nicoll 2012). Total dust emissions for the duration of the event across GSLD were close to 0 (Fig. 6c), despite the fact that modeled winds were in excess of 15–20 m s−1 for an extended period of time.
MODIS satellite overpasses early in the afternoon indicated that large dust plumes were originating from Lake Sevier and the Milford Flats regions, and blowing to the NNE toward SLC as a result of the strong prefrontal southerly winds (Fig. 7). This verifies that the modeled dust emissions across central Utah were likely real; however, dense cloud cover over northern Utah prevented the images from verifying dust emissions over the GSLD.
With no obvious transport errors in the modeled wind fields, and inconclusive results when comparing MODIS satellite images with modeled emissions, we decided to investigate the soil categories applied to the GSLD, which are obtained directly from WRF. Recent work by Massey et al. (2014) concluded that the WRF default USGS database is insufficient for providing realistic land cover information across much of western Utah due to the USGS’s lack of a dry lake basin (playa) soil-type classification for the GSLD. As a result, most of the GSLD was classified as silty–clay in the original WRF simulations, with those settings carrying over to the FENGSHA dust emission model (Fig. 8a). To rectify this issue, we updated WRF’s land cover and soil-type classification using the National Land Cover 2006 Database (NLCD 2006), which has a more realistic representation of land cover and soil types across western Utah (Massey et al. 2014) (Fig. 8b). Furthermore, we added playa to the FENGSHA dust emission model, which previously had no playa classification within it. Recent literature has highlighted the complexities of playa, whose physical characteristics are seasonally dependent. Playa often is capped by a hard crust during the summertime, which significantly increases u*t and results in minimal summer dust production, whereas the playa in wintertime and early spring is often associated with softer crusts resulting in a lower u*t and higher dust production (Gillette et al. 2001; Reynolds et al. 2007). As our analysis focuses on the spring, we opted to treat the western playa across Utah as a softer crust with a u*t set to 0.3 m s−1 for barren surfaces. This is significantly lower than the original u*t for this region, which was set as 0.6 m s−1 as a result of the previous silty–clay classification.
c. 30 March 2010 after NLCD update
After the updates to the WRF land cover and soil-type data, along with FENGSHA, the dust simulations for SLC greatly improved (Fig. 9a). In addition, the land cover and soil-type update also improved the timing of the frontal passage within WRF when comparing wind directions seen in Figs. 6a and 9a. During the peak of the event, modeled concentrations of PM2.5 spiked to ~270 μg m−3, which was comparable to the observed PM2.5 (251 μg m−3). As expected, significant dust contributions now originated from the GSLD after the cold frontal passage at 2100–2200 UTC with the GSLD seeing a large increase in total dust emissions (Figs. 9b,c). This is likely the result of decreased u*t across GLSD resulting from the updated playa soil-type classification.
It is worth noting that the prefrontal dust (1800–2200 UTC) within the model was still underestimated, even after the fix. We suspect that this issue is potentially related to underestimations in dust emission sources across the southern GSLD, and/or transport errors. Overall, the majority of dust contributions for the 30 March 2010 event originated from the GSLD (74%) while western and central Utah dust sources accounted for 15%. Sources outside of Utah accounted for the remaining 11% (Fig. 10).
Large improvements were also seen in modeled simulations of PM2.5 across the five other observation sites along the Wasatch Front (Table 2). After the land cover and soil-type update, average concentrations of PM2.5 increased significantly at all six sites, which ultimately compared more favorably to the observed PM2.5. In addition, correlation coefficients between the observations and modeled value of PM2.5 improved significantly (Table 2).
d. 30 March 2010 forecast simulation
While the dust modeling framework presented above performed adequately using simulations driven by reanalysis data, simulations that utilized forecast boundary conditions were necessary for determining the model’s capabilities in an operational forecast setting. To determine the predictability of the 30 March 2010 event, we cold started our WRF simulations 24 h before the onset of dust at 1200 UTC 30 March 2010 using the North American Mesoscale Forecast System (NAM) forecast (http://nomads.ncep.noaa.gov/). For example, if the forecast was 12 h into the simulation, WRF would be forced laterally by the NAM 12-h forecast.
Similar to the hindcast runs, the forecast simulations showed a large increase in PM2.5 concentrations during the event (Fig. 11a). While the duration and magnitude of the forecast was similar to that of the hindcast simulation, the timing of peak PM2.5 concentrations was delayed by approximately 2 h. This delay in enhanced PM2.5 can likely be attributed to the lag in wind shifts in the forecast simulation between 2000 and 0000 UTC. Comparisons between the WRF modeled wind fields for the hindcast (Fig. 11b) and forecast simulation (Fig. 11c) suggest that the cold front in the forecast simulation was significantly delayed. The slower progression of the cold front delayed the STILT footprints from veering toward the north over the GSLD, thus slowing the timing of the maximum PM2.5 (not shown). The failure of WRF to capture the timing of the cold frontal passage is reflected in the modeled PM2.5 uncertainties (±1σ), which are significant throughout the duration of the forecast, indicating low model confidence in the timing of the event (Fig. 11a).
As a result of the timing issue, correlation coefficients were <0.2; however, the magnitude and average of the event compared favorably to the observations (32 vs 39 μg m−3). This case study highlights the sensitivity of dust forecasts to wind direction and how it is imperative that the model forecasts correctly capture the timing of major wind shifts, which are often associated with frontal boundaries; otherwise, significant errors can result. In cases where the timing of the front is poorly resolved, model uncertainties can exceed 500%. However, the case study presented here also demonstrates the ability of STILT to indicate periods of significant uncertainty, which could be useful information for an operational weather forecaster.
e. 27–28 April 2010
The second wind-blown dust event occurred on 27–28 April 2010, and was also associated with a storm system moving in from the northwest. Winds out ahead of the storm system were out of the southwest (Fig. 12a), with wind speeds increasing to >15 m s−1 at 1100–2300 UTC 27 April. Similar trends were also observed at other sites across central and northern Utah. Winds remained elevated through 0800 UTC 28 April before subsiding during the early morning after the frontal passage at 0900 UTC across SLC. Observed PM2.5 in SLC started to increase after 2200 UTC 27 April, spiking to 237 μg m−3 by 0500 UTC the following day (Fig. 12a). Shortly after reaching peak value, the observed PM2.5 sharply dropped off to ~50 μg m−3, and remained elevated through 1200 UTC before returning to near background levels.
STILT-simulated PM2.5 showed reasonable agreement with the observed concentrations of PM2.5 at SLC with the magnitude and timing of the event well resolved by the model (Fig. 12a). In addition, the event-averaged PM2.5 compared favorably to the observations, with a correlation coefficient of R = 0.5. In contrast with the 30 March 2010 event, most of the modeled dust contributions originated from dust source regions across central Utah (Figs. 12b,c). This was likely the consequence of the predominantly southerly flow throughout the duration of the event (Fig. 12a). As a result, SLC saw limited changes in average modeled PM2.5 for the April event before (37.2 μg m−3, R = 0.5) and after (44.7 μg m−3, R = 0.5) the land cover fix, since most of the changes occurred over the GSLD and select dust hot spots across western Utah.
Prior to 0600 UTC, both the model and observations were in agreement, with PM2.5 averaging between 50 and 70 μg m−3. However, after 0600 UTC, modeled concentrations of PM2.5 started to diverge significantly from the observations, with overestimations exceeding +100%. During this time, the observed PM2.5 concentrations fell outside of the range of modeled uncertainty, indicating an issue that may not necessarily be limited to transport errors. Finally, by 1100–1200 UTC, both observed and modeled PM2.5 concentrations dropped significantly after the passage of the cold front (Fig. 12a).
According to the model (Figs. 12b,c), the majority of PM2.5 measured at SLC originated from the Milford Flats region (52%). Delta, which is the other central Utah dust source region, accounted for 6% of PM2.5 contributions. The western Utah dust source regions (Lake Sevier, Escalante Desert, Tule, GSLD), accounted for the remaining 35% of PM2.5 dust contributions when winds started to veer from the south-southwest to west between 0600 and 1200 UTC (Fig. 12a). It was during this period of time that the model exhibited the most significant errors in PM2.5, which could not be entirely explained by uncertainties as a result of model transport errors (Fig. 12a). This result suggests there may be errors associated with weather-dependent emission values and/or soil profiles over this region.
Model simulations for the other sites performed reasonably well, with most of the modeled event averages comparing well to the observations (Table 2). Correlation coefficients across most sites were reasonable, with the exception of Provo, which had a correlation of R = −0.1 (Table 2), due to a 2-h timing mismatch between the model and observations (not shown). While the correlation coefficients for Tooele were acceptable (R = 0.5), there was a considerable overestimation (+260%) of PM2.5 (Table 2; Fig. 13a). Furthermore, simulated values of PM2.5 prior to the land cover and soil-type update actually compared more favorably to the observations, with an overestimation of only +90%. This suggests that the land cover and soil-type update applied to the dust model may not necessarily be a perfect remedy for all emission sources.
The initial onset of the dust event at Tooele was well captured by the model. However, by 0300 UTC, modeled PM2.5 started to significantly diverge from observed PM2.5 concentrations (Fig. 13a). Throughout the remainder of the event, modeled concentrations remained extremely elevated, with concentrations > 100 μg m−3 until 1200 UTC the next day. Modeled uncertainties for PM2.5 suggest that some of these errors can be attributed to transport errors within WRF-STILT; however, it is possible that there are also issues with the dust emission model.
Similar to SLC, the majority of modeled PM2.5 contributions for Tooele originated from central and western Utah (Figs. 13b,c). However, a smaller percentage of dust originated from Milford Flats while a greater percentage originated from dust sources located farther to the west (Lake Sevier, Tule, Lake Sevier, Escalante Desert, and GSLD). Model errors for both SLC and Tooele were most pronounced when PM2.5 contributions originated from western Utah, with model errors being well correlated with PM2.5 originating from western Utah (Fig. 13a; R > 0.8). Wind errors at observational sites across western Utah were relatively small, as the model and observations were in good agreement throughout much of the event with wind speed differences rarely exceeding 2.5 m s−1 and wind directional errors <15°. This wind speed analysis suggests that the model errors are less likely the result of emission-dependent variables such as u* and more likely the result of errors associated with modeled soil characteristics such as the threshold friction velocity (u*t). However, it is plausible that transport errors are responsible for some of the model error since the shaded uncertainties in Fig. 13 only accounts for ±1σ.
Brigham City saw an improvement in simulated PM2.5 during the 27–28 April 2010 event (Table 2). Overall, PM2.5 concentrations compared reasonably to observations, albeit with a slight overestimation of +30% (Fig. 14a). However, the timing of the event, especially during the latter half, was offset by approximately 2 h. Unfortunately, the meteorological station near Brigham City was missing wind observations during the frontal passage, rendering it difficult to determine whether this error was wind related or caused by something else. Similar to Tooele, the majority of dust arriving at Brigham City came from western Utah (Figs. 14b,c); however, Brigham City did see a larger contribution from the GSLD (40%) relative to Tooele (20%). From this, we are left to hypothesize that the land cover and soil update may have improved dust emissions across the GSLD for the 27–28 April 2010 event (along with 30–31 April). However, the Tooele analysis suggests that the Tule and Lake Sevier dust emission hot spots are still problematic for the dust emission model, despite the land cover and soil-type update.
4. Summary and discussion
To the best of our knowledge, the modeling framework presented here is the first implementation of a backward LPDM for dust modeling applications. A major benefit of the backward-LPDM approach is that this method is a natural tool for identifying the source regions of tracers and quantifying their impacts on downwind locations. In addition, the backward-LPDM approach presented here has the unique ability of estimating the impacts of wind errors to determine uncertainties in modeled dust concentrations.
To demonstrate the capabilities of the LPDM approach presented here, the model was used to simulate dust events across Utah throughout the spring of 2010. During this time, two significant wind-blown dust events occurred on 30 March and 27–28 April 2010. Initial results showed that the model framework initially struggled with resolving a significant dust event on 30 March 2010, with considerable underestimations occurring at multiple locations along the Wasatch Front (−90%). By using the receptor-oriented approach with transport error capabilities, we were able to isolate model errors for the 30 March 2010 event to poorly characterized soil profiles across the GSLD, which led to a successive update to WRF’s land cover and soil-type categories. Additional updates were then made to the FENGSHA dust emission model to account for playa soil types, which have a lower u*t than the previous silty–clay classification over the GSLD.
Ultimately, the land cover and soil-type update resulted in significant increases in modeled PM2.5 concentrations across all sites along the Wasatch Front as a result of increased dust production across the GSLD. Furthermore, this update also resulted in the better timing of the cold frontal passage at 2200–2300 UTC, which further improved the dust simulations. We suspect that the better timing of the frontal passage was the result of improved interactions between the land surface and atmosphere, which is consistent with results from other studies across this region (Massey et al. 2014, 2017). From this example, we are able to demonstrate the backward LPDM’s ability in identifying upwind source regions, and how that knowledge can be used to make refinements to emission inventories based upon observed model errors. This case study also demonstrates the necessity of using up-to-date land cover and soil-type data, especially for dust modeling applications.
To assess the model’s veracity in the forecast setting, a 24-h forecast simulation was performed, which was driven by forecast boundary conditions for the 30 March 2010 dust event. The model was able to capture the duration and magnitude of the event. However, a lag in the timing of maximum PM2.5 concentrations existed within the forecast simulation, which can be attributed to wind direction errors as a result of a delayed frontal passage. Through this forecast case study, we are able to demonstrate the potential for backward LPDMs to be used as a forecasting tool for predicting dust events. Furthermore, STILT’s ability to quantify periods of significant uncertainty in a dust forecast could be a useful tool for the operational weather forecaster, as demonstrated in the 30 March 2010 case study.
With that said, there are some caveats that need to be considered when identifying the correct modeling approach to deploy for a particular forecast application. The backward-LPDM approach can be a computationally cheaper alternative when compared with forward-trajectory and Eulerian modeling methods. However, this is only true when the number of sources exceeds the number of receptor locations. For the 30 March and 27–28 April 2010 dust events, we identified over 15 individual dust sources in Utah, some of which exhibited significant spatial heterogeneity in emission profiles. For this particular application, the receptor-oriented approach would be computationally cheaper by a factor of 2.5 when compared with forward-oriented approaches since there were six receptors and approximately 15 source regions. However, for cases where many more receptor locations are needed—either for forecasts covering a larger geographic region and/or for a forecast requiring dust simulations to be resolved at much finer scales—the backward-Lagrangian modeling approach becomes a less effective tool relative to Eulerian and Lagrangian forward approaches. Finally, additional considerations are needed on how to correctly implement model transport errors into forecast simulations. For this study, we had the benefit of using observations to calculate the forecast transport errors. However, in a real operational forecast setting, these observations would be unavailable. Thus, some type of model validation study would likely be needed to estimate the average (climatological) transport error associated with major wind shifts such as during a frontal passage. Some of this is discussed in greater length by Lin et al. (2007), who used STILT for planning regional-scale trace gas field experiments.
While the dust simulations for PM2.5 performed adequately for the case studies presented here, further improvements can likely be made to dust modeling systems, as indicated by the large model errors observed at Tooele during the 27–28 April 2010 event. Despite using updated land cover and soil-type data resolved at a fine grid resolution, significant overestimations (+250%) were still reported for certain receptor locations. Using the backward-LPDM approach at multiple receptor locations, combined with a transport error analysis, we speculate that there may be some remaining issues with the dust emission model over specific hot spots over western Utah. For example, the land cover update appeared to significantly improve emissions over the GSLD, while potentially overestimating dust emissions for Tule Lake and Lake Sevier. As a result, we hypothesize that a better representation of soil types and/or u*t is likely needed over some regions, although it is plausible that some of this error was due to model transport errors. In its current formulation, u*t is determined by the soil type assigned to a particular grid cell. However, similar soil types may have different u*t as a result of surface disturbances such as fires, human activity, or previous wind events (Okin et al. 2011; Draxler et al. 2001, Gillette et al. 1980).
Dust sources across western Utah are further complicated by the fact that a large number of these sources consist of playa, which can have a physically complex emission profile. For example, u*t can also vary significantly for playa soil types, which can emit more dust when they are soft (Gillette et al. 2001). Here, we assumed that the crust was soft, which is often the case for the GSLD during the spring months when precipitation is maximized across this region. However, this may not be a good assumption for other seasons, particularly for the summer and early fall. In addition, different dry lake beds across the globe may experience different seasonality patterns in terms of the hardening and softening of playa crusts. Furthermore, dry lake beds can sometimes consist of standing water, which can significantly inhibit dust production. This requires the need to implement some type of high-resolution product such as the MODIS shortwave infrared band, which has been previously used for surface water mapping applications during relatively cloud-free days (Ticehart et al. 2014). Finally, another assumption made in this study was that the PM2.5/PM10 ratio was fixed for all emissions sources. However, this ratio can be affected by soil type, texture, and wind speed (Fu et al. 2014).
With the global population set to increase over the next century in addition to the impacts of climate change on land surfaces (Munson et al. 2011; Okin et al. 2011), dust events will likely become more frequent, requiring the need for advanced atmospheric transport models that can accurately simulate wind-blown dust. The GSL is set to see decreasing lake levels due to increased water consumption that will result in more exposed erodible lake bed, which will further deteriorate air quality along the Wasatch Front (Wurtsbaugh et al. 2016). Additionally, other water bodies across the globe, such as the Aral Sea, Dead Sea, and Lake Urmia (Smithsonian Institution 2014), are all experiencing rapid decreases in lake levels, which may lead to significant lake bed exposures and dust emissions. The modeling framework presented here can be used to identify dust emission hot spots where mitigation techniques can be deployed to reduce the severity and frequency of future dust events, as well as serving as a forecasting tool for predicting dust events for urban centers such as SLC and beyond.
This study was made possible by funding from the Utah Department of Environmental Quality, the Friends of Great Salt Lake, and the University of Utah’s Global Change and Sustainability Center. The authors also thank Joshua Benmergui for providing scripts that calculate depositional loss and Erik Crosman, Jeff Massey, and Maura Hahnenberger for their valuable input. The support and resources from the Center for High Performance Computing at the University of Utah are also gratefully acknowledged. The data used to produce the results in this study are available from DVM upon request. Finally, we thank the anonymous reviewers for their suggestions, which led to significant improvements to the text.
Measurement Sites Used in This Study
List of air quality and meteorological observation sites used in this study, in addition to their geographical location (Table A1).
Dry Deposition Calculation
The dry deposition velocity used within STILT is size segregated with a dependence on land cover based on the work of Zhang et al. (2001). Using this scheme, we estimate the mass survival rate M as
where t is the model time step (2 min), h is the PBL height, and z is the STILT trajectory height above ground level. The dry deposition velocity Vd can be calculated as
where Vg is the gravitational settling velocity, while Ra and Rs are the aerodynamic and surface vegetation resistances, respectively. We can define Vg as
where ρa is the density of the particle (1.3 μg m−3) and D is the mean diameter of the aerosol, which was assumed to be 0.25 μm for this study. The gravitational acceleration constant g = 9.8 m s−2; n is the dynamic viscosity of air and equal to 1.8 × 10−5 kg m−1 s−1; and Cc is the Cunningham slip correction, which can be written as
λmfp is the mean molecular free path,
where λstp = 6.53 × 10−8 m and ρstp = 1.275 kg m−3 are the mean free path of air molecules and the density at standard temperature and pressure, respectively. The Ra, which is used in the calculation for Vg, can be defined as
where u* is the friction velocity, k = 0.4, is the von Kármán constant, while zo is the roughness length of the underlying surface, and ψh is the Businger profile estimate of a dimensionless temperature gradient. The Rs is defined as
where εo = 3 is an empirical constant for surface resistance (Zhang et al. 2001), R1 is the correction factor that represents the fraction of particles that stick to the surface. In addition, EB is the efficiency of surface collection, impaction EIM, and interception EIN, which can be formulated as
where Sc is the Schmidt number, ξ ranges from ½ to ⅔ and is a function of the surface roughness, β = 2, α is dependent on the land cover type, and AR is the characteristic radius. The Stokes number St for smooth (Stsmo) and vegetated surfaces (Stveg) can be given as
where υ is the kinematic viscosity of air and equal to 1.51 × 10−5 m2 s−1. We can simply write R1 as
Land cover categorizations used for the dry deposition script were obtained from the Land Cover Type Climate Modeling Grid product (MCD12C1 type 1) provided by the USGS Land Processes Distributed Active Archive Center. Zhang et al. (2001) further translated this product to include seasonal patterns in land cover changes.
27–28 April 2010 Forecast Case Study
The 27–28 April 2010 event was also simulated in forecast mode, similar to the 30 March 2010 event. Our WRF run was initialized 24 h prior to the onset of enhanced PM2.5 at SLC at 1800 UTC 26 April 2010. Unlike the hindcast simulation, which realistically captured dust effects, the forecast significantly overpredicted PM2.5 concentrations (+200%) at SLC throughout the duration of the event (Fig. C1) while the correlation coefficient for this event was poor: R < 0.2. At 0100 UTC, the STILT footprint was oriented out of the south-southwest (Fig. C1b), which was fairly similar to the reanalysis footprint shown in (Fig. C1a). However, large contributions of dust were occurring approximately 100 km to the south-southwest of SLC near Delta, Utah, which was responsible for the large discrepancy in PM2.5 concentrations between the hindcast and forecast simulation at 2300–0200 UTC (~380 vs 60 μg m−3) (Fig. C1b). Delta is located next to the Little Sahara Sand Dunes, which is classified as “sand” within the FENGSHA dust emission model. Wind speeds across that area were 4 m s−1 higher in the forecast than in the hindcast simulation, which can likely explain the differences between the two PM2.5 simulations.
A similar overestimation was also observed during the second half of the event between 0800 and 1100 UTC. During this time, however, differences between the forecast and observed PM2.5 concentrations can be attributed not only to the overestimation in wind speed during the first half of the event but also to a delay in the veering of forecast wind directions, which resulted in the source region being erroneously placed farther to the south and east (Figs. C2c,d).