2.1 MERRA-2 aerosol reanalysis data

MERRA-2 is the latest atmospheric reanalysis version for the modern satellite era provided by the NASA Global Modeling and Assimilation Office (Gelaro et al., 2017), using the Goddard Earth Observing System, version 5 (GEOS-5), Earth system model (Molod et al., 2012, 2015), which includes atmospheric circulation and composition, ocean circulation and land surface processes, and biogeochemistry. Note that in MERRA-2, in addition to providing assimilation of traditional meteorological observations, a series of AOD observation datasets, including bias-corrected AODs, retrieved from the Advanced Very High Resolution Radiometer (AVHRR) instrument over the oceans (Heidinger et al., 2014) and MODIS (onboard both the Terra and Aqua satellites) (Levy et al., 2010; Remer et al., 2005), and non-bias-corrected AODs retrieved from MISR (Kahn et al., 2005) over bright surfaces and ground-based AERONET observations (Holben et al., 1998), were also assimilated within the GEOS-5 Earth system model. An overview of the MERRA-2 modeling system and a more detailed description of aerosols in the MERRA-2 system can be found in Gelaro et al. (2017) and Buchard et al. (2017), respectively. In this study, the 3-hourly MERRA-2 analyzed AOD fields, at a resolution of 0.5^∘ latitude by 0.625^∘ longitude, were used for evaluation, while the monthly mean AOD values were used for climate analysis.

2.2 Satellite aerosol data

Two AOD datasets during 2001–2016 retrieved from MODIS and MISR, both onboard the Terra platform, were used in this study. The MODIS sensor onboard the Terra satellite observes the Earth at multiple wavelengths (range: 410–1450 nm; 36 bands) with a 2330 km swath, which has provided near-daily global coverage since 2000 (King et al., 2003; Levy et al., 2015). This study employed the combined Dark Target–Deep Blue (DTB) AOD algorithm at 550 nm, with a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution, from the level 3 (L3) monthly global aerosol dataset for MODIS Terra, Collection 6.1. The average mean absolute error (RMSE) of the level 3 MODIS/Terra DTB monthly AOD data have been estimated to be about 0.075 (0.120) over land (Wei et al., 2019). Note that MODIS/Aqua L3 was not used because it started late (June 2002). In addition, compared with the linear trend in MODIS/Aqua AOD during the overlapping period (2003–2016), MODIS/Terra AOD shows a similar performance worldwide (including spatial–temporal consistency and distribution patterns of trend values) (Fig. S1), although the Terra sensor has been documented to suffer from degradation issues. The similar performance between MODIS/Terra and MODIS/Aqua is mainly attributed to a new calibration approach in the C6 version, which can remove major non-polarimetric calibration trends from the MODIS data (Levy et al., 2013, 2015; De Leeuw et al., 2018).

Total column AOD observations from the MISR sensor onboard the Terra satellite, which provides observations of the Earth's atmosphere with nine different along-track viewing zenith angles at four different spectral bands (440–866 nm) (Diner et al., 1998), were utilized. It should be noted that, although MISR has a much narrower swath (∼ 360 km) compared with MODIS, the multi-angle observation from MISR provides the capability for retrieving a more reliable AOD over bright surfaces, such as desert areas (Diner et al., 1998; Kahn et al., 2010). The AOD retrieval in the 555 nm channel from monthly global aerosol datasets at a spatial resolution of $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ were used in this study. The uncertainty of the MISR level 2 (L2) AOD data over land and ocean has been estimated to be ±0.05 or $\pm (\mathrm{0.2}\times \text{AOD})$ (Kahn et al., 2005). Note that the wavelength of AOD (555 nm) reported by MISR is different from that of the MERRA-2 and MODIS/Terra datasets (550 nm); however, this slight wavelength difference is not expected to affect our analysis and conclusions regarding AOD annual and seasonal trends.

2.3 Ground-based reference data: AERONET and CARSNET

Owing to the accuracy of ground-based AOD observations, long-term instantaneous AOD observation records from two independent operational networks – AERONET and CARSNET – were used to validate the 3-hourly MERRA-2 AOD values. Since there are not enough long-term AERONET observations in China, it was necessary to examine the performance of the MERRA-2 analyzed AOD fields using additional AOD observations from CARSNET. CARSNET is a ground-based network for monitoring aerosol optical properties that was first established by the China Meteorological Administration in 2002 (Che et al., 2009). Both AERONET and CARSNET use the same types of sun photometers, which can observe direct solar and sky radiances at seven wavelengths (typically 340, 380, 440, 500, 670, 870 and 1020 nm) within a 1.2^∘ full field of view at intervals of about 15 min (Holben et al., 1998; Che et al., 2009). For CARSNET, operating instruments are calibrated and standardized using CARSNET reference instruments, which in turn are regularly calibrated at Izaña, Tenerife, Spain, together with the AERONET program (Che et al., 2009, 2018). The cloud-screened AOD (based on the work of Smirnov et al., 2000) in CARSNET has the same accuracy as AERONET, with an estimated uncertainty of 0.01–0.02 (Eck et al., 1999; Che et al., 2009).

Figure 1Geographical locations of the AERONET (yellow dots) and CARSNET sites (magenta dots) used in this work. The red boxes represent the 12 regions of interest selected in this study: northeastern Asia (NEA), northern China (NC), southern China (SC), southeastern Asia (SEA), northwestern China (NWC), southern Asia (SA), Middle East (ME), western Europe (WEU), Sahara desert (SD), central Africa (CF), eastern United States (EUS) and Amazon zone (AMZ).

In this work, we collected ground-based AOD observations (more than 1 year of data) from 468 AERONET sites worldwide and 37 CARSNET sites in China. The locations of these ground-based sites are shown in Fig. 1. Detailed information about these AERONET and CARSNET sites is given in Tables S4 and S5. The combined instantaneous AOD data collected by AERONET (quality-assured and cloud-screened L2 data) during 1993–2016 and CARSNET (cloud-screened L2 data) during 2002–2014 were used. Moreover, to ensure the reliability of AOD evaluation, the AOD measurements in two adjacent channels (i.e., 440 and 675 nm) from AERONET and CARSNET were subsequently interpolated to 550 nm for MERRA-2, using a second-order polynomial fit to ln (AOD) vs. ln (wavelength) (Eck et al., 1999).

