On the seasonal variation in observed size distributions in northern Europe and their changes with decreasing anthropogenic emissions in Europe: climatology and trend analysis based on 17 years of data from Aspvreten, Sweden

Size-resolved aerosol trends were investigated based on a 17-year data set (2000–2017) from the rural background site Aspvreten located in southern Sweden (58.8^∘ N, 17.4^∘ E). Cluster analysis of the size distributions was performed to aid in the interpretation of the data. The results confirm previous findings of decreasing aerosol mass and number during the last decades as a result of reduced anthropogenic emissions in Europe. We show that both particle modal number concentration and size have substantially been reduced during the last 17 years. Negative trends in particle number concentration of about 10 cm⁻³ yr⁻¹ are present for nuclei, Aitken, and accumulation modes. In total, integral particle number concentration has decreased by 30 %, from 1860 to ca. 1300 cm⁻³. The reduction in modal number concentration is accompanied by a decrease in modal size, and this decrease is largest for the accumulation mode (2 nm yr⁻¹ or about 17 % for the whole period). These reductions have resulted in a decrease in submicron particle mass (< 390 nm) by more than 50 % over the period 2000–2017. These decreases are similar to observations found at other stations in northern Europe.

Although all size classes show a downward trend as annual averages, we also show that observed trends are not evenly distributed over the year and that a rather complex picture emerges where both sign and magnitude of trends vary with season and size. The strongest negative trends are present during spring (accumulation mode) and autumn (Aitken mode). The strongest positive trends are present during summer months (Aitken mode). The combined trajectory and data analyses do not present evidence for an increase in new particle formation formed locally, although some evidence of increased new particle formation some distance away from the receptor is present. Observed aerosol size distribution data, together with an adiabatic cloud parcel model, were further used to estimate the change in cloud droplet concentration for various assumptions of updraught velocities and aerosol chemical composition. The results indicate a substantial increase in the atmospheric brightening effect due to a reduction in cloud reflectivity corresponding to 10 %–12 % reduction in cloud albedo over the period 2000–2017.

1 Introduction

Turbidity of the atmosphere is a result of light being scattered and absorbed by particles suspended in the air. For most ambient conditions there is a strong proportionality between the light scattered and the total particle volume in the accumulation mode range between about 0.1 and 1 µm diameter (Willeke and Brockmann, 1977). This size dependence arrives from the fact that the wavelength of the peak in solar radiation is in the same size range as the accumulation mode particles. Hence, variations in the number of accumulation mode particles directly affect the turbidity of the atmosphere.

An increase in aerosol loading tends to scatter more of the incoming solar radiation back to space, the so-called “dimming effect”, which leads to less radiation reaching the Earth's surface and a net cooling of the climate system (Myhre et al., 2013). Observations of solar radiation incident at Earth's surface indicate that global dimming increased up to about 1990, but after that it started to decrease, the so-called “brightening effect” (Wild et al., 2006). Streets et al. (2006) suggested that this trend in global dimming–brightening was due to the combined effect of economic growth and the recent decrease in emission of aerosol particles and their gaseous precursors as a result of legislative measures.

The most important aerosol precursor is sulfur dioxide, which has both natural and anthropogenic sources. Sulfur dioxide is oxidized to sulfuric acid vapour and later transforms into sulfate aerosol mass, fully or partially neutralized by ammonia. Sulfuric acid vapour both participates in new particle formation and contributes to aerosol mass concentration by condensation on already existing particles. The above-mentioned air quality policies mainly targeted sulfur emissions, which have led to a general decrease in the anthropogenic emissions of SO₂ since the 1970s (Smith et al., 2011; Klimont et al., 2013). With perhaps the exception of the Indian subcontinent, this decrease appears to also continue through the last decade (Li et al., 2017; United States Environmental Protection Agency, 2018).

Zhao et al. (2017) showed that the decrease in SO₂ and other important anthropogenic emissions, such as nitrogen oxides and carbonaceous material, is consistent with the decrease in the trends in dimming observed from space-borne instruments. They studied three regions, the eastern United States, western Europe, and eastern and central China, and concluded that the trends in aerosol optical depth (AOD) are consistent with the trends in emissions for these three regions.

Collaud Coen et al. (2013) studied trends in observed aerosol optical properties measured at ground-based stations mainly located in North America and Europe. They concluded that, even if the trends are not homogeneously distributed geographically, the decreasing trend observed for most stations in North America is related to the decrease in anthropogenic emissions of particles and their gaseous precursors. The pattern of trends over Europe was not as clear as over North America, and the Arctic and Antarctic stations did not show significant trends at all.

On the other hand, Turnock et al. (2015) found that the reduced emissions have led to a substantial decrease in the concentration of aerosols over Europe. They used a global environmental model to consolidate the observed trends in aerosol loading over Europe and changes in emissions. They also explored the trends in the microphysical properties of the aerosol by comparing the model to observed number densities of different size ranges. They concluded that the model under-predicts the observations in general, and in particularly this was the case for the two Arctic stations included in the study. The authors did not look at the trends for different particle size ranges per se but focused on how the model reproduced the data at the various sites.

By studying long-term observations of aerosol microphysical data available through the Global Atmospheric Watch (GAW) programme, organized by the World Metrological Organization (WMO) and the European Union infrastructure project ACTRIS, Asmi et al. (2013) explicitly studied the decadal trends in total aerosol number concentrations. They came to a similar conclusion as the studies above, that the most likely cause for the generally decreasing trends in aerosol number concentrations in the Northern Hemisphere is the decrease in anthropogenic emissions of SO₂ or other co-emitted species. To address the question about possible trend in cloud condensation nuclei (CCN) particles, the authors divided the observed size distribution into two number densities of particles larger than 20 nm (N20) and 100 nm (N100). Here, N100 was a proxy for CCN particles. However, due to claimed lack of stations with more than 10 years of data, only five stations were available for their study. From this analysis they made two interesting observations. The first observation was that the trends in both N20 and N100 were consistent at the four sites (three stations in Finland and one in Sweden), which led the authors to conclude that the trends are uniform over the size distribution. The second observation was the lack of a diurnal difference in trends. This led the authors to hypothesize that either the characteristics of the size distributions are mainly controlled by long-range transport or the frequency or intensity of new particle formation has stayed rather constant over the investigated time period.

The statistics of the occurrence of nucleation events over an 8-year period were explicitly studied by Dal Maso et al. (2005). This study was conducted for one of the Finnish stations participating in the study by Asmi et al. (2013), namely the Hyytiälä SMEAR II station. The authors do not explore any possible trends in the occurrence of new particle formation, but their Fig. 6 suggests that between 1996 and 2003 the number of event days (days with clear new particle formation) and undefined days (cases with small particles but not a clear event day) increases over time, while clear non-event days suggest little change. How this potentially small increase in events is distributed over the year is not possible to assess from their study.

