Developing and diagnosing climate change indicators of regional aerosol optical properties

Given the importance of aerosol particles to radiative transfer via aerosol-radiation interactions, a methodology for tracking and diagnosing causes of temporal changes in regional-scale aerosol populations is illustrated. The aerosol optical properties tracked include estimates of total columnar burden (aerosol optical depth, AOD), dominant size mode (Ångström exponent, AE), and relative magnitude of radiation scattering versus absorption (single scattering albedo, SSA), along with metrics of the structure of the spatial field of these properties. Over well-defined regions of North America, there are generally negative temporal trends in mean and extreme AOD, and SSA. These are consistent with lower aerosol burdens and transition towards a relatively absorbing aerosol, driven primarily by declining sulfur dioxide emissions. Conversely, more remote regions are characterized by increasing mean and extreme AOD that is attributed to increased local wildfire emissions and long-range (transcontinental) transport. Regional and national reductions in anthropogenic emissions of aerosol precursors are leading to declining spatial autocorrelation in the aerosol fields and increased importance of local anthropogenic emissions in dictating aerosol burdens. However, synoptic types associated with high aerosol burdens are intensifying (becoming more warm and humid), and thus changes in synoptic meteorology may be offsetting aerosol burden reductions associated with emissions legislation.

The direction of change and the presence of significant trends are consistent for mean and extreme (P90) AOD in all regions, but the magnitude of the change is larger for extreme AOD, indicating a narrowing of the AOD probability distributions (Fig. 3a). Significant regional AOD trends are ~1% year −1 , while the magnitude of the extreme AOD trends are 1.2-1.4% year −1 in regions of decreasing AOD and 1.9% year −1 for the NW (Figs 2, 3, and S2). There is marked seasonality in some regions in terms of both the magnitude of and temporal trends in the aerosol-CIs. For example, extreme (P90) AOD significantly decreased in summer (the season of highest historical values), spring, and fall in NE, summer and fall in SE and MW (p-value = 0.06 for MW summer), and during fall in GPl. Conversely P90 AOD increased in summer and fall in NW (Fig. 3d).
The key utility of including two indices of spatial structure of the fields is illustrated by the divergent trends in these two aerosol-CIs. All regions exhibit decreased AOD spatial autocorrelation (AOD-I), but increased AOD spatial coherence (AOD-SC) is observed over the NW, GPl, MW, SE, and AK, and decreased AOD-SC is observed in the SW (p-value = 0.07), GPu (p-value = 0.15), and NE (Figures 2 and S3). Causes of these differences and the inter-annual variability in the aerosol-CI trends are discussed below.
Mean AE significantly increased across all eight regions, indicating a decrease in mean particle size (Figures 2  and S2). This shift to higher AE is observed across the probability distribution, implying a shift in fine mode aerosols to smaller sizes, as opposed to a relative increase in fine versus coarse mode aerosols (Fig. 3b). However, trends in the spatial metrics of AE are not uniform across the regions. Significant negative trends in AE-I are observed in NW, SW, GPu, MW, and AK (Figures 2 and S3), but only two regions exhibited significant changes in AE-SC and they showed different signs (increased in SW and decreased in NE). Thus, there is evidence that as the aerosol populations are, on average, decreasing in diameter at the regional scale, but there remain sub-regions within many of the NCA regions with high coarse mode concentrations (e.g., across all days, 50% of grid cells have AE ≤ 1.2 in the NW, SW, and GPu; Fig. 3b), possibly due to wind-blown dust events 24 .
Mean SSA and SSA-SC decreases are observed in all eight regions (Fig. 2). There are also decreases in SSA-I for all regions except SE where there were significant increases in SSA-I, although the significance of the trend is lower in GPl (p-value = 0.06) and AK (p-value = 0.16). It is noted that SSA is determined by the aerosol composition and the dynamic range of SSA in MERRA-2 is lower than observations 17,25 ( Figure S1); therefore the aerosol-CIs that relate to SSA must be viewed with caution in the current reanalysis product. However, these trends are consistent with a tendency towards a relatively more absorbing aerosol, thus reducing the net cooling SCIENTIFIC RePoRTS | (2017) 7:18093 | DOI:10.1038/s41598-017-18402-x from aerosols. Further, the trends in SSA-I and SSA-SC imply aerosol populations are becoming more spatially heterogeneous in terms of the relative contribution of absorption to total radiative extinction.
When applied to the U.S. NCA regions, the aerosol-CIs thus indicate substantial evolution of aerosol populations through time in ways that are relevant to regional climate forcing. Overall aerosol burdens have declined (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and on average aerosol populations have changed to become more dominated by smaller diameter and more absorbing aerosols. They are also evolving in a way that causes a decrease in spatial autocorrelation, but increases in spatial coherence.
Attribution of temporal trends in the aerosol-CIs. Attribution of observed trends in the aerosol-CIs, particularly deconvoluting changes resulting from changing anthropogenic emissions, natural emissions, and atmospheric conditions is critical to demonstrating the effectiveness of emission reduction policies, exploring and prioritizing potential climate change mitigation strategies, and making projections of possible future values of the aerosol-CIs. Thus, the aerosol-CIs for the NCA regions are examined below in the context of these key drivers of aerosol populations.
Aerosol-climate interactions are reciprocal. Aerosols are a major driver of climate variability and change, but equally changes in climate alter aerosol concentrations and composition [26][27][28] . Further, previous research has  (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). Upward and downward facing triangles indicate significant positive and negative trends as determined using Kendall's tau-b, while square markers indicate no significant trend (at α = 0.05). (c) Percentage change per year in the CIs estimated using a linear regression fit (shown in Figures S2 and S3). The middle circles denote the normalized regression slopes (i.e. trends), and the inner and outer circles are the lower and upper bounds, respectively, of the 95% confidence intervals of these slopes. Black circles indicate trends that are not significant at α = 0.05. SCIENTIFIC RePoRTS | (2017) 7:18093 | DOI:10.1038/s41598-017-18402-x illustrated a key role of synoptic scale meteorological conditions in determining regional aerosol concentrations under the current 3,29-31 and possible future climate 32,33 . Consistent with that research, in each of the NCA regions, a number of synoptic types (i.e. repeated meteorological patterns) derived in a PCA of MERRA-2 meteorological output are associated with 10-20% AOD anomalies (positive and negative from the mean) (Fig. 4). The link to meteorological conditions at the synoptic scale is less pronounced for AE (the anomalies are <10%) and it appears SSA is relatively insensitive of the prevailing meteorological conditions (no synoptic type had a regionally average SSA anomaly of >2%). This finding re-emphasizes the complexity of aerosol populations and their related climate forcing, and highlights the importance of having multiple aerosol-CIs in order to fully characterize changes in climate-relevant aerosol properties.
Over all regions, synoptic types characterized by cooler (or milder) and drier conditions are associated with lower AOD. Conversely, anomalously high AOD is associated with warm and/or humid synoptic types, consistent with enhanced AOD under stagnant flow 29 and aerosol growth by water uptake 34 . Over the northern and western regions of the contiguous U.S. (NW, SW, GPu, MW) southwesterly geostrophic flow is typically associated with positive anomalies in both mean and extreme AOD, while northwesterly flow is associated with negative anomalies in mean and extreme AOD (Fig. 4). Anomalously low AE in virtually all regions is often associated with cool, dry synoptic conditions, consistent with an increase in dust loading during dry conditions 24 . Conversely, high AE is associated with warm, humid conditions at the synoptic scale consistent with predominance of hygroscopic secondary aerosols.
Consistent with prior research that has indicated changes in global and regional temperature and humidity are likely to result in changing characteristics of the synoptic types 29,35 , the majority of synoptic types associated with large positive AOD anomalies in each region exhibit a significant positive trend in PC scores. Conversely, synoptic types associated with negative AOD anomalies exhibited trends that are divided between increasing and decreasing trends (Fig. 4). While there is evidence that some cool, dry days are also becoming cooler and drier, the dominant signal in this analysis is thus that synoptic types associated with elevated AOD are evolving to become more intense, i.e. warm, humid days becoming warmer and more humid. These changes in the synoptic-scale  26,28 . While the intensity of the synoptic types has changed, the frequencies of individual synoptic types over each region do not exhibit significant temporal trends over the period 2000-2015.
Consistent with policy enacted under the U.S. Clean Air Act that has resulted in declining anthropogenic pollutant emissions over the study period, regionally integrated emissions of key aerosol precursor species, sulfur dioxide (SO 2 ) and nitrogen oxides (NO x ), exhibit a significant negative trend for all eight NCA regions over the period 2000-2015. Further ammonia (NH 3 ) emissions exhibit a negative trend in all regions except the MW and NE, and volatile organic compounds (VOC) emissions exhibit a negative trend in all regions except the NW and SE ( Figure 5) 36 . Consistent with this, mean and extreme AOD significantly decreased in GPl, MW, SE, and NE, and seasonal extreme AOD decreased in the fall in GPl, summer and fall in MW and SE, and spring, summer, and fall in NE. The overall tendencies in aerosol-CIs, including the significant decrease in mean and extreme AOD over GPl, MW, SE, and NE, are thus consistent with a decrease in sulfate aerosol abundance due to the reduction in SO 2 emissions (e.g., correlation coefficients between annual SO 2 emissions and extreme summer (except GPl) and fall AOD are >0.57 over these regions). Congruent with this decline in SO 2 emissions, the annual deviations from the overall cumulative distribution functions (cdf) imply that almost the entire probability distribution of AOD has shown a shift towards lower values (Fig. 3a). Further, because sulfate has a high SSA (near unity) 37 , Figure 4. Mean synoptic conditions for synoptic types associated with anomalously low and high AOD for each region (locations shown in Fig. 1). The mean temperature at 700 hPa (in K) are shown by the background colors, the solid black lines depict the 500 hPa geopotential isoheights (in m), and the red, magenta, cyan, and blue stippling represent 700 hPa specific humidity anomalies −2, −1, +1, and +2 standard deviations from the mean for all days. The arrows beside the panels indicate the presence and direction of significant trends in the PC scores associated with these synoptic types. The abscissa and ordinate axes are longitude (degrees East) and latitude (degrees North), respectively.
SCIENTIFIC RePoRTS | (2017) 7:18093 | DOI:10.1038/s41598-017-18402-x a reduction in secondary sulfate aerosol would also contribute to the observed decline in regionally-averaged SSA. Reduced production of sulfuric acid may also lead to a reduction in mean aerosol diameter, implied by the increase in AE, due to a reduction in condensational growth. While historic trends in black carbon (BC) emissions are highly uncertain (e.g., from biomass burning), it is estimated emissions from mobile sources, the largest BC source in the U.S., decreased by 32% from 1990-2005 (ref. 38 ). Further, BC only contributes to ~4% of global AOD 18 . Thus changes in SSA are likely not due to changes in anthropogenic BC emissions. Secondary organic aerosols are also a substantial component of aerosol mass and AOD over much of the eastern U.S. 39 . Thus an additional contributory factor to declining AOD in these regions is the reduction in anthropogenic VOC emissions and secondary organic aerosol formation. Accordingly, the correlation coefficients between annual VOC emissions and extreme summer and fall AOD in the NE and MW are >0.61. Thus, consistent with prior research, historical temporal trends of AOD across much of the contiguous U.S. are strongly responsive to emission reductions associated with the Clean Air Act.
Despite reductions in anthropogenic aerosol precursor gas emissions, it is worthy of note that primary aerosol emissions exhibit a significant trend only in the NW, GPu, and MW (Fig. 5), and that biogenic VOC, dust, and wildfire emissions exert a substantial impact on aerosol burdens and optical properties 40,41 . For example, there is a clear peak in extreme AOD in the spring of 2011 in the GPl, MW, and SE when wildfire burned area in the GPl was approximately four times greater than any other year (Figs 3 and 5). In the GPl, the lack of association (i.e. lower correlation coefficients) between annual anthropogenic emissions and extreme AOD in three of the four climatological seasons and the observed decreased SSA may also be in part due to increased abundance of dust aerosols, consistent with remote sensing measurements that indicate increased dust-related absorption aerosol optical depth (AAOD) over the central U.S. 24 . The declining trend in AOD in AK is also not very strongly linked to changes in anthropogenic emissions, but there is a significant positive association between extreme summer AOD and wildfire burned area (r = 0.96). This is clearly evident in 2004, 2009, and 2015, when positive excursions in monthly burned area (Fig. 5) coincide with spikes in summer extreme AOD (Fig. 2).
Only the NW region exhibits a significant positive trend in annual mean AOD, with extreme AOD increasing in the summer and fall (Figs 2 and 3). This is despite declines in regional anthropogenic emissions (Fig. 5), and may reflect confounding influences from increased wildfires (seasonal burned area and extreme AOD in summer and fall exhibit co-variability with r = 0.53 and 0.75, respectively) and long-range transport. For example, Siberian fires in the summer of 2012 impacted air quality in the Pacific NW 41 , and are evident in high P90 AOD during the 2012 summer and fall (Fig. 2).
The decrease in the spatial autocorrelation in AOD (Figures 2 and S3) along with the decreased anthropogenic aerosol precursor emissions in each region (Fig. 5) indicates an increasing influence of local sources on sub-regional aerosol concentrations and thus increased grid cell-to-grid cell variability in aerosol populations. Conversely, scales of spatial coherence (distance at which grid cells become independent) are increasing, which may be linked to changes in synoptic scale conditions (Fig. 4). High and low AOD are generally associated with warm, humid and cool, dry conditions, respectively. The positive trend in PC scores for synoptic types associated with high positive AOD anomalies indicate a tendency towards intensification of meteorological conditions associated with large direct aerosol radiative forcing that may be offsetting some of the effects of emission controls. As climate conditions continue to evolve, this highlights the critical need to better understand the feedbacks between climate and aerosol populations.

Discussion
Use of climate indicators to represent key components of the climate system is an increasing focus of the U.S. NCA. For this reason, we advocate that aerosol-CIs are urgently needed to track a key aspect of the radiation balance of Earth, air quality, and biogeochemical cycles, and that aerosol-CIs should be generated and interpreted at the regional scale. The guidance for developing CIs is that they should be relatively straightforward to compute and readily evaluated in both the contemporary and possible future climate. Thus, the aerosol-CIs we propose can be readily derived for any gridded data set and therefore can be applied to any region using current and future generation reanalysis products and/or output from regional and climate models.
The aerosol-CIs presented herein are designed to be useful in tracking changes in climate relevant aspects of the aerosol population and to assist in diagnosing the causes of changes in aerosol populations at the regional scale. Their utility in the former regard is illustrated by application to the NCA regions, and specifically the finding that mean and extreme AOD and SSA is declining and AE is increasing over most of the U.S. consistent with a tendency towards lower aerosol burdens that are increasingly dominated by smaller diameter and relatively more absorbing aerosols. This implies a decline in the degree to which aerosols have offset greenhouse gas related warming of the climate over much of the contiguous U.S.
The aerosol-CIs are also defined using two geospatial metrics: Spatial correlation and spatial coherence. The former (Moran's I) characterizes normalized co-variability and is a measure of the degree to which daily fields of AOD, AE, and SSA exhibit spatial clustering. The latter is a measure of the distance (range in the semivariogram) at which spatial fields become independent, and thus the extent to which the aerosol forcing can impact regional climate. The utility of these two spatial metrics in terms of diagnosing causes of changes in aerosol populations at the regional level is also indicated by the presence of divergent trends in AOD-I and AOD-SC in the NCA regions. These findings imply a tendency towards more grid cell-to-grid cell variability in aerosol populations, due to declining regional precursor and aerosol emissions leading to an increase in the relative importance of local emissions, within larger areas of increased spatial coherence (i.e. large range values from the semivariograms) in part due to an increase in the intensity of the predominant modes of synoptic scale meteorology.
Future work is needed to examine aerosol trends in global regions outside of the U.S. that are characterized by markedly different emissions and climate trends. Additionally, analyses of reanalysis products are only as good as the assimilation data and model used to develop the product. Thus the CIs should be applied to future reanalysis SCIENTIFIC RePoRTS | (2017) 7:18093 | DOI:10.1038/s41598-017-18402-x products that assimilate improved bias-correction assimilated data, data from additional, recently launched sensors, and more sophisticated model frameworks with improve aerosol treatment and emissions inventories.

MERRA-2.
MERRA-2 is derived using assimilation of both meteorological and aerosol observations every 6 and 3 hours, respectively, into the Goddard Earth Observing System, version 5 (GEOS-5) model 18 . It provides hourly, global gridded output of meteorological variables and aerosol optical properties including AOD, AE, and aerosol scattering extinction at 0.625° by 0.5° resolution. The aerosol characteristics are constrained using a wide suite of remote sensing products. For example, AOD at 550 nm is derived from Moderate Resolution Imaging Spectroradiometer (MODIS) measurements on both the Terra and Aqua satellites (Collection 5) 42 of reflectances, solar and instrument geometry, cloud cover, and surface features 18 using a neural network retrieval (NNR) trained using AERONET measurements. A similar approach is used to assimilate Advanced Very High Resolution Radiometer (AVHRR) 43 measurements of radiances, total precipitable water, wind speed, and solar and instrument geometry trained to the MODIS NNR. MISR AOD is assimilated only over bright surfaces 44 , and ground-based AOD measurements from the AERONET 21 are assimilated after 1999. As the density of assimilated aerosol optical properties and meteorological measurements increases greatly after 2000 (refs 18,45 ), the analysis presented here is limited to 2000-2015.
MERRA-2 output includes surface short-and longwave radiation fluxes, with and without clouds, and with and without aerosols, which could be used to estimate aerosol radiative forcing. However these properties are dependent on the radiative transfer model and treatment of aerosol optical properties within the reanalysis model. Thus, herein we only use observable variables that are more closely tied to the assimilated data.
MERRA-2 aerosol properties that are not directly assimilated have been compared to, and found to be in reasonable agreement with, satellite-based radiometric measurements. For example, monthly mean biases relative to the Ozone Monitoring Instrument (OMI) retrieved absorption aerosol optical depth (AAOD) are typically < |0.02| over the NCA regions, and MERRA-2 reproduces the aerosol vertical profile (e.g., height of peak attenuation backscatter) retrieved from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) over the continental U.S. (CONUS) 17 . MERRA-2 has also been evaluated relative to near-surface measurements of PM 2.5 . Again the results indicate a relatively high degree of consistency with independent observations. For most months across the CONUS, MERRA-2 PM 2.5 is within one standard deviation of the in situ measurements, although there is an underestimation of winter PM 2.5 concentrations over the northwest and northeast U.S., potentially due to lack of nitrate aerosols in MERRA-2 (ref. 17 ).
Our analysis of the joint probabilities of AOD, and AE and SSA from MERRA-2 relative to AERONET, indicate good agreement, although MERRA-2 underestimates the dynamic range of AE and SSA ( Figure S1). Such underestimation is common when comparing gridded aerosol datasets that represent area means (~2,500 m 2 for MERRA-2) versus in situ observations such as the pseudo-point measurements from AERONET. MERRA-2 reproduces the observed region-to-region variability in aerosol radiative properties and the MERRA-2 versus AERONET differences tend to be smaller than region-to-region differences ( Figure S1).
Physical variables from MERRA-2 used here within the synoptic-scale meteorological classification have also been extensively evaluated in the previous MERRA release. For example, the mean residual between MERRA and observations is <0.5 hPa for Northern hemisphere surface pressures and ~<1 K for temperature through the depth of the atmosphere relative to radiosonde measurements 46 . Since the original MERRA reanalysis, the GEOS model has been further updated to reduce erroneous trends and discontinuities deriving from breaks in assimilated measurements, and to reduce biases in the water cycle. For all regions in the CONUS, MERRA-2 mean summer precipitation is within ~0.5 mm day −1 (~0.1-0.2 mm day −1 averaged across the CONUS) of surface rain gauge measurements and exhibits an anomaly correlation of ~0.9 for 1980-2011 (ref. 47 ).
The advantages of using the MERRA-2 product for development of aerosol-CIs are manifold. These include use of a consistent data assimilation system for the entire period of record. However, any reanalysis system is subject to inherent uncertainties due either to assimilated variables and/or the model system. For example, an artificial trend exists in Terra radiances assimilated into MERRA-2, which may confound the trend analysis presented herein. Thus trends identified here should be further validated with future MERRA releases in which this trend is corrected and/or with other aerosol reanalysis products as they become available.
Wildfire and anthropogenic emissions. Estimates of wildfire occurrence and spatial extent used herein to diagnose trends in the aerosol-CIs derive from the Global Fire Emissions Database (GFED4) monthly burned area product. GFED4 provides monthly estimates of hectares of burned area on a 0.25° grid derived from the MODIS (Collection 5.1) monthly burned area product 48 .
Annual estimates of anthropogenic emissions of carbon monoxide (CO), NH 3 , NO x , PM 10 , PM 2.5 , SO 2 , and VOCs are also used in attribution of changes in the aerosol-CIs. These estimates are accumulated for all states within each of the NCA regions and derive from the EPA's state level National Emissions Inventory (NEI) 36 . It is noted that there is inherent uncertainty in emissions estimates due to spatiotemporal variability in emission sources, measurement and sampling errors, and the simplification of modeled emissions processes. For example, SO 2 emissions rely on the sulfur content of the combustible material, biogenic emissions vary with environmental conditions, and NH 3 emissions lack wide-spread regulatory restrictions and ambient NH 3 measurements are scarce 49,50 . Additionally, MERRA-2 aerosol speciation depends, in part, on the magnitude of prescribed emissions, which do not evolve (i.e. persistency is assumed) during the later years of the study period 18 . Despite these uncertainties, measurements of species important for secondary aerosol formation, e.g. SO 2 , suggest that trends in emissions are robust 13,51 . Statistical methods used to derive and interpret the aerosol-CIs. The aerosol-CIs we propose quantify the regionally-averaged mean AOD, AE, and SSA; extreme (90 th percentile) AOD; and two geostatistical metrics of spatial autocorrelation and spatial coherence of AOD, AE, and SSA. The regionally averaged mean and P90 values are computed from hourly output that are aggregated in space and time to generate daily mean values for each property that then comprise each CI. While a spatial mean is used here, previous work indicates that spatiotemporal averages are sensitive to averaging methodology 52 , particularly for variables such as AE 53 . The spatial autocorrelation (AOD-I, AE-I, SSA-I) and spatial coherence (AOD-SC, AE-SC, SSA-SC) statistics are computed from the daily mean of the hourly output for each grid cell.
The global spatial autocorrelation for each region and aerosol parameter is computed at the daily timescale and quantified using Moran's I 22 : SCIENTIFIC RePoRTS | (2017) 7:18093 | DOI:10.1038/s41598-017-18402-x where N is the number of grid cells, w ij is the weight for grid cells i and j, X i is the daily mean value (AOD, AE, or SSA) at grid cell i, X is the mean of the daily means for all grid cells, and D ij is the great circle distance between the centroid of grid cell i and j. Values approaching 1 and −1 indicate positive and negative spatial autocorrelation, respectively, while 0 indicates a random spatial field. Significance for rejecting the null hypothesis of no spatial autocorrelation is determined by calculating a z-score for each I: The spatial coherence of each variable in each region is computed using semivariograms which describe the semivariance as a function of separation distance between all grid cell pairs 23 : where N(h) is the number of grid cell pairs that are separated by a great circle distance of h, X i and X j are the daily mean values (AOD, AE, or SSA) at grid cells i and j, respectively, h is a bin range of separation distances, and Q is the set of all grid cells not within three grid cells of the domain border. The empirical semivariogram fit, γ(h), is binned in 100 km increments (i.e. -γ( 1 100 km) includes all grid cell pairs separated by 1-100 km). An exponential fit is used to model γ(h) assuming an exponential decay in correlation with distance and for physical interpretability of the model 53,54 .
where γ′(h) is the exponential model fit; C n is the nugget describing the semivariance at zero spatial lag, resulting from variability at scales below data resolution 54 ; C p is the partial sill, where the sill, C n + C p , is the semivariance as h → ∞; and a is the range or distance at which 95% of the sill is reached, indicating the distance at which two locations are no longer correlated. γ(h) is calculated for each day, and γ′(h)is fit to the mean γ(h) for all days in each season 53 . For the CIs to be tracked through time, a single daily quantity is required. Thus, the daily "scale of SCIENTIFIC RePoRTS | (2017) 7:18093 | DOI:10.1038/s41598-017-18402-x spatial coherence", SC, is herein defined as the minimum h where γ(h)> 0.75 × C p (a s ), where C p (a s ) is the partial sill for that season. While the spatial structure of the AOD and SSA fields is well represented by an exponential model, within the spatial extent of the individual regions AE semivariance tends to increase linearly with distance leading to higher uncertainty in a range determined using the exponential semivariogram model. Temporal trends in the aerosol-CIs are quantified and the significance determined using Kendall's tau-b (τ b ) rank coefficient 55 . τ b is calculated by comparing all pairs of observations, {(t i , X i ), (t j , X j )} where X i and X j are the variable (AOD, AE, SSA) at time t i and t j , respectively: where N is the number of observations. τ b > 0 indicates a positive trend and τ b < 0 indicates a negative trend. The significance of the trend is quantified using z-scores 56 : The slope of the trends, in terms of percentage change per year, is estimated to be the slope of a linear regression fit to the CIs' time series.
It is hypothesized that changes in anthropogenic and natural precursor and primary aerosol emissions will be associated with changes in the aerosol populations. The significance of this association is quantified using the Pearson's r correlation coefficient.
Prior research indicates that synoptic meteorological conditions are also a key control of aerosol concentrations 29,30 . Thus, PCA is used to derive a daily synoptic classification for all days in the study period and investigate the interaction between synoptic conditions and aerosol properties, and to determine the impact of meteorology on the CIs trends. Predictors used in the PCA are air temperature and specific humidity at 700 hPa plus 500 hPa geopotential heights from MERRA-2. The number of PCs to retain for each region was determined using a scree test 57 and the retained factors are rotated using a Varimax rotation 58 . Between six and nine components (i.e. unique synoptic types) were retained for each of the eight NCA regions. The PC scores for each day (i.e. similarity to the major modes of variability as characterized by the PCs) are used to track changes in the frequency of each synoptic type (i.e. counts of days with highest similarity to each of the modes) and the intensity of the types (i.e. the magnitude of the scores for each PC). The mean anomaly of each aerosol-CI on all days classified by each synoptic type, calculated relative to the mean aerosol-CI computed for all days, is used to illustrate the importance of meteorological conditions at the synoptic (regional) scale in determining aerosol properties. Data availability. MERRA-2 data is available from the Goddard Earth Science Data and Information Services Center (https://disc.sci.gsfc.nasa.gov/), AERONET data is available from https://aeronet.gsfc.nasa.gov/, GFED4 is available from http://www.globalfiredata.org/, and NEI is available from https://www.epa.gov/sites/ production/files/2016-12/state_tier1_90-16.xls.