Application of a satellite-retrieved sheltering parameterization (v1.0) for dust event simulation with WRF-Chem v4.1

Roughness features (e.g., rocks, vegetation, furrows) that shelter or attenuate wind flow over the soil surface can considerably affect the magnitude and spatial distribution of sediment transport in active aeolian environments. Existing dust and sediment transport models often rely on vegetation attributes derived from static land use datasets or remotely sensed greenness indicators to incorporate sheltering effects on simulated particle mobilization. However, these overly simplistic approaches do not represent the three-dimensional nature or spatiotemporal changes of roughness element sheltering. They also ignore the sheltering contribution of non-vegetation roughness features and photosynthetically inactive (i.e., brown) vegetation common to dryland environments. Here, we explore the use of a novel albedo-based sheltering parameterization in a dust transport modeling application of the Weather Research and Forecasting model with Chemistry (WRF-Chem). The albedo method estimates sheltering effects on surface wind friction speeds and dust entrainment from the shadows cast by subgrid-scale roughness elements. For this study, we applied the albedo-derived drag partition to the Air Force Weather Agency (AFWA) dust emission module and conducted a sensitivity study on simulated PM₁₀ concentrations using the Georgia Institute of Technology–Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) model as implemented in WRF-Chem v4.1. Our analysis focused on a convective dust event case study from 3–4 July 2014 for the southwestern United States desert region discussed by other published works. Previous studies have found that WRF-Chem simulations grossly overestimated the dust transport associated with this event. Our results show that removing the default erodibility map and adding the drag parameterization to the AFWA dust module markedly improved the overall magnitude and spatial pattern of simulated dust conditions for this event. Simulated PM₁₀ values near the leading edge of the storm substantially decreased in magnitude (e.g., maximum PM₁₀ values were reduced from 17 151 to 8539 µg m⁻³), bringing the simulated results into alignment with the observed PM₁₀ measurements. Furthermore, the addition of the drag partition restricted the erroneous widespread dust emission of the original model configuration. We also show that similar model improvements can be achieved by replacing the wind friction speed parameter in the original dust emission module with globally scaled surface wind speeds, suggesting that a well-tuned constant could be used as a substitute for the albedo-based product for short-duration simulations in which surface roughness is not expected to change and for landscapes wherein roughness is constant over years to months. Though this alternative scaling method requires less processing, knowing how to best tune the model winds a priori could be a considerable challenge. Overall, our results demonstrate how dust transport simulation and forecasting with the AFWA dust module can be improved in vegetated drylands by calculating the dust emission flux with surface wind friction speed from a drag partition treatment.

1 Introduction

Surface roughness features such as rocks, vegetation, and soil ridges created by tillage attenuate wind flow over the soil surface, considerably affecting aeolian sediment transport and dust emission patterns (see reviews by Mayaud and Webb, 2017; Shao et al., 2015). Representing aerodynamic roughness effects and their spatiotemporal dynamics effectively in numerical dust models is therefore critical for estimating and forecasting the spatial patterns, timing, magnitude, and frequency of mineral dust emission accurately (Evans et al., 2016; Fu, 2019; Ito and Kok, 2017; Li et al., 2013; Tegen et al., 2002; Webb et al., 2014a, b). Several methods are available for characterizing vegetation effects in dust emission models (e.g., Evans et al., 2016; Ginoux et al., 2001; Ito and Kok, 2017; Kim et al., 2013; King et al., 2005; Koven and Fung, 2008; Marticorena and Bergametti, 1995; Okin, 2008; Raupach et al., 1993; Tegen et al., 2002). However, many of these techniques require in situ data that are challenging to obtain over large spatial footprints or incorporate parameterizations with dependencies on land use or vegetation datasets that are often out of date or fail to represent the complex three-dimensional heterogeneity of landscape roughness elements (e.g., Pierre et al., 2014; Raupach and Lu, 2004).

Tuning dust models to satellite- and ground-based observations of dust concentrations and dust optical depth has generally been a sufficient approach for making dust models useful for large-scale climate and air quality applications (e.g., Cakmur et al., 2006). These particular applications tend to focus on total atmospheric dust aerosol loading trends instead of individual dust events. If the goal is to capture the general magnitude and frequency of large-scale dust aerosol patterns, tuning practices are often a viable solution to simulated dust emission errors. However, the relatively poor representation of surface roughness effects on wind erosivity in current dust models has limited our capacity to forecast individual dust events accurately in desert regions with vegetation. Furthermore, poor roughness effect representation limits our ability to accurately simulate the influence of land management practices and desertification on spatial and temporal patterns of dust emission (Webb and Pierre, 2018).

Sediment mobilization schemes are often represented in terms of wind friction speed, u_*, a scalar parameter commonly used to describe processes related to wind shear stress (τ; note that $\mathit{\tau }=\mathit{\rho }(u_{*}{)}^{\mathrm{2}}$, and ρ is air density). Near the land surface, u_* represents the total wind shear stress (τ) acting on both the horizontal soil surface (τ_s) and roughness elements (τ_r) (i.e., $\mathit{\tau }={\mathit{\tau }}_{\mathrm{s}}+{\mathit{\tau }}_{\mathrm{r}}$; see Raupach, 1992; Raupach et al., 1993). This process is typically termed drag partitioning and is often expressed in terms of u_* rather than τ. Since τ is proportional to u_*, we can similarly divide u_* into soil surface (u_s*) and roughness (u_r*) components (i.e., $u_{*}=u_{\mathrm{s}*}+u_{\mathrm{r}*}$). The wind shear stress that reaches the immediate soil surface (i.e., τ_s) governs particle mobilization, so dust emission models driven by u_s* (or wind erosivity) instead of u_* may produce better outcomes (e.g., Darmenova et al., 2009; Okin, 2008; Webb et al., 2020). However, this assumption is worth testing for individual modeling applications since more physically sophisticated methods do not always lead to better simulation outcomes, especially when the more advanced modeling approaches are limited by uncertainties in their required input parameters (e.g., LeGrand et al., 2019).

Though roughness elements affecting u_s* could comprise any ground feature obstructing airflow over the land surface, traditional aeolian transport models generally equate roughness to vegetation cover, with varying degrees of sophistication in their approach. For example, some models restrict or reduce dust fluxes based on static prescribed vegetation or land use attributes (e.g., Ginoux et al., 2001; LeGrand et al., 2019; Woodward, 2001). Others implement dynamic masks or spatially varying dust flux scaling factors based on real-time or climatological datasets derived from greenness fraction, leaf area index, or normalized difference vegetation index satellite data (Asadov and Kerimov, 2019; Collins et al., 2011; Evans et al., 2016; Ito and Kok, 2017; Kim et al., 2013; Kok et al., 2014; Solomos et al., 2019; Spyrou et al., 2022; Tegen et al., 2002; Vukovic et al., 2014). While some of these techniques have proven useful in certain modeling applications, they depend on datasets that primarily highlight plant productivity phases (e.g., Yu et al., 2016), which may not translate well to soil surface sheltering (Okin, 2010).

Chappell and Webb (2016) proposed a method for parameterizing u_s* in aeolian transport models using remotely sensed surface albedo (ω). Their technique infers the drag partition from shadows (1−ω) cast by roughness elements based on the assumption that shadows can serve as a proxy for the sheltered surface area extent (Chappell et al., 2010). As described in Ziegler et al. (2020), a generalizable form of the Chappell and Webb (2016) approach equates to

$$\begin{array}{}\text{(1)}& u_{\mathrm{s}*}=u_{\mathrm{ns}*}{U}_{\mathrm{10}\mathrm{m}},\end{array}$$

where U_10 m is the wind speed 10 m a.g.l. (above ground level), and u_ns* is a normalized u_s* parameter (unitless) derived from albedo. By using an albedo-based calculation for each pixel area, the Chappell and Webb (2016) method provides an areal estimate of u_s* for both the sheltered and non-sheltered zones of the soil surface within each grid pixel.

Here, we explore the use of the Chappell and Webb (2016) albedo-based drag partition within the Weather Research and Forecasting model with Chemistry (WRF-Chem). WRF-Chem is a physically based Earth system model that simulates atmospheric motion on a non-hydrostatic, Eulerian grid in addition to the emission, transport, and mixing of gases and aerosols simultaneously with the meteorology (Fast et al., 2006; Grell et al., 2005; Peckham et al., 2017; Skamarock et al., 2019). The WRF-Chem framework is configurable and can be run with a variety of aerosol and atmospheric chemistry parameterizations, depending on a user's specific interests and computational resources. For this effort, we chose to use the AFWA dust emission module (LeGrand et al., 2019) with the Georgia Institute of Technology–Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) model (Chin et al., 2000; Ginoux et al., 2001) as implemented in WRF-Chem v4.1.