Clearly there is a lack of long-term aerosol microphysical data from which trend analysis has been reported or even can be reported because of too short measurement periods or interrupted data sets. Trend analysis requires long-term uninterrupted time series to be meaningful. This does not mean that data are completely missing. Through collaborations within European Supersites for Atmospheric Aerosol research (EUSAAR, http://www.eusaar.net/, last access: 5 April 2019), long-term standardized observations of aerosol number size distribution observations are available from several European sites, some of which are long enough to support time series analysis and determination of trends (e.g. Asmi et al., 2013). The characteristics of the general properties of aerosols are better known as a large number of studies address the different aspects of the aerosol number size distribution observed in Fennoscandia. In this region, several long-term monitoring stations exist, covering vegetation zones from nemoral to polar sites north of the treeline, e.g. the southern Sweden nemoral site Vavihill (56.01^∘ N, 13.09^∘, 172 m a.s.l.; Tunved et al., 2003; Kristensson et al., 2008), the southern boreal site Hyytiälä (61.85^∘ N, 24.29^∘, 179 m a.s.l.; Dal Maso et al., 2005; Kulmala et al., 2008), and the northerly located stations Värriö (67.76^∘ N, 29.61^∘, 390 m a.s.l.) and Pallas (67.97^∘ N, 24.12^∘, 560 m a.s.l.; Laakso et al., 2003; Tunved et al., 2003; Dal Maso et al., 2007). Available reports show that both the shape and number concentration over Scandinavia are highly variable. At northerly located stations such as Värriö and Pallas (∼67^∘ N), reported annual average number concentration of nuclei mode particles is around 50–100 cm⁻³, Aitken mode particles ca. 200 cm⁻³ and accumulation mode particle concentration ∼150 cm⁻³ (Tunved et al., 2003; Komppula et al., 2006; Dal Maso et al., 2007). Observations at Hyytiälä (61^∘51^′ N) yield averages of ca. 300 cm⁻³ nucleation mode size particles and 750 and 290 cm⁻³ in the Aitken and accumulation mode size ranges, respectively. In southern Sweden, Kristensson et al. (2008) report accumulation mode concentration around 300 cm⁻³ and Aitken mode concentration around 800 cm⁻³. At all stations above, there is a strong seasonality in aerosol size distribution properties, with a typically bimodal wintertime size distribution with, on average, low number concentration. This is contrasted by a near unimodal summertime size distribution with high total number concentration. It has further been demonstrated in studies of air mass transport across Scandinavia that the Boreal region acts as a net source of particle number (and mass) during southerly transport of marine clean air (Tunved et al., 2006a, b; Väänänen et al., 2013), while acting as a relative sink when polluted air is transported northwards. During the southerly transport, it is believed that the forest is playing an active role in particle formation and growth, via nucleation and condensation. During northerly transport, transport of air masses with a comparable high aerosol condensation sink leads to rarer nucleation events and sink processes; in particular wet deposition leads to reduction in number and mass.

How these mean aerosol characteristics change over time may contain surprises, as pointed out by Julin et al. (2018). Using a chemical transport model, the authors showed that aerosol properties in the Nordic countries may have a very different temporal evolution compared to central Europe with respect to the total number density of particles. The model simulations suggest that changes in aerosol properties may not only follow the reduction in emissions, but the number density of particles may even increase as a result of legislative actions. The explanation for this is an already relatively low emission of precursor gases together with a general decrease in the condensational and coagulation sink by pre-existing particles in neighbouring countries, which sets the stage for new particle formation (Kulmala et al., 2001; Kerminen and Kulmala, 2002).

Based on the literature referenced above, several locations exhibit a decrease in aerosol mass observed in overpopulated areas in the last decades and this decrease is related to the decrease in emissions of anthropogenic primary particles and precursor gases. There is also a suggestion that this decrease may not be uniformly distributed over the whole size distribution. In our opinion it is important to establish how this reduction is distributed in order to understand the implications of the reduction in emissions, which means that the whole size distribution must be studied. There might be a general decrease in aerosol loading that is shared by all particle sizes, but if one type of aerosol size distribution becomes less frequent, this must be compensated for by some other type of size distribution (e.g. size distributions dominated by accumulation mode, Aitken mode or nucleation mode).

We have chosen to study hourly averaged size distributions observed at Aspvreten, Sweden. The data extend over a 17-year period from June 2000 through November 2017. The measurement site at Aspvreten is characterized as a rural background station. Our aim is to study the aerosol size-resolved trends and try to establish how the characteristics of the aerosol size distributions have changed over time as a potential result of reduction in anthropogenic emissions. We will explore seasonal differences and possible changes in transport patterns. To test how potential changes in aerosol properties at Aspvreten may influence cloud droplet number concentrations, a simple parcel model was used to simulate cloud droplet activation. The model was initiated using observed size distributions for a range of updraught velocities and characteristic chemical aerosol properties.

2 Methods

2.1 Aerosol number size distribution

Aspvreten observation station (58.8^∘ N, 17.4^∘ E, 25 m a.s.l.) is located in the county of Södermanland, about 80 km south of Stockholm. The station is situated close to the Baltic Sea, and the coastline is located a few kilometres to the east and about 1 km to the south of the station. The surroundings are dominated by deciduous and conifer forests, combined with farmlands. The station represents typical continental rural background conditions with few local pollution sources or densely trafficked roads. Climatologically, the surroundings represent the boreal–nemoral zone, a transition region between the temperate southern nemoral zone (mainly deciduous broadleaf forest) and the boreal zone.

This study focuses on aerosol number size distribution observations collected at the Aspvreten observation station using a differential mobility particle sizer (DMPS). The station has been housing aerosol number size distributions since the year 2000 and Stockholm University has performed observations of aerosol number size distribution properties with nearly identical instrumental setup for almost 2 decades. In 2018 the measurements were moved to Norunda, which is about 90 km north of Stockholm, and these data are not included in current analyses. Nevertheless, the DMPS data set from Aspvreten is one of the longest in Europe. The instrument setup consists of a medium differential mobility analyser (DMA) with a TSI 3010 condensation particle counter (CPC), and the system is capable of monitoring the aerosol number size distribution between 10 and roughly 400 nm with a temporal resolution of around 15 min. The actual length of each scan varied slightly over the years, but our base data will be hourly averages, which makes this a minor issue.

More details on the instrumental setup and station characteristics may be found in Tunved et al. (2003, 2004, 2005) The station, being a part of the European Supersites for Atmospheric Aerosol Research project (EUSAAR, http://www.eusaar.net/, last access: 5 April 2019), has undergone audits to assure the application of harmonized measurement and sampling protocols regarding aerosol number size distribution.

In this study, we utilize all the data recorded at Aspvreten, from 2000 to 2017. After screening for instrument failures and maintenance periods, the data coverage is typically better than 75 % on annual basis. The only years when less than 75 % data coverage occurred was for the starting year 2000 (the year of the measurement start-up, with 25 % data coverage), for 2002 (70 % data coverage) and for 2010 (50 % data coverage). In this study, we have made use of hourly average size distributions, and a total of 132 182 hourly data points have been calculated, equivalent to more than 5500 d of data, or 83 % coverage for the whole period. Since the actual bin sizes of the instrument have changed slightly over the years, all data have been harmonized to a fixed logarithmic equidistant size grid distributed over 34 bins between 10 and 390 nm. Previous studies using aerosol number size distribution observations from Aspvreten include Tunved et al. (2003, 2004, 2005) and Dal Maso et al. (2007, 2008). These studies have mainly focused on new particle formation, seasonality, life cycle and transport characteristics of the Fennoscandian aerosol. The current study is the first decadal trend analysis based on aerosol size distribution data from Aspvreten.

The availability of a long time series (17 years) that has used nearly identical instrument setup during the full operational time makes this data set most suitable for the study of long-term trends of aerosol size distribution properties, on par with other long-term archives such as the approximately 20-year aerosol size distribution data set from Hyytiälä (e.g. Asmi et al., 2013).

During the entire period of interest, the instrument has been maintained by the same technician. This by itself assures consistent calibration routines and standard operating procedures. Further, the DMPS system has been part of the European Supersites for Atmospheric Aerosol Research (EUSAAR, http://www.eusaar.net/, last access: 5 April 2019) intercalibration workshop, and the setup follows the recommendations made by EUSAAR regarding sampling and inversion. This certifies that any error introduced in the trend analysis by changes in sampling routine quality assurance is minimized.

2.2 Log-normal fitting

All hourly aerosol size distributions were fitted assuming that the aerosol can be described by the sum of log-normally distributed modes. We decided a priori to fit the distribution over three modes after concluding that three modes in most cases were sufficient to capture the size distribution properties over the available size range. Equation (1) gives the general form of a log-normal mode,

$$\begin{array}{}\text{(1)}& n_N\left(\log D_{\mathrm{p}}\right)={\displaystyle \frac{N}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\mathrm{1}/\mathrm{2}}}}\exp \left(-{\displaystyle \frac{{\left(\log D_{\mathrm{p}}-\log {\stackrel{\mathrm{‾}}{D}}_{\mathrm{pg}}\right)}^{\mathrm{2}}}{\mathrm{2}{\log }^{\mathrm{2}}{\mathit{\sigma }}_g}}\right),\end{array}$$

where n_N is the number size distribution function, N is the number concentration of the mode, D_p is the size of the particle diameter, ${\stackrel{\mathrm{‾}}{D}}_{\mathrm{pg}}$ is the geometrical mean diameter of the mode and σ_g is the geometric standard deviation, henceforth denoted GSD, giving an estimate of the width of each log-normal mode centred around ${\stackrel{\mathrm{‾}}{D}}_{\mathrm{pg}}$.

The fitting was performed utilizing the MATLAB function fmincon.m to perform a constrained fit of the three modes over the size range covered by the instrument. The number concentration of each mode is allowed to vary from 0 to infinity, and the GSD is allowed to vary between 1.1 and 2.5. When fitting, we do not force the modes into predefined size ranges, i.e. typical nuclei, Aitken or accumulation mode size ranges, but instead allow the algorithm to find the best numerical solution over the whole size range covered by the instrument. After fitting has been performed, the modal parameters are arranged according to size (${\stackrel{\mathrm{‾}}{D}}_{\mathrm{pg}}$), resulting in modes 1–3, with mode 1 being the smallest and mode 3 being the largest. The modal sizes derived in this way roughly correspond to nuclei, Aitken and accumulation modes.

2.3 Cluster analysis

Clustering of aerosol size distribution has proven to be a very useful tool in studies of aerosol life cycle and different process studies, as demonstrated in several studies following Tunved et al. (2004), including Beddows et al. (2009), Wegner et al. (2012), Freud et al. (2017) and Varanda Rizzo et al. (2018) in rural, urban and Arctic environments.

Contrary to standard averaging of number size distributions, cluster analysis and associated centroids can conserve the shape of the aerosol size distribution. Hence, size distribution clusters represent “signature distributions” that reflect contribution from members that are likely to have undergone similar processing in the atmosphere prior to observations. Thus, clustering size distribution and combining the cluster analysis with auxiliary parameters, such as trajectory-derived source areas, temporal distribution of members and parameters related to sink processes (e.g. precipitation), can provide a deeper insight into the multitude of factors defining the aerosol over time.

In this study we have applied a MATLAB version (Statistics and Machine Learning Toolbox) of k-means clustering (kmeans.m) to perform a clustering of the hourly size distributions. Naturally, there are several ways to determine the optimum number of clusters to use in order to maximize the inter-cluster variability while at the same time minimizing the intra-cluster variability. Various functions are available in MATLAB for this, e.g. the silhouette.m function. Although similar cluster evaluation algorithms rely on sound mathematical foundations, in the end a subjective decision must be made in order to decide the criteria used, and this criterion should reflect the needs with respect to level of detail of the project at hand. In Tunved et al. (2004), eight clusters were used to characterize 1 year of data from Aspvreten. For the trend analysis in this study, we increased the number of clusters to 12, in order to better capture changes between clusters over time of the much larger data set. Tests were performed with both more and fewer clusters. The best balance between information content and data amount to be presented was found to be reasonable around 12 clusters. A balance of information content refers to being able to follow a logical context between different clusters, but also that the whole set of clusters can characterize the domain of different size distributions observed at Aspvreten.

The clustering was performed on hourly averaged data, using “max iterations” of 10 000 and “number of replicates” set to 10 in MATLAB. The distance function applied was squared Euclidean distance, assuming that the difference is calculated from the centroids defined as the mean of the points in the clusters,

$$\begin{array}{}\text{(2)}& d\left(x,c\right)=\left(x-c\right){\left(x-c\right)}^{\prime },\end{array}$$

where x is an observation (i.e. the size distribution vector) and c is the centroid.

2.4 Trajectory analysis

In order to explore source dependence of aerosol properties and to allow the study of the effect different meteorological parameters have on the observed size distribution characteristics, air mass back trajectories were calculated using the HYSPLIT4 model (Draxier and Hess, 1998). Trajectories were calculated on an hourly basis, with a total trajectory length of 240 h. Arriving altitude was set to 100 m a.g.l., and model-derived relative humidity, precipitation amount, mixing layer height and temperature were saved for each trajectory point. The meteorological input used by the model is the FNL archives for the time period 2000–2005 and the GDAS1deg for the time period 2006–2017. The main difference between these two data sets lies in the resolution; GDAS1deg is given with a 1 ^∘ resolution, while the FNL data set uses a resolution of about 191 km. The FNL archive is no longer updated and was replaced by the 1 ^∘ GDAS data set in 2005. Another difference is that the FNL data are provided as hemispheric data, two files per month, while the GDAS1deg data are provided as weekly files with global coverage (for details cf. https://www.ready.noaa.gov/archives.php, last access: 1 December 2019).

3 Results

3.1 General seasonal size distribution properties

As discussed above, the aerosol over the Fennoscandian region is highly variable, and observations show a strong seasonality in aerosol size distribution properties. Figure 1 presents the average seasonal variation in all validated number size distributions for the period 2000–2017. Integral number is superimposed on the surface plot as median and 25th–75th percentile range. As can be seen, a strong seasonality is present. The lowest number concentration is found in November and December whilst the highest number concentration is found during summer. It is also evident that the changes of aerosol number size distribution occur gradually and the transition between the seasons is smooth, quite opposite to the seasonal characteristics of, for instance, the Arctic site Zeppelin (Tunved et al., 2013). The seasonally averaged hourly number size distributions are presented in Fig. 2 as a complement to Fig. 1. Figure 2 shows seasonal median number size distribution for the whole period 2000–2017 including the 25th–75th percentile range for each size distribution bin. Four seasons are considered: spring (March–May), summer (June–August), autumn (September–November) and winter (December–February). During spring the median integral concentration N_tot is around 1300 cm⁻³. The median size distribution has a bimodal shape and is dominated by an Aitken mode centred on 50 nm and an accumulation mode around 150–200 nm. During the summer period, the median size distribution is rather mono-modal and has a less pronounced accumulation mode and a larger Aitken mode around 70 nm. The integral concentration is around 2000 cm⁻³. The autumn period is fairly similar compared to the spring period, although the number concentration is slightly higher (∼1900 cm⁻³). The median winter size distribution is distinctly bimodal, with one dominating Aitken mode around 40–50 nm and an accumulation mode around 150–200 nm. Median integral number concentration is the lowest of all seasons, around 900 cm⁻³.

Figure 1 Seasonal variation in the aerosol number size distribution between 10 and 390 nm presented as daily median aerosol number size distribution for the whole study period, 2000–2017. Superimposed on the surface plot are the median and quartile ranges of integral number concentration.

Figure 2 Seasonally average size distributions observed at Aspvreten 2000–2017. The shaded blue area gives the 25th–75th percentile range, and the dashed line shows the median size distribution. Shown in the sub-frames are the same data but now in log-log scale. Spring: March–May; summer: June–August; autumn: September–November; winter: December–February.

It should be noted that seasonal averages of daily mean aerosol number size distributions do not preserve the signature of new particle formation events (for details regarding new particle formation, see Kulmala et al., 2004). The lack of a distinct nuclei mode in Fig. 2 does however not imply that nucleation is absent in the data set, but rather it suggests that the intermittent behaviour and short lifetime of the nuclei mode under conditions characteristic for the Aspvreten station lead to a masking of these features when performing long-term averaging. Thus, in Supplement Fig. S1, we show time-of-day seasonal mean size distributions. As evident, the signature of new particle formation events is present for all seasons except wintertime. The 10 nm particles are typically observed around noon but grow rapidly into larger size classes after a couple of hours.

Table 1 Statistics of modal fits per season. Table shows statistics derived from fitted hourly number size distributions. Indicated in the table are Nuclei, Aitken and accumulation mode parameters as the median and 25th–75th percentiles. GSD represents the geometric standard deviation and D_g the geometric mean diameter of each of the log-normal modes.

Seasonal averages of the three log-normal modes (referred to as nuclei, Aitken and accumulation mode) fitted to the hourly size distributions are presented in Table 1. The accumulation mode is largest during summer and smallest during winter, while autumn and spring have approximately the same number concentration and modal diameter as the accumulation mode (around 270 cm⁻³ and 160 nm, respectively). The winter period has the largest accumulation mode size (173 nm) and the lowest number concentration (∼200 cm⁻³). Aitken mode is largest during summer, regarding both number concentration and size (800 cm⁻³ and 76 nm, respectively). Spring and autumn have quite similar parameters for Aitken mode (average modal number and size of around 600 cm⁻³ and 60–70 nm, respectively), while winter has the lowest concentration. The nuclei mode is smallest in size, but at the same time largest in number during spring and autumn. This suggests more recent formation or slower growth of newly formed particles compared to the summer period when the size of the nuclei mode is slightly larger. Wintertime data have the lowest concentration of nuclei mode particles.

In general, the winter size distribution shares features with aged aerosol, and the low number concentration suggests comparably efficient removal. The distinct minima between Aitken and accumulation mode are indicative of cloud processing (Hoppel and Frick, 1990). Absence of small particles suggests a reduced amount of nucleating species likely due to reduction in photochemical activity and/or seasonal changes in primary sources. Sulfuric acid vapour has a short residence time and is rapidly consumed by either gas-particle formation or deposition.

Summertime distributions suggest ageing by condensation growth of both nucleated particles and primary emissions, and with a lesser influence from cloud processing and wet removal. This is in line with more photochemistry, less frequent precipitation and larger emissions from vegetation during the summer months, as pointed out earlier in Tunved et al. (2004). The fact that both spring and autumn have smaller modal size of the Aitken mode suggests less efficient growth by condensation, which may result from less effective photochemistry and lesser source strength of biogenic volatile organic compounds such as monoterpenes, sesquiterpenes and isoprene from vegetation. The larger fraction of smaller nuclei mode particles during these seasons is well in line with studies of nucleation events over the boreal region that suggest a maximum in nucleation occurring during spring and autumn (Dal Maso et al., 2005). On average, size distribution properties from the 17-year data record fall in between those values presented for Hyytiälä and Vavihill (e.g. Tunved et al., 2003, 2005; Kristensson et al., 2008).

3.2 Cluster analysis of size distributions

This section describes the results from the cluster analysis. As stated under Sect. 3.2, the clustering was performed on hourly means, comprising roughly 130 000 size distributions. This approach captures signature size distributions in different stages of the aerosol life cycle, including aerosol number size distribution types that originate from new particle formation events (e.g. Kulmala et al., 2004). For each of the 12 clusters (see Sect. 2.3) a number-weighted mean diameter was calculated, and the clusters were organized in order of increasing diameter except for what became the starting cluster 1 (this cluster stands out by having significantly lower number concentration compared to other clusters). This organization was performed to better follow the ageing of an aerosol population, where cluster 1 represents a cyclic point for the aerosol evolution (both a starting and end point). A similar approach was used by Tunved et al. (2004) in order to place the different clusters in the context of the aerosol life cycle. The median and quartile distributions for each cluster are presented in Fig. 3 and the modal parameters are provided in Table 2.

Table 2 Fitted modal parameters per cluster shown in Fig. 2. GSD represents the geometric standard deviation and D_g the geometric mean diameter of each of the log-normal modes.

Figure 3 Aerosol number size distribution clusters for the 2000–2017 data record from Aspvreten, 58.8^∘ N, 17.4^∘ E. (a) Total of 132 182 hourly size distributions included in the clustering. Number of cluster members is indicated in each frame. Clusters are ordered from left to right and up to down. The red shaded area indicates the 25th–75th percentile range of the cluster members and (b) the dashed lines indicate median number size distribution of that cluster. Indicated in blue shading is the 25th–75th percentile. GSD represents the geometric standard deviation and D_g the geometric mean diameter of each of the log-normal modes.

The seasonal and diurnal frequencies of occurrence of each cluster are presented in Figs. 4 and 5, respectively. Cluster 1 represents the most common aerosol size distribution type, a clearly bimodal size distribution with a dominating Aitken mode. This size distribution is observed about 26 % of the time, and both shape and number concentration suggest extensive processing by clouds. This cluster is observed mainly during the winter period as can be seen in Fig. 4. The average number concentration is about 700 cm⁻³, which is a little more than one-third of the annual average concentration (around 1600 cm⁻³). It is interesting to note that observations of size distributions belonging to cluster 1 exhibit a quite pronounced diurnal pattern (Fig. 5), with a maximum around noon. As this cluster is not related to new particle formation (clear absence of freshly formed particles), we hypothesize that cluster 1 has contributions from air entrained from the free troposphere, resulting in the midday enhancement in occurrence.

While cluster 1 in some respects represents the terminal size distribution, aged under the influence of clouds and precipitation, it also sets the stage for new particles to form due to the low associated condensation sink. In this cycle, cluster 2 represents the youngest form of aerosol size distribution observed. The size distributions belonging to this cluster are evidently shifted towards the smaller size classes with an open-ended size distribution towards particles smaller than instrumental limit of 10 nm. Cluster 2 has the highest average integral number concentration, around 7500 cm⁻³. The number of larger sized particles is however very small, indicating transport from cleaner regions or alternatively recent removal by wet deposition. At the same time, this is also the least frequent cluster. The mean diameter of the smallest mode is 15 nm, which also controls the total aerosol number. From Figs. 4 and 5 it is also clear that this is mainly a spring and autumn phenomenon occurring around noon. As the lifetime of such small particles is typically short, this cluster must represent recently formed aerosols. This is consistent with the timing of occurrence, i.e. the time of most efficient photochemical production of aerosol precursors.

Figure 4 Seasonal distribution of clusters 1–12. Panels indicate the number of cluster members observed for a given month (UTC).

Figure 5 Diurnal variability of clusters 1–12. Panels indicate the number of cluster members observed at a given time of day (UTC).

Clusters 3, 4, 5 and 6 all share a similar seasonal pattern, with maxima during spring and autumn. The different clusters have varying number concentration, but as pointed out above they are organized according to increasingly larger modal diameter. We understand these clusters as representing different stages of growth following nucleation (i.e. cluster 2). This view is supported by the diurnal distribution of the four clusters; cluster 2 peaks around noon, cluster 3 peaks slightly later around 15:00 UTC and cluster 4 is shifted yet later into the afternoon. Cluster 5 is mainly observed around 22:00 UTC, and cluster 6 is interpreted as growth of the aerosol extending into the following day. New particle formation and subsequent growth observed at boreal sites in northern Europe can typically be traced back in time for about 12 to 24 h. The peak frequency of occurrence and the increase in modal sizes of clusters 2–6 are consistent with the evolution observed for individual events. The cluster 2 nuclei mode appears at 15 nm around noon, and cluster 6 is observed 12 h later with a shift in nuclei modal diameter of ∼40 nm. This change in size would correspond to an average growth rate of the nuclei mode of about 2 nm per hour, which is in the lower range of reported values from Hyytiälä (see Kulmala et al., 2004, and references therein).

Clusters 7, 8 and 9 show less distinct features and are more difficult to place in the evolution timeline, although the most frequent occurrence of cluster 7 between midnight and 10:00 UTC would still fit growth that extends into the following day. The three clusters do however present larger accumulation mode sizes than clusters 2–6. Moreover, clusters 7 and 9 are mostly encountered during summer, while cluster 8 is observed mainly during spring and autumn. It is also noted that cluster 7 and cluster 9 have higher-than-average number concentration, with cluster 9 having about 2 times the average number concentrations. We view these clusters as various stages of aged aerosols, but less so than the last three clusters.

Clusters 10, 11 and 12 represent the most aged aerosol size distributions with the largest modal size centred around a dominating accumulation mode of about 170–200 nm. These clusters are mainly observed during summer months. Integrated amount of precipitation along the trajectories shows that clusters 10 to 12 have the lowest amounts, with cluster 12 having the least precipitation of all clusters (data not shown). Thus, these clusters represent aerosol size distributions that undergo ageing by condensation and coagulation with limited removal from the atmosphere. The size distributions in these clusters share features with aged continental aerosol observations, and as they represent the aerosol observations with the highest submicron mass, we will refer to these clusters as polluted.

Much of the reasoning above is based on knowledge of the tropospheric aerosol. Some additional support to the analysis was given by the precipitation history of trajectories belonging to the four major cluster groups. The precipitation intensity along trajectories belonging to each one of the major cluster groups washout, nucleation, intermediate and polluted was averaged and the result was a bar graph used to illustrate the resulting average precipitation during the last 10 d of transport. It can be clearly seen in Supplement Fig. S2 that the integrated precipitation is largest for washout and nucleation clusters and smallest for intermediate and polluted clusters. This suggest that the washout-type cluster is indeed more likely to have experienced a higher precipitation amount en route to Aspvreten. The nucleation-type clusters exhibit an interesting pattern: precipitation rate is on average high up to some 20 h prior to arrival after which average intensity decreases. We hypothesize that washout followed by a clearing of skies just before arrival paves the way for new particle formation, which in turn highlights the need for both relatively clean air (low condensation sink) together with high photochemical activity to create favourable conditions for new particle formation in the size range of a few nanometres. Hence, the epithet names given to the cluster types are consistent with the expected relation to the evolution of precipitation along the trajectories (Tunved et al., 2004).

Based on the analysis above, we can reduce the clusters into four broader categories: washout (cluster 1), which characterizes wintertime continental boundary layer (CBL) aerosol size distributions; nucleation (clusters 2 through 6), which characterizes new particle formation and growth; intermediate (clusters 7 through 9), which characterizes further growth and processing; polluted (clusters 10 through 12), which characterizes summertime aged clusters with a low influence of precipitation and strong influence from long-range transport of continental aerosol.

We will use these four groupings in linking clusters and their corresponding trajectory source regions. We will also use them in the evaluation of how cluster members have been shifted between the four types during the studied period.

3.3 Source regions of dominating cluster types

Dominating transport paths of air mass trajectories belonging to each of the major cluster groups are depicted in Fig. 6. As can be seen, the main difference in air mass source regions is found comparing the polluted and the nucleation cluster groups. Nucleation-type size distributions are commonly associated with northerly air mass transport between Greenland and Svalbard, while polluted-type clusters show pronounced southwesterly transport on average, with major pathways over Great Britain and northern parts of continental Europe. This is in agreement with previous findings (Tunved et al., 2003, 2005; Sogacheva et al., 2005, 2008), showing that nucleation observed over the Nordic region predominantly occurs in clean marine Arctic air masses. The transport characteristics for the intermediate-type size distribution clusters show resemblance to the transport pattern observed for the nucleation-type clusters, suggesting this type of size distributions indeed fits the picture of an aged nuclei mode resulting from nucleation taking place upstream of the receptor. The washout cluster type shares source regions with the polluted-type cluster, but with a high degree of wet removal masking the influence of sources apparent in the polluted cluster type.

Figure 6 Dominant source regions of air mass trajectories belonging to cluster groups denoted washout (cluster 1), nucleation (clusters 2 through 6), intermediate (clusters 7 to 9) and polluted (clusters 10–12). Aspvreten shown by the red dot.

4 Decadal trends in aerosol number size distribution properties

4.1 General trends of the aerosol size distribution

Figure 7 shows the estimated annual linear trend (and 95 % confidence interval, CI) for each of 16 size bins expressed as Theil–Sen estimator slope together with the grand average of the number size distribution with 25th–75th percentile ranges indicated. The Theil–Sen estimator represents the median of all slopes derived from pairs of values (i,j) from the data set, $(y_j-y_i)/(x_j-x_i)$, and constitutes a robust measure of trends for data series that are not normally distributed. For the statistical calculations we have made use of the ktaub statistical MATLAB package (Burkley, 2006) and references therein. Prior to trend analysis, all size distribution data were remapped to 16 size bins over the observed size range, using a bin width of dlog D_p equal to 0.1. The trend is expressed as the change in size distribution as d(dN∕dlog D_p) per annum. It is clear from Fig. 7 that a negative trend is present for all size classes. The nuclei mode particles exhibit only a modest downward trend (> −5 cm⁻³ yr⁻¹) and remain about the same up to 40 nm. Then the decreasing trend become gradually larger with increasing size, reaching peak values of −35 cm⁻³ yr⁻¹ around 80 nm. Hereafter the negative trend decreases in magnitude up to around 120 nm. Between 120 and 200 nm, the negative trend again becomes larger, reaching a second mode peaking at a downward trend of −28 cm⁻³ yr⁻¹. After this peak, the negative trend decreases with increasing size to the upper limit of the instrument at 400 nm. It is noted from Fig. 7 that the trend does not correlate with the average concentration of particles. Figure 7 also suggest two modes of the decreasing trend, one around 80 nm and one centred at 200 nm. The gradient in trend is further very sharp between 40 and 80 nm. The presence of the two modes could be suggestive of mainly two different processes acting independently on the size distribution. Sulfur emissions in Europe have seen a dramatic decrease during the last decades (Smith et al., 2011; Aas et al., 2019). Being one of the prime precursors for formation of aerosol number and mass formation, we attribute the strong negative trends around 80 and 200 nm to a reduction in aerosol sulfate. The first mode likely reflects the reduced gas-particle production occurring far from the receptor through condensation of sulfates. The second mode however is likely the result of less in-cloud oxidation of SO₂, a well-known production pathway for accumulation mode sulfate mass (e.g. Langner et al., 1992). The lack of any major trends present for particles smaller than 40 nm suggests that the period 2000–2017 has seen few changes in particle production close to the site. The two modes at 80 and 200 nm are however consistent with changes in processes occurring far away from the measurement site in air masses where the aerosol has had time to age from both condensation and in-cloud processes.

Figure 7 Calculated annual average trend per size bin for the studied period 2000–2017 (right y axis). The figure also shows median number size distribution and 25th–75th percentile interval for reference (left y axis).

Figure 8 Contour plot depicting the calculated trends (as Theil–Sen slopes) for different size classes between 10 and 390 nm. Colours indicate calculated linear trend for binned particle number concentration at Aspvreten as particles per cubic centimetre per year for the time period 2000–2017. Daily averaged data were used for calculation of the Theil–Sen slope. Areas enclosed by the dashed red line represent pairs of month–size bin where tests for significance failed.

The trend analysis was further extended to investigate if any seasonal variation could be found in the size-dependent trends. Figure 8 shows the seasonal size-resolved trend for the period. As can be seen, a rather complex picture emerges, where trends vary with both size and season. Some generalizations can however be made. During December, January and February all size classes show a decreasing trend. This is followed by the most substantial negative trend, found in March through May for the size range 50–300 nm. The trend, d(dN∕dlog D_p) per annum, is in excess of −40 cm⁻³ yr⁻¹. During the same time of the year there is a concurrent moderate increasing trend for particles smaller than about 40 nm (∼10 cm⁻³ yr⁻¹). During the summer period between mid-May and August, the negative trend for particles larger than about 70 nm is less distinct compared to the spring period. However, particles of 25–70 nm present a strong positive trend, and this trend is similar in magnitude but opposite in sign compared to the accumulation mode trend in spring (−30 cm⁻³ yr⁻¹). In autumn the positive trend for the small particles is abruptly changed to a negative trend, and essentially the whole aerosol size range displays decreasing concentration with time. The negative trend is strong in the size range of about 40 and 70 nm (around −40 cm⁻³ yr⁻¹). Evidently, strong gradients exist in estimates of trends over the observed size range and season, and same-sized aerosols can show both positive and negative trends, depending on season. In Fig. S4 the 95 % confidence interval of Theil–Sen slopes for lower (left) and upper (right) confidence intervals are shown. Colours indicate the calculated linear trend for binned particle number concentration at Aspvreten as particles per cubic centimetre per year for the time period 2000–2017. Areas bounded by the dashed red line represent pairs of month–size bin where the test for significance was below the 95 % threshold.

4.2 Trends of fitted size distribution properties

In order to attribute the trend information depicted in Fig. 8 to the different modes described in Sect. 2.2, trends were derived for each mode, now on an annual basis. From the information in Fig. 8, we are aware that the nuclei and Aitken modes will be a mix of positive and negative trends on an annual basis. Nevertheless, this procedure will help in establishing a confidence in the observed trends over the data set as a general feature rather than a step change due to changes in infrastructure.

Figure 9 Trends of fitted log-normal parameters 2000–2017 of nuclei, Aitken and accumulation modal diameter and modal concentration. Panels (a)–(c) show trends for fitted geometric modal diameter, and panels (d)–(f) show trends for fitted modal number concentration. Whiskers indicate the 25th–75th percentiles of data. Fitted trend (as Theil–Sen estimator) slope is given in each title. Upper and lower 95 % confidence intervals of slopes are indicated in figures. The CI of slopes is calculated following Hollander and Wolfe (1973).

Figure 9 shows the annual medians of modal diameter (Fig. 9a–c) in micrometres for modes 1–3 and the corresponding modal number concentration (Fig. 9d–f) in reciprocal cubic centimetres. Also indicated in the figure title is the Theil–Sen estimator of the slope of the geometric mean diameter and modal number concentration. All trends were tested using the seasonal Kendall test on the hourly data (95 % CI), and Table 3 summarizes the results from calculations of Theil–Sen estimator slopes together with confidence intervals of derived slopes for the whole year.

Table 3 Theil–Sen estimator slope of fitted variables D_g and N and associated upper and lower confidence intervals during the time period 2000–2017. Bold values indicate significance on the 0.05 level according to the seasonal Kendall test (based on 12 months). Calculated σ from seasonal Kendall test is indicated. Trend of D_g for nuclei mode is not significant on the 0.05 significance level. Upper and lower CIs calculated following Hollander and Wolfe (1973).

As can be seen in Table 3, the median trend for the nuclei mode on an annual basis is −10.9 cm⁻³ yr⁻¹ (95 % CI −8.5 to −12.2 cm⁻³ yr⁻¹). Aitken mode number concentration trend was found to be −9.2 cm⁻³ yr⁻¹ (95 % CI: −7.1 to −11.4 cm⁻³ yr⁻¹), and accumulation mode concentration trend was calculated to be −8.2 cm⁻³ yr⁻¹ (95 % CI: −7.4 to −9.1 cm⁻³ yr⁻¹). All trends of number concentrations were found significant when applying a seasonal Kendall test following the method proposed by Hirsch and Slack (1984). The choice of a seasonal Kendall test was motivated by the serial dependencies apparently present in the data set (i.e. seasonal variation). The upper and lower 95 % confidence intervals of slopes were calculated following Hollander and Wolfe (1973).

Regarding the trends of modal geometric diameters, calculations of the Theil–Sen estimator slopes gave only a slight negative trend for the nuclei mode of −0.07 nm yr⁻¹ (95 % CI: −0.00 to −0.14 nm yr⁻¹) and a little stronger decrease for the Aitken mode geometric mean diameter (−0.31 nm yr⁻¹, 95 % CI: −0.20 to −0.44 nm yr⁻¹). The most marked decrease was found for the accumulation mode geometric mean diameter. The average trend over the studied period was found to be −2.01 nm yr⁻¹ (95 % CI: −1.76 to −2.24 nm yr⁻¹). This means that since observations started in 2000, the annual mean modal size of the accumulation mode has decreased by almost 40 nm, conjunct with a number decrease of around 150 cm⁻³. These numeric values are very large compared to the mean seasonal values listed in Table 1 and indicative of major changes in both size and number concentration of the accumulation mode.

Figure 10 Theil–Sen estimator slopes for annual trend per month for fitted number concentration of nuclei, Aitken and accumulation modes, January–December 2000–2017. Indicated by error bars are upper and lower confidence intervals of slopes at σ(0.05). All trends are tested with the Mann–Kendall test and shown to be significant for σ=0.05 except the Aitken mode trend during June–July and November. Upper and lower CIs calculated following Hollander and Wolfe (1973).

A seasonal trend analysis for the modal number concentration is presented in Fig. 10, and a similar seasonal variation as depicted in Fig. 7 is apparent. The accumulation mode is negative for all months with the strongest decrease present during April, −13 cm⁻³ yr⁻¹. Over the rest of the year the decrease is in the range between −5 and −10 cm⁻³ yr⁻¹. Numerical values for each mode and month are summarized in Table 4.

Table 4 Theil–Sen estimator trend slopes (cm⁻³ yr⁻¹) and upper and lower confidence intervals for monthly inter-annual trends in modal number concentration 2000–2017.

Different from the accumulation mode, the nuclei and the Aitken mode trends are weakly positive in June and July. From a positive trend in summer, the nuclei mode presents the overall strongest monthly decreasing trend of the three modes in October, with about −20 cm⁻³ yr⁻¹. Except for the months of June, July and November regarding the Aitken mode, all other trends in Fig. 10 and Table 4 are significant according to the Mann–Kendall test at the 95 % confidence level.

Figure 11 Box–whisker plot of daily mean mass concentration 10–390 nm assuming density of 1 g cm⁻³ for Aspvreten 2000–2017. Theil–Sen estimator slope and associated 95 % CI of slope is indicated in the figure by solid black lines. Note that outliers are excluded in the figure for improved clarity. Upper and lower bounds of boxes show the 25th–75th percentile range, and whiskers show the 5th–95th percentile range of data. Red lines show the median.

As an additional analysis, the temporal evolution of aerosol submicron mass was investigated. The aerosol number size distribution was recalculated to mass size distribution assuming a density of 1 g cm⁻³, followed by integration over the instrumental size range. In this way mass was calculated as daily means (µg m⁻³) for the period 2000–2017. Theil–Sen estimator slope was calculated for this data set. Figure 11 presents the result of annual mass concentration as a box–whisker plot and the estimated Sen slope. As can be seen, the reduction of mass follows the general trend of aerosol accumulation mode number concentration closely. A Theil–Sen slope of 0.077 µg m⁻³ yr⁻¹ (95 % CI −0.085 to −0.071 µg m⁻³ yr⁻¹) was estimated. Based on this calculation, we estimate that mass was reduced by roughly 52 % during the studied period. Tørseth et al. (2012) estimated an average decrease in PM_2.5 of 27 % for 2000–2009 for several European stations. In Figs. 5 and 6 of Tørseth et al. (2012), average PM_2.5 mass and an estimated trend are visible for Aspvreten. PM_2.5 average concentration at Aspvreten is in the range of 4–7 µg m⁻³, and our calculated mass (assuming spherical particles and unit density) is about 2.2 µg m⁻³. This suggests that at least 50 % of PM_2.5 is controlled by particles smaller than approximately half a micron (this is a lower range estimate since we use density of 1 g cm⁻³, so doubling our assumption of density will put sub-390 nm mass in the range of 4.5 µg m⁻³). Our calculated decrease in mass for 2000–2009 is roughly 30 %, i.e. slightly less decrease in PM_2.5 than estimated by Tørseth et al. (2012), in excess of 40 %). This suggests that the negative trend of particles > 390 nm is slightly larger for particles above the size range covered by our instrument. Nevertheless, our estimated trend is still in good agreement with observed PM_2.5 trends and it is evident that a major part of PM_2.5 reductions is controlled by particles smaller than 400 nm. Furthermore, Glantz et al. (2019) found a significant decrease in aerosol optical thickness (AOT) over the western Gotland basin area of 1.5 %, 1.1 % and 1.6 % per annum derived from MODIS c051, MODIS c061 and AERONET Gustav Dalen Tower (GDT), respectively. This would translate to an overall decrease of about 25 % over the 17-year period covered by our study, which agrees well with our derived submicron aerosol mass decrease of ∼30 %.

4.3 Trends in cluster groups

Each size distribution observed at the station is attributed one of the 12 cluster types and further grouped into one of four categories (washout, nucleation, intermediate and polluted), as described in Sect. 3.2. For each year, the sum of the different categories represents 100 %. If the relative contribution of a category remains the same over the time period, there will be no significant trend for that category. If a category increases or decreases significantly, this must be reflected in an opposite change in one or more of the other categories. The trend analysis of the cluster groups is presented in Fig. 12.

Figure 12 Inter-annual trends of major cluster groups washout (cluster 1), nucleation (clusters 2 through 6), intermediate (clusters 7 to 9) and polluted (clusters 10–12). Fitted linear trend and annual decrease are indicated in figure.

The two clearest changes over time are the increase in the washout type (cluster 1) and a corresponding decrease in the polluted type (clusters 10–12), both at a rate of about 1 % units per year. A closer look at the polluted category shows that cluster 12 almost completely disappears over the measurement period. The two other categories show smaller changes. The nucleation presents a small decreasing trend of 0.1 % units per year, and the intermediate category shows an increase of about the same magnitude, 0.2 % units per year.

Whereas the washout category almost doubled over the measuring period, this was almost entirely compensated for by the decrease in the polluted category. These changes combined could conceivably lead to a lowering of the condensational sink (CS) by the aerosol (see Table 3), and this would tentatively promote new particle formation. However, based on the data from Aspvreten, the negative trend in the category nucleation does not support a simple direct link between CS and new particle formation at the station. The slight increase in the intermediate category could indicate that there was still a small enhancement in new particle formation, but that it must have occurred at some distance away from Aspvreten and had time to age several days.

4.4 Trends in potential CCN concentration

In order to study the potential effect of the observed changes presented above on availability of cloud condensation nuclei (CCN), each observed size distribution was used to initialize a simple adiabatic cloud parcel model. The model, which constitutes a module in the aerosol transport model CALM (Tunved et al., 2010), calculates cloud droplet activation in warm clouds (i.e. without ice component) under adiabatic conditions.

This is a highly generalized assumption that neglects the influence of competitive growth in mixed-phase clouds. The introduction of ice crystals into a pure liquid phase cloud causes the ice crystals to grow at the expense of the liquid phase due to the relatively lower saturation vapour pressure over ice compared to liquid droplets. In clouds that are warmer than −36 ^∘C, formation of ice crystal occurs on so-called ice nuclei, typically mineral dust and soot particles of larger diameters (DeMott et al., 2010), although other particle types may serve as ice nuclei as well (e.g. Murray et al., 2012). In Kanitz et al. (2011) Figure 3 presents the observed fraction of ice-containing clouds in different environments as a function of cloud-top temperature. As can be seen, there is a high degree of variation in the fraction of ice-containing clouds depending on environment studied. For the Leipzig data set, 70 % of the clouds contain ice at −19 ^∘C, while a similar ice fraction is only achieved at −34 ^∘C in the Punta Arenas data set.

If we assume our calculated estimate is restricted to low-level clouds, and given that the average surface temperature at Oxelösund (a few kilometres southwest of Aspvreten) is below 0 ^∘C during winter months only (Data source: Surface temperature, hourly values, from Swedish Meteorological and Hydrological Institute, SMHI, https://www.smhi.se/q/Stockholm/2673730, last access: 4 June 2019, 1961–2018), the problem arising from mixed-phase clouds is likely contained to the winter months December–January. Thus, our assumption regarding liquid-phase-only cloud droplet activation and growth is likely prone to errors during this season and is likely to overestimate the CCN number in our idealized simulations.

A simplified chemical system of internally mixed aerosols was assumed, consisting of ammonium bisulfate and some hypothetical organic compound with a solubility of 10 % and a van 't Hoff factor of 1. The van 't Hoff factor for NH₄HSO₄ is set to 2.4, and the solubility is assumed to be 100 %. Water accommodation coefficient is set to 1.

The evolution of the droplet spectra for each cloud simulation is calculated until peak supersaturation (S) is reached in the parcel, i.e. for as long as $\mathrm{d}S/\mathrm{d}t>\mathrm{0}$, to assure that no more activation can take place. Once calculations are terminated, the number of activated particles is integrated and stored as the CCN value for that particular size distribution. This approach will allow for a slight variation in the altitude where peak supersaturation is reached, and this in turn depends on updraught velocity and chemical makeup.

Each aerosol size distribution was tested using three different chemical characteristics based on NH₄HSO₄:organic ratios of 50:50, 90:10 and 10:90. For each chemical mixture three different updraught velocities were used, 0.1, 0.5, and 1 m s⁻¹, respectively. The nine different scenarios are summarized in Table 5.

Table 5 Median (and 25th–75th percentile) range (cm⁻³) of calculated CCNs for different updraught and chemical composition scenarios. See text for details. In the calculation of ΔN∕3N the year 2000 has been used as the reference year. Value used is the intercept at the year 2000 of the fitted line. ΔN is calculated as the Theil–Sen slope multiplied with 18 years.

The resulting annual median CCN for the 50:50 mixture and three different updraughts is presented in Fig. 13. All three updraughts give a similar impression where the CCN drops decreased by more than one-third of the original concentration during the observation period. The linear trends, calculated as Theil–Sen slopes, provide numerical values of about −5, −10 and −13 cm⁻³ yr⁻¹ for the low, medium and high updraught scenarios. The two other cases, with 90:10 and 10:90 mixtures, follow the negative trend as presented in Fig. 13, but with the difference of being parallel shifted, that is, stronger trends for the more soluble mixtures as opposed to weaker trends for the less soluble scenarios. All nine slopes are summarized in Table 5. The difference in trend is largest for the high-updraught scenarios comparing the different solubility assumptions. The trends are similar for the low-updraught case, −5.0 (−5 to −5), $-\mathrm{4}.$3($-\mathrm{4}.$8 to −3.9) and −3.5(−3.8 to −3.2) cm⁻³ yr⁻¹ for compositions of 50:50, 90:10 and 10:90, respectively. For the high-updraught case corresponding values are −13.3(−15.0 to −11.7), −15.0(−17.0 to −13.0) and 9.3(−10.3 to −8.4) cm⁻³ yr⁻¹.

Figure 13 Trend of calculated CCN for different updraughts assuming a chemical composition ratio NH₄SO₄:organic of 50:50. The figure shows the median calculated CCN per year and corresponding 25th–75th percentile interval. The calculated linear trend as Theil–Sen slope and associated 95 % confidence intervals are indicated by solid and dashed lines, respectively. Slope indicated in figure title. See text for further details.

It is possible to make simple estimates of changes in cloud optical thicknesses based on these trends. If the cloud water content in our clouds is unaffected by the decreasing trend in CCN, we can simply write

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\mathrm{\Delta }\mathit{\tau }}{\mathit{\tau }}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{\Delta }N}{N}},\end{array}$$

