Effects of Non-Photosynthetic Vegetation on Dust Emissions

Mineral dust is among the top contributors to global aerosol loads. Ability of non-photosynthetic vegetation (NPV) to suppress dust emission has been widely acknowledged but a realistic representation of NPV has not been tested with regional-to-global scale models. In this study, we implemented a satellite-based total vegetation data set, which included NPV, into a regional atmospheric chemistry model and conducted simulations for the year 2016 over the conterminous United States. To test the response of dust simulations to the NPV coverage, we conducted a control simulation incorporating only the photosynthetic vegetation (PV). Simulated dust emissions decrease by 10%–70% over most of the southwestern US from spring to autumn due to NPV. Reductions in dust concentrations are the largest in spring, which attenuate the overpredictions of fine soil concentrations, but accentuate the underpredictions in summer. Overall, the mean errors and correlations of annual simulations are slightly improved with NPV. NPV modulates dust emissions mainly by sheltering the surface and increasing the threshold velocity through drag partitioning. Moreover, we investigated the effect of vegetation height and addressed its uncertainties through a series of sensitivity tests. We observed that a 50% variation in predefined vegetation heights results in small changes in soil concentrations over majority of southwestern US, but causes up to 30% changes at local hotspots. This study highlights the significance of including NPV into the dust model and points out the importance of validation of total vegetation datasets as well as more realistic representation of vegetation heights and seasonality.


