Comparison of Aerosol Optical Depth from MODIS Product Collection 6.1 and AERONET in the Western United States

: Recent observations reveal that dust storms are increasing in the western USA, posing imminent risks to public health, safety, and the economy. Much of the observational evidence has been obtained from ground-based platforms and the visual interpretation of satellite imagery from limited regions. Comprehensive satellite-based observations of long-term aerosol records are still lacking. In an effort to develop such a satellite aerosol dataset, we compared and evaluated the Aerosol Optical Depth (AOD) from Deep Blue (DB) and Dark Target (DT) product collection 6.1 with the Aerosol Robotic Network (AERONET) program in the western USA. We examined the seasonal and monthly average number of Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua DB AOD retrievals per 0.1 ◦ × 0.1 ◦ from January 2003 to December 2017 across the region’s different topographic, climatic, and land cover conditions. The number of retrievals in the southwest United States was on average greater than 37 days per 90 days for all seasons except summer. Springtime saw the highest number of AOD retrievals across the southwest, consistent with the peak season for synoptic-scale dust events. The majority of Arizona, New Mexico, and western Texas showed the lowest number of retrievals during the monsoon season. The majority of collocating domains of AOD from the Aqua sensor showed a better correlation with AERONET AOD than AOD from Terra, and the correlation coefﬁcients exhibited large regional variability across the study area. The correlation coefﬁcient between the couplings Aqua DB AOD-AERONET AOD and Terra DB AODAERONET AOD ranges from 0.1 to 0.94 and 0.001 to 0.94, respectively. In the majority of the sites that exhibited less than a 0.6 correlation coefﬁcient and few matched data points at the nearest single pixel, the correlations gradually improved when the spatial domain increased to a 50 km × 50 km box averaging domain. In general, the majority of the stations revealed signiﬁcant correlation between MODIS and AERONET AOD at all spatiotemporal aggregating domains, although MODIS generally overestimated AOD compared to AERONET. However, the correlation coefﬁcient in the southwest United States was the lowest and in stations from a higher latitude was the highest. The difference in the brightness of the land surface and the latitudinal differences in the aerosol inputs from the forest ﬁres and solar zenith angles are some of the factors that manifested the latitudinal correlation differences. recommend the application of medium spatiotemporal collocating resolutions (temporally ± 15 min and spatially 30 km × 30 km) for the estimation of ground-based air quality parameters. Incorporating aerosol products from geostationary satellites like NASA’s Geostationary Operational Environmental-16 (GOES-16) Advanced Baseline Imager (ABI) into the dataset used in this study can signiﬁcantly improve the characterization of the aerosol loading and air quality assessment rather than using MODIS- and AERONET-derived AOD only. Aerosol products from GOES-16 ABI possess a higher temporal resolution than MODIS and they can mitigate the sparseness of aerosol observations from MODIS-AERONET couplings. More investigations are needed to understand the aerosol impacts of long-term land cover and land use changes around AERONET sites as these surface types are very sensitive to climate change and variability, and this parameter in turn is crucial in determining the response of the land-surface-to-satellite signals.


