Introduction

Aerosols are one of the most important, yet uncertain, factors in determining the surface energy balance in the Arctic1. They can directly affect the surface energy budget through scattering and absorption of light (direct radiative effects) and indirectly through the modification of cloud microphysical properties (indirect radiative effects)2,3. However, the processes of natural Arctic aerosols and their effects on the radiative balance are poorly characterized in global climate models, thus a more complete understanding of the factors affecting the aerosol burden is pertinent for the prediction of the future Arctic climate4.

In the Arctic, there are strong seasonal variations in aerosol particle number size distributions and concentrations, with long-range transport of accumulation (Acc) mode particles (particles with a diameter (Dp) > 100 nm) in late winter and early spring and local production of ultrafine (UF) mode particles (Dp < 100 nm) in the summer and autumn seasons5,6,7,8,9,10. Expansion of the polar dome, limited turbulent transfer due to stably stratified atmospheric layers, and inefficient wet removal in the late winter and early spring allow for the long-range transport and accumulation of anthropogenic pollution (e.g., sulfate and black carbon) in the lower Arctic11,12,13,14,15,16, resulting in the Arctic haze phenomenon17.

During the summer, contraction of the polar dome, increase in wet removal, and retreating sea ice cause regional sources of aerosols and precursor gases to become more dominant, which is manifested in an increase in nucleation and Aitken mode particles through new particle formation (NPF) and subsequent growth8,18,19,20,21,22. Marine and cryosphere emissions of dimethyl sulfide (DMS) and iodine-containing compounds can photochemically oxidize into low-volatility vapors (e.g., sulfuric acid and iodic acid), which can initiate new particle formation and undergo further condensational growth to sizes that can affect cloud condensation nuclei (CCN) concentrations23 and thus the atmospheric radiative budget indirectly24,25,26. These emissions have been shown to proceed under biotic (marine biological activity in the marginal ice zone) and abiotic (photochemically induced snowpack emissions) conditions27,28. Biogenic sources can also contribute to the atmospheric burden of aerosol particles through the direct emission of primary particles, especially in the summer months29,30,31,32.

Changes in the atmospheric burden of aerosol particles have been observed within the Arctic. Recently, Heslin-Rees et al.33 showed that coarse mode aerosols (most likely of sea spray origin) have increased over the past 2 decades (1999–2016) at Zeppelin Observatory on Svalbard. This was mainly due to changing air mass circulation patterns, with air masses more frequently arriving from the Northern Atlantic, with sea ice decrease having only a marginal influence. Collaud Coen et al.34 recently demonstrated that aerosol radiative properties in the polar regions showed a mix of increasing and decreasing trends, which were seasonally and spatially inhomogeneous. Moffett et al.35 analyzed summertime particulate methanesulfonic acid (MSA) and non-sea-salt sulfate (nns-SO42−) concentrations at Utqiaġvik (formerly known as Barrow), AK from 1998–2017, an increase of 2.5 and 2.1% yr−1 was observed for MSA and nns-SO42−, respectively. Sharma et al.36 showed that during January–March at Alert, CA, anthropogenic aerosol constituents (SO42−, H+, NH4+, K+, Cu, Ni, Pb, Zn, non-soil V, non-soil Mn, and equivalent black carbon) have decreased between 23 and 80% between 1980 and 2013 due to emission reductions in Eurasia. Subsequently, NO3 has increased by 20% while aerosol acidity (defined as the concentration of H+)36 has decreased by 76%. Sea salt components Na+ and Cl have increased by 19 and 43%, respectively, during November–March. Schmale et al.37 analyzed nine aerosol chemical species at 10 Arctic sites, they found the majority of significantly decreasing trends are occurring during the winter (due to reductions in anthropogenic emissions) while the summer period displayed a mix of negative and positive trends. Interestingly, during June–August at Alert, CA, MSA decreased by 0.28 ng yr−1 from 1981–1999, and no significant change after 2000 was observed, while at Kevo, Finland, MSA increased by 0.326 ng yr−1 from 1965–201037. These studies highlight how aerosol properties, both microphysical and chemical, have changed in recent decades, both in different regions and periods, and emphasize the importance of long-term measurements of aerosol microphysical and chemical properties in the High Arctic.

With the Arctic region experiencing drastic changes in the recent decades, illuminating the effects of sea ice retreat, temperature increase, atmospheric transport patterns, and changing biogeochemical cycles, on the formation, growth, and removal of natural aerosols remains a significant hurdle for experimentalists and climate modelers22. There remain key knowledge gaps concerning Arctic aerosols, which need to be addressed to understand the role of anthropogenic and natural aerosols in aerosol–radiation–cloud interactions, and how this relates to Arctic climate change. Therefore, a detailed investigation of the changes in aerosol particle properties (number concentrations and types) and air mass history in the High Arctic was undertaken. This study aims to accomplish the following objectives: (1) determine how the aerosol particle number concentration has changed during the past decade, (2) examine which aerosol types are responsible for this change, and (3) investigate which processes are driving this change.

Results and discussion

Seasonality of PNSD, clusters, and air mass history

The particle number size distribution (PNSD) and number concentration of total, UF mode, and Acc mode aerosols at Villum experience a strong seasonality, as evident from Fig. 1a and b, for the years 2010–2018 (Supplementary Fig. 1 displays the data capture for this period). During the winter, the total number concentration is low, although, during late winter/early spring, Acc mode particles start to build up (Fig. 1a). This is due to inefficient wet removal (Supplementary Fig. 2) and expansion of the polar dome, which enables the transport and subsequent accumulation of anthropogenic emissions from the mid-latitudes to the Arctic38. Acc mode particles peak in early spring, which drives the local maximum of total number concentration in April5,6,7,8,11,17, although UF mode particles make a small contribution due to the return of photochemical activity8,21,26. During May, when the polar dome contracts, long-range transport loses its importance, while wet deposition becomes more prominent (due to the transition from less-efficient ice-phase cloud scavenging to more-efficient warm cloud scavenging39), and thus accumulation mode particles significantly diminish. During the summer, in the near absence of local sources of Acc mode particles, UF mode particles increase in concentration. The increase occurs in air masses confined to the High Arctic, where there is little to no influence of anthropogenic sources. An increase in biogenic emissions of DMS and organics from the freshly exposed open ocean and marginal ice zones, which are oxidized in the atmosphere, drive NPF and growth19,40,41,42,43. During the autumn, precipitation along the trajectories reaches a maximum at Villum (Supplementary Fig. 2), which causes significant removal of aerosol particles which create extremely pristine conditions9.