of 22
per unit horizontal area (Raupach et al., 1993). A classic work by Hagen (1996) distinguished the erosion-control mechanisms of flat and standing residues and developed equations for modified dust flux with both residues based on compiled wind tunnel data. These equations used coverage of flat residue, as well as height and frontal area index of standing residue to determine emission coefficient, raise threshold wind velocity, and adjust wind profile and subsequently friction velocity. In addition, empirical relations between dust generation and residue density were established, such as the exponential dependence of soil loss measured in a wind tunnel on the frontal area index of three types of crop stubbles (Pi et al., 2020). However, these previous studies on the effect of NPV on dust emissions were confined to wind tunnel experiments or small-scale field measurements over erodible lands with limited vegetation species. Evaluation of these findings at larger scales and over regions with more heterogeneous plant community is of research interest.
The challenge to extend this knowledge to regional-to-global scales lays in providing accurate information about the temporally and spatially variant vegetation to dust emission models. In practice, the parameterization of vegetation in windblown dust schemes implemented in chemical transport models or general circulation models relies on data for vegetation fractional coverage (Duncan Fairlie et al., 2007;LeGrand et al., 2019). Currently, the vegetation fractional maps used by most dust models are based on readily available vegetation indices from satellite retrievals which only represent green or photosynthetic vegetation (PV) (e.g., fraction of absorbed photosynthetically active radiation (FPAR) or normalized difference vegetation index (NDVI)). Most recently, the total vegetation (sum of PV and NPV) data was made available in the Multiscale Online Non-Hydrostatic AtmospheRe CHemistry (MONARCH) model version 2.0 (Klose et al., 2021), but analysis on how NPV modulates dust dynamics using this data set are yet to be conducted.
To the authors' knowledge, there has been only one attempt to quantify the changes in regional dust emissions after NPV was included in the vegetation map (Kang et al., 2014). The researchers estimated the NPV fractions in East Asia by assuming that the NPV fractions follow a linear decrease from the maximum fractions of green vegetation from last year. The conclusion that the approximated NPV fractions improved the simulation for a dust event can be further validated with a more realistic representation of NPV coverage that accounts for the non-linear growth-decay cycle of plants in different environments.
Remote sensing techniques have an advantage in capturing the heterogeneity of vegetation at a larger scale compared to field measurements or ecosystem modeling (Mougin et al., 1995). These techniques can identify NPV elements based on their different reflectance spectra in the visible light to short-wave inferred regions, due to their lower pigments and water contents than PV, and higher cellulose and lignin contents than soils (Z. Li & Guo, 2015). Multispectral imagery is more widely used than hyperspectral imagery at large scales due to availability (Z. Li & Guo, 2015). Several vegetation indices that represent NPV coverage calculated from selected bands of multispectral reflectance have been developed, such as the Normalized Difference Senescent Vegetation Index (NDSVI) (Qi & Wallace, 2002), the Soil Adjusted Total Vegetation Index (SATVI) (Marsett et al., 2006), and the Dead Fuel Index (DFI) (Cao et al., 2010), but they are generally considered to be site-specific (X. Li et al., 2016). Another technique to acquire NPV fractions is the spectral mixture analysis (SMA) which uses all bands of the reflectance spectrum (Asner & Heidebrecht, 2002). The SMA method assumes that the surface reflectance is a combination of the reference reflectance of certain surface components or endmembers, and then resolves the fractions of all endmembers given their reference spectra. Variations of SMA methods (X. Li et al., 2016;Okin et al., 2013) mostly differ in the selection of reference spectra and comparisons among these variations showed that allowing the spectra for a certain endmember to vary among pixels improved the estimation of fractions. Overall, the SMA approach is promising to determine NPV coverage on the surface.
In addition to vegetation coverage, some dust models also use vegetation height to calculate aerodynamic roughness length, which modifies the wind profiles Klose et al., 2021). Although near-global maps for vegetation height based on the Geoscience Laser Altimeter System (GLAS) data were developed and evaluated (Los et al., 2012;Nie et al., 2015), the implementation of vegetation height in dust models remain quite simple. The MONARH model version 2.0 (Klose et al., 2021) preselects a maximum vegetation height based on field studies in Sahel, and then assumes that the height scales with leaf area index (LAI). The Community Multiscale Air Quality (CMAQ) model version 5.3 (Appel et al., 2021; uses a look-up table for monthly vegetation heights which aims to represent their seasonal dynamics and differences among land types. The objectives of this study are to analyze and understand the responses of seasonal dust emissions to NPV coverage and to test the sensitivity of dust emission to vegetation heights over the southwestern United States through a modeling approach. Influence of including NPV coverage was analyzed by implementing satellite-based maps for both PV and NPV fractions derived using an SMA method into the windblown dust scheme in the CMAQ model. Simulations were conducted for the entire year 2016 over the conterminous United States (CONUS) (hereafter the TOTAL run). As a comparison, a control run was conducted which used PV-only coverage data represented by the FPAR from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument (hereafter the FPAR run). The vegetation heights were predefined and incorporated into the dust model, and responses of dust simulations to varying vegetation heights were investigated and quantified through a series of sensitivity tests. Section 2 describes the methods for generating both the vegetation coverage and height data, the parameterization of vegetation in our dust model, the design of sensitivity tests with adjusted vegetation heights, and the evaluation methodology. In Section 3, dust simulations from two runs are compared to determine impacts of NPV coverage and the underlying mechanisms are discussed. Results of the sensitivity of dust model to varying vegetation heights are also analyzed. Finally, summaries of our major findings and proposed future improvements are included in Section 4.

Vegetation Coverage Data
In this study, we used two datasets for vegetation fractional coverage. The total vegetation data was obtained from the MODIS Nadir BRDF-Adjusted Reflectance (MCD43A4) product collection 5, which has a spatial and a temporal resolution of 500 m and every 16 days, respectively. Details about the spectral mixture analysis (SMA) method used to develop this data set have been described in Guerschman et al. (2015) and Scarth et al. (2011). In brief, a linear unmixing method was performed to calculate the fractions of three pure surface components, namely PV, NPV, and bare soils based on the observed surface reflectance and the synthetic reference reflectance of the three components. Synthetic reference reflectance was derived from field measurements of vegetation fractional cover and satellite imagery using multiple linear regression. Seven bands of the reflectance from the MCD43A4 product were used to perform the unmixing, as well as their log transforms and interactive terms to account for the non-linear spectral mixing. To avoid overfitting, a subspace truncation method was applied to control the number of reflectance terms used for unmixing and that number was determined with a 100-fold cross-validation method. During the unmixing, the three fractions in each pixel were constrained to be non-negative and must add up to 100%. The resulting data set includes monthly averaged vegetation fractions at 5 km resolution. It was re-gridded to 12 km over the study domain using the nearest-neighbor space-filling method. The processed monthly data was then transformed into daily data using linear interpolation, and meanwhile, some missing values were replaced using values from adjacent months. For a small amount of grid cells with missing values throughout the year, the total vegetation fractions were set to be 1. The rationale was that the missing values were likely due to snow cover because most of those grids were in the north of the study domain, and complete coverage of vegetation would eliminate the dust emissions from these cells just as snow cover.
The other data set uses an index for green vegetation, that is, the fraction of absorbed photosynthetically active radiation (FPAR). The FPAR data was retrieved from the MOD15A2GFS product (Tan et al., 2009) with 1 km resolution and every 8 days and then re-gridded and interpolated to a daily 12 km-resolution data set. This approach and the resulting data set have been previously used in Weather Research and Forecasting (WRF)-CMAQ simulations by Ran et al. (2015) and . The spatial distribution and seasonal variation of the vegetation fractional coverage from these two datasets are discussed in Section 3.1.

Vegetation Height Data
The CMAQ dust model contains a landuse-specific look-up table for vegetation heights  which reflects the growth-decay cycle of general vegetation and was adopted as the height for PV (h PV ) in this study (Table 1). Since the phenological and geometric characteristics of dead plant and litter are different from those of PV, we, for the first time, defined the new parameter of NPV height and incorporated it into the dust model in CMAQ along with the implementation of NPV coverage (Table 1). The magnitude of assigned NPV heights were based on the PV heights and the relations between PV and NPV heights which were inferred from simulations and measurements of biomass and coverage. We assumed that the biomass roughly scales with the product of coverage and height. Then, considering that the modeled and observed vegetation coverage of annual grasses does not change significantly during the transition of the maximum green biomass into the maximum NPV biomass (Pierre et al., 2015), we confined the maximum NPV height to be around half of the maximum height of PV over each landuse following the quantitative relation between the interannual peaks of PV and NPV biomass (Pierre et al., 2015). As for the seasonality, the assigned NPV heights peak in September or October, based on the evidence that the biomass of senescent plants in grassland peaked between September and Novem ber and that the interannual averaged senescence period for typical steppe plants was between September and October (Nandintsetseg & Shinoda, 2015;Shinoda et al., 2011).
We realize that the predefined vegetation heights are not sufficient to represent the realistic plant community and acknowledge this as a source for uncertainties. To be more specific, first, the vegetation height can vary by a factor of two to four within a hundred meters in dust source regions in North America (Sankey et al., 2015;Weltz et al., 1994), so it is challenging to accurately represent the height with the current model resolution. Second, while the vegetation community is very diverse in the North American deserts (Warner, 2004), the magnitude of our NPV heights are based on biomass information of merely annual grass. Ideally, additional types of vegetation (e.g., shrubs) should be considered to represent the heterogeneity of vegetation in North America. Third, the values for PV heights were developed for Asian deserts (Xi & Sokolik, 2015) and the NPV heights were estimated using information about African and Asian vegetation (Nandintsetseg & Shinoda, 2015;Pierre et al., 2015), which would introduce uncertainties to our study on the Southwestern US. For example, the vegetation in North America and in central Asia show different seasonal trends (compare Figure 3 and Figure S1 in Supporting Information S1).
To better understand and quantify the effects of these uncertainties in vegetation heights, particularly those associated with the seasonality discrepancies between African and American plant species, we tested the model's sensitivity to both PV and NPV heights by conducting two sets of sensitivity runs over selected time span in different seasons. Details about the design of these sensitivity tests are described in Section 2.4.

Parameterization of Vegetation in the Dust Model
The windblown dust scheme used in CMAQ is a physics-based model described in detail by . Here, we present a short overview of the scheme and focus on the representation of vegetation in the model. Saltation bombardment is deemed as the main mechanism of aeolian dust emissions. The module calculates the bulk vertical dust emission and assigns the total mass into fine and coarse modes. The total mass of vertical dust flux is determined based on a total horizontal flux and a vertical-to-horizontal flux ratio (Equation S1 in Supporting Information S1). The ratio is dependent on soil properties and scales with the friction velocity (Lu & Shao, 1999 Note. All units are in cm. in Supporting Information S1). The size-resolved horizontal flux from each landuse, F H (L,D) is determined using the following expression: where L annotates the erodible landuse type, and D is the particle diameter. C is a constant of proportionality set to 1.0, ρ a is the air density, g is the gravitational acceleration, u * (L) is the friction velocity, and u *,t (L,D) is the threshold friction velocity.
The threshold friction velocity governs the initiation of saltation. It is modeled as an ideal threshold friction velocity corrected with two factors for soil moisture and roughness elements. * , ( , ) = * , 0( ) ( ) Here, u *,t0 (D) is the ideal threshold velocity for particles with diameter D on dry and smooth surfaces. The ideal threshold velocity, u *,t0 (D), is calculated based on Equation S3 in Supporting Information S1 (Shao & Lu, 2000). The f m and f r (L) are correction factors for soil moisture and surface roughness, respectively, both of which are equal or greater than 1.0. The soil moisture factor is determined according to a model by Fécan et al. (1998).
The roughness factor is determined using a double drag partitioning concept to take both the solid elements and the vegetation into account (Darmenova et al., 2009;Raupach et al., 1993).
Here, σ V and σ S are the basal-to-frontal area ratios of vegetation and solid elements, β V and β S are ratios of drag coefficients on vegetation and solid elements to the drag coefficient on bare surface, m V and m S account for the differences between average surface stress and maximum surface stress. The values for above parameters are from Darmenova et al. (2009), that is, σ V = 1.45, σ S = 1.0, β V = 202, β S = 90, m V = 0.16, and m S = 0.5. A V is the fractional coverage of total vegetation, and λ V and λ S (L) are surface roughness densities of vegetation and solid elements. Values of λ S (L) for each landuse type were adopted from Xi and Sokolik (2015) and Darmenova et al. (2009).
The λ V is calculated from vegetation coverage, A V following a relation proposed by Shao et al. (1996).
where κ = 0.4 is the von Kármán constant, U 10 is the 10-m wind speed, and z 0 is the surface roughness length relevant to dust emission processes (cf. the aerodynamic roughness used by WRF). Equation 5 uses the Monin-Obukhov similarity theory and assumes a neutrally stable boundary layer. The error from this assumption is likely to be insignificant during dust storms, which are usually associated with strong winds (Darmenova et al., 2009). Comparisons of the original and corrected friction velocities are shown in Figure S2 in Supporting Information S1. In general, the corrected friction velocity used in the dust model is lower than the mesoscale friction velocity calculated in WRF.
The aeolian roughness length (z 0 in Equation 5) scales with the physical height of roughness elements on the surface. To determine z 0 , we adopted the empirical relation developed by  10.1029/2021JD035243 6 of 22 where λ is the total roughness density, and it is defined as the sum of roughness density for solid elements and total vegetation (λ = λ S +λ V ). The h is the total effective height of roughness elements. In this study, the effective heights of roughness components were updated with the inclusion of NPV. It was calculated as the weighted average of roughness heights based on the roughness density: In general, it can be seen that the vegetation coverage and heights are important components of the dust scheme. The vegetation coverage not only determines the fractions of available surface to wind erosion, but also alters the friction velocity u * (via surface roughness length, Equations 5-6) and its threshold value u *,t (via surface roughness factor, Equation 3). The vegetation heights contribute to roughness height and modifies the roughness length and thus the friction velocity u * (Equations 5-7).

Sensitivity Tests on Vegetation Heights
As described above, the impacts of NPV coverage on dust emissions can be compounded by the uncertainties associated with vegetation heights. Therefore, we designed two sets of sensitivity tests to quantify the influence of uncertainties in PV and NPV heights on dust simulations. We particularly considered the uncertainty related to varied seasonality of vegetation heights among different continents as described in Section 2.2 and adjusted vegetation heights in the sensitivity runs so that their annual trends are more representative for North American plants. In order to make these sensitivity tests more physics-based and systematic, we determined the perturbation in vegetation height based on the average vegetation coverage over the western US as observed by the satellite (Figure 3). Specifically, the averaged NPV coverage gradually grows from September and stays at the highest throughout winter and early spring, so we accordingly lowered the NPV height in September, and raised it in late January to early February and March in our sensitivity tests. Similarly, for the PV height, we raised its value in July and September and lowered it in March to represent a shifted peak from spring toward summer based on the trend of PV coverage over the western US. The levels of variation were set to ± 50% of the predefined PV and NPV heights in Table 1. Simulations were run for at least 15 days for each selected season with spin-up time for no less than 6 days. Results of the sensitivity simulations will be presented later in Section 3.5.

Model Setup
The CMAQ model version 5.3 (Appel et al., 2021) was used in this study. The model domain consisting of the CONUS, as well as parts of Mexico and Canada ( Figure 1) was discretized using a 12-km horizontal grid and 35 vertical layers. Simulations were performed for the entire year 2016 with a clean initial condition and 10-day spin-up time. The meteorological inputs to CMAQ were generated by a WRF version 3.8 simulation and then processed with the Meteorology-Chemistry Interface Processor (MCIP) version 4.3. Anthropogenic emissions input data were provided by the emissions modeling platform run by United States Environmental Protection Agency and biogenic emissions were calculated in-line. The boundary conditions were derived from hemispheric simulations of CMAQ version 5.3. The Biogenic Emission Landcover Database version 3 (BELD3) was used in the dust scheme and three landuse types were considered as erodible lands, namely USGS_shrubland, USGS_shrubgrass, and USGS_sprsbarren. The total fractions of these three types of erodible lands are shown in Figure 1. The soil type information was based on US State Soil Geographic (STATSGO) soil database (Miller & White, 1998) and four soil textures (clay, silt, fine-to-medium sand, and coarse sand) were identified for each soil type following Tegen et al. (2002). Ammonia bi-directional flux and an updated M3dry model were used for deposition. The deposition processes considered include turbulent dry deposition, gravitational settling, Brownian diffusion, surface impaction, and wet deposition by cloud (Byun & Ching, 1999;Pleim & Ran, 2011). CB06r3 chemical mechanism and AERO7 aerosol model were used for atmospheric chemistry. Details on all other settings as well as the model evaluation can be found in Appel et al. (2021).

Evaluation Methodology
The windblown dust emissions in CMAQ are confined to erodible lands, so the evaluation of dust simulations was focused on the western states. We selected six states of Wyoming, Nevada, Utah, Colorado, Arizona, and New Mexico, which cover most of the active dust source regions in the US, to perform model evaluations. Dust events are essential sources for minerals in the air. Therefore, we used "fine soil" (hereafter simply soil) concentrations as defined by the Interagency Monitoring of Protected Visual Environments (IMPROVE) sites (http://vista.cira. colostate.edu/Improve/) to represent dust emissions and used the IMPROVE observations to evaluate the simulations. Modeling results suggest that the fine soil only makes up less than 40% of total fine particulate matter (PM 2.5 ) ( Figure S3 in Supporting Information S1), so the particulate matter (PM) concentrations are not as suitable as soil concentration for evaluating the model performance on dust simulations. The IMPROVE network was designated to monitor the visibility in national parks and its sites concentrate in the western US. The observatory data were available throughout 2016 every three days. The outputs from the CMAQ model were post-processed and the soil concentrations were calculated following the equation for soil adopted by the IMPROVE sites.
This equation considered the chemical composition of the oxides for predominant elements in soil (Malm, 1994).
The mean bias (MB), the normalized mean bias (NMB), the mean error (ME), the normalized mean error (NME), and the Pearson correlation coefficient between simulations and observations were used to access the simulated results. These statistics were computed based on all pairs of observations and simulations from all IMPROVE sites in selected states during the selected time span, so both the temporal and spatial variations were accounted for.

Vegetation Coverage
We present snapshots of the MODIS FPAR, the fractions of PV, NPV, and total vegetation in the middle of each season, as well as the ratio of NPV to total vegetation fractions in Figure 2. Clear contrasts in coverage and seasonal variation between the eastern and the western parts of the region are seen in maps for green and brown vegetation (Figures 2a-2c), which are most apparent in summer. Therefore, we divided the study domain into the eastern and the western parts to conduct quantitative analysis. The separation was chosen to be along 97°W with consideration of both our vegetation maps and the shifted 100th meridian. The 100th meridian is a historical divide between the humid eastern and the arid western America, and it was found to shift eastwards due to climate change and human activities over past centuries . Erodible lands, which are identified based on Biogenic Emission Landcover Database version 3 (BELD3) (Figure 1) are basically located in the resulting western part. We present the average vegetation fractions from two datasets over the east and the west throughout 2016 in Figure 3.
The vegetation fractions over Australia derived using the same SMA method were evaluated, and the root-meansquare errors for PV, NPV, and bare soil fractions were 0.13, 0.18, and 0.16, respectively (Guerschman et al., 2015). the Mogollon Rim in the western US, and they decrease northwestwards from the southeast coast of the US. The averaged MODIS FPAR and SMA derived PV fractions differ by less than 9% in the east and less than 2% in the west throughout 2016. The comparable values of PV fractions and MODIS FPAR demonstrate that the SMA method is generally good at resolving green vegetation over the study domain.
The NPV coverage generally has reversed spatial and temporal patterns relative to the green vegetation coverage. The average NPV fraction over the western domain ranges between 26% and 40%, which is at the minimum in late-June to mid-July, starts to increase in August, and reaches a maximum between December and February (Figure 3a). These seasonal trends resemble the accumulation and decay of senescent plant materials. Note that, the maximum average coverage of NPV (40%) in 2016 exceeds that of PV (31%) in the west. This suggests that there may be other sources of NPV in the winter in addition to withered green vegetation from the same year, probably perennial dead biomass. This might also suggest some overestimations of NPV coverage using the SMA method. In the southwestern US, the average NPV fractions are high over most arid or semi-arid areas including the Chihuahuan Desert, the Colorado Plateau, and the Great Basin. The NPV coverage varies around 40%-50% across the year in these areas. In forested areas in the southwest where PV fractions are relatively high, the NPV coverage is relatively low and varies around 30% throughout 2016.
Maps for the ratios of NPV to total vegetation coverage or the relative NPV ratios (Figure 2e) highlight areas where NPV is the dominant component of vegetation. The relative NPV ratios have similar spatial distributions as the NPV fractions, but with stronger contrasts. In the western US, the regions with high relative NPV ratios (nearly 100% in winter) greatly overlap with those with the erodible landuse ( Figure 1). Most of these areas also have high NPV coverage (>40%in winter), except for the Mojave Desert in the southmost corner of California where the NPV fractions are relatively low (at around 30% in winter). The relative ratios of NPV are constantly 100% in some southwestern desert lands (e.g., southern Nevada and eastern California) all year round due to no detection of PV.
The total vegetation coverage seems quite stable over the year (Figure 2d), as a result of the opposite spatial and seasonal trend of PV and NPV coverage. The average total vegetation fraction varies slightly around 56% in the west over the year. As expected, the coverage is relatively low in desert lands in the southwest. The annual average of total coverage is around 30% in the Great Basin, the Red Desert, the Sonoran Desert, the Colorado Plateau, and the Chihuahuan Desert, and it gradually increases toward the edges of these drylands, reaching around 60% at the rims. The annual standard deviation of total vegetation coverage is within 10% over most of the erodible lands, and reaches 20% in the Salt Lake Desert and some areas in northern Nevada ( Figure S4 in Supporting Information S1). In the west, at places where NPV is abundant (both fractions and relative ratios), the PV is relatively sparse, and the total vegetation coverage is lower than other areas. These NPV-rich regions highly overlap with regions with erodible landuse (Figure 1). Therefore, NPV is expected to have an impact on dust emissions from these source regions.

Response of Seasonal Dust Emission Fluxes to the NPV Coverage
We compared the seasonally averaged simulation results of dust emission flux from the FPAR run and the TOTAL run in Figures  Similarly, the NPV coverage in eastern New Mexico and northern Arizona is also high (above 40% in the middle of the spring). Though the dust source regions in the TOTAL case become very limited, the maximum dust fluxes still reach very high values (>20 mg/km 2 /s). Figure S5 in Supporting Information S1 further compares the dust flux from two simulations during three spring dust events. As shown in the three cases, the areas of dust sources in northwestern Nevada, in Colorado Plateau, and in the Salt Lake desert greatly shrink after introducing NPV.

Response of Seasonal Dust Concentrations to the NPV Coverage
With dust emission fluxes being highly sporadic both in location and time, we focused on dust concentrations in our spatiotemporal averaging and seasonal analysis. We presented the seasonal averages of simulated soil concentrations from the TOTAL run and the FPAR run in Figures 5a and 5b, respectively. Comparisons of soil concentration among all seasons reveal that spring has the highest soil concentrations for both cases. In spring, simulations from the TOTAL run show that the average soil concentrations are higher (above 1.5 μg/m 3 ) in the south, including most areas in Arizona and New Mexico, southern California, and northern Mexico, than other areas. As for the FPAR simulations, soil concentrations are consistently above 1.5 μg/m 3 over the southwestern US during spring and the most pronounced soil concentrations (above 3 μg/m 3 ) are seen in several sub-regions with desert lands, including the Salt Lake Desert, the Colorado Plateau, the Sonoran Desert, and the Chihuahuan Desert. In both cases, the magnitude and spatial distribution of soil concentrations during summer and autumn are comparable, but more soils are emitted from the Great Basin in summer than in autumn from the FPAR run. The soil concentrations in winter are similar to those in autumn, but are less in the northern states, probably due to snow cover. In general, simulations from both runs are able to capture the seasonal variation of dust emissions and highlight locations with high dust emissions. Note that the discussion above does not apply to Canada because high soil concentrations in that region are unlikely caused by windblown dust emissions due to the low fractions of erodible land (Figure 1), and instead are probably sourced from agricultural activities or unpaved roads.
Then, we analyzed the differences in simulated soil concentrations between the two runs to understand how NPV modulates windblown dust emissions. Figure 5c shows the changes in soil concentrations as percentage of the simulations from the FPAR run. Spring witnesses the most dramatic changes in soil concentrations in terms of affected areas. With the total vegetation data, the seasonal averaged soil concentrations reduce the most (above 50%) in dust source regions that generate the most soils in the FPAR run ( Figure 5b). The transport of emitted dust aerosols results in above 30% of changes in soil concentrations over most areas in the southwestern US. During summer, averaged soil concentrations reduce by over 50% in the Sonoran Desert and the Salt Lake Desert, and by over 20% also in the Mojave Desert, the Great Basin in Nevada, Wyoming, and western Colorado when using the total vegetation data. Changes in soil concentrations during autumn are less significant and major changes (over 50%) mainly occur over small areas in the Salt Lake Desert, northwestern Nevada, Wyoming, and eastern New Mexico. In winter, the most pronounced suppression on soil emission due to NPV occurs in the same regions as those for spring except for the Salt Lake Desert, and in eastern Montana additionally. The reductions in soil emissions, nevertheless, influence smaller areas downwind. Areas including southern Nevada, western Utah, eastern Arizona, and western New Mexico experience less than 10% of changes in soil concentrations, likely due to the shorter distances traveled by the dust plumes.
Most areas in the southwestern US. experience 10%-70% of reductions in soil concentrations after replacing the MODIS FPAR data with the total vegetation data during all seasons but winter. The highest percentage (above 50%) of changes in soil concentrations are seen in dust source regions with high NPV fractions (>40% in winter). Regions with high percentage of changes across multiple seasons highlight places where the dust emissions are most susceptible to NPV. These regions include the Sonoran Desert in Baja California, Mexico, the Chihuahuan Desert in New Mexico, the Great Basin in northwestern Nevada, and the Colorado Plateau in southeastern Utah, where the seasonal averaged soil concentrations are significantly suppressed by NPV throughout the year. Furthermore, the reductions in soil emissions induced by NPV are significant in the Salt Lake Desert from spring to autumn.
To shed light on the performances of the two model runs, we evaluated both simulations with ground observations from the IMPROVE sites. The annual and seasonal statistics of the models' performances over selected western states are presented in Table 2. During spring, the NMB of simulated soil concentrations from the TOTAL run (3.6%) is much smaller than that from the FPAR run (72.8%). The correlation coefficient between simulations and observations increases from 0.44 to 0.52 when replacing MODIS FPAR with total vegetation. Similar results are seen in winter, when the NMB decreases by more than a half with total vegetation. In general, including NPV in the model attenuates the overestimations of dust emissions during spring and winter. As discussed above, the effect of NPV on reducing dust overpredictions is more pronounced in spring than in winter, even though the NPV fractions are the highest in winter.
Both simulations, however, underpredict dust emissions during summer and the accuracy of simulations from the TOTAL run drop compared to those from the FPAR run, with the mean error increasing by 0.05 μg/m 3 and the mean bias increasing by around 0.1 μg/m 3 . The underpredictions in summer are likely due to the inability of the model to capture small-scale convective storms, as discussed by several previous studies (Anisimov et

Table 2 Seasonal Statistics for Simulated Soil Concentrations From the Fraction of Absorbed Photosynthetically Active Radiation Run and the TOTAL Run
assimilation and sub-grid wind distribution in the CMAQ dust model to simulate convective storms, but these modifications were not included in this study. The effects of NPV on dust emissions after the convective storms are included and the systematic negative biases are addressed need to be further analyzed. In autumn, the implementation of NPV slightly reduces the NME by 1.3%, but increases the magnitude of NMB by 9.8% and reduces the correlation coefficient by 0.06.
Overall, introducing NPV in these simulations changes the overprediction of soil concentration to underprediction, with the annual simulations from the TOTAL run underpredicting the soil concentration by 0.2 μg/m 3 and those from the FPAR run overpredicting by 0.05 μg/m 3 . The ME for the TOTAL simulations decreases by 0.1 μg/m 3 while the absolute value of MB for the TOTAL simulation increases compared to the FPAR simulations, indicating that the TOTAL simulations are less scattered but have a higher systematic bias. The correlation from the TOTAL run slightly increases by 0.03.
We further investigated the spatial distribution of differences in model biases between the two runs. The seasonal averaged changes in biases after using total vegetation data are shown in Figure 6. During spring, the averaged biases in soil concentrations decrease in the TOTAL run at 93% of the IMPROVE sites in the western study domain as defined in Section 2.1, which are reductions in overpredictions as informed by statistics from  −/+ 0.2 μg/m 3 with no apparent spatial characteristics. Overall, comparisons of annual biases from the two runs ( Figure S6 in Supporting Information S1) show that the annual averages are improved by over 0.1 μg/m 3 in New Mexico, southern Utah, southern Colorado, but underpredictions are worsened in southern Arizona after implementing NPV. The spatial distribution of annual correlations at each site for the two simulations ( Figure S7 in Supporting Information S1) shows that the correlation increases over the northwestern states while decreasing over the southwestern states.

Mechanisms for Modulating Dust Emissions by the NPV
Vegetation modulates windblown dust emissions through multiple pathways. In-depth comparisons between the intermediate parameters from the FPAR run and the TOTAL run simulations allow us to understand different mechanisms by which NPV impacts windblown dust emissions.
First, the vegetation coverage directly controls the fraction of land susceptible to wind erosion. In the model, the total flux of dust in each grid cell is calculated as the weighted sum of fluxes from three erodible landuse (see Figure 1) multiplied by the fraction of land that is not covered with vegetation (Equation S2 in Supporting Information S1). The mechanism considered here is that vegetation covering the bare soil prevents the saltating particles from impacting the surface and smaller dust particles from being ejected into the atmosphere. Figure 7 presents the fractions of vegetation-free land available for dust emissions from the two runs. The decrease in the fractions of land susceptible to dust emission is significant in all seasons after the MODIS FPAR is replaced with the total vegetation. The reduction is the largest in autumn, followed by summer, and is comparable in spring and winter. Places with the most changes are the Salt Lake Desert, the northern Great Basin located in Nevada, Oregon, and Idaho, northeastern Wyoming, and Montana, where the vegetation-free erodible lands decrease by up to 50% in autumn. In most grid cells containing dust source regions, the fractions of vegetation-free erodible lands decrease by at least 0.2 in all seasons. Because the total emission of dust from each grid cell is proportional to the fractions of vegetation-free erodible land, these significant reductions suggest that one strong mechanism for NPV to reduce dust emission is to prevent dust particles from leaving the surface. Nevertheless, potential underestimation of dust emissions when using total vegetation is possible due to two reasons. First, the actual fractions of land available to dust emission in a grid cell may be underestimated because the vegetation on non-erodible lands should not affect the exposed areas with erodible landuse but these sub-grid details were not resolved. Higher resolutions of landuse and vegetation inputs should help to address this issue. Second, this mechanism of dust suppression is only true for vegetation that is closer to the ground. Some NPV, such as standing dead trees, are detected by the satellite but cannot prevent soil erosion through this pathway.
Second, the initiation of dust generation is jointly controlled by the friction velocity and its threshold value (see Equation 1). Emissions of windblown dust can only occur and sustain when the friction velocity exceeds the threshold velocity. Dust fluxes are eliminated from some regions (e.g., Wyoming) due to NPV (Figure 4), which may be, in theory, induced by either increased threshold velocity or decreased friction velocity, or both. To pinpoint the exact pathway for emission suppression, we investigated the differences in modeled threshold velocity between the two simulations. Vegetation can increase the threshold velocity by extracting the wind stress exerted on the ground surface through drag partitioning process . In our model, the threshold velocity was calculated as the ideal threshold velocity modified by correction factors for soil moisture and surface roughness (see Equation 2). The inclusion of NPV would not change the ideal threshold velocity or soil moisture (vegetation can affect soil moisture by intercepting rainfall and uptake soil water (Teuling, 2005), but these interactions were not implemented in our model), so the change in the correction factor for surface roughness can reflect the change in the magnitude of threshold velocity. We illustrate the seasonal averages of the correction factor for roughness calculated as the average of correction factors for three erodible landuse types (f r (L)) weighted by the fractions of each type (A L ) in Figure 8. In practice, the model actually used the roughness correction factor separately over each landuse type, and the values shown here were not used directly in the model. However, they serve as a good representation of the magnitude of overall roughness correction factor.
Seasonal averaged roughness factor from the FPAR run is below 3 in most places across the year, whereas for the TOTAL run, roughness factor ranges between 2 and 12 from spring to autumn and between 1 and 7 in winter over most erodible lands. The significant changes in roughness factor induced by NPV can be explained by Equations 3 and 4, where a higher vegetation coverage in the case of total vegetation results in a higher roughness factor. Largest increases, over 100%, in seasonal averaged roughness factor are seen in Montana, Idaho, northeastern Wyoming, the northern Great Basin, eastern Colorado, and eastern New Mexico from spring to autumn. These areas generally align with regions where the contribution of NPV to total vegetation is higher as shown in Figure 2. In general, the NPV induced increases in roughness factor would raise the threshold velocity by over 50% in most areas and over 200% in certain areas such as northeastern Wyoming, and hence lessen the potential for dust emissions in these regions.
Third, vegetation can impact the friction velocity. Friction velocity is not only associated with the initiation of dust events, but also essential in determining the dust flux. The horizontal flux of dust is proportional to the third power of the friction velocity (Equation 1) and is thus very sensitive to its value. Vegetation elements that protrude the ground can alter the flow pattern and thus change the friction velocity exerted on the surface. The effect of vegetation on wind profile is accounted for using roughness length, which has a non-linear relation with the roughness density or fractional cover of vegetation as well as the vegetation height (Equations 5 and 6).
We present the seasonal averages of the weighted mean friction velocity in Figure 9, which were calculated using the same method as for the roughness factor. In general, friction velocities make non-monotonic responses to NPV. During spring and summer, the friction velocity increases by less than 0.02 m/s (below 10%) in most places, but by 20%-30% in the Salt Lake Desert, and decreases in northeastern Wyoming and southeastern Colorado with the inclusion of NPV. In autumn and winter, the increases in friction velocity due to NPV are more significant, reaching 0.03-0.05 m/s (10%-20%) on most erodible lands. Considering that the friction velocity increases in most of the erodible lands with inclusion of NPV, which would facilitate the initiation of saltation, we can conclude that the observed preventions of dust events are mainly attributed to the increases in threshold velocity instead. However, a sensitivity analysis (Darmenova et al., 2009) showed that the percentage of increase in friction velocity could lead to up to two orders of magnitude increase in the dust flux when the wind velocity is low (less than 0.05 m/s), which is our case. As a result, the increases in friction velocity caused by NPV could potentially amplify the dust flux by significant amounts.

Sensitivity of Seasonal Dust Emissions to Vegetation Heights
The responses of simulated dust emissions to the adjustment of PV and NPV heights (as described in Section 2.4) are presented here. Figure 10 shows the percentage changes in soil concentration from the six sensitivity runs as compared to the TOTAL run. Perturbing NPV heights based on the seasonal trend of NPV coverage leads to more dust emissions in winter and spring, and less emissions in autumn. The changes in NPV heights pose very slight changes on soil concentration over the majority of the southwestern US. (within ± 2%), but the local effect over certain areas are stronger. In spring, around 20% of increases in soil concentrations at the centers of hotspots in northwestern Nevada and northern Mexico leads to elevated soil concentration by up to 6% downwind over the Salt Lake. The most sensitive regions in autumn are the same as those in spring. Up to 30% of decreases in soil concentrations in these hotspots cause around 5%-10% of drop in soil concentration near the lower eastern border of California. In winter, up to 30% of increases in soil concentrations occur in hotspots in the southeastern corner of California and northern Mexico, and the soil concentrations along the western and southern borders of Arizona increase by 5%-10%. As for sensitivity tests on PV heights, 2%-6% of increases in soil concentrations occur over most of Arizona, southern New Mexico, and the northwestern Utah, after increasing PV heights by 50% in summer. In addition, a hotspot with elevated average soil concentration by around 20% inside Utah leads to a shallow plume toward the northeast direction. The changes in spring and autumn are negligible. In all the sensitivity runs, increased vegetation heights, while keeping the vegetation coverage fixed, lead to higher soil concentration because of increased friction velocity ( Figure S8 in Supporting Information S1) resulted from higher roughness length (Equation 5), and vice versa. This can be physically explained by that taller NPV (i.e., higher roughness Figure 10. The ratios of simulated soil concentration from six sensitivity runs to that from the TOTAL run. The first row presents results from three sensitivity runs where heights of non-photosynthetic vegetation (NPV) were raised by 50% in spring (16-30 March) and winter (27 January-10 February), and reduced by 50% in autumn (16-30 September) in the model. The second row presents results from the other three sensitivity runs where heights of photosynthetic vegetation (PV) were reduced by 50% in spring (15-31 March) and raised by 50% in summer (1-14 July) and autumn (14-28 September). The adjusted PV and NPV heights are intended to reflect the seasonality of the average coverage of the corresponding vegetation components over the western US.
length) induces stronger turbulence near the ground surface. Overall, 50% of changes in NPV or PV heights can result in up to 30% of changes in soil concentration over limited areas and up to about 10% of changes regionally.
Further quantitative analysis on the sensitivity of model performance to the changing NPV heights is performed based on statistics of dust simulations from the test runs and the original TOTAL run (Table 3). In spring, NMB is slightly improved from −14.5% to −13.7% and NME is improved from 50.5% to 50.1%. In autumn, the NMB and NME get worse by 1.6% and 0.5%, respectively. Improvements in NMB (by 0.3%) and NME (by 0.2%) are very mild during winter. The Pearson correlation coefficient is slightly improved in spring and winter and is almost constant in autumn with the adjustment in NPV heights. Overall, 50% of change in NPV heights leads to within ±5% of change in MB and ME in all three seasons and does not flip the pattern of over/underpredictions.

Conclusions
In this study, we incorporated a total vegetation coverage data set, which includes both photosynthetic (PV) and non-photosynthetic vegetation (NPV), into the dust module in the CMAQ model version 5.3. The objective was to investigate the responses of simulated dust emissions to the NPV coverage, as well as to the changes in vegetation heights. We conducted a pair of simulations (namely the TOTAL run and the FPAR run) using the total vegetation (PV + NPV) data set and a PV data set represented by the MODIS FPAR data over the conterminous United States for the entire year 2016. The fractional coverage of total vegetation was derived from the MODIS surface reflectance data using a spectral mixture analysis (SMA) method. A qualitative evaluation of the total vegetation fractional cover data set showed that the average PV fractions derived from the SMA approach and the average MODIS FPAR over western United States agree very well. The NPV coverage is around 40%-50% over most of the arid and semi-arid areas in the southwestern US throughout the year, which likely contains perennial dry biomass. The areas with high NPV-to-total-vegetation ratios highly overlap with the areas covered by erodible landuse, suggesting that the consideration of NPV is important for dust emissions from these source regions.
We first compared dust emission fluxes from two model runs for analyzing the effects of NPV on dust emissions. Seasonal averaged dust fluxes are significantly suppressed over most of erodible lands but the maximum dust fluxes can still be high with the inclusion of NPV. We then focused on dust concentration and analyzed the averaged percentage of differences in simulated soil concentrations between the two simulations for all seasons. Simulated soil concentrations decrease by 10%-70% due to NPV in most areas in southwestern US from spring to autumn, and these affected areas are more confined in winter. NPV induced reductions in soil concentrations are most significant in spring. Note. The simulations are compared with soil concentration observations from the Interagency Monitoring of Protected Visual Environments network. The area of concern includes states of Wyoming, Nevada, Utah, Colorado, Arizona, and New Mexico. Observation, simulation, mean bias, and mean error are in μg/m 3 , and normalized mean bias and normalized mean error are in percent. (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30),  attenuates the overpredictions at 93% of the sites in the western study domain during spring, except for those near southern Arizona. In summer, however, the underpredictions in soil concentrations are accentuated, especially in Utah, Colorado, and Wyoming. This suggested that addressing other issues of the model, such as implementing convective storm simulation, could be more effective in enhancing model performance during summer. Overall, the normalized mean error and the correlation of annual simulations are slightly improved (reduces from 79% to 67% and increases from 0.32 to 0.35, respectively) which are mainly due to the reduced overpredictions in spring.

Table 3 Statistics of Soil Concentration Simulations From the Three Sensitivity Tests on Non-Photosynthetic Vegetation Heights and the TOTAL Run During Spring
Analyses of changes in intermediate parameters in the dust model after the inclusion of NPV provide insights into the mechanisms by which NPV modulates dust emissions. The fraction of land susceptible to wind erosion is reduced by 20% in most grid cells due to NPV, indicating that NPV effectively prevents dust particles from being ejected from the ground by covering the land surface. The prevention of several dust events resulted from NPV is associated with the increases in threshold velocity and, in limited places, the decreases in friction velocity. The friction velocity experiences bi-directional changes. But on most erodible lands, it is increased by up to 10% and 10%-20% from spring to summer and from autumn to winter, respectively, which could potentially amplify the dust flux.
We then focused on the vegetation heights, which directly affect the calculation of roughness length in our model. We recognized several sources of uncertainties in the predefined NPV and PV heights and conducted sensitivity analysis on modeled dust emissions to the variation of vegetation heights. Specifically, we adjusted the original heights of NPV and PV by ± 50% in three seasons (spring, autumn, and winter for NPV; spring, summer, and autumn for PV) so that their seasonal dynamics were more consistent with their average coverage over the western US. Test simulations were conducted for at least 15 days in each of the six cases. The perturbation of NPV heights poses insignificant changes to modeled soil concentration over the majority of southwestern US, but has stronger local effects over certain areas. The hotspots with biggest changes in dust emission during spring and autumn locate in northwestern Nevada and near the southern border of California, where soil concentrations rise by about 20% and drop by nearly 30%, respectively. In winter, hotspots concentrate in southeastern corner of California and experience nearly 30% of decrease in soil concentration. Results of the sensitivity tests on PV heights show negligible changes in average soil concentration during spring and winter. In summer, a hotspot with around 20% of elevated soil concentration is seen in Utah, and regional mild increases (2%-6%) in soil concentration cover most of Arizona, southern New Mexico, and northwestern Utah. Overall, 50% of changes in vegetation heights lead to up to 30% of changes in soil concentrations over small hotspots areas, and up to around 10% of changes in the nearby downwind regions.
This paper showed that dust emissions from a large portion of erodible lands in the southwestern US are modulated by the NPV. The improvement on dust model performance with the inclusion of NPV was most prominent in spring but was not consistent across all seasons. Therefore, replacing the fractional coverage data of green vegetation currently used in many dust models with the total vegetation coverage data derived from satellite-based surface reflectance, while carefully adjusting vegetation height and seasonality, is a promising step toward better representation of land-atmosphere interactions in atmospheric and air quality models. More evaluation on the uncertainties associated with the SMA-derived total vegetation fractions will magnify the benefits of including NPV into windblown dust models. Another aspect of the representation of NPV in dust models relates to the vegetation heights. The existing implementation of vegetation heights in dust models are simplistic and subjected to large uncertainties. The moderate sensitivity of dust emission to vegetation heights found in this work suggested that incorporating a more realistic data set of vegetation heights could improve the estimation of surface roughness and thus dust emissions in future studies.
coverage of PV and NPV is described and distributed by Guerschman (2014) and alternative formats (i.e., geotiff and png files) of the coverage data can be accessed at https://eo-data.csiro.au/remotesensing/v310/global_5km/. CMAQ output data from the TOTAL run, the FPAR run, and the six sensitivity runs used in the analysis are available at the Virginia Tech research repository (Huang & Foroutan, 2022). Source codes and model input data for CMAQ version 5.3 are provided by the United States Environmental Protection Agency (2019a, 2019b) and are freely available to the public.