Introduction
Aerosols affect the climate system through direct radiative forcing [1,2], indirect radiative forcing [3], and modulation of heterogenous chemistry [4]. They adversely impact human health [5,6], the environment [7,8], and the economy [9]. Linked to its major deserts in the southwest and highly vegetated mountainous ecoprovinces, the western United States is a source of many types of aerosols including both direct and chemically transformed emissions injected from natural and anthropogenic sources (e.g., dust, biomass burning, power plants, biogenic aerosols, and mobile source emissions) [10,11]. Considerable evidence indicates that during the past few decades the western United States has been experiencing an increase in the frequency and intensity of dust events and wildfires; the dust dominated aerosols mainly originate from the deserts and plains, and wildfire aerosols from forest, woodland, and grassland ecoprovinces. Tong et al. [12] reported a rapid intensification of dust storm activity over American deserts in the time range 1988-2011. Brahney et al. [13] revealed that the American West is becoming increasingly dusty based on the Ca +2 composition of precipitation from 1994 to 2009. Using iron as a proxy of mineral dust, Hand et al. [14] showed an increasing trend of dust in large parts of the western United States with +5% per year change from 1995 to 2015 during the normally dusty month of March, an effective doubling of March dust concentration during the two decades. Wildfires have also been contributing to large amounts of aerosols into the atmosphere over the American West [15][16][17][18]. McClure & Jaffe [18] showed that since the mid 1980s, the total area of the United States that has been burned by wildfires has been increasing with fires in the northwest United States accounting for around 50-60% of that increase. Climate change and variability in different climate parameters significantly contributes to the increase in wildfires [19]. Littell et al. [16] documented that most of the vegetation zones in the western United States exhibited strong year-of-fire relationships with low precipitation, intense drought, and high temperature.
In recent decades, aerosol products from satellite remote sensing, ground-based observations, and models have been widely synthesized to monitor the temporal and spatial distribution of aerosol loading on global and regional scales [20][21][22][23][24][25]. Remote sensing imagery products from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor aboard the Aqua and Terra satellites provide retrieved data on aerosol properties at a near-daily global coverage and a wide range of spatial and spectral resolutions. However, due to the potential influences of retrieval algorithms, surface albedo, topography, shadow, and image resolutions, the satellite products should be continuously validated to minimize the uncertainties [26][27][28]. Aerosol data acquired from ground stations, more specifically ground-based remote sensing approaches such as the Aerosol Robotic Network (AERONET), are crucial resources in ground-truthing satellite products as they have a higher precision of measurement at a low spatial reach [29]. The AERONET is a direct columnar measurement acquired by sunphotometer and aggregated in time at the point location. However, the satellite Aerosol Optical Depth (AOD) is a columnar measurement retrieved through radiative-transfer models based on multiple assumptions about aerosols and surface characteristics leading to various sources of uncertainty [30,31]. Satellite AOD (including MODIS) represents the AOD concentration for a fraction of a minute when a satellite is on orbit over a region [32].
Several researchers have investigated the relationship between satellite remote-sensed and ground-based aerosol products by collating both data products to characterize aerosol loading to the atmosphere in space and time [33][34][35][36]. Bréon et al. [33] evaluated the accuracy of aerosol products retrieved from the measurements of multiple satellite remote sensing instruments by comparing with AOD from the AERONET program on a global scale. Saide et al. [35] examined the performance of the aerosol data matrix from the assimilation of AOD retrieval from Geostationary Ocean Color Imager (GOCI) and MODIS instruments by comparing against the ground measurement of AERONET AOD and PM 10 in northeast Asia. Vijayakumar et al. [36] investigated dust-dominated aerosol characteristics during a specific dust event in western India by augmenting aerosol data from MODIS and AERONET, and with meteorological data of model simulations from the European Center for Medium-Range Weather Forecasts Reanalysis-Interim (ERA-Interim). Similarly, many studies utilized satellite, model simulations, and ground-based data to improve air quality modeling and forecasting in the United States [37][38][39][40], mainly focusing on the eastern United States. For example, Li et al. [39] used AOD data processed from MODIS, Multi-angle Imaging SpectroRadiometer (MISR), Sea-viewing Wide Field-of-view Sensor (SeaWiFS), and Ozone Monitoring Instrument (OMI) image products and the results were verified against AERONET AOD measurements.
Although air quality in the western United States has improved during the past few decades [41], different studies have revealed that air pollution has increased nationally since 2016 [42]. For example, according to a report from the United States Environmental Protection Agency [43], the southwest has experienced an abrupt increase of PM 10 pollution from 2015 primarily linked to mineral aerosol loading and the northwest showed a spike in PM 2.5 pollution contributed by wildfires. In addition, Tong et al. [12] documented that in the American deserts, the frequency of windblown dust storms has increased by 240% from the 1990s to the 2000s. In an effort to develop a more comprehensive and long-term fusion of ground-based and satellite-based observations, and to facilitate the validation of satellite remote sensing imagery products, this study compares and evaluates the AOD from MODIS Deep Blue (DB) and Dark Target (DT) product collection 6.1 with AERONET AOD in the western United States bounded by the −92 • meridian in the east.

Study Area
The study area encompasses the western part of the United States bounded by the −92 • meridian in the east, including parts of 22 states (Figure 1). This region is a land of extremes. It contains several major ecological regions ( Figure 1a): Warm and cold deserts, coastal and inland forests, coastal and alpine environments, shrublands and grasslands, distributed from −105 to 4302 m elevation (Figure 1b). Mainly influenced by differences in elevation and the climatological influence of the Pacific Ocean, the climate of the western United States is very diverse in nature [44]. Due to this variation, the region is subjected to a diverse spectrum of natural and anthropogenic aerosols and impacts. The black ( Figure 1a) and blue (Figure 1b) dots show the AERONET stations over the western United States, located in a wide range of ecological regions and elevations. Depending on their proximity to different ecosystems and land uses, these stations can be exposed to different aerosol sizes and types.