The AFWA dust emission equations assume that wind-driven dust entrainment primarily occurs through a process called saltation bombardment, in which wind-lofted sand-sized particles (∼ 50–2000 µm diameter) too heavy to remain suspended in the air collide with the land surface and eject smaller dust-sized particles (generally <20 µm diameter) upon impact (e.g., Gillette, 1977; Kok et al., 2012). Several studies have investigated the use of the AFWA dust emission module for a variety of dust modeling applications (e.g., Aragnou et al., 2021; Francis et al., 2022; Hamzeh et al., 2021; Karumuri et al., 2022; Kuchera et al., 2021; Mesbahzadeh et al., 2020; Miller et al., 2021; Mohebbi et al., 2019; Péré et al., 2018; Saidou Chaibou et al., 2020; Solomos et al., 2018; Teixeira et al., 2016; Tsarpalis et al., 2018, 2020; Uzan et al., 2016; Xu et al., 2017; Zhang et al., 2022; Zhou et al., 2019). In general, these studies highlight useful applications of the model. However, several authors noted the need for improved dust source and land surface characterizations in regions with heterogeneous terrain (Cremades et al., 2017; Fountoukis et al., 2016; Hyde et al., 2018; Kim et al., 2021; Ma et al., 2019; Mohebbi et al., 2020; Nabavi et al., 2017; Nguyen et al., 2019; Nikfal et al., 2018; Parajuli et al., 2019, 2020; Parajuli and Zender, 2018; Rizza et al., 2016, 2021; Spyrou et al., 2022; Su and Fung, 2015; Yuan et al., 2019; Zhao et al., 2020).

Here, we aim to evaluate the sensitivity of dust transport simulated by WRF-Chem to the Chappell and Webb (2016) albedo-based drag partition. Our analysis focused on a convective dust event case study from 3–4 July 2014 for the southwestern US desert region previously discussed in other published works (e.g., Hyde et al., 2018; Yu and Yang, 2016). The results from Hyde et al. (2018), in particular, motivated our case study choice. They found that the AFWA dust emission scheme used in WRF-Chem grossly overestimated dust transport for this event. We hypothesized that better geographic sheltering-effect representation from the Chappell and Webb (2016) drag partition would improve the overall magnitude and spatial pattern of the emitted dust and ultimately result in an improved simulation of transported dust. This hypothesis is further strengthened by recent findings from Bukowski and van den Heever (2022), who found that accurate roughness effect characterizations are critical for predicting dust patterns associated with cold pool events similar to our chosen case study.

2 Methodology

2.1 Dust emission flux calculation

The AFWA dust emission module (LeGrand et al., 2019) is an adaptation of the dust emission scheme originally described by Marticorena and Bergametti (1995). Equations comprising the AFWA code are in terms of u_* and primarily include the threshold friction speed required for particle entrainment over an idealized smooth surface (u_*ts), the horizontal saltation flux (Q), the bulk vertical dust flux (F_B), and the size-resolved emitted dust flux. For a detailed overview of the dust emission module equations, see LeGrand et al. (2019). The following two subsections summarize the main components of the AFWA dust emission model and the modifications required to incorporate the Chappell and Webb (2016) drag partition.

2.1.1 The original AFWA module configuration

In the AFWA dust emission calculation, soil particles are divided into a predetermined number of bins based on their effective particle size (referred to as size bins). Tables 1–2 provide attributes associated with the nine size bins used for saltation-based processes and the five size bins used for emitted dust. Here, we denote effective diameters for saltation and dust particles by D_s,p and D_d,p, respectively.

Table 1 Saltation particle size bins and their associated effective diameters (adapted from LeGrand et al., 2019, Table 1). Particle diameters are presented here in micrometers but handled in units of centimeters within the model.

Table 2 Dust particle size bins and their associated attributes (adapted from LeGrand et al., 2019, Table 2). Particle diameters are presented here in micrometers but handled in units of centimeters within the model.

As the simulation evolves, saltation for a given size bin initiates and ceases as u_* exceeds or falls below size-resolved values of u_*ts, respectively. First, the module generates a semi-empirical u_*ts estimate (in units of centimeters per second) for each saltation size bin, $u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}}\right)$, assuming air-dry soil conditions:

$$\begin{array}{}\text{(2)}& u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}}\right)=\mathrm{0.13}{\displaystyle \frac{{\left(\frac{{\mathit{\rho }}_{\mathrm{s},\mathrm{p}}g{D}_{\mathrm{s},\mathrm{p}}}{{\mathit{\rho }}_{\mathrm{a}}}\right)}^{\mathrm{0.5}}{\left(\mathrm{1}+\frac{\mathrm{0.006}}{{\mathit{\rho }}_{\mathrm{s},\mathrm{p}}g{D_{\mathrm{s},\mathrm{p}}}^{\mathrm{2.5}}}\right)}^{\mathrm{0.5}}}{{\left[\mathrm{1.928}{\left(a_{\mathrm{mb}}{\left(D_{\mathrm{s},\mathrm{p}}\right)}^x+b_{\mathrm{mb}}\right)}^{\mathrm{0.092}}-\mathrm{1}\right]}^{\mathrm{0.5}}}},\end{array}$$

where g=981 cm s⁻² is the gravitational constant, ρ_a is the spatiotemporally varying air density from the lowest model atmospheric layer, ρ_s,p is the particle density of saltation bin p, x=1.56, a_mb=1331 cm^−x, and b_mb=0.38. The code then applies a correction function, f(θ), that incorporates soil water content and clay composition to account for the effects of soil moisture on particle cohesion based on the approach of Fécan et al. (1999):

$$\begin{array}{}\text{(3)}& u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)=u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}}\right)f\left(\mathit{\theta }\right).\end{array}$$

Next, the module diagnoses size-resolved saltation flux (Q(D_s,p): in units of grams per centimeter per second) values for each saltation size bin following

$$\begin{array}{}\text{(4)}& Q\left(D_{\mathrm{s},\mathrm{p}}\right)=\left\{\begin{array}{l}C\frac{{\mathit{\rho }}_{\mathrm{a}}}{g}u{{}_{*}}^{\mathrm{3}}\left(\mathrm{1}+\frac{u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)}{u_{*}}\right)\left(\mathrm{1}-\frac{u_{*\mathrm{ts}}{\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)}^{\mathrm{2}}}{{u_{*}}^{\mathrm{2}}}\right),\\ u_{*}>u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)\\ \mathrm{0},u_{*}\le u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)\end{array}\right.,\end{array}$$

where C is an empirical proportionality constant set to 1.0. The module then multiplies Q(D_s,p) by bin-specific weighting factors (dS_rel) determined from prescribed mass distribution assumptions and spatially varying sand, silt, and clay mass fractions. The integrated sum of the weighted Q(D_s,p) values determines the total streamwise saltation flux, Q, associated with each model grid cell:

$$\begin{array}{}\text{(5)}& Q=\sum _{s,p=\mathrm{1}}^{\mathrm{9}}Q\left(D_{\mathrm{s},\mathrm{p}}\right)\mathrm{d}{S}_{\mathrm{rel}}\left(D_{\mathrm{s},\mathrm{p}}\right).\end{array}$$

After determining Q, the module generates the bulk vertical dust emission flux (F_B; in units of grams per centimeter squared per second) triggered by saltation by multiplying the Q field by a topographic-based dust source strength parameter (S) and a sandblasting efficiency factor (β; in units of value per centimeter) derived from the soil clay fraction. The F_B calculation also includes an aerodynamic roughness length (z₀) conditional to limit dust emission to non-urban regions with relatively sparse vegetation coverage:

$$\begin{array}{}\text{(6)}& F_{\mathrm{B}}=\left\{\begin{array}{ll}QS\mathit{\beta },& z_{\mathrm{0}}\le \mathrm{20}\mathrm{cm}\\ \mathrm{0},& z_{\mathrm{0}}>\mathrm{20}\mathrm{cm}\end{array}\right.,\end{array}$$

where

$$\begin{array}{}\text{(7)}& \mathit{\beta }=\left\{\begin{array}{ll}{\mathrm{10}}^{\mathrm{0.136}\left(m_{\mathrm{clay}}\right)-\mathrm{6}},& m_{\mathrm{clay}}<\mathrm{0.2}\\ \mathrm{1.06}\times {\mathrm{10}}^{-\mathrm{6}},& m_{\mathrm{clay}}\ge \mathrm{0.2}\end{array}\right.,\end{array}$$

and m_clay is the soil clay mass fraction. Here, z₀ represents the theoretical height at which the mean wind speed near the surface falls to zero due to surface-induced drag (e.g., Stull, 1988; Zobeck et al., 2003). Importantly, this z₀ parameter is not part of the new drag partition treatment, but rather a value from the parent WRF model used for a variety of land surface and airflow processes and diagnostics.

Effectively, the m_clay conditional in Eq. (7) limits the sandblasting efficiency parameter to 1.06 × 10⁻⁶ cm⁻¹ when the soil composition exceeds 20 % clay content. However, as noted by LeGrand et al. (2019), the β function primarily serves as a dimensional scaling factor. The overall effects of soil clay content variability on F_B are relatively small.

The S parameter (originally described by Ginoux et al., 2001) represents the availability of loose erodible soil material at a given location based on the degree of topographic relief of the surrounding area. This approach assumes that soil composition remains consistent over time, and the simulated land surface will neither run out of dust material nor acquire new dust material through fluvial or atmospheric deposition as the simulation evolves. Some papers refer to this spatially varying S field as an erodibility or dust source map. However, both labels provide an inaccurate description of how S functions in the AFWA module. Accordingly, we will not use that terminology here. Essentially, S is a spatially varying tuning parameter ranging from 0 to 1 that assumes that erodible material accumulates in low points in the terrain, determined by

$$\begin{array}{}\text{(8)}& S={\left({\displaystyle \frac{z_{\max }-z_i}{z_{\max }-z_{\min }}}\right)}^{\mathrm{5}},\end{array}$$

where z_i is the elevation of the cell, and z_max and z_min are the maximum and minimum terrain elevation in the surrounding 10^∘ × 10^∘ area, respectively.

Of note is that the AFWA module uses interpolated values of S initially derived from a $\mathrm{1}/\mathrm{4}{}^{\circ }$ elevation dataset. In addition, the S field incorporates a static vegetation mask that blocks dust emission (i.e., S=0) from vegetated areas derived from a 1^∘ × 1^∘ resolution land cover dataset. While these settings may be appropriate for some modeling applications, the coarse nature of these input datasets likely limits the spatial viability of S at mesoscale and convective-permitting model resolutions (e.g., Saleeby et al., 2019; Vukovic et al., 2014; Walker et al., 2009).

Next, the module applies a prescribed particle size distribution derived using the Kok (2011) brittle fragmentation theory to obtain size-resolved dust emission fluxes:

$$\begin{array}{}\text{(9)}& F\left(D_{\mathrm{d},\mathrm{p}}\right)=F_{\mathrm{B}}\mathit{\kappa }\left(D_{\mathrm{d},\mathrm{p}}\right),\end{array}$$

where F(D_d,p) is the size-resolved emitted dust flux (in units of grams per centimeter squared per second) and κ_d,p is the distribution fraction of the suspended dust size bin (see Table 2). Lastly, the AFWA module uses F(D_d,p) values to determine the size-resolved dust concentrations injected into the first model atmospheric level for transport during each model time step. From this point forward, functions from other modules in the GOCART suite take over processing the fate and transport of the airborne dust aerosols.

2.1.2 Incorporation of drag partitioning

To obtain gridded values of u_ns* for use in WRF-Chem, we estimated surface shadowing using data from the snow-masked 500 m Moderate bidirectional reflectance distribution function (BRDF) albedo daily product (Hall and Riggs, 2016; Schaaf and Wang, 2021a). Following Chappell and Webb (2016) and Ziegler et al. (2020), we derived the normalized proportion of shadow (represented using a normalized albedo, ω_n) by

$$\begin{array}{}\text{(10)}& {\mathit{\omega }}_{\mathrm{n}}={\displaystyle \frac{\mathrm{1}-{\mathit{\omega }}_{\mathrm{dir}}\left(\mathrm{0}{}^{\circ }\right)}{f_{\mathrm{iso}}}},\end{array}$$

where ω_dir(0^∘) is the daily nadir “black-sky” albedo for MODIS band 1 (620–670 nm wavelength; Schaaf and Wang, 2021b) and f_iso is the band 1 BRDF isotropic weighting parameter. In essence, Eq. (10) integrates the albedo across viewing angles for a single illumination angle (solar noon), producing an areal shadow estimate that, in theory, represents non-erodible roughness element conditions for an integrated (500 m pixel) area (e.g., Fig. 1 from Ziegler et al., 2020). We then determined daily u_ns* by

$$\begin{array}{}\text{(11)}& u_{\mathrm{ns}*}=\mathrm{0.0311}\left(\exp {\displaystyle \frac{-{{\mathit{\omega }}_{\mathrm{ns}}}^{\mathrm{1.131}}}{\mathrm{0.016}}}+\mathrm{0.007}\right),\end{array}$$

where ω_ns represents an empirically scaled proportion of shadow obtained via

$$\begin{array}{}\text{(12)}& {\mathit{\omega }}_{\mathrm{ns}}={\displaystyle \frac{\left(a-b\right)\left({\mathit{\omega }}_{\mathrm{n}}-\mathrm{35}\right)}{-\mathrm{35}}}+b,\end{array}$$

using the scaling factors a=0.0001 and b=0.1 to align ω_ns with the ray casting performed on the reconstructed roughness element configurations in the wind tunnel study by Marshall (1971).

Figure 1 A summary overview of the atmospheric forcing conditions associated with a convective dust storm that occurred 3–4 July 2014. The event was characterized by persistent broad high pressure (blue H), clockwise upper-air circulation (black streamlines and arrows), mid- to low-level moisture transport (blue arrows) from the Gulf of California and the Gulf of Mexico, and surface wind vectors (purple arrows) converging downslope of the Mogollon Rim. These conditions initiated and sustained several convective systems that merged near and around Phoenix, Arizona (denoted with a black ×). The shaded overlay shows the national radar reflectivity composite imagery at the storm's peak intensity on 4 July 2014 at 01:00 UTC. Shortly after, the gust front at the leading edge of the convective cell south of Phoenix (dashed black line) moved northwest (storm motion denoted by thin black arrows) over the greater Phoenix area. A wall of thick dust associated with the storm lofted and transported along the gust front.

To incorporate u_s* into the AFWA dust emission module, we configured WRF-Chem to ingest daily MODIS-derived u_ns* fields (Eq. 11) that had been interpolated to the model grid domain into the WRF-Chem framework through an auxiliary channel at model runtime and modified the dust emission equations to use u_s* in place of u_*. We then converted daily u_ns* values to instantaneous u_s* estimates following Eq. (1) using the WRF-Chem simulated wind speed at 10 m a.g.l. that updates each model time step as U_10 m. Finally, we replaced u_* in the saltation function from Eq. (4) with u_s*:

$$\begin{array}{}\text{(13)}& Q\left(D_{\mathrm{s},\mathrm{p}}\right)=\left\{\begin{array}{l}C\frac{{\mathit{\rho }}_{\mathrm{a}}}{g}u{{}_{\mathrm{s}*}}^{\mathrm{3}}\left(\mathrm{1}+\frac{u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)}{u_{\mathrm{s}*}}\right)\left(\mathrm{1}-\frac{u_{*\mathrm{ts}}{\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)}^{\mathrm{2}}}{{u_{\mathrm{s}*}}^{\mathrm{2}}}\right),\\ u_{\mathrm{s}*}>u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)\\ \mathrm{0},u_{\mathrm{s}*}\le u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}},\mathit{\theta }\right)\end{array}\right..\end{array}$$

Note that estimated dust emissions are particularly sensitive to small changes in wind friction speed due to the cubed component of the saltation equation (e.g., Tegen et al., 2002). This replacement of u_* with u_s* in Eq. (13) makes the saltation equation consistent with the physics of aeolian transport in the presence of roughness, where excess wind friction speed at the soil surface (i.e., friction speed above the threshold required for particle mobilization) governs the saltation mass flux (Webb et al., 2020). In addition, we added model runtime configuration options to disable the z₀ conditional and remove the S factor from the F_B function (Eq. 6) since both of these features incorporate broad-scale vegetation masks. Readers are encouraged to review the report by Michaels et al. (2022) for a detailed overview of the code implementation process and model runtime instructions.

The use of the snow-masked version of the MODIS product is essential for this particular modeling application. Although the AFWA module restricts dust emissions from snow-covered areas, the snow-covered grid spaces in the model may not align with snow-covered pixels in the MODIS product. Accordingly, our model implementation process assumes $u_{\mathrm{s}*}=\mathrm{0}$ m s⁻¹ for grid spaces with missing data in the MODIS product. This assumption effectively blocks dust emissions from MODIS-detected snow-covered regions and water bodies. We note the potential for missing data problems if the MODIS retrievals have poor spatial coverage. If this scenario occurs, we recommend that users consider filling the gaps through interpolation techniques (assuming the missing data gaps are minimal), using seasonal or monthly climatological u_ns* values, or using a u_ns* input dataset from a previous period.

2.2 Description of case study event

As described in detail by Gallagher et al. (2022), our case study dust event was driven by a convective outflow boundary associated with a thunderstorm cluster that developed and evolved over central and southern Arizona following convective weather patterns characteristic of the North American monsoon's summer phase (e.g., Adams and Comrie, 1997). The event began at 18:00 UTC on 3 July 2014, peaked at around 00:00 UTC on 4 July 2014, and concluded mainly by 12:00 UTC on 4 July 2014. At the peak of the event, convective cells organized into an extensive north–south line, pivoting clockwise about its northern end. By 4 July 2014 at 02:00 UTC, the convective line collapsed into three components, with each portion progressing in different directions as they diverged. The gust front associated with the central cell just south of Phoenix, Arizona, generated a thick wall of dust that preceded the storm as it continued west-northwest over the metropolitan area. Figure 1 provides a conceptual overview of the general environmental forcing conditions associated with the dust event. For a more in-depth review of the storm evolution, including synoptic, mesoscale, and local condition assessments using a broad collection of analysis fields, radar composites, and observations for support, we encourage readers to review the Gallagher et al. (2022) report.

Our simulation analysis primarily focused on dust produced by the cell that affected Phoenix. However, we also chose this particular case study because dust concentrations outside the immediate Phoenix area were relatively low in magnitude or negligible. The model's ability to simulate minimal dust loading and dust-free (i.e., clear sky) conditions accurately is also important for it to be a reliable prediction tool. For the simulated results to be considered successful for this case study, the model should be able to both capture the progression of the dust wall over Phoenix and limit dust production from the rest of the broader model domain.

2.3 Model configuration