2.4 Emissions inventory and meteorological data

The anthropogenic emissions inventories used in this study were obtained from the Peking University (PKU) website (http://inventory.pku.edu.cn/, last access: 31 March 2019), including total suspended particles (TSP) (Huang et al., 2014), SO₂ (Su et al., 2011), BC (Wang et al., 2014) and organic carbon (OC) (Huang et al., 2015), with a spatial resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ and spanning the period 1980–2014. The emissions were calculated using a bottom-up approach based on fuel consumption and an emissions factor database. Huang et al. (2015) showed that the PKU emissions inventories are broadly similar to those of EDGAR version 4.2 (Edgar, 2011). Monthly meteorological fields from the MERRA-2 global reanalysis were also utilized, including total surface precipitation, surface wind speed, surface relative humidity (RH), mean sea level pressure, etc. These data have a spatial resolution of $\mathrm{0.5}{}^{\circ }\times \mathrm{0.625}{}^{\circ }$ and span the period 1980–2016 (Gelaro et al., 2017). For more detailed information on the selected meteorological parameters, see Table 1.

Table 1Prediction variables used in the stepwise MLR models.

^a Units: g km⁻² (TSP, SO₂, BC, OC); kg m⁻² s⁻¹ (Pre); m (PBLH, GH); 1 (SM, CLF); Pa (SLP); K (T, dT); % (RH); pa s⁻¹ (Ome) and m s⁻¹ (U,V, WS, VWS_500−850). ^b VWS_500−850 was calculated as $\sqrt{(U_{\mathrm{500}}-U_{\mathrm{850}}{)}^{\mathrm{2}}+(V_{\mathrm{500}}-V_{\mathrm{850}}{)}^{\mathrm{2}}}$.

Download Print Version | Download XLSX

2.5 Regions of Interest (ROIs)

In this study, 12 regions of interest (ROIs) dominated by different aerosol types were selected to study the long-term trends in regional aerosol loading and how they are related to local emission changes as well as the variation in meteorological variables. These 12 ROIs included three mineral-dust-dominant regions, SD (3–25^∘ N, 17^∘ W–20^∘ E), ME (14–33^∘ N, 38–56^∘ E), and northwestern China (NWC; 35–47^∘ N, 73–94^∘ E); three biomass-burning-dominant regions, the Amazon zone (AMZ; 1–22^∘ S, 46–60^∘ W), central Africa (CF; 2–18^∘ S, 12–33^∘ E), and southeastern Asia (SEA; 8^∘ S–18^∘ N, 96–127^∘ E); and six anthropogenic-aerosol-dominant regions, EUS (29–45^∘ N, 73–94^∘ W), western Europe (WEU; 37–59^∘ N, 10^∘ W–18^∘ E), southern Asia (SA; 10–30^∘ N, 72–90^∘ E), northern China (NC; 30–40^∘ N, 108–120^∘ E), southern China (SC; 20–30^∘ N, 108–120^∘ E) and northeastern Asia (NEA; 30–41^∘ N, 125–145^∘ E). The geographical boundaries of these ROIs are shown in Fig. 1.

Figure 2Flowchart with the procedure followed for (a) the evaluation of MERRA-2 global AOD, using the AERONET and CARSNET ground-based reference dataset, and (b) the evaluation of global and regional AOD trends.

Download

2.6 Statistical analysis

2.6.2 Trend analysis and stepwise MLR model

Long-term trend analysis of the AOD from MERRA-2, MODIS/Terra and MISR was performed on monthly time series data, using ordinary least-squares linear regression – a technique widely employed for trend analysis of aerosol data (Hsu et al., 2012; Pozzer et al., 2015; Klingmüller et al., 2016; Ma et al., 2016; Hammer et al., 2018). Prior to regression, these data were first deseasonalized by subtracting the monthly mean for different study periods for each grid cell to eliminate the large influence of the annual cycle. To better compare the results of the trend analysis, the MERRA-2 and MISR datasets at high spatial resolution ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.625}{}^{\circ }$ and $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$, respectively) were bilinear interpolation to the MODIS/Terra resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ (see Fig. 2b for the whole procedure). Incomplete sampling from the satellite instruments may introduce biases in long-term trend analysis. Thus, to ensure the reliability of the trend analysis, each grid cell for the MISR and MODIS/Terra AODs was required to have valid data for at least 60 % of the time period before regression was performed. Two-tailed Student's t tests were used to assess the robustness of each trend estimate, and the criterion for statistical significance was set at the 95 % confidence level.

Pearson's R was used to measure the strength of the relationship between AOD, anthropogenic emissions and meteorological parameters. MLR models of monthly MERRA-2 AODs were built for the 12 ROIs using emission factors, meteorological parameters or both as predictors. A total of 4 emission factors and 32 meteorological parameters were considered in the MLR models (Table 1). For each ROI, the MLR model could be expressed as

$$\begin{array}{}\text{(8)}& y={\mathit{\beta }}_{\mathrm{0}}+\sum _{i=\mathrm{1}}^n{\mathit{\beta }}_i{x}_i+\mathit{\epsilon },\end{array}$$

where y is the standardized monthly AOD and (x₁, …, x_n) is the ensemble of standardized monthly explanatory variables. The standardized regression coefficient β_i was determined by the least-squares method, and ε is an error term.

In each step of the MLR model, a variable is considered to be moved or removed from the set of explanatory variables using the stepwise regression method to obtain the best model fit. In other words, for each step the model adds a significant (P<0.05) explanatory variable to the model, it can be removed only if it is insignificant (P>0.1) after adding or removing another variable. A similar model has been widely used to investigate the relationship between aerosols and meteorology (e.g., Yang et al., 2016; Lu et al., 2016; Zhai et al., 2019; Tai et al., 2010).

Although the most important explanatory variables were obtained via the above stepwise MLR model, there might be multiple collinearities among different explanatory variables. In that situation, the standardized regression coefficient as an explanation of relative importance is unstable and misleading. To eliminate the influence of multi-collinearity, the variance inflation factor (VIF) (Altland et al., 2006) was used to test whether there was a multi-collinearity problem among the variables. VIF is often regarded as a measure of collinearity between each variable and another variable in the model. VIF can be calculated from the following relationship:

$$\begin{array}{}\text{(9)}& \mathrm{VIF}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-R_i^{\mathrm{2}}}},\end{array}$$