AERONET
The ground-based data used for this study were acquired from the AERONET network program. AERONET is a federated instrument network of sunphotometers and data archive for aerosol characterization under the administration of the NASA Goddard Space Flight Center [45]. The AERONET program uses sunphotometers that measure direct solar irradiance and sky radiance at the Earth's surface [46]. The sunphotometers are directed at the sun (or the moon if the sensitivity is high enough) and the output delivers a direct measurement of the extinction of light between the top of the atmosphere and the Earth's surface, which is directly related to AOD. In addition, using this measurement enabled other columnar optically effective aerosol properties to be acquired including volume size distribution, complex index of refraction, Angstrom exponent, and single scattering albedo. The spatial representation of AERONET sites vary between a 200-and 500-km radius depending on the specific locations influenced by topographic and climatic conditions [47]. The measurement is acquired every 15 min at eight spectral bands (340, 380, 440, 500, 675, 870, and 1020 nm). Three versions of data are available for three data quality levels: Level 1.0 (unscreened), Level 1.5 (cloud-screened and quality controlled), and Level 2 (quality assured).
We used Level 2 data of 23 AERONET stations ( Figure 1 and Table 1) over the western United States. These sites were selected for analysis based on the availability of more than five years of data record and the geographic location west of 90 • west longitude. Since AERONET AOD data were not available at a 550-nm wavelength, equivalent AERONET AOD at 550 nm (τ 550nm ) values were determined by interpolating AERONET AODs reported at 500 and 675 nm spectral channels and Angstrom exponents at 440-675 nm (α 440−675 ) using the Angstrom exponent formula in (Equation (1)) [48,49].