where the term on the left is the relative change in cloud optical thickness and the term to the right is our estimated relative change in CCN based on the analysis above divided by 3. The relative changes in cloud optical thicknesses are summarized in Table 5.

For all our scenarios, the expected relative increase in cloud optical thickness is greater than 9 %. The largest relative change in cloud albedo of 16 % between 2000 and 2017 occurs for the assumption of least-soluble particles and the smallest updraught velocity. The high sensitivity for these conditions is expected since both low updraught and less soluble particles will generate few cloud droplets and thus high sensitivity to even moderate changes in cloud droplet number concentration, all in agreement with Twomey (1991) and Platnick and Twomey (1994). If we assume that the most realistic composition is somewhere between the 50:50 and 90:10 chemical makeup, this would yield a cloud albedo change in the range of 10 %–12 %.

Figure 14 Calculated Theil–Sen slopes for the nine different cloud activation scenarios.

The annual trends provide a very robust picture of a decrease in potential CCN and the likely decrease that follows in cloud optical thickness. However, our analysis above also emphasizes that the calculated trends in size distribution properties are not the same over the whole year. We therefore repeat or model calculations on a monthly basis. Figure 14 combines all nine combinations of updraught and chemical mixture and the resulting trends in CCN as a function of month.

Figure 14 points to some interesting general features. Firstly, the strongest contributors to the overall decrease in CCN are March and April. The strong decrease in accumulation mode particles cannot be compensated for by the small local enhancement in nuclei mode (see Fig. 18). Secondly, despite the summer increasing trend in nuclei and Aitken modes, the net effect on CCN is negative. Thirdly, the strongest decreasing trend presented by the nucleation mode in September and October is not reflected in the trend of CCN in the same period, suggesting that CCN is controlled by larger particles transported to Aspvreten.