where $R_i^{\mathrm{2}}$ is the coefficient of determination of linear regression between the ith independent variable and other independent variables in the model. The present study used a VIF threshold of 10, as is widely recommended in the literature (e.g., Hair et al., 2007; Barnett et al., 2006; Field, 2005), to represent the maximum acceptability of collinearity.

Finally, to better quantify the relative contributions of each independent explanatory variable, which were obtained from the stepwise MLR model, to AOD variability, the LMG method (Bi, 2012; Grömping, 2006; Lindeman et al., 2014) was applied. This approach is one of the most advanced methods for determining the relative importance of explanatory variables in a linear model and provides a decomposition of the fraction of model-explained contributions (i.e., R²) into nonnegative contributions using semi-partial R values. The LMG measure for the ith regressor x_i can be expressed as

$$\begin{array}{}\text{(10)}& \mathrm{LMG}\left(x_i\right)={\displaystyle \frac{\mathrm{1}}{p\mathrm{!}}}\sum _{r\text{permutation}}\text{seq}{R}^{\mathrm{2}}\left(\left\{x_i\right\}\mathrm{|}r\right),\end{array}$$

where r represents the rth permutation (r=1, 2, …, p!) and seqR²({x_i}|r) represents the sequential sum of squares for the regressor x_i in the ordering of the regressors in the rth permutation.

Figure 3Evaluation of the 3-hourly MERRA-2 AOD against the (a) AERONET and (b) CARSNET AODs. The color-coded dots indicate the number of samples. The solid red line is the line of best fit and the dashed black line is the 1 : 1 line. For descriptions of statistical metrics, see the comparison methods section.

Download

Figure 4Comparison of the 3-hourly MERRA-2 AOD datasets with AOD observations of 468 AERONET sites worldwide and 37 CARSNET sites in China: site performance maps for the (a) correlation coefficient (R), (b) mean absolute error (MAE), root-mean-square error (RMSE), (c) relative mean bias (RMB), (d) mean fractional error (MFE), (e) fractional gross error (FGE), and (f) the index of agreement (IOA) between MERRA-2 AOD and ground-based AOD observations. The size of the circles in (b) represents the RMSE and their inner color represents the MAE. The bars in the lower-left inset in each panel represent the frequency distribution histograms for the R, MAE, RMSE, RMB, MFE, FGE and IOA between MERRA-2 and all ground-based observations incorporating AERONET and CARSNET, respectively. Note that all sites within each region of interest (ROI) are integrated to assess the accuracy of the MERRA-2 AOD dataset in that area. The performance of the MERRA-2 AOD dataset in each ROI is illustrated in Figs. S2 and S3.

For a detailed introduction to and description of the calculation process of the LMG measure, refer to Grömping (2006). For all variables (including the AODs from MERRA-2, MISR, and MODIS/Terra; the meteorological variables from MERRA-2; and the emission estimates from PKU), the regional mean was calculated by averaging valid variable values over all grids within the 12 ROIs. For the seasonal analysis, the four seasons were considered as follows: spring (March–April–May), summer (June–July–August), autumn (September–October–November) and winter (December–January–February).

3.1 Assessing the performance of the MERRA-2 AOD datasets on the global scale

Although the official documentation points out that a large number of AOD observations have been assimilated into the system (Buchard et al., 2017), the global performance of MERRA-2 AOD is still unknown. In addition, since MERRA-2 assimilates a variety of AOD datasets from different observation periods (such as AVHRR before 1999, AERONET since 1999 and Earth Observing System (EOS) era satellites after 2000) (Buchard et al., 2017), it is difficult to disentangle the influence of each assimilated dataset alone on the overall accuracy of MERRA-2. Strictly speaking, we need to point out that MERRA-2 and AERONET are not independent of each other (after 1999). Nevertheless, we hope that this assessment will still provide some reference for other studies using the MERRA-2 AOD dataset. Therefore, we first use AERONET to evaluate the overall performance of MERRA-2 AOD on the global scale and then use CARSNET to independently examine the performance of the MERRA-2 analyzed AOD field in China.

Table 2Statistical measures of the 3-hourly MERRA-2 AOD vs. AERONET and CARSNET AODs over the 12 regions of interest.

^* Indicates the statistical results for CARSNET sites.

Download Print Version | Download XLSX

3.1.3 Global AOD distribution and interannual evolution of regional AOD

Figure S6 shows the global annual- and seasonal-mean AOD distribution calculated from the MERRA-2 AOD products during 1980–2016. Furthermore, the distributional characteristics of the global annual-mean AOD from MERRA-2, MODIS and MISR during the same period (2001–2016) are also compared in the figure. The comparison shows that, although MISR underestimated the AOD (e.g., in SA and eastern China), as expected because of insufficient sampling (Mehta et al., 2016; Kahn et al., 2009), the three AOD products were generally closely consistent on the global scale (also see Fig. S5). Generally, high AOD loading was mainly observed in areas of high anthropogenic and industrial emissions, such as in eastern China and India, and major source areas of natural mineral dust – particularly the Saharan, Arabian and Taklimakan deserts.

Figure 5Temporal evolution of regional monthly averaged AOD for the 12 regions of interest. Each year is represented by an irregular ring with 12 directions. Each direction of the ring represents a specific month, the distance from the center of the ring represents the regional monthly mean AOD value and the color of the ring represents the year. A special cyan-colored ring represents the monthly mean AOD for the period 1980–2016.

Download

Due to the seasonal variation in the atmospheric circulation driven by solar radiation and the intensity of human activities in different regions, the global distribution of AOD also shows obvious seasonal differences, with global aerosol loading reaching its maximum in spring and summer. On the one hand, this can mainly be attributed to the enhanced circulation in spring and summer, which increases the likelihood of natural mineral dust from several major dust sources in the Northern Hemisphere (i.e., the Sahara and Sahel, the Arabian Peninsula, Central Asia, and the Taklimakan and Gobi deserts) being brought into the atmosphere; additionally, along the westerly belt airflow dust can be transmitted to surrounding sea areas (such as the strip of the northern tropical Atlantic stretching between West Africa and the Caribbean, the Caribbean, the Arabian Sea, and the Bay of Bengal) and more remote areas (such as South America; the Indo-Gangetic Plain; and the eastern coastal areas of China, Korea, and Japan) (Mao et al., 2014). On the other hand, higher temperatures and damp air in summer can create favorable conditions for the hygroscopic growth and secondary formation of aerosols (Minguillón et al., 2015; Zhao et al., 2018), which raises the AOD in some areas, such as NC and northern India, dominated by anthropogenic aerosol emissions in summer. Moreover, frequent local biomass-burning aerosol emissions in central Africa during summer are the main cause of high AOD in the region (Tummon et al., 2010).

