Improving the simulation of convective dust storms in regional-to-global models.

Convective dust storms have significant impacts on atmospheric conditions and air quality and are a major source of dust uplift in summertime. However, regional-to-global models generally do not accurately simulate these storms, a limitation that can be attributed to (1) using a single mean value for wind speed per grid box, i.e., not accounting for subgrid wind variability and (2) using convective parametrizations that poorly simulate cold pool outflows. This study aims to improve the simulation of convective dust storms by tackling these two issues. Specifically, we incorporate a probability distribution function for surface wind in each grid box to account for subgrid wind variability due to dry and moist convection. Furthermore, we use lightning assimilation to increase the accuracy of the convective parameterization and simulated cold pool outflows. This updated model framework is used to simulate a massive convective dust storm that hit Phoenix, AZ, on 6 July 2011. The results show that lightning assimilation provides a more realistic simulation of precipitation features, including timing and location, and the resulting cold pool outflows that generated the dust storm. When those results are combined with a dust model that accounts for subgrid wind variability, the prediction of dust uplift and concentrations are considerably improved compared to the default model results. This modeling framework could potentially improve the simulation of convective dust storms in global models, regional climate simulations, and retrospective air quality studies.

1. Introduction

Windblown dust has a major impact on weather, climate, air quality, and human health. Therefore, it is important to include its effects in regional-to-global models. However, this is difficult because windblown dust emissions are typically parameterized using a tight, but highly nonlinear relationship with the meteorological fields [Darmenova ]. In particular, dust emissions rate is shown to be a power function of surface wind speed with the power ranging from 3 to 5 [Shao, 2008]. Surface winds not only dictate the onset of dust uplift but also impact the dust emissions flux [Kok ]. Hence, representing the meteorological processes that generate strong surface winds is important to correctly predict windblown dust emissions.

Several processes can generate strong surface winds and windblown dust, including large-scale monsoon type flows, synoptic-scale systems, dry convection, and moist convective storms [Knippertz, 2014]. Dust emissions generated by cold pool outflows (also known as haboobs) from moist convective storms are of particular interest because coarse-resolution (>10 km) atmospheric models that use grid-averaged surface wind speed to parameterize dust emissions typically underestimate them [Largeron ]. This shortcoming can be attributed to at least two issues: first, most windblown dust parameterizations in regional-to-global models use a threshold wind speed (or threshold friction velocity) value for determining whether emissions occur or not. However, use of such a threshold might bias coarse-resolution model simulations that incorporate a single mean value of the wind speed (or friction velocity) for the entire grid box [Ridley ]. For instance, no dust will be generated if the mean modeled wind speed is slightly lower than the threshold value. However, in reality, if the mean wind is near the threshold value, it is likely that the wind exceeds the threshold value at some locations within the grid box, which would result in dust emissions.

The second issue is that regional-to-global models often use convective parameterizations that generally misrepresent moist convective storms and their resulting cold pool outflows. For example, through analysis of rain and lightning observations, Pope concluded that “haboobs (cold pool outflows from moist convection) are an important dust source in reality but are badly handled by the model’s convection scheme.” Heinold made simulations over West Africa with 40, 12, and 4 km grid spacing. The domains using a convective parameterization (40 and 12 km) were unable to realistically simulate the afternoon peak in convectively driven dust emissions. This was a significant shortage because approximately 40% of the total dust emissions in that region are generated by moist convection. Heinold concluded that the convective parameterization could not realistically simulate the timing and location of storms and their associated cold pool outflows, leading to a deficiency in convectively driven dust emissions. Many other studies [e.g., Marsham ; Garcia-Carreras ; Sodemann ; Pantillon ] have reached similar conclusions regarding convective parameterizations and dust emissions. Although it is recognized that the use of higher resolution alleviates some of these issues [e.g., Marsham ], many long-term and global studies will rely on convective parameterizations in the foreseeable future. Therefore, the inability of convective parameterizations to realistically generate cold pool outflows and their related strong surface winds presents a clear limitation for windblown dust modeling in regional-to-global models.

