Evaluation of climate model aerosol trends with ground-based observations over the last 2 decades – an AeroCom and CMIP6 analysis

Abstract

This study presents a multiparameter analysis of aerosol trends over the last 2 decades at regional and global scales. Regional time series have been computed for a set of nine optical, chemical-composition and mass aerosol properties by using the observations from several ground-based networks. From these regional time series the aerosol trends have been derived for the different regions of the world. Most of the properties related to aerosol loading exhibit negative trends, both at the surface and in the total atmospheric column. Significant decreases in aerosol optical depth (AOD) are found in Europe, North America, South America, North Africa and Asia, ranging from $-\mathrm{1.2}$ % yr${}^{-\mathrm{1}}$ to $-\mathrm{3.1}$ % yr${}^{-\mathrm{1}}$. An error and representativity analysis of the spatially and temporally limited observational data has been performed using model data subsets in order to investigate how much the observed trends represent the actual trends happening in the regions over the full study period from 2000 to 2014. This analysis reveals that significant uncertainty is associated with some of the regional trends due to time and space sampling deficiencies. The set of observed regional trends has then been used for the evaluation of 10 models (6 AeroCom phase III models and 4 CMIP6 models) and the CAMS reanalysis dataset and of their skills in reproducing the aerosol trends. Model performance is found to vary depending on the parameters and the regions of the world. The models tend to capture trends in AOD, the column Ångström exponent, sulfate and particulate matter well (except in North Africa), but they show larger discrepancies for coarse-mode AOD. The rather good agreement of the trends, across different aerosol parameters between models and observations, when co-locating them in time and space, implies that global model trends, including those in poorly monitored regions, are likely correct. The models can help to provide a global picture of the aerosol trends by filling the gaps in regions not covered by observations. The calculation of aerosol trends at a global scale reveals a different picture from that depicted by solely relying on ground-based observations. Using a model with complete diagnostics (NorESM2), we find a global increase in AOD of about 0.2 % yr${}^{-\mathrm{1}}$ between 2000 and 2014, primarily caused by an increase in the loads of organic aerosols, sulfate and black carbon.

1 Introduction

As one of the key gears involved in the climate mechanism and as a predominant component of air quality that affects human health , aerosols have been increasingly subject to observation over the last 2 decades, both from ground- and space-based platforms . Aerosols are also recognized to have an important role in the fertilization of the Amazon forest , and in other socioeconomic fields such as solar energy production .

Through their direct, semidirect and indirect effects , aerosol particles are crucial for the estimation of radiative forcing. Currently, the overall estimate of aerosol radiative forcing is associated with high uncertainties . Some of the reasons for these uncertainties reside in the heterogeneity of atmospheric particles, in terms of both their microphysical and their optical properties, as well as in terms of the high variability in these aerosols in space and time. Different regions of the world exhibit contrasting aerosol properties , which can vary depending on the seasons, from year to year, and possibly exhibit interannual trends . In addition to natural emissions such as sea salt and dust, anthropogenic sources of aerosols add another layer of complexity. The second industrial revolution, which relied on the use of fossil fuel energy, has had a significant impact on the aerosol load on a global scale and on the local air quality, resulting in severe pollution episodes, such as the famous 1952 smog event in London that caused the deaths of thousands of people within a few days. Starting in the 1970s in the USA and in the 1990s in Europe, mitigation measures were implemented to limit the emission of particles and other pollutants , resulting in significant improvements in terms of air quality and particle concentration levels . In recent decades there has been a shift in anthropogenic emissions from Europe and North America to the developing nations, which are now facing, in varying degrees, the air quality issues that affected Europe and North America 40 years ago .

In order to provide realistic radiative-forcing estimates and projections, it is important for the atmospheric models to be able to capture the long-term aerosol trends caused by both natural and anthropogenic variations.

Assessing and improving the modeling of aerosols in global Earth system models is the main objective of the AeroCom project. Specific experiments are conducted within this initiative with a focus on individual aerosol species, such as dust or organic aerosols , while dedicated control experiments aim to enable an assessment of global aerosol modeling. Both and, more recently, present evaluations of simulations of global aerosol optical properties by focusing on AeroCom control experiment data for a specific year.

This study presents an overview of the aerosol trends for multiple aerosol parameters (optical and chemical) over the last 2 decades using ground-based observation network data as a reference for the evaluation of the models' skills in reproducing those trends.

To serve that purpose, this study addresses the following two questions:

• What are the observed aerosol trends over the last 2 decades in the different regions of the world (Sect. )?

• Can the climate models reproduce these observed trends (Sect. )?

Then, having developed an understanding of the models' skills in reproducing the observed aerosol trends, the last section of this study aims to answer the following question:

• What are the global aerosol trends derived from the model data (Sect. )?

The CAMS reanalysis dataset and output from six AeroCom models and four CMIP6 models (both model groups performed historical experiments) are evaluated in this work. CMIP6 (Coupled Model Intercomparison Project, Phase 6) is an intercomparison project organized by the WCRP (World Climate Research Programme). Participating models will contribute to the assessment of climate change in the upcoming 2024 IPCC (Intergovernmental Panel on Climate Change) report.

Figure  presents the time series of modeled global AOD (aerosol optical depth) between 1850 and 2014. All of the climate models appear to exhibit a large increase in AOD, especially between 1950 and 1990 , followed by more stable conditions up to the present. While the models show some diversity in absolute values, the trends (focus of this paper) seem to be consistent among the different models at a global scale. Long-term monitoring of many optical and chemical parameters was initiated in the late 1990s e.g.,, providing a high-quality dataset for the investigation of aerosol trends over the last 2 decades e.g.,. These observational datasets also offer an opportunity to validate the modeled trends in this period. Since 2014 is the last year available from the CMIP6 historical runs, we focus this study on the aerosol trends in the 2000–2014 period.

Figure 1

Global AOD computed from model historical runs (Oslo CTM3, GFDL-AM4, CanESM5, CESM2, IPSL-CM6A, ECHAM-HAM) at monthly (gray lines) and yearly (black lines) resolutions, overlaid with the number of active observation sites in the AERONET sun photometer network.

[Figure omitted. See PDF]

2 Datasets

A set of nine columns and in situ surface aerosol datasets are used in this study. The observation networks and the models providing output for these parameters are reported in Table .

Table 1

List of observations and model datasets used in this study (see text for explanation).

Parameter	Type	Observation networks	Models
AOD	Column	AERONET${}^{\mathrm{1}}$	ECMWF-Rean, NorESM2, SPRINTARS, ECHAM-HAM, GEOS, Oslo CTM3, GFDL-AM4, BCC-CUACE, CanESM5, CESM2, IPSL-CM6A
${\mathrm{AOD}}_{\mathrm{f}}$	Column	AERONET	NorESM2, SPRINTARS, ECHAM-HAM, GEOS, Oslo CTM3, GFDL-AM4
${\mathrm{AOD}}_{\mathrm{c}}$	Column	AERONET	ECMWF-Rean, NorESM2, SPRINTARS, ECHAM-HAM, Oslo CTM3, GFDL-AM4
AE	Column	AERONET	ECMWF-Rean, NorESM2, SPRINTARS, ECHAM-HAM, GEOS, Oslo CTM3, GFDL-AM4
${\mathrm{PM}}_{\mathrm{2.5}}$	Surface	EMEP${}^{\mathrm{2}}$, IMPROVE${}^{\mathrm{3}}$	ECMWF-Rean, NorESM2, SPRINTARS, ECHAM-HAM, GEOS
${\mathrm{PM}}_{\mathrm{10}}$	Surface	EMEP, IMPROVE	ECMWF-Rean, NorESM2, SPRINTARS, ECHAM-HAM, GEOS
${\mathrm{SO}}_{\mathrm{4}}$	Surface	EMEP, IMPROVE, CASTNET${}^{\mathrm{4}}$, CAPMoN${}^{\mathrm{5}}$, EANET${}^{\mathrm{6}}$	ECMWF-Rean, NorESM2, SPRINTARS, ECHAM-HAM, GEOS, Oslo CTM3, BCC-CUACE
${\mathit{\sigma }}_{\mathrm{sp}}$	Surface	GAW-WDCA${}^{\mathrm{7}}$ (incl. IMPROVE, NOAA FAN${}^{\mathrm{8}}$, ACTRIS${}^{\mathrm{9}}$, EMEP)	NorESM2, Oslo CTM3
${\mathit{\sigma }}_{\mathrm{ap}}$	Surface	GAW-WDCA (incl. NOAA FAN, ACTRIS, EMEP)	NorESM2, Oslo CTM3

${}^{\mathrm{1}}$ Aerosol Robotic Network. ${}^{\mathrm{2}}$ European Monitoring and Evaluation Programme. ${}^{\mathrm{3}}$ Interagency Monitoring of Protected Visual Environments. ${}^{\mathrm{4}}$ Clean Air Status and Trends Network. ${}^{\mathrm{5}}$ Canadian Air and Precipitation Monitoring Network. ${}^{\mathrm{6}}$ Acid Deposition Network in East Asia. ${}^{\mathrm{7}}$ Global Atmosphere Watch – World Data Centre for Aerosols. ${}^{\mathrm{8}}$ National Oceanic and Atmospheric Administration Federated Aerosol Network. ${}^{\mathrm{9}}$ Aerosol, Clouds, and Trace Gases Research Infrastructure.