In contrast, global aerosol loading is relatively low in autumn and winter. The atmosphere in autumn and winter is generally more stable and vertical mixing is weaker, and thus it is difficult for more aerosols – particularly natural mineral dust – to be brought into the atmosphere, which leads to lower AOD in autumn and winter (Zhao et al., 2018). Nevertheless, the AOD in autumn in South America, SEA, SC and CF is clearly high, which is mainly attributable to the emission of large amounts of fine aerosol particles (i.e., BC and OC) from frequent biomass burning in these regions (Thornhill et al., 2018; Ikemori et al., 2018; Chen et al., 2017). Notably, fine particulate matter composed of sulfate–nitrate–ammonium aerosols, which is produced by high-intensity anthropogenic activities in autumn and winter, is still the main contributor to high AOD in eastern China and India (Gao et al., 2018; David et al., 2018).

To better characterize the temporal evolution of regional AOD, the monthly mean AODs over the 12 ROIs from 1980 to 2016 were calculated. As illustrated in Fig. 5, the monthly regional AOD had large seasonal variability, in addition to varying degrees of fluctuation in different periods. In areas dominated by smoke aerosols from biomass burning (i.e., AMZ, CF and SEA), biomass-burning events tend to occur in the warm season (May–October), leading to a more prominent monthly AOD at this time of the year compared with the cold season (November–April). It is noteworthy that MERRA-2 also captured several well-known forest-fire events, such as those in Indonesia in 1983 and 1997, which have been proven to be mainly related to climatic drying caused by El Niño and large-scale deforestation (Page et al., 2002; Goldammer, 2007). In the CF region, the monthly mean maximum AOD experienced a transformation process – that is, the monthly maximum AOD often occurred in June and July before 2000, whereas after 2000 it occurred more frequently in August and September. This shift may be attributed to the fact that MERRA-2 did not assimilate any land-based AOD observations before 1999, which made it difficult for the model to simulate the monthly variation in regional AOD (Gelaro et al., 2017; Buchard et al., 2017). In the AMZ and SEA regions, September and October seem to be the two most frequent months for the occurrence of high AOD values, but the magnitude of AOD values has decreased in recent years, which may be related to changes in meteorological conditions (Torres et al., 2010).

Figure 6Spatial distributions of the linear trends in annual and seasonal MERRA-2 AOD calculated from the time series value of the deseasonalized monthly anomaly during (a) 1980–2016, (b) 1980–1997 and (c) 1998–2016. Only trend values with statistical significance at the 95 % confidence level are shown.

In areas dominated by natural mineral dust aerosol (i.e., the SD, the ME and NWC), the monthly maximum AOD mainly occurred in March–August. Before 2000, there were many anomalies in the AOD monthly maximum, which also implied frequent sandstorms. In contrast, the frequency of monthly AOD anomalies decreased after 2000, which may be attributable to the reduced surface wind speed and increased vegetation cover (Kim et al., 2017; Wang et al., 2018; An et al., 2018). Compared with the areas dominated by smoke and dust aerosols, the seasonal differences of AOD in the areas dominated by anthropogenic aerosol emissions appear to be smaller, but their temporal evolution is more pronounced. In NEA, the monthly maximum AOD often occurred in March–June, possibly related to the long-distance transportation of sand and dust in the China–Mongolia deserts (Taklimakan and Gobi). However, as the frequency of sandstorms has decreased in the past 10 years (An et al., 2018), the monthly maximum AOD has also shown a downward trend. In NC and SA, the monthly AOD has gradually expanded outward since 1980, indicating that AOD has experienced a gradual increase. Monthly AOD had large seasonal variability in the SC region, reaching its maximum in February–April. The increased aerosol emissions from biomass burning in spring seem to be one of the main reasons for high AOD in the SC region (Chen et al., 2017). For the EUS and WEU regions, the characteristics of the monthly variation in AOD were similar, i.e., large values of AOD occurred in summer. With time, the monthly AOD showed a tendency to gradually shrink inwards, suggesting AOD has experienced a significant decline over the past few decades in the EUS and WEU. The main drivers of the interannual variability of AOD over each ROI are discussed in detail in Sect. 3.5 and 3.6.

3.2 Global AOD trend maps

Annual and seasonal linear trends of the MERRA-2 AOD anomaly were separately calculated for each $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell for the whole of 1980–2016 period (period 1) and for the first 18 years (1980–1997, period 2) and last 19 years (1998–2016, period 3). Figure 6 shows the spatial distribution of these trends on the global scale. Throughout period 1, the regions where annual AOD showed a significant upward trend (p<0.05) were mainly located in eastern China, SA, the ME, northern South America and the southern coastal areas of Africa, whereas some significant downward trends were observed in the whole of Europe and the EUS. However, compared with the annual trends, the seasonal AOD trends had obvious regional differences in terms of their spatial distribution. For instance, a strong positive trend throughout East Asia, including Korea and Japan, was found in spring. In summer, there was a significant upward and downward AOD trend in Siberia (and other areas of Asian Russia) and the Amazon basin, respectively. In contrast, winter AOD had a significant downward trend in the area north of 40^∘ N. These differences in seasonal trends are closely related to the seasonal variations in anthropogenic aerosols generated by local emissions and natural aerosols driven by meteorological conditions (De Meij et al., 2012; Chin et al., 2014).

Figure 7Spatial distributions of annual and seasonal trends in AOD calculated from the time series value of the deseasonalized monthly anomaly from (a) MERRA-2, (b) MODIS/Terra, and (c) MISR between 2001 and 2016. Only trend values with statistical significance at the 95 % confidence level are shown.