Mansell implemented lightning assimilation into the Kain-Fritsch (KF) [Kain, 2004] cumulus parameterization and found that “the lightning assimilation was successful in generating cold pools that were present in the surface observations.” The lightning assimilation technique used a simple approach: force deep convection where lightning is observed and suppress deep convection where it is not. Forecasts that used lightning assimilation during spin-up had more skill because of the more realistic production of cold pools. Heath adapted and modified the same lightning assimilation technique for retrospective applications and found substantial improvements in simulated rainfall during July 2012 and July 2013. Their runs using lightning assimilation also better simulated localized heavy rainfall events (see their Figure 4) that are typically missed by convective parameterizations [e.g., Stephens ].

Figure 4

Time series of spatial correlation of 1 h accumulated rainfall for WRF-CTRL (black) and WRF-LTGA (blue). The calculation was performed over the subdomain shown in Figure 2 The numbers next to the legend in the top right are the average values from 00:00 to 06:00 UTC 6 July 2011.

The results of these previous studies are encouraging and motivate the following question: could the combination of lightning assimilation and a dust model that accounts for subgrid wind variability improve convectively driven dust emissions in models that use a convective parameterization (e.g., regional-to-global models)? To answer this, here we develop a new dust model that includes the variability of the wind speed within a grid box and also implement lightning assimilation to improve the simulation of moist convection. This new modeling framework is implemented in the Weather Research and Forecasting (WRF) and Community Multiscale Air Quality (CMAQ) modeling systems and evaluated for a haboob dust storm event that occurred near Phoenix, AZ, on 5 July 2011. The next section describes the model updates that were designed to improve convective dust storm simulation. Section 3 briefly reviews the 5 July 2011 haboob dust storm and its resulting air quality impacts. The model configurations used to simulate the event are described in section 4. Section 5 examines the impacts of the modeling updates on the simulation of rainfall and surface wind, and the resulting effects on particulate matter (PM) concentrations (used as a proxy for dust) during the haboob event. A summary of the research and its key conclusions are presented in section 6.

2. Model Updates

2.1. Windblown Dust Emissions Parameterization

The windblown dust emissions parameterization developed by Foroutan is used in this study, which formulates the dust emission due to saltation bombardment (sandblasting). The physics of saltation include the horizontal movement of sand particles due to wind shear, the impact of these particles to the surface, and the vertical emission of dust particles [Kok ]. These processes were mathematically formulated in Foroutan as follows. First, the friction velocity is calculated using the wind speed and a surface roughness length relevant to the physics of dust generation [see Foroutan , sections 2.2 and 2.3]. Second, if the friction velocity exceeds a threshold value (which depends on the size of particles, soil moisture content, and the presence of nonerodible roughness elements), the saltation (horizontal movement) flux is obtained. Finally, the vertical flux of the dust is calculated based on a sandblasting efficiency formulation (vertical-to-horizontal dust flux ratio).

A distinct feature of the windblown dust model used here [Foroutan ] is that it incorporates a newly developed dynamic relation for the surface roughness length, which is important in correctly predicting both the friction velocity and its threshold value used in the dust emission model. Using this new relation, the effects of solid (such as pebbles) and vegetation nonerodible elements in local wind acceleration, drag partitioning, and protective coverage are formulated in a consistent manner. Additionally, the fraction of absorbed photosynthetically available radiation (fPAR) from the Moderate Resolution Imaging Spectroradiometer (MODIS) is used to represent more realistic time-varying vegetation coverage in this model. By performing an annual 2011 simulation over the continental United States (CONUS), Foroutan showed that the modeled soil concentrations agreed quite well with observations in the springtime. However, there was a considerable underestimation of soil concentration in the summertime when convective dust storms are most frequent, which was part of the motivation for the current study.

2.2. Subgrid Variability of the Surface Wind