Fig. 1: Overview of aerosol particle properties and air mass history.
figure 1

Annual cycle of (a) daily median PNSD expressed as a 3-D contour plot and (b) daily median number concentration of total (red), UF mode (blue), and Acc mode (black) shown as the solid lines, along with the 25th/75th percentiles displayed as the shaded area, (c) PNSD for each cluster with the red line representing the median and the blue shading the 25th/75th percentiles, (d) the monthly occurrence of each cluster, (e) the monthly occurrence of each surface footprint type, and (f) the occurrence of each surface footprint type for each cluster. The data period for each figure spans 2010–2018.

The k-means clustering algorithm was applied to investigate the different sources and processes responsible for the observed PNSDs at Villum. The median PNSD for each cluster and their monthly occurrence are displayed in Fig. 1c and d, respectively. Eight clusters were found to describe the full range of aerosols types observed at Villum (Fig. 1c): five that peak in the UF mode range (Pristine, Nucleation, Bursting, Nascent, UF Bimodal) and three that peak in the Acc mode range (Acc Bimodal, Haze, and Aged). The UF mode clusters occur most frequently during the summer and autumn months, while the Acc mode clusters occur primarily during winter and spring (Fig. 1d). These are the same clusters found in previous studies at Villum albeit over shorter periods and using different temporal averaging (daily versus hourly in this study)21,44,45. For a more in-depth description of these clusters refer to Supplementary Note 1.

The air mass history for each cluster, and their seasonality, can help clarify the origin of the PNSD clusters. The monthly occurrence of each surface footprint type and the occurrence of each cluster is displayed in Fig. 1e and f, respectively. For a schematic of the air mass history calculation see Supplementary Fig. 3. During winter and early spring, air masses are above the mixed layer (AML) most of the time, when air masses reside below the mixed layer, they traverse mainly over sea ice and snow with minor contributions from the sea. This minor contribution from the sea likely reflects the longer distances traveled by trajectories during the winter months38, while the majority contribution from above the mixed layer reflects the expansion of the polar dome and accompanying long-range transport. During early spring (March and April), similar transport patterns are observed compared to winter, as the polar dome reaches its maximum. In May, when the polar dome starts to recede, air masses spend considerably less time above the mixed layer and more time over sea ice. During the summer, air masses spend most of the time over sea ice and snow with increased contributions from the sea and land, likely due to retreating sea ice and melting snow, as well as air masses being confined to the Arctic region. During autumn, air masses experience the maximum time over the sea (reflecting the minimum sea ice extent) and the time above the mixed layer starts to increase (reflecting the expansion of the polar dome).

The particles associated with the Acc mode clusters (Acc Bimodal, Haze, and Aged) and the Pristine cluster spend most of the time above the mixed layer, with contributions from sea ice and snow when below the mixed layer. For the Acc mode clusters, this reflects the influence of long-range transport on the origin of these clusters, while for the Pristine cluster, this shows the influence of cloud processing (indicated by the Hoppel minimum at 60 nm,) and wet removal (indicated by the low number concentration for this cluster, Fig. 1c). The particles associated with the UF mode clusters (Nucleation, Bursting, Nascent, and UF Bimodal) spend less time above the mixed layer compared to the Acc mode clusters (although they still experience substantial time above the mixed layer, ~50–59%), and increased time over sea ice, sea, and land. These clusters are associated with surface emissions from marine (sea) and cryosphere environments (snow on land and sea ice) with a small contribution from terrestrial sources (land without snow), supporting previous observations21,28,46. While these clusters are mainly connected with surface processes, no one single surface type dominates and air masses from the free troposphere still make a sizable contribution. This indicates that several sources and processes are contributing to the formation and growth of aerosol particles at Villum, such as marine and cryosphere emissions26,28,47 and possibly entrainment or cloud-mediated transport from the free troposphere48,49, albeit unlikely50,51. The distribution of each surface footprint type for the PNSD clusters largely reflects the seasonality in the cluster occurrence (Fig. 1d) and the surface footprint types (Fig. 1e), with Acc mode clusters being long-range transported and UF mode clusters being affected by local processes within the Arctic.

Trends in aerosol particle number concentration, PNSD clusters occurrence, and air mass history

To address research objective 1 (determine how the aerosol particle number concentration has changed during the past decade), we explored changes in the total (9–915 nm), UF (Dp < 100 nm), and Acc (Dp > 100 nm) mode number concentration. When referring to a direction and magnitude of a trend, it is implied that said trend is statistically significant. The total, UF, and Acc hourly number concentration was temporally aggregated to daily resolution using percentiles from 1 to 100. These daily percentile concentrations were segmented into the meteorological seasons and analyzed for trends using the 3PW algorithm52. The results are shown in Fig. 2a–c. The use of percentiles as an aggregation method allows us to inspect changes in the distribution of daily number concentrations (e.g., is there any trend for the lower or upper concentration regimes of the daily number concentration), which otherwise would not be evident using a single aggregation method (i.e., daily median). The lower and upper limits of the 95th% confidence intervals for these trends are displayed in Supplementary Figs. 4 and 5, respectively.

Fig. 2: Seasonal trends of particle number concentrations.
figure 2

Trends in the different percentiles of the daily aerosol particle number concentrations in different meteorological seasons for (a) total, (b) UF mode, (c), and Acc mode number concentrations. The slope is displayed by the color bar and percentiles with significant trends are indicated by the black and white dots.