Overall, it is clear that a substantial reduction in CCN number is likely to have taken place during the studied period. In fact, if we assume the scenario 50:50, with an average updraught of 0.5 ms⁻¹, the 17-year period has brought about a change in CCN of as much as 360 cm⁻³ in the month of April. Of course, other changes different from the decreasing trend in aerosol concentration may also have taken place, which could have altered the potential availability of CCN.

4.5 Trends in general transport patterns

The indisputably and generally significant trends in the aerosol size distribution on an annual and monthly basis throughout the measurement period are expected as a response to a reduction in anthropogenic emissions of aerosols and their precursor gases, mainly sulfur dioxide. However, these analyses must be complemented with an analysis of possible variations in transport patterns over the same time period. To assess the contribution of observed trends from changes in air mass transport patterns, a polar coordinate grid system was arranged around Aspvreten spaced in 2.5^∘ increments in the lateral direction and 10^∘ in the longitudinal direction.

Each time a trajectory passed over a grid box was recorded as an event. For each year, the sum of all used trajectories for that year was used to normalize each of the grid boxes. This results in a new value associated with the grid box representing a yearly relative probability that a trajectory passed a given grid box en route to Aspvreten (compare with Fig. 6). For each grid box a 1st-order fit was performed to represent the trend using the polyval.m function in MATLAB. The result is presented in Fig. 15.