In the two different historical periods (i.e., period 2 and 3), these trends seem to have experienced a remarkable shift. During period 2, the annual AOD had a significant upward trend throughout the Southern Hemisphere, and similar upward trends also existed in eastern and northwestern China. This upward trend in the Southern Hemisphere, which was most likely associated with two giant volcano eruption events in the early 1980s (El Chichón, Hofmann and Rosen, 1983) and early 1990s (Pinatubo volcanoes, Stenchikov et al., 1998; Bluth et al., 1992; Kirchner et al., 1999), is also reflected in the regional annual mean AOD time series shown in Fig. S4. The eruptions led to a strong increase in volcanic ash and SO₂ emissions, consequently increasing AODs from place to place via airflow transport, which was captured accurately by MERRA-2. Meanwhile, AOD had a significant downward trend throughout Europe and the EUS, which appears to be related to the reduction of TSP and SO₂ emissions (see Sect. 3.5). Seasonally, a significant upward trend seems to be prevalent in all seasons in the Southern Hemisphere. Compared with other seasons, the decline of AOD was more obvious in Europe and America. In winter, except for the positive trend that still existed in the marine area of the Southern Hemisphere, the fluctuations in other regions were smaller and relatively stable.

During period 3, AOD began to show a significant upward trend in most regions, especially in SA, SEA, the ME, central Russia, the western United States and northern South America, while still maintaining an upward trend in eastern China with greater intensity. These upward trends over SA, the ME and eastern China are in good agreement with the results of Hsu et al. (2012), who used SeaWiFS AOD records from 1997 to 2010. It is worth noting that the trends for the whole of Europe shifted from significantly positive to statistically insignificant, while the region that had shown a significant downward trend before 1997 in the EUS was also shrinking. Furthermore, the region showing a positive trend, prevailing in the Southern Hemisphere, shrunk dramatically. Similarly, the spatial distribution of the trend also had significant differences in different seasons of this period. In spring and winter, only significant upward trends could be observed on a global scale, mainly in eastern China, SA, the ME and South America. Conversely, significant downward trends were apparent in the EUS, northwestern Africa and central South America in summer. Additionally, it was also found that the region with a significant downward trend in Africa shifted from the northwest in summer to the southwest in autumn. The joint effect of the changes in local emissions and meteorological conditions determined these trends in these regions. See Sect. 3.5 for a more detailed explanation.

Figure 8Intercomparisons of global and regional annual trends in AOD calculated from the time series value of the deseasonalized monthly anomaly of MERRA-2, MODIS/Terra and MISR, during the four periods of 1980–2016, 1980–1997, 1998–2016 and 2001–2016. Error bars represent the uncertainty associated with the calculated trend. The hatched trend bars indicate statistical significance at the 95 % confidence level.

Download

Ensuring the accuracy of AOD trends calculated by MERRA-2 is critical for quantifying the contribution of local emissions and meteorological factors to the inter-decadal variation in AOD in different regions. For comparison, the resulting annual and seasonal trends of the MERRA-2, MODIS/Terra and MISR AOD anomaly over the whole globe were derived, using the same method, between 2001 and 2016; the results are shown in Fig. 7. This comparison shows that the AOD trends during 2001–2016 calculated by MERRA-2 in most regions of the world agreed well with the results of MODIS and MISR, on both annual and seasonal timescales. Although MERRA-2 assimilates MODIS and MISR at the same time, the relatively small difference between MERRA-2 and MISR may be mainly due to the insufficient sample size of MISR (MODIS produces 3 to 4 times more data than MISR) (De Meij et al., 2012).

For the annual trend, the significant upward trend observed by MODIS/Terra and MISR in SA and the ME and the significant downward trend observed in the EUS, WEU and central South America were consistent with the results of the MERRA-2 trend. Similar trends were reported in a previous study based on 14 years (2001–2014) of observational records (Mehta et al., 2016). Similarly, upward trends also existed in spring, autumn and winter, while downward trends were also apparent in spring, summer and autumn. It should be noted that the trend signals calculated from MERRA-2 and MODIS/Terra were opposite in SC. The difference in sign associated with trends during 2001–2016 could mainly be due to the larger deviation between MERRA-2 and MODIS/Terra between 2001 and 2004 (Fig. S4c). The large deviation directly led to a reversal of trend throughout the period 2001–2016. This deviation may be related to the use of different versions of MODIS data: in the MERRA-2 AOD observing system, MERRA-2 assimilated the bias-corrected AOD derived from MODIS radiances, Collection 5 (C5; Buchard et al., 2017), and the MODIS data used in this study was the latest collection (Collection 6.1, C6). Different versions mean differences in algorithms (Fan et al., 2017), which may affect the statistical error.

3.3 Regional AOD trends

To examine the spatial and temporal changes in more detail, the annual trend over the globe and in the 12 ROIs, derived based on MERRA-2 during periods 1, 2 and 3, were calculated. In addition, for comparison purposes, the regional trends in AODs from MERRA-2, MODIS and MISR during 2001–2016 were also estimated. It should be noted that different datasets may have a certain effect on global and regional trend assessment due to their own uncertainties. Nevertheless, we include them for completeness but exercise caution when interpreting the differences in trend values between different datasets. The comparisons of the magnitudes of global annual trends with these regional trends are summarized in Fig. 8 and Table S1. In general, the annual trends derived from different datasets were small on the global scale. As indicated by the results in Fig. 8 and Table S1, the trend values were −0.00068 yr⁻¹ for the globe during period 1, with statistical significance at the 95 % confidence level. In contrast, no statistically significant trend was detected at the global scale for period 2 (0.00050 yr⁻¹) or 3 (0.00038 yr⁻¹). Analyzing the global AOD trends during 2001–2016 from MERRA-2 and the two satellite datasets, it was found that the MERRA-2 trends were negligible, whereas significant positive (negative) trends were found for MODIS (MISR).

However, the trends could be considerable on regional scales. For example, over the anthropogenic-aerosol-dominant regions for periods 1, 2 and 3, strong positive trends were apparent over NEA, NC, SC and SA, while strong and statistically significant negative trends were found over WEU and EUS. For biomass-burning regions (SEA, CF and AMZ but not CF, which had a negligible and insignificant trend), there was a positive trend during periods 1, 2 and 3. For the mineral-dust-dominant regions, although there seemed to be an upward trend over the ME, the estimated trends were not statistically significant for other areas, such as NWC and the SD. During 2001–2016, the estimated MERRA-2 AOD trend in most ROIs (i.e., NEA, SA, ME, WEU, EUS and AMZ) was comparable to and had the same sign as the trend from both the MODIS and MISR sensors. However, it was opposite in sign to the MISR data over NC, NWC and the SD and to the MODIS data over SC, SEA and CF during overlapping years. These differences in global trends between MERRA-2 and satellites may be related to several aspects, including the difference in sample number, data accuracy, different measurement methods, etc. (De Meij et al., 2012).