MODIS
MODIS sensors aboard NASA's Aqua and Terra satellites acquire near-global daily measurements using 36 spectral bands in a wide spectral range (410 to 1500 nm) [30]. MODIS on the Aqua platform was launched in December 1999 and on the Terra platform in late May 2002 [50]. Retrieved from the calibrated radiance/reflectance measurements, the MODIS Level 2 aerosol products (MYD04 from Aqua and MOD04 from Terra) monitor the ambient columnar aerosol loading and other aerosol properties over cloud-free and snow/ice-free land and ocean surfaces. Aerosol products for land are retrieved using Deep Blue (DB) and Dark Target (DT) algorithms at a nominal spatial resolution of 10 km × 10 km at nadir. These algorithms synthesize different assumptions about the Earth's surface and the expected aerosol types above these surfaces [51]. The standard aerosol product from Dark Target method works best over dark vegetated regions and does not work well over bright land surfaces. Extraction of aerosol properties from satellite imagery over bright surfaces (such as prevalent in arid and semiarid lands) is a challenging issue at many wavelengths of visible light. The DB algorithm is powerful in retrieving aerosol products over bright desert surfaces because it utilizes two blue channels (412 and 470 nm) to discern the reflectance between the land surface and aerosols [52][53][54], but it can also retrieve over highly vegetated targets. In these channels, aerosols tend to be bright and land surface features dark. The DB algorithm reports AOD at 550, 412, 470, and 670 nm, Angstrom exponent at 412-470 nm, and single scattering albedo at 412, 470, and 670 nm. In the present study, we used the latest Deep Blue Level 2 MODIS aerosol products collection 6.1 from Aqua and Terra satellites, MYD04_L2 and MOD04_L2, respectively. The uncertainty of the retrieved Terra MODIS DB AOD collection 6.1 is 0.03 ± 21% and 0.03 ± 18% with respect to the ground-truth measurement of AERONET AOD (0.03 ± 0.21 × AODAERONET and 0.03 ± 0.18 × AODAERONET) for arid and vegetated path retrievals, respectively [55]. Thus, collection 6.1 has a slightly better performance than collection 6.0. The DB aerosol product in collection 6.1 is generated by an Enhanced DB algorithm [54] by improving collection 6.0. The improvements include heavy smoke detection, artefact reduction over heterogenous terrain, improved surface modeling for elevated terrain, and bug fixes and updated regional/seasonal aerosol optical models [56]. The MYD04_L2 and MOD04_L2 data products were ordered and downloaded from the Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC) website (https://ladsweb.modaps.eosdis.nasa.gov (accessed on 1 June 2018)) covering the western United States from July 2002 to June 2018. We used daily Deep Blue AOD at 550 nm and Angstrom exponent for land at a 10 km × 10 km spatial resolution with all the available quality flags. The DB AOD is represented by the data field name "Deep_Blue_Aerosol_Optical_ Depth_550_Land". The DT AOD is represented by the field name "Optical_Depth_Land_ And_Ocean" which is aerosol optical thickness at 550 nm for both ocean (best) and land (corrected). The interpolation of AERONET and MODIS AOD, the extraction of MODIS AOD for the AERONET coordinate locations, and the statistical analysis were batch processed using MATLAB and Python software packages.

Assessment of the Number of Pixels of MODIS DB AOD Retrieval
We investigated the long-term monthly and seasonal average distribution of MODIS Aqua DB AOD retrieval per the interpolation grid cell per season and month averaged from 2003 to 2017 for the western United States. This analysis can give us information about the applicability of the aerosol product in different regions, frequency of days impacted by poor air quality associated with aerosols, and impact of cloud and snow on the number of retrievals. For this analysis we used the interpolated images created from the Level 2 granules overlapping the study area. This methodology was applied to facilitate the assessment by utilizing images with regular latitude and longitude grids rather than on irregularly spaced pixels along the orbit track (Level 2) [57]. This is performed by linearly interpolating the Level 2 granules from each day using an interpolation function from MATLAB. All pixels covering the study area were arranged into one daily file. Then, linear interpolation was applied at every fixed matrix point (0.1 • × 0.1 • ) weighing the surrounding Level 2 values by the fraction of overlapping area between the two grids. As a result, a total of 5786 files were produced extending from 4 July 2002 to 22 May 2018. Aqua products were used based on the higher correlation with AERONET AOD, discerned during the research process.

Comparison of MODIS and AERONET AOD
Aerosol optical depth retrieved from MODIS DB and DT product collection 6.1 were compared with AOD data from 23 AERONET sites. We also investigated the correlation of AOD between the couplings of MODIS Aqua-AERONET and MODIS Terra-AERONET. For this analysis, we used the original (irregular grid) Level 2 product with the data field name "Deep_Blue_Aerosol_Optical_Depth_550_Land" from both Aqua (MYD04_L2) and Terra satellites (MOD04_L2). The time series of MODIS AOD were extracted for the AERONET station locations by considering different combinations of temporal and spatial domains.
Spatially, we considered three different scenarios. The MODIS AOD pixels nearest to the AERONET sites were identified using the Euclidean distance method. The determined pixels were denoted as single pixels or a centroid representing AERONET sites, by limiting the distance between the center of pixels to a radius ≤ 10 km. In addition, we examined the impact of regions surrounding the centroids by considering two different spatial averaging domains of AOD values, the nearest nine and twenty-five pixels. For the nine pixels scenario (30 km × 30 km box or 3 × 3 pixels), AOD from the identified central points and its nearest eight surrounding pixels were averaged by setting the threshold radius distance of 20 km or less. Similarly, for the 25-pixel scenario (50 km × 50 km box or 5 × 5 pixels), AOD from the identified central points and its nearest 24 pixels were averaged by limiting the radius distance to 30 km or less.
Likewise, we considered three averaging temporal domains, nearest data, 30-min and 3-h averaging domains. The closest single MODIS AOD to the AERONET AOD were identified by limiting the difference to less than or equal to 10 min, otherwise the data were not included in the correlation analysis. For the 30-min and 3-h averaging domains, all of those AERONET AOD measurements within ±15 min and ±90 min of the MODIS overpass time were aggregated. These three temporal scenarios were applied to the AERONET AOD corresponding to the MODIS AOD retrievals averaged at the nearest single pixel and within 30 km × 30 km and 50 km × 50 km boxes. Using the mean as a statistical measure of central tendency is appropriate as it saves the mode of the quality flags from each pixel within the averaging domain [58]. After running the correlation analysis between all the corresponding collocated MODIS DB AOD and AERONET AOD measurements at all spatiotemporal scenarios, the results were compared with the correlation between MODIS Dark Target (DT) AOD and AERONET AOD at similar spatiotemporal settings. Here, we used the scientific data set named "Optical_Depth_Land_And_Ocean" at 550 nm from the Dark Target (DT) algorithm with all quality assurance flags.

Assessment of the Number of Pixels of MODIS DB AOD Retrievals
The long term seasonal and monthly average numbers of MODIS DB AOD retrievals per 0.1 • × 0.1 • grid averaged from January 2003 to December 2017 are presented in Figure 2 and 3, respectively. The results revealed that the southwestern United States showed a significant number of retrievals per long-term average seasonal and monthly distribution during the same period. During the summer months (June, July, and August) the majority of Arizona, New Mexico, and western Texas, which experience a wet monsoon during that time, showed the lowest number of AOD retrievals. As the climate is very dry, the number of retrievals in the southwestern United States was on average greater than 14.8 days (≈50%) per 30 days for the majority of months and more than 37 days (≈41%) per 90 days for the three seasons other than June, July, and August. On the other hand, the states of California and Oregon exhibited the highest retrieval during summer reaching from 65 to 85 days per 90 days. During winter (December, January, and February), the northern part of the study area experienced either no or the lowest number of retrievals due to the presence of snow and cloud cover. Due to their high reflectivity, snow and clouds interfere with the AOD retrieval by making the separation of the radiance contribution from the surface and aerosols very difficult. In addition, some regions did not show any retrieval possibly due to the exclusion of pixels due to the contribution of reflectance from water bodies.

Evaluation of MODIS AOD Using AERONET AOD
We evaluated the robustness of the MODIS AOD product against the in-situ records of AERONET AOD. To validate the respective MODIS DB (land) and DT (land and ocean) AOD from Aqua and Terra satellites against a ground-based AERONET AOD, we fitted a linear regression model at all spatiotemporal scales. After running the model, we also extracted the Pearson correlation coefficient (R), p-value (P), number of samples (N), and root mean square error (RMSE) to measure the performance of the fitted models. Spatially, we considered the nearest single pixel to the AERONET station, 3 pixels by 3 pixels (corresponding to 30 km × 30 km), and 5 pixels by 5 pixels box (corresponding to 50 km × 50 km) as aggregating domains of MODIS AOD retrievals using the AERONET locations as the center of the boxes. Temporally, we used the nearest single record to the MODIS AOD, 30 min (15 min before and after) and 3 h (90 min before and after) for averaging AERONET AOD using the MODIS AOD as a center time. The results for correlation analysis between MODIS (Aqua and Terra) DB AOD and AERONET AOD at different spatiotemporal scenarios are presented in Figure 4-6, and Table A1. The findings revealed strong and statistically significant (p-value < 0.05) correlation between MODIS DB and AERONET AOD at all spatial and temporal scales in the majority of stations, except the Red Mountain Pass station which showed a very weak and insignificant correlation at limited spatiotemporal scenarios. This is mainly linked to high cloud and snow cover in the high altitude area of the Rocky Mountains ( Figure 5 and Table A1). The Pearson correlation coefficient between the couplings Aqua AOD-AERONET AOD and Terra AOD-AERONET AOD ranges from 0.1 to 0.94 and 0.001 to 0.94, respectively, with 19 stations scoring greater than 0.4 as the correlation coefficient. The Missoula station exhibited the strongest correlation ( Figure 6) and the Red Mountain Pass the weakest correlation at all averaging domains. Overall, at the majority of collocating domains AOD from Aqua sensor showed a better correlation with AERONET AOD as compared to that of AOD from Terra. Therefore, in this region, Aqua AOD in combination with ground-based AOD and fine particulate matter measurements may be a better choice for air quality assessment, unless there is a diurnal component to the AOD-PM relationship influenced by changes in the boundary layer, or changes in relative humidity, etc. For ex-ample, a combination of these data can be employed for estimating the ground-level fine particulate matter from columnar aerosol loading using different statistical approaches including machine learning and deep learning methods. The relationship between Aqua MODIS DB AOD and AERONET AOD showed a combination of decreasing and increasing correlation with the increase in spatiotemporal averaging domains (Figure 4 and 5). The stations that showed a strong correlation (greater than 0.6 except the University of Houston and Rimrock sites) and a slight increase in the number of matched samples, experienced a gradual decrease in the correlation coefficient with an increase in the spatial domain from the nearest single pixel to 50 km × 50 km aggregating domain. The University of Houston and Rimrock sites showed very similar concentrations in the immediate and up-to-50-km radius region. However, in the majority of the sites that exhibited less than a 0.6 correlation coefficient and very few matched data points at the nearest single pixel, the correlations gradually improved when the spatial domain increased to the 50 km × 50 km box averaging domain ( Figure 5). The correlation analysis between Terra DB AOD and AERONET AOD showed a similar pattern to the Aqua DB AOD-AERONET couplings at all stations, except at CART Site, Monterey, Sevilleta, Sioux Falls, and Table Mountain stations that reversed the correlation trend by changing from the lowest to highest spatial domains. This reversal could be mainly due to the changes on the availability of matched AOD measurements linked to the overpass time of Aqua and Terra satellites, as the overpass time affects the number of matched samples. When the temporal collocation of AERONET AOD to the nearest ±15 min and ±90 min of Aqua and Terra MODIS overpass time was performed, the relationship was changed in terms of correlation coefficient, p value, and number of samples ( Figure 5 and Table A1). From this analysis, the correlation between Aqua MODIS DB AOD from the nearest pixel and AERONET AOD increased at 14 stations (representing a wide range of ecoregions), when the temporal scenario of AERONET AOD increased from the nearest single pixel to ±15 min ( Figure 5). A total of 11 of these 14 stations then exhibited a drop in correlation when the temporal aggregating domain increased from ±15 to ±90 min. However, the other nine stations showed a decreasing correlation trend from the nearest single pixel to ±15 min averaging extent ( Figure 5). Five out of the nine stations then further showed a decrease in correlation when the AERONET AOD averaged increased to ±90 min of Aqua MODIS overpass time. During the aforementioned series of temporal changes, the number of samples also changed significantly (Table A1). The number of samples increased with an increase in temporal scales. However, the rate of change was different from site to site. The rate of change from the nearest single pixel to ±15 min ranges from 4% to 11.4% (only the White Sands HELSTF site showed a change of 115%) and from ±15 to ±90 min ranges from 9.3% to 34%. In the scenario from the nearest single pixel to ±15 min, the lowest and highest percentage of changes happened at Goldstone station and Frenchman station, respectively. Increasing from ±15 to ±90 min, the lowest and highest percentage of change in the number of samples was shown at Railroad Valley station and White Sands HELSTF, respectively. These changes could be partly contributed by the difference in the availability of AERONET AOD records at different temporal averaging domains. The change in the number of samples either weakened or strengthened the correlation between the variables. It is worth noting that the gaps in AOD measurements between the AERONET sites might impact the correlation spectrum. Consequently, it is difficult to conclude how the rate of change in the number of matched samples influenced the relationship, as it may depend on the geographical location, aerosol type, meteorology, etc. In addition, the increase in the temporal scale caused the disappearance of the only non-significant correlation (Table A1). Perhaps due to the difference in satellite overpass time and associated differences in the matched number of samples with AERONET AOD measurements, Terra DB AOD from the nearest single pixel (spatially) showed a significant difference with Aqua DB AOD when the temporal domain increased from the nearest single pixel to ±15 min and then to ±90 min. Twelve stations (a different set of stations as compared to Aqua) experienced an increase in the Pearson correlation coefficient when the temporal distance increased from the nearest single pixel to the ±15 min averaging extent. In eight of these stations, the correlation further increased as the temporal averaging domain increased from ±15 to ±90 min, while the remaining four stations showed a decrease in correlation. On the other hand, 11 of the stations exhibited a decreasing correlation when the temporal aggregating domain increased from the nearest single pixel to ±15 min. Extending the temporal averaging domain from ±15 to ±90 min led to nine out of the eleven stations to show further decreases in correlation, while only two stations bounced back to an increase in correlation.
Similar to the Aqua DB AOD-AERONET coupling, Terra DB AOD-AERONET coupling revealed an increase in the number of matched samples and change in the significance of correlation with an increase in temporal resolution. The percentage change of the number of samples extends from 3.6% to 12.9% and 7.4% to 31.7% when the temporal domain increases from the nearest single pixel to ±15 min and from ±15 min to ±90 min, respectively. From the nearest single pixel to ±15 min, the lowest and highest percentage of change was observed at Goldstone station and Frenchman station, respectively, stations representing the Mojave Desert ecoregion. The dynamics of these changes could have emerged from the impact of the bright desert surface on the MODIS AOD retrieval. Changing from the ±15 to ±90 min domain, the lowest and highest percentage of change occurred at Red Mountain Pass station and Rimrock station, respectively.
From the investigation of the Pearson correlation coefficient between Aqua DB AOD aggregated using a 30 km × 30 km box and AERONET AOD, 14 stations showed an increase in correlation when the temporal extent increased from the nearest single pixel to ±15 min ( Figure 5). These stations represent a wide range of ecoregions. This finding suggests that over these 14 stations a major fraction of the aerosols may have been in suspension for a long period of time or the aerosol sources were continually injecting aerosols into the atmosphere. The other nine stations showed a decreasing correlation. This decrease in correlation could be linked to the dispersion or settlement of aerosol particles from the atmosphere. Extending temporal domains from ±15 min to ±90 min resulted in a decline of correlation over 16 stations. This observation indicates that as the spatiotemporal averaging domains for AERONET AOD increase, the relationship between both variables starts to diminish. When the time interval increased from the nearest single pixel to ±15 min, the lowest and highest percentage change of the number of samples were observed at the Bozeman station with 4.6% and the Frenchman Flat station with 12.1%, respectively. As the temporal aggregation is associated with the AERONET AOD data, this spectrum in the rate of change could be linked to the availability of data before and after the satellite overpass time. Increasing the temporal extent from ±15 min to ±90 min, the number of samples quadrupled compromising the strength of correlation coefficients. The lowest and highest correlation was observed at the Goldstone station with 14.5% and the Rimrock station with 46.5%, respectively (Table A1).
Terra MODIS DB AOD aggregated using a 30 km × 30 km box showed an increase in correlation at 14 stations when compared with AERONET AOD when the temporal extent increased from the nearest single pixel to ±15 min. The remaining nine stations experienced a decrease in correlation. With an increase in the temporal averaging range from ±15 to ±90 min, 13 stations exhibited a decrease in correlation, while the rest of the stations showed an increase in correlation. From the nearest single pixel to ±15 min, the lowest and highest percentage of change occurred for the number of samples at the Goldstone station with 4.8% and at the Frenchman Flat station with 13.6%, respectively. Changing from ±15 min to ±90 min, the lowest and highest percentage of change observed in the number of samples at the Goldstone station with 15.9% and at the Ames station with 42.4%, respectively.
Aqua MODIS DB AOD averaged using a 50 km × 50 km box (corresponding to 5 pixels × 5 pixels) was also compared with AERONET AOD processed at the nearest single pixel, ±15 min and ±90 min to the Aqua MODIS overpass time. As in the previous scenarios, the ±15 min for aggregating AERONET AOD showed a better performance in improving the correlation coefficient. Fourteen stations displayed an increase in correlation when the aggregation of AERONET AOD increased from the nearest single pixel to ±15 min. A further increase to ±90 min led to a decrease in correlation in 13 stations and an increase in 15 stations. This finding indicated that with coarser spatiotemporal domains (50 km × 50 km box and ±90 min) the correlation strength was diminished by the AOD values acquired from distant pixels in time and space. In other words, with a further increase in the spatiotemporal scale the pixels start to decrease their representation with the center of the box. The Bozeman station (4.9%) showed the lowest percentage of change in the number of samples and the University of Houston station (13.6%) the highest when increased from the nearest single pixel to ±15 min. The heterogeneity of geographical location, land cover/land use setting, and sources of aerosols in both these stations could be attributed to the variability. However, when the temporal domain increased from ±15 min to ±90 min, the Goldstone station saw the lowest (16.8%) and the Rimrock station the highest (52%).
Furthermore, the correlation coefficient between Aqua DB AOD and AERONET AOD exhibited large regional variability across the western United States. In particular, the correlation coefficient values in the stations from the southwestern United States were relatively lower than those in other regions of the study area across all spatiotemporal scenarios, with the exception of the Red Mountain Pass station ( Figure 6). The stations in the northern part of the study area displayed the strongest correlation. The correlation between Terra DB AOD and AERONET AOD showed the same distribution pattern of correlation as the Aqua MODIS DB AOD-AERONET AOD coupling, however, slightly lower in strength. The Railroad Valley and Red Mountain stations showed the lowest correlation at all temporal and spatial domains.
The latitudinal differences in the correlation between MODIS and AERONET AOD in the western United States may be primarily a manifestation of uncertainty in the satelliteretrieved AOD from MODIS. Given its areas of bright surface reflectance where it is difficult for the satellite to separate the surface from the atmospheric signal, the desert regions of southwest North America provide challenges for accurate satellite retrievals. This would result in reduced correlations in the southwestern regions of the United States (Thomas F. Eck, personal communication, 2021). Additionally, the prevalence of large forest fires in the northwestern United States result in high AOD values, and the correlation often improves when the AOD data range is larger. This factor likely increases the correlation in the northern latitudes (Thomas F. Eck, personal communication, 2021). Furthermore, the accuracy of the AERONET-measured AOD is very high at around 0.01 (around 0.02 in ultraviolet) for the overhead sun [59,60]. This accuracy is even better at larger zenith angles, therefore there is some improvement in the accuracy of AERONET-measured AOD at higher latitudes.
The AOD measurement at 550 nm from the MODIS DT aboard Aqua and Terra satellites was evaluated against AERONET AOD. The results from this comparison were also compared with the correlation coefficient from the MODIS DB AOD-AERONET AOD coupling, with the objective to identify the product that better characterizes the atmospheric aerosol loading across the western United States. We used the AOD at 550 nm retrieved from land and ocean through the DT algorithm with the scientific data name of "Optical_Depth_Land_And_Ocean". The results for the evaluation of the MODIS DT AOD are presented in Figure 7-9 and Table A2.  Figure 7). However, in the case of Aqua DB AOD and AERONET AOD comparisons it recorded a correlation coefficient of 0.941 (Figure 4). Similar patterns of correlation differences were observed across all the remaining spatiotemporal settings. The University of Houston station, in a coastal and urban environment showed a percentage change increase of 55% going from Aqua DB AOD (0.612) to Aqua DT AOD (0.952) from the nearest single pixel in space and time. The inclusion of land and ocean pixels in the DT algorithm could have partly influenced the strength of the correlation. In addition, an increase in correlation from DB to DT in the coastal stations of Monterey and UCSB to those from central California may be partly contributed by the presence of more vegetated surface and sea salt aerosol presence that was commensurate with the merits of the DT algorithm [61]. The Maricopa site showed the lowest correlation between Aqua DT AOD acquired from the nearest pixel in space and AERONET AOD aggregated from the nearest pixel, ±15 min, and ±90 min (Figure 9). This is due to an insufficient number of matched data points (Table A2). Sayer et al. [50] suggested that the difference in correlation between DB AOD-AERONET AOD and DT AOD-AERONET AOD couplings was linked to the persistently low AOD measurements, tricky surface conditions, and algorithmic assumptions. According to their global-scale study, southwest North America was one of the few regions that showed the largest differences in seasonal mean AOD between DB AOD and DT AOD.