2.1 Observations

For each of the parameters used in this study, data of the highest quality level provided by the different observation networks were used. Mountain sites, corresponding to an elevation above 1000 m, were excluded, mainly because global models have problems simulating the aerosol distribution in complex terrain .

2.1.1 Columnar aerosol optical properties

The Aerosol Robotic Network (AERONET) is a network that was established by NASA (National Aeronautics and Space Administration) and has been expanded by national and international collaborations. AERONET operates aerosol ground-based measurements in the different regions of the world . Observation of the columnar aerosol properties is performed by standardized and calibrated solar-powered Cimel Electronique sun photometers. These instruments measure the solar radiation reaching the surface of the Earth at different wavelengths and for different optical geometries. A new version of the sun photometer (CE318-T) is also able to perform nighttime measurements using the moon as a light source . The direct measurements (aiming at the light source) allow for the derivation of the aerosol optical depth (AOD) and the Ångström exponent (AE), which are related to the number and size of the particles, respectively. The spectral information can be further utilized to derive the AOD for the fine and the coarse particles, split by diameters less than or greater than 1.2 $\mathrm{\mu }$m . Three different data quality levels are available depending on the application of cloud filtering and correction for instrument calibration derivations . The level 2.0 version 3 daily data, which provide automatic instrument anomaly quality controls , are used in this study for four different parameters: AOD (calculated at 550 nm); AE (calculated using 440 and 870 nm channels); and ${\mathrm{AOD}}_{\mathrm{f}}$ (or fine AOD) and ${\mathrm{AOD}}_{\mathrm{c}}$ (or coarse AOD), corresponding to the AOD of the particles whose diameter is less than and greater than 1.2 $\mathrm{\mu }$m, respectively.

2.1.2 Particulate matter concentrations

The particulate matter (PM) measurements are from EMEP (covering Europe) and IMPROVE (for North America). The PM data have been made available either via the EBAS database infrastructure (http://ebas.nilu.no, last access: 26 October 2020) or in the original IMPROVE data to be found in the VIEWS database (http://views.cira.colostate.edu, last access: 26 October 2020). Both ${\mathrm{PM}}_{\mathrm{10}}$ and ${\mathrm{PM}}_{\mathrm{2.5}}$ (in $\mathrm{\mu }$g m${}^{-\mathrm{3}}$) are used in this study. Note that the ${\mathrm{PM}}_{\mathrm{10}}$ size fraction of particles below 10 $\mathrm{\mu }$m encompasses the ${\mathrm{PM}}_{\mathrm{2.5}}$ aerosol mass below 2.5 $\mathrm{\mu }$m.

The first PM measurements in EMEP started in 1996, and the number of sites increased steadily in the following decade . Most of the sites use the gravimetric method for both size fractions, though some have used automated monitors, i.e., tapered element oscillating microbalance (TEOM) and filter dynamics measurement system (FDMS) or beta attenuation. The EMEP monitoring complies with European standards, i.e EN 12341:2014 for the gravimetric methods and EN 16450:2017 for the automatic methods.

The IMPROVE network has been operating since 1988 at predominantly remote and rural sites across the United States. This ensures a good representativity of the measurements as some chemical species contributing to PM observations (i.e., organic carbon) can exhibit different seasonality and spatial variability. IMPROVE uses four separate modules to collect samples for speciated ${\mathrm{PM}}_{\mathrm{2.5}}$ analysis and gravimetric ${\mathrm{PM}}_{\mathrm{2.5}}$ and ${\mathrm{PM}}_{\mathrm{10}}$ bulk mass measurements. Samples are collected every third day for 24 h and reported at local conditions. ${\mathrm{PM}}_{\mathrm{2.5}}$ and ${\mathrm{PM}}_{\mathrm{10}}$ mass concentrations are determined from Teflon filters from two separate modules sampling with ${\mathrm{PM}}_{\mathrm{2.5}}$ and ${\mathrm{PM}}_{\mathrm{10}}$ inlets, respectively. The gravimetric mass measurements are not performed at controlled relative humidity and temperature, and a laboratory relocation in 2011 resulted in unstable weighing conditions. Therefore, gravimetric mass measurements from 2011 to 2018 were subject to potentially high-relative-humidity conditions and likely contain particle-bound water on the filters that could bias trends .

2.1.3 Sulfate aerosol concentration

The sulfate aerosol (${\mathrm{SO}}_{\mathrm{4}}$) concentration dataset is a subset of the global data presented in and is based on measurements obtained in different regional networks as described in Table . The sulfate aerosol concentrations are obtained from analysis of aerosol filters. In the EMEP, CASTNET, CAPMoN and EANET networks, these are sampled with either a ${\mathrm{PM}}_{\mathrm{10}}$ inlet or a total aerosol inlet, with no specific-size cutoff effective, using a filter pack sampler. In the IMPROVE network, sulfate measurements are carried out using a filter pack sampler with a ${\mathrm{PM}}_{\mathrm{2.5}}$ inlet. The filters are typically analyzed by ion chromatography after water extraction of the aerosol filter.

The data have been screened to be of satisfactory quality. Urban sites are not included and nor are sites where the surroundings have changed considerably in the period in question. In the data were averaged to monthly means. When the data have a lower sampling frequency than daily, samples are weighted prior to averaging in accordance with how many days were sampled in a given month.

2.1.4 Scattering and absorption coefficients

For the surface in situ PM measurements, the scattering and absorption coefficient measurements were accessed through EBAS database infrastructure. The level 2 data (quality controlled, hourly averaged, reported at standard temperature and pressure – STP – conditions) were used. Detailed information on the quality assurance and quality control procedures for GAW aerosol in situ data is available in . The difference in measurement conditions (i.e., observations being made at STP versus models simulating at ambient conditions) was not expected to impact the calculated trends, so no adjustment was made to account for this.

Scattering and absorption coefficients are measured by different instruments:

• Scattering coefficients (${\mathit{\sigma }}_{\mathrm{sp}}$; in ${\mathrm{Mm}}^{-\mathrm{1}}$) were measured by integrating nephelometers. For better consistency in the model comparisons (model data for these parameters are reported for RH $=$ 0 %), only the measurement data obtained when the relative humidity in the instrument was lower than 40 % were utilized .

• Absorption coefficients (${\mathit{\sigma }}_{\mathrm{ap}}$; in ${\mathrm{Mm}}^{-\mathrm{1}}$) were obtained from filter-based absorption photometers.

Due to the scarcity of stations available for long-term trend analysis (only 28), the presence of regionally nonrepresentative stations (e.g., stations located near roads, in cities), difficult to capture by global models, can have large effects on the computation of the regional-average time series. The urban stations have therefore been removed from this analysis.

Altogether the same data selection procedures (exclusion of stations, removal of outliers) and corrections (conversion to coefficients at 550 nm wavelength) were applied as in , which describes the AeroCom evaluation of the 2019 control experiment, analyzing AeroCom simulations of the year 2010 in detail.

2.2 Models

A set of 10 climate and aerosol models and an aerosol reanalysis dataset are used in this study. Their main characteristics are reported in Table . These models can be separated into three main groups.

Table 2

Information on models used in this study (CAMS-Rean: CAMS reanalysis, AP3: AeroCom phase III, CMIP6: historical experiments from CMIP6).

Model	Group	Natural	Anthropogenic	Meteorology	Resolution	References
		interactive emissions	emissions		(degrees)
ECMWF-Rean	CAMS-Rean	D, SS	MACCity	RA	0.7 $\times$ 0.7
SPRINTARS	AP3	D, SS, DMS, Oce VOCs	C: ${\mathrm{SO}}_{\mathrm{2}}$, BC, OC	N	0.56 $\times$ 0.56
ECHAM-HAM	AP3	D, SS, DMS	C: ${\mathrm{SO}}_{\mathrm{2}}$, BC, OC	fSST	1.875 $\times$ 1.875
GEOS	AP3	D, SS, DMS, Oce VOCs	O: ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{SO}}_{\mathrm{4}}$, BC, OC, ${\mathrm{NH}}_{\mathrm{3}}$	N	1.00 $\times$ 1.00
Oslo CTM3	AP3	D, SS	C: ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{SO}}_{\mathrm{4}}$, BC, OC, ${\mathrm{NH}}_{\mathrm{3}}$	S	2.25 $\times$ 2.25
GFDL-AM4	AP3	D, SS, DMS, Oce & Veg OC	C: ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{SO}}_{\mathrm{4}}$, BC, OC	fSST&N	1. $\times$ 1.25
BCC-CUACE	AP3	D, SS, DMS	C: ${\mathrm{SO}}_{\mathrm{2}}$, BC, OC	F	2.8 $\times$ 2.8
NorESM2	CMIP6	D, SS, DMS, MSA, VOCs	C: ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{SO}}_{\mathrm{4}}$, OC, BC	F	1.89 $\times$ 2.50
CanESM5	CMIP6	D, SS, DMS	C: ${\mathrm{SO}}_{\mathrm{2}}$, OC, BC	F	2.77 $\times$ 2.81
CESM2	CMIP6	D, SS, DMSclim	C: ${\mathrm{SO}}_{\mathrm{2}}$, OC, BC	F	0.94 $\times$ 1.25
IPSL-CM6A	CMIP6	D, SS, DMSclim	*C: ${\mathrm{SO}}_{\mathrm{2}}$, BC, OC, ${\mathrm{NH}}_{\mathrm{3}}$	fSST	2.50 $\times$ 1.27

