Multidecadal trend analysis of in situ aerosol radiative properties around the world

In order to assess the evolution of aerosol parameters affecting climate change, a long-term trend analysis of aerosol optical properties was performed on time series from 52 stations situated across five continents. The time series of measured scattering, backscattering and absorption coefficients as well as the derived single scattering albedo, backscattering fraction, scattering and absorption Ångström exponents covered at least 10 years and up to 40 years for some stations. The non-parametric seasonal Mann–Kendall (MK) statistical test associated with several pre-whitening methods and with Sen's slope was used as the main trend analysis method. Comparisons with general least mean square associated with autoregressive bootstrap (GLS/ARB) and with standard least mean square analysis (LMS) enabled confirmation of the detected MK statistically significant trends and the assessment of advantages and limitations of each method. Currently, scattering and backscattering coefficient trends are mostly decreasing in Europe and North America and are not statistically significant in Asia, while polar stations exhibit a mix of increasing and decreasing trends. A few increasing trends are also found at some stations in North America and Australia. Absorption coefficient time series also exhibit primarily decreasing trends. For single scattering albedo, 52 % of the sites exhibit statistically significant positive trends, mostly in Asia, eastern/northern Europe and the Arctic, 22 % of sites exhibit statistically significant negative trends, mostly in central Europe and central North America, while the remaining 26 % of sites have trends which are not statistically significant. In addition to evaluating trends for the overall time series, the evolution of the trends in sequential 10-year segments was also analyzed. For scattering and backscattering, statistically significant increasing 10-year trends are primarily found for earlier periods (10-year trends ending in 2010–2015) for polar stations and Mauna Loa. For most of the stations, the present-day statistically significant decreasing 10-year trends of the single scattering albedo were preceded by not statistically significant and statistically significant increasing 10-year trends. The effect of air pollution abatement policies in continental North America is very obvious in the 10-year trends of the scattering coefficient – there is a shift to statistically significant negative trends in 2009–2012 for all stations in the eastern and central USA. This long-term trend analysis of aerosol radiative properties with a broad spatial coverage provides insight into potential aerosol effects on climate changes.

1 Introduction

Climate change has been considered a premier global problem in the scientific community for decades. Thirty years ago, the community organized to produce the first Intergovernmental Panel on Climate Change (IPCC) report (IPCC, 1990) about the state of scientific, technical and socio-economic knowledge on climate change, its impacts and future risks, and options for reducing the rate at which climate change was taking place. Aerosols have been recognized as an important active climate forcing agent since the 1970s and, in the last IPCC report (IPCC, 2013), the impact of aerosols on the atmosphere was still considered to be one of the most significant and uncertain aspects of climate change projections and, for the first time, decadal trend analysis of in situ aerosol optical properties around the world was reported.