WRF-Chem is a version of the Weather Research and Forecasting (WRF) model by Skamarock et al. (2019) with additional modules for atmospheric chemistry processes and feedbacks (Fast et al., 2006; Grell et al., 2005). Like WRF, WRF-Chem is a fully compressible finite-difference model that simulates atmospheric motion on the Arakawa C-grid and incorporates a variety of parameterizations for simulating subgrid atmospheric motion, cloud microphysics, radiation, and terrain processes. We used WRF-Chem v4.1 for our test case simulation with WRF parent model configuration settings suggested by Gallagher et al. (2022) and chemistry settings from LeGrand et al. (2019) and Letcher and LeGrand (2018). The study by Gallagher et al. investigated the sensitivity of WRF-simulated atmospheric forcing conditions for the dust event studied here. In particular, they focused on the effects of model initialization (spin-up) time, initial atmospheric conditions, horizontal and vertical model resolutions, and several WRF physics package settings to determine the optimal model configuration that minimized environmental forcing condition errors in the dust simulation.

Table 3 provides a brief overview of the model chemistry and physics configuration. However, complete pre-processor and runtime configuration files (referred to as the namelist.wps and namelist.input files, respectively) for this effort are available from the code repository by Letcher et al. (2022) and in the report by Michaels et al. (2022). The model vertical grid used the default spacing distribution with 40 levels following a stretched hybrid sigma–pressure vertical coordinate that favors higher resolution near the ground. We note that the text and Table 2 from Gallagher et al. (2022) erroneously list their model configuration as having 42 vertical levels instead of 40. However, we confirmed that their study simulations were generated with the same vertical level configuration used for our assessment. Figure 2 displays the three telescoping model domains (D01, D02, and D03, hereafter) with grid resolutions of 18, 6, and 2 km, respectively. In Figs. 3–4, we show key terrain attributes associated with the AFWA dust emission functions for D02 and D03, which primarily encompass the southwestern US desert region.

(Han and Pan, 2011)(Ek et al., 2003)(Nakanishi and Niino, 2004)(Nakanishi and Niino, 2004)(Iacono et al., 2008)(Thompson et al., 2008)(LeGrand et al., 2019)

Table 3 WRF-Chem physics and chemistry parameterizations. Note that the option numbers may be different in later model versions.

Figure 2 Telescoping model domains with shading showing terrain elevation in meters. The black × marks the location of Phoenix, Arizona. The elevated terrain feature northeast of Phoenix marked by a line of hollow circles is the Mogollon Rim.

Figure 3 Prescribed land use categories assumed for the model domain. The black × marks the location of Phoenix, Arizona.

Figure 4 Important terrain attribute features considered in the AFWA dust emission flux calculation, including (a) soil sand mass fraction, (b) soil clay mass fraction (m_clay), (c) the clay mass fraction threshold affecting dust emission, (d) dust source strength (S; gray areas are masked), (e) the threshold aerodynamic roughness length (z₀) required for dust emission, and (f) sandblasting efficiency (β). Note that β (panel f) is relatively homogenous due to the clay-rich soil content of the domain and is capped at 1.06 × 10⁻⁶ cm⁻¹, where the soil composition exceeds 20 % clay content. The black × marks the location of Phoenix, Arizona.

We set the WRF-Chem initial and lateral boundary forcing conditions using 3-hourly analysis fields from the Rapid Refresh Model (Benjamin et al., 2016) supplemented with 6-hourly soil moisture information from the National Centers for Environmental Prediction North American Model analysis product. Though these data are obtainable from multiple resources, we acquired all of our forcing datasets from the National Centers for Environmental Information data portal: https://www.ncei.noaa.gov/products/weather-climate-models/numerical-weather-prediction (last access: 2 February 2023).

We performed our simulations for a 2 d period from 3 July 2014 at 00:00 UTC to 5 July 2014 at 00:00 UTC (2 July 2014 at 17:00 LT to 4 July 2014 at 17:00 LT for Arizona). The first 12 h of the simulation (00:00–12:00 UTC on 3 July 2014) were disregarded as spin-up to allow the model to adjust to the initial and lateral boundary conditions. Atmospheric dust was initialized using a “cold start” approach, which assumes an initial atmospheric dust concentration of zero. Settings engaged by the GOCART simple option in WRF-Chem generated atmospheric fields for other non-dust aerosols, including sea salt, black carbon, organic carbon, and dimethyl sulfide. This particular configuration assumes background sea salt emissions based on the lowest model level wind speeds over the oceans (Gong, 2003) and emissions for the other three aerosol species using climatological emission datasets and the WRF-Chem PREP-CHEM-SRC preprocessing software (Freitas et al., 2011). All dust aerosols originated from the local model domain, meaning that no additional dust aerosols entered the D01 lateral boundaries as the simulations evolved. We consider this a reasonable assumption given the relatively short duration of the case study and because the dust event was generated from localized, convection-induced conditions well within the D01 boundaries.

To better understand the sensitivity of the model to the drag partition treatment, we assessed five test configurations, including a simulation produced using the original AFWA code configuration without a drag partition as a control (CTRL) and four alternative versions using the modified size-resolved saltation treatment from Eq. (13) (see Table 4). The first alternative configuration (ALT1) applied the u_s*-based Q(D_s,p) estimates and the original F_B calculation (Eq. 6), which restricts dust emission from grid cells classified by the model as areas with z₀ greater than 20 cm. Because we used the model default 21-class MODIS International Geosphere–Biosphere Programme (IGBP) dataset to characterize land use, this conditional effectively limits dust production to areas classified by the model as barren, grassland, shrubland, savanna, or cropland. The second alternative configuration (ALT2) is identical to ALT1 but with the z₀ conditional removed. Our third alternative configuration (ALT3) further simplifies ALT2 by removing the preferential source strength term from the F_B calculation (i.e., S=1). Lastly, the fourth alternative configuration (ALT4) estimates u_s* by applying a global scaling factor (C_s) to U_10 m:

$$\begin{array}{}\text{(14)}& u_{\mathrm{s}*}=C_{\mathrm{s}}{U}_{\mathrm{10}\mathrm{m}},\end{array}$$

with C_s set to a value within the general range of u_ns* values estimated for the model domain by Eq. (11). Similar to ALT3, the ALT4 configuration also ignores the z₀ conditional and assumes S=1.

Table 4 Size-resolved saltation and bulk vertical dust flux treatments associated with each test case.

The ALT3 configuration tested the ability of the albedo-based approach to accurately represent the spatially varying surface aerodynamic sheltering afforded by non-erodible roughness elements on the land surface. However, we specifically included the ALT4 configuration to test if the model is able to achieve similar results for this case study event by using globally scaled surface winds without the additional input data acquisition and preprocessing requirements inherent to the Chappell and Webb (2016) method. In other words, ALT4 explores if a model user could achieve the same outcome as the Chappell and Webb (2016) drag partition scheme by simply tuning down the winds if the dust source areas in their domain of interest were generally prone to coverage by vegetation or other roughness elements. Effectively, ALT3 and ALT4 are the only model configurations without some form of vegetation mask built into the dust emission code (see Fig. 5).

Figure 5 Erodible and non-erodible areas associated with the prescribed roughness and S parameter masks. Erodible area masks in the plot on the left (a) largely overlap; maroon indicates regions masked by both the z₀ conditional and S parameter, the darker blue represents the area masked by only the S parameter, and orange implies the masking is isolated to the z₀ conditional. Light blue regions in (a) have no masking. The middle panel (b) shows the combined z₀ and S parameter masks applied in the CTRL and ALT1 configurations, and panel (c) shows the S parameter mask applied in ALT2. No masks are applied in the ALT3 or ALT4 configurations. The black × marks the location of Phoenix, Arizona.

Figure 6 provides a comprehensive schematic summary of the five test configurations and their required input parameters. This multi-tiered testing approach enabled us to explore the drag partition's influence on simulation outcomes systematically. The z₀ dependency removal is relatively straightforward since it primarily serves as a vegetation mask. However, omitting the S parameter also removes the spatially varying available sediment supply tuning parameter from Eq. (6). Since the drag partition provides no direct information about sediment supply, it is not clear if removing S is a universally valid approach. We note, however, that previous studies have found that the default WRF-Chem S parameter provides a poor representation of dust source strength in our case study area (e.g., Vukovic et al., 2014; Parajuli and Zender, 2018). Setting S to 1 forces the model to assume that all areas are functionally equally erodible and that dust emissions are solely a result of saltation and sandblasting efficiency. Replacing the built-in WRF-Chem S field without the static vegetation mask, though, is beyond the scope of this effort.

Figure 6 Schematic of the original and modified components of the dust emission module and their required inputs. A black diamond indicates that the parameter varies spatially and temporally. A black circle indicates that the parameter varies spatially but not temporarily, and a hollow diamond means the term is related to a particle size bin attribute table. The hollow circle indicates that the field is derived from MODIS data. These input parameters are provided to the AFWA dust emission module by the WRF-Chem model. Color boxes imply which module components are used in the five different test configurations. See the comprehensive variable list in Table A1 for variable definitions.

Note that β, in this case, is homogenous in the study area (e.g., Fig. 4f) due to the relatively clay-rich soil content of the region (e.g., Fig. 4b) and is capped at 1.06 × 10⁻⁶ cm⁻¹ everywhere the soil composition exceeds 20 % clay content (e.g., Fig. 4c). Accordingly, we can assume that the resultant simulated dust emission patterns produced by our test configurations are a function of Q, S, the z₀ mask, or some combination of these parameters.