Anthropogenic emissions (C $=$ CMIP6 CEDS; O $=$ other; *C $=$ CMIP6 modified); natural interactive emissions (D $=$ dust; SS $=$ sea salt; O $=$ biogenic organic; V $=$ volcanic; Oce $=$ oceanic; Veg $=$ vegetation; DMS $=$ dimethyl sulfide; DMS_clim $=$ dimethyl sulfide from climatology; VOCs $=$ volatile organic compounds, MSA $=$ methane sulfonate); meteorology (S $=$ prescribed varying meteorology CTM (chemical transport model); N $=$ GCM (general circulation model) nudged to analyzed meteorology; fSST (fixed sea surface temperature) $=$ fixed SST (sea surface temperature) $-$ SIC (sea ice concentration) monthly fields, GCM not nudged; F $=$ free-running coupled GCM; RA $=$ combined reanalysis of meteorology and composition).

2.2.1 CAMS reanalysis

The CAMS reanalysis, which is the successor to the MACC (Monitoring Atmospheric Composition and Climate) reanalysis, is the latest global reanalysis dataset of atmospheric composition produced by the Copernicus Atmosphere Monitoring Service . It is produced using 4D-Var data assimilation in the CY42R1 model cycle of the ECMWF (European Centre for Medium-Range Weather Forecasts) Integrated Forecast System (IFS), with 60 hybrid sigma–pressure vertical levels. The model used in the CAMS reanalysis includes several updates to the aerosol and chemistry modules on top of the standard CY42R1 release. The IFS model assimilates several satellite products, from aerosols (AOD) to greenhouse gases (CO${}_{\mathrm{2}}$, CH${}_{\mathrm{4}}$; ), whereof most relevant for aerosol trends are data from both MODIS sensors and AATSR and ATSR-2. Daily data, from the ECMWF data archive (MARS), were used in this study. The CAMS reanalysis dataset covers the period January 2003 to near real time. The first 3 years of this study period (2000–2002) are missing for this model.

2.2.2 AeroCom phase III