From Fig. 2a, it is evident that the total aerosol particle number concentration is significantly increasing across the majority of percentiles for the spring and summer and only the top few percentiles during autumn (75th–88th and 91st–94th), although these trends for the different seasons are not of a similar magnitude. The slopes of the total number concentration trends are highest during summer, followed by spring and autumn, and the slopes are largest at the highest percentiles for all three seasons. This indicates that the upper concentration regime of the daily aerosol particle number concentration is increasing at a higher rate than the middle or lower concentration regimes. Given the similarity in the percentiles that show significant trends and similarities in the magnitude of the slopes during spring (total: 2.2–6.4 cm−3 yr−1, UF: 0.7–3.9 cm−3 yr−1) and summer (total: 3.8–18.1 cm−3 yr−1, UF: 2.7–15.5 cm−3 yr−1), the increase in total number concentration is likely driven by increases in the UF mode number concentration (Fig. 2b). Interestingly, during the autumn, there are more significant trends in the lower and middle percentiles for the UF mode (10th–22nd and 56th–95th) than for the total number concentration (75th–88th and 91st–94th). Total number concentration slopes (1.33–2.20 cm−3 yr−1) and UF mode number concentration slopes (0.6–1.8 cm−3 yr−1) are both low during the autumn. For the Acc mode number concentration, there are only a few significant trends in the upper percentiles during spring (86th–88th) and summer (91st–100th) and the magnitude of the slope is lower compared to the other size ranges (spring: 1.17–1.23 cm−3 yr−1, summer: 1.15–1.23 cm−3 yr−1) (Fig. 2c). The increase in the upper regimes for the total and UF number concentrations suggests an increase either in the frequency or intensity of NPF events, which contribute greatly to the number concentration at Villum8.

While it is clear that the total and UF number concentrations are increasing during spring and summer and this increase is strongest in the upper regimes, the mechanisms behind these increases remain unclear. Therefore, to address research objective 2 (examine which aerosol types are responsible for this change), we examined the changes in the temporal occurrence of the PNSD clusters. The results are presented in Fig. 3. These trends were calculated using the monthly contribution of each cluster (as described in the Methods section), as the daily contribution produced too many zero values for a reliable trend calculation. The only clusters to exhibit significant trends are Nucleation during summer (1.51 [0.04–2.31] % yr−1) and Pristine during autumn (1.41 [−1.68–3.36] % yr−1). It should be noted that the k-means clustering algorithm normalizes the PNSD to the vector length53, thus each cluster is defined by the shape and not the magnitude, i.e. total particle concentration of the PNSD cluster. Therefore, this trend analysis informs us about which PNSD clusters are occurring more or less often and not the magnitude of the PNSD in each cluster. The increasing trend for Nucleation during summer suggests the increase in the upper regimes of the total/UF mode number concentrations is likely due to an increase in the frequency of NPF events, although further research is required to elucidate whether the intensity of NPF events is increasing or not.

Fig. 3: Seasonal trends of PNSD cluster occurrence.
figure 3

Trends for the monthly occurrence of each cluster for each meteorological season. Boxes represent the Theil Sen Slope and red error bars represent the 95th% confidence intervals. Blue boxes signify statistically significant trends on the 90th% confidence level and black boxes are not statistically significant.

The eight-year average surface footprint analysis for each cluster provided us with information regarding the air mass history for each cluster. Similarly, trends in each surface footprint type can help elucidate which processes might drive the trends for each cluster. The results of the trend analysis for the surface footprint types are presented in Fig. 4. Analogous to the trends in the PNSD clusters, these trends were calculated using the monthly contribution. The time air masses spend over the sea is increasing during spring (0.35 [−0.002–0.93] % yr−1) and summer (0.69 [0.15–1.23] % yr−1). The geospatial coverage of these surface footprint types and their trends will be discussed below. The time air masses spent above the mixed layer is decreasing during autumn (−1.58 [−3.53–1.17] % yr−1). The trend in accumulated precipitation was calculated using the daily median as well as for each percentile (same methodology as used for the trends in number concentration). The daily median of accumulated precipitation is increasing during autumn (0.29 [0.14–0.45] mm yr−1) (Fig. 4). Trends and their associated confidence intervals for the different percentiles of accumulated precipitation are shown in Supplementary Fig. 6. Accumulated precipitation is increasing across all percentiles (0.06–0.56 mm yr−1) during autumn, with the largest increase seen in the upper regime (Supplementary Fig. 6). Similar to the increase in the upper regimes of the UF/total number concentrations, this suggests that precipitation along the trajectory path is either becoming more intense or more frequent. This increase in precipitation during the autumn has also been observed at Zeppelin Observatory (~5% yr−1) during autumn from 1999 to 2016 by Heslin-Rees et al.33, who applied a similar approach to the one utilized in this study (integrating precipitation along HYSPLIT produced trajectories). While Heslin-Rees et al.33 used a different meteorological field to drive the HYSPLIT model than this study (GDAS1 vs NCAR/NCEP), this gives confidence that the increase in accumulated precipitation is realistic, it also highlights the need to produce accurate meteorological fields which will aid researchers in future atmospheric studies. This increase in Arctic precipitation during the autumn is also predicted by state-of-the-art climate models54. The summary statistics (slope, confidence intervals, and the number of data points used in each calculation) for each variable and season exhibiting a statistically significant trend are shown in Table 1.

Fig. 4: Seasonal trends of air mass history.
figure 4

Trends for the monthly occurrence of each surface footprint type when air masses were below the mixed layer, above the mixed layer, and the daily median of accumulated precipitation for each meteorological season. Boxes represent the Theil Sen Slope and red error bars represent the 95th% confidence intervals. Blue boxes signify statistically significant trends on the 90th% confidence level and black are not statistically significant.

Table 1 Summary of the statistics for statistically significant seasonal trends.

Relationships between aerosol particle number concentration, PNSD clusters occurrence, and air mass history

To examine the effects of air mass history (time spent over the sea, accumulated precipitation, and time spent above the mixed layer) on the PNSD, the median total/UF mode number concentration and PNSD was calculated for percentiles of the air mass history parameters (Fig. 5), a correlation analysis was performed (see Supplementary Note 2), and boxplots for deciles of air mass history parameters were calculated.