To overcome the issue related to the use of the grid-averaged wind speed in regional-to-global models, we represent the wind speed with a probability distribution function and also account for the subgrid processes contributing to the variability of the wind speed within a grid box. The impact of including the subgrid variability of wind in large-scale atmospheric models has been investigated by several previous studies. Cakmur represented the components of the wind velocity by a bivariate normal distribution and derived a probability density function (PDF) for the wind speed. Marcella and Eltahir [2010] followed the methodology of Cakmur and considered a wind speed that follows a Gaussian distribution. More recently, Zhang derived a Weibull probability distribution to represent subgrid wind speed variability, which is the approach we adopt in this study, and is described below.

2.2.1. Weibull Distribution of the Wind Speed

As a first step in including the subgrid processes within each grid box, we make a distinction between the mean surface wind speed and the resolved mean surface wind speed where u and v are the two horizontal components of the wind vector and the overbar indicates the average value. The latter is usually available in an atmospheric model since ū and v̄ are being calculated. To obtain the total mean surface wind speed in a model, the effects of subgrid variability should be included, i.e., where u is the subgrid standard deviation of surface wind speed (also known as wind gustiness) and needs to be parameterized in the model.

Considering the mean surface wind speed U to be a random variable, a probability density function (PDF) can be used to represent its behavior. The Weibull distribution has been used in several previous wind studies [Justus ] and has specifically been shown to improve the dust emission prediction when included in general circulation models [Grini ; Ridley ]. The PDF of U, as a Weibull random variable, can be written as

In equation (4), k > 0 is the shape parameter and λ > 0 is the scale parameter of the distribution. They can be obtained based on the mean surface wind speed and its standard deviation [Justus ] using where Γ is the Gamma function.

2.2.2. Parameterization of the Wind Gustiness

The wind gustiness u is a key parameter in formulating the subgrid variability of surface wind. It represents subgrid scale processes causing the inhomogeneity of wind speed within a grid box in the model. In this study, we consider two physical processes responsible for the subgrid scale wind gustiness: dry and moist convection.

2.2.2.1. Dry Convection

In an unstable convective boundary layer, buoyancy results in updrafts and turbulent mixing, which brings higher momentum down from aloft and generates surface wind gustiness. Schumann [1988] used a large eddy simulation (LES) of a convective boundary layer for zero mean wind velocity to show that free convection can result in the horizontal wind variability. The level of gustiness depends upon the intensity of thermals and plumes, i.e., the magnitude of the surface buoyancy flux. Therefore, the contribution of dry convection to subgrid wind gustiness u can be written as [Godfrey and Beljaars, 1991] where w* is the convective velocity scale, which depends on the gravity g, boundary layer height H, virtual potential temperature θ, and the surface buoyancy flux . Several values for the scaling factor β have been proposed in the literature, ranging from 0.6 to 1.25 [Redelsperger ]. Here we use the value of 0.65 obtained by the LES study of Mondon and Redelsperger [1998] and recommended by Redelsperger to be used in mesoscale and general circulation models.

2.2.2.2. Moist Convection

The impacts of moist convection on subgrid wind variability are especially important for our case study of the haboob dust storm that occurred during July 2011, which was triggered by thunderstorm cold pool out-flows (section 3). In this study, we employ the parameterization of Zeng for the wind gustiness due to moist convection (precipitation and cloudiness). Through the analysis of the Goddard cloud-ensemble (GCE) model output for the cases of the Atlantic tropical experiment (GATE) and the coupled ocean-atmosphere response experiments (COARE1 and COARE2), Zeng developed the following relation for the contribution of moist convection to subgrid wind gustiness: Here R is the precipitation rate (mm h−1) and f is considered to be the grid cloud fraction not exceeding 0.5 (f=min [f, 1−f]). This choice of f is because the effect of clouds on grid-averaged variables is already accounted for in regional-to-global models, and for example, a cloud with 100% coverage does not expect to affect wind gustiness [Zeng ]. It remains unclear whether equation (8), which gives very similar results to those of Redelsperger , is appropriate over all regions of the globe. Nevertheless, with the lack of a globally derived relation, we have chosen equation (8) as a first step in our model development.