2.4 Observation data

We evaluated model performance using in situ observations from Automated Surface Observing Stations (ASOS). ASOS observations are available in standard Meteorological Aerodrome Report (METAR) format (e.g., WMO, 2019) every hour and in the less conventional 1 min or 5 min format, depending on the station. Gallagher et al. (2022) previously assessed model wind speed for our simulation configuration in the vicinity of the main event convection using 1 min frequency ASOS data. Thus, our analysis primarily focused on comparing simulated data against ASOS present weather and visibility observations as well as wind speed measurements for ambient conditions away from the main convective dust event using ASOS data accessed through the Meteorological Assimilation Data Ingest System (MADIS) portal (https://madis.ncep.noaa.gov, last access: 12 November 2021). We also compared simulated results against observed concentrations of particulate matter up to 10 µm in diameter (PM₁₀) using hourly PM₁₀ in situ measurements obtained from the US Environmental Protection Agency's Air Quality System database (retrieved via an application programming interface at https://www.epa.gov/outdoor-air-quality-data, last access: 1 July 2020). We did not consult additional satellite-retrieved products like false color dust-enhanced imagery and aerosol optical depth to supplement our assessment due to cloud obscuration over the area of interest.

To evaluate simulated storm structure, location, and timing, we analyzed observed Next Generation Weather Radar (NEXRAD) composite imagery (available at https://mesonet.agron.iastate.edu/docs/nexrad_mosaic/, last access: 2 February 2023) to qualitatively compare and verify the location, structure, and intensity of the simulated convective storms. The radar composite images comprise a blend of radar reflectivity from individual locations with overlapping coverage to generate a cohesive radar map across the continental United States. This mosaic approach has the advantage of filling in blockages from individual radar sites with neighboring ones, providing a more complete snapshot of convective activity across the radar network coverage area.

3 Results

Experimental results are reviewed as follows. We begin with an overview of the u_ns* field, then review the environmental forcing conditions and dynamic components of the dust emission scheme simulated by the model. Lastly, we assess the resultant dust-related parameters produced by each test configuration outlined in Table 4. This holistic component-based approach helps to break down how the many factors affecting modeled dust conditions contributed to each test simulation outcome.

3.1 Albedo-derived normalized surface friction speed

Figure 7 shows the normalized surface wind friction speed for the case study period diagnosed from MODIS BRDF data and Eq. (11) interpolated to the model domain. Although we consider u_ns* to be a dynamic parameter, u_ns* here is a static field (i.e., u_ns* only varies spatially, not temporally) due to the relatively short duration of our case study simulation (48 h). Note that the u_ns* values in the areas characterized by the model as shrubland range from 0.0225 to 0.0325, barren areas are around 0.03 to 0.035, croplands are 0.0125 to 0.02, and forested and urban areas generally range from 0.005 to 0.0225. The upper limits of u_ns* associated with forested and urban areas are likely higher than those of croplands due to roughness element height consistency and density. For example, plants in individual crop fields generally grow at relatively uniform rates and spacing, giving the crop canopy a “smoother” aerodynamic roughness than urban centers and forested regions with more heterogeneity in roughness element height, shape, and density.

Figure 7 Normalized surface friction speed (u_ns*) diagnosed from MODIS data. Due to the relatively short duration of our case study event, the u_ns* field remains fixed throughout the simulation. The black × marks the location of Phoenix, Arizona.

3.2 Environmental forcing conditions

We conducted a rigorous evaluation of the simulated environmental forcing conditions so we could account for any wind flow errors in our assessment of simulated PM₁₀ against hourly station-based PM₁₀ measurements. The model was able to reproduce the storm's general structure and timing, including the formation of the initial quasi-linear convective system and the collapse of the convective line into individual cells (results consistent with findings by Gallagher et al., 2022). Furthermore, the simulated near-surface wind speeds were in good agreement with wind speeds observed at ASOS. However, simulated wind speeds peaked 1 to 2 h early in some locations with slightly higher (about +1 m s⁻¹) intensity. According to Gallagher et al. (2022), these minor wind speed errors may be partly due to erroneous land use characterization, particularly in the higher terrain elevation areas where the storm initiated. Gallagher et al. (2022) also performed a full statistical analysis of simulated surface wind speeds against all available ASOS wind speed data from the innermost domain (D03). The average wind speed bias for the entire forecast period was +0.59 m s⁻¹. However, a large portion of this overestimation occurred during non-convective nocturnal periods (3 July at 05:00–15:00 UTC and 4 July at 08:00–16:00 UTC).

To better assess the potential influence of these minor errors in wind flow on simulated dust concentrations, we conducted a more in-depth review of the location and timing of the primary cell that produced the main dust event. Specifically, we examined locational differences between a central point along the storm's leading edge in both the observed and modeled radar reflectivity over the dust-producing cell's lifetime (Fig. 8). Overall, the simulated convective activity developed and propagated more northwest than the patterns we observed in the composite radar images. This divergence began with a slight difference in where the cells initiated, followed by a more prevalent westward motion in the simulated conditions than the storm's actual southwestern track. As the storm's evolution continued, the simulated event maintained its path and velocity, while the observed event pivoted and accelerated. Though variation between the simulated and actual storm tracks occurred, this distance became smaller over time (Fig. 8c).

Figure 8 Observed and simulated radar reflectivity for 4 July 2014 at 01:00 UTC, including an (a) observed composite radar image and (b) the model-simulated radar image. The black × marks the location of Phoenix, Arizona. Triangle and square markers indicate the location of a central point along the leading edge of the convective cell that affected Phoenix in the observed and simulated storms, respectively. The plot on the right (c) shows how these leading-edge central point locations evolved over time. Marker colors denote time in UTC.

Due to these discrepancies, we must limit our simulated dust assessment to a qualitative evaluation of the overall storm system rather than comparing hourly simulated time series of PM₁₀ to observed values associated with an exact geospatial point. This approach, however, is not inconsistent with convective weather analyses since simulating the precise placement and timing of individual convective cells is generally considered more luck than skill due to irreducible limitations in nonlinear model physics, numerical processing, and initialization accuracy. Capturing the general spatial pattern and timing of the dust event is the higher priority as model users can often correct persistent magnitude errors through global model tuning parameters.

Simulated time series of U_10 m, u_*, and albedo-derived u_s* values associated with the peak of the dust event (4 July 2014, 00:00–03:00 UTC) are provided in Fig. 9. Overall, the U_10 m patterns align well with the ASOS observations and areas of convective activity seen in the radar composite imagery (see Fig. S1 in the Supplement). We note that the simulated u_* values are generally an order of magnitude stronger than their associated u_s* partition. This outcome is likely due to the drag partitioning scheme interpreting the relatively “dark” terrain surfaces of the domain landscape (presumably caused by the prevalence of shrubs and grasses or trees) as areas with substantial roughness element coverage.

Figure 9 Times series of simulated (a–d) 10 m wind speed (U_10 m; m s⁻¹), (e–h) total wind friction speed (u_*; m s⁻¹), (i–l) surface friction speed determined from the albedo-based drag parameterization (u_s*; m s⁻¹), (m–p) threshold friction speed for a 69 µm diameter particle for an idealized smooth, dry, and barren soil surface ($u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm); m s⁻¹), (q–t) the soil-moisture-based correction function (f(θ); dimensionless), (u–x) and soil-wetness-corrected threshold friction speed for a 69 µm diameter particle ($u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ); m s⁻¹) from 00:00 to 03:00 UTC on 4 July 2014. The u_{*ts_dry} and u_{*ts_wet} symbols used in the figure labeling represent $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm) and ($u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ), respectively. The wetness-corrected threshold is the product of the threshold friction speed and the soil-moisture-based correction function. Note the differences in color bar scaling. The black × marks the location of Phoenix, Arizona.

Figure 9 also shows time series of the threshold friction speed for air-dry soil conditions, the soil-moisture-based correction function, and the soil-wetness-corrected threshold friction speed for a particle with a diameter of 69 µm ($u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm), f(θ), and $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ), respectively). This particular particle size is the effective diameter for saltation size bin 7 (see Table 1). As described in a review paper by Kok et al. (2012, and references within), the u_*ts parameter represents the minimum surface wind shear stress needed to overcome the inertial, cohesive, and adhesive forces holding the particle to the soil bed. Interparticle forces (e.g., capillary, Van der Waals, electrostatic, chemical binding, and Coulomb forces) keeping the particles fixed to the ground have a stronger effect on smaller particles due to their high ratio of surface area to volume. However, larger particles require more energy from the wind to mobilize because inertial force effects become stronger with increasing particle mass. The $u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}}\right)$ curve produced by Eq. (2) has a minimum of about 74 µm (assuming standard atmospheric pressure with 1.225 kg m⁻³ air density) because fine sand grains are less subject to interparticle forces than smaller clay- and silt-sized particles but are still relatively lightweight. Hence, fine sand-sized particles tend to be the first particles to mobilize from direct wind shear stress. Accordingly, particles associated with saltation size bin 7 ($D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm) in the AFWA module will be the first to mobilize. Thus, we consider spatial patterns related to saltation size bin 7 to be good diagnostic indicators of overall dust emission model behavior.