Fig. 5: Relationships between PNSDs and air mass history.
figure 5

Median 2-D number size distribution for percentiles of (a) time spent over the sea during summer, (b) accumulated precipitation during autumn, and (c) time spent above the mixed layer during autumn. The 3-D size distributions from (a–c) are displayed in (d–f). The median Nucleation cluster size distribution during summer is shown in red in (a) and the median Pristine cluster size distribution during autumn is shown in gray in (b) and (c). The median UF (black line) and total (black line with diamonds) number concentration is displayed on the right axis for (df).

During spring, a correlation analysis yielded a significantly positive relationship between the time air masses spent over the sea and total (R = 0.56, p < 0.005) and UF mode (R = 0.69, p < 0.001) particle concentrations (Supplementary Fig. 7). To examine this relationship in more detail, deciles of the time air masses spent over the sea were calculated and the total and UF mode particle concentrations were binned into these deciles. Supplementary Figures 8 and 9 show that as the time air masses spent over the sea reaches the higher deciles, the distribution of total and UF mode number concentration increases during the spring and summer, respectively. Thus, the increase in total and UF mode number concentrations during spring can likely be attributed to the increase in time air masses spent over the sea.

During summer, a correlation analysis between the monthly occurrence of the Nucleation cluster and time spent over the sea exhibits a moderately positive correlation (R = 0.46, p < 0.05). The Nucleation cluster is associated with high amounts of total and UF mode aerosol particle number concentrations both in spring and summer (Supplementary Figs. 10 and 11, respectively). As the percentiles of time spent over the sea increase during summer, the median PNSD experiences a decrease in concentration for particles >100 nm and an increase in particles <100 nm (Fig. 5a). The median Nucleation cluster PNSD during summer is presented for comparison. Air masses associated with our Nucleation cluster spend most of the time over the sea and below the mixed layer during summer compared to other UF mode clusters (Fig. 1f and Supplementary Fig. 12). While the PNSDs for the upper percentiles of the time spent over the sea closely resemble that of the Nucleation cluster >100 nm, the PNSDs in the UF mode deviate from the Nucleation cluster <100 nm. There is a large increase in the concentration of particles between 30 and 40 nm for the upper percentiles of time spent over the sea as well as a slight increase for particles >20 nm (Fig. 5a). This indicates that as the time spent over the sea increases so does the number of ultrafine, and thus total, particles (Fig. 5d). It also indicates that the time air masses spent over the sea results in a lowering of the condensation sink (CS) (this relationship is also displayed in a 3-D format in Fig. 5d). The lowering of the CS for air masses below the mixed layer during summer is likely caused by a combination of diminished sources and enhanced sinks. As air masses spend more time below the mixed layer, they are less likely to be influenced by long-range transport of accumulation mode particles (diminished sources) and are more likely to be subject to dry and especially wet deposition (enhanced sinks)10. Wet deposition occurs through in-cloud scavenging (cloud droplet activation, in-cloud impaction, and diffusion scavenging) and below-cloud impaction scavenging55. Tunved et al.18 showed precipitation to be a major sink of Acc mode particles and a source of UF mode particles at Zeppelin Observatory during the summer. This results from Acc mode preferentially forming CCN (in-cloud scavenging) and their removal creates conditions conducive for NPF as well as the inefficient in-cloud scavenging of UF mode aerosols56,57. Freud et al.10, showed that during the summer an increase in the time spent within a cloud decreases the Acc mode number concentration and that an increase in time spent within the boundary layer resulted in a decrease in Acc mode number concentration. These results are valid not only for Villum but for several other Pan-Arctic stations. Supplementary Figure 3 shows the seasonal cycle of accumulated precipitation and the percentage of time the trajectory RH was above 95% (a proxy for time spent in-cloud), respectively. While accumulated precipitation peaks during the autumn, the in-cloud proxy peaks during the summer. In-cloud scavenging has been shown to be the more efficient removal mechanism for Acc mode particles over UF mode particles39,58 while below cloud scavenging is most efficient for UF mode particles compared to Acc mode particles59. This likely points toward in-cloud scavenging being the dominant sink, although we cannot draw robust conclusions at this time from our dataset and further research is required to fully elucidate the removal mechanisms of aerosols during the High Arctic summer. The reduction in the CS as air masses spend time below the mixed layer and over the sea is therefore likely a result of diminished sources and enhanced sinks, notably in-cloud scavenging60.

While a reduction in the CS sets the stage for NPF, other environmental parameters are likely required to reproduce the Nucleation cluster PNSD. Other factors that have been shown to influence the aerosol particle number concentrations in the Arctic include increased emission strength of precursor vapors, or the rate of oxidation of said vapors25,26, increased occurrence of low-level inversions50, or decreased cloud cover over the Arctic19. Interestingly, emissions of DMS, north of 70 °N are increasing in recent decades61 and downwelling shortwave radiation has been shown to be increasing on Svalbard62. Both of these trends were calculated outside of the period in this study, so while they are expected to affect aerosol formation as well as aerosol–cloud–radiation interactions their effect on aerosol particle properties at Villum remains to be seen during this period (2010–2018). Further analysis is required to elucidate the effects of the above-mentioned phenomena on total/UF mode number concentrations as well as the frequency and intensity of new particle formation and growth at Villum.

Figure 5b shows the effects of accumulated precipitation on the PNSD during autumn. It can be seen that as accumulated precipitation reaches the upper percentiles, the PNSD begins to closely resemble the median Pristine cluster PNSD during autumn, although less so for the sub-30 nm size range. The monthly Pristine cluster occurrence and monthly mean accumulated precipitation display a moderately positive correlation (R = 0.59, p = 0.013, Supplementary Fig. 7) and the period of peak Pristine cluster occurrence and peak accumulated precipitation are both during the autumn (Fig. 2d and Supplementary Fig. 2). Given the positive correlation and similarity between the PNSD for the upper percentiles of accumulated precipitation and the median Pristine cluster PNSD, the increase in the occurrence of the Pristine cluster is most likely due to the increase in accumulated precipitation during autumn.