Figure 15 Trend in probability of transport across different source regions expressed as a linear fit of the trajectory probability function over the 18-year period. Units expressed as per mille.

The figure shows the sign (positive or negative) and magnitude of the trend of the probability of a random sampled trajectory to cross a given grid cell. A positive trend thus means that more trajectories are crossing that grid cell today compared to 17 years ago. All in all, given the length of the trajectories this reflects changes in the mesoscale circulation.

As can be seen in Fig. 15, the likelihood of direct northerly transport has decreased, while the southwesterly component has increased, resulting in more frequent air mass transport from the North Atlantic region, including areas around Greenland and Iceland. There is only very minor reduction in the amount of air transported from continental regions.

Despite the systematic change in the transport pattern presented in Fig. 15, the trends are small in comparison to the large changes in aerosol properties over the same time period. We cannot however completely exclude any possibility of transport-related changes in observed aerosol concentrations at Aspvreten. We do however believe that contributions from observed small changes in transport patterns to the very large observed trends are small.

We further attempted to relate the increased relative frequency of washout-type clusters to trends in precipitation amount along the trajectories. The analysis did not result in any clear time-dependent trends. In order further investigate the relative roles of wet removal and source strength, we revisited the data and made an analysis of average precipitation intensity along cluster 1 (i.e. washout-type cluster), creating hourly resolved bar graphs for pre-2009 and post-2008 periods, to see if there are any differences. If post-2008 data contain trajectories with on average more precipitation, enhanced wet removal may play an important role. It is clear that integral precipitation on average is about the same for cluster 1, but slightly different patterns in timing of precipitation can be noted. This tells us that at least for the last ∼2 d, we have, if anything, less precipitation influence. This in turn does suggest that the relative increase in size distribution observations belonging to cluster 1 is not likely to be the result of increase in precipitation, and if we assume that the sink strength would be the same, the redistribution of cluster members in favour of cluster 1 must be dominated by emission reductions and not enhanced sink processes.