Air density is the only spatiotemporally varying parameter in the calculation of $u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}}\right)$. Though we can discern a slight reduction in $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm) over time immediately under the convective line (Fig. 9m–p), the overall effect of air density on $u_{*\mathrm{ts}}\left(D_{\mathrm{s},\mathrm{p}}\right)$ for this case is relatively negligible. Under air-dry soil conditions, $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm) ranges between 0.17 and 0.19 m s⁻¹ across most of the domain. These results also align with findings by Darmenova et al. (2009) in their assessment of the sensitivity of the Marticorena and Bergametti (1995) dust emission scheme to uncertainties in its required input parameters.

We also see relatively little change in the f(θ) field during the dust event, except for the area associated with a line of precipitation that occurred within the convective cell behind the main wall of dust (Fig. 9q–t). The $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ) maxima over the Mogollon Rim align with isolated areas of convective precipitation that occurred earlier in the simulation (Fig. 9u–x). Note, however, that the $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ) values along the Mogollon Rim adjacent to these maxima are around 0.2 to 0.3 m s⁻¹, comparable to $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ) values in the southwestern Arizona region where the dust event occurred and, for the most part, well below the simulated values of u_*. This particular aspect of the model behavior is important because the only component of the AFWA module preventing dust emissions from these drier forested areas in the ALT3 and ALT4 simulations is the drag partition treatment.

Figure 10 shows time series of the excess friction speed, u_*ex, for saltation size bin 7 associated with different treatments for particle mobilization during the peak of the event, where u_*ex is the model simulated friction speed driving the saltation and dust emission equations (i.e., u_* or u_s*) minus $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ). The six examples provided in Fig. 10 include u_*ex calculated using the model-generated u_* values, u_s* values derived using the Chappell and Webb (2016) method (Eq. 1), and four instances of u_s* values estimated from 10 m wind speeds scaled by a global tuning constant (C_s; Eq. 14) as the wind forcing parameter. We considered values of C_s=0.02, 0.025, 0.03, and 0.035 based on the general range of u_ns* values diagnosed for the area where the dust event occurred (e.g., Fig. 7; note that u_ns* has a theoretical maximum of 0.0381). Effectively, positive u_*ex values indicate that wind friction speeds are strong enough in the model to initiate particle mobilization, while negative u_*ex values imply that the wind friction speeds are too weak. Immediately, we see that practically the entire domain is capable of lofting dust in the u_*-driven simulations if sufficient dust and saltator material are present (Fig. 10a–d). Conversely, dust production in the u_s*-driven simulations is restricted to isolated areas, including areas near the leading edges of individual storm cells (e.g., Fig. 10e–x). The u_*ex results generated with globally scaled U_10 m and C_s=0.025 (Fig. 10e–h) are markedly similar to the u_*ex values calculated with U_10 m scaled by u_ns* (Fig. 10q–t). In particular, the excess friction speeds for these two instances match rather closely over Arizona, with minor spatial extent differences for positive u_*ex values in California and southern Nevada. Thus, we chose C_s=0.025 as the U_10 m global scaling constant for our ALT4 test configuration. Interestingly, areas in the u_ns* field equal to 0.025 (Fig. 7) tend to align with grid cells characterized as closed shrubland (Fig. 3) in the model domain where there were strong winds, though this may be a coincidence.

Figure 10 Times series of excess friction speed (u_*ex; (u_*ex; m s⁻¹) calculated by subtracting the friction speed for a 69 µm diameter particle, $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ), from the wind forcing parameter driving dust production from 00:00 to 03:00 UTC on 4 July 2014. The u_{*ts_wet} symbol used in the figure labeling represents $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ). Plots provided include (a–d) u_*ex calculated using total friction speed, (e–h) surface friction speed determined from the albedo-based drag parameterization, and (i–x) four instances of 10 m wind speeds scaled using a global tuning constant as the wind forcing parameter. The black × marks the location of Phoenix, Arizona.

To further explore how our code modifications influenced simulated outcomes, we diagnosed time series of the horizontal saltation flux for the peak of the dust event associated with each test configuration (Fig. 11). Though the non-erodible area masking components of the AFWA dust emission code are not applied in the AFWA module until the F_B calculation (Eq. 6), we included the appropriate non-erodible area masks for the CTRL, ALT1, and ALT2 configurations (see Fig. 5) in the Fig. 11 Q diagnostics. Strong saltation occurred everywhere in the CTRL configuration plots unless the grid spaces were masked out (e.g., Fig. 11a–d). The four alternative model configurations showed Q patterns similar to each other, suggesting that the masking had little influence on the resultant Q flux in AFWA module configurations that accounted for drag partition. The ALT1 and ALT2 Q outcomes were practically identical because the z₀>20 cm mask largely coincides with the areas where S=0 for this domain (e.g., Fig. 5). With the exception of a small area of saltation to the east of Phoenix, the ALT3 Q results are also similar to ALT1 and ALT2, which makes sense as developed and urban areas as well as higher-elevation forested regions generally yield low values of u_ns*. Lastly, the ALT4 Q spatial patterns mimic ALT3. However, the magnitude of the ALT4 Q flux in the areas outside where the main dust event occurred is lower than the ALT1, ALT2, and ALT3 results.

Figure 11 Times series of horizontal saltation flux (Q; µg m⁻¹ s⁻¹) for each test configuration (CTRL, ALT1, ALT2, ALT3, and ALT4) with non-erodible area masking applied from 00:00 to 03:00 UTC on 4 July 2014. The black × marks the location of Phoenix, Arizona.

3.3 Simulated dust conditions

We compared our simulation results against observed PM₁₀ concentrations (Figs. 12–14). To ensure our simulated PM₁₀ conditions could primarily be attributed to dust transport, we also reviewed the general contribution of each aerosol species to the overall simulated PM₁₀ patterns (not pictured). For this case, most of the simulated PM₁₀ came from either dust or sea salt transport. Sea salt and dimethyl sulfide distributions were relatively isolated to areas over the ocean and the immediate coastlines, while dust was the primary source of PM₁₀ over land. Contributions from black carbon, organic carbon, and dimethyl sulfide to simulated PM₁₀ totals were negligible, except in dense urban areas along the southern California coastline.

Figure 12 Times series of observed and simulated PM₁₀ concentration for each test configuration (CTRL, ALT1, ALT3, and ALT4) (µg m⁻³) from 00:00 to 03:00 UTC on 4 July 2014. The black × marks the location of Phoenix, Arizona. The ALT2 configuration results are not included in this figure since they are nearly identical to results from the ALT1 configuration.

Figure 13 (a) Time series of observed (OBS) and simulated (CTRL, ALT1, ALT3, and ALT4) maximum PM₁₀ concentration (µg m⁻³) for the combined Maricopa and Pinal counties in Arizona. Panel (b) highlights how most EPA monitoring stations in Arizona are in Maricopa county and Pinal county, clustered around the Phoenix metropolitan area (marked by the black ×). We used the hourly maximum PM₁₀ concentrations from the combined Maricopa and Pinal county areas for our assessment due to the shifted nature of the main dust-producing outflow boundary recorded by the EPA stations situated around Phoenix, Arizona, in the simulations. The map zoom in (c) of the Maricopa county and Pinal county boundaries shows the EPA station locations relative to the hourly PM₁₀ maximum values sampled for each test configuration and the observations on 4 July 2014 from 00:00–03:00 UTC. The ALT2 configuration results are not included in this figure since they are nearly identical to results from the ALT1 configuration. Note that the CTRL, ALT1, and ALT3 PM₁₀ maximum value sites overlap on 4 July 2014 at 01:00 UTC, and the ALT1 and ALT3 sites are co-located on 4 July 2014 at 02:00 UTC.

Figure 14 Simulated and observed (a) surface wind speed and (c) PM₁₀ values for ASOS and EPA station locations in the Mojave Desert area associated with the erroneously simulated PM₁₀ plume. The spatial plot (b) shows the station locations relative to the area of maximum wind speed simulated by the model over the Mojave Desert on 4 July 2014 at 02:00 UTC, when observed PM₁₀ values for the general location peaked.

Note that PM₁₀ monitoring stations in the United States are generally situated around populated areas, so most of the PM₁₀ sites associated with our model domain were tightly clustered around Phoenix (e.g., Fig. 13b). The dust wall crossed over the Phoenix metropolitan area between 00:00 and 04:00 UTC on 4 July 2014. Importantly, the air quality conditions rapidly recovered once the main dust wall passed, and the areas outside the immediate dust event were relatively clear. According to surface weather station visibility, present weather records, and observer remarks for Arizona and southern California (not pictured), it is unlikely the strong dust obscurations extended much beyond the immediate convective cell outflow boundary area.

The original CTRL configuration produced a strong dust signal along the convective outflow boundary, but it also generated widespread dust concentrations in areas that in reality were clear (e.g., Fig. 12). Similar to the Q flux patterns, simulated PM₁₀ concentrations produced by the ALT1 and ALT2 configurations were nearly identical because dust in the u_s*-driven simulations did not originate near a masked area (see Fig. 5). Due to this artifact, we opted to limit our PM₁₀ assessment to the CTRL, ALT1, ALT3, and ALT4 simulations. The ALT1 configuration produced much lower PM₁₀ concentrations than the other test settings, with some dust emerging near the gust front. The ALT1 “dust wall”, however, was an order of magnitude lower than the observed conditions. In contrast, the magnitude and spatial pattern of the simulated PM₁₀ values associated with ALT3 and ALT4 aligned better with reported observations.