The effect of time spent above the mixed layer on the PNSD and UF/total number concentrations during autumn is displayed in Fig. 5c and f. As the time spent above the mixed layer reaches the upper percentiles during autumn, the PNSD strengthens in the Acc mode, although the number concentration generally decreases (except for an increase around the ~85th percentile). Conversely, as the time spent above the mixed layer decreases, the UF mode number concentration begins to increase (which is driving a similar pattern for the total number concentration), indicating a surface source of UF mode particles (Fig. 5f). The distribution of UF mode number concentration binned by deciles of the time spent above the mixed layer supports this posture, with increasing UF mode concentrations with decreasing time spent above the mixed layer during the autumn (Supplementary Fig. 13). This observation, combined with the negative correlation displayed in Supplementary Fig. 7, suggests that decreased time spent above the mixed layer is responsible for the increase in UF mode particles during autumn. While total/UF mode number concentrations both slightly increase with increasing percentiles of accumulated precipitation (Fig. 5e and Supplementary Fig. 14) concentrations were consistently low for all percentiles. Furthermore, no correlation was found for these parameters during autumn, therefore, there appears to be little effect of the amount of accumulated precipitation on the total/UF mode number concentrations. Increased precipitation will remove Acc mode aerosols, thus lowering the condensation sink, therefore creating conditions conducive for secondary particle formation18. However, this process appears to be insignificant for the UF mode number concentration at Villum during autumn, although it strongly affects UF mode concentrations during summer (Supplementary Fig. 14).

The decrease in the time spent above the mixed layer in autumn could also be contributing to the increase in the sub-30 nm size range for the Pristine cluster, indicating that a combination of these processes might be driving the increase in the occurrence of the Pristine cluster during autumn. This is, however, unclear as the median size distributions for the lower percentiles of time above the mixed layer deviate from the median size distribution for the Pristine cluster during autumn.

Changing transport patterns or retreating sea ice?

During spring and summer, the UF/total number concentration, and the time air masses spend over the sea is increasing, and the Nucleation cluster occurrence is increasing during summer. A logical hypothesis to the outcome of research objective 3 (investigate which processes are driving this change) would be retreating sea ice and the exposure of air masses to the open ocean, where biogenic emissions are more likely to escape to the atmosphere if transport patterns are similar between years. However, there was no statistically significant trend for the time spent over the sea ice during any season (Fig. 4), therefore, we hypothesize that changing transport patterns are driving the increased time spent over the sea during the summer. To explore this hypothesis, the changes in air mass transport patterns were examined in detail.

During spring, the trend in trajectory frequency for only trajectories corresponding to available SMPS measurements and only trajectory steps over the sea was calculated. The results are displayed in Fig. 6. During spring, air masses below the mixed layer mainly resided over the southern part of the Greenland Ice Sheet and the central Arctic Ocean (Fig. 6a). Freud et al.10 analyzed transport patterns for 240-hour back trajectories, regardless of altitude, arriving at Villum from 2010–2013, and found a similar transport pattern during March and April, although with a greater geographical extent due to the extended trajectory length. For trajectory steps over the sea during spring, the main source regions are the Greenland Sea, the Denmark Strait, and the Davis Strait (Fig. 6b). During the spring, trajectories arriving only when SMPS measurements are available are more frequently arriving from the Greenland and Barents Sea and less frequently from the central Arctic Ocean (Fig. 6c). Thus, air masses are experiencing a shift from a northerly origin to a more southerly one. Trajectory steps over the sea are more frequently arriving from the Greenland and Barents Sea, although with much fewer significant trends due to the scarce occurrence of the open ocean during the spring (Fig. 6d). The shifting transport patterns during spring are supported by a decrease in the frequency of air masses arriving from sea ice-covered regions in the central Arctic Ocean and an increase in the frequency from sea ice-covered regions on the east coast of Greenland and in the Barents Sea, north and west of Svalbard (Supplementary Fig. 15). This indicates that the changing transport patterns are occurring below the mixed layer and not that air masses are spending less time above the mixed layer and more below. The lower and upper confidence intervals for trends in Fig. 6 and Supplementary Fig. 15 are displayed in Supplementary Fig. 16, respectively.

Fig. 6: Trajectory frequencies and trends during spring.
figure 6

Frequencies and trends in air mass transport patterns for air masses that were below the mixed layer during the spring. Frequencies and trends in transport patterns for (a and c) all trajectories occurring when SMPS measurements are available and (b and d) for trajectory steps over the sea. The Theil Sen Slope is indicated by the color bar and statistically significant trends on the 90th% confidence level are outlined in black.

During summer, this method was applied to only trajectories corresponding to available SMPS measurements, only Nucleation cluster trajectories, and only trajectory steps over the sea. The results of this analysis are presented in Fig. 7. During the summer, air masses below the mixed layer mainly arrived from the southern part of the Greenland Ice Sheet, Baffin Bay, the Canadian Archipelago, and the central Arctic Ocean (Fig. 7a). Freud et al.10 observed a similar transport pattern during June and July at Villum from 2010–2013. A similar pattern is observed for Nucleation cluster trajectories that were below the mixed layer during summer (Fig. 7b). Previous analysis of PNSD clustering and their geographical extent yielded similar source regions for a cluster associated with Aitken mode particles10. For trajectory steps over the sea during summer, Baffin Bay and the Greenland Sea are the major source regions, with minor contributions from the Beaufort Sea and Laptev Sea (Fig. 7c).

Fig. 7: Trajectory frequencies and trends during summer.
figure 7

Frequencies and trends in air mass transport patterns for air masses that were below the mixed layer during the summer. Frequencies and trends in transport patterns for (a and d) all trajectories occurring when SMPS measurements are available, (b and e) for only Nucleation cluster trajectories, and (c and f) for trajectory steps over the sea. The Theil Sen Slope is indicated by the color bar and statistically significant trends on the 90th% confidence level are outlined in black.