Initiated in 2000, the AeroCom project (https://aerocom.met.no, last access: 26 October 2020) is an open international initiative of scientists interested in the advancement of the understanding of global aerosols and their impact on climate . Different model experiments have been conducted during the third phase of this project, started in 2015, in order to investigate not only specific topics (e.g., dust, volcanic aerosols, aerosol absorption, hygroscopicity) but also the general modeling of the aerosols (control experiments). The model versions and parametrizations used in AeroCom are closely linked to the versions used for CMIP6 and, for instance, AerChemMIP (Aerosol Chemistry Model Intercomparison Project) climate experiments.

In this study, we use the model outputs from the historical AeroCom experiment, whose main aim is to understand the regional trends in aerosol distribution from 1850 to 2015 and to quantify the aerosol forcing with a main emphasis on the direct aerosol effect. The models were also run in various configurations, providing different degrees of constraints on the evolving meteorological conditions, such as using monthly fixed sea-surface temperature (SST); historically evolving SSTs; and basic meteorology fields, e.g., wind for a given year.

2.2.3 CMIP6

The upcoming 2024 IPCC sixth assessment report (AR6) will feature new state-of-the-art CMIP6 models with model runs in higher resolution and with new physical processes. An overview of the experimental design and organization can be found in . In this study, we use a preliminary extract of the data of four CMIP6 models from the historical experiment, as available on ESGF (Earth System Grid Federation) nodes (https://esgf.llnl.gov, last access: 26 October 2020), which provided output from 1850 to 2014. The last year of the study period of the analysis presented here was selected as 2014.

3 Methods

3.1 Regional time series

Due to the nature of the processes involved in the emission and the deposition of aerosols, one can expect different trends in different regions of the world. Instead of investigating the trends obtained at each individual observation station in a given region, we resort here to the analysis of average regional time series as computed by assembling all measurements at stations in each region. One advantage of this method is that a single trend can be computed in a given region, with an associated significance and uncertainty. It is difficult, apart from in a diversity analysis, to define such an uncertainty when combining individual trends. Also, with our aggregation method, even a station that has not provided a sufficient number of data for computing a trend at its location can still contribute to the computation of a regional time series. The computation of such aggregated regional time series makes most sense in regions exhibiting similar seasonal patterns.

3.1.1 Region definition and observation coverage

Seven regions are considered in this study. The definition of these regions has been determined in a pragmatic way to limit the number of geographic areas investigated, but altogether it also provides global coverage when considering the ensemble of all regions. The Americas and Africa have been separated into northern and southern sections. In order to assemble the sites most affected by Saharan dust, the North Africa region has been extended to the north beyond the continent itself. Stations located in the south of Spain, Cyprus and Greece contribute to the regional time series in the region we are calling North Africa. The regions' coordinates can be found in the Supplement.

As seen in Fig. , the regions do not have similar coverage in terms of observations. North America and Europe have the highest concentrations of instruments monitoring aerosol trends.

• AERONET is the most important network in terms of number of instruments. More than 1000 observation points, with more or less long time series, are found across the globe. The highest density of instruments is in Europe and in the central part of North America (USA). The lowest densities are found in southern Africa and Australia.

• For particulate matter, 227 instruments are used in this study and are spread mostly over Europe and North America.

• For ${\mathrm{SO}}_{\mathrm{4}}$, 346 instruments altogether have been operating, mostly in North America and Europe. A few stations are also located in Asia and North Africa.

• For both parameters ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$ combined, approximately 50 stations are spread over North America, Europe, North Africa and Asia. Due to time coverage issues (2005 is the first year where in situ optical data are available in the European time series), the data from 2000 up to the year 2018 were used to compute the regional time series of these two parameters.

Figure 2

Distribution of the observations within the different regions considered in this study. The numbers reported within each region correspond to the maximum number of stations given for the observation networks, corresponding to the five observation types found in the legend.

[Figure omitted. See PDF]

3.1.2 Time series aggregation requirements

The regional time series are computed by combining, for each month, the valid data of all the stations in the corresponding region. In order to construct consistent and robust regional time series, some additional criteria are required to be met to provide a valid point (a station with valid measurements) for the regional time series. Stations with very short time series (e.g AERONET DRAGON campaign stations) are eliminated by requiring a minimum of 300 valid daily measurements in the whole period from 2000 to 2014, which reduces, as an illustration, the number of AERONET stations from 1015 to 437. A minimum of three valid stations is required to be present in the overall regional time series to produce a valid point. In other words, if the available time resolution is daily, at least three stations need to provide valid data for a certain day in order to produce a valid regional mean for that day. The list of the station names contributing to the computation of the regional time series can be found in the Supplement.

When all criteria are fulfilled for a given month in the regional time series, the median and the first and third quartiles are computed from all valid data points available. The quartiles provide an indication of the intraregional variability. An example of regional time series is shown in Fig.  for AOD.

Figure 3

Regional time series of AOD. The dark blue line corresponds to the median, and the light blue envelope is bound by the first and third quartiles of all valid points at the corresponding month, respectively. The blue dots correspond to the yearly averages which are used to compute the linear trend. The latter is displayed as a continuous line when the trend is significant and as a dashed line when it is not. Trend values, an error estimate and a significance value are given in each panel.

[Figure omitted. See PDF]

3.2 Trends calculation

3.2.1 Yearly regional time series

For all of the parameters, the trends are computed based on the yearly averages of the regional time series. Using the yearly averages eliminates any issues caused by the seasonal cycles (observed for most of the aerosol parameters used in this study) during the calculation of the trend slope. In order to ensure the statistical robustness of these yearly averages, the time averaging is performed step-by-step with specific time constraints. By starting at the finest time resolution available in the data, monthly, seasonal and then yearly averages are computed when the following criteria are fulfilled:

• at least 5 d per month are available (when daily observations are available)

• at least 1 month with data per season is present (seasons defined as JFM, AMJ, JAS, OND)

• all four seasons are available for a given year.

These temporal constraints offer a reasonable compromise between the availability and robustness of the yearly statistics.

3.2.2 Trends computation

We use the same methodology as described by to derive the trends in the regional time series. The significance of the trends is tested with the Mann–Kendall test . The related $p$ value is used to determine if the trend is significant or not within a confidence interval of 95 %. The slope is calculated with the Theil–Sen estimator, which is less sensitive to outliers than standard least-squares methods . At least seven valid yearly regional averages (50 % of time coverage) are required in the regional time series for the computation of a slope.

An uncertainty is provided for each trend by combining the error in the slope calculation itself with the error in the residuals:

1$$\mathrm{Uncertainty}=\sqrt{{\left({\displaystyle \frac{\mathrm{\Delta }m}{y\left(\mathrm{2000}\right)}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{m\cdot \mathrm{\Delta }r}{y(\mathrm{2000}{)}^{\mathrm{2}}}}\right)}^{\mathrm{2}}},$$ where $\mathrm{\Delta }m$ is the Theil–Sen estimator 95 % confidence interval, $y\left(\mathrm{2000}\right)$ is the value of the regression line at the year 2000, $m$ is the value of the Theil–Sen slope, and $\mathrm{\Delta }r$ is the averaged error in the residuals computed based on the difference between the linear trend and the yearly mean values of the regional time series.

The trend is provided as a relative trend (% yr${}^{-\mathrm{1}}$) with respect to the first year of the time period (2000).

3.3 Representativity of the trends

The number of available points used to compute the regional time series is not constant in time. For a given observation station, the number of points available might vary in time due to the nature of the measurements. For instance, classic sun photometers only measure in the daytime and in cloud-free conditions. Due to seasonal daylight and cloud condition variations, clear seasonal cycles are observed in the number of observations of AOD. The density of the different observation networks can also change with time. The early development of the different observation networks usually coincided with an increase in the number of observation stations. More recently, primarily for funding reasons, some networks have reduced the number of stations. This variation in the number of available measurements raises the question of time representativity for the computation of the trends.

Associated with this time representativity issue is the space representativity issue. The data coverage is uneven across the different regions. Moreover, within a single region, the observation stations might be located in contrasting environments. Stations located in environments that are more urban, more rural or mostly affected by natural particles might have trends differing from the trend associated with the whole region.

Some studies have focused on the representativity of the observation stations by investigating the biases of different optical properties . The analysis here is dedicated to characterizing the representativity of the observation networks specifically for the purpose of computing the trends. These two perspectives on representativity might give different results, since a station associated with a bias could still have a representative tendency in time. In order to evaluate the effect of the partial space and time sampling of the observations for the evaluation of the trends, two sensitivity studies, focusing on the time sampling and the space sampling, have been conducted using NorESM2 model data subsets. For each of these studies, the trends are computed for one reference (Ref) and one experiment (Exp) dataset and compared with each other. The reference dataset corresponds to the model data co-located with the available observations, while the experiment dataset uses all model points.

• Time representativity study.

  • –

    For Ref${}_{\mathrm{time}}$, model data are co-located in space and time with available observations.

  • –

    For Exp${}_{\mathrm{time}}$, model data are co-located in space with available observations using the complete model time series.

• Space representativity study.

  • –

    For Ref${}_{\mathrm{space}}$, model data are co-located in space with available observations using the complete model time series ($=$ Exp${}_{\mathrm{time}}$).

  • –

    For Exp${}_{\mathrm{space}}$, all of the model grid points in the region use the complete model time series.

The difference between the relative trends is computed for each parameter and region. In order to summarize the representativity, those differences are then converted into a score ($\mathrm{\%}$) by using a mapping function which has been defined based on a normal distribution. The choice of the parameters describing this function leads to a representativity score of 100 % when there is no difference in the trends computed for a reference and an experiment dataset, while a difference of 0.5 $\mathrm{\%}{\mathrm{yr}}^{-\mathrm{1}}$ obtained with these two datasets would indicate a representativity score of 50 %. Finally, the total score is computed as the mean of the time and the space representativities.

An example of the calculation is presented in Fig.  for AOD in Europe and North America. In both regions, the Ref${}_{\mathrm{time}}$ dataset, corresponding to the available observations, reveals strong seasonal cycles when considering the number of points used to compute the regional time series. These cycles are observed with most of the sun photometer datasets since the instruments only operate during daytime and cloud-free conditions and the amount of daylight and cloud varies with the season. Together with this seasonal cycle, one observes an increase in the number of points with time, which reflects the increasing number of stations over these two regions.

Figure 4

Three regional AOD time series and respective trends, constructed from model data (NorESM2) for the investigation of the representativity of the observational data. Panels (a) and (b) correspond to the number of points used to compute the regional time series for the three different datasets. Panels (c) and (d) show the time series, the trends and the resulting representativity value (black, bold). The blue color (Ref${}_{\mathrm{time}}$) corresponds to the model output co-located in space and time with the available observations. Panels (a) and (b) show an overall increase in the number of available observations (more stations) combined with a seasonal cycle (less AOD available in wintertime). The orange color (Exp${}_{\mathrm{time}}/{\mathrm{Ref}}_{\mathrm{space}}$) corresponds to the model output co-located in space with the stations providing measurements, using the complete time series from 2000 to 2014. The green color (Exp${}_{\mathrm{space}}$) corresponds to the model output in the whole geographic region (see Fig. ), using all of the grid boxes without any co-location with the observations.

[Figure omitted. See PDF]

The trends in Europe show similar values for the time study, which means that the trend is not greatly affected by the variation in the available measurements in time. The difference is larger when considering all the grid boxes of the domain, but the overall difference in the two studies corresponds to a representativity of 76 %. In North America, the difference in the three trends is larger, particularly for the space study trend. This means that the trend obtained in the whole region is significantly different from the trend obtained when considering only the grid points where observation stations are located. It should be mentioned that the ocean grid points are not filtered out when computing the trends over the whole domain. For this reason, the regions containing a greater proportion of ocean grid points, where the trends are most likely to differ from those observed over land, will tend to have a lower spatial representativity, such as in North America.

This representativity study illustrates that the partial coverage in time and space of the observations leads, in some cases, to artificial trends. The representativity scores are discussed for each parameter in the following section together with the trend estimate results.

4 Results

4.1 Trends in observations

This section presents the trends in the observations computed for the different parameters and over the predefined regions. In order to compare the trends observed for the set of nine aerosol parameters in a consistent manner, we focus on the relative trends, with the reference set to the year 2000 as the first year of the study period. The means for the year 2000, reported in Table , reveal a large interregional variability.

Table 3

Observational mean values for the year 2000, the reference year used for computing relative trends. Each value is extracted as the intercept of the linear trend computed in the 2000–2014 period, except for ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$, where the trends have been computed over 2000–2018. Because the required minimum number of yearly averages was set to seven, no trend could be computed in the southern African region.

	Europe	North America	South America	North Africa	Asia	Australia
AOD	0.17	0.10	0.15	0.26	0.37	0.10
${\mathrm{AOD}}_{\mathrm{f}}$	0.15	0.08	0.11	0.11	0.22	0.04
${\mathrm{AOD}}_{\mathrm{c}}$	0.03	0.02	0.03	0.10	0.09	0.03
AE	1.43	1.48	1.26	0.70	1.16	1.00
${\mathrm{PM}}_{\mathrm{2.5}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$)	12.3	6.9	–	9.0	–	–
${\mathrm{PM}}_{\mathrm{10}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$)	16.8	12.4	–	19.7	–	–
${\mathrm{SO}}_{\mathrm{4}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$)	2.01	1.45	–	2.98	1.97	–
${\mathit{\sigma }}_{\mathrm{sp}}$ (${\mathrm{Mm}}^{-\mathrm{1}}$)	34.4	23.4	–	–	–	–
${\mathit{\sigma }}_{\mathrm{ap}}$ (${\mathrm{Mm}}^{-\mathrm{1}}$)	6.3	2.6	–	–	–	–

The AOD is more than 3 times higher in Asia (AOD $=$ 0.37) than in North America and Australia (AOD $=$ 0.10). Intermediate AOD values are found in Europe and South Africa, while the second-highest load is found in North Africa (AOD $=$ 0.26). In most regions, the AOD is largely dominated by its fine-mode fraction (${\mathrm{AOD}}_{\mathrm{f}}$), but this is not the case in North Africa (or Australia), where the persistent presence of desert dust makes the coarse-mode (${\mathrm{AOD}}_{\mathrm{c}}$) contribution to the total AOD similar in size to the fine-mode contribution. This predominance of coarse particles is reflected in the AE values which exhibit lower values in North Africa (AE $=$ 0.70) and Australia (AE $=$ 1.00).

The PM observations are primarily available from Europe and North America. ${\mathrm{PM}}_{\mathrm{10}}$ observations are also available in the North Africa region as defined in this analysis, but these stations are located in the northern part of the region, i.e., in southern Europe, which is less affected by the dust sources than the AERONET stations, which cover the whole region including the surrounding deserts. Both ${\mathrm{PM}}_{\mathrm{10}}$ and ${\mathrm{PM}}_{\mathrm{2.5}}$ are larger in Europe than in North America, with different relative proportions. In Europe, ${\mathrm{PM}}_{\mathrm{2.5}}$ represents 76 % of the ${\mathrm{PM}}_{\mathrm{10}}$, as compared to only 57 % in North America. This difference in the relative proportion of fine particles against coarse particles in Europe and North America may be due in part to our definition of regions. presented a phenomenology of PM data in Europe showing coarse aerosols tended to be highest in southern Europe which in our study is part of the North Africa region. The discrepancy in the relative proportions of coarse and fine aerosols in Europe and North America may be exacerbated both by a decrease in North America of the fine-particle concentration due to pollution mitigation strategies and by the growth of the coarse mass due to increasing contributions of natural and agricultural sources, particularly in the western half of the USA .

${\mathrm{SO}}_{\mathrm{4}}$ means (surface mass concentrations) for the year 2000 range between 1.45 and 2.98 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$ with the lowest value occurring in North America and the highest value in North Africa (sites in southern Europe). Similar means are found in Europe and Asia, around 2 $\mathrm{\mu }$g m${}^{-\mathrm{3}}$, though one should bear in mind that there are relatively few sites in Asia and they are not located in the most polluted areas in China and India .

Analogous to the surface ${\mathrm{PM}}_{\mathrm{10}}$ measurements, ${\mathit{\sigma }}_{\mathrm{sp}}$ is higher in Europe (34 ${\mathrm{Mm}}^{-\mathrm{1}}$) than in North America (23 ${\mathrm{Mm}}^{-\mathrm{1}}$). The same feature is found for ${\mathit{\sigma }}_{\mathrm{ap}}$, which also has higher values in Europe than North America.

The relative trends for the 2000–2014 period are shown in Fig. . The heatmap is dominated by blue colors, which indicate mostly negative trends, especially when considering the parameters related to aerosol burden (i.e., the extensive parameters). Usually, the lowest $p$ values ($\mathit{<}$ 0.05) are associated with the lowest uncertainties, although these are not shown in the same figure. The largest circles (highest significance of trend) are more confidently associated with a decrease or increase in the aerosol property in the time period 2000–2014 since the value of the trend is greater than the uncertainty. The uncertainties are presented in Fig. . Among the 38 computed trends, 22 are associated with a representativity score higher than 50 % and 24 are significant at a 95 % confidence level.

In Europe, both columnar and surface parameters reveal statistically significant decreases. With the exception of ${\mathit{\sigma }}_{\mathrm{ap}}$, for which the associated uncertainty in the trend exceeds the trend itself, the trends computed for other parameters are associated with uncertainties lower than the trends values. A decrease in AOD ($-\mathrm{2.8}$ % yr${}^{-\mathrm{1}}$) is found for both fine- and coarse-mode particles. This is consistent with the negative trends found at some individual stations in this region . The fine mode decreases more than the coarse mode, which is consistent with the decrease observed for the AE. The same shift in aerosol size is found at the surface since ${\mathrm{PM}}_{\mathrm{2.5}}$ has decreased by a factor of 2 relative to ${\mathrm{PM}}_{\mathrm{10}}$. These trends could result from the mitigation measures aimed at reducing anthropogenic aerosol emissions. This is more directly observed in the decrease in ${\mathrm{SO}}_{\mathrm{4}}$ ($-\mathrm{1.5}$ % yr${}^{-\mathrm{1}}$). We find a somewhat lower trend than what was reported in ($-\mathrm{2.67}$ % yr${}^{-\mathrm{1}}$), but that could be explained by the differences in the methodology (trends computed from the regional time series, in this study, against a statistical average of the trends computed at the individual stations) and/or the definition of the region. The stations in the Mediterranean Basin, where a larger decrease is found ($-\mathrm{4.3}$ % yr${}^{-\mathrm{1}}$), are attributed to the North African region in this study.

The representativity study reveals that the observed trends are actually representative of the whole period and region for all of the parameters. A good agreement is found with the trends obtained at individual stations reported by , who found decreases of $-\mathrm{2.92}$ % yr${}^{-\mathrm{1}}$ for ${\mathit{\sigma }}_{\mathrm{sp}}$ and $-\mathrm{4.2}$ % yr${}^{-\mathrm{1}}$ for ${\mathit{\sigma }}_{\mathrm{ap}}$, as compared to $-\mathrm{2.5}$ % yr${}^{-\mathrm{1}}$ and $-\mathrm{2.0}$ % yr${}^{-\mathrm{1}}$ in this study.

In North America, similar trends to those in Europe are found for the columnar properties as were found for Europe. AOD decreases at a rate of 1.3 % yr${}^{-\mathrm{1}}$, a 55 % percent smaller trend than observed in Europe, but the North America reference value in 2000 is 40 % lower than the reference value in Europe. One can note that the representativity scores are higher for the AE than for AOD, while these two parameters have the same number of data. This means that the trends in the AE are probably more homogeneous in space and time, which makes the same number of available observations more representative in the case of the AE. The decreases observed for both ${\mathrm{PM}}_{\mathrm{2.5}}$ ($-\mathrm{2.0}$ % yr${}^{-\mathrm{1}}$) and ${\mathrm{PM}}_{\mathrm{10}}$ ($-\mathrm{1.2}$ % yr${}^{-\mathrm{1}}$) are significant and in the same range of values as the trends found in Europe. However, the actual trends in ${\mathrm{PM}}_{\mathrm{10}}$ and ${\mathrm{PM}}_{\mathrm{2.5}}$ are probably somewhat more negative than those found here. The possible bias is caused by increased relative humidity during weighing, resulting in more particle-bound water and thus higher mass, after the relocation of the laboratory in 2011. reported that the decrease in ${\mathrm{PM}}_{\mathrm{2.5}}$ from 2005 through 2016 was $-\mathrm{2.6}$ % yr${}^{-\mathrm{1}}$, while it was $-\mathrm{3.9}$ % yr${}^{-\mathrm{1}}$ for the reconstructed fine mass, correcting for the possible bias in the measurements. ${\mathrm{SO}}_{\mathrm{4}}$ decreases by about 3 % yr${}^{-\mathrm{1}}$, which is twice as large as the decrease observed in Europe, where the reference value is however larger than in North America. The sulfate trend is similar to the trend reported by in this region ($-\mathrm{3.15}$ % yr${}^{-\mathrm{1}}$). The regional time series extend farther back in time for ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$ in North America than in Europe. However, no significant trends are found for these datasets. This is in contrast to , who found a large decrease for ${\mathit{\sigma }}_{\mathrm{sp}}$ ($-\mathrm{2.66}$ % yr${}^{-\mathrm{1}}$). Our study used regionally averaged time series to calculate the trend rather than regionally averaged trends as used by . This probably illustrates the difference in methodology which consists of computing the mean of station trends in one case and the trend in a regional time series in the other case, especially when only a few measurements are available. However, as shown by the representativity study (Fig. ), the nonsignificant increase of $+\mathrm{0.0}$ % yr${}^{-\mathrm{1}}$ found in this study with the observations is similar to the trend derived over the whole region and using complete time series of the NorESM2 model data. Similar values are found in this study and by for ${\mathit{\sigma }}_{\mathrm{ap}}$ ($-\mathrm{1.85}$ % yr${}^{-\mathrm{1}}$) although the trend is, here, not significant. The IMPROVE network also measures filter absorption using a hybrid integrating plate and sphere (HIPS) system . These data are not included in this study, but report a significant decrease ($-\mathrm{2.7}$ % y${}^{-\mathrm{1}}$) in the light absorption coefficients from 2005 to 2015.

All of the columnar properties, except for ${\mathrm{AOD}}_{\mathrm{c}}$, show significant decreasing trends in South America. As shown in the regional time series in Fig. , the observed decrease in AOD coincides with a global diminution of the intensity of the seasonal peaks happening around September and resulting from the Amazonian forest fires . These peaks are highly variable from year to year and could greatly affect the trend when considering another time period. With a rate of $-\mathrm{2.0}$ % yr${}^{-\mathrm{1}}$, the largest decrease in the AE is found in this region. While no significant trend is found for ${\mathrm{AOD}}_{\mathrm{c}}$, the tendency towards increasing coarse particles is probably due to the production of local dust as a result of the increasing deforestation .

In North Africa, while significant decreases are found for all AOD parameters, an increase in the AE ($+\mathrm{1.2}$ % yr${}^{-\mathrm{1}}$) is observed, which indicates an increase in the proportion of fine particles with time. This is consistent when considering the AOD of the fine and coarse modes, which reveals a larger decrease for ${\mathrm{AOD}}_{\mathrm{c}}$. also found a decrease in dust in the Sahara and Sahel in the time period 1980–2009 due to reduced 10 m wind speed, possibly caused by an increase in sea-surface temperature (SST) in the North Atlantic.

The AE also increases in Asia as a combination of a (not significant) increase in ${\mathrm{AOD}}_{\mathrm{f}}$ and a significant decrease in ${\mathrm{AOD}}_{\mathrm{c}}$. The increase in the AE is likely tied to increases in anthropogenic emissions which are associated with fine-mode aerosols. This result is consistent with the trend reported by at some individual stations. At the same time, we observe an increase in ${\mathrm{SO}}_{\mathrm{4}}$ of 3.8 % yr${}^{-\mathrm{1}}$, which is consistent with the trend reported in . This increase is associated with a large uncertainty ($\pm$ 4 % yr${}^{-\mathrm{1}}$) due to a drop in the already small number of stations available in the region, especially between 2010 and 2012. Indeed, with a maximum of 12 stations, a few stations missing can greatly affect the computation of the regional time series. This is reflected by the representativity study which reveals a score lower than 40 % for this parameter.

No significant trends are found in Australia, although the representativity is greater than 50 % for ${\mathrm{AOD}}_{\mathrm{f}}$.

Figure 5

Regional trends in the aerosol properties computed with the observation datasets. The color of the circles corresponds to the slope, while the radius indicates the $p$ value. The largest circles represent the trends that are significant with a confidence level of 95 %. The circles outlined in black indicate the trends associated with a representativity greater than 50 %.

[Figure omitted. See PDF]

This multiparameter trend analysis reveals a decrease in most of the parameters relating to aerosol burden, both in the total column and at the surface level. In Asia, the trends in ${\mathrm{AOD}}_{\mathrm{f}}$, the AE and ${\mathrm{SO}}_{\mathrm{4}}$ suggest an increase in the proportion of the finer particles. While differences might be expected when comparing regional trends with trends computed at individual stations, the trends are usually consistent with those previously reported in the literature. The work of focused on regional AOD trends in the 2000–2009 period; despite the differences in the study periods and the methodologies involved, trends consistent with those found in this study are found for most of the regions.

4.2 Evaluation of the models trends against observations

In order to evaluate the trends from the models, the regional time series have been computed with the model output co-located in space and time with the available observations at the station level. The model trends are computed in a similar manner to the trends for the observation datasets. The results, shown in Fig. , reveal (a) the differing abilities of the various models to reproduce the observed trends and (b) that the model performance depends on the parameters and the regions.

• For AOD, the models show trends with the same sign as the observed trends over all the regions except in Asia, where the associated uncertainties are usually larger than the trend values. Some differences among the three model groups are observed when investigating the different regions:

  • –

    Europe. All the groups underestimate the observed decrease in AOD. With an average decrease of $-\mathrm{1.0}$ % yr${}^{-\mathrm{1}}$, the CMIP6 models exhibit the largest underestimation, while the best performance is obtained with CAMS reanalysis (CAMS-Rean; $-\mathrm{2.0}$ % yr${}^{-\mathrm{1}}$). The AeroCom phase III (AP3) models' trends range from $-\mathrm{1.3}$ % yr${}^{-\mathrm{1}}$ to $-\mathrm{2.3}$ % yr${}^{-\mathrm{1}}$.

  • –

    North America. In contrast to the results for Europe, on average, all of the models overestimate the observed AOD decrease in North America even though two models of the AP3 group simulate lower trends than are found for observations. The consistency in the trends is very high within the CMIP6 group over this region.

  • –

    South America. CAMS-Rean slightly overestimates the observed AOD decrease, while almost all of the models in the two other groups underestimate this decrease. A few of the models simulate positive trends, but these are associated with large uncertainties.

  • –

    North Africa. All the models capture the observed decreasing AOD tendency. With a trend of $-\mathrm{3.0}$ % yr${}^{-\mathrm{1}}$, CAMS-Rean is the closest to the observed trend ($-\mathrm{2.8}$ % yr${}^{-\mathrm{1}}$). AP3 and CMIP6 multimodel trend averages are both equal to $-\mathrm{2.1}$ % yr${}^{-\mathrm{1}}$.

  • –

    Asia. A large intermodel variability is found in this region where the uncertainty is also important. The means of the trends in each group are close to 0 % yr${}^{-\mathrm{1}}$.

• For ${\mathrm{AOD}}_{\mathrm{f}}$, usually, the same patterns are found as for AOD. The models that underestimated the AOD underestimate ${\mathrm{AOD}}_{\mathrm{f}}$ and vice versa. For ${\mathrm{AOD}}_{\mathrm{f}}$ and the following parameters, only NorESM2 provides data for the CMIP6 group.

  • –

    Europe. The underestimation of the decrease in ${\mathrm{AOD}}_{\mathrm{f}}$ captured by the models is larger than the underestimation of AOD.

  • –

    Asia. As for AOD, the trends are associated with large uncertainties and show a large intermodel variability.

• For ${\mathrm{AOD}}_{\mathrm{c}}$, the performance of the models is not as good as for ${\mathrm{AOD}}_{\mathrm{f}}$. This is also observed when evaluating the models for a single year . The intermodel variability is also higher since some models simulate ${\mathrm{AOD}}_{\mathrm{c}}$ trends in opposite directions in some regions.

  • –

    Europe. While the observations exhibit a significant decrease, CAMS-Rean and all of the AP3 models exhibit increasing values for ${\mathrm{AOD}}_{\mathrm{c}}$. NorESM2 from CMIP6 simulates a decrease consistent with the observations.

  • –

    South America. All of the models simulate large increases, from $+\mathrm{1.2}$ % yr${}^{-\mathrm{1}}$ up to $+\mathrm{8.5}$ % yr${}^{-\mathrm{1}}$, which are not visible in the observations ($+\mathrm{0.1}$ % yr${}^{-\mathrm{1}}$).

  • –

    North Africa. The models reproduce the observed decrease of 3.1 % yr${}^{-\mathrm{1}}$ to some extent (from $-\mathrm{0.6}$ % yr${}^{-\mathrm{1}}$ to $-\mathrm{2.7}$ % yr${}^{-\mathrm{1}}$). The fact that some models with a fixed SST (e.g ECHAM-HAM) reproduce this decrease does not support the hypothesis of the SST changes impacting dust emissions .

  • –

    Asia. CAMS-Rean captures a similar negative trend to that computed with the observation dataset. As with ${\mathrm{AOD}}_{\mathrm{f}}$, no certain trend can be identified in this region with the NorESM2 CMIP6 model.

• For the AE, the trends are usually smaller than for AOD in the respective regions. This can mean that the number of particles is more subject to variations than the size (type) of these particles but could also illustrate that the AE is less sensitive to the change in a relative sense. This feature is visible with both observations and models.

  • –

    Europe and North America. One model in the AP3 group (ECHAM-HAM) simulates a significant positive trend in the AE, while negative tendencies are found in the observation and with the other models.

  • –

    South America. All of the models simulate negative AE trends, most of them significant, in agreement with the observations. CAMS-Rean and the AP3 models tend to underestimate the decrease, while NorESM2 CMIP6 model tends to overestimate it.

  • –

    North Africa. CAMS-Rean does an excellent job of reproducing the observed AE increase ($+\mathrm{1.3}$ % yr${}^{-\mathrm{1}}$ versus $+\mathrm{1.2}$ % yr${}^{-\mathrm{1}}$). The significant trends in the AP3 models range from $-\mathrm{0.4}$ % yr${}^{-\mathrm{1}}$ to $+\mathrm{2.0}$ % yr${}^{-\mathrm{1}}$. The increase in the AE supports the theory of enhanced scavenging of dust by anthropogenic aerosols.

  • –

    Asia. The AP3 models and the NorESM2 CMIP6 model exhibit significant positive trends, which is also the case for the observations. CAMS-Rean does not capture any significant trend in this region.

• For ${\mathrm{PM}}_{\mathrm{2.5}}$, almost all the models simulate significant decreases over Europe and North America, in good agreement with the observations. The CMIP6 model performs better however in North America, while it underestimates the extent of the decrease in Europe. Further analysis reveals that, despite the fact that it does a good job reproducing the ${\mathrm{PM}}_{\mathrm{2.5}}$ trend in North America, CAMS-Rean exhibits a large positive bias in this region when considering the absolute values ($+\mathrm{100}$ %). In North Africa, both CAMS-Rean and AP3 models capture the significant decrease seen in the observations.

• For ${\mathrm{PM}}_{\mathrm{10}}$, in North Africa, only CAMS-Rean reproduces the observed significant decrease. Positive trends are found for all the models in the AP3 (except GEOS) and CMIP6 groups. As for ${\mathrm{PM}}_{\mathrm{2.5}}$, NorESM2 has better performance in North America. CAMS-Rean produces a trend that is twice as high as the observed trends over both Europe and North America.

• For ${\mathrm{SO}}_{\mathrm{4}}$, the AP3 and CMIP6 models perform quite well for the ${\mathrm{SO}}_{\mathrm{4}}$ surface concentration. The magnitude of the model trends is however higher than the observed trends in all the regions except North Africa.

• For ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$, as mentioned in the previous section, the observation trends have been computed for these two parameters using data up to 2018. The two models providing output for these parameters are NorESM2 and Oslo CTM3. NorESM2 provides data until 2014, so the NorESM2 trends correspond to the period (2000–2014), while Oslo CTM3 provides data until 2017 and the respective trends correspond to (2000–2017).

  • –

    Europe. A significant decrease is found in the observations for both ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$ but this is not captured by the models for which the calculated trends are associated with large uncertainties.

  • –

    North America. A significant decrease of $-\mathrm{1.3}$ % yr${}^{-\mathrm{1}}$ is found with both NorESM2 and Oslo CTM3 for ${\mathit{\sigma }}_{\mathrm{sp}}$ which is not seen in the observations. For ${\mathit{\sigma }}_{\mathrm{ap}}$, the two models capture a similar trend as derived from the observations ($-\mathrm{1.5}$ % yr${}^{-\mathrm{1}}$).

Figure 6

Regional trends in the aerosol properties computed with observations and models co-located in space and time with the observations. The error bars correspond to the uncertainty in the trend as calculated using both the uncertainty in the Theil–Sen slope and the residuals. The bold font indicates that the trends are significant at a confidence level of 95 % ($p$ value $\mathit{<}$ 0.05).

[Figure omitted. See PDF]

This model trends evaluation reveals some key points. First, CAMS-Rean, which assimilates AOD, performs the best for capturing the trends of this parameter. Second, a large intermodel variability is generally found over Asia, where the observed trends are also the most uncertain. Considering the total column, the models usually perform rather well for AOD, ${\mathrm{AOD}}_{\mathrm{f}}$ and the AE, but they show lower skill for ${\mathrm{AOD}}_{\mathrm{c}}$. At the ground level, the models perform well for both ${\mathrm{SO}}_{\mathrm{4}}$ concentration and PM. The trends in ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$ computed from regional time series are associated with large uncertainties due to the limited number of stations. This is exacerbated by the fact that data were only available from two models for these parameters.

4.3 Trends in models

4.3.1 Global trends

As discussed previously, the regional trends found are probably not always representative of the trends in the extended regions and over the whole study period. The reason for this is the partial spatial and temporal coverage of the ground-based observations. Moreover, the observation stations are obviously located on land. This does not allow for a depiction of global aerosol trends and is unfortunate as sea salt particles are among the most predominant aerosols on Earth .

Unlike observations, models provide data at a global scale and for the entire study period. The completeness of these model datasets offers the opportunity to derive global aerosol trends. In order to provide an assessment of the aerosol trends at a global scale, we present, in this section, the trends computed with NorESM2 (CMIP6 group), which provides data for all of the nine parameters considered in this study. The calculation of the global trend is made by averaging the absolute trends computed at each grid point of the model and using all timestamps in the study period. In order to provide a relative trend, this absolute trend is normalized to the global average of the considered parameter for the year 2000. The global trends are reported for the nine aerosol parameters in Table . The global maps, shown in Fig. , enable investigation of the spatial variability in these trends.

Table 4

Global means and trends of aerosol parameters using NorESM2 model data. The values in parentheses are obtained by aggregating only grid points where observation stations are located while using the complete model time series. The relative trends are calculated by averaging the absolute trends within the considered grid points and normalizing this average to the global mean for the year 2000.

	Mean${}_{\mathrm{2000}}$	Trend (% yr${}^{-\mathrm{1}}$)
AOD	(0.16) 0.14	($+\mathrm{0.1}$) $+\mathrm{0.2}$
${\mathrm{AOD}}_{\mathrm{f}}$	(0.09) 0.05	($+\mathrm{0.4}$) $+\mathrm{0.6}$
${\mathrm{AOD}}_{\mathrm{c}}$	(0.06) 0.09	($-\mathrm{0.2}$) $+\mathrm{0.1}$
AE	(0.78) 0.43	($+\mathrm{0.2}$) $+\mathrm{0.3}$
${\mathrm{PM}}_{\mathrm{2.5}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$)	(6.3) 9.1	($-\mathrm{1.1}$) $+\mathrm{0.2}$
${\mathrm{PM}}_{\mathrm{10}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$)	(11.7) 18.7	($-\mathrm{0.6}$) $+\mathrm{0.1}$
${\mathrm{SO}}_{\mathrm{4}}$ ($\mathrm{\mu }$g m${}^{-\mathrm{3}}$)	(1.86) 0.64	($-\mathrm{3.9}$) $+\mathrm{0.4}$
${\mathit{\sigma }}_{\mathrm{sp}}$ (${\mathrm{Mm}}^{-\mathrm{1}}$)	(13.9) 21.2	($-\mathrm{0.3}$) $+\mathrm{0.2}$
${\mathit{\sigma }}_{\mathrm{ap}}$ (${\mathrm{Mm}}^{-\mathrm{1}}$)	(2.0) 0.9	($+\mathrm{1.4}$) $+\mathrm{1.5}$

Figure 7

Global trends in aerosol properties using NorESM2 data regridded at a 5${}^{\circ }$ $\times$ 5${}^{\circ }$ resolution. The blue and red dots indicate significant negative and positive trends, respectively.

[Figure omitted. See PDF]

While the observed trends in the three AOD parameters show a decrease in most of the regions of the world, the global AOD trend is actually positive ($+\mathrm{0.2}$ % yr${}^{-\mathrm{1}}$). This global increase is also found with other models. Averages of the models from the CAMS-Rean and the AP3 groups simulate global trends of about $+\mathrm{0.2}$ % yr${}^{-\mathrm{1}}$ and $+\mathrm{0.3}$ % yr${}^{-\mathrm{1}}$, respectively. Within the CMIP6 group, IPSL and CESM2 also exhibit positive trends ($+\mathrm{0.7}$ % yr${}^{-\mathrm{1}}$ and $+\mathrm{0.3}$ % yr${}^{-\mathrm{1}}$), consistent with NorESM2, while CanESM simulates a negative trend ($-\mathrm{0.8}$ % yr${}^{-\mathrm{1}}$). The relative increase of 0.2 % yr${}^{-\mathrm{1}}$ found with NorESM2 corresponds to an absolute rate of $+\mathrm{0.0028}$ per decade, which is in excellent agreement with the global trend (over the oceans) of $+\mathrm{0.003}$ per decade reported by using MODIS data. The increase in AOD is observed to be larger for the fine fraction, with an increase of about $+\mathrm{0.6}$ % yr${}^{-\mathrm{1}}$, as compared to $+\mathrm{0.1}$ % yr${}^{-\mathrm{1}}$ for ${\mathrm{AOD}}_{\mathrm{c}}$. As seen in Fig. , similar geographical patterns are found for the three AODs: an increase in South Africa and East Asia and decrease in Europe and in the USA. The increasing AOD observed in Canada is dominated by an increase in ${\mathrm{AOD}}_{\mathrm{f}}$ in this region. The prominent increase in AOD in Indonesia seems to be linked to a large increase in ${\mathrm{AOD}}_{\mathrm{c}}$. The tropical Pacific Ocean, off the west coast of South America, has significant positive modeled trends in both AOD and ${\mathrm{AOD}}_{\mathrm{f}}$. Almost no significant trend is found south of 60${}^{\circ }$ S.

The model also simulates an increase for the AE on a global scale, with a rate of $+\mathrm{0.3}$ % yr${}^{-\mathrm{1}}$. This suggests a shift towards smaller particles. The largest increases are found over Canada, Greenland, Siberia and the Pacific Ocean. There are some distinct outliers around 60${}^{\circ }$ S. In the Atlantic, we find a decrease in the AE off the east coast of the USA, which is consistent with the decrease in ${\mathrm{AOD}}_{\mathrm{f}}$ in the same region.

The trends in both ${\mathrm{PM}}_{\mathrm{2.5}}$ and ${\mathrm{PM}}_{\mathrm{10}}$ exhibit similar geographical features to those observed for AOD. In addition, one finds large and significantly increasing trends in the high Arctic that could be explained by a change in the air mass circulation pattern or by the increase in open sea, which might contribute to a higher production of sea salt aerosols . The global averages show that ${\mathrm{PM}}_{\mathrm{2.5}}$ is increasing faster than ${\mathrm{PM}}_{\mathrm{10}}$ ($+\mathrm{0.2}$ % yr${}^{-\mathrm{1}}$ vs. $+\mathrm{0.1}$ % yr${}^{-\mathrm{1}}$), which is consistent with the increasing AE, suggesting a relatively higher fraction of fine particles with time.

The surface ${\mathrm{SO}}_{\mathrm{4}}$ concentration trends map reveals two large contrasting regions. Significant decreases are found over North America and Europe, while significant increases are found over southern and eastern Asia and southern to central parts of Africa. This illustrates the shift in polluting activities from the developed countries to the developing countries during the last 2 decades. With an overall increase of $+\mathrm{0.4}$ % yr${}^{-\mathrm{1}}$, the global trend is positive.

The ${\mathit{\sigma }}_{\mathrm{sp}}$ trends are very similar to those observed for both ${\mathrm{PM}}_{\mathrm{2.5}}$ and ${\mathrm{PM}}_{\mathrm{10}}$. The same geographical patterns are found, and a similar global average trend amounts to an increase of 0.2 % yr${}^{-\mathrm{1}}$ over the study period.

${\mathit{\sigma }}_{\mathrm{ap}}$ reveals increasing tendencies over most of the grid boxes of the model, except in Europe, the eastern part of the USA and Australia. This explains why a large positive global trend is obtained for this parameter, with an average of $+\mathrm{1.5}$ % yr${}^{-\mathrm{1}}$. Further analysis shows a good spatial correlation with the black carbon OD (optical depth), which exhibits a strong global positive trend of $+\mathrm{2.3}$ % yr${}^{-\mathrm{1}}$, as discussed below.

Table  also contains the trends computed for the different aerosol parameters when combining only the grid points where an observation station is located, whether measurements are available or not. Significant differences in so-called global trends can be found when observations are not provided over some regions. This is most obvious for ${\mathrm{SO}}_{\mathrm{4}}$, for which the observation stations are located mostly in Europe and North America and which exhibits decreasing values, while only a few stations are located in the regions associated with increasing values. In this case, the computation of the trends by considering only observation station grid boxes leads to a global decrease of $-\mathrm{3.9}$ % yr${}^{-\mathrm{1}}$, while consideration of all of the grid boxes of the model leads to a global increase of $+\mathrm{0.4}$ % yr${}^{-\mathrm{1}}$.

4.3.2 Contribution of main aerosol species to the AOD trends

The averaged global trend computed by NorESM2 indicates an increase in AOD in the 2000–2014 period with a rate of about 0.2 % yr${}^{-\mathrm{1}}$. The trends in the AE, ${\mathrm{AOD}}_{\mathrm{f}}$ and ${\mathrm{AOD}}_{\mathrm{c}}$ indicate that fine-mode particles are primarily responsible for this increase in the atmospheric column.

In this section, we investigate the trends in the major aerosol species simulated by NorESM2. For that purpose, the absolute trends in the individual contribution of these species to the AOD were computed, as well as the trends in the loads and the emissions. The trends in OD and loads are shown in Fig. . In this version, NorESM2 simulates a large proportion of sea salt. This is the result of model tuning used for reaching climate equilibrium. While the model attributes too much OD to sea salt (SS), the trends should not be affected by this tuning.

Figure 8

Absolute trends in OD and emissions of the main aerosol species computed with NorESM2. The $y$ axis of the trends in OD and the emissions is given according to the power of 10 indicated at the top left corner of each of the subplots.

[Figure omitted. See PDF]

The relative increase in AOD of $+\mathrm{0.2}$ % yr${}^{-\mathrm{1}}$ corresponds to an absolute increase in AOD of $+\mathrm{3.1}\times {\mathrm{10}}^{-\mathrm{4}}$ yr${}^{-\mathrm{1}}$. This positive trend is dominated by an increase in the species-specific ODs of the organic aerosols (OAs), ${\mathrm{SO}}_{\mathrm{4}}$ and black carbon (BC), which are responsible for an increase in the OD of about $+\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{4}}$ yr${}^{-\mathrm{1}}$, $+\mathrm{0.7}\times {\mathrm{10}}^{-\mathrm{4}}$ yr${}^{-\mathrm{1}}$ and $+\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{4}}$ yr${}^{-\mathrm{1}}$, respectively. The relative OD trends give a different ranking since the highest increase is found for BC ($+\mathrm{2.5}$ % yr${}^{-\mathrm{1}}$), followed by OAs ($+\mathrm{0.5}$ % yr${}^{-\mathrm{1}}$). On average, the trends for dust and sea salt OD are slightly negative ($-\mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{4}}$ yr${}^{-\mathrm{1}}$). Note – these species trends include any associated water which can change as a function of relative humidity.