Though the general PM₁₀ patterns of the main dust event are similar between ALT3 and ALT4, the areas of maximum PM₁₀ extend further south along the gust front in the ALT3 simulation than in the ALT4 version (Fig. 12). We also note higher PM₁₀ values for ALT3 in southern California over the Mojave Desert and near the northeastern corner of Arizona compared to ALT4. However, the maximum PM₁₀ values associated with these two areas occurred in observation-limited regions. PM₁₀ observations from an EPA station in Barstow, California (site ID: 06-071-0001; 34.8939080^∘ N 117.024804^∘ W), near the simulated area of maximum PM₁₀ over the Mojave Desert suggest that simulated PM₁₀ values in both ALT3 and ALT4 were higher than observed at this location (Fig. 14c). For example, simulated PM₁₀ values for ALT3 and ALT4 at the Barstow site peaked at 570 and 346 µg m⁻³, respectively, at 00:00 UTC on 4 July 2014, while observed PM₁₀ values ranged between 5 and 47 µg m⁻³ during the 4 July 2014 00:00–04:00 UTC period. However, the PM₁₀ concentrations in the simulations changed rapidly over short distances near the Barstow EPA station as the simulations evolved. Thus, small shifts in the simulated dust position could greatly affect the apparent skill of the simulated output.

Surface wind speed observations from Edwards Air Force Base (KEDW, K9L2), which is primarily upstream and west of the simulated Mojave Desert PM₁₀ plume, and Barstow (KDAG) suggest that there were strong wind speeds over the area in question (Fig. 14). Observed wind speeds at the Barstow ASOS station closest to the simulated PM₁₀ plume ranged from ∼ 5 to 10 m s⁻¹ between 00:00 and 03:00 UTC on 4 July 2014, which aligns relatively well with the modeled 10 m wind speeds. Though the simulated wind speed ramp-up for the Barstow location was slightly early on 3 July 2014, the observed strong winds persisted longer over a similar time span before decaying. Accordingly, we can conclude that the model reproduced a realistic wind evolution pattern over the Mojave Desert. While the actual dust conditions for the area are less clear, the lack of commentary about dust in the recorded ASOS observer remarks leads us to believe that the simulated PM₁₀ values for both the ALT3 and ALT4 configurations in the Mojave Desert region were higher than what occurred. However, the cause of this discrepancy is unclear. The erroneous PM₁₀ feature could be due to issues with the drag partition treatment, unresolved sediment supply influences, a combination of the two factors, or some other source of error within the model.

Though we are unable to do point-specific quantitative evaluations for the Phoenix metropolitan area due to the shifted location of the main dust-producing outflow boundary in the simulations, we can see from our spatial PM₁₀ time series analysis (e.g., Fig. 12) that the maximum PM₁₀ value during the peak of the dust event was generally co-located with the storm outflow boundary. Thus, we compared the simulated maximum PM₁₀ values for the area surrounding Phoenix with the maximum PM₁₀ values observed by the EPA stations as the dust wall passed over the city. Figure 13 shows time series of observed and simulated maximum PM₁₀ concentration for the combined Maricopa and Pinal counties in Arizona. This combined county area generally encompasses the footprint of both the simulated and observed track of the storm. Note that the gust front in the simulation advanced into the combined county area earlier than the recorded event (see Fig. 8), and the main dust wall moved away from the EPA stations around Phoenix after 03:00 UTC. Therefore, we expect the simulated combined county maximum PM₁₀ values to be higher than the observed PM₁₀ conditions before and after the observed PM₁₀ peak at 02:00 UTC on 4 July 2014. During the 00:00 to 03:00 UTC period, the observed maximum PM₁₀ was 9132 µg m⁻³, while the CTRL, ALT1, ALT3, and ALT4 configurations produced maximum PM₁₀ values for the combined county area of 17 151, 419, 8539, and 7240 µg m⁻³, respectively. If we emphasize the order of magnitude over the exact value of PM₁₀ from these results to account for observed and simulated gust front location differences, both ALT3 and ALT4 produced reasonable results.

4 Discussion

Findings from this case study show that incorporating a drag partition treatment into the AFWA dust emission module can improve simulated dust transport for vegetated dryland regions. The CTRL simulation results with no drag partitioning align with previous findings that the default AFWA module settings produce excessive dust conditions in the southwestern United States. While the albedo-based drag partition method improved the results with respect to PM₁₀ spatial patterns and concentration magnitudes, we also achieved similar outcomes with a global wind scaling parameter.

We argue that the ALT4 globally scaled wind speed approach, at least how it is applied in our case study, should also be considered a crude form of drag partitioning. However, unlike the Chappell and Webb (2016) method, this global drag partitioning approach assumes that the effects of the variability in vegetation cover over the dust-producing areas within the model domain are relatively negligible for the dust entrainment process. Our case study results imply that the spatial variability of u_ns* was not a key factor in the evolution of the dust storm simulation, but this outcome may simply be an artifact of the conditions associated with our chosen case study. This particular dust event occurred in a confined, localized area due to its convective nature over a shrub-dominated landscape with relatively homogenous topography (e.g., Figs. 2–3). The spatial variability of u_ns* in the broader region may have more influence on dust events generated by synoptically driven wind forcing conditions (i.e., horizontal length scales on the order of several hundred kilometers to 1000+ km). Furthermore, the ALT4 configuration does not provide a viable modeling solution for users interested in exploring the effects of land management or land cover change on sediment transport and air quality and may be limited in its use to short-term atmospheric forecasting applications.

It appears that drag partition treatments effectively suppressed the broad-scale erroneous dust emissions in southern California and southwestern Arizona inherent to the CTRL configuration without restricting the dust emissions that generated the Phoenix haboob. However, the relatively low PM₁₀ concentrations produced by the ALT1 simulation suggest that the S parameter was masking the absence of drag representation in the original AFWA code. Although the model performance for this case study is arguably better with the drag partition when S=1, this may not be the case for all regions. Further investigation of the role of the S parameter in areas with more heterogeneous sediment supply is necessary to resolve this issue.

In order to better isolate the role of the u_ns* spatial heterogeneity in our results from the effects of soil moisture on the dust emission process, we reproduced the diagnostic excess friction speed time series from Fig. 10 but used the threshold friction speed for air-dry soil conditions for saltation size bin 7 instead of the moisture-corrected threshold value (u_{*ex_dry}; Fig. 15). The resultant u_{*ex_dry} values for all wind forcing treatments are, at least to some degree, largely higher in magnitude than their moisture-corrected counterparts in Fig. 10. However, both the u_*ex and u_{*ex_dry} values for the forested areas along the Mogollon Rim are generally low, even when using values for C_s that align with the albedo-derived u_ns* estimates associated with the barren and shrubland areas. For example, the u_ns* values for the forested Mogollon Rim region in Fig. 7 are primarily less than 0.02 throughout the peak of the dust event. This suggests that model sensitivity for this domain to the albedo-based drag partition may vary for alternate wind flow conditions and that the albedo-based drag partition scheme could potentially benefit from additional validation studies.

Figure 15 Times series of excess friction speed (u_*ex; m s⁻¹) calculated by subtracting the friction speed for a 69 µm diameter particle, $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ), from the wind forcing parameter driving dust production from 00:00 to 03:00 UTC on 4 July 2014. The u_{*ts_wet} symbol used in the figure labeling represents $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm, θ). Plots provided include (a–d) u_*ex calculated using total friction speed, (e–h) surface friction speed determined from the albedo-based drag parameterization, and (i–x) four instances of 10 m wind speeds scaled using a global tuning constant as the wind forcing parameter. The black × marks the location of Phoenix, Arizona. Excess friction speed (u_{*ex_dry}; m s⁻¹) calculated by subtracting the dry threshold friction speed for a 69 µm diameter particle, $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm), from the wind forcing parameter driving dust production from 00:00 to 03:00 UTC on 4 July 2014. The u_{*ex_dry} symbol used in the figure labeling represents $u_{*\mathrm{ts}}(D_{\mathrm{s},\mathrm{p}}=\mathrm{69}$ µm). Plots provided include (a–d) excess friction speed calculated using total friction speed, (e–h) surface friction speed determined from the albedo-based drag parameterization, and (i–x) four instances of 10 m wind speeds scaled using a global tuning constant as the wind forcing parameter. The black × marks the location of Phoenix, Arizona.

We note, however, that modeled 10 m wind speeds will be reduced over forested areas due to the aerodynamic drag from the trees on simulated wind speeds and that these particular aerodynamic roughness effects vary as a function of the z₀ settings in the parent WRF model. The drag partition treatment, at least with respect to how it is configured in our modeling setup, has no influence on simulated winds outside the dust emission calculation. Even if the u_ns* values were high for a model grid cell with forest cover, the simulated winds would likely be reduced by the internal WRF model physics, potentially mitigating (or masking) a shortfall in the drag partition correction.

These peculiar results over the forested area also highlight the need for consistency of aerodynamic roughness treatments between the dust emission module and other components of the parent weather model, including the land surface and planetary boundary layer schemes, so that u_* and u_s* diagnosed for the dust emission treatment are consistent with the other model elements governing simulated wind flow. Findings by Webb et al. (2014b) demonstrate how these kinds of inconsistencies can lead to erroneously modeled increases in sediment transport with increasing vegetation coverage. As dust emission models like the AFWA module are further refined, careful consideration should be taken when configuring model components that incorporate land surface conditions on aerodynamic processes.