Trajectories arriving only when SMPS measurements are available and that are below the mixed layer are more frequently arriving from Baffin Bay during the summer, and they are less frequently arriving from the Greenland Ice Sheet and the Central Arctic Ocean (Fig. 7d). Trajectories arriving when the Nucleation cluster was observed are also more frequently arriving from Baffin Bay and less frequently from the Greenland Ice Sheet and the Central Arctic Ocean (Fig. 7e), although with a reduced number of significant trends due to the reduced number of trajectories. Figure 7f shows that trajectory steps over the sea are more frequently occurring at Baffin Bay and the Greenland Sea. Lower and upper confidence intervals for trends in Fig. 7 are displayed in Supplementary Fig. 17. A similar analysis was performed for trajectory steps over the snow and sea ice during summer, with the trajectory frequency and trends displayed in Supplementary Fig. 18. During the summer, the main source region for trajectory steps over the snow is the Greenland Ice Sheet with minor contributions from the Canadian Archipelago (Supplementary Fig. 18). The trend analysis revealed air masses are arriving less frequently from snow-covered areas on the Greenland Ice Sheet and more frequently from the Canadian Archipelago adjacent to Baffin Bay (Supplementary Fig. 18). For trajectory steps over sea ice during summer, the main source regions are the central Arctic Ocean although with contributions from the Davis and Denmark Strait (Supplementary Fig. 18). The trend analysis showed that air masses are arriving less frequently from the central Arctic Ocean (Supplementary Fig. 18). Lower and upper confidence intervals for trends in Supplementary Fig. 18 are displayed in Supplementary Fig. 19.

This analysis shows that, during the summer, trajectories below the mixed layer are more frequently arriving from Baffin Bay for all trajectories when SMPS data are available, and less frequently from the Greenland Ice Sheet and the central Arctic Ocean during the summer. Therefore, we postulate that the increase in time air masses spent below the mixed layer and over the sea is due to changing transport patterns. Nucleation cluster trajectories are also arriving more frequently from Baffin Bay and less frequently from the Greenland Ice Sheet and the central Arctic Ocean. This observation further supports the postulate of changing transport patterns as the reason for the increased time spent over the sea. Trajectory steps over the sea are more frequently occurring in Baffin Bay and the Greenland Sea. Conversely, trajectory steps over the snow on the Greenland Ice Sheet and sea ice in the central Arctic Ocean are less frequently occurring during the summer. While trajectory steps over the sea are occurring more frequently from the Greenland Sea and Baffin Bay, the increase in the frequency of occurrence of Nucleation cluster air masses from Baffin Bay indicates the source region mainly associated with this increase is Baffin Bay. This gives compelling evidence that the cause of the increase in Nucleation cluster occurrence, total and UF mode number concentration, and time spent over the sea, during the investigated period, is primarily due to changing transport patterns to areas of the open sea during the summer.

Trends in trajectory frequency for autumn showed very few significant trends, thus changes in the UF mode number concentration and the Pristine cluster occurrence can be ascribed to changes in other air mass history parameters (as discussed above) and not changing transport patterns.