In addition to the annual trend, the seasonal trend of AOD for different datasets in different ROIs and different historical periods was also studied (Fig. S7 and Table S1). Globally, negative trends were observed throughout the four seasons during period 1, especially during summer, autumn and winter (−0.00078, −0.00092 and −0.00097 yr⁻¹, respectively; statistically significant at the 95 % confidence level). In contrast, there was a negative trend in period 2, although it was not significant. In the subsequent period, period 3, the trend values shifted from negative to positive. The positive trend was more significant in spring and autumn (0.00053 and 0.00070 yr⁻¹). Regionally, strong positive trends were apparent over both NC and SC throughout the four seasons during periods 1, 2 and 3. Strong upward trends were also found over SA. These upward trends were most likely associated with an increase in urban and industrial pollution in China and India. Meanwhile, some similar but relatively moderate upward trends also existed over NEA in spring. In contrast, strong negative trends were observed over the WEU and EUS regions, especially during spring, summer and autumn. The negative trends over WEU and the EUS may partly have been due to a decrease in polluting aerosols associated with emission control measures (De Meij et al., 2012; Li et al., 2014). A statistically significant upward trend was also found over the SD, NWC and the ME in spring during periods 1, 2 and 3 (0.00252, 0.00300 and 0.00463 yr⁻¹), respectively. In contrast to the strong downward trends over AMZ in summer during periods 1, 2 and 3, there appeared to be upward trends in spring over AMZ and in winter over CF and AMZ. When compared with the regional trends during 2001–2016 calculated by the two satellite datasets, we found that the seasonal trends of MERRA-2 were highly consistent with the satellite results in almost all regions, especially in spring and autumn. It is worth noting that the trend differences among the three different datasets in all four seasons still existed in NC and SC, and the differences had different seasonal characteristics. For example, over NC, the most significant difference occurred in spring and summer, whereas it occurred in summer and winter over SC. Seasonal differences in trends are mainly due to insufficient accuracy of MERRA-2 in China (see Sect. 3.1.2).

Figure 9Sliding-window trend analyses of the annual mean MERRA-2 AOD from 1980 to 2016 over the 12 ROIs (see Fig. 1 for names and locations of regions), with at least 10 years used to calculate trends. The x axis and y axis indicate the start year and the length of the time series to calculate the trend, respectively. The colors of the rectangles represent the intensity of the trend (units: yr⁻¹) and those with black “x” signs indicate linear trends above the 95 % significance level.

Since the sign of a trend value often varies with the span of the calculation period, it was necessary to evaluate the sliding trend of different periods to help examine the time node of the changes more precisely. Therefore, sliding trend analyses were used to present a more comprehensive analysis of annual trends over the 12 ROIs during different historical periods (Fig. 9). These trends were calculated for all periods, starting each year from 1980 to 2007 and ending in 2016 with increments of at least 10 years. As shown in Fig. 9, in the EUS and WEU the AOD experienced a large decline up until the 1981–1990 period, then the trend reversed moderately from 1984 to 1986, declined sharply from 1989 after a short increase from 1996 to 1999 and sustained a moderate downward trend in the last 17 years. A similar pattern was found for NWC, the SD and AMZ, although there was a stronger upward trend and relatively weaker downward trend in the corresponding period. In SC and NC, the AOD experienced a slight increase in the 1980s and a short-term decline around the 1990s, and then showed its largest positive trend since 1995, before reversing sharply over the last 10 years (Sun et al., 2019). A similar evolution also existed in NEA and the ME, although the intensities of the trends were relatively weak. In addition to the negligible downward trend in the 1980s and 1990s, SA showed overall positive trends throughout the period, corresponding to increasing anthropogenic emissions (Fig. 11). Furthermore, in CF, a moderate increasing trend was detected from 1983 to 1985 and then again in 1990, and the trends became relatively stable but unexpectedly showed sharp increases after 1993, followed by a significant decline in the 2000s and reversal in the last 10 years. The trends for SEA were much smaller and relatively stable. Also, note that around 1985 and 1990 two distinct opposite trend signs were found in all regions. These two unexpected trends indicated that large volcanic eruptions not only greatly affect short-term changes in local aerosols but also impose different degrees of disturbance in long-term trends of aerosols in different regions of the world (Hofmann and Rosen, 1983; Stenchikov et al., 1998; Kirchner et al., 1999).

Figure 10Temporal evolution of sliding decadal trends in the annual and seasonal mean AOD from MERRA-2, MODIS/Terra and MISR over the 12 ROIs. The trends were calculated for each 10-year interval from 1980 to 2007 for MERRA-2 and from 2001 to 2007 for MODIS/Terra and MISR. The colors of the rectangles represent the intensity of the decadal trend (units: yr⁻¹) and those with black “x” signs indicate linear trends above the 95 % significance level.

Furthermore, considering that aerosol concentration and composition usually have strong seasonal cycles (Li et al., 2018), the trends for each season were also calculated separately and compared with the MODIS and MISR trends in the period of overlap (2001–2016). Note that Fig. 10 only shows the evolution of seasonal and annual trends for every 10-year period starting from 1980 to 2007 for MERRA-2 and from 2001 to 2007 for MODIS and MISR; refer to Figs. S8–11 for a fuller presentation of the regional seasonal trend. For all regions, the trends for all seasons, except autumn in SEA, CF and AMZ and spring in the SD, were in phase with the annual trend (also see Fig. S12). In general, autumn trends over SEA, CF and AMZ were larger and often out of phase, possibly attributable to the sudden increase in aerosol concentration caused by biomass-burning events. Similarly, the spring trend over the SD was also larger and more asynchronous than in other seasons. This phenomenon can mainly be attributed to active spring dust events (Liu et al., 2001). In addition, compared with the annual and seasonal regional trends during 2001–2016 (Figs. 8 and S7), the decadal trends of MERRA-2 agreed better with the trend results from MODIS and MISR. This implies that the trends can change relatively quickly with time (Li et al., 2018). Supporting evidence was also found from the strongest trends on both annual and seasonal scales being mostly concentrated in the lower y axis values (Figs. 9 and S8–S11). These results also highlight the importance of evaluating temporal shifts or decadal AOD trends.