Even though we were able to largely replicate the PM₁₀ case study simulations produced with the albedo-based drag partition by using a global scaling factor, it is important to remember that our C_s value selection was informed by the albedo-derived u_ns* values and that the relatively short duration of our case study event meant that the u_ns* values associated with our study were static throughout the event life cycle. In any event, choosing an appropriate C_s value a priori could be challenging when using the globally scaled surface wind speed approach for operational forecasting, climatological applications, or a larger domain with more variability in vegetation coverage over dust-producing areas.

Spatial u_ns* patterns derived from MODIS data will likely evolve over longer time spans from processes like annual vegetation growth cycles, drought conditions, anomalous wet periods, and land cover change. To investigate, we reviewed the temporal coefficient of variation (CV) of 500 m gridded monthly mean snow-masked u_ns* values in our case study area over the 2001–2021 MODIS record (Fig. 16) and found that CV values generally ranged between 1 % and 12 % for grid pixels in the dust-producing regions of southwestern Arizona, with a few clusters of pixels ranging between 15 % and 20 % over active agricultural areas. By comparison, with the global C_s drag partitioning (ALT4), we saw markedly different spatial patterns in areas of the model domain capable of initiating saltation simply by varying the C_s parameter by increments of 0.005 (e.g., Figs. 10 and 15), or $C_{\mathrm{s}}=\mathrm{0.025}\pm \sim \mathrm{20}$ %. Given that u_ns* and C_s are mathematically equivalent with respect to the model implementation (i.e., in terms of being multiplied by U_10 m to generate u_s* in Eqs. 1 and 14, respectively), our MODIS-derived u_ns* CV analysis suggests that C_s=0.025 may be a reasonable setting for this domain in general. Nonetheless, though these CV values may appear to indicate a relatively low amount of variability in MODIS-derived u_ns* over the majority of the dust-producing region, even modest changes in u_ns* may result in simulated u_s* values exceeding or falling below u_*ts.

Figure 16 The spatiotemporal variability of u_ns* for the 2001–2021 MODIS record shown as the coefficient of variation of monthly u_ns* as a percentage. The black × marks the location of Phoenix, Arizona. Note that the color bar tick marks increase in increments of 5 % after the 30 % mark.

Thus, the degree to which the observed variability in u_ns* matters from a mesoscale dust modeling context is still unclear. Okin (2005) explored the influence of surface heterogeneity on saltation flux following the original method of Raupach et al. (1993), whereby the drag partition adjusts u_*ts directly rather than modifying the u_* value driving the saltation flux equation. When reviewing the daily saltation flux as a function of the CV of u_*ts, Okin (2005) found that CV values of 10 % or less produced minimal differences in saltation flux. Additional modeling case studies focused on seasonal and interannual land cover changes are needed to determine the sensitivity of the AFWA dust emission module to u_ns* variability.

Based on these findings, we argue that the Chappell and Webb (2016) drag partition approach may have the potential to improve dust transport model performance in vegetated dryland regions. Still, it is worth discussing the nuances associated with its use. On the one hand, the technique removes the dependence on greenness-related remote sensing indicators. By using albedo in a way that is interpreted as shadow to discern roughness elements and their associated aerodynamic properties, the Chappell and Webb (2016) technique may capture the drag induced by plants as well as drag generated by engineered and natural structures like fences, tillage furrows, and subgrid-scale topographic relief. In theory, this method incorporates subgrid-scale heterogeneity and temporal variability in its calculation.

Another aspect to consider is the use of 10 m wind speed in Eq. (1), which is not physically based. The original derivation for u_s* multiplied u_ns* by a parameter Chappell and Webb (2016) referred to as “the free-stream wind speed” or U_f (also referred to by some studies as U_h). Yet this U_f parameter is not associated with wind speed in the free atmosphere (i.e., wind speed above the atmospheric boundary layer). Instead, the value stems from the free-stream flow associated with the wind tunnel data collected by Marshall (1971) that Chappell and Webb (2016) used to develop their equations, which has no direct replacement. Additional research is needed to assess the degree of error introduced by the U_f=U_10 m assumption.

The Chappell and Webb (2016) scheme also assumes that all observed shadows emanate from solid subgrid-scale roughness objects within the flow field. This assumption could be obviated if shadows cast by tall objects extend into neighboring grid cells, creating a sensitivity to wind direction. While this is not likely to be an issue at the mesoscale or coarser modeling scales with u_ns* values derived from 500 m MODIS data, additional lower-bound grid-resolution limitation studies may be necessary before applying the Chappell and Webb (2016) drag partition method to finer-scale dust transport models (e.g., large eddy simulation and field-scale models) with u_ns* initialized from high-fidelity data sources.

Like many other parameterizations, the Chappell and Webb (2016) drag partition is empirical. It incorporates multiple rounds of normalization and draws conclusions about nature from idealized “black-sky” (i.e., 100 % direct illumination) scenarios that do not exist in real-world settings. Furthermore, the Marshall (1971) wind tunnel experiments that Chappell and Webb (2016) used to establish their drag partition equations used non-porous, rigid cylinders to imitate natural roughness elements. Multiple studies have since demonstrated that the cylinder method is often a poor representation the aerodynamic behavior of live vegetation (e.g., Walter et al., 2012). Arguably, the effects of uncertainties in the rescaling parameters associated with the albedo approach on simulated dust entrainment outcomes may be comparable to errors created by uncertainties in inputs and configuration parameters required for other drag partition methods. Additional studies comparing the many available drag partitioning approaches in mesoscale and global-scale dust transport models over extended periods of vegetation green-up and senescence are needed to fully resolve this issue.

Given the relative simplicity of implementing the Chappell and Webb (2016) drag partition into dust emission modeling code, it would be worthwhile to continue investigating use of the scheme in convection-permitting and regional-scale dust transport models. At these scales, many of the aforementioned issues would likely be small or negligible compared to errors introduced by other parts of the modeling framework that govern u_* and erodible sediment supply. In fact, some US agencies have already begun to incorporate the Chappell and Webb (2016) drag partition into their operational air quality models (e.g., Tong et al., 2020).

Finally, though the MODIS data preprocessing steps required by the Chappell and Webb (2016) drag partition are comparable to the processing requirements of other air quality and weather model input datasets, operational environmental prediction agencies may find the additional processing burden prohibitive. The computational resource allocations of production systems often operate near capacity in order to meet a multitude of simulation parameter output requirements and data dissemination schedules, making it challenging to justify added data ingest and preprocessing steps for a single model parameter. Accordingly, future efforts could explore the feasibility of using monthly or seasonal climatologies of u_ns* for weather forecasting applications.

5 Conclusions

This study explored the use of an albedo-based drag partitioning scheme introduced by Chappell and Webb (2016) to incorporate roughness effects into a widely used dust emission and transport model within WRF-Chem. Our results for a convective dust event case study for the desert region around Phoenix, Arizona, support previous findings that the original WRF-Chem/GOCART/AFWA model configuration is prone to producing widespread erroneous dust emissions in the southwestern United States. Furthermore, our findings suggest that incorporating a drag partition treatment into the AFWA dust emission module may improve dust transport simulation and forecasting in vegetated dryland regions through better sheltering-effect representation on dust emissions. This code adaptation could take the form of a formal drag partition treatment or a wind speed tuning approach if the roughness element coverage of the dust-producing areas associated with the model domain is somewhat homogenous and the simulation is run over a short enough period that the surface roughness configuration does not change. We note, however, that we cannot generalize model performance from a single case study event. Instead, we view this study as a preliminary demonstration of capability and motivation to further explore the use of drag partitioning in dust transport models to account for roughness effects on the dust entrainment process. Additional case studies and seasonal analyses are necessary to determine if the Chappell and Webb (2016) drag partition method can reliably improve dust simulations used for various weather prediction, land management, and climate change modeling applications.

The benefits of using a drag partition methodology in the AFWA module are manifest in the vast improvements we see with ALT3 and ALT4 over the initial CTRL configuration. Follow-on studies are still necessary to explore if these benefits persist over longer simulation periods. However, we anticipate that satellite-derived roughness information could markedly improve investigations of the role of short- and long-term changes in vegetation in dust emission patterns. This, in turn, could be of benefit to model users interested in drought hazard, climate change, land management, land use and land cover change, and post-wildfire condition modeling applications.

Appendix A: Variable list

Table A1 provides the symbol, name, units, and value or description of variables referred to throughout this paper. Values of prescribed constants are listed. Variable arrays are described as “variable” for size-bin-related settings, “spatially varying parameter” for static fields, and “spatiotemporally varying parameter” for temporally dynamic fields. See Fig. 6 for a schematic overview of how these parameters are used in the AFWA module calculations.

Table A1 Variable list. Units are provided in terms of mass [M], length [L], and time [T]. In the Type column, “C” implies the parameter is a prescribed constant or global scaling factor, “V” indicates the parameter is a prescribed one-dimensional array of size-bin-related settings, “SV” parameters are spatially varying static fields, and “STV” parameters are spatiotemporally varying dynamic fields.