The trends in OD do not necessarily represent the trends in the aerosol loads, which do not include associated water. The different species have different global mass extinction coefficients (calculated in this study as OD per load; dust – 1.8 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$, SS – 4.3 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$, OA – 5.6 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$, ${\mathrm{SO}}_{\mathrm{4}}$ – 5.3 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$, BC – 7.6 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{g}}^{-\mathrm{1}}$). For sea salt, opposite trends are observed for the sea salt OD (positive trend) and the sea salt load (negative trend). The analysis of the global maps (not shown in this study) reveals that the largest increases in the sea salt load happen in Indonesia and near the North Pole and result in a relatively larger increase in OD in these areas. These localized increases in sea salt OD drive the global sea salt OD trend and are due, at least in part, to the higher relative humidity at these latitudes which makes the sea salt, which is very hygroscopic, more efficient at light extinction.

5 Conclusions

The main findings of this multiparameter trend analysis are listed below:

• The observations exhibit mostly negative trends regarding the extensive parameters in the different regions of the world. Significant decreases are found in Europe, North America, South America, North Africa and Asia. In Asia, the AE increases in time and is consistent with increases in ${\mathrm{AOD}}_{\mathrm{f}}$ and ${\mathrm{SO}}_{\mathrm{4}}$, which reflects the regional increase in the anthropogenic aerosols in that region in the overall study period from 2000 to 2014.