3.4 Response of inter-decadal variation in regional AOD to local emissions and meteorological parameters

Previous studies have shown that the interannual variations in regional AOD are mainly controlled by changes in emissions and meteorological factors (De Meij et al., 2012; Pozzer et al., 2015; Itahashi et al., 2012; Zhao et al., 2017; Lee et al., 2016; Chin et al., 2014). First, the trends of the four emission factors (i.e., TSP, SO₂, BC and OC) and their correlations with AOD were calculated for the whole study period (1980–2014), as well as for two individual periods (i.e., 1980–1997 and 1998–2014). Note that the PKU global emissions inventories were only available for 1980–2014, which limited our research to a relatively short period. Figures 11 and S13 show the linear trends in emissions and their relationships with MERRA-2 AOD during 1980–2014, respectively. The decreasing AOD trends over Europe and the EUS (see Fig. 6) coincided with substantial reductions in the emissions of primary anthropogenic aerosols (TSP and BC) and precursor gases (SO₂), corresponding to pollution controls (Hammer et al., 2018; De Meij et al., 2012). This was also supported by significant positive correlation between AOD and emissions in most regions of Europe and the EUS (Fig. S13).

Positive trends in TSP and SO₂ were present over India and eastern China, which explained the significant upward trend of AOD in these two regions. In addition, eastern China and India experienced a shift in the emissions trend during the two periods (Figs. S14 and S16). In 1980–1997, a significant upward trend existed in both regions (Huang et al., 2014). In contrast, in 1998–2014, India at least maintained this upward trend for all four emission factors, with it sometimes being even stronger, while the positive trends in emissions of TSP and SO₂ over eastern China were interspersed with negative trends. More importantly, the trend of BC and OC in eastern China reversed completely. The shift in these emission trends in eastern China can mainly be attributed to the implementation of multiple emission reduction policies (Zheng et al., 2018). The reductions in emissions were at least partly responsible for the decreasing trend of AOD in the NC and SC regions in the last 10 years (see Fig. 9). The trends in primary BC emissions followed a similar pattern to the trends in OC emissions, except there were positive trends over northeastern China and the positive (negative) trends over CF, AMZ and SEA (WEU and SC) were lower in magnitude, reflecting regional changes in fire activity. There were positive AOD trends in areas dominated by biomass burning (especially in CF and SEA), in response to increased BC and OC emissions. Because human activities are scarce in desert areas, there was no direct relationship between AOD and emissions, as expected. Therefore, this highlights the importance of studying how natural factors (here, this refers to meteorological parameters) control the interannual variation in AOD in different desert areas. Furthermore, it is worth noting that in the two short periods (especially 1998–2014), these regions with significant positive correlation shrunk and were no longer significant (Figs. S15 and S17), suggesting other factors such as meteorological parameters might be driving the interannual trend of regional AOD.