Aerosol optical properties are the relevant parameters that determine the radiative forcing of particulate matter. While some of these optical properties are currently measured by satellite (Choi et al., 2019), airborne and ground-based remote sensing (REM) technologies (https://aeronet.gsfc.nasa.gov/, last access: 20 July 2020, https://www.earlinet.org/, last access: 20 July 2020), the ground-based, in situ measurements represent some of the longest time series, allowing assessment of the long-term time evolution of aerosol radiative properties in the lower troposphere.

The first in situ measurement network began in the mid 1970s at several remote locations (Bodhaine et al., 1995). Through national and international programs and/or on an individual organization's initiatives, the number of stations with systematic aerosol monitoring activities in regional background locations has continued to increase since the 1990s. As of 2017 absorption has been measured for at least 1 year (yr) at 50 sites, for 5 years at 37 sites and for 10 years at 20 sites, while scattering has been measured for at least 1 year at 56 sites, for 5 years at 45 sites and for 10 years at 30 sites. The companion paper (Laj et al., 2020) provides a historical view and a complete description of the present networks for aerosol measurements. The longest datasets cover up to 40 years of measurements, BRW (40 years), SPO (40 years), and MLO (31 years) (see Table 1 for station acronyms), whereas some stations with long time series recently closed or moved (THD, SGP, MUK, CPT). The spatial and temporal variability of aerosol properties is extremely high due to the short lifetime of aerosol particles (on the order of days to weeks), the wide variety of sources, as well as the chemical and microphysical processing occurring in the atmosphere; a dense network of stations is consequently required to obtain a global view of aerosol changes. The growing number of stations with long-term (> 10 years) time series of aerosol particle optical properties – 24 in 2010 (Collaud Coen et al., 2013, hereafter referred to as CC2013) and now 52 in 2016–2018 – is a positive factor. Detracting from that growth is the continued lack of sites in South America, Africa, Oceania and Asia.

Table 1 List of observatories included in this study, arranged alphabetically by Global Atmospheric Watch (GAW) acronyms, including their names, countries, coordinates and elevation, site environmental characteristic (geographical category and footprint), size cut, type of nephelometer and absorption filter photometer deployed, time period used, and nephelometer RH percentiles.

¹ Geographical category: Mountain: Mt, Polar: P, Continental: Con, Coastal: coast. Footprint: Rural background: RB, Forest: F, Desert: DE, (Sub-)Urban: U, Pristine: P Mixed: Mix.
² the mention of two size cuts separated by “–” corresponds to a modification of inlet during the tim series, whereas the “&” corresponds to measurements at two size cuts.
³ T: TSI nephelometer, O: Optec nephelometer, R: Radiance Research nephelometer; E3: Ecotech nephelometer Aurora 3000, E4: Ecotech nephelometer Aurora 4000.
⁴ AE16(C_ref= 1.8) 120∕ AE22 (C_ref= 1.8) ∕ AE31 (C_ref= 3.5) ∕ AE33 (C_ref= 3.5): Aethalometer, P1: one-wavelength or P3: three-wavelength PSAP, M: MAAP, C: NOAA CLAP, ET: ES95L Thermo 5012$==$M.

Long-term measurements are the only possible approach for detecting change in atmospheric composition resulting from either changes in natural or anthropogenic emissions and/or changes in atmospheric processes and sinks. However, detecting long-term trends of aerosol optical properties remains a challenge, due to their high natural variability, uncertainties caused by changes and biases in measurement methodology, the ill-defined statistical distribution of the parameters, the presence of high autocorrelation in aerosol parameters, as well as the occasional issues regarding traceability of historic operating procedures. Trend analysis can only be performed on time series without breakpoints or on homogenized time series that account for changes in measurement conditions (e.g., relocations, instrument calibration/repair/upgrades, inlet changes) (CC2013). Once homogenized datasets are available, appropriate techniques must be used to identify potential trends. The trend analysis methodology must take into account the non-normal distribution of most aerosol parameters, the high autocorrelation of the parameters, and the presence of gaps and negatives in the datasets.

In this current analysis, a considerable effort was made to detect time series breakpoints, to find explanations for them in the logbooks and station history and, if possible, to correct or homogenize the time series. These homogenized time series were then subjected to an array of statistical tests to identify trends. These tests include (1) the non-parametric seasonal Mann–Kendall test (hereafter referred to as the MK test) associated with Sen's slope. The applied MK test is however applied with a new pre-whitening method (Collaud Coen et al., 2020a), (2) a generalized least squares (GLS) method associated with a Monte Carlo bootstrap algorithm and (3) the least mean squares fit (LMS). While the MK test with pre-whitening was considered the most robust method, the other tests were included to allow a comparison between various simple and frequently used methods.

The first long-term trend analyses of aerosol optical properties, number concentration and particle size distribution (CC2013; Asmi et al., 2013) covered 2001–2010 as the shortest period and longer periods if data were available. The main observations were (1) a general statistically significant (ss) – at 95 % confidence level – decrease in number concentration, scattering and absorption coefficients in North America, (2) a ss decrease in number concentration in northeastern Europe, (3) no ss trends in central Europe for any of the parameters and (4) no ss scattering coefficient trends but increasing 10-year absorption coefficient and number concentration trends in the polar regions. These trends were related to the decrease in anthropogenic primary aerosol emissions and in precursors of secondary aerosol formation. The high-altitude station Mauna Loa (MLO) in the Pacific was unique in exhibiting increasing optical property trends that were mostly attributed to long-range transport from Asia. The results in CC2013 are in line with the 1996–2013 trend analysis at the BND and SGP stations in North America (Sherman et al., 2015) showing a decreasing scattering coefficient and a sub-micron scattering fraction and increasing backscattering fraction. More recently, Pandolfi et al. (2018) presented the long-term trends of in situ surface aerosol particle optical properties (scattering) measured in Europe until 2015. The ss decreasing trends of aerosol particle scattering observed in Europe at around 40 % of the stations (mostly in Nordic and Baltic countries and southwestern Europe) were attributed to the implementation of continental to local emission mitigation strategies. Pandolfi et al. (2018) also reported that the scattering Ångström exponent decreased at around 20 % for the European stations included in their study (at remote Nordic and Baltic locations and at two mountain sites in central and eastern Europe), whereas an increase was observed at 15 % of the stations (one urban site in southwestern Europe and one in central Europe). In the same study, the backscattering fraction was observed to increase. Trends in horizontal visibility synoptic observations over 1929–2013 from 4000 stations over the USA, Europe and Asia (Li et al., 2016) generally agreed with extinction coefficient trends, with a significant decrease in all regions but with different evolutions of the trends. Hand et al. (2014) also found a significant drop of the ambient light extinction coefficient at all IMPROVE (Interagency Monitoring of Protected Visual Environment, http://vista.cira.colostate.edu/Improve/, last access: 20 July 2020) stations over the 1990s through 2011, with a larger decrease in the eastern USA. To our knowledge, no further trend analyses of surface in situ aerosol optical properties involving a network of stations or several stations have been published up to now.

This study is part of the SARGAN (in-Situ AeRosol GAW observing Network) initiative (see the companion paper by Laj et al., 2020) with the objective of supporting a global aerosol monitoring network to become a GCOS (Global Climate Observing System) associated network. This trend analysis is intended to answer the following questions.

1. Are there homogeneous long-term trends of in situ aerosol optical properties over the covered regions of the world? Do they differ as a function of the length of the data series? How do the trends evolve with time?

2. Are there regional similarities or differences in the observed trends among stations? Are there similarities or differences in trends among aerosol parameters at regional and continental scales?

3. How do the observed optical property trends compare with trends in other aerosol and gaseous properties reported in the literature?

The results of this study provide the best representation of change in surface aerosol optical properties considering the available in situ aerosol optical property datasets and highlight the possible side-effects of air pollution control policies on radiative forcing.

Figure 1 Map of stations with their GAW acronyms.

2 Experimental

2.1 Measurement sites

The long-term trend analysis presented in this study analyzes in situ aerosol time series from 52 observatories worldwide shown in Fig. 1 with site information listed in Table 1. The network, which is a subset of the station network described in Laj et al. (2020), comprises 16 stations in Europe, 21 in North America, 5 in Asia, 2 in Africa, 6 in the polar regions and 2 in the southwestern Pacific. The stations included in this study are primarily located in rural or remote areas and are expected to exhibit regional- to large-scale representativeness (e.g., Wang et al., 2018). Apart from MUK, all the stations are regional or global GAW (Global Atmospheric Watch, https://gawsis.meteoswiss.ch/GAWSIS//index.html#/, last access: 20 July 2020) sites or IMPROVE stations. The GAW aerosol data are archived at and available from the World Data Centre for Aerosol (WDCA, http://www.gaw-wdca.org, last access: 20 July 2020) located at the Norwegian Institute for Air Research (NILU). The WDCA data repository is the EBAS database (http://ebas.nilu.no, last access: 20 July 2020), an e-infrastructure shared with other frameworks targeting atmospheric aerosol properties, such as the Co-operative Programme for Monitoring and Evaluation of the Long-range Transmission of Air pollutants in Europe (EMEP) and the European Aerosols, Clouds, and Trace gases Research InfraStructure Network (ACTRIS). The IMPROVE data are available from the IMPROVE website (http://vista.cira.colostate.edu/Improve/data-page/, last access: 20 July 2020) and from the WDCA. To ensure that the long-term trend analysis was performed on homogeneous time series, a substantial effort of quality control, rupture detection and homogenization (see Sect. 2.4) was performed in close collaboration with each station's operator on the data. As has been noted in previous papers, it is critical to have outside review of data to improve the quality of long-term time series (CC2013; Asmi et al., 2013). The final time series used in this analysis are available from the following DOI: https://doi.org/10.21336/c4dy-yw57.

The stations' environments were classified into four types (continental, coastal, mountain, or polar) that are represented by 22, 8, 16 and 7 time series, respectively. The type of measured aerosol at each site is further characterized by their footprints comprising six types (rural background, forest, desert, (sub)-urban, pristine and mixed). While the environments of Europe, North America and the polar regions are fairly well represented, the number of long-term stations in the rest of the world is currently quite low, resulting in a lack of information from the largest deserts (e.g., Sahara, Gobi, Australian, Arabian, Atacama), from many mountain ranges (e.g., Himalaya, Andes, Southern Great Escarpment, Great Dividing Range, Urals) and from whole continents (South America (no site), Africa (one island in the Atlantic and one coastal site), and Australia (one coastal site)). Some stations from these underrepresented areas currently have 4 to 7 years of measurements available and will potentially be used for trend analyses in the future (see Table 3 in Laj et al., 2020).

Sites were chosen based on the following criteria: (1) availability of at least 10 years of continuous data (two sites with 9 years and one site with 8 years of data for at least one parameter have also been included to improve spatial coverage (CPT, EGB and GSN, respectively)); (2) continuous measurements without ruptures in the aerosol light scattering and/or absorption measurement; (3) submission of quality-assured data to the WMO WDCA data repository; and (4) responsiveness of site operators to questions concerning data quality and homogeneity.

The longest time series with 40 years of measurements are the Arctic and Antarctic stations of BRW and SPO, followed by the high-altitude MLO station (31 years). During the 1990s NOAA began extending their network (Andrews et al., 2019), the IMPROVE network installed numerous stations in the USA (Malm et al., 1994), and the first long-term measurements in Europe, JFJ (Bukowiecki et al., 2016) and HPB, began in 1995. To have the largest representativity and to minimize the number of stations with less than 10 years of measurement, the current long-term trends were computed from time series ending in 2016, 2017 or 2018 (whichever year was most recently available). To obtain an overview of the long-term trend evolution in the past 40 years, all stations with at least 10 years of measurements were considered (see results in Sect. 3.2).

2.2 Instruments

The relevant instruments operating at each site are listed in Table 1 and further instrument details are given in the Supplement (Table S1). Some particular instrumental features that could influence the trend analysis or comparison between stations are briefly discussed below.

Nephelometers measure aerosol light scattering over a truncated angular range (Müller et al., 2009, and references therein), leading to non-idealities often called “truncation error”. The truncation adjustment accounts for scattering over the angles outside the measurement range and non-ideality of the light source. All TSI nephelometer scattering and backscattering sets were adjusted for truncation and instrument non-idealities using the Anderson and Ogren (1998) correction. Thus, for times when enhanced amounts of large diameter (D_p>1 µm) particles are present, the measured scattering will be lower than true scattering by a substantial amount since the truncation correction increases with particle size (Anderson and Ogren, 1998; Molenar et al., 1997). The Radiance Research nephelometer has similar truncation characteristics to the TSI nephelometer (Müller et al., 2009). The Optec nephelometer measures over a wider angular range (Molenar, 1997) than the other nephelometers and, like the Radiance Research measurements, the scattering has not been corrected for truncation in this study. The Optec nephelometers measure at ambient conditions with no size cut (they are open-air instruments) so they can sample the very large particles present due to both hygroscopic growth at high humidities and/or the occurrence of precipitation, fog, dust, pollen, etc. The Ecotech nephelometers have a similar angular range to the TSI nephelometers, and the measurements are corrected for truncation errors using the Müller correction (Müller et al., 2011b), adapted from the Anderson and Ogren correction.

For better comparability of aerosol properties amongst sites and to minimize the confounding effects of water associated with the aerosol, GAW recommends drying the sample air to RH < 40 % (WMO/GAW report 227, 2016). While most of the nephelometer scattering time series are accompanied by sample RH measurements, this was not the case for all stations and for the entire measurement period. The calculated RH trends are therefore not always complete. Many breakpoints were detected in sample RH data and exchanges with the individual station operators revealed that humidity sensors often suffer from artifacts, offsets, and modifications that were not considered problematic. These sensor problems were often not resolved due to the secondary status of this housekeeping diagnostic, leading to problematic time series. Nonetheless, apart from the IMPROVE network, the majority of nephelometers appeared to have sampled at RH < 40 %. The IMPROVE scattering measurements were analyzed at the measurement conditions with some constraints on acceptable scattering values, although the IMPROVE network recommends screening the data when RH > 90 % (Prenni et al., 2019). For this study and according to CC2013, the IMPROVE scattering coefficient was restricted to σ_sp values lower than 500 M m⁻¹ for stations in the eastern USA (ACA, GSM, MCN and SHN) and lower than 100 M m⁻¹ for stations in the western USA to minimize the influence of rain, fog, snow and ice. These screening constraints minimized the issues associated with high RH but do not correspond to a screening based on RH.

Measurement of the absorption coefficient was always performed by some type of filter-based photometer but relied on a variety of instruments. These instruments include Multi-Angle Absorption Photometers (MAAPs), Particle Soot Absorption Photometers (PSAPs) and Continuous Light Absorption Photometers (CLAPs), as well as various models of the Aethalometer (AE16, AE21, AE31 and AE33). All these instruments suffer from various artifacts, from which the loading effect can influence the wavelength dependence. However, the largest uncertainty in filter-based photometer measurements lies in the effect of the multiple scattering of light into the filter matrix, leading to over-prediction of absorption aerosol (e.g., Bond et al., 1999; Lack et al., 2008; Müller et al., 2011a; Collaud Coen et al., 2010; Bernardoni et al., 2019). This artifact is roughly corrected by the multiple scattering constant C_ref and is probably largest for the Aethalometer and smallest for the MAAP.

The ACTRIS community has suggested that Level 2 AE31 data submitted to EBAS utilize a multiple scattering constant C_ref=3.5; most of the analyzed AE31 time series were corrected with this new rule. The AE33 adds a simultaneous measurement of the light transmission through a second filter spot sampling the same air at a different flow rate associated with a real-time compensation algorithm. This two-spot technique allows for correction of the filter loading artifact. This improvement, however, has no effect on the largest artifact (multiple scattering artifact) and, as of yet, there is no agreed upon correction for the AE33 by the aerosol community. Previous AE models used a white light diode (AE10 and AE16) and a C_ref=1.6 is usually applied. At FKL, the AE21 used a C_ref=1.8 and the AE33 a C_ref=3.0. The various versions of the Aethalometer require then different corrections, whereas the real C_ref value depends on the filter and on the aerosol type. For background rural aerosol, the real C_ref value is between 2.5 and 4.5 (Collaud Coen et al., 2010; Bernardoni et al., 2019), the Asian plume has a relatively high C_ref between 4 and 5.5 (Kim et al., 2018), in the Arctic C_ref is suggested to be 3.45 (Backman et al., 2017), whereas pure mineral dust leads to a lower C_ref of 1.75–2.56 (Di Biagio et al., 2017).

The MAAP measures not only the light transmission through the filter, but also the light backscattered at two different angles. This design takes into account the scattering and multiple scattering artifacts (see Collaud Coen et al., 2010), which are two of the most significant artifacts for filter-based absorption photometers so that no correction is needed (C_ref=1). The MAAP measured absorption coefficient is consequently more reliable.

The CLAP was developed by NOAA as a replacement for the PSAP (Ogren et al., 2017). The CLAP was designed to have the same optical characteristics as the PSAP so that either the Bond et al. (1999) correction along with the Ogren (2010) update for wavelength and spot size correction or the Virkkula et al. (2005, 2010) corrections can be applied to account for scattering artifacts at multiple wavelengths as well as other instrument non-idealities (e.g., filter-loading artifacts, variability in spot size and flow calibrations). These correction algorithms rely on co-located scattering measurements from a nephelometer and may have issues in the presence of large, primarily scattering aerosol such as sea salt or dust (e.g., Bond et al., 1999) and also may not work well when organic aerosol is abundant (e.g., Lack et al., 2008).

The differences in instrumentation, measurement conditions and post-processing data treatment do not allow the absolute values of aerosol optical parameters for all sites to be compared; however, because there was consistency of data treatment for each individual time series, the trends across the different sites can be compared.

2.3 Aerosol optical properties

The data used in this paper consist of hourly-averaged, quality-checked, spectral light scattering (σ_sp), backscattering (σ_bsp) and absorption (σ_ap) measurements. The quality checks correspond to the Level 2 requirements of EBAS (Laj et al., 2020). After further visual quality control by the authors, the hourly data were aggregated into daily medians with the requirement that at least 25 % of the daily data be valid. The median was chosen to minimize the effect of extreme values on the average since the measured parameters are strongly not normally distributed and most of the calculated parameters also do not follow a normal distribution. Such a low requirement for data coverage was chosen since six hourly measurements a day corresponds to half of the potential data coverage at many of the NOAA stations, where the operation mode consists of alternating between PM₁ and PM₁₀ size cutoff on a sub-hourly basis (Andrews et al., 2019).

All the nephelometers and the multi-wavelength absorption photometers measure at a green wavelength (∼ 525–550 nm), which is the channel for which the parameters are reported. For the AE31 and AE33 models, the 520 nm channel was chosen. At several sites, the light absorption was measured by white light (∼ 840–880 nm) Aethalometers (AE16), two-channel Aethalometers (AE21) using 370 and 880 nm or by MAAPs (Multi-Angle Absorption Photometers) at 637 nm (Müller et al., 2011a), requiring the use of another wavelength, typically a red wavelength. In some cases, the blue or red wavelength was preferred due to inhomogeneities or gaps in the green data. Since the trend analysis is not sensitive to the multiplication by a constant, the data series used to determine scattering and absorption trends were not adjusted to 550 nm.

In addition to the measured parameters, the following parameters were computed when the appropriate measurements were available:

• backscatter fraction, $b={\mathit{\sigma }}_{\text{bsp}}/{\mathit{\sigma }}_{\text{sp}}$,

• scattering Ångström exponent,
  ${\mathit{å}}_{\text{sp}}=-\ln ({\mathit{\sigma }}_{\text{sp,1}}/{\mathit{\sigma }}_{\text{sp,2}})/\ln ({\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}})$,

• absorption Ångström exponent,
  ${\mathit{å}}_{\text{ap}}=-\ln ({\mathit{\sigma }}_{\text{ap,1}}/{\mathit{\sigma }}_{\text{ap,2}})/\ln ({\mathit{\lambda }}_{\mathrm{1}}/{\mathit{\lambda }}_{\mathrm{2}})$, or by a linear fit between the logarithm of the seven absorption coefficients as a function of the logarithm of the seven wavelengths of the Aethalometers (AE31 and AE33), and

• single scattering albedo, ${\mathit{\omega }}_{\mathrm{0}}={\mathit{\sigma }}_{\text{sp}}/({\mathit{\sigma }}_{\text{sp}}+{\mathit{\sigma }}_{\text{ap}})$,

where σ_sp,i is the scattering coefficient at wavelength i, λ_i is the wavelength i, σ_bsp is the hemispheric backscattering coefficient, and σ_ap is the absorption coefficient.

å_sp and å_ap were usually computed from the blue (∼ 450 nm) and green wavelengths, because the red channel of the nephelometers was frequently less stable and more prone to rupture in the time series due to calibrations or instrument changes. However, in some cases, other wavelength pairs were used to utilize the longest time series. å_ap computed from AE31 and AE33 is always more homogeneous if fitted on the seven wavelengths, so that the fitted å_ap was always chosen for these two instruments.

The single scattering albedo was computed from σ_sp and σ_ap after σ_ap was adjusted to match the nephelometer green wavelength with an assumed absorption Ångström exponent of one (i.e., 1∕λ dependence). In order to maintain similar data treatment for absorption instruments with single or multiple wavelengths, the measured absorption Ångström exponents were not used for the wavelength adjustment for the ω₀ calculation.

It should be recalled that all parameters calculated using ratios of the σ_sp, σ_bsp and/or σ_ap may have higher uncertainties for two reasons: (1) the ratio of two similar values has a larger uncertainty than the σ_sp, σ_bsp or σ_ap uncertainties and (2) the σ_sp difference between the wavelengths depends on the nephelometer calibration that is performed independently for each wavelength. These uncertainties are particularly enhanced for clean locations with low aerosol loading.

2.4 Discontinuities, data consistency and homogenization

Long-term climate analyses require homogeneous time series to be accurate. A homogeneous climate time series is defined as one where variations are caused only by variations in weather and climate (Conrad and Pollak, 1950) and in emissions of aerosol particles and their precursor gases. Long-term climatological time series can be affected by a number of non-climatic factors called breakpoints (e.g., relocation, instrument upgrades, inlet changes, calibrations, nearby pollution sources) that mask the real climate variations. The breakpoints can be detected either by subjective visual inspection or by objective statistical methods (Peterson et al., 1998; Beaulieu et al., 2007) and must correspond to an event recorded in logbooks describing the station/instrumental history. Many statistical methods are only suitable for normally distributed data and cannot therefore be applied to aerosol optical property measurement without data transformation (Lindau and Venema, 2018). Moreover, they are often applied not only to the data, but also to ratios or differences between various time series that are not systematically available at all the measuring sites of this study.

Visual inspection was used to detect breakpoints and to assess the validity of the time series to be used for climatic trend analysis. For this study, each measured and calculated (see Sect. 2.3) parameter at all wavelengths, as well as all the possible ratios between measured parameters (including the number concentration if available), at each station were visually inspected in linear and logarithmic time series plots. The treatment of minimum and maximum values, outliers and negatives along with the consistency of seasonal cycles were looked at closely when inspecting the time series plots. In addition, the data owners responded to a questionnaire about potential breakpoints, providing metadata that could be used to confirm/dismiss possible breakpoints or to accurately locate them. The identified breakpoints were discussed with the data owners, leading to corrections, homogenization, invalidations or splitting of the time series into two parts. In one case (absorption data from SUM measured by AE16 and CLAP), the two time series were homogenized by multiplying the AE16 data by the median of the ratio between both datasets during the 10.5 months of simultaneous measurements. Only datasets considered homogeneous by the authors and the data owners were analyzed in this study.

In the older networks, several modifications likely lead to inhomogeneities that occurred at sites in the network around the same time. Some of these include the following.

1. Two of the longest running NOAA stations changed their TSP (total suspended particle) inlets for PM₁₀ size cuts in the middle of the multi-decade time series (MLO: 2000, BRW: 1997). Some other stations outside the NOAA network also modified the measurement size cuts over their long-term measurement period. Usually this change in size cut (TSP to PM₁₀) did not generate a breakpoint for aerosol optical properties so that the time series could be considered homogeneous. A differentiation between periods of sampling inside or outside of clouds was not made, even though TSP and PM₁₀ could respond differently in these situations. In contrast, the modification of TSP or PM₁₀ size cuts to PM_2.5 or PM₁ cutoffs usually led to visible breakpoints. PAL is the only station where changes between PM₁₀, PM₅ and PM_2.5 did not induce a visually obvious breakpoint, likely due to the minimal presence of supermicron particles at this site.

2. The NOAA stations used the single green wavelength PSAP until the years 2005–2007, when they replaced them with a three-wavelength (3w) PSAP (see Table S1). This instrumental change usually did not induce a visually obvious breakpoint.

3. A further instrument change for the absorption coefficient at NOAA sites occurred in 2013–2015 through the introduction of the 3w CLAP. The 3w PSAP to 3w CLAP change usually induced no breakpoint in the green absorption coefficient. The red channel sometimes exhibited a visible breakpoint (APP and BND), resulting in breakpoints in the absorption Ångström exponent. In those cases, calculation of the absorption Ångström exponent with the blue and green channels was preferred.

4. The long time series from MLO and JFJ were subject to the removal of negative values during the first years of measurements until 2000 and 1999, respectively. The raw data prior to these years were not archived by the data providers for either site. This change in minimal values does not seem to produce a clear breakpoint in the sense that the computed trends were not affected strongly enough to modify the climatic trends.

To compare long-term trends between stations from various networks, instruments and operators, instrumentation, measurement conditions and data treatment consistency is critical, but some lenience amongst stations was deemed acceptable. Specifically, some discretion was allowed, including whether the datasets had the same corrections applied (e.g., truncation or not), how the sites dealt with sample RH and very low aerosol amounts, and inlet size cuts. Table 1 includes columns indicating information about the size cuts and RH conditions at the various sites. No screening or analysis as a function of cloud amount/clear-sky conditions was done since these criteria/flagging were not available at all stations. Below, the impacts of sample RH, size cut and general instrument conditions and corrections on trend evaluation are briefly discussed.

1. Humidity: one important factor affecting all aerosol measurements is the relative humidity (RH) at which the measurements are made. For σ_sp, measurements at controlled RH enable minimization of the confounding effects of aerosol hygroscopic growth, resulting in increases in the amount of scattering aerosol (Nessler et al., 2005; Fierz-Schmidhauser et al., 2010; Burgos et al., 2019). The disadvantage of making measurements at low RH is that aerosol hygroscopic properties must be measured or assumed in order to adjust the aerosol optical properties to ambient conditions. As noted above (see Sect. 2.2), within the GAW program, recommendations have been given to measure σ_sp at low (RH < 40 %) humidities. Apart from the IMPROVE and CPR nephelometers, the instruments typically operated at RH < 50 %, with only six stations having a RH 95th percentile value larger than 50 % (AMY, CMN, EGB, GSN, IPR and SGP) but with a median clearly much lower than 50 %. In contrast, the IMPROVE network instruments measure at near-ambient conditions (Malm et al., 1996). The scattering restriction method (see Sect. 2.2) was chosen in order to maintain the highest data coverage – simply removing scattering values associated with RH > 50 % from the ambient IMPROVE dataset would have eliminated most of the summertime measurements, particularly for the eastern USA locations. For all stations with some contribution of scattering made at RH values larger than 50 %, the dry scattering and backscattering coefficients were calculated by removing values corresponding to hourly RH median > 50 %.

   Ensuring a low humidity in the nephelometer reduces but does not suppress the potential influence of the hygroscopic growth on nephelometer measurements (Zieger et al., 2013). Therefore, if RH data were available, the RH long-term trends were also computed and their potential effect on the trend of σ_sp, σ_bsp, b and å_sp was evaluated (see Sect. 4.1).

   The filter-based absorption photometers are also sensitive to rapid RH changes (e.g., Anderson et al., 2003), but daily absorption averages are usually not biased by such rapid fluctuations (Bernardoni et al., 2019). Very high sample RH could lead to higher uncertainties, but absorption measurements at GAW stations are usually connected to inlets with some sort of conditioning intended to reduce sample RH (e.g., diffusion or membrane dryers, dilution with dry air and, in some cases, heating). Additionally, CLAPs are gently heated to ∼ 37 ^∘C to minimize RH effects. In this study, stations with high sample RH in the nephelometer sample (Table 1) are also the most likely to have issues with high sample RH in the collocated absorption photometer.

2. Size cut: as described in Table 1, the size cuts differ amongst the stations, but most of the sites measure TSP or PM₁₀. The GAW program generally recommends a PM₁₀ size cut, except for stations in extreme environments (clouds, etc.), where a whole air inlet is recommended (WMO/GAW report 227, 2016; GAW/WCCAP recommendations https://www.wmo-gaw-wcc-aerosol-physics.org/files/WCCAP-recommendation-for-aerosol-inlets-and-sampling-tubes.pdf, last access: 20 July 2020). Many stations in the NOAA Federated Aerosol Network measure at a second size cut (PM₁) as well. PAY and SUM are the only stations that have no measurement of coarse-mode aerosol, with only a PM_2.5 inlet. As reported previously, the amount of aerosol particles larger than 10 µm is usually sufficiently low to enable consideration of TSP and PM₁₀ results as being in the same category. Moreover, the trend results of PM₁₀ and PM₁ sampling are found to be quite similar for all stations with both size cuts, so that the results of TSP/PM₁₀ size cut will be presented in this study and, if not specified, PM₁ results can be assumed to be similar to those of the larger size cut (PM₁₀ or TSP).

3. Absorption filter photometer artifacts: the first main point to consider is that all filter-based absorption photometers suffer from various measurement artifacts and that continuous reference measurements to assess the absolute σ_ap values are not available at long-term monitoring sites. If the variability and the long-term trends of absorption coefficients are to be analyzed with high confidence, the σ_ap absolute value is necessary to compute the ω₀. As stated in Sect. 2.2, the real C_ref values can potentially vary by a factor of 4 (1.5 to 5.5). Using an erroneous C_ref value can influence the magnitude of the ω₀ trends. Similarly, an applied correction depending on the wavelengths can affect the absorption Ångström exponent calculation and its trends. Both ω₀ and å_ap long-term trends therefore must be interpreted with greater care.

4. Nephelometer truncation correction artifacts: as explained in Sect. 2.2, the various types of nephelometers measure at different truncated angular ranges that were corrected by several algorithms or even not corrected. The absence of truncation correction leads to lower scattering and backscattering coefficients than the true values and the correction algorithm effects are known to increase with particle size. The most important requirement that was verified for this trend analysis is the coherent treatment of nephelometer data for each time series. The bias leading to a higher contribution of Aitken and accumulation modes than the coarse mode is difficult to estimate, but the minimal differences in PM₁ and PM₁₀ results (see Sect. 4.2) suggest this artifact is small. The effect of the humidity on the nephelometer measurements is regarded as the most significant artifact.

Finally, in order to minimize the potential artifacts in the determination of the long-term trends in the case of large seasonal variability (de Jong and de Bruin, 2012), only full start and end years of the time series, that is, without gaps in the data, were considered. For some stations, we did allow gaps of up to 4–6 weeks without measurements after checking that the removal of the whole year led to similar trend results.

The differences in instrumentation, measurement conditions, and post-processing data treatment do not allow the absolute values for all sites to be compared; however, because there was consistency of data treatment for individual sites, the trends can be compared.

2.5 Trend analyses

The aerosol extensive parameters (σ_sp, σ_bsp and σ_abs) are not normally distributed and they exhibit varying degrees of autocorrelation. They can be represented approximately by a lognormal distribution but are usually better fitted by a distribution in the Johnson distribution family (Johnson, 1949). The intensive parameters (b, å_sp, å_ap and ω₀) also exhibit distributions that differ to varying degrees from the normal distribution. We chose, therefore, to rely mostly on the non-parametric seasonal Mann–Kendall (MK) test associated with Sen's slope. The MK test does not require normally distributed data. Additionally, as described under Sect. 2.5.1, the MK test was adapted to correctly handle autocorrelated datasets. To allow a comparison with other studies, the trends were also computed with the generalized least squares analysis associated with the autoregressive or block bootstrap confidence intervals (GLS) and the least-mean square (LMS) fit applied to the data logarithms.

2.5.1 Mann–Kendall test and Sen's slope estimator

This non-parametric method based on rank (Gilbert, 1987; Sirois, 1998) is the most appropriate test to compute optical property trends because it can be applied regardless of missing values, statistical distribution and presence of negatives or below detection limit values in the dataset. The MK test determines whether a monotonic increasing or decreasing long-term trend exists; the slope and the confidence limits are then computed by Sen's slope estimator that is based on the median of the slopes calculated from all possible data pairs. For this study, the MK test was applied on daily medians.

The MK test is designed for serially independent data and is, consequently, influenced by autocorrelation in the time series, leading to inflated type 1 error; that is, there is increased probability of rejecting the no-trend hypothesis (i.e., a false positive). Several correction schemes for the MK test were proposed to correctly handle autocorrelated datasets, and the problems induced by autocorrelation and its various corrections have been clearly described (Wang and Swail, 2001; Yue et al., 2002; Zhang and Zwiers, 2004; Bayazit and Önöz, 2007; Blain, 2013; Wang et al., 2015). A new method (Collaud Coen et al., 2020) has been used for this study that tends to minimize the type 1 and 2 errors (type 2 error is non-rejection of a false null hypothesis, i.e., a false negative) as well as issues with the modification of the slope due to pre-whitening procedures by the application of three pre-whitening (PW) methods. The standard pre-whitening by removing the first lag autocorrelation (von Storch, 1995) has a very low type 1 error but also a low test power, whereas the so-called trend-free pre-whitening procedure published by Yue et al. (2002) (called TFPW-Y in Collaud Coen et al., 2020a) restores the test power at the expense of the type 1 error. Both these pre-whitening procedures were applied prior to the MK test to assess the statistical significance of the trend. A trend was then considered to be ss only if both PW and TFPW-Y were ss at the 95 % confidence level or if PW is ss but not TFPW-Y (false negative). Among the trends of all parameters at all stations calculated for this paper, none was ss for the PW but not for the TFPW-Y, meaning that the PW procedure was always powerful enough. In contrast, many trends were not ss when PW was applied, but were ss with the TFPW-Y procedure, leading to false positives and showing that the TFPW-Y rejection rate of the no-trend hypothesis is too high.

After having determined the statistical significance, a third pre-whitening procedure, the variance-corrected trend-free pre-whitening procedure (VCTFPW) allowing an increase in the slope accuracy (Wang et al., 2015), was applied prior to Sen's slope estimation. The confidence limits of Sen's slope were computed at the 90 % confidence level.

Since many of the time series exhibited clear seasonal cycles, the modified seasonal MK test (Hirsch et al., 1982) was always applied to the four meteorological seasons. The annual trends were considered only if the slopes of the four seasons were homogeneous at the 90 % confidence level (Gilbert, 1987; Sirois, 1998).

Figure 2 Seasonal MK results for the σ_sp trend for three stations with long time series: JFJ, MRN and MLO. The trends are plotted for the last 10-year period (2009–2018) as well as for all possible longer periods (15 years = 2004–2018 to 30 years = 1989–2018). The seasons correspond to meteorological seasons (MAM: March–April–May, JJA: June–July–August, SON: September–October–November and DJF: December–January–February). The dots correspond to the slope, large dots being ss at the 95 % confidence level, whereas small dots are not ss trends. The cyan triangles correspond to false positive trends (with type 1 error). Red squares correspond to annual trends where the seasonal results are homogeneous.

Figure 2 presents three examples of seasonal MK results and Sen's slopes of σ_sp. At JFJ, σ_sp has ss negative annual trends for all of the analyzed periods, with the most recent 10-year period with a larger negative slope than the longer periods. Spring and fall are the seasons at JFJ with the strongest ss trends; winter has tiny ss negative trends. MRN also exhibits σ_sp annual negative trends for all of the analyzed periods, but only the 15-, 20- and 25-year trends are ss, and their slopes are more negative for the longest periods. At MRN, summer and fall are the seasons with the largest trends, and that is true for all the trend periods (10–25 years), while spring and winter have more scattered and less significant slopes. Finally, MLO has annual trends that are ss negative for the last 10 years, not ss for the last 15 years, and ss positive for the longest periods (20, 25 and 30 years). The spring season at MLO exhibits a not ss negative trend for the last 10 years and positive trends for the longest periods, with only 25- and 30-year trends being ss.

2.5.2 Least mean square analysis (LMS)

Following the Weatherhead procedure (Weatherhead et al., 2000), the trend is estimated by fitting the following frequently used statistical model for monthly data with an LMS approximation:

$$\begin{array}{}\text{(1)}& Y_t=m+C_t+\mathit{\rho }\cdot (t/\mathrm{12})+M_t,t=\mathrm{1}\mathrm{\dots }n,\end{array}$$

where m is a constant term, C_t is a seasonal component, and ρ is the magnitude of the trend per year. The unexplained noise term M_t is modeled as an [AR(1)] process $M_t=\mathit{\varphi }\cdot M_{t-\mathrm{1}}+\mathit{ϵ}$, where ϕ is the autocorrelation coefficient of the data noise. For this study, either the logarithm of the monthly medians or the monthly medians were taken for all the parameters. Due to the non-normal distribution of the studied parameters, the LMS method applied on the logarithm is considered the standard method according to previous trend analyses (CC2013 and Asmi et al., 2013). A trend is considered to be ss at the 95 % confidence level if $|\mathit{\rho }/{\mathit{\sigma }}_p|>\mathrm{2}$, σ_p being the standard deviation of the slope. Figure 3a and c show the LMS trends and statistics for MLO σ_sp, respectively. The LMS results are similar to the MK analysis, the last 10-year trend is negative but ss at only the 90 % confidence level, the 15- and 20-year trends are not ss and the 25- and 30-year trends are ss positive. The normal probability plot of the residue (Fig. 3c) shows that the use of the logarithm of the data results in normally distributed residues as required by this statistical tool.

2.5.3 Generalized least squares associated with the autoregressive bootstrapping method (GLS/ARB)

A similar GLS method based on the minimization of the least square errors similar to ordinary least squares fitting (including similar sensitivity to outliers), but taking into account the autocorrelation in the covariance matrix, was also used in this study. The GLS uses an autoregressive bootstrapping algorithm (ARB) to evaluate the potential differences in the GLS trends arising from the noise terms (Asmi et al., 2013). The ARB methodology was used to produce 1000 realizations of the original time series, with randomized noise terms, and the resulting set of trends was used to determine the 5th to 95th percentile confidence intervals (ARB CLs) of the GLS trends. If the ARB CLs did not include a zero trend, we considered the GLS trend to be ss. The GLS and ARB methodologies were adapted from Mudelsee (2010) and applied to both daily and monthly medians. The previous trend analyses (CC2013 and Asmi et al., 2013) used daily medians.

Figure 3 LMS and GLS/ARB results of MLO σ_sp: (a) logarithm of the monthly medians (blue circles), LMS fit (red) and the 10- to 30-year slopes (ss slopes are plotted in black and not ss slopes in grey), (b) daily medians (light blue dots) and their GLS/ARB fit (orange line), monthly medians (blue circles) and their GLS/ARB fit (red) and the 10- to 30-year slopes, (c) normplot of LMS residues, (d) monthly medians of the GLS/ARB residues, (e) cumulative summation of monthly median GLS/ARB residues and (f) normplot of GLS/ARB residues for daily medians (light blue crosses and orange line) and monthly medians (blue crosses and red line).

Figure 3b and f show the GLS/ARB results for MLO σ_sp for daily and monthly medians. Here again the results are similar to the MK analysis, where the 10-year trend is positive but not ss, the 15-year trend is not ss, while the longer periods exhibit ss positive trends. As with many other stations included in this study, the use of daily or monthly medians did not result in normally distributed residues (Fig. 3f), and, in fact, the residues of the daily and monthly medians appeared to represent different types of distributions. It is also obvious that the seasonality fits (fits from monthly and daily medians in red and orange in Fig. 3b) are different for the two time granularities, with similar shape but higher absolute trend values if fitted from daily medians. The timing of the winter minima is also more precisely defined with the daily data.

3 Results

3.1 Long-term trends ending in the present day (2016–2018)

To assess the aerosol optical property long-term trends, the largest number of stations around the world was included in this study. This overview takes into account the 10-year (or longer) trends ending in 2016, 2017 or 2018. The results shown here comprise not only the 10-year trends, but also the longer periods for 15 to 40 years in 5-year increments also ending in 2016–2018. The results are presented for the MK analysis and a comparison between the trend analysis methods will follow in Sect. 3.3. Complete results for all the other methods can be found in the Supplement.

3.1.1 Total scattering and hemispheric backscattering coefficients

Long-term trend analysis of σ_sp has been performed on 37 datasets. Since some nephelometers only measure σ_sp (Optec and Radiance Research nephelometers) and σ_bsp was determined to be unusable for several other sites due to various discontinuities (see Sect. 2.4), the hemispheric backscattering coefficient trends were computed on only 28 datasets. The detailed results of MK trend analyses are given in Table 2, while the overall picture for σ_sp is presented in Fig. 4. The results for σ_bsp are very similar to those for σ_sp for sites where both measurements existed; corresponding figures for σ_bsp can be found in the Supplement (Figs. S1, S2 and S7).

Table 2 MK trends for all parameters in units per year for the last 10, 15 and 20 years of measurements ending in 2016–2018. The ss trends are given in bold. Results in % yr⁻¹ are given in Table S1.

Figure 4 MK trend results for the scattering coefficient. Black symbols correspond to stations with no significant trends. Green and orange symbols correspond to ss negative and positive trends, respectively. The magnitude of the trends (slope) is given by the colors as stipulated in the legend. The size of the circles is proportional to the length of the datasets, with the central dots representing the most recent 10-year trend ending in 2016, 2017 or 2018. If possible, trends for longer time periods were calculated and the larger circles denote the trends for 15 to 40 years in 5-year increments.

The σ_sp ss trends are predominantly negative: 20 stations have ss negative 10-year trends, 5 stations ss positive trends and 12 stations no ss trends dispersed across all continents. Eight (nine) stations with time series longer than 10 years have ss negative 15-year (20-year) trends and none (two) of the 15-year (20-year) trends are ss positive. The MK slopes range between −2.45 and +0.39 Mm⁻¹ yr⁻¹ with a mean of −2.19 Mm⁻¹ yr⁻¹. The main results are as follows.

• Over North America, all the σ_sp trends for periods longer than 10 years are ss negative, and the most recent 10-year trends are generally ss negative. Three stations do not have ss trends. (1) EGB's 9-year time series does not allow for a ss trend (too short), but was included as one of only two Canadian sites. (2) MRN is an IMPROVE station on the western coast of the USA with very high humidity, leading to condensation that can disturb the humidity measurement. This makes it difficult to know whether the ss positive RH 10-year trend (Table S4) is real or due to measurement artifacts and uncertainties. If the ss positive RH trend is real, it could mask a decreasing σ_sp trend, resulting in a not ss trend. The time coverage for the dry σ_sp (σ_sp restricted to RH < 50 %) for MRN is too low to be representative for trend analysis. It should be mentioned that the 10-year trends for MRN ending in 2014–2018 are all not ss (see Sect. 3.2.1), so that the absence of σ_sp trends seems to be a real phenomenon. (3) GLR is also an IMPROVE station with high humidity. The RH trends at GLR are also not ss, and the dry σ_sp has a ss negative trend, similar to other stations in its vicinity.

  In the previous decadal trend paper (CC2013), the trends in scattering for the arid state of Arizona were not consistent (ss positive: IBB, ss negative: SIA, PAZ, not ss: HGC, SCN). Four of the five Arizona sites (IBB, PAZ, SIA and SCB) were closed in 2010 and HGC now exhibits a ss decreasing scattering trend. MZW, the other IMPROVE station with ss positive scattering trends in 2010, also closed in 2010.

• Most (7 out of 11) of the European sites have present-day ss decreasing σ_sp trends. The other four stations do not have ss trends: (1) one urban station also influenced by Saharan dust (UGR), (2) two sites in eastern European countries (KPS and BEO) and (3) a high-altitude station in the Central Range in France (PUY). The ss negative scattering trends of the Scandinavian stations have lower absolute slopes than in central Europe. PAL, the northernmost station, has a ss positive trend. PAL is geographically situated in Europe, but it can be climatologically considered an Arctic station (Schmeisser et al., 2018). PAL (slope = 0.06 Mm⁻¹ yr⁻¹) has a similar trend to ZEP (slope = 0.05 Mm⁻¹ yr⁻¹), the nearest Arctic station, with the largest ss trend in summer (JJA) when PAL is largely influenced by Arctic air masses. The increasing trend at PAL may be due to increasing biogenic secondary organic aerosol formation related to emissions from the surrounding boreal forest (Lihavainen et al., 2015a), changes in circulation patterns or a larger influence of open water with increasing concentration of sea salt aerosol.

• Sites in the polar regions exhibit two ss positive σ_sp trends. In addition to ZEP and PAL, SPO also has a ss positive present-day 10-year trend but with lower slope, whereas no ss trend is found for the other Antarctic site (NMY). BRW and ALT both exhibit ss negative 10-year trends. The BRW 15-year σ_sp trend is not ss, whereas longer periods up to 40 years lead to ss negative trends. SPO also has very long time series but with alternating trend slopes, from ss positive for the shortest periods (10–25 years) to ss negative for the longest periods (35–40 years), with some not ss trends in between. The aerosol load is very low at BRW and SPO, leading to scattering coefficients near the instrumental detection limits, so that the measurement uncertainties are proportionally larger than for middle-latitude stations.

• CPR, a site on the Caribbean island of Puerto Rico, has a ss positive σ_sp trend. At CPR, the largest scattering trend is found in summer and the scattering trend of the PM₁₀ trend is 5 times larger than the PM₁ trend. The most probable explanation is increased Saharan dust transport over the Atlantic Ocean; more dust transport has been reported at an IMPROVE site in the Caribbean (Hand et al., 2017, 2019).

• The only two stations representing the Pacific region are MLO and CGO. The recent MLO 10-year σ_sp trend is ss decreasing, the σ_sp 15-year trend is not ss, whereas the trends for the longer time periods (20–30 years) are ss positive (see Fig. 2). In the previous decadal trend paper (CC2013), MLO exhibited a ss positive trend for the 10-year period ending in 2010. MLO σ_sp trends changed from previously ss positive to currently ss negative trends. The recent 10-year trend at CGO is found to be positive and quite homogeneous with the seasons, with fall being the only season without a ss trend.

• The σ_sp trends are mostly (70 %) not ss for stations at middle to high altitudes. From the 10 stations higher than 1100 m a.s.l., only SPO in Antarctica has a present-day ss positive 10-year trend and only JFJ in the European Alps, HGC in Arizona and GBN in Nevada exhibit ss negative 10-year trends. In contrast, only 26 % of the stations lower than 1100 m a.s.l. do not have ss trends. New particle formation (NPF) and growth are favored at high altitudes (> 1000 and up to 5000 m) due to low temperatures, high solar radiation and low pre-existing particle concentrations, leading to limited condensational sinks for nucleation precursor gases (Sellegri et al., 2019). This higher frequency of nucleation at high altitude leads to a high contribution of secondary particles to the total number concentration that largely contributes to the total scattering coefficient. The decreasing σ_sp trends from anthropogenic pollution in the planetary boundary layer can, consequently, be masked by the presence of NFP at high-altitude stations.

Figure 5 Seasonal results of the MK trend of the scattering coefficient. Other details same as Fig. 4.

Figure 6 MK trend results for the absorption coefficient. Other details same as Fig. 4.

The seasonal MK results for σ_sp are presented in Fig. 5. Spring is the season with the largest number of ss decreasing trends and winter the season with the lowest. ZEP and PAL exhibit ss positive trends only in summer and BRW has ss negative trends only between December and May. The SPO annual trend is ss positive, whereas it is not for NMY. Both Antarctic stations exhibit, however, a coherent seasonality with ss positive trends only in spring. While the 25- and 30-year trends at MLO are all ss positive, with the largest slope in spring when MLO is influenced by Asian long-range transport (CC2013), the most recent 10–20-year trends are not ss for the individual seasons.

3.1.2 Absorption coefficient

The analysis of σ_ap long-term trends has been performed on 33 datasets (see Fig. 6 and Tables 2 and 3). The long-term trends are ss decreasing (21 stations) or not ss (12 stations) for all stations around the world, leading to a mean decreasing trend of −3.05 Mm⁻¹ yr⁻¹. No ss σ_ap positive trends are measured for any of the stations. The other main results are the following.

• In North America the number of σ_ap datasets is much lower than the number of σ_sp datasets (IMPROVE sites do measure aerosol absorption, but with a different instrumental setup; White et al., 2016). From the five sites with long-term aerosol absorption, APP and BND, two continental rural sites, and the marine Caribbean island (CPR) station have ss negative trends. The other three stations representing the continental rural USA (SGP, EGB) and marine western coast of the USA (THD) exhibit not ss trends in σ_ap.

• In Europe, most (12 stations) of the 10-year σ_ap trends are ss negative. Only three stations, one Scandinavian (PAL), one eastern rural continental (KPS) and one coastal Mediterranean (FKL) station exhibit no ss trends. The 15-year σ_ap trends at JFJ and FKL are ss negative.

• In Asia, both the high-altitude stations of LLN in Taiwan and WLG in China exhibit annual ss decreasing σ_ap trends. The South Korean coastal station of AMY has no ss annual trend.

• For the polar regions, the Antarctica site of NMY, the American Arctic site of BRW and the Russian Arctic site of TIK have slight ss negative σ_ap trends, whereas SUM, ALT and ZEP have no ss trends. Thus, there is no common clear σ_ap trend in the polar regions.

• In the southwestern Pacific, the high-altitude station of MLO has a ss decreasing trend for the last 10 years but no ss trend for the last 15 years, whereas the coastal station of CGO in Australia exhibits not ss σ_ap trends.

• In contrast to the σ_sp trends, σ_ap trends at high-altitude stations (> 1100 m a.s.l) are mostly (6 out of 8) ss decreasing; the trends at the other two high-altitude stations are not ss.

The seasonal trends are more strongly negative and more ss in spring than in summer (see Fig. S3). Winter is the season with the smallest number of ss decreasing trends in Europe (only 2∕15) and with the only ss positive trend (ZSF), the others being not ss, whereas fall seems to be the season with the lowest ss trend in North America.

Figure 7 MK trend results for the single scattering albedo. Other details same as Fig. 4.

Figure 8 Seasonal results of the MK trend of the single scattering albedo. Other details same as Fig. 4.

3.1.3 Single scattering albedo

As described under Sect. 2.4, ω₀ trends have to be considered with greater caution since the σ_ap absolute values suffer from a certain uncertainty related to filter-based absorption photometer artifacts.

The ω₀ trends depend directly on both the magnitude and the sign of the σ_sp and σ_ap trends. If expressed in % yr⁻¹, a σ_ap trend larger (smaller) than the σ_sp trend will result in an increasing (decreasing) ω₀ trend, respectively (see Fig. S8 and the related estimation of ω₀ uncertainty due to measurement and C_ref errors). The ω₀ trends are consequently much more diverse than the σ_sp and σ_ap trends with 52 % of ss positive (relatively more scattering), 22 % of ss negative (relatively more absorption) and 26 % not ss trends (see Fig. 7 and Table 2). One peculiarity is that all ω₀ ss negative trends are found between latitudes 30 and 50, but this is perhaps due to the low spatial coverage outside of North America and Europe. The main results are the following.

Figure 9 MK trend results for the backscattering fraction. Other details same as Fig. 4.

• The ω₀ is decreasing at three stations in North America (BND, SGP and THD), whereas APP and CPR exhibit ss positive ω₀ trends. The CPR ω₀ increasing trend can perhaps be related to increased Saharan dust load. The seasonal ω₀ trends at CPR are, however, ss not only in summer when Saharan influence is greatest, but for every season except spring (Fig. 8). EGB has no ss trend.

• European stations exhibit ss increasing ω₀ trends at the urban station of UGR and at most eastern and Scandinavian stations (KPS, SMR, PAL) and at the mid-altitude station of HPB. These ss positive ω₀ trends in eastern and northern Europe are strongest in summer (Fig. 8), when MEL and BIR are also ss positive, and weakest in winter when only PAL is ss positive (possibly related to increased particle formation from biogenic emissions, as mentioned above). In central Europe, JFJ, IPR and MSY have ss negative ω₀ trends for the entire year as well as for all seasons. PUY, a station at 1465 m in France's Central Range, has a ss positive annual trend due to strong positive trends in fall and winter, even if a strong ss negative trend is found in summer. Because the site is located at a mid-range elevation (1465 m a.s.l.), PUY has a large probability of being influenced by different air masses as a function of the season, with a large impact of the planetary boundary layer in summer (Collaud Coen et al., 2018; Hervo, 2013).

• The high-altitude stations of LLN and WLG in Asia have strong and weak ss positive annual ω₀ trends, respectively. This pattern is also observed for all seasonal trends at LLN, but only in fall at WLG. The coastal station of AMY has a ss decreasing annual ω₀ trend that is due to decreasing trends in MAM, SON and DJF. AMY is located in an agricultural and touristic region that is influenced not only by these regional aerosol sources (e.g., traffic, field burning), but also by long-range transported plumes with high aerosol load.

• The Arctic stations of ALT and ZEP have ss positive ω₀ annual trends, which are due to ss positive trends from December to August for ALT and from December to May for ZEP. The two polar stations (BRW and NMY) exhibit no ss ω₀ annual trends, although there is a ss positive trend in summer at BRW for the most recent 10-year time series.

Figure 10 MK trend results for the scattering Ångström exponent. Other details same as Fig. 4.

Figure 11 MK trend results for the absorption Ångström exponent. Other details same as Fig. 4.

3.1.4 Backscattering fraction and scattering Ångström exponent

The present-day trends for the backscatter fraction b are mostly ss positive (65 %) across all regions (Fig. 9 and Table 2). This suggests a shift in the size distribution towards smaller accumulation-mode aerosol. The two stations with ss negative trends are CPR in Puerto Rico and BEO, located on a summit in the Balkan range. Not ss trends are mostly found in eastern and northern Europe (KPS, MEL, BIR and PAL), in Antarctica (NMY), as well as at BND and AMY for the last 10 years. The Arctic sites (ALT, BRW and ZEP) all exhibit ss positive b trends. CPR's seasonal trend is ss negative only in fall; trends in b for the other seasons at CPR are not ss (see Fig. S4). Similarly, the BEO b seasonal trend is ss negative only in summer, and not ss otherwise. PAL has ss positive b trends in spring and summer and ss negative b trends for fall, leading to an annual not ss trend.

The scattering Ångström exponent (å_sp) trends exhibit a higher variability than the trends in other parameters, with 33 % of ss positive trends, 37 % of ss negative trends and 30 % of not ss trends (Fig. 10 and Table 2). There are ss positive and negative trends in North America, Europe and the polar regions, and the various trends cannot be attributed to specific regions or environments. It should be recalled, however, that å_sp is affected by higher uncertainties (see Sect. 2.3) that may contribute to the larger observed variability. The seasonal results also exhibit high variability, with summer being the season with the least number of ss å_sp trends (10 out of 26 sites), while spring and fall are the seasons with the largest number of ss positive and negative trends in å_sp (8 out of 26 sites), respectively (see Fig. S5).

3.1.5 Absorption Ångström exponent

The number of stations with long-term å_ap measurement is low, with only 14 time series available. Seven stations situated in various geographical regimes exhibit ss positive trends (Fig. 11 and Table 2): polar regions (ALT and ZEP), a Caribbean coastal station (CPR), a rural continental North America station (SGP), and high-altitude stations in the remote Pacific (MLO), in continental (WLG) and in coastal (LLN) Asia. SMR and JFJ, two stations in Europe but with very different environmental footprints and altitudes, exhibit ss decreasing å_ap trends; the six other stations, consisting of three coastal and two continental sites, have no ss trends.

While CPR and SGP å_ap trends are ss positive and JFJ ss negative for all four seasons, the other stations exhibit higher variability as a function of the meteorological seasons (see Fig. S6). The absorption Ångström exponent is principally a function of the particle chemical composition and material properties, but its assignment to an aerosol type is not uniquely defined and also depends on the particle size, with larger particles corresponding to lower å_ap values (Liu et al., 2016; Schmeisser et al., 2017). For example, å_ap>2 corresponds to mineral dust in the case of big particles and to brown carbon in the case of small particles. In contrast, å_ap<1 corresponds to large particles with small absorption like sea-salt-dominated aerosol in the case of big particles and to black carbon (BC)-dominated aerosol in the case of small particles. Following these observational constraints, the JFJ and SMR aerosol tends to represent the category “mixed BC/BrC” according to Schmeisser et al. (2018). CPR absorption has a strong contribution from mineral dust and sea salt, whereas at MLO, SGP, ALT and ZEP contributions to absorption are from mixed sources, including various light-absorbing carbon species and dust. Ideally, direct chemical composition measurements would provide more precise information on the aerosol type, but the necessary chemical composition measurements are not yet readily available at many sites.

Figure 12 Time series of sequential 10-year σ_sp trends as a function of station longitude. Stations in the South Pacific and in the polar regions were grouped for clarity. The red arrow indicates the end of the time periods covered in CC2013.

3.2 Time evolution of 10-year trends

The previous section describes the present-day trends for different periods extending from 10 to 40 years. Another interesting analysis is to follow the evolution of the trends in time and space. For this purpose, all the possible 10-year trends were computed and plotted as a timeline for each station. In what follows each point on the timeline represents a 10-year trend ending in the year it is located on the graph. For example, in Fig. 12 the two black points for AMY represent the 10-year trends covering the periods 2008–2017 and 2009–2018, respectively. These timelines can be presented as a function of the latitude, longitude, altitude or environment. Depending on the results, the most interesting representation has been chosen for each parameter.

3.2.1 Scattering and backscattering coefficients

Figure 12 presents the σ_sp 10-year trend timelines as a function of the longitude of the stations. The σ_bsp 10-year trend timelines are similar to the σ_sp trend timelines (see Fig. S7). The polar stations (in both the Arctic and Antarctic) have been gathered to the bottom of the figure just after the two Pacific stations of MLO and CGO. The main result is that the sites in eastern and central North America (longitude between −68 and −112^∘) have ss negative σ_sp 10-year trends ending after 2009–2012 regardless of their altitude (200–2200 m a.s.l.) and their environments. This is a clear signature of continental-scale modification due to air quality regulations, and this very clear feature relates to the sulfate-dominated aerosol in the eastern USA and to large SO₄ reductions in power plant emissions (Hand et al., 2014; McClure and Jaffe, 2018). Almost all the σ_sp 10-year trends in the southwestern USA (MZW, SCN and HGC) ending before 2011 are ss positive, as published in the previous trend analysis (CC2013). MLO also exhibits ss positive trends for the same period. These four stations are also high-altitude sites (2000–3400 m), so that it is possible that all of them were influenced by long-range transport of highly polluted air masses from Asia (CC2013). It is further interesting to note that the high-altitude site JFJ (3580 m) in Europe also exhibited a ss positive 10-year trend ending in 2005–2008.

The evolution of the European σ_sp 10-year trends does not show a clear time for trend modification as is seen in North America, probably due to variable timing in implementation of abatement policies in each individual country. Apart from PAL (which can be considered, to some extent, to be a polar station), the σ_sp 10-year trends in Europe ending after 2008 are all ss negative or not ss. The four stations in Asia do not have ss trends ending in the last 5 years. The two African stations exhibit no ss trend. For polar sites, BRW, ALT and NMY have mostly not ss trends, whereas SPO exhibits alternating ss positive and negative trends, with the oldest 10-year trends being not ss. In contrast, ZEP exhibits positive trends for all three 10-year periods, which is similar to the 10-year trends at PAL for the same time periods. Due to the very low aerosol concentrations at these sites and, thus, larger measurement uncertainty, it is difficult to interpret the evolving σ_sp polar trends. They could be related to increased influence from the boreal forest and/or changed circulation patterns modifying the sea/ice influence.

The GSN dataset only covers 8 years, with some missing periods due to the destruction of the station by a typhoon. Due to the very low number of long-term measurements in Asia, GSN was included in this study. While GSN σ_sp summer trends are not reliable (low data coverage and issues in humidity control), the ss negative winter–spring trends corresponding to the dry season are valid and in line with the PM₁₀ decreasing trends in Korea (Kim and Lee, 2018; Nam et al., 2018).

Figure 13 Time series of sequential 10-year σ_ap trends as a function of station longitude. The red arrow indicates the end of the time period covered in CC2013.

3.2.2 Absorption coefficient

The lengths of the σ_ap time series are much shorter than for σ_sp (Fig. 13). This means that the oldest 10-year trends cover the period 1998–2007 (BRW and BND), followed by MLO (2001–2010) and JFJ (2002–2011). For these four stations, the most recent 10-year trends are either not ss or ss negative. The present-day (i.e., trends covering 2009–2018) ss negative 10-year trends (JFJ, BND, MLO and BRW) are preceded by not ss trends. The σ_ap 10-year trend evolution of each station is usually homogeneous with either ss negative or not ss 10-year trends in Asia, Europe, Africa and North America. ALT σ_ap 10 years ending in 2017 is ss positive. The BRW polar station and MLO high-altitude station exhibit also some ss positive 10-year trends ending between 2010 and 2014. Unfortunately, only MLO has a long enough σ_ap time series to compare with the ss positive σ_sp 10-year trends at high-altitude sites (Fig. 12). At MLO, the series of ss positive σ_sp 10-year trends ended in 2008, while the series of ss positive σ_ap 10-year trends occurred for the period ending 2009–2013.

Figures 12 and 13 suggest that mid-latitude σ_sp and σ_ap sequential 10-year trends were ss positive for some periods between 2000 and 2013, followed by not ss trends and ending in the present day with ss negative trends. The evolution from increasing to decreasing σ_sp and σ_ap trends appears to be not simultaneous, with the σ_sp inflection points occurring some years before those for the σ_ap trends. The sparse number of stations with long enough time series does not allow generalization of this result.

Figure 14 Time series of sequential 10-year ω₀ trends as a function of station longitude.

Figure 15 Time series of sequential 10-year b and å_sp trends as a function of station longitude. The red arrows indicate the end of the time period covered in CC2013.

3.2.3 Single scattering albedo

Because it is limited by the length of σ_ap time series, the ω₀ 10-year trend evolution also only covers the last decade. The following results can be seen in Fig. 14.

• All stations at longitude > 10^∘ have ss positive ω₀ 10-year trends except for AMY, which exhibits a ss negative 10-year trend ending in 2018, and MUK with a not ss 10-year trend ending in 2013. For European sites, ss positive ω₀ 10-year trends exist for all stations at latitude > 46.8^∘, apart from BIR, which has a not ss trend ending in 2018. This suggests that the decreasing σ_ap trends in Asia and in eastern and northern Europe are proportionally larger than the decreasing σ_sp trends.

• The central and western European sites exhibit mostly ss negative or not ss ω₀ 10-year trends. At JFJ and IPR, a shift from not ss to ss negative 10-year trends occurred in 2013–2014. The JFJ time series is moreover long enough to monitor a ss positive ω₀ 10-year trend ending in previous years (2010). The urban station of UGR in Spain exhibits an increasing trend in ω₀ (decrease in contribution of absorbing aerosol) for the most recent 10-year period (2009–2018), possibly related to long-term effects of the 2008 financial crisis (e.g., Lyamani et al., 2011).

• In North America, the ω₀ 10-year trends ending after 2013 are ss negative or not ss, apart from CPR in Puerto Rico and APP. The sites with the longest series of ω₀ 10-year trends (BND and THD) exhibit ss positive trends followed by not ss and ss negative trends. In contrast, MLO ω₀ 10-year trends shifted from ss negative trends (10-year trends ending in 2010–2015) to not ss trends (10-year trends ending 2017) to a ss positive trend in 2018. This is consistent with the observed increase in σ_ap 10-year trends ending in 2010–2012.

• The two polar sites of BRW and NMY exhibit mostly not ss ω₀ 10-year trends, whereas ALT and ZEP, similar to northern European sites, exhibit ss positive ω₀ for all the 10-year trends.

3.2.4 Backscattering fraction and scattering Ångström exponent

Both the b and å_sp 10-year trends in Asia and Africa exhibit similar 10-year trend patterns that are either ss positive or not ss (Fig. 15). In this context, similar means that the b and å_sp trends are never ss when opposite signs of the slope are observed. These results suggest that particle average size tends to decrease at the Asian and African sites. In Europe, b and å_sp 10-year trends have a majority of ss negative or not ss trends in the northeast (longitude > 10^∘). At lower European longitudes, there is a discrepancy between b and å_sp 10-year trends, with IPR, JFJ and PUY having opposite ss trends for b and å_sp, b trends being often ss positive and å_sp trends often ss negative. The discrepancy in the signs of the trends for b and å_sp may be related to shifts in both the fine and coarse modes of the aerosol size distribution – this is discussed more below (see Sect. 4.2).

In North America, the b 10-year trends ending after 2012 are almost all ss positive, whereas the previous 10-year trends are ss negative at BND and MLO. As in western America, one can see discrepancies in the sign of the slope for b and å_sp 10-year trends, with APP, CPR, MLO, SGP and THD having at least one 10-year period with opposite signed ss trends. Also, as in western Europe, the b trends are usually ss positive and the å_sp trends ss negative. BND records four 10-year periods with opposite ss trends. In contrast to the other stations, BND å_sp trends are ss positive, while BND b exhibits ss negative trends.

In the polar regions, the b 10-year ss trends ending after 2014 are all ss positive, whereas the older trends are primarily ss negative. Here again, the discrepancy between b and å_sp 10-year trends is large, with all the 10-year b and å_sp trends for ZEP and ALT being ss with opposite signs. In contrast, BRW exhibits trends with the same sign for both parameters that can be interpreted as an increase in average particle size for early years followed by a decrease after 2014.

Figure 16 Time series of sequential 10-year å_ap trends as a function of station latitude.

Table 3 Number of trends analyzed for each parameter, of ss cases for each trend analysis method, of trends with similar statistical significance in MK, GLS/day and LMS/log, of trends with similar statistical significance for all five methods, of trends with at least GLS/day or LMS/log ss similar to MK, and of trends with no agreement between these two methods and MK ss.

3.2.5 Absorption Ångström exponent

The å_ap time series are not very long, because the first generation of absorption photometers used either white light or only one wavelength (Fig. 16). The longest time series of å_ap begins in 2002 at JFJ and exhibits a continuous ss å_ap decrease. Similar to the results for JFJ, most of the stations have consistently ss negative (JFJ, SMR), ss positive (ALT, ZEP, IPR, SGP, WLG, LLN, MLO and CPR) or not ss 10-year trends (TIK, BRW, THD, APP, GSN, MUK and CPT). The not ss trends of TIK and GSN may be due to datasets shorter than 10 years.

3.3 Comparison of the trends among methods

As described under Sect. 2.5, the long-term trends were computed with three methods (MK, GLS and LMS), where GLS was used on both daily and monthly medians and LMS with and without taking the logarithm of the monthly medians. These methods are thereafter called GLS/day, GLS/month, LMS/log and LMS/lin, respectively. Tables S2 and S3 give the GLS/day and LMS/log results for all parameters and stations. Table 3 presents an overview of the number of present-day 10-year trends that are ss with each method. For the reasons described in Sect. 2.5.1, MK is considered to be the most appropriate method for aerosol optical parameters. The agreement between the three methods used in CC2013 (MK, GLS/day and LMS/log) and between all five methods are then also reported in Table 3, as well as the number of cases with either GLS/day or LMS/log agreement with MK and with both GLS/day and LMS/log disagreement with MK. The following conclusions can be derived.

• Generally, the trends computed by the various methods agree very well with one another. Among all parameters, all stations and all periods, none of the present-day trends presents ss results with opposite slope for different methods. In all cases, the differences among methods relate either to the degree of the ss or to the sign of the slopes for not ss trends. This implies that the main conclusions of this study would not have been fundamentally different if the other methods were used.

• GLS applied on daily medians is the method that has the largest number of ss trends for all parameters.

• The three methods applied on monthly data have a lower number of ss trends for all the computed parameters (ω₀, b, å_sp and å_ap).

• The three methods used in 2013 have similar statistical significance (comprising cases with no ss trend) in 44 % to 86 % of the cases, whereas the five methods used here exhibit consistency in 37 % to 82 % of the cases. The measured parameters, which are less uncertain than the calculated parameters, always exhibit the largest agreements amongst the methods (> 69 % for the three methods used in 2013 and > 63 % for the five methods utilized here). ω₀ is always the parameter with the largest dissimilarity among the methods and σ_bsp the parameter with the largest similarity among methods.

• The MK statistical significance is similar to at least one of the methods applied in 2013 in more than 90 % of the cases for all of the parameters apart from ω₀ (88 %) and å_ap (78 %). This lower level of agreement can be explained by the fact that ω₀ and å_ap are almost normally distributed, so that the use of the LMS/log is not appropriate.

Figure 17 Median, interquartile ranges and whiskers of the slopes in % yr⁻¹ computed by the five methods for the scattering coefficient, the absorption coefficient, the single scattering albedo, the backscattering fraction and the scattering Angström exponent. The outliers are not always visible in the figure for the purpose of clarity.

The boxplots of the slopes computed by the various methods (Fig. 17) show first that the application of the logarithm to transform to a normal distribution for ω₀ and å_ap (not shown) is not suitable and leads to very large interquartile ranges. While the measured parameters are clearly not normally distributed, the derived parameters usually have distributions that more closely approximate normal distributions. No systematic rule could be deduced, since the distributions of each computed parameter largely depend on the individual stations. It seems however that ω₀ and the Ångström exponents are closer to the normal distribution than to the lognormal distribution.

The Sen slope estimator applied to the variance-corrected pre-whitening (Wang et al., 2015) leads, in almost all cases, to a median of the slope nearer to zero than the other methods. The VCTFPW method was developed specifically to get rid of the falsely increased slope by the trend-free pre-whitening process (Collaud Coen et al., 2020a). The LMS/log method sometimes results in lower absolute slope medians, and this effect is probably due to the almost normal distribution of the data (= log of the monthly median). Both the GLS (GLS/day and GLS/month) and LMS/lin methods lead to higher absolute slopes, probably due to misuse of statistical methods developed for normally distributed data.

The GLS/day method leads to a broader range of slopes than the GLS/month method. This larger variance may be due to (1) the larger variability of daily data, leading to a less distinct seasonal cycle and, consequently, to a worse fit of the seasonal variation and (2) a higher autocorrelation in the daily time series with, possibly, an autocorrelation order larger than one.

4 Discussion

4.1 Considerations related to measurement humidity

As explained in the instrumental section, GAW protocol suggests that the σ_sp and σ_bsp be measured at low and controlled humidity, and that is the case for almost all stations considered here, except for those in the IMPROVE network which measure at ambient conditions due to their different monitoring goals. Temporal cycles and variations of RH with time are observed in a number of datasets. There are also some clear breakpoints in measurement RH that have been identified at several stations (e.g., an insulating jacket was installed on the nephelometer at THD in late 2012, resulting in a clear decrease in sample RH due to warmer nephelometer temperatures). It is evident that high RH will enhance particle diameters and, consequently, increase σ_sp, σ_bsp and ω₀ while resulting in decreased b and å_sp. This particle diameter enhancement depends not only on the RH values, but also on the particle hygroscopicity, which is a function of the aerosol size distribution and chemical composition.

Similarly to the previous aerosol optical property trend study (CC2013), dry σ_sp was calculated by removing data when measurement RH was higher than 50 % in order to minimize the impact of aerosol hygroscopicity on the scattering trends. However, hygroscopic growth can occur for RH < 50 %; for example, for sea salt aerosol, up to 25 % of the scattering could be due to water at RH = 40 % (e.g., Fig. 5 in Zieger et al., 2013). The confounding effects of aerosol water impact the reported scattering values and, hence, the trends presented here to a greater or lesser extent. The effect of hygroscopic growth at RH < 50 % on the reported trends would depend on the temporal variability in sample RH, composition and size; investigating the interactions amongst those parameters is beyond the scope of this study.

For this study, if RH was frequently larger than 50 % at a station, relationships between RH and aerosol parameter trends were analyzed as follows. In the case of the RH trend being not ss, the aerosol parameter trends were considered to be independent of the RH variation. In the case where a ss RH trend was detected (see Table S4 in the Supplement), an attempt was made to try to determine the influence of RH trend on each aerosol parameter by considering the following situations: (1) if all aerosol trends follow the RH trends, (2) if σ_sp at all measurement RH and dry σ_sp trends are similar, and, finally, (3) the features of σ_ap trends, which are less likely to be influenced by long-term RH variation. The distinct patterns exhibited by the evolution of the 10-year trends were very helpful in this analysis. Below we describe the assumed implications for scattering trends at sites where trends in RH were observed for several cases.

• Trends in RH are the opposite of both σ_sp and σ_bsp trends: this implies that the aerosol optical property trends are real and not influenced by humidity (SMR, SHN, MRN).

• Trends in RH are ss, but trends in σ_sp are not ss: this implies that the absence of statistical significance for the σ_sp trends is real if the slopes of the RH and σ_sp trends have the same sign (IZO, LLN) or can be partially induced by the RH trend if the slopes have opposite signs (EGB, PUY, UGR).

• Trends in RH and σ_sp are similar, the overall and dry (RH < 50 %) σ_sp trends are similar, and σ_sp and σ_ap exhibit similar trends: this implies that the σ_sp trends are probably influenced by RH but also have an intrinsic aerosol trend (APP, BIR, MZW, SGP).

• Trends in RH and σ_sp are similar, but the dry (RH < 50 %) and overall σ_sp trends are dissimilar and the trends in σ_sp and σ_ap are also dissimilar: this implies that the RH influence is major (THD, CPR). THD and CPR are coastal stations with a dominant influence of sea salt. At THD on the North America western coast, RH, σ_sp and σ_bsp trends are ss decreasing, whereas b and å_sp trends are ss increasing and σ_ap trends are also ss decreasing but with lower slope and ss than σ_sp. Further, the PM₁ trends were less ss and exhibited much lower slopes, suggesting that the large sea salt particles are probably sensitive to the RH decrease, leading to the decreasing σ_sp trend. The 10-year trends show that RH decreasing trends are particularly important until 2015 and likely explain the ω₀ ss positive 10-year trends ending in 2012 and 2013. At CPR, a coastal site in the Caribbean, RH, σ_sp and σ_bsp trends are ss increasing, but the σ_ap 10-year trends do not have the same shape or statistical significance as the σ_sp trends. As observed at THD, the PM₁ trends at CPR are less ss and have much lower slopes than PM₁₀ trends.

As mentioned in the instrumental section, RH trends measured by the nephelometers have to be considered with caution. Because the measurement RH is only a secondary parameter, the instrument humidity sensors are typically not maintained or calibrated with the same care as the scattered light detectors. The influence of humidity variations on the optical property trends presented here can generally be considered to be low, apart from the cases of very hygroscopic particles like sea salt (e.g., at THD and CPR). A better knowledge of the particle hygroscopic growth at low RH (< 40 %) would be valuable in order to interpret σ_sp and σ_bsp trends as well as trends in ω₀, b and å_sp.

4.2 Particle size trends

Both the scattering Ångström exponent and the backscattering fraction are indicators of the particle's average size, with the general interpretation that lower values of b and å_sp correspond to the presence of larger particles albeit at different parts of the aerosol size distribution (Collaud Coen et al., 2007). However, the relation between b and å_sp is not uniquely defined for several reasons. First, the scattering efficiency has an oscillating response to particle size rather than a constant increase. Second, the measured particle size distribution is usually composed of several modes. Since the sensitivity of scattering to the mode depends on the size parameter (proportional to the ratio diameter–wavelength), b (here usually taken at 550 nm) and å_sp (here usually computed with the 450–550 nm pair) do not always exhibit similar sensitivity to the various size modes. Further, the extinction Ångström exponent (analogous to å_sp) was found to be more sensitive to fine-mode volume fraction if computed from long wavelengths and to fine-mode effective radius if computed from short wavelengths (Schuster et al., 2006). Lastly, the relation between b and å_sp also depends on the refractive index and consequently on the absorption coefficient (Hervo, 2013): for a constant particle diameter, an increase in the refractive index real part will decrease å_sp but increase b.

In this analysis, some stations exhibit b and å_sp trends with the same sign (BRW, CPT, SMR, LLN, SPO, UGR), while for other stations b and å_sp trends are in opposite directions (ALT, APP, BEO, BND, EGB, JFJ, MLO, MSY, PAL, PUY, SGP, SPO, ZEP). The plots showing the evolution of the 10-year trends (Fig. 15) demonstrate that b and å_sp can exhibit either similar or opposite trends depending on the considered periods (CPR, IPR, MLO, THD). The plots showing the evolution of the 10-year trends also suggest that the variations of the 10-year slopes are often identical in sign but with different magnitude (e.g., shifted towards larger trend values for b; see for example MLO, Fig. 18).

Figure 18 10-year slopes of b and å_sp at MLO. Ss negative and positive trends are plotted in blue and red, respectively. Not ss trends are plotted in black. The dots correspond to PM₁₀ and the upwards triangles to PM₁.

We can attribute both b and å_sp ss positive trends (ALT, BRW, SMR, LLN, THD, UGR) to a shift in the accumulation mode towards smaller sizes and a decrease in the coarse-mode particle concentration. In contrast, ss negative trends (BEO, CPR) for both b and å_sp suggest a shift to bigger sizes, specifically an increase in the coarse-mode particle concentration and perhaps also a shift towards larger diameters of the accumulation mode. At a boreal forest site in northern Europe (SMR), size distribution data suggest that seasonal variation of b and å_sp was caused by a shift in the accumulation mode and not by changes in the coarse-mode fraction (Luoma et al., 2019). Trends towards smaller particle sizes might be due to an increase in near-anthropogenic sources of pollution, to an increase in new particle formation, to a decrease in long-range transport of anthropogenic pollution, to increased scavenging of larger particles due to changes in atmospheric conditions, to a modification of atmospheric chemistry (Banzhaf et al., 2015) or to a change in both primary and secondary natural aerosol (e.g., an increase in biogenic secondary aerosols and their precursors as demonstrated by Ciarelli et al., 2019). Trends towards bigger particles can relate to a decrease in near-anthropogenic emissions, to larger influence of mineral dust caused by variation in desert emissions or dust transport, to changes in agricultural activities or to an increase in humidity.

For stations with opposite b and å_sp trends, the chemical composition may play an important role in identifying reasons for the changing trends. It is however out of the scope of this paper to study these kinds of dependencies.

Table 4 Overview of the aerosol optical property decadal MK trends ending between 2016 and 2018 for all the stations and per continent/region of the world.

4.3 Single scattering albedo trends

The single scattering albedo is the most important variable determining the direct radiative impact of aerosol, so that its trend analysis – derived for the first time for a large number of stations – has a high relevance. The filter-based absorption photometer artifacts lead to uncertain absorption absolute values that have no effect on σ_ap trends but impart higher uncertainties to ω₀ trends. The results of ω₀ trends depend directly on the relative values of σ_sp and σ_ap trends. The global picture is nuanced, with about half of ss positive, one-fifth of ss negative and one-fourth of not ss trends leading to an annual positive median trend of 0.02 % yr⁻¹ (Table 4). The medians of ss trends are increasing in Asia, in the Arctic and in the Pacific but decreasing in Europe and North America. The largest median slopes are found in Asia and in the Pacific (+0.13 and 0.14 % yr⁻¹, respectively), whereas the decreasing median slopes in other regions are relatively small (< 0.01 % yr⁻¹). The beginning of the decrease in the aerosol burden varies with region; the earliest decrease is found in Europe in the 1980s (Tørseth et al., 2012), followed by North America in the 1990s (Bodhaine and Dutton, 1993; Hand et al., 2012) and by Asia some 10–15 years ago (Sogacheva et al., 2020; Zhao et al., 2019; Paulot et al., 2018). The median slope of the ω₀ trends seems to be proportional to the length of the mitigation efforts, which for some relevant pollutants (e.g., black carbon, SO₄ and NO_x) are still ongoing. In Europe, the diversity of the timing of abatement policies with earlier impact in western Europe than in eastern Europe (Vestreng et al., 2007; Crippa et al., 2016; Huang et al., 2017) is also directly visible in the decreasing and increasing ω₀ trends (Figs. 7 and 14), respectively.

These results suggest that policy regulations induced first a ω₀ increase (cooling effect) and, in a second phase, a ω₀ decrease (warming effect). The Emission Database for Global Atmospheric Research (EDGAR V4.3.2) (easy accessible via https://eccad.aeris-data.fr/, last access: 20 July 2020) shows that both the BC and SO₄ emissions decreased rapidly during the 1990s and that currently emission reductions of SO₄ are larger than the reductions for BC. From this we conclude that the reduction of primary particles, such as BC, leads first to the ω₀ increase, whereas the reduction of SO₄, a precursor of secondary particle formation, tends to result in a ω₀ decrease. Moreover, emission changes can lead to modification of the atmosphere chemistry. Banzhaf et al. (2015) show, for example, that sulfate and nitrate formation increased in efficiency by factors between 20 % and 25 % between 1990 and 2009. The decrease in sulfate and total nitrate concentrations is consequently smaller than expected (non-linear response), leading to lower trends than the trends in precursor emissions and concentrations. This different timing and evolution in primary and secondary aerosol concentrations could explain the evolution of the 10-year ω₀ trend at IPR, JFJ, BND and THD (Fig. 14), but the time series are not long enough to properly assess this change.

These observed ω₀ trends are in line with the modeled impact of aerosol on climate (Zhao et al., 2019). They found that a global cooling effect of −0.41 K due to growth of aerosol burden caused by an increase in energy use in the Northern Hemisphere (particularly in Asia) is counterbalanced by a global warming of +0.10 K caused by the decreased aerosol emissions due to technology advances particularly in North America and Europe. This illustrates the complex nexus of environmental pollution regulations which have positive effects for health and the environment (air pollution is a primary cause of premature deaths in much of the world, Landrigan et al., 2018), but may have an adverse effect on efforts to reduce climate change. Ideally, abatement policy aimed at decreasing atmospheric pollutant levels would take into account both climate and health impacts.

4.4 Comparison with other trends and causality

The current study has focused on surface in situ aerosol optical properties at point locations, primarily in North America and Europe, but also in Asia and the polar regions. Comparison with reported trends from other long-term measurements of aerosol properties (e.g., surface aerosol mass concentrations, surface chemical mass concentrations, ground-based and satellite column optical properties) can provide a more holistic and global view of changes in the atmospheric aerosol. Model simulations of aerosol trends can also supply insight into global impacts of emission changes. We, thus, present a (non-exhaustive) comparison of the trend results from this study with some other relevant aerosol trend studies in the literature. The Supplement of Li et al. (2017) includes a summary of trends reported in the literature for aerosol optical depth (AOD), PM_2.5 and several aerosol constituents (e.g., sulfate, BC).

There are some important caveats to keep in mind when comparing aerosol trends across platforms and instruments. First, they represent different aspects of the aerosol (chemical, physical, or optical), at different conditions (dry or ambient), different wavelengths (300–1100 nm), different techniques (in situ, REM) and different locations (ground-based, airborne or satellite). Second, there are differences in the statistical methodologies, both in terms of methods used and data treatment. Third, the periods covered often overlap, but are not the same. Further, some REM measurements can only be made under certain conditions (e.g., daylight and cloud-free conditions versus continuous sampling, over land versus over ocean), meaning temporal coverage may be quite different. Because of all these differences, we only discuss general tendencies rather than absolute values when comparing trends from different studies. Below we first compare our results with trends from other surface in situ measurements and REM observations. Finally, we discuss causes of the observed trends and speculate specifically on some of the trends in intensive aerosol properties, which have received less attention in the literature than properties related to aerosol loading.

4.4.1 Comparison with other surface, in situ aerosol trends

A comparison of the present-day trends derived here to our previous trend ending in 2010 (CC2013) demonstrates that the larger number of stations, particularly in Europe, permits a more detailed view of regional trends. The current wide coverage across continental Europe shows decreasing present-day trends. Decreasing σ_sp, σ_bsp and σ_ap trends were confirmed for individual stations (e.g., SMR, Luoma et al., 2019; PAL, Lihavainen et al., 2015b; ARN, Sorribas et al., 2019), as well as at ACTRIS sites including JFJ, HPB, IPR, IZO, PAL, PUY, SMR and UGR (Pandolfi et al., 2018). There are some discrepancies in the trends between our current study and Pandolfi et al. (2018) that seem to be principally due to differences in the analyzed periods. Three additional years of data were included in this study and some older periods included in Pandolfi et al. (2018) were invalidated following the evaluations described in Sect. 2.4. The European b and å_sp trends computed by Pandolfi et al. (2018) are similar to the results of this study for most of the stations, in that they also found a general ss increase in b and variable å_sp trends. In North America the ss decreasing trends in aerosol extensive properties observed in CC2013 are found to have continued in this work with the extended datasets. These results are confirmed by the two other trend studies for in situ aerosol optical properties in North America. While the methodology and time period of Sherman et al. (2015) were different, the sign and ss of their σ_sp, b, and å_sp trends for BND and SGP were the same as reported here. White et al. (2016) found a decreasing trend in absorption coefficient (estimated from light transmittance measurements on 24 h filter samples) at 110 IMPROVE stations for the 2003–2014 period. SPO σ_sp, b and å_sp trends for the 1979–2014 period (Sheridan et al., 2016) do agree with CC2013 results, whereas the 1979–2018 trends reported in this study suggest an evolution towards more ss positive trends. The very low aerosol concentrations in Antarctica and the difference in the MK algorithm could however also explain the differences amongst these three analyses.

There have been multiple trend studies on carbon species (also referred to as BC, elemental carbon, equivalent black carbon, brown carbon or other terms) which are closely related to aerosol absorption. A decreasing trend in BC concentration is found in Europe (Singh et al., 2018; Kutzner et al., 2018; Grange et al., 2020), related primarily to traffic emission decreases rather than changes in wood burning and/or industrial emissions. Similarly, Lyamani et al. (2011) noted a decrease in BC in southern Spain due to the 2008 economic crisis. In contrast, Davuliene et al. (2019) reported an increasing trend in equivalent black carbon (eBC) for the Arctic site of TIK. In North America, White et al. (2016) found that the decreasing elemental carbon trend at IMPROVE sites was larger than the aerosol absorption trend at the same sites due to the impact of Fe content in mineral dust. BC trends in the Arctic have been extensively studied (e.g., AMAP, 2015; Sharma et al., 2019, and references therein) and suggest a decreasing trend. This is consistent with our general trend in absorption for the polar regions (Table 4), although for individual stations most trends were statistically insignificant.

Particulate mass (PM) and visibility are other metrics for atmospheric aerosol loading that can be most readily compared with our trends in aerosol scattering. Tørseth et al. (2012) detailed decreases in PM across Europe, while Hand et al. (2014, 2019) report significant decreases in PM_2.5 mass across the USA, with larger trends in the eastern than in the western USA. Both these trends were also confirmed by the PM trend analysis in Mortier et al. (2020) and are consistent with our reported scattering trends. Li et al. (2016) used visibility to assess trends in atmospheric haze and aerosol extinction coefficient around the world. The time delay in when the trends switch sign between North America (late 1970s), Europe (early 1980s) and China (mid 2000s) correlates with SO₄ trends, and the trend differences between the eastern and western parts of the USA and Europe are consistent with what is presented in our study.

Many atmospheric aerosols are formed in the atmosphere rather than being directly emitted, so understanding trends in aerosol precursors is also relevant for understanding changes in the atmospheric aerosol. Our study found similar results for scattering to those found for sulfate trends (Aas et al., 2019), i.e., decreasing sulfate trends across Europe and the USA, albeit with the sulfate decrease in Europe beginning before the decrease was observed in the USA. Aas et al. (2019) also describe potential increases in sulfate in India and increases followed by decreases in SE Asia. Vestreng et al. (2007) monitored the sulfur dioxide emission reduction in Europe and concluded that SO₄ emission reductions were largest in the 1990s, with a first decrease in western Europe in the 1980s followed by a large decrease in eastern Europe in the 1990s. Similarly, Crippa et al. (2016) simulated a larger impact of policy reduction in western Europe than in eastern Europe for NO_x, CO, PM₁₀ and BC between 1970 and 2010. Likewise, Huang et al. (2017) simulated the non-methane volatile organic compound emissions and found a rapid decrease in Europe and in North America since the 1990s, whereas the emissions of Africa and Asia clearly increased between 1970 and 2012.

4.4.2 Comparison with remote sensing trends

A significant advantage of many REM platforms is their global coverage. Satellites often provide coverage over both land and ocean, and the major ground-based REM network AERONET (Holben et al., 1998) is more globally representative than the sites used in this study. However, there are some inherent limitations in comparing aerosol optical property trends from REM retrievals with surface in situ trends. Our study used aerosol optical measurements made at low RH (typically RH < 40 %) at the surface, while column aerosol optical retrievals are made at ambient conditions and represent the atmospheric column including layers aloft. Only in the situation of a well-mixed atmosphere will it be reasonable to compare trends in surface in situ optical properties with those obtained by ground-based or satellite retrievals. It has also to be mentioned that satellite measurements are less sensitive to the near-ground layers containing the greatest aerosol load. Thus, while our trends can be compared with those for column aerosol properties, there is no reason to expect them to be in complete agreement. Below we discuss trends in PM, AOD, column σ_ap and column SSA.

Satellites have been used to assess the decreasing PM trends in North America and Europe and also to estimate PM trends in other regions with sparse surface measurements. For example, Nam et al. (2017) evaluated the trend in satellite-derived PM₁₀ over Asia and reported mixed annual trend values depending on the subregion they looked at. Li et al. (2017) found satellite-derived PM_2.5 to continuously increase in some parts of Asia (e.g., in India) for the 1989–2013 period – we also find an increasing trend (for aerosol absorption) at the one site we studied in India (MUK). For China, Li et al. (2017) report that the PM_2.5 trend transitions from an increasing to a decreasing trend, with the transition occurring in the 2006–2008 time period similarly to the sulfate trend pattern reported by Aas et al. (2019). The in situ measurements from China (WLG) and Taiwan (LLN) used in our study are not long enough to detect this transition.

Multiple ground-based REM studies (e.g., Yoon et al., 2016; Wei et al., 2019; Mortier et al., 2020) report decreasing trends in AOD over the USA and Europe, with larger decreasing trends over Europe than over the USA, which is the case in our study (see Table 4) as well. The lack of measurements in many regions, similar to the lack of representativeness in the surface in situ aerosol sites discussed in this study (Asia, Africa, South America, etc.), is also emphasized. Ningombam et al. (2019) analyze AOD 1995–2018 trends from 53 remote and high-altitude sites, of which 21 had ss negative trends. Regionally, Ningombam found primarily negative trends at sites in the USA, Europe and the polar regions. Their findings for sites in China and India suggested mixed trends, with some being positive and some negative in those regions. Some of the sites in Ningombam et al. (2019) were also involved in our study. The trends they find for AOD at LLN and MLO are similar to ours (i.e., not ss trends) at SPO (i.e., ss increasing) and at SGP (i.e., decreasing (note: they refer to SGP as “car”)). Their results are different for IZO (we found no ss trends for scattering, while they reported ss decreasing AOD) and for BRW and BIR (we found ss decreasing scattering trends, but they found not ss AOD trends).

Satellite retrievals can offer an even more global picture of aerosol trends than the surface-based REM data. Various satellite trend analyses present a picture of trends in aerosol optical depth for different regions of the world that is quite consistent across satellite (and ground-based) AOD datasets. For example, for the satellite literature that we surveyed, all found decreases in AOD over the USA and Europe (e.g., Hsu et al., 2012; Mehta et al., 2016; Zhao et al., 2017; Alfaro-Contreras et al., 2017; Wei et al., 2019) consistent with what we have reported for the AOD from ground-based, REM instruments. As we note above, this is also consistent with surface in situ scattering trends. There are some discrepancies in the various satellite-derived AOD trends over Asia that are likely due to differences in time period of analysis, trend methodology, regional definitions and/or perhaps satellite data product. Nam et al. (2017) found that AOD trends varied depending on what part of Asia was being evaluated. Zhao et al. (2017) reported an increasing then decreasing trend over China, which was also suggested by others (e.g., Sogachova et al., 2019; Alfaro-Contreras et al., 2017). Wei et al. (2019) found a slightly negative but statistically insignificant AOD trend for China. Our study found statistically insignificant trends in aerosol loading for both the high-altitude surface site in China (WLG) and in Taiwan (LLN), perhaps because measurements at both these sites span the AOD increase/decrease periods mentioned by Zhao et al. (2017). Over India, increasing trends in satellite AOD were reported by all the literature we surveyed (e.g., Wei et al., 2019; Mehta et al., 2016; Hsu et al., 2012; Alfara-Contreras et al., 2017). This is consistent with our finding of an increasing trend for aerosol absorption for the one Indian site (MUK) in our study.

The satellite measurements also enable evaluation of aerosol loading changes in regions with few to no long-term surface in situ aerosol optical property measurements. The Middle East exhibited an increasing trend in AOD, while South America exhibited variable trends (e.g., Wei et al., 2019; Metha et al., 2016; Hsu et al., 2012; Alfaro-Contreras et al., 2017). Wei et al. (2019) found a statistically insignificant trend in South America and suggested it was due to complex and changing aerosol sources. Mehta et al. (2016) looked specifically at Brazil and found a decreasing annual AOD trend but an increasing AOD trend in springtime. Decreasing AOD trends were found over central Africa (Wei et al., 2019), over the African deserts (Metha et al., 2016) and on African coasts (Alfaro-Contreras et al., 2017) regardless of whether they are dominated by smoke aerosols (southwest) or dust (northwest).

In addition to AOD, trends for other column aerosol property such as column σ_ap and column SSA can be considered. While there appear to be many investigations focusing on trends in column aerosol properties other than AOD at individual sites, there are only a few papers that take a more global, multi-site approach (e.g., Li et al., 2014; Zhao et al., 2017; Mortier et al., 2020). There have been several studies related to changes in column σ_ap using AERONET REM retrievals. For example, Li et al. (2014) suggest an increase in column σ_ap over the USA and a decrease over Europe and at most sites in Asia. More recently, Mortier et al. (2020) found ss decreasing σ_ap trends in Europe and North America and ss increasing σ_ap trends in Asia and Africa. Zhao et al. (2017) used satellite retrievals and reported decreasing trends of column σ_ap over both the USA and Europe and a not ss column σ_ap trend over China. Nam et al. (2018) suggested there was an increasing trend in the column extinction Ångström exponent over Asia based on satellite observations. These findings are mostly consistent with our results (Table 4) which indicated decreasing å_sp trends in the USA and Europe but perhaps an increasing trend in Asia.

Comparisons of in situ and column ω₀ trends are more fraught, because, in addition to the above-mentioned caveats related to comparing surface and column measurements, column ω₀ can only be obtained from REM techniques under higher aerosol loading conditions. For example, Kahn and Gaitley (2015) indicate that MISR SSA retrieval requires AOD > 0.15–0.2. Similarly, AERONET retrievals require AOD (at 440 nm) > 0.4 (Dubovik et al., 2000). This limits the sites for which column SSA can be retrieved. Andrews et al. (2017) present a plot derived from global model simulations suggesting more than 80 % of the globe has annual AOD values below 0.2, and, indeed, many of the surface in situ sites discussed here are in remote locations with annual AOD consistently below 0.2. Andrews et al. (2017) also suggest there is a systematic variability of SSA with loading that might result in column SSA biases if retrievals are constrained to higher levels of AOD. With these caveats in mind, we can compare our surface ω₀ trend results with satellite column ω₀ trends. Li et al. (2014) studied 2000–2013 trends in column ω₀ at select AERONET sites. Their findings suggest that column ω₀ is increasing in the USA, Europe and Asia. However, they noted that the uncertainty in these trends is high because they used level 1.5 data (AOD < 0.4) in order to have enough data points for their analysis. Zhao et al. (2017) utilized satellite retrievals and reported decreasing trends in column ω₀ over the eastern USA and Europe and a not statistically significant trend over China for the 2001–2015 period. Their results over the USA and western Europe are consistent with the overall regional ω₀ trends reported in this study (i.e., Table 4), although Fig. 7 suggests there is a fair amount of variability in the surface ω₀ trends at the individual sites in these two regions. Our study found an increasing trend in ω₀ at the surface in Asia (based on three sites), which is consistent with Li et al.'s column ω₀ trend but not with the lack of trend in column ω₀ over China suggested by Zhao et al. (2017). But, as noted above, remote sensing retrievals of column SSA should be considered with caution and, clearly, further effort in column SSA trend analysis is warranted.

4.4.3 Causality

While it is beyond the scope of this effort to explore in depth the causes of the observed trends of aerosol optical properties, some general comments can be made. First, tendencies in regional trends for variables representing aerosol loading (e.g., surface in situ aerosol scattering, PM, and AOD) are generally consistent across multiple datasets. Overall, the main cause of observed decreasing trends in loading is likely a strong reduction of both primary aerosols and precursors of secondary aerosol formation connected to mitigation strategies on regional to continental scales (e.g., Huang et al., 2017; Crippa et al., 2016; Pandolfi et al., 2016; Vestreng et al., 2007). Detailed analysis of PM reductions and composition changes in Europe and the USA has enabled attribution of the trends to changes in source types and emission levels (e.g., Hand et al., 2019; Pandolfi et al., 2016; Ealo et al., 2018). The explanations of the trends based on long-term measurements are supported by modeling efforts. Like many satellite retrievals, model simulations also provide global coverage and, in addition, can be used to investigate reasons for observed changes in aerosol. Model simulations described in Li et al. (2017) suggested that the decrease in PM_2.5 in western and central Europe is principally due to sulfate, whereas in eastern Europe decreases in organic aerosol also play a role. The EMEP status report (2019) notes that the difference in emission trends between western and eastern Europe has become more significant since 2010. Further, the EMEP status report suggests that estimated increasing emissions of all pollutants since 2000 in eastern Europe are mainly influenced by emission estimates for the remaining Asian areas in the EMEP modeled domain. Similarly, Zhao et al. (2019) used a model to attribute the AOD, ω₀ and å_sp decreases in North America and Europe to considerable emission reductions in all major pollutants except in mineral dust and ammonia.

For Asia, modeling by Li et al. (2017) suggests aerosol changes are principally related to increases in organic aerosol and secondary inorganic aerosol, whereas the increases in BC, nitrate and ammonium are comparably moderate. Yoon et al. (2016) use a model to ascribe the observed increases in AOD over India to increases in BC and water-soluble materials – both related to anthropogenic emissions. Over China, Yoon et al. (2016) observe a disconnect between the model chemical composition and the measured AOD, which they explain by noting that the measurement sites they rely on in the region are far from the population centers where most of the emissions occur. Zhao et al. (2019) use a model to attribute the increase in AOD followed by a decrease in AOD to emission increases induced by rapid economic development until 2008–2009 followed by decreases in both anthropogenic primary aerosols and aerosol precursor gases.

Zhao et al. (2017) suggest that the larger reductions in aerosol precursors (e.g., SO₄ and NO_x emissions) rather than primary aerosols, including mineral dust and black carbon, can explain the decreases in ω₀ and å_sp observed over Europe and the USA. This is because the secondary aerosols formed from such precursors tend to be primarily scattering, so less secondary aerosol would change the relative balance between scattering and absorption, driving ω₀ down. Similarly, secondary aerosol particles tend to be small, so a decreasing trend in secondary aerosol would change the relative contribution of small to large particles in the aerosol size distribution and lead to a decreasing trend in å_sp. In contrast, in Asia simultaneous increases in aerosol precursors and BC before 2006 and a simultaneous decrease after 2011 explains the trends ω₀ and å_sp they observed there. Modifications in emissions of aerosol precursors also impact the atmospheric chemistry, leading to non-linear response of the formation of secondary inorganic aerosol (Banzhaf et al., 2015).

While regional changes in emissions are one driving factor in trends, because of long-range transport, out-of-region changes in sources also have the potential to affect trends. For example, Saharan dust impacts CPR, IZO and UGR (e.g., Denjean et al., 2016; Rodriguez et al., 2011; Garcia et al., 2017; Lyamani et al., 2008), and its emissions may change (decrease) in a warmer world (Evan et al., 2016). Other examples of sites clearly impacted by long-range transport include IZO impacted by northern African pollution due to developing industries (Rodriguez et al., 2011) and the high-altitude station of MLO which is impacted by Asian pollution (e.g., Perry et al., 1999). Mountainous stations can also be affected by modifications of the planetary boundary layer or of the continuous aerosol layer heights responding to ground temperature or mesoscale synoptic weather changes (e.g., Collaud Coen et al., 2018, and references therein). The oscillation in trend sign for several variables at the Arctic sites is potentially caused by the very low aerosol loading, but the Arctic region is changing rapidly and the impact of evolving transport patterns, atmospheric removal processes or local sources cannot be excluded (e.g., Willis et al., 2018) and requires closer study. While both increasing and decreasing levels of aerosol due to changes in anthropogenic emissions have been observed, the role of non-anthropogenic sources may become more important in the future. For example, climate change also affects soil drought, and the positive feedback between drought and wildfires can also affect aerosol optical properties (Hallar et al., 2017; McClure and Jaffe, 2018). The number and intensity of wildfires are increasing in several regions (e.g., Moreira et al., 2020; Turco et al., 2018; Hand et al., 2014). McClure and Jaffe (2018) confirmed an increasing trend of PM_2.5 98th percentiles in the northwestern USA due to an increase in wildfires superimposed on the global decrease in anthropogenic emissions. Yoon et al. (2016) also note an increase in extreme AOD events in the western USA, which they hypothesize could be due to wildfires. Another example of potential changes in natural aerosol may take place in the Arctic, where decreases in sea ice coverage might play a role in natural aerosol increases in the region (e.g., Willis et al., 2018) (decreases in sea ice coverage may also lead to enhanced anthropogenic emissions due to increased human activity (e.g., Aliabadi et al., 2015)). Whether such changes in natural aerosol emissions lead to observable changes in overall aerosol trends or trends at the extremes of aerosol loading is something to look for in future trend analyses.

Detailed studies at each station are necessary to discriminate between direct causes like changes in anthropogenic emissions and indirect causes related to general climate changes such as drought, changes in surface albedo, biogenic aerosol concentration, atmospheric chemistry, sea ice coverage or atmospheric circulation patterns. The availability of the homogenized dataset from this study will provide a useful tool for these types of analyses.

In order to get a truly global overview of aerosol trends, surface in situ measurements need to be paired with model simulations and satellite observations. This will enable evaluation of the uncertainty in regional and global trends based on deficiencies in spatial and/or temporal coverage. Satellites and models are able to fill the gaps in coverage from ground-based measurements, but both rely on surface measurements for ground truth.

5 Conclusion and recommendations

This second long-term trend analysis of in situ aerosol measurements derived from stations with large spatial representation leads to a more coherent picture of aerosol radiative properties around the world. Results from this study provide evidence that the aerosol load has significantly decreased over the last 2 decades in North America and Europe. The low number of stations on the other continents means global tendencies cannot be assessed and the results are more variable. The mean extensive property trends are decreasing for all parameters (σ_sp, σ_bsp and σ_ap) and all regions apart from the σ_sp trend in the South Pacific and in the polar regions (see Table 4). These decreases in aerosol burden are assumed to be a direct consequence of decreases in primary particles and particulate precursors such as SO₄ and NO_x due to pollution abatement policies. This assumption is supported by trend results for the USA, where the inflection point between not ss and ss decreasing σ_sp 10-year trends consistently occurred over the same time period (2009–2012) for all central and eastern stations. While the annual σ_ap decrease (−2.5 to −5 % yr⁻¹ for the ss trends in all regions) is larger than that for σ_sp, the σ_ap time series are not long enough to detect the beginning of σ_ap decreasing 10-year trends.

The single scattering albedo trend analysis – derived for the first time from a large number of stations – has the greatest climatic relevance. The uncertainty of the ω₀ trend is higher than for the other aerosol parameters due to uncertainties in absorption coefficient absolute value. The general picture is nuanced, with ss positive trends mostly in Asia and eastern Europe and ss negative in western Europe and North America, leading to an annual positive median trend of 0.02 % yr⁻¹. It appears that the historical abatement policies for gaseous species and primary aerosol particles (e.g., in western Europe in the 1980s) have resulted in present-day decreasing ω₀ trends in the Western Hemisphere, whereas more recent regulations (Asia) are leading to increasing ω₀ trends. Again, this suggests it is necessary to consider how regulatory policies designed to improve health and environmental outcomes impact climate change efforts and vice versa.

The backscattering fraction and scattering Ångström exponent trends relate mostly to the average particle size distribution and to the relative concentrations in the accumulation and coarse modes of the size distribution, but the mean refractive index also plays a role. The interpretation of the results for these parameters is less straightforward as, depending on the site, the trends for b and å_sp may have the same or opposite signs. The causes of particle size change encompass not only the primary aerosol emission, but also the emission of secondary aerosol precursors, the particle chemistry and condensation rate, the hygroscopic growth and the humidity condition during the measurement. In general, the interpretation of b, å_sp and å_ap trends is more difficult, and the effects of global climate change on aridity, wildfire frequency and intensity, planetary boundary layer height, transportation patterns or natural oscillation must also be investigated in order to find the causality of aerosol optical property changes.

This study was limited by the lack of information from many WMO regions. Since 2010, the number of stations with time series longer than 10 years has doubled (24 in 2010, 52 currently), so that the spatial coverage is improved and various additional environments are covered in Europe, North America and the polar regions. A first result of this study is that, while aerosol exhibits a very weak spatial and temporal homogeneity, general features can be deduced with the present station density in Europe and North America, while the picture in the polar regions is less clear. The few stations in Asia, Africa, South America and the Oceania/Pacific region cannot, however, be considered representative of their continents/regions, first, because of their small number and, second, because mountainous and coastal environments are overrepresented relative to the continental environment with rural, forest or desert footprints. According to information from the GAWSIS metadata base, more stations located in underrepresented regions are now in operation, which promises a better spatial coverage in a few years; however, sustaining these operations is still an open issue (the longest time series in India closed in 2016), and not all stations are actually providing their data in open access with the proper associated metadata. Even in developed countries, the financial resources needed to operate long-term monitoring are not always secure, leading to the closing of stations, to a decrease in time series quality and/or to a delay in data submission to the international data banks.

In this study, a number of datasets were not used or were only partially used due to the occurrence of breakpoints following instrumental or inlet changes or even calibration shifts. High-quality data rely on attention to international recommendations for measurements, on a regular maintenance schedule, on participation in inter-comparison efforts and on high-level quality control. The existence of metadata, logbooks and a station's history is crucial for determining causes of any detected breakpoints and is necessary to enable the generation of a final homogenized time series for trend analysis. This homogenization process provides us with an important finding: a critical review of the data by others outside the measurement network is very important in improving the quality of the reported data. This study has resulted in a large improvement to the EBAS database and in the quality of the reported datasets.

Based on the results of this study and with a view toward future trend analyses, the following recommendations concerning the improvement of aerosol optical time series are raised.

• The station history, metadata and logbooks have to be detailed and handled with great care, since they are absolutely necessary to evaluate long-term trends on homogenized time series.

• Time series are affected not only by the instrument type or inlet changes, but also by replacement by instruments of the same type and by shifts in calibrations.

• A rotation between instruments in a network (e.g., to enable repairs) will decrease potential missed data losses but has a potential to increase breakpoints in the time series, particularly in the wavelength dependence of the parameters.

• The scattering and backscattering coefficients, the backscattering fraction and the scattering Ångström exponents are very sensitive to the humidity conditions in the nephelometers due to the hygroscopic growth of particles even at low RH. The nephelometer humidity sensors should be better checked and characterized in order to assess long-term trends of dry particles.

• Long-term trend analysis should not be computed on time series shorter than 10 years, since short datasets lead to a larger probability of false trend detection because of the low number of elements in the time series.

• Stations with long-term records have to be sustained and their funding should be assured in order to study the future impact of aerosol on climate change. Station maintenance as well as new station creation in regions with a low spatial coverage (Africa, South America, Asia and Oceania) should be particularly encouraged.