• Some observation networks allow for the derivation of representative trends over the whole study period. In other cases, the limited temporal and spatial coverage of the observations can induce artificial and/or highly uncertain trends when using regional time series. Among the 38 computed trends with observation data, 22 are considered as representative of the actual trends occurring in the whole region and study period.

• The models tend to capture observed AOD, AE, ${\mathrm{SO}}_{\mathrm{4}}$ and PM trends but show larger discrepancies regarding ${\mathrm{AOD}}_{\mathrm{c}}$. The smaller number of data available for establishing ${\mathit{\sigma }}_{\mathrm{sp}}$ and ${\mathit{\sigma }}_{\mathrm{ap}}$ trends makes the validation of the modeled trends more uncertain.

• The rather good agreement of the trends across different aerosol parameters between models and observations, when co-locating them in time and space, implies that global model trends, including those in poorly monitored regions, are likely correct.

• The global trends computed with NorESM2 (CMIP6 group) model data give a different picture than the trends obtained when using only ground-based observations. Global positive trends are found for all of the parameters related to aerosol loading. The trends in AOD are dominated by the increase in the fine particles both in the column and at the surface. This tendency towards finer particles is consistent with the positive trend in the AE. This increase appears to be dominated by organic aerosols, of which the emissions have increased in the study period, and by ${\mathrm{SO}}_{\mathrm{4}}$ aerosols, whose sources were shifted from Europe and North America to Africa and East Asia and for which a global positive ${\mathrm{SO}}_{\mathrm{4}}$ trend is found.