Figure 11Spatial distributions of linear trends (units: kg km⁻² yr⁻¹) in total anthropogenic emissions of total suspended particles (TSP), SO₂, black carbon (BC) and organic carbon (OC) during 1980–2014 derived from the Peking University emissions inventory (http://inventory.pku.edu.cn/, last access: 31 March 2019) (Huang et al., 2014). Only linear trend values with statistical significance at the 95 % confidence level are shown.

To investigate the roles of meteorological parameters in the decadal variation in AOD, Pearson's R values between AOD and meteorological parameters (a total of 32; see Table 1) and over the 12 ROIs for the three periods (i.e., 1980–2014, 1980–1997 and 1998–2014) were calculated. Some of these meteorological variables, such as surface precipitation, surface wind speed, wind velocity, RH and surface wetness, have been shown before to be correlated with regional AOD (Klingmüller et al., 2016; Pozzer et al., 2015; Chin et al., 2014; He et al., 2016). Correlation analysis showed similar correlation patterns between AOD and meteorological parameters for the three different periods over all ROIs. During the period 1998–2014, the correlation was generally stronger than in the other two periods (see Fig. S18), suggesting meteorological factors may have played a more important role in this period. In addition, these correlations seemed to be similar in regions dominated by the same aerosol type. For example, in the mineral-dust-dominated regions (i.e., NWC, the ME and the SD), AOD had a significant positive (negative) correlation with near-surface wind speed (soil moisture), suggesting that surface wind speed and soil moisture may be the main factors controlling the dust cycle, which is consistent with previous studies in the ME (Klingmüller et al., 2016). In the biomass-burning-dominated regions (i.e., SEA, CF and AMZ), AOD had a significant negative correlation with humidity-related meteorological parameters (such as surface precipitation, RH and soil moisture), implying that ambient humidity (including the atmosphere and soil) may be a direct correlation factor in controlling the frequency of biomass-burning events (Torres et al., 2010). In contrast, in the regions dominated by anthropogenic aerosols, the correlation was regionally dependent and their signs differed from place to place.

Correlation analysis cannot directly identify the main factors affecting the inter-decadal change of AOD in different regions. Here, MLR models were used to diagnose the influences of local anthropogenic emissions and other meteorological parameters on the inter-decadal variation in AOD over the 12 ROIs. Figure 12 shows the time series of monthly mean MERRA-2 and MLR model-predicted normalized AOD anomalies, which used the emission factors, meteorological parameters or both as input predictors over the 12 ROIs for the whole study period (1980–2014). Similar comparisons for the two individual periods (i.e., 1980–1997 and 1998–2014) are also presented in Figs. S19 and S20, respectively. Table S2 summarizes the predictors included in the MLR models and their performance for the three different periods over each ROI. The MLR models with both emissions and meteorological parameters as predictors generally reproduced the AOD values in most regions during 1980–2014, except for high AOD values (Fig. 12), which are discussed below. For all the ROIs, the MLR models explained most of the MERRA-2 AOD variability (R²=0.42–0.76). However, when meteorology and emissions alone were used as predictors, there were considerable differences in different ROIs. When emission factors alone were used as the predictor, it could account for more than 35 % of the AOD variability in regions dominated by anthropogenic aerosols and biomass burning (except NEA, 14 %), with the largest explanation occurring in NC (58 %). In contrast, in the mineral-dust-dominated regions (the SD and ME), emission factors contributed little (< 0.05 %) to the interannual variation in AOD (Fig. 11g and i). Moreover, emission factors contributed 37 % of the AOD variability in NWC, which is mainly because of the strong anthropogenic emission sources in northern Xinjiang (mainly encompassing Ürümqi, Korla, Kashgar, etc.). However, compared with meteorological factors, emissions were not the main factors driving the interannual change of AOD (Fig. 12e).

Figure 12Time series of MERRA-2 (in black) and modeled AOD monthly normalized anomalies from 1980 to 2014 over the 12 regions of interest. The coefficient of determination (R²) of the regression fit of the stepwise MLR model with emission factors (in blue), meteorology (in green), and both emissions and meteorology (in red) as predictors are given in the top-right of each panel.

Download

On the other hand, when meteorological factors were used as predictors in the MLR models, it was surprising that they explained a larger proportion of the AOD changes in all ROIs except NC and SEA, where emission factors accounted for slightly lower AOD changes of 42 % and 33 %, respectively. Further analysis indicated that this difference in contribution between emissions and meteorology seemed to be greater for the two shorter periods of 1980–1997 and 1998–2017 (see Figs. S19 and S20). Besides, it should also be noted that the total explained variances of the MLR model for 1980–1997 were generally lower than those of the MLR model for 1998–2014, in all ROIs. The difference can be explained by two reasons: (1) a greater number of high AOD anomaly values occurred during the period 1980–1997 (Figs. 12 and S19), especially in relation to the two volcanic eruption events in the 1980s and 1990s, which directly reduced the total explained variances of the MLR model, because the model only considers the inter-decadal variations in local emissions and meteorological factors and the large-scale transport of pollutants is not considered, and (2) meteorology and emissions were confirmed to explain more AOD changes during the period 1998–2014.

Figure 13The LMG method-estimated relative contributions (%) of total variances in the stepwise MLR model explained by the local emission factors (left-hand bars) and meteorological variables (right-hand bars) over the 12 regions of interest during three periods: (a) 1980–1997, (b) 1998–2014 and (c) 1980–2014. Note that meteorological parameters were combined as follows: temperature, T (T_s, T₈₅₀, T₇₀₀, T₅₀₀, dT_900−s, dT_850−s); geopotential height, GH (GH₈₅₀, GH₇₀₀, GH₅₀₀); relative humidity, RH (RH_s, RH₈₅₀, RH₇₀₀, RH₅₀₀); vertical velocity, Ome (Ome₈₅₀, Ome₇₀₀, Ome₅₀₀); and wind speed, WS (U₈₅₀, U₇₀₀, U₅₀₀, V₈₅₀, V₇₀₀, V₅₀₀, WS_s, WS₈₅₀, WS₇₀₀, WS₅₀₀, VWS_500−850). Refer to Table S3 for the detailed relative contributions of each variable in the stepwise MLR models.

Download

3.5 Relative contributions of local emissions and meteorological parameters to inter-decadal variations in regional AOD

Application of the LMG method (see Sect. 2) to the MLR model allowed the relative contributions of each anthropogenic emission type and meteorological factor to the inter-decadal variations or trend of regional AOD to be quantified. Figure 13 shows the relative contributions of the local emissions and meteorological factors to the changes in regional AOD for the period 1980–2014, as well as for 1980–1997 and 1998–2014, using both emissions and meteorology as predictors in the MLR model. During the period 1980–2014, over the anthropogenic-aerosol-dominant regions, SO₂ was the dominant emissions-driving factor, explaining 24.9 %, 15.2 %, 32.6 %, 21.7 % and 12.7 % of the variance of AOD over NC, SC, SA, WEU and the EUS, respectively (also see Table S3). The above results also confirm that particulate sulfate is the main contributor to fine-mode AOD in anthropogenic-aerosol-dominant regions (Itahashi et al., 2012; David et al., 2018). Meanwhile, wind speed (including surface and upper wind speed) was the dominant meteorological driving factor, explaining 11.4 %, 14.2 % and 17.9 % of the variance of AOD over NC, SC and the EUS, respectively. In addition, planetary boundary layer height, temperature (including surface temperature, upper temperature, and the temperature difference between the surface and upper atmosphere) and RH (including surface and upper RH) were the strongest meteorological driving factors over NEA, SA and WEU, contributing 30.2 %, 15.9 % and 21.5 %, respectively.

In contrast, over the biomass-burning-dominant regions, BC (OC) was the dominant emissions-driving factor over SEA (AMZ), explaining 27.7 % (24.0 %) of the variance of AOD. Meanwhile, soil moisture and RH were the top meteorological driving factors over SEA, AMZ and CF, contributing 11.7 %, 35.5 % and 28.5 %, respectively. Furthermore, over the dust-dominant regions, wind speed was the strongest meteorological driving factor, explaining 30.3 % and 29.8 % of the variance in AOD over NWC and the SD, respectively. In contrast to wind speed being the primary meteorological driving factor over NWC and the SD, it was the second-most important factor over the ME, while sea level pressure was the primary driving factor, accounting for 60.9 % of the variation in AOD. This large variance, explained by sea level pressure and significant anticorrelations of the AOD with it (see Fig. S18c), further confirms the previous studies' findings that frequent sandstorms over the ME often correspond to a large horizontal pressure gradient differences caused by the enhanced high-pressure system across the eastern Mediterranean Sea and the enhanced low-pressure system across Iran and Afghanistan (Hamidi et al., 2013; Yu et al., 2016).

By comparing the estimated results of the two independent study periods (i.e., 1980–1997 and 1998–2014), it was found that in almost all ROIs (except NC and AMZ), meteorological factors contributed a larger explained proportion of AOD changes during 1998–2014, which indicates that meteorological factors seem to be becoming increasingly more important in dominating the inter-decadal change of regional AOD. It is worth noting that, in addition to the increased explained proportion of SO₂ and BC, among these meteorological factors, the role of diffusion-related parameters (such as horizontal and vertical wind speed, representing horizontal and vertical diffusion, respectively) seems to be the most prominent. This is consistent with the findings of Gui et al. (2019), who found wind speed to be the dominant meteorological driver for decadal changes in fine particulate matter over SC, based on a 19-year record of satellite-retrieved fine particulate matter data (1998–2016).