Discussion
In summary, this study demonstrated that the augmentation of aerosol optical depth data products from satellite remote sensing and ground-based platforms can improve the characterization of aerosol loading in space and time across the western United States. Monthly and seasonal maps exhibited a wide range of MODIS aerosol optical parameter retrievals and identified the time periods impacted by the data gaps. These maps can be used as a reference for spotting the temporal and spatial spaces during the process of location selection for the deployment of ground-based instruments to augment satellite remote sensing measurements. We evaluated the MODIS Level 2 AOD from DB and DT retrieval algorithms using a ground-based remote sensing product of AOD from AERONET by utilizing different spatiotemporal aggregation domains using the coordinate information of the 23 AERONET sites as a centroid for aggregation. The findings from this analysis can help in the calibration of MODIS aerosol retrieval algorithms in the locations where their correlation with AERONET AOD were very low, such as the southwestern United States. In addition, it can provide information about the performance of AERONET sites.
Although MODIS Level 2 (10 km × 10 km) aerosol products are available at a better spatial resolution than Level MODIS Level 3 (1 • × 1 • ) almost all the studies similar to the current study focusing in the United States utilized used MODIS Level 3 gridded aerosol data. This is because the quantitative use of Level 2 aerosol data requires a significant level of expertise in aerosol remote sensing, computer programming, and computing resources [62]. However, the current study is the first of its kind to utilize MODIS Level 2 aerosol data products at finer and multiple spatiotemporal aggregating domains involving 23 ground stations distributed at different topographic and climatic conditions of the western United States. In addition, no single study has been conducted to compare MODIS and AERONET aerosol products that are apparently focusing over the sub-continental western United States. In the previous studies, the western United States has been considered under coarser scales, either a continental [63] or global scale [62,64]. A continental study by Drury et al. [63] reported a correlation coefficient of 0.32 (R = 0.32) between MODIS and AERONET AOD in the western and central United States using a daily data for the period of 1 July-15 August 2004; although their correlation increased to 0.64 with a weekly averaging domain. Our study utilized the available resource of a finer temporal resolution of both data products in order to capture the dynamic nature of aerosol loading influenced by topographic and climatic conditions. On a regional scale, Loria-Salazar et al. [61] evaluated the MODIS Collection 6.0 aerosol product using an AERONET product in semi-arid Nevada and California, United States, during the summer of 2012. Similar to our results, their study found a weak correlation between both data products, a correlation coefficient of 0.316-0.447 during non-fire periods, and 0-0.557 during fire periods. However, our study covers 17 years of aerosol data records.