Some elements were not considered in this study, and they could be investigated in order to complete the aerosol trends picture:

• Some regions are associated with strong seasonal cycles. In South America, the regional time series shows high peaks in AOD, associated with forest fires in the late summer, whose intensity greatly varies from year to year. In Africa, a strong seasonal contrast is also found due to the transport of desert dust at altitude in the summer months . The computation of the seasonal trends would allow characterization of the tendencies in such extreme or synoptic aerosol events.

• This study shows that the trends computed from the ground-based observation networks are not representative of the global aerosol trends due to the inhomogeneities in data spatial coverage. The satellites providing a global Earth observation could be utilized for the evaluation of the model trends in the regions lacking observations and over the oceans .

• The trends in the meteorological parameters could be investigated in parallel with the aerosol trends because they affect the aerosol life cycle and their optical properties . Hypothetical trends in wind velocity could produce trends in the loads of sea salt and dust and, as seen in the last section, trends in OD could also be enhanced by relative humidity changes. Changes in temperature could impact the magnitude of the biogenic emissions. Indeed, increasing temperatures, associated with changes in land use and high atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentrations, have been shown to lead to an increase in the VOC emissions . Finally, trends in precipitation that are responsible for aerosol wet scavenging would directly impact trends in aerosol loads.

• Several studies have linked the trends in anthropogenic aerosols to radiative-forcing variations while investigating sources of global dimming and brightening . It could be of interest to evaluate how much the modeled trend deviations, as compared to the observations, are affecting the calculation of the radiative forcing, in the different regions of the world and at a global scale.