It should also be added that there is a substantial uncertainty in calculation of trajectories, especially in regions where possibilities of meteorological model validation are poor like in the Arctic region. The results presented here do however indicate that if changes in transport pattern in fact take place, the changes over the studied period seem small or non-linear.

5 Discussion

In the presented study we have made use of almost 20 years of continuous data comprising aerosol number size distribution data between 10 and 390 nm. Data availability during the studied period is above 80 %. The data have been analysed from several different aspects in order to shed light on the degree of change and when and where changes in emissions may have had an impact on the observed aerosol size distribution. Results show generally negative trends for annual mean concentrations and sizes of the three modes, referred to as nuclei, Aitken and accumulation modes. Based on the combined cluster and air mass analysis, this decreasing trend, on the order of 10 particles per mode and year, is consistent with reduced emissions of sulfur dioxide and other anthropogenic primary particles as well as aerosol precursor gases. Aerosol integral number has been reduced by roughly 30 %, from ca. 1800 cm⁻³ in the beginning of the period to about 1300 cm⁻³ at the end of the studied period. The general trend found in this study is in good agreement with the findings presented by Asmi et al. (2013), and extrapolation of their derived trends (N20, Asmi et al., 2013, Table 2) over the same time period presented in this study would give results in the range of 23 %, 20 % and 40 % at Vavihill, Hyytiälä and Värriö and an increase of 2 % at Pallas. A simple estimate of mass trends based on the presented trends for N100 yields values similar to our estimates of submicron mass trends found in the current study. In addition, cluster analysis of hourly number size distributions has been demonstrated to be a useful tool in trend studies. The method has been applied to study how the aerosol observations are distributed over 12 signature distributions and further applied to investigate how the aerosol has been re-distributed between these 12 dominating cluster types during the period 2000–2017. The method seems well suited for studying trends in new particle formation events. It is however clear that there is more to the story than simply an overall decrease over all seasons and over the entire size range. We have shown that the trends can be negative or positive depending on what season and over which size range the trend is calculated. This is indicative of the complex interplay between aerosol sources and sinks, their seasonal variation, and the dynamics of new particle formation.