The total wind gustiness due to both dry and moist convection can then be obtained by

2.2.3. Dust Emissions Calculation With Subgrid Variability of Wind

The windblown dust emission parameterization is further enhanced by including the effect of the wind speed variability. The procedure can be summarized as follows (also see supporting information Figure S1):

For each grid box in the model, calculate the total mean surface wind speed U using equation (3) with U being calculated by the model, and u being calculated by equations (7)–(9).

Calculate the shape k and scale λ parameters of the Weibull distribution using equations (5) and (6).

Determine a Weibull distribution for the surface wind speed U using equation (4).

Equally divide the central 99% of the Weibull distribution into 50 bins, with each bin i having a sample wind speed U and a weight factor ω associated with it. This means that we consider the lowest wind speed in the grid box being less than 99% of the winds and the highest wind speed in the grid box being more than 99% of the winds. It should be noted that we performed a series of sensitivity analyses using 20 and 100 bins and obtained similar results.

For each bin i, calculate the friction velocity and consequently windblown dust emission using U as the wind speed [see Foroutan for details].

Calculate the total dust emission for the grid box by integrating over 50 bins considering the weight factor ω for each bin (resulting in the desired Weibull PDF).

2.3. Lightning Assimilation

Although the windblown dust model described above accounts for the effects of moist convection on subgrid wind variability (equation (8)), the parameterization is still dependent on the accuracy of the moist convection within the meteorological model itself. As previously stated in the introduction, coarse-resolution models that use convective parameterizations typically cannot accurately simulate many aspects of moist convection (e.g., timing, location, and distribution of rainfall intensity), which all impact resulting cold pool outflows. To mitigate this issue, we use the lightning assimilation technique of Heath , which builds off of the methods of Mansell , Lagouvardos , and Giannaros . The assimilation technique uses a straightforward approach of activating KF deep convection where lightning is observed and suppressing deep convection where lightning is absent. Particularly, Heath showed that triggering deep convection where lightning was observed and only allowing shallow convection where lightning was not observed was the optimal assimilation technique for retrospective applications. Below, we briefly describe the assimilation technique.

The KF scheme is called at 10 min intervals during the WRF simulation. During each call to KF, the code searches for lightning within −10 and +20 min of the current time at the current grid point. If lightning is present, the scheme first goes through its normal updraft calculations using the layer with the highest moist static energy as its updraft source layer (USL). If the resulting cloud does not meet the criteria for deep convection, water vapor and temperature perturbations are added to the USL using 0.1 g kg−1 and 0.1 K increments, respectively, until deep convection is forced. If deep convection cannot be achieved after 1 g kg−1 and 1 K have been added to the USL, no further action is taken and KF is not activated. Note that the method used by Heath added only moisture to the USL, but here we use a combination of moisture and temperature perturbations to increase the likelihood of activating deep convection at lightning grid points. Finally, if lightning is not observed at a grid point, only KF shallow convection is allowed to occur.

Cloud-to-ground lightning observations from the National Lightning Detection Network (NLDN) [Orville, 2008] were used for the assimilation described above. NLDN has a detection efficiency of 90%–95% and a location accuracy of approximately 500 m.

3. Test Case Overview

A major haboob dust storm occurred on 6 July 2011 (late afternoon/evening of 5 July local time) near Phoenix, AZ. The haboob was generated by late afternoon severe thunderstorms occurring near Tucson, AZ. Cold pool outflows associated with this region of large storms moved northwest toward Phoenix, bringing with them a wall of dust extending approximately 160 km wide and 1.5–1.8 km high (Figure 1). The dust storm brought near-zero visibility and wind gusts greater than 22 m s−1, making it important from both air quality and meteorological perspectives [Raman ; Vukovic ]. In fact, part of the motivation for choosing this case is that standard CMAQ simulations could not accurately simulate the haboob [Foroutan ], leading to large errors in PM predictions. The National Weather Service (NWS) reports that approximately 1–3 dust storms occur near Phoenix every year and that more than 100 were reported across all of Arizona between 2001 and 2011. Thus, it is important that we adequately simulate these storms and accurately assess their impacts on air quality and regional climate. For a complete description of the haboob studied here, please refer to the NWS report, which is available online at http://www.wrh.noaa.gov/psr/pns/2011/July/DustStorm.php.