Conclusions
The current study contributes a comprehensive approach to the evaluation of remotely sensed satellite aerosol products using ground-based aerosol products by utilizing different aggregating domains covering a very large study area, different land cover/land use settings, and geographical locations. Overall, this study demonstrated that the augmentation of aerosol data products from satellite remote sensing and ground-based platforms can improve the characterization of aerosol loading in space and time across the western United States. It explored deeper into the finer temporal and spatial scenarios than previous studies by optimizing the aggregating approaches performing the long-term inventory of the aerosol data records from both the satellite and ground-based platforms. In order to assess the variability of aerosol concentrations in space and time from the AERONET sites, the study integrated high aggregating domains. Based on these in-depth approaches, this study identified stations and regions with a strong potential for evaluating and calibrating satellite products, and spotted stations and regions that showed avery low correlation with the satellite observation and stations with many data gaps and few data records.
In general, the majority of the stations revealed a significant correlation between MODIS and AERONET AOD at all spatiotemporal aggregating domains. Specifically, the stations from the higher latitude sites showed a stronger correlation, making the sites suitable for evaluating satellite retrievals. However, the stations from the southwestern United States desert showed the lowest correlation coefficients. In view of this region being a dominant source of mineral dust, these findings suggest that the MODIS AOD should be further improved and locally calibrated. Aerosol sources, such as windblown dust and anthropogenic fugitive dust (e.g., road dust from unpaved road), are short-lived and heterogeneous over space and time. Such source characteristics, along with bright, reflective desert surface, make it particularly challenging for satellite retrievals in this region. We recommend the application of medium spatiotemporal collocating resolutions (temporally ±15 min and spatially 30 km × 30 km) for the estimation of ground-based air quality parameters. Incorporating aerosol products from geostationary satellites like NASA's Geostationary Operational Environmental-16 (GOES-16) Advanced Baseline Imager (ABI) into the dataset used in this study can significantly improve the characterization of the aerosol loading and air quality assessment rather than using MODIS-and AERONET-derived AOD only. Aerosol products from GOES-16 ABI possess a higher temporal resolution than MODIS and they can mitigate the sparseness of aerosol observations from MODIS-AERONET couplings. More investigations are needed to understand the aerosol impacts of long-term land cover and land use changes around AERONET sites as these surface types are very sensitive to climate change and variability, and this parameter in turn is crucial in determining the response of the land-surface-to-satellite signals.

Acknowledgments:
The authors would like to thank the National Aeronautics and Space Administration (NASA) for providing MODIS Level 2 Aerosol Products Collection 6.1 from Terra and Aqua satellites, the Aerosol Robotic Network (AERONET), and the principal investigators of the 23 AERONET sites for allowing us to use ground-based Level 2 aerosol data products, the United States Environmental Protection Agency (EPA) for the ecological classification map of the United States, and the United States Geological Survey (USGS) for the digital elevation model (DEM). The authors acknowledge Carlos Montana of the Department of Earth, Environmental, and Resource Sciences from the University of Texas at El Paso. We also acknowledge Thomas F. Eck, Santiago Gasso, and Robert Levy for their helpful comments.

Conflicts of Interest:
The authors declare no conflict of interest. USDA is an equal opportunity employer and provider. Mention of trade names is for information purposes only and does not infer endorsement by or exclusion of equivalent products by USDA.

Abbreviations
The following abbreviations are used in this manuscript:   St. 6 n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a St. 13 n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a n/a