In the case of the accumulation mode, fitted modal diameter has decreased by 2 nm yr⁻¹ (in total, a 17 % decrease or ∼35 nm absolute decrease; see Fig. 9). This decrease, together with the accumulation mode number decrease, has been associated with a major reduction in submicron mass. We show that integrated particle mass of 10–390 nm has decreased by more than 50 % during the studied period (as annual averages; see Fig. 11). This decrease is in good agreement with observed trends of PM_2.5. We further show that a substantial part of PM_2.5 mass is made up by particles smaller than 400 nm. This in turn suggests that observed PM_2.5 changes to a large degree are controlled by particles smaller than about 0.5 µm. We also show results suggesting that two main sulfate production pathways have decreased in magnitude, namely gas-phase oxidation followed by condensation and wet-phase sulfate production in non-precipitating clouds.

More complex is the behaviour of the Aitken and nuclei modes. Since the nuclei mode reflects recent formation of particles, the general reduction in nucleation mode particles must reflect reduced particle formation in the relative vicinity of the station. However, as particle formation is largely controlled by the condensation sink (CS) of pre-existing particles (Boy et al., 2005; Dal Maso et al., 2008; Kulmala et al., 2001), and the accumulation mode mainly controls CS, a reduction in accumulation mode number and size will also decrease the condensation sink. Therefore, one would expect a reduction in CS to promote new particle formation. As shown by the calculated trends, a general increase in nuclei mode is not observed. This indicates that the reduction in CS is accompanied by a concurrent reduction of precursors available to nucleate new particles, which counterbalance the effect of a reduced CS. As sulfuric acid is one of the main key candidates for new particle formation, we conclude that the reductions of sulfur emissions thus feed back through both reduced CS and reduced nucleation. Moreover, there is likely no other nucleating species that appears to have replaced sulfuric acid with respect to new particle formation in the studied region. In that case, we would have expected a stronger response in the nuclei mode due to the strong changes in CS. Another possibility that we cannot rule out is that changes in precipitation during the studied period have resulted in more efficient removal of both accumulation mode particles and nucleating species. Recalling Fig. S3 in the Supplement, any trend as indicated by trajectory-derived precipitation includes a decrease or redistribution along the transport paths comparing pre-2009 to post-2008 periods. As we have not explicitly studied cloudiness over the site, we also can not rule out the possibility that on average increased cloudiness has reduced the photolysis rate close to Aspvreten. The role of both precipitation and cloudiness in changes of new particle formation would likely deserve a dedicated study. In Fig. 8 we show that a positive trend for Aitken mode size particles is present for the summer months June–July. This increase may reflect improved conditions from a new particle formation perspective, and this may reflect favourable relation between H₂SO₄ and CS or some other nucleating species. As this increasing trend is absent for the smallest size classes, we conclude that if this increase is the result of enhanced nucleation, we argue that it must take place upstream of the station rather than in the immediate vicinity in order to reach the size range where we observe the positive trend. In essence, nucleation occurs earlier during the marine–continental transition and this in turn slows further nucleation closer to the station. We have also performed simulations to assess the potential influence of observed trends on cloud droplet number concentration. We show that a negative trend is present for all months, but the calculated decrease is largest during spring. The springtime reduction ranges from −5 to $\sim -\mathrm{30}$ cm⁻³ yr⁻¹ depending on assumptions regarding updraught and chemical composition. We estimate that the reduction of potential CCNs corresponds to a reduction in cloud albedo of about 10 %.