Figure 1

Phoenix, AZ, dust storm on 5 July 2011. Photo by Daniel Bryant.

4. Data and Methods

4.1. WRF and CMAQ Model Configurations

The WRF model version 3.8 [Skamarock and Klemp, 2008] was used for this study. The simulations used 35 vertical layers and 12 km grid spacing covering the entire CONUS (472 × 312 cells). WRF runs were performed for all of 2011 as part of a separate U.S. Environmental Protection Agency (EPA) project [e.g., see Appel ]. Here we focus on the WRF output in a subdomain (31°N–35°N, 110°W–114°W) including southwest Arizona, and from July 2011, when the haboob event occurred near Phoenix, AZ (see Figure 2).

Figure 2

The modeled CONUS domain and the subdomain investigated in this study. The black dot shows the location of Phoenix, AZ, and background is the percentage of CMAQ grid cell covered by erodible land indicating the potential for windblown dust generation.

Two WRF runs were performed: a run without lightning assimilation (WRF-CTRL) and a run with lightning assimilation (WRF-LTGA). Initial and boundary conditions came from the North American Mesoscale (NAM) model 12 km analyses at 3 h intervals [Environmental Modeling Center ]. All simulations used a 50 hPa model top. In WRF-CTRL, the KF scheme was called at every model time step (60 s), whereas it was called every 10 min in WRF-LTGA (section 2.3). Cloud microphysics were parameterized using the Morrison two-moment scheme [Morrison ]. Short and long wave radiations were calculated using the RRTMG schemes [Iacono ]. Radiation interacted with the KF subgrid convection using the formulation in Alapaty and Herwehe . The Asymmetric Convection Model 2 (ACM2) [Pleim, 2007] boundary layer scheme was used with the Pleim-Xiu surface layer [Pleim, 2006] and Pleim-Xiu land surface model (PX LSM) [Pleim and Xiu, 1995; Xiu and Pleim, 2001; Pleim and Xiu, 2003]. Model output was saved every hour to drive CMAQ.

The two WRF simulations used indirect soil nudging in the PX LSM to adjust the soil moisture and temperature, which has been shown to reduce biases in 2 m temperature and 2 m mixing ratio [Pleim and Xiu, 2003; Pleim and Gilliam, 2009]. Surface observations were blended into the NAM analyses using WRF’s OBSGRID tool for this indirect nudging. Both simulations also used four-dimensional data assimilation (FDDA) [Stauffer and Seaman, 1994; Liu ], i.e., grid nudging, using the NAM analyses. Temperature, humidity, and horizontal winds were nudged toward the NAM analyses in grid cells above the boundary layer. The nudging coefficient for temperature and winds was 1×10−4s−1 and the coefficient for moisture was 1×10−5s−1.

The meteorological fields from these two WRF runs are processed by running version 4.2 of the Meteorology-Chemistry Interface Processor (MCIP) [Otte and Pleim, 2010] to provide input for CMAQ simulations. Here we use CMAQ model version 5.1 [Appel ] with a modified windblown dust scheme (section 2). Three sets of CMAQ runs are performed: (1) a run based on the meteorological fields of the WRF-CTRL run with the standard CMAQ dust model [Foroutan ] (CTRL), (2) a run based on the meteorological fields of the WRF-CTRL run with a dust model described in section 2 that incorporates the subgrid wind variability (SGWV), and (3) a run based on the meteorological fields of the WRF-LTGA run and the dust model that incorporates the subgrid wind variability (SGWV + LTGA). All the CMAQ simulations use the same 12 km horizontal resolution (with a domain slightly smaller than that of WRF, i.e., 459 × 299 cells) and 35 vertical layers extending up to 50 hPa as in WRF. Additional details of CMAQ run configurations can be found in Foroutan .