To further test this hypothesis and investigate other factors that could affect aerosol number concentrations and the occurrence of PNSD clusters on a regional Arctic scale, a trend analysis of daily sea ice concentration, sea ice parameters (sea ice extent, first-year ice (FYI), ice free–marginal ice zone (IF–MIZ), and sea ice change (growth/melt) and teleconnections (North Atlantic Oscillation and Arctic Oscillation) was performed. Trends in daily sea ice concentrations are displayed in Supplementary Fig. 20 with the lower and upper confidence intervals displayed in Supplementary Figs. 21 and 22, respectively. Supplementary Fig. 23 displays trends for the sea ice parameters and teleconnections. For details about these parameters and their data sources, please refer to Supplementary Notes3 and 4, respectively. For sea ice concentrations, there are no significant trends in the Greenland Sea during spring nor in Baffin Bay during summer and none of the other parameters listed above exhibit significant trends on this time scale. While sea ice is declining on long time-scales63,64,65,66 the period analyzed in this study is likely too short to reflect large-scale significant changes in these parameters. The Arctic Oscillation has exhibited a shift towards more negative values in recent years67, which indicates increased meridional transport and supports observations in this study, although no significant trends were observed during this period (2010–2018).

This analysis shows that on short time scales at a single site, changes in aerosol number concentrations and PNSDs are affected by transport and precipitation patterns. During the spring and summer, changing transport patterns, to areas of open ocean, are leading to increased total and UF mode aerosol number concentrations, and the increased occurrence of the Nucleation cluster. Interestingly, the spring and summer seasons show contrasting source regions, with the Greenland Sea for spring and Baffin Bay for summer. Whilst during the autumn, increases in precipitation are likely driving an increase in the occurrence of the Pristine cluster (representing clean, natural Arctic background conditions), and changes in atmospheric dynamics (air masses residing more often below the mixed layer) result in the increase of UF mode particles. The synopsis of this study is that atmospheric dynamics (transport patterns) and sinks (precipitation) influence the aerosol population, however, it remains unclear how the impact of source strengths has changed over this period. The contribution of changing source strengths (concentrations of chlorophyll-a along with individual phytoplankton functional types, emissions of DMS and iodine, and the conversion efficiency of DMS to MSA and sulfuric acid) in Baffin Bay and the Greenland Sea to changes in the aerosol population require further investigation.

The changing total/UF mode aerosol number concentrations will have important implications for aerosol–cloud–radiation interactions during the summer. These changes will alter the short and long-wave radiation budget through increased scattering and absorption (depending on the chemical composition), however, the implications for cloud formation remain to be seen. This will only occur if the aerosols grow to sizes where they can act as CCN, which has been demonstrated at Zeppelin Observatory26 and the central Arctic Ocean25, although not currently for Villum26. Resolving the effect of increased aerosol concentrations on cloud formation will have important implications for the radiation budget during the summer, when the cloud-albedo effect results in negative radiative forcing, whereas for the other seasons it is positive68. Similarly, for the autumn season, if the observed trends persist (increase in the occurrence of the Pristine cluster and accumulated precipitation, decrease in time air masses spend above the mixed layer) then this could have implications for the CCN burden as well as the CS. With the autumn freeze-up expected to occur later in the year69,70, a decrease in the CS (through increased precipitation) and an increase in air masses residing in the mixed layer could extend the period of NPF and growth into the autumn. This will have impacts on cloud formation and the radiative budget, as the autumn is often a CCN-limited regime, therefore increases in CCN active aerosols could increase the downward longwave emissivity71. The increase in the occurrence of the Pristine cluster could also increase the supersaturation (through fewer Acc mode particles), and coupled with the increase in UF mode particles could allow sub-micrometer size aerosols to activate more easily into CCN72. If the observed increase in precipitation continues in the future, there will also be drastic consequences for the Arctic environment. Increased precipitation will affect ice sheet mass balance73,74, surface water density75,76, and river discharge77,78 (which will possibly affect the Atlantic meridional overturning circulation); numerous feedback mechanisms will be affected (water vapor, cloud, and ice-ocean/snow albedo)79; biological and ecological cycles may be altered80; and there will be socio-economic implications for Arctic communities81. Overall, an altered hydrological cycle will have positive and negative feedback on the Arctic climate, although both the sign and magnitude of these changes are yet to be determined82.

When coupling this study with the recent work of Heslin-Res et al.33, it appears the properties of aerosols both in the sub- and super-micron size range are undergoing changes due to changing transport patterns, and while sea ice is changing, this does not seem to be substantially affecting aerosol on these time scales. This work demonstrates how the state of the Arctic is changing fast, as stressed by Schmale et al.22. These changes will likely drive changes in the aerosol burden for the Arctic region as a whole further impacting the climate locally and by teleconnections globally. This calls for a need to analyze the aerosol microphysical properties from a pan-Arctic viewpoint37, and highlights the importance of long-term aerosol in-situ measurements, especially in the Arctic regions.

Methods

Measurement site

Villum Research Station (Villum) is located on the Danish military base, Station Nord, in Northeastern Greenland (81° 36’ N, 16° 40’ W, 24 m above mean sea level). Monitoring of atmospheric constituents took place about 2.5 km southwest of the main facilities of the Station Nord military camp at Flyger’s hut (2010–2014) and in the Air Observatory (2014–2018). The Air Observatory and Flyger’s hut are separated by ~200 m’ distance and they are both upwind >95% of the time from local pollution sources at the military base and any indications of local pollution affecting atmospheric measurements were removed.

Villum is situated in a region with a dry and cold climate with a mean temperature of −13.4 °C. The dominating wind direction is from the southwest, except for the summer which experiences a larger contribution from easterly winds, with an average wind speed of 3.7 m s−1. Villum experiences polar night from October to February, polar twilight from February to April and September to October, and polar day from April to September. Temperatures routinely reach below −25 °C during the winter months (December, January, and February) and only exceed 0 °C during the summer months (June, July, and August). Wind speed exhibits two annual maxima, one in February where the median wind speed is consistently elevated, and another one in the summer months when maximum wind speeds are experienced.

PNSD measurements & k-means clustering analysis

PNSDs were measured at Villum using a scanning mobility particle sizer (SMPS), consisting of a medium Vienna type differential mobility analyzer (DMA) and a condensation particle counter (CPC), in the size range 9–915 nm. A description of the instrument can be found in Wiedensohler et al.83. The amount of data captured for each month during 2010–2018, divided into the meteorological seasons (winter: December, January, February; spring: March, April, May; summer: June, July, August; autumn: September, October, November), is displayed in Supplementary Fig. 1. These data were averaged to hourly time resolution, using the arithmetic mean, which allowed for collocation with other environmental variables. Total, ultrafine (UF, Dp < 100 nm), and accumulation (Acc, Dp > 100 nm) mode number concentrations were calculated by integrating the entire PNSD, below, and above 100 nm, respectively. PNSDs were grouped using the Hartigan–Wong k-means clustering algorithm53 which notably clustered the PNSDs normalized to the vector length. For the number of clusters, only 8 and 9 cluster solutions were considered with the 8-cluster solution being preferred to identify 8 a priori groups, namely: Pristine; Nucleation; Bursting; Nascent; UF Bimodal; Acc Bimodal; Haze and Aged. This approach was in part motivated by the long dataset which prevented the random-access memory (RAM) limited cluster statistics from being carried out on the whole dataset. Calculated for each cluster number, the Silhouette Widths and the Dunn Index define a range of cluster numbers, specified from a lower cluster number that results in under clustering of the data to a higher cluster number that results in over clustering of the data. This range is useful when investigating data collected at new sites, where a strategy is used, of recombining groups of PNSDs for an over-clustered solution until the optimum description is given by the analysis. For an established research station such as Villum, such an approach was deemed unnecessary in this case and only 8 clusters were required with a 9-cluster solution used as a check for hidden detail. This method has been used before to describe different aerosol types at Villum and Svalbard20,21,44,45,84. The occurrence of a cluster for a given month was calculated by summing the total number of hourly occurrences of each cluster in each month and dividing by the total number of hourly measurements for that month multiplied by 100%.

Trend analysis

A trend analysis was performed on total, UF, and Acc mode number concentrations, PNSD cluster occurrence, time spent over different surfaces, air mass history, and trajectory frequencies. The data were segmented into the different meteorological seasons. A threshold of 50% data coverage for each day (for aerosol number concentration and accumulated precipitation) and month (for PNSD cluster and surface footprint types) was applied to ensure proper statistical representation of each temporal aggregation. The Mann–Kendall test was used to determine the presence of a significant trend at the 90th% confidence level85,86 and the Theil Sen slope estimator was used to calculate the magnitude of the trend slope87,88, via the 3PW algorithm from Collaud Coen et al.52. The upper and lower 95th% confidence intervals contain the middle 95% of the pairwise slopes, which are normally distributed89. These tests require the data to be serially independent and homogeneously distributed, however, these measurements experience significant lag-1 autocorrelation and exhibit distinct seasonal patterns. Therefore, to satisfy the requirements of the Mann–Kendall test, the data were prewhitened (to remove the lag-1 autocorrelation) and the seasonal Mann–Kendall test was applied (to address the seasonality present in the data)90. The strong seasonality in the Arctic results in different processes and sources dominating at different times of the year (anthropogenic pollution from long-range transport in winter/spring and natural emissions from regional sources in summer/autumn), therefore investigating trends on a seasonal basis investigates the different processes/sources occurring throughout the year. An alternative approach to satisfying the assumptions of the Mann–Kendall test is averaging to coarser time granularities (monthly, seasonal, or annual averages), although this might not necessarily remove the autocorrelation and will significantly reduce the number of data points, which will decrease the statistical robustness of the results, and gives no seasonal information. Therefore, using the seasonal Mann–Kendall test and prewhitening over coarsely time-resolved data will increase the number of data points (thus increasing the statistical robustness) and determine the season in which these changes are occurring (which, due to the strong seasonality of solar radiation and atmospheric constituents, can have large implications for the climate forcing). The 3PW algorithm utilizes two prewhitening methods prior to testing for statistical significance. The first method from Kulkarni and Storch91 removes the lag-1 autocorrelation; this method has a low rate of false positives but lowers the power of the Mann–Kendall test. The second prewhitening method from Yue and Wang92 removes the trend from the data before prewhitening, this method restores the power of the Mann–Kendall test although increases the rate of false positives. Trends must be statistically significant using both methods to be considered significant in the 3PW algorithm and the significance of the trends is taken as the higher one between the two methods. If a significant trend is present, then the slope is calculated using the variance corrected trend free prewhitening (VCTFPW) method from Wang et al.93, which gives an unbiased estimate of the Theil Sen slope. This method maximizes the advantages of these prewhitening methods and minimizes their disadvantages. It should be noted that the slope and confidence intervals are determined from the Theil Sen slope estimator and the statistical significance is determined from the Mann–Kendall test, therefore trends can be statistically significant while the confidence intervals of the slope can cross zero52.

The occurrence of these clusters and time spent over different surface types is based on their relative contribution to the total amount of observations expressed as a percentage, thus their contributions sum to 100%. Normally any changes in the occurrence of one will result in an equal in magnitude and opposite in sign change in the occurrence of another over the period analyzed. Although due to the prewhitening procedure, a zero-sum of the contribution change is not achieved, meaning an increase in occurrence isn’t necessarily accompanied by a decrease in another occurrence, as would be expected through a linear regression analysis.

When referring to a direction and magnitude of a trend, it is implied that said trend is statistically significant.

Back trajectory analysis

Air mass history was examined by use of the HYSPLIT trajectory model94,95,96. Air mass back-trajectories of 168-hour length were calculated arriving at 50 m above ground level for every hour from 2010 to 2018. The trajectory starting height of 50 m was selected as a compromise between capturing air masses that are representative of our sampling site and avoiding trajectories intercepting the surface, which can produce unrepresentative trajectories97. The trajectory length of 168 h was selected to capture the lifetime of aerosols in the atmosphere and assess the geographical extent of air masses but also avoid the uncertainty associated with extremely long trajectory calculations. Trajectories were calculated based on meteorological files from the NCEP/NCAR Reanalysis Data which has a resolution of 2.5° latitude/longitude98. The HYSPLIT model output included meteorological variables along the trajectory path such as precipitation and mixed layer height. Precipitation was integrated along each trajectory to calculate the amount of accumulated precipitation. The percentage of time the trajectory RH (produced from the HYSPLIT model) was above 95%, was used as a proxy for the time air masses spent within a cloud, following Freud et al.10. Only trajectories corresponding temporally to available SMPS measurements were used in this study.

To analyze trends in the trajectory frequency for a given season, a concentric grid centered around Villum, consisting of 2° × 8° grid cells, was constructed. The trajectory frequency in each grid cell was calculated by summing the number of trajectory steps that were below the mixed layer and intersecting each grid cell and normalizing this by the total number of trajectory steps that were below the mixed layer, multiplied by 100%. This method was applied for each spring, summer, and autumn over the measurement period (2010–2018). The 3PW algorithm was used to calculate these trends, although cells rarely exhibited autocorrelation, therefore a simple trend analysis (no prewhitening, only Mann–Kendall and Theil-Sen Slope on the 90th% confidence level) was performed for each grid cell that exhibited no autocorrelation. This method analyzes only trajectories that are below the mixed layer, therefore, changes in a given grid cell could either be due to a shift in trajectory frequency horizontally (i.e., changing circulation patterns below the mixed layer) or vertically (i.e., air masses residing below the mixed layer more often in a given grid cell). Since autumn was the only season to exhibit a change in time spent above the mixed layer, any changes in the trajectory frequency during spring and summer can be interpreted as changing transport patterns.

Surface type footprint analysis

A surface type footprint analysis was performed for each trajectory, as depicted in Supplementary Fig. 3. For each step along the trajectory, the altitude was compared to the height of the mixed layer. If the trajectory altitude was above this height, that step was classified as being above the mixed layer (AML). If the trajectory altitude was below the mixed layer, then the underlying surface footprint type (land without snow, sea, sea ice, or snow on land) was recorded using a polar stereographic map of the Northern Hemisphere classified into 1024 × 1024 24 km grid cells. The snow and ice coverage values used for the surface footprint type analysis were produced by the National Oceanic and Atmospheric Association/National Environmental Satellite, Data, and Information Service (NOAA/NESDIS) Interactive Multisensor Snow and Ice Mapping System (IMS) developed under the direction of the Interactive Processing Branch (IPB) of the Satellite Services Division (SSD). The time spent over different surfaces is expressed as a percentage of the total trajectory length.

For clarity, when referring to air masses over a surface type throughout the text, it is implied that said air mass is below the mixed layer.