• While the mountain sites were excluded from this study, it could be of interest to investigate the trends at higher altitude (which may be related to changes in long-range transport) by including the in situ and remote sensing stations higher than 1000 m (Jungfraujoch, Mauna Loa Observatory, etc.). Similarly, it may also be of interest to look at trends in smaller regions (e.g., split North America into several subregions which are more internally consistent in terms of climate and environment than the large North America region defined here or consider southern Europe as its own region rather than combining it with the North Africa region as was done here).

Code availability

The observation and model data were read and co-located with the pyaerocom Python library, developed by MET Norway (https://doi.org/10.5281/zenodo.4139879, version 0.9.1) .

The supplement related to this article is available online at: https://doi.org/10.5194/acp-20-13355-2020-supplement.

Author contributions

AM coordinated the study, was responsible for the statistical calculation and analysis, and wrote the paper. JG is the main developer of the pyaerocom library. MS provided feedback on the methods and the manuscript. WA, EA, JH and PL provided in situ data, contributed to the writing of the observation dataset section and provided feedback on the manuscript. BH is the principal investigator of AERONET. HB, MC, PG, HZ, ZK, AK, TL, GM, DN, DO, KvS, TT, RBS and ST provided model output data and feedback on the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “The Aerosol Chemistry Model Intercomparison Project (AerChemMIP)”. It is not associated with a conference.

Acknowledgements

Data providers from all the regional and global networks are greatly acknowledged for sharing and submitting their data to be used. The ECHAM-HAMMOZ model is developed by a consortium composed of ETH Zurich, Max Planck Institut für Meteorologie, Forschungszentrum Jülich, University of Oxford, the Finnish Meteorological Institute and the Leibniz Institute for Tropospheric Research and is managed by the Center for Climate Systems Modeling (C2SM) at ETH Zurich. Computing and data storage resources, including the Cheyenne supercomputer (10.5065/D6RX99HX), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. All simulations were carried out on the Cheyenne high-performance computing platform https://www2.cisl.ucar.edu/user-support/acknowledging-ncarcisl (last access: 20 November 2020), and are available to the community via the Earth System Grid. High-performance computing and storage resources were provided by MET and the Norwegian infrastructure for computational science (through projects NN2345K, NN9252K, NS2345K and NS9252K).

Financial support

This research has been supported by the European Commission, H2020 Research Infrastructures (FORCeS (grant no. 821205), CRESCENDO (Coordinated Research in Earth Systems and Climate: Experiments, Knowledge, Dissemination and Outreach; grant no. 641816)), the National Science Foundation (NSF; grant no. 1852977), the Research Council of Norway (grant no. 295046) and the Copernicus Atmospheric Monitoring Service (grant no. CAMS84).

Review statement

This paper was edited by Pedro Jimenez-Guerrero and reviewed by two anonymous referees.