4.2. Observational Data

Precipitation observations come from the National Centers for Environmental Prediction (NCEP) stage-IV data set [Lin and Mitchell, 2005]. These data are a combination of radar-derived and rain gauge observations that have undergone human quality control at the NWS River Forecast Centers. Here we use the 1 h rainfall accumulation data that are provided on a 4 km grid, which we interpolated to our 12 km WRF grid for comparisons. The interpolation was done using the Gridded Analysis and Display System (GrADS) (documentation available at http://cola.gmu.edu/grads/gadoc/gadoc.php).

Particulate matter (PM) concentration observations are obtained from EPA’s Air Quality System (AQS; https://www.epa.gov/aqs) network. Specifically, we focus on concentrations of PM with diameter less than 10 µm (PM10) as a proxy for dust, which are measured hourly.

5. Results and Discussion

We first evaluate the ability of the WRF simulations to capture the convective event that triggered the haboob on 6 July 2011. Figure 3 shows snapshots of the 1 h accumulated rainfall from 02:00 to 06:00 UTC for the stage-IV observations (left column), WRF-CTRL (middle column), and WRF-LTGA (right column). For reference, the black dot in Figure 3 denotes Phoenix, AZ. At 02:00 UTC (top row), observations show an east-west oriented band of convection south of Phoenix. The WRF-CTRL simulation generates a more concentrated, circular area of rainfall that is displaced to the southeast when compared to observations. WRF-LTGA, on the other hand, generates an east-west band of rainfall that agrees better with observations. However, both models seem to overestimate the area of intense rainfall (>20 mm/h) that can affect dust emissions through changing cold pool outflow as well as soil moisture. An hour later at 03:00 UTC (second row), the spatial distribution of WRF-LTGA rainfall again matches better with observations than does WRF-CTRL. Over the next three hours (bottom three rows), the convective system moves toward the northwest and the resulting precipitation is more widespread over the western portion of the domain. The WRF-CTRL simulation shows the general movement toward the northwest but does not adequately capture the rainfall spatial distribution and instead continues to generate a more isolated, circular area of convection. WRF-LTGA shows better agreement with observations, generating more widespread precipitation over the western portion of the domain. With 12 km grid spacing, we cannot expect the model to perfectly simulate this transient convective event; however, it is encouraging that WRF-LTGA matches the spatial distribution of precipitation fairly well from 02:00 to 06:00 UTC.

Figure 3

Snapshots of 1 h accumulated rainfall (mm) from 02:00 to 06:00 UTC 6 July 2011 (19:00–23:00 local time on 5 July 2011) from the (left column) stage IV observations, (middle column) WRF-CTRL, and (right column) WRF-LTGA

To further evaluate the ability of the simulations to generate convection in the right location, which will subsequently be important for realistic dust generation, Figure 4 shows a time series of the simulated rainfall spatial correlation with stage-IV observations for WRF-CTRL (black) and WRF-LTGA (blue) from 00:00 to 06:00 UTC. As was qualitatively seen in the spatial snapshots of Figure 3, the WRF-CTRL run produces poor spatial correlations, with the maximum values only approaching 0.2. Conversely, WRF-LTGA produces higher spatial correlations that range between 0.4 and 0.6. Averaged over the 00:00–06:00 UTC period, the spatial correlation is considerably improved when lightning assimilation is used, increasing from 0.09 in WRF-CTRL to 0.46 in WRF-LTGA (numbers next to the legend in Figure 4).

Ultimately, cold pool outflows and low-level wind changes associated with the deep convection analyzed above caused the haboob dust event on 6 July 2011 (section 3). To further examine the observed and simulated cold pools, Figure 5 shows the correlation coefficient from the Phoenix, AZ, WSR-88D radar (Figures 5a and 5d) and the 2 m potential temperature anomaly and 10 m wind vectors for WRF-CTRL (Figures 5b and 5e) and WRF-LTGA (Figures 5c and 5f) at 02:00 and 03:00 UTC. The radar correlation coefficient (Figures 5a and 5d) shows the similarity between the horizontally and vertically oriented radar pulses in a pulse volume. Hydrometeors typically have a correlation coefficient near one whereas values lower than one indicate nonmeteorological features, such as insects, chaff, or for our case, dust. Thus, correlation coefficient can be used to help identify the leading edge of the haboob [e.g., Dempsey, 2014] and subsequently the observed location of the triggering cold pool outflow. The black dashed line in Figure 5 is the approximate location of the leading edge of the haboob based on the observed correlation coefficient (see also http://www.wrh.noaa.gov/psr/pns/2011/July/DustStorm.php). At 02:00 UTC (top row), the cold pool outflow in the WRF-CTRL simulation (Figure 5b) is farther southeast than is observed. The WRF-LTGA run (Figure 5c), with its more accurately placed convection (Figures 3 and 4), shows the leading edge of the cold pool outflow to be very similar to observations, and much more realistic than the WRF-CTRL run. At 03:00 UTC (bottom row), the leading edge of the haboob has moved past Phoenix (Figure 5d). The WRF-CTRL run cold pool outflow is again too far south/southeast. WRF-LTGA’s placement is again more realistic, but the progression of its cold pool outflow is still slower than was observed. Overall, as will be shown below, the more accurate location of WRF-LTGA’s cold pool outflow contributes to a more realistic simulation of the haboob dust event when subgrid wind variability is accounted for.

Figure 5

(left column) Correlation coefficient from the Phoenix, AZ, WSR-88D radar at (a) 02:00 UTC and (d) 03:00 UTC. (middle and right columns) Simulated 2 m potential temperature anomaly (K; shaded) and 10 m wind vectors (m s−1) from (b, e) WRF-CTRL and (c, f) WRF-LTGA. The radar correlation coefficient shows the similarity between the horizontally and vertically oriented radar pulses in a pulse volume and can be used to help identify the leading edge of the haboob. The dashed black line on each plot shows the approximate location of the gust front (based on the observed correlation coefficient) that generated the haboob dust storm.

Figure 6 shows hourly PM10 concentrations obtained from CMAQ simulations between 02:00 and 06:00 UTC, overlaid with the observations from 11 AQS sites located in the region (mainly in the Maricopa County and concentrated around the Phoenix metropolitan area). Observations from AQS sites show the episodic and rapid nature of this severe dust storm—PM10 concentrations intensely change in the Phoenix area within a few hours. The highest concentrations are observed at 03:00 UTC (20:00 local time, 5 July), when PM10 concentration reaches several thousands of µg m−3, while three hours later at 06:00 UTC, concentrations are at least one order of magnitude lower. Simulation results in Figure 6 are shown for three sets of runs: CTRL (left column), SGWV (middle column), and SGWV + LTGA (right column). It is clearly seen that CTRL, which uses default windblown dust model [Foroutan ], is unable to capture any signal related to the haboob dust storm. The downdraft governed by the convective parameterization in this relative coarse-grid simulation is not strong enough to provide the wind speed required for triggering dust generation. This, in general, results in underprediction of dust concentration in summertime as reported by Foroutan . When subgrid wind variability is introduced within the model, the friction velocity (calculated based on the surface wind) is allowed to exceed the threshold value at some locations within the grid box resulting in dust generation, as seen for SGWV (middle) and SGWV + LTGA (right). The intensity and spatial distribution of dust emission is closely related to the wind gustiness u, which in turn is determined by the strength of the simulated convection (cloudiness and precipitation). For the SGWV run, which uses the WRF CTRL run without lightning assimilation, dust generation is mainly isolated in a circular area that follows the pattern of the simulated convection (see Figure 3, middle column). Introducing lightning assimilation results in more scattered areas of dust generation which agree (at least qualitatively) better with AQS observations. Specifically, SGWV + LTGA is able to capture the dust front that passed through the Phoenix metropolitan area, while SGWV is not able to properly locate the dust storm. This is especially important from a practical point of view when the model is being used as a decision making tool.

Figure 6

Simulated hourly PM10 surface concentrations (µg m−3) between 02:00 and 06:00 UTC on 6 July 2011 (19:00–23:00 local time on 5 July 2011) from three runs (CTRL, SGWV, and SGWV + LTGA), overlaid with the observations of 11 AQS sites in the region.

The improvements in SGWV + LTGA are due to the change in both the resolved wind speed U and the wind gustiness u, as shown in Figure 7. As discussed above, the case with lightning assimilation better captures the location of the cold pool, which manifested itself as more realistic predictions of U. This improvement alone, however, is not enough to capture the dust outbreak in present simulations with a relatively coarse grid. In fact, our tests (see supporting information Figure S2) reveal that a run with only LTGA (and no SGWV) is not able to generate any dust. The main factor in enhancing the friction velocity above its threshold value somewhere within a grid box is SGVW, while LTGA has a major role in locating the cold pool (and consequently the dust uplift) in a correct place. When the lightning assimilation is combined with subgrid wind variability, an improved u is realized through better modeling of the precipitation, which allows the friction velocity distribution to pass its threshold value and generate dust in locations that agree better with observations. This combined effect can be clearly seen in Figure 7 at 02:00 and 03:00 UTC, when dust concentrations peak in the Phoenix area.

Figure 7

Simulated hourly values of resolved mean wind speed U and wind gustiness u (m s−1) between 02:00 and 06:00 UTC on 6 July 2011 (19:00–23:00 local time on 5 July 2011) from two runs (SGWV and SGWV + LTGA).

It should be noted, however, that we do not expect to capture all of the details with the 12 km grid spacing used in these simulations. For instance, both SGWV and SGWV + LTGA runs underestimate the level of PM10 concentrations, which may be attributed to predicting too much rainfall right before the dust storm (see Figure 3), resulting in increased soil moisture and threshold friction velocity, and consequently suppressed dust generation. Nevertheless, the considerable improvements obtained through the implementation of subgrid wind variability and lightning assimilation is encouraging for retrospective applications and potentially for operational modeling efforts that use relatively coarse resolution.

6. Summary and Conclusion

Windblown dust affects weather, climate, air quality, and human health, making its impacts an important component of regional-to-global atmospheric models. In particular, moist convective storms that generate windblown dust, also known as haboobs, are difficult to simulate because regional-to-global models (1) typically use convective parameterizations and (2) do not account for subgrid wind variability in their windblown dust parameterizations. In this study, we set out to improve the simulation of haboobs in regional-to-global models by using lightning assimilation to increase the accuracy of the convective parameterization and by developing a dust parameterization that accounts for subgrid wind variability (through a Weibull probability distribution of surface wind) and is linked to output from the convective scheme.

Lightning assimilation was applied in WRF and a new dust model was developed for CMAQ to simulate a haboob that occurred near Phoenix, AZ, on 5 July 2011. The use of lightning assimilation substantially improved the location of simulated precipitation when compared to stage-IV observations (Figures 3 and 4). Moreover, when compared to radar observations, the improved rainfall led to a more realistic simulation of the cold pool outflows (and resulting strong surface winds) that generated the haboob (Figure 5). Major improvements were obtained in the simulated PM10 concentrations associated with the haboob (Figure 6) when the improved meteorology was combined with the new dust model in CMAQ. Lastly, we found that the more accurate rainfall patterns generated by the lightning assimilation resulted in more widespread subgrid wind gustiness in the new dust model (Figure 7), explaining the more realistic simulation of the dust storm.

In conclusion, we found that the combination of lightning assimilation and a new dust parameterization with subgrid variability of surface wind improved the simulation of a convective dust storm in our relatively coarse-resolution model runs. These results are especially important for regions like West Africa, where convective dust storms can account for up to 40% of the total summertime dust emissions [Heinold ]. Therefore, when lightning data are available, this modeling framework could potentially improve dust simulations in global models, regional climate simulations, and retrospective air quality studies.