We also studied the potential influence on observed trends resulting from changes in transport pattern. We do however conclude that these are small and likely of less importance for the observed reduction in particle number and mass. Instead we suggest that emission reductions of anthropogenic emissions, notably SO₂, are the main driving candidates for the observed change. These changes have the largest effect during the spring, and this period is most sensitive to changes in anthropogenic emissions in absolute terms. Further, the sensitivity is also largest for the accumulation mode number concentration and size, which in turn represents the most efficient light-scattering and cloud-forming aerosol particles.

6 Summary and conclusions

In summary, our main findings are presented as bullet points below.

• As revealed by the aerosol size distribution clustering, it is evident that the cluster representing clean, cloud-processed aerosol is increasing at the expense of the polluted-type monomodal size distribution. There is only a marginal increasing trend of cluster members belonging to clusters showing signs of recent new particle formation. At the same time, the most polluted cluster 12 has been reduced from around 5 % of observations to around 1 % during the period 2000–2017. We have shown that cluster analysis can be successfully used to study aerosol trends. Further we have shown that cluster analysis also can be used for studying trends of intermittent processes such as new particle formation events. The method has clear advantages compared to standard time-averaging techniques as it preserves the shape and number concentration of the aerosol number size distributions within the clusters, which otherwise can easily get lost in the usually applied time averaging such as mean and medians.

• There is on average a negative trend of 10 particles cm⁻³ yr⁻¹ for nuclei, Aitken and accumulation modes. During the 2000–2017 period, annual average integral number concentration decreased by 20 %–30 %.

• Integral number concentration decreased from about 1800 to about 1300 cm⁻³, or by almost 30 %

• Annual mean accumulation mode size decreased by 35 nm or 17 % during the 2000–2017 period. This suggests that the reduced emissions result not just in fewer particles, but also that they are, on average, smaller.

• Aerosol mass (10–390 nm) decreased by more than 50 % during the 2000–2017 period. This trend is found to be in good agreement with trends derived for PM_2.5. The results further suggest that a large part of previously presented PM_2.5 trends are controlled by particles smaller than 0.5 µm.

• There is a pronounced seasonality in the trends of aerosol number concentration. The strongest reduction is found during spring in the accumulation mode and during autumn for the Aitken mode. In spring, the annual trend in accumulation mode number concentration is found to be close to 15 cm⁻³ yr⁻¹ and the reduction in nuclei mode concentration during autumn is calculated to be around 20 cm⁻³ yr⁻¹ during September and October.

• Although we do not find evidence for increased nucleation close to the measurement site, we do find a positive trend for Aitken mode sized particles during summer months. Since the accumulation mode concentration is decreasing, we interpret the increase in Aitken particle concentrations as a response to the decrease in condensation sink, leading to more frequent nucleation upstream of the station (see Sect. 5).

• A calculation of potential CCN using an adiabatic cloud parcel model for various assumptions regarding updraught velocity and chemical composition showed that the number of CCN likely substantially decreased during the 2000–2017 period. This in turn yields as a best guess about a 10 %–12 % decrease in cloud optical thickness. The negative trend is present for all months and for all assumptions regarding updraught and chemical composition, but the largest negative trend is found during spring, ranging from −5 to almost −35 CCN cm⁻³ yr⁻¹, depending on assumptions applied regarding chemical composition and updraught velocity.

• We suggest that main drivers for the observed reduction in aerosol number and mass are reductions of anthropogenic emissions in general and reduced SO₂ emissions in particular.