Characterization of size-segregated particles' turbulent flux and deposition velocity by eddy correlation method at an Arctic site

Estimating aerosol depositions on snow and ice surfaces and assessing the aerosol lifecycle in the Arctic region is challenged by the scarce measurement data available for particle surface fluxes. This work aims at assessing the deposition velocity of atmospheric particles at an Arctic site (Ny-Ålesund, Svalbard islands) over snow, during the melting season, and over dry tundra. The measurements were performed using the eddy covariance method from March to August 2021. The measurement system was based on a condensation particle counter (CPC) for ultrafine particle (UFP; < 0.25 µm) fluxes and an optical particle counter (OPC) for evaluating particle size fluxes in the accumulation mode (ACC; 0.25 < d_p < 0.7 µm) and quasi-coarse mode (CRS; 0.8 < d_p < 3 µm). Turbulent fluxes in the ultrafine particle size range were prevalently downward, especially in summertime. In contrast, particle fluxes in the accumulation and quasi-coarse mode were more frequently positive, especially during the colder months, pointing to surface sources of particles from, for example, sea spray, snow sublimation, or local pollution. The overall median deposition velocity ($V_{\mathrm{d}}^{+}$) values were 0.90, 0.62, and 4.42 mm s⁻¹ for UFP, ACC, and CRS, respectively. Deposition velocities were smaller, on average, over the snowpack, with median values of 0.73, 0.42, and 3.50 mm s⁻¹. The observed velocities differ by less than 50 % with respect to the previous literature in analogous environments (i.e. ice/snow) for particles in the size range 0.01–1 µm. At the same time, an agreement with the results of predictive models was found for only a few parameterizations, in particular with Slinn (1982), while large biases were found with other models, especially in the range 0.3–10 µm, of particle diameters. Our observations show a better fit with the models predicting a minimum deposition velocity for small-accumulation-mode particle sizes (0.1–0.3 µm) rather than for larger ones (about 1 µm), which could result from an efficient interception of particles over snow surfaces which are rougher and stickier than the idealized ones. Finally, a polynomial fit was investigated (for the ACC-CRS size range) to describe the deposition velocity observations which properly represents their size dependence and magnitude. Even if this numerical fit is driven purely by the data and not by the underlying chemical–physical processes, it could be very useful for future model parameterizations.

1 Introduction

The Arctic region is experiencing rapid climate change in response to the increase in greenhouse gases, aerosols, and other climate drivers, and it is warming 2 to 3 times faster than the global average (Stjern et al., 2019), as indicated both by observations and climate models (Cowtan and Way, 2014; Hartmann et al., 2013). Although this phenomenon, known as Arctic amplification (Schmale et al., 2021), is mainly driven by changes in anthropogenic greenhouse gases and climate feedbacks, short-lived climate forcers such as methane, tropospheric ozone, and aerosols contribute to the observed environmental change (Arnold et al., 2016; Law et al., 2014; Quinn et al., 2008; Sand et al., 2015; AMAP 2015a, b). The most important factors controlling aerosol climate forcing in the Arctic are the long-range transport of aerosols from mid-latitudes, local sources from both terrestrial and marine emissions, the surface energy budget at high latitudes, the annual cycle of the cryosphere, and atmospheric depositions (Quinn et al., 2008; Willis et al., 2018).

Aerosol depositions on uncontaminated snow and ice can be a key factor affecting melting processes (Skiles et al., 2018) because particles can decrease snow and/or ice albedo directly (as particles can contain light-absorbing materials) or indirectly by affecting ice metamorphism. This phenomenon causes a reinforcing feedback that melts the snow and ice, thereby exposing a darker underground and leading to further surface warming (Di Mauro, 2020; Abbatt et al., 2019). The aerosols' lifetime and hence the radiative forcing are strongly influenced by wet and dry deposition processes (e.g. Johnson et al., 2018). Dry deposition is a complex process that is influenced by the microphysical properties of aerosols and their sources, meteorological conditions, and surface morphological characteristics (Donateo and Contini, 2014; Urgnani et al., 2022). The particle exchange between the atmosphere and the surface is also controlled by frictional drag and terrain-induced flow modification (Giorgi, 1986; Stull, 1988). Knowing the factors controlling dry deposition allows us to estimate the residence time of particles in the atmosphere that governs their transport distance and potential climate effects (Nemitz et al., 2002; Pryor et al., 2008). The accurate quantification of particle deposition rates is a necessary prerequisite for the modelling of aerosol cycles, particle size distribution, long-range transport, and radiative forcing potential (Wang et al., 2011; Liu et al., 2011; Zhou et al., 2012; Browse et al., 2012; Menegoz et al., 2012; Lee et al., 2013; Eckhardt et al., 2015; Qi et al., 2017;). Despite its importance, however, the size dependence of the deposition rates is poorly understood, and therefore, the aerosols' atmospheric lifetime remains uncertain (Farmer et al., 2021; Emerson et al., 2020).

Many aerosols' properties have been investigated in the Arctic with regards to chemical composition (Quinn et al., 2009; Köllner et al., 2021), total number and mass concentrations (Croft et al., 2016), aerosol optical properties (Ferrero et al., 2019), ability to act as cloud condensation nuclei (Bulatovic et al., 2021), and number and size distribution (Lupi et al., 2016; Song et al., 2022). However, relatively few cases of aerosol deposition measurements on snow or iced surfaces exist, especially for those using the direct eddy covariance (EC) method (Farmer et al., 2021). These limitations are mainly due to the logistical challenges of collecting continuous data sets in remote areas (Abbatt et al., 2019) and micrometeorology-based measurement techniques, which require fast and sensitive detectors for EC (Burba, 2022). Dry deposition is typically described by the deposition velocity V_d = $-F/C$, where F is a flux, and C is the concentration of the species of interest. V_d provides a particularly useful metric for comparing the results across sites and for modelling particle removal because it is independent of the ambient concentration (Farmer et al., 2021). Furthermore, the current understanding of Arctic amplification is limited by a lack of robust model representations of the regional Arctic feedback processes, with major challenges in representing the aerosol sources and sinks processes (Schmale et al., 2021). In particular, to reduce model uncertainty, a deeper knowledge of the dry deposition velocity as a function of the particles' size is required. Moreover, an assessment of the current deposition models against an observational data set of aerosol fluxes over the cryosphere is clearly needed (Farmer et al., 2021). Global model skills in simulating the Arctic aerosol behaviour have improved in the last few years (Eckhardt et al., 2015; Arnold et al., 2016); however, significant discrepancies remain in our understanding of Arctic aerosol deposition and removal (Saylor et al., 2019; Emerson et al., 2020; IPCC, 2021). Global models generally make use of an aerosol deposition module with a particle-size-dependent resistance approach developed for specific deposition surfaces (Slinn, 1982; Wesely and Hicks, 2000). The model developed by Zhang et al. (2001) expanded on the Slinn (1982) approach by incorporating simple empirical parameterizations for dry deposition processes. Zhang et al. (2001) also expanded on the application of the resistance approach to also include ice and/or snow surfaces. The deposition rates described by Zhang et al. (2001) have been compared to some observations over vegetated surfaces; however, the parameters used to tune modern deposition models over the cryosphere have not been tested against observations (Khan and Perlinger, 2017).

The principal aim of this work is to measure the particle number fluxes and the related dry deposition velocities for size-segregated particles (from ultrafine to quasi-coarse range) at an Arctic site located in the Svalbard archipelago (Norway). In aiming to characterize the effect of surface properties on dry deposition, we performed continuous observations from the coldest months (on snow surface) to the snow melting period, and all through the early summer (snow-free surface). For these three conditions, a characterization of the deposition velocity as a function of particle diameters will be provided. Finally, a data fit parameterization of the deposition velocity on micrometeorological friction velocity has also been obtained.

2 Methodology

2.1 Measurement site

Aerosol fluxes were measured at the Gruvebadet laboratory located southwest of the village of Ny-Ålesund (78^∘55^′ N, 11^∘56^′ E; Fig. 1) in the archipelago of Svalbard (Norway). The measurement campaign started on 15 March 2021 and lasted until 15 August, for a total of 5 months. The site is characterized by small hills and depressions (height differences below 10 m), whereas the land cover is characterized by dry tundra or bare soil (Magnani et al., 2022) during summer months when the snowpack disappears. Ancillary meteorological measurements of air temperature and relative humidity were collected at 5 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ (above ground level) at the Amundsen–Nobile Climate Change Tower (CCT), which is located about 1 km northwest of the Gruvebadet laboratory (Mazzola et al., 2016).

Figure 1 (a) Location map of the study site at Ny-Ålesund (Svalbard, Norway). Purple and gold points indicate the Amundsen–Nobile Climate Change Tower and the Gruvebadet laboratory, respectively. Height contours (above sea level) are also presented. Figure credit: ©Norwegian Polar Institute (https://www.npolar.no, last access: 5 September 2022). In panel (b), the flux footprint of the EC system is represented (see Sect. 3.2) with the fractions of the total flux originating within the respective contours.

2.2 Instruments

The measurement system was located on the rooftop of the Gruvebadet laboratory at 9.6 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ on a pneumatic mast (Fig. S1a and b in the Supplement). Specifically, the EC station included an ultrasonic anemometer (R3-100; Gill Instruments Limited, UK) with a 100 Hz acquisition frequency, a condensation particle counter (CPC; TSI 3756; TSI Incorporated, USA) measuring the total particle number concentration, and an optical particle counter (OPC; 11-D; GRIMM AEROSOL TECHNIK, Germany) measuring in 16 size channels (from 0.25 to 3 µm). To measure the particle concentration and fluxes from both the CPC and OPC, a maximum acquisition rate of 1 Hz was used. Air was sampled by a unique inlet system for CPC and OPC (Fig. S2). It was placed at a small distance from the sampling volume of the anemometer (about 25 cm). The air sample was driven into the laboratory by a silicon conductive tube that is 6 m long (internal diameter 26 mm), ending with a (three-way) steel flow splitter that is 0.15 m long and 25 mm in diameter, to supply the air sample to both the CPC and the OPC. A nominal flow rate of 60 L min⁻¹ was applied to the flow splitter for the sampling operation, giving a turbulent flow (Reynolds number equal to 4073) and minimizing the temporal distortion of concentration fluctuations. The CPC was connected to the flow splitter by a 0.82 m long silicon conductive tube (6 mm internal diameter) operating at a flow rate of 1.5 L min⁻¹, whereas a portion of 1.2 L min⁻¹ was aspirated by the OPC through a 0.80 m long (4 mm internal diameter) tube of the same type. The particle losses of the sampling system were calculated according to the formulation by Hinds (2012) for turbulent flow in the larger section of the tube (6 m) and according to Kulkarni et al. (2011) for the laminar flow inside the two narrow tubes of the inlet system (beyond the flow splitter). The total particle losses amount to 9 % (on average; mainly diffusional for smaller particles) for CPC and about 2.3 % (on average; mainly gravitational for greater particles) for OPC (Fig. S3a and b). Flux is affected by particle losses in the same way as the particle concentration; consequently, the deposition velocity should not be affected or marginally overestimated by the particle losses being calculated from the ratio of the flux on the concentration. According to TSI Incorporated, the D₅₀ cut-off diameter (at 50 % efficiency) of the TSI 3756 CPC is 2 nm under normal laboratory conditions. The particle penetration curve through the whole sampling system, calculated as the product of the penetration factors in the two tubing sections for CPC (Mordas et al., 2008; Kupc et al., 2013), shows that D₅₀ is about 5 nm (Fig. S3). Therefore, the data analysis for concentration, fluxes, and deposition velocity refers to particles with a diameter starting from 5 nm, which causes large particle losses for smaller diameters. Both particle instruments were running together through a specific software, developed by our research group, which is able to synchronize measured data with the anemometer output. Air temperature and relative humidity were measured at CCT (5 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) by a conventional thermo-hygrometer (Vaisala, Finland; model HMP45AC). The 10 min average physical size distributions above the canopy were measured using a Scanning Mobility Particle Sampler (SMPS; TSI Incorporated; 3034) sampling from the base of the roof (5.5 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$), corrected for particle-size-specific tubing losses, and then used to compute the particle number size distributions for each half-hour in the size range from 10 to 470 nm. More details of the particle size distribution measurements which were conducted at Gruvebadet laboratory are given in Lupi et al. (2016). Precipitation data were measured during the campaign by a laser optical disdrometer, Parsivel² (OTT HydroMet, Germany), installed on the roof of Gruvebadet. The Parsivel² can measure the size and fall speed of the hydrometeors for a comprehensive measurement of all precipitation types (rain, snow, and hail). The Parsivel² disdrometer can measure droplet sizes from 0.25 mm to about 25 mm, with 32 classes of varying diameter intervals. The velocity categories range from 0 to 22.4 m s⁻¹, with 32 classes of varying intervals. Details of the instrument and the measurement technique used to determine the size and velocity of hydrometeors can be found in the literature (e.g. Löffler-Mang and Joss, 2000; Battaglia et al., 2010; Tapiador et al., 2010). An automatic snow meter station was located close to Gruvebadet laboratory at about 100 m southeast during the campaign period. It automatically provides continuous snow data, including near-infrared images of the snow cover area, snow depth, internal snow temperature, and liquid water content profiles at different depths, with a time resolution of 10 min.

2.3 Eddy covariance data analysis

The eddy covariance micrometeorological technique was applied in this work to quantify the atmosphere–surface particle exchange in a size-dependent mode. The flux is computed as the mean cross-product of the fluctuations in a scalar concentration (s^′) and the vertical component of wind velocity (w^′) as F_s = $\stackrel{\mathrm{‾}}{w^{\prime }{s}^{\prime }}$, where the overbar indicates the average values (Kaimal and Finnigan, 1994; Stull, 1988). Turbulent fluxes were calculated on a 30 min basis, using a homemade code developed in MATLAB^® 2018b. In this work, the micrometeorological convention was used, according to which upward fluxes are positive, thus corresponding to emissions, while downward fluxes (towards the surface) are negative. EC measurements allow us to calculate the following main turbulence characteristics: the virtual sensible heat flux, which is $(H=c_p\mathit{\rho }\stackrel{\mathrm{‾}}{w^{\prime }{T}_{\mathrm{s}}^{\prime }}$; hereafter sensible heat flux), where T_s represents the sonic temperature, c_p = 1005 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ is the specific heat at constant pressure, and ρ is the air density; the turbulent kinetic energy, which is $\mathrm{TKE}=\mathrm{1}/\mathrm{2}({\mathit{\sigma }}_{\mathrm{u}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{w}}^{\mathrm{2}}$), where σ_u,v,w is the standard deviation of the wind components; and the friction velocity, which is $u_{\ast }=({\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{4}}$, where $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ and $\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}$ are the horizontal momentum fluxes. Using CPC and OPC measurements together, a complete characterization of the particle number concentration can be obtained from ultrafine to coarse particle size ranges across 17 size bins. Particle fluxes from the particle number concentration were calculated according to $F_{\text{Ni}}=\stackrel{\mathrm{‾}}{w^{\prime }{c}_i^{\prime }}$ for each particle size bin (index i). Size-resolved exchange velocities (V_ex i) of aerosol particles were defined according to the following:

$$\begin{array}{}\text{(1)}& V_{\text{ex}i}=-F_{\text{Ni}}/N_i,\end{array}$$

which shows the turbulent flux of each stage normalized by the respective particle number concentration. The minus sign is used to define the positive values of the deposition velocity ($V_{\mathrm{d}}^{+}$) as a transport toward the surface (deposition) and negative values ($V_{\mathrm{d}}^{-}$) as a transport into the atmosphere (emission), respectively. Positive and negative fluxes will be treated as two separate processes following the uneven distribution of positive and negative fluxes around zero. The separation of upward and downward flux periods or systematic removal periods of positive flux based on other parameters is a common practice in deposition studies (e.g. Nilsson and Rannik, 2001; Vong et al., 2004; Pryor et al., 2013; Lavi et al., 2013; Emerson et al., 2020). Atmospheric stability is a crucial quantity for the characterization of the turbulent pollutant dispersion in the environment, and it is considered in the parameterizations of turbulent characteristics and wind profiles (Nordbo et al., 2013). To classify the different atmospheric conditions over the entire measurement period, the atmospheric stability parameter is defined as ζ = $z/L$, where z is the measurement height, and L represents the Obukhov length, which is computed as follows:

$$\begin{array}{}\text{(2)}& L={\displaystyle \frac{-u_{\ast }^{\mathrm{3}}\stackrel{\mathrm{‾}}{T_{\mathrm{s}}}}{\mathit{\kappa }g\stackrel{\mathrm{‾}}{w^{\prime }{T}_{\mathrm{s}}^{\prime }}}},\end{array}$$

where T_s is the sonic temperature measured by the anemometer. Five stability classes are estimated, according to Nordbo et al. (2013), as follows: very unstable (ζ < −1) and unstable (ζ < −0.01), neutral (−0.01 $\le \mathit{\zeta }\le$ 0.01), stable (ζ > 0.01), and very stable (ζ > 1).

Figure 2 (a) Normalized median cospectra of the kinematic heat flux (blue), ultrafine particle (UFP; red), and accumulation mode and quasi-coarse mode (ACC-CRS) particle flux (green) as a function of f_n. The binned median cospectra (about 1500) were computed from a continuous period of 1 h. The dashed black line represents the slope theoretically predicted in the inertial subrange (Kolmogorov, 1941). (b) Median normalized ogives for kinematic heat flux as a function of natural frequency (Hz) for the three classes of atmospheric stability.

2.4 Data pre-processing

Half-hour periods were hard flagged for dropouts, discontinuities caused by power loss, or values outside the absolute limits and discarded from the data set. These events resulted in 25 % and 24 % of data being rejected for the CPC and OPC fluxes, respectively. The OPC measurements were discontinued on 29 July due to technical issues. The raw data were pre-processed by first applying a despiking procedure to eliminate spikes caused, for example, by electronics issues. Spikes in the 100 Hz (anemometer) and 1 Hz (CPC and OPC) time series were removed from the data set and replaced by linear interpolation, using a procedure described by Vickers and Mahrt (1997). A spike-detection algorithm was applied to the raw high-frequency data defining spikes as absolute deviations from the mean of a threshold value that is 6-fold σ (where σ is the variance of the 10 min sub-interval). Using a closed-path instrument (i.e. CPC or OPC) with a sampling tube, the air sampled by the sonic anemometer takes many seconds to reach the scalar sensor, leading to a time lag between the vertical wind component fluctuations w^′ and concentration fluctuations c^′. Without correcting for the time delay, the vertical wind component fluctuations do not correlate properly with concentration fluctuations, leading to incorrect flux estimation. Such a time lag was estimated by means of a cross-correlation analysis, by moving the time series forward, defining the maximum covariance between the vertical component of wind speed and particle number concentration, and considering a time window between 3 and 5 s. This time window was estimated using the flow rate, tube length, and response time (Deventer et al., 2015). A mean time lag of 4.14 s was calculated for CPC and 3.77 s for OPC (for all size channels), respectively. The wind data were rotated along the streamline reference system (McMillen, 1988) via three rotations. The first two rotations set a reference system that, for each averaging period, aligns the streamline velocity component with the mean wind velocity vector. The third rotation was performed, except when the absolute value of the angle of attack (McMillen, 1988) was greater than 15^∘ (about 3 % of total cases). In neutral or stable atmospheric conditions with low wind speed and weak and intermittent turbulence (Sun et al., 2012; Schiavon et al., 2019), the submeso motions do not follow the surface layer similarity, but they can still contribute to observed statistics (Vickers and Mahrt, 2006; Liang et al., 2014; Mortarini et al., 2016). According to some authors (e.g. Vickers and Mahrt, 2003), similarity relationships should be evaluated only after filtering out the contribution from these motions. In this work, the energy contributions related to non-turbulent submeso motions, with timescales often longer than the investigated time window, were removed by a recursive digital filter both for heat and particles fluxes (Falocchi et al., 2018; Pappaccogli et al., 2022). The recursive digital filter worked on a different timescale, depending on the atmospheric stability conditions. Ogive analysis (Fig. 2b) was carried out in order to estimate a proper timescale (Metzger and Holmes, 2008). A detailed description of the methodology used for spectral and ogive analysis is reported in Pappaccogli et al. (2022). The timescale for unstable atmosphere is 522 s, whereas the value decreases to 350 and 340 s for neutral and stable conditions, respectively. This filter does not introduce any phase shift or signal amplitude attenuation in the filtered time series. The filtering procedure used two data buffers (1800 s long) before and after the considered 30 min period of investigation. A fundamental assumption of the EC method is that fluctuations are statistically stationary during the chosen averaging time to ensure the calculation of an ensemble average. Stationarity tests, as reported by Mahrt (1998), were carried out on sonic temperature and the particles' concentration fluctuations (Cava et al., 2014; Věcenaj and De Wekker, 2015). A lower detection limit for the fluxes in the sampling system was computed using the method proposed by Langford et al. (2015) and defined at 2.8 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the CPC and 0.3 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ for the OPC. To ensure the study of particle dry deposition only, all data corresponding to a precipitation intensity greater than 0.1 mm h⁻¹ for a period greater than 5 min (on the averaging period of 30 min) were also rejected. Errors associated with the random and limited statistical counting (relative error; %) were estimated through the approach reported in Deventer et al. (2015) for the particle number concentration δ(N) and fluxes $\mathit{\delta }\left(w^{\prime }{N}^{\prime }\right)$. The method reported in Fairall (1984) was used for the deposition velocity δ(V_d) for each size range (Table 1). Obviously, uncertainties due to discrete counting δ(N) are negligibly small for all stages, as they increase from ultrafine to coarse particles (greater than 100 %). Hence, the relative flux uncertainty $\mathit{\delta }\left(w^{\prime }{N}^{\prime }\right)$, due to limited counting statistics, is also moderate (on average 9 %), and the average exceeds 25 % only for the 11th OPC channel. If the counting errors in the deposition velocity δ(V_d) are considered, then on the first size channel (CPC) it was very low (< 5 %). The same error for the first 11 channels of OPC (0.25–0.80 µm) was on average about 85 %, while for the remaining channels (1–3 µm) it was on average greater than 100 % (on average 109 %). To lower the associated statistical counting error, especially on deposition velocity, the first nine channels have been pooled together, as have the remaining seven channels (Whitehead et al., 2012; Conte et al., 2018; Donateo et al., 2019).

Table 1 The relative error associated with the random and limited statistical counting for the particle concentration, fluxes, and deposition velocity is shown. Aerodynamic (D₅₀) cut points and the respective geometric mean diameter (D_gm) are also reported. Note that UFP is for ultrafine particle, ACC is for accumulation mode, and CRS is for quasi-coarse mode.

The previous considerations confirm that it is possible to aggregate the particles into the three size ranges mentioned below, namely ultrafine (UFP; 5 nm < d_p < 0.25 µm), accumulation (ACC; 0.25 < d_p < 0.7 µm), and quasi-coarse (CRS; 0.8 < d_p < 3 µm) mode, with the last mode indicating a particle size range between large accumulation mode and small coarse particles. It is worth noting that UFP particle concentration has been obtained as the difference between the total number concentration (CPC measurement) and the OPC total concentration in the size range 0.25–1 µm. The relative counting errors in V_d for these two groups are approximately 29 % for ACC and 41 % for CRS (Table 2). This size aggregation is also supported by the correlation analysis of the concentration time series associated with the different size classes. The resulting correlation shows that all classes in the first and second group have a good temporal correlation (Pearson coefficient > 0.6) with each other. The UFP fraction was obtained from the difference between the sum of the first 12 OPC channels (up to 1 µm) and the total number concentration of the CPC. This subdivision will be useful for studying the characteristics and trends of the particles and their possible correlation with meteorological and micrometeorological parameters.

Table 2 The relative error associated with the random and limited statistical counting for the particle concentration, fluxes, and deposition velocity is shown. Aerodynamic (D₅₀) cut points and the respective geometric mean diameter (D_gm) are also reported.

2.5 Spectral analysis and corrections

The cospectra were computed on the vertical wind speed component and the analysed scalars (sonic temperature, ultrafine, accumulation, and coarse particle concentration) by means of the fast Fourier transform (FFT). Figure 2a shows normalized median cospectra (about 1500) of the kinematic heat flux, UFP, and ACC-CRS particle number flux as a function of a normalized frequency $f_{\mathrm{n}}=fz/u$, where u is the mean wind velocity, z is the measurement height, and f is the sampling frequency. The cospectra were calculated over 1 h periods. ACC-CRS cospectra were calculated based on the sum of all OPC size channels. The cospectra ${\text{Co}}_{w^{\prime }{x}^{\prime }}$ between the vertical wind velocity and the scalar x are normalized with the correlation $\stackrel{\mathrm{‾}}{w^{\prime }{x}^{\prime }}$. Cospectra calculated for all variables display similar patterns, showing a $-\mathrm{7}/\mathrm{3}$ decay for frequencies above 10⁻¹, which is within the universal equilibrium range.

Specifically, the kinematic heat flux cospectrum shows a slight slope at the high frequencies, which tend to level off in the low-frequency range. Results indicate that the contribution to the correlation between w and particle concentration follow a kinematic heat flux trend up to normalized frequencies in the interval 0.5–0.8. Therefore, the system used was generally able to measure concentration fluctuations at frequencies that made a substantial contribution to the vertical turbulent fluxes of particles. Figure 2b shows the normalized median ogives of kinematic heat flux. A comparison of the ogives for the different atmospheric stability classes (defined below in Sect. 3.1) shows that the higher frequencies contribution is larger under stable conditions. Indeed, no relevant contribution to kinematic heat flux occurs above 1 Hz (1 s timescale) under unstable conditions, whereas a small turbulence contribution is still present up to about 2 Hz (0.5 s) and 5 Hz (0.2 s) under neutral and stable conditions, respectively. Ogives obtained for neutral conditions display some irregularities due to the low number of available cases. In Fig. S4, the ogives for the CPC and OPC data for the whole data set (with no distinction for stability conditions) were shown. Spectral corrections were applied to the UFP, ACC, and CRS particle number fluxes, using the theoretical approach described in Aubinet et al. (2012). The first-order time constant of the CPC and OPC measurement systems was determined by estimating the time response (at the first order) to a concentration step in a laboratory experiment. The results were τ_CPC = 0.77 ± 0.01 s and τ_OPC = 0.40 ± 0.03 s (identical for each size channel). High frequency under the sampling of the EC system, due to the relative low-frequency response of the CPC and OPC, leads to an error (underestimation) of the measurement of the particles' fluxes. It means that the full atmospheric cospectra between the vertical wind velocity and the particle concentration can be under sampled at high frequency (Fig. 2a). High-frequency losses were corrected following the parametric/in situ approach developed by Horst (1997), and these losses have been quantified on average to 30 % for CPC and 21 % for OPC. A spectral correction of the attenuation of the fluctuations during the transport of the aerosol sample in the inlet tube was obtained according to the method reported by Massman and Ibrom (2008). This spectral correction for our measurement system resulted in a very low associated error (lower than 0.01 %).

Figure 3 Time series of the principal meteorological variables measured during the campaign in Ny-Ålesund. (a) Air temperature (^∘C), (b) relative humidity (%), (c) wind velocity (m s⁻¹), (d) snow depth (m), (e) daily cumulative precipitation (mm), and (f) wind direction (^∘ N). In panel (d), three colour bands were used to separate the period characterized by snow cover (dark cyan) from the melting snow period (cyan) and, finally, the snow-free phase (light cyan).

3 Results and discussion

3.1 Site meteorology and micrometeorology

During March and April, the air temperature was on average −7 ^∘C, while during summertime (June–August) the temperature was above 5 ^∘C (Fig. 3a).

Positive temperature values occur approximately at the end of May (29 May), leading to a gradual melting of the snowpack, which reached a maximum height of 1.20 m on 21 May (Fig. 3d). Specifically, the snowpack had an average height of 1 m until 30 May (hereafter named the snow cover period or SC). Snow depth decreased gradually for about 8 d due to its compacting processes. Snowpack was on average 0.54 m from 30 May for about a month until 3 July. In this period (hereafter named the snow melting period or SM), snowpack went through a complete melting phase until the summer period, with the presence of widespread ponds and snow patches in the footprint of the EC system. This period was characterized by positive daily average temperatures. The last period of the measurement campaign, from 3 July to 15 August, will be referred to as the snow-free period (SF), with the surface covered by dry tundra. Relative humidity in general increases from a minimum of 37 % during a severe storm in the month of March up to 82 % in July (Fig. 3b).

A direct correlation is observed between snow melting and increased relative humidity, even if the latter variable is also strongly influenced by the presence of the Kongsfjorden close to the site. Two prevailing wind directions can be distinguished, with one characterized by winds coming from the northwesterly sector for 27 % of the cases and one with winds coming from the southeast for 42 % of the cases (Fig. 3f). Indeed, the prevailing wind directions are essentially along the direction of the Kongsfjorden. On the other hand, from mid-to-end June onwards, the wind field shows two components from west–southwest and northeast (with low winds). The highest wind speeds are measured in the northwesterly direction, with an average of 4.9 m s⁻¹, while in the southeasterly direction the average wind speed is about 4.0 m s⁻¹ (Fig. 3c). In the measurement period, 22 rain days were observed starting from 1 June, with a cumulative precipitation above 2 mm for 9 d. In this period, the maximum amount of cumulative rain on daily basis was of 23 mm on 4 June. In the first part of the campaign, until 30 May, 16 d of precipitation were observed, with a peak of snow height (1.20 m) on 20 May.

Figure 4 A time series of the principal micrometeorological variables measured during the campaign carried in Ny-Ålesund. (a) Sensible heat flux, (b) momentum turbulent flux, and (c) friction velocity u^∗. In each panel, three colour bands were used to separate the period characterized by snow cover (dark cyan) from the melting snow period (cyan) and, finally, the snow-free phase (light cyan).

Figure 5 Percentage of atmospheric stability over the daily cycle for the (a) snow cover, (b) melting, and (c) snow-free periods. Stability classes based on threshold values of the stability parameter ζ were defined as very unstable (dark red; ζ < −1) and unstable (red; ζ < −0.01), neutral (green; −0.01 $\le \mathit{\zeta }\le$ 0.01), stable (blue; ζ > 0.01), and very stable (dark blue; ζ > 1).

The results from this work are completely in line with those reported in (Mazzola et al., 2016) for the data collected at CCT from 2010 to 2016. In the measurement period, the sensible heat flux was on average negative (−0.8 W m⁻²; Fig. 4a). The results show the presence of a long period with negative energy fluxes (with an average value of −10.7 W m⁻²) and a minimum around −121 W m⁻² both while snow cover is present and during the snow melting phase, when the atmosphere is warmer than the surface. As the snow melts, positive sensible heat fluxes begun to appear, with values of up to 191 W m⁻² (on a half-hourly base) in the month of July. This behaviour has also been observed in other measurement campaigns in the Antarctica (Van den Broeke, 2005; Van As et al., 2005) and in the Arctic (Kral et al., 2014), where the snowpack acts as a sink of heat during the winter and spring months. The kinematic momentum flux was typically downward (Fig. 4b), with an average of −0.034 m² s⁻². However, for about 16 % of the data, the direction of the momentum flux showed positive values. A particular event can be observed (Fig. S5) on 2 June from 14:00 to 18:30 LT, with large positive values of momentum flux (on average 0.40 m² s⁻²). This event was characterized by air masses coming from the glacier Brøggerbreen, which is located west (250–270^∘) of our measurement site, with relatively high velocity (on average 8 m s⁻¹). More frequent positive momentum flux events with lower intensity occurred, coinciding with air masses coming from east–southeast (from the Kongsfjorden), most likely due to the effect of a very shallow sea breeze swell (Grachev and Fairall, 2001). Figure 4c reports the time series of friction velocity, with a mean value of 0.17 m s⁻¹ over the whole measurement period. No specific differences can be noted in the friction velocity behaviour due to the changing in snowpack characteristics. As it can be noted in Fig. 5, stable conditions prevailed during the snow cover period (Fig. 5a), with 62 % of cases (and 26 % very stable), while unstable stratification was more frequent in the summer months (Fig. 5c; 45 % and 43 % for unstable and very unstable, respectively). Neutral stratification occurred with a maximum frequency of 7 % during the melting period and with no diurnal variability. It is worth noting the difference in a typical daily pattern at the mid-latitudes, characterized by a much greater contrast between day and night planetary boundary layer (PBL) temperature profile structures, compared to more unstable conditions during the day and more frequent stable conditions at night.

3.2 Footprint analysis

The ratio between wind velocity and friction velocity in atmospheric neutral conditions was used to evaluate the average roughness length z₀ for the measurement site, using a parameterization based on similarity theory (Toda and Sugita, 2003). The results gave z₀ = 0.005 m if the whole measurement period is considered. After separating the period with snow coverage from that without snow, the calculated values were z₀ = 0.002 m and z₀ = 0.004 m, respectively. No significant statistical differences were found when exploring the roughness length for different wind sectors. Besides, an evaluation, using the method reported in Toda and Sugita (2003), gave a null displacement height when considering again all wind directions together. The source area for the scalar fluxes have been evaluated using a Lagrangian footprint model by Kljun et al. (2015). The results of the flux footprint analysis of the EC system are shown in Fig. 1b, with fractions of the total flux originating within the respective contours (until 80 %). The particle fluxes measured at the rooftop site represented a surface area of about 0.48 km² (with respect to the 80 % contour line), with a maximum distance of 400 m both in the northwesterly and southeasterly directions (otherwise 300 m in northeasterly and southwesterly sector). The flux peak contribution was in the wind direction sectors at about 40 m (± 20 m; Fig. 1). However, the source areas were quite homogeneous for all wind direction sectors around the measurement site, with 100 % of snow coverage for the first period of the campaign; from then onwards, the land cover was largely dry tundra and bedrocks, with about 1 % water surfaces (an arctic lake) in the flux footprint (Fig. 1).

Table 3 Statistical quantities of UFP, ACC, and CRS for the particle number concentration and their turbulent fluxes for the different phases considered in this work (snow cover, snow melting, and snow-free period). Note that SD is for the standard deviation, and SE is for the standard error.

Figure 6 Upper panels show the monthly box plots for (a) N_UFP, (b) N_ACC, and (c) N_CRS. Bottom panels show the monthly box plots for (d) F_UFP, (e) F_ACC, and (f) F_CRS. Continuous red lines represent median values. value. Boxes represent the 25th and 75th percentiles. Whiskers correspond to ± 2.7σ and 99.3 % data coverage. In black (square) are the average values for the measured variables.

3.3 Particle concentration and turbulent fluxes

Average particle number concentration, over the whole measurement period, for the UFP (hereafter N_UFP) was 595 cm⁻³ (median 238 cm⁻³) and ACC mode (hereafter N_ACC) was 25 cm⁻³ (median 19 cm⁻³) and for the CRS mode (hereafter N_CRS) was 0.47 cm⁻³ (median 0.30 cm⁻³; Table 3). N_UFP increases from March (on average 166 cm⁻³) to August 2021, reaching its maximum in July (951 cm⁻³; Fig. 6a). On the other hand, N_ACC and N_CRS show a different, anti-correlated behaviour with respect to the N_UFP (Fig. 6b and c). N_ACC and N_CRS reach the maximum concentration in March with, respectively, 40 and 0.67 cm⁻³, and then they decrease down to a minimum concentration in June (18 cm⁻³ for N_ACC and 0.47 cm⁻³ for N_CRS). The same behaviour was observed by Croft et al. (2016) at the Zeppelin Observatory on the Zeppelin Mountain. The complete annual cycle is remarkably similar to that observed at Mount Zeppelin over an earlier 10-year period from 2000 to 2010 (Tunved et al., 2013). The particle size distribution was continuously monitored, combining the SMPS and OPC measurements, from the nanoparticle (4 nm) up to the micro-particle (2.5 µm) size range. The total size distribution plot (Fig. 7; black markers) refers to the whole measurement period. This plot shows a two-mode size distribution, with a high concentration of particles in the 0.02–0.04 µm range, with a peak at 0.03 µm (ultrafine mode), and 0.1–0.2 µm range, with a peak at 0.15 µm (accumulation mode). Considering the snow cover period, the particle size distribution (red markers) also presented two modes, in which one peaked at 0.03 µm and the other at 0.15 µm; with respect to the whole period of analysis, the dominant mode was the accumulation one.

Figure 7 Particle size distribution for the measurement period.

Proceeding towards the melting period and throughout the snow-free period, the two-mode size distribution tends to peak in the ultrafine range (0.02–0.03 µm), with an increasing concentration. In particular, during the snow melting, the particle concentration in the accumulation mode decreased, while the UFP mode significantly increases in this period. Previous studies at the mountain station of Zeppelin Observatory (south of Ny-Ålesund) showed that new particle formation events seem to be a rather common phenomenon during the summer season, and this is the result of both the photochemical production of nucleating or condensing species and a low condensation sink (Tunved et al., 2013; Ferrero et al., 2016). The particle size distribution pattern from the snow cover to the snow-free period can be explained as follows: during March–May, there is sunlight and also a large deposition aerosol surface due to the presence of long-range transported aerosol (Giardi et al., 2016). On the other hand, during June–August, atmospheric conditions are quite different, with high solar insolation and more daylight hours bringing new particle production and supporting their growth to a large size (accumulation mode). Furthermore, the seasonal effect of nucleation precursors is also related to sources because marine biota (the phytoplankton blooms) grow in late spring and early summer and emit methane sulfonic acid (MSA), which gives rise to nucleation (Beck et al., 2021).

The turbulent fluxes of ultrafine (F_UFP), accumulation (F_ACC), and quasi-coarse (F_CRS) particle statistics are listed in Table 3 and consider the snow cover, snow melting, and snow-free periods. The average F_UFP was −16.99 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (median −8.20 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) over the entire measurement period; specifically, the measurement site behaved, on average, like a deposition area (negative fluxes) for particles in that size range, both with and without snowpack (Table 3). The average value of F_ACC for the whole period was 0.64 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (median −0.15 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), but from snow melting onwards, the fluxes became negative in median values, showing a sink behaviour for the measurement site for this size range of particles. Finally, the mean of F_CRS was 0.07 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (median 0.03 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) for the whole period. Overall, in median values, quasi-coarse mode fluxes are positive for all periods, even if, on average, these fluxes were very small. However, this net deposition is not a steady feature. In fact, total and size-resolved turbulent fluxes are highly variable both in magnitude and direction (emission and deposition). Turbulent fluxes in the ultrafine mode are negative for 63 % of the quality-assured cases (Fig. 6d), establishing a deposition behaviour in this size range for the measurement site. Similarly, F_ACC is negative in 52 % of validated cases, while F_CRS is negative for 46 % of available good flux cases. The values of F_UFP have a positive median (0.59 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) at the beginning of the measurement period when ultrafine particle concentration is small and then turn significantly negative (median −17.3 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) in the summer when concentrations sharply increase (Fig. 6a). If the high summertime concentrations of UFP are due to new-particle-formation processes, then these do not occur near the surface and propagate upward but rather take place in the full mixing layer, topping at ≤ 500 m above the ground (Ferrero et al., 2016), with UFP mixing down and being deposited onto the surface. F_ACC and F_CRS generally have a very similar monthly behaviour (Fig. 6e and f), with a maximum value in the month of April (median 0.50 and 0.09 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively) and a minimum in the month of May (median −0.35 and −0.02 ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively). Although remote Arctic environments are normally considered to be receptors for pollution transported from the mid-latitudes, and the seasonal cycle of atmospheric aerosol concentrations witnesses it, at the surface level, small deposition fluxes can overlap with – and being exceeded by – small emission fluxes of the particles. In Ny-Ålesund, the emission fluxes of natural aerosols can originate from soil erosion (in the summer) and snow resuspension and sublimation; besides, a contribution from the village itself and from shipping (Eckhardt et al., 2013) has been occasionally identified but not quantified from our measurements. Finally, it is unclear whether a residual sea spray emission can lead to detectable emission fluxes of sea salt particles at ∼ 1 km inland from the coast. Our measurements – relying on fast particle counting – do not provide direct insights to the chemical nature, properties, and actual origin of the aerosols associated with either emissions or depositions.

3.4 Particle deposition velocity V_d

Exchange velocity (V_ex) calculated from normalization of the turbulent fluxes (Eq. 1) on the number concentration, in general, as seen in the Sect. 3.3, presents both positive and negative values, representing the emission ($V_{\mathrm{d}}^{-}$) and deposition events ($V_{\mathrm{d}}^{+}$). V_ex calculated for the measurement campaign resulted in median values of 0.34 mm s⁻¹ (IQR₁₀₋₉₀ −2.00 to 2.81 mm s⁻¹) for the ultrafine size range, 0.10 mm s⁻¹ (IQR₁₀₋₉₀ −2.02 to 1.53 mm s⁻¹) for the accumulation mode, and −1.36 mm s⁻¹ (IQR₁₀₋₉₀ −11.25 to 8.48 mm s⁻¹) for the coarse range (Table 4). IQR₁₀₋₉₀ represents the interquartile range from the 10th and 90th percentile. Thus, the measurement site acts, on average, as a deposition area for the ultrafine and accumulation particles, while it was a source area for a greater size range (quasi-coarse mode). To explore the behaviour of the deposition events on the measurement site, only the positive values of the deposition velocity will be taken in consideration by studying $V_{\mathrm{d}}^{+}$. The median $V_{\mathrm{d}}^{+}$ value for the UFP range was 0.90 mm s⁻¹ (IQR₁₀₋₉₀ 0.18 to 3.89 mm s⁻¹); in the ACC mode, the median $V_{\mathrm{d}}^{+}$ was 0.62 mm s⁻¹ (IQR₁₀₋₉₀ 0.16 to 2.38 mm s⁻¹), and in the CRS mode, $V_{\mathrm{d}}^{+}$ was 4.42 mm s⁻¹ (IQR₁₀₋₉₀ 1.57 to 13.09 mm s⁻¹). Also, in this case, V_ex and $V_{\mathrm{d}}^{+}$ analyses were separated into three periods (with or without snow and melting). In particular, $V_{\mathrm{d}}^{+}$ increases from the snow surface to the snow-free conditions in all size ranges. It increased from a median value of 0.73 mm s⁻¹ (SC) to 1.14 mm s⁻¹ (SF) for the UFP particles. In the same way, $V_{\mathrm{d}}^{+}$ increased from the SC to the SF phase, with a maximum during SM period for the ACC and CRS size ranges (Table 5).

Table 4 Statistical parameters for exchange (V_ex), deposition ($V_{\mathrm{d}}^{+}$), emission velocity ($V_{\mathrm{d}}^{-}$), and deposition velocity normalized for the friction velocity ($V_{\mathrm{d}}^{+}/u_{\ast }$) calculated on the whole measurement period, separating the UFP, ACC, and CRS particle size ranges.

Table 5 Statistical parameters for the exchange (V_ex), deposition ($V_{\mathrm{d}}^{+}$), and emission velocity ($V_{\mathrm{d}}^{-}$) and deposition velocity normalized for the friction velocity ($V_{\mathrm{d}}^{+}/u_{\ast }$) calculated on snow (SC), snow melting (SM), and snow-free (SF) periods, separating the UFP, ACC, and CRS particle size ranges.

Figure 8 Median $V_{\mathrm{d}}^{+}$ as a function of the geometric mean diameter measured during (a) snow cover, (b) melting period, and (c) snow-free phases. A comparison with the model predictions is also reported in each panel. Observations are shown in symbols and models in lines. Error bars and shaded areas represent the interquartile range. In panel (a), the median values measured $V_{\mathrm{d}}^{+}$ are compared to previous measurements of the deposition velocities over snow. (d) Functional fit for $V_{\mathrm{d}}^{+}$, depending on the geometric mean particle diameter for the different period analysed in this work.

The observed deposition or exchange velocity is in good agreement with previous measurements on snow and ice in the polar regions; however, a comparison with the measured deposition velocities in the literature is a challenge because work on ice and/or snow surfaces is very seldom reported, and often different instruments and methods (such as passive sampling, gradient method, and eddy covariance) are used. Measurements reported by Duann et al. (1988) on snow gave an average $V_{\mathrm{d}}^{+}$ = 0.34 ± 0.14 mm s⁻¹ for particles in the size range between 0.15 and 0.5 µm (Fig. 8a). In the High Arctic, Nilsson and Rannik (2001) report a mean deposition velocity $V_{\mathrm{d}}^{+}$ = 1.4 mm s⁻¹ over ice in the nucleation mode and $V_{\mathrm{d}}^{+}$ = 0.51 mm s⁻¹ in the Aitken mode (with a size distribution peaking at about 60–70 nm).

Gronlund et al. (2002) reported a median $V_{\mathrm{d}}^{+}$ of about 3.3 mm s⁻¹ (and an interval between 0.8 and 18.9 mm s⁻¹) over snow for total particles larger than 10 nm. The values reported by this last work are significantly larger than the values observed in the present data set. Held et al. (2011a) observed median $V_{\mathrm{d}}^{+}$ values ranging from 0.27 to 0.68 mm s⁻¹ during deposition-dominated periods. In Held et al. (2011b), a deposition velocity on the snow surface in the Arctic pack ice ranged from 0.28 to 0.58 (± 0.4) mm s⁻¹ by eddy covariance. There was a fair agreement with the observations from Contini et al. (2010) at the Nansen Ice Sheet (Antarctica) over ice and/or snow, where median V_ex = 0.19 mm s⁻¹ was observed for total particle larger than 10 nm. Furthermore, Contini et al. (2010) observed a median deposition velocity of $V_{\mathrm{d}}^{+}$ = 0.65 mm s⁻¹. In the Hudson Bay (Canada) area, Whitehead et al. (2012) measured a mean V_ex of 0.12 ± 0.11 mm s⁻¹ on sea ice for particles greater than 2.5 nm. Seasonal differences in the deposition have also been observed in the cryosphere, with greater values during warmer months (Macdonald et al., 2017) when collecting snow samples on average every 4 d at Alert, Nunavut (Arctic Canada). In Ibrahim et al. (1983), a deposition velocity between 0.39 and 0.96 mm s⁻¹ over snow surface was observed for 0.7 µm diameter particles using artificial collectors for the liquid scintillation counting technique. Even if the particle sizes range were not exactly the same, the V_ex (1 mm s⁻¹) values reported by Tschiersch et al. (1991) are still comparable with values measured at the closest size range (Fig. 8a). For a more detailed and comprehensive review of the deposition velocities, see Pryor et al. (2008), Whitehead et al. (2012), Saylor et al. (2019), Emerson et al. (2020), and Farmer et al. (2021). From the graphs in Fig. 8, $V_{\mathrm{d}}^{+}$ (median values) increases with particle diameter, and $V_{\mathrm{d}}^{+}$ is greater in all size bins during the melting phase with respect to the other periods. This could probably be due to the increase in the roughness length z₀ as the snow melts (and the consequent shallow pond formation), which leads to an increased capture of the particles (in all size ranges) by the interception and impaction mechanisms. $V_{\mathrm{d}}^{+}$ has the lowest results for the snow surface coverage. To put our observations into context, deposition velocity observations from previous studies have been plotted in Fig. 8, alongside four models for dry deposition velocity developed by Slinn (1982; hereafter S82), Zhang et al. (2001; hereafter Z01), Pleim and Ran (2011; hereafter PR11), and Emerson et al. (2020; hereafter E20). Dry deposition models were based on a custom code (DepoBoxTool) translated in MATLAB (version 2018b) from the original form in Python (https://doi.org/10.5281/zenodo.4749548; Shu et al., 2021; Shu, 2021). Figure 8a shows the deposition velocity measured and predicted during the snow period, whereas Fig. 8b and c display the results during the melting and snow-free periods, respectively. Furthermore, it seems that the $V_{\mathrm{d}}^{+}$ values observed in this work follow the shape of the S82 model output. As reported by Saylor et al. (2019) for smooth surfaces such as water or snow and/or ice, models provide a minimum deposition rate for particles in the accumulation mode. Data suggest a minimum in the 0.08–0.15 µm range, with a strong increase in $V_{\mathrm{d}}^{+}$ above 0.25 µm, rather than the gradual decrease predicted by some of the algorithms (i.e. Z01 and PR11). Our results are in line with the observations of E20, as they revisit the importance of interception processes with respect to Brownian diffusion. Apparently, interception can be a significant term of particle loss, even over ice and/or snow surfaces, whose actual structures and fine irregularities contrast with the idealized smoothness traditionally assumed for this land type.

Figure 9 (a) Relation between median $V_{\mathrm{d}}^{+}$ and friction velocity for the different size ranges. Vertical bars represent the interquartile range of the deposition velocity within the specific interval of u_∗. Friction velocity intervals were selected to optimize the number of data points within each interval and, hence, provide a statistically reliable median deposition velocity. (b) Functional fit for V_n as a function of the geometric mean particle diameter (excluding UFP) for the different periods analysed in this work.

The selected deposition velocity models take as input our meteorological and micrometeorological observations (temperature, pressure, relative humidity, L, and u_∗), site-specific surface properties (roughness length, canopy height, displacement height, and land use category), and particle properties (density and diameter). Note that most model outputs are highly sensitive to canopy structure parameters and friction velocity. It is worth considering, for example, the high sensitivity of PR11 to the value of the convective velocity scale, which can vary significantly with the surface roughness defined for each land use type (Saylor et al., 2019). On the other hand, Z01 produces the widest differences in the predicted deposition velocities among land use types, showing weak deposition velocities compared to other models, especially in the accumulation and coarse mode. The comparison of predictive models with our observations shows a good agreement, especially with S82, which seems to better fit the data over a broader range of particle diameters and properly predict a minimum in the range 0.1–0.3 µm in addition to the deposition velocities for accumulation and coarse mode. On the contrary, the Z01 and PR11 models show the minimum as being around 1 µm, providing a broad underestimate of the deposition velocity for particles above 0.3 µm. Specifically, during the snow period (Fig. 8a), the considered models predict deposition rates for particles below 0.1 µm well, except for Z01, which overestimates $V_{\mathrm{d}}^{+}$ in that range. In contrast, a slight and broad underestimation of $V_{\mathrm{d}}^{+}$ above 0.3 µm was observed for the E20, Z01, and PR11 parameterizations. During the melting period (Fig. 8b), a good agreement was observed in the 0.03–0.3 µm size range, although slight underprediction was observed for S82 and E20. For particle diameters above 0.3 µm, a broad underestimation occurs for all parameterizations, except S82, which shows a good agreement, especially in the particle range above 1 µm. Finally, during the snow-free period, an overestimation of the deposition velocity for ultrafine particles was observed in Z01 and PR11, which remained in agreement with measurements up to 0.3 µm particles. Although the S82 and E20 models show an underestimation of $V_{\mathrm{d}}^{+}$ for ultrafine particles, a good agreement was observed for accumulation and coarse modes, especially for S82. In view of these results, most analysed parameterizations show an underestimation of the deposition velocity that is greater than an order of magnitude and in the range of 0.2–0.3 µm, except for S82. In general, due to the sensitivity of the models to input parameters, the differences between model outputs and our data should be considered with caution. However, our data clearly show that the minimum deposition velocity is reached at around 0.1 µm particle diameter so that $V_{\mathrm{d}}^{+}$ values monotonously increase over the particle size interval from 0.1 to 10 µm, and our data support the parameterizations in the high range of $V_{\mathrm{d}}^{+}$ for this particle size interval.

3.5 Particle deposition velocity parameterization

Regarding the trend of the deposition velocity as a function of particle diameter, $V_{\mathrm{d}}^{+}$ can be described by a third-order polynomial function of the mean geometric diameter in a typical saddle-shaped graph, with an increase in magnitude toward the lowest and highest side of the size spectrum and a minimum in the Aitken mode. A robust (bi-square) polynomial fit to our observations (excluding UFP) was computed for the three considered periods, with a resulting goodness of fit R² and a root mean square error (RMSE) as reported in Table 6.

$$\begin{array}{}\text{(3)}& V_{\mathrm{d}}^{+}\left(D_{\text{gm}}\right)=p_{\mathrm{1}}(D_{\text{gm}}{)}^{\mathrm{3}}+p_{\mathrm{2}}(D_{\text{gm}}{)}^{\mathrm{2}}+p_{\mathrm{3}}\left(D_{\text{gm}}\right)+p_{\mathrm{4}},\end{array}$$

where the p values in Eq. (3) denote the polynomial fitting coefficients.

Table 6 Polynomial fit parameters and goodness of fit for each observed curve.

Deposition velocity increased with friction velocity in all size ranges, and several studies (Nilsson and Rannik, 2001; Pryor et al., 2008; Contini et al., 2010) found a positive correlation between $V_{\mathrm{d}}^{+}$ and the friction velocity. In particular, for u_∗ > 0.1 m s⁻¹, the relation can be described linearly. Our data confirm a linear correlation (see Fig. 9a), but the correlation coefficient varies considerably with particle size, with an angular coefficient m = 1.27 (R² = 0.9) for UFP, m = 0.78 (R² = 0.67) for ACC, and, finally, m = 7.15 (R² = 0.86) for CRS. In this context, a reasonable normalization of the deposition velocity is the ratio $V_{\mathrm{n}}=V_{\mathrm{d}}^{+}/u_{\ast }$ to obtain better comparability between different studies (Table 5). An analogous fit for the normalized V_n resulted in a similar goodness of fit and the respective fitting coefficients reported in Table 6, which have been segregated for the presence (or not) of the snowpack.

4 Conclusions

Aerosol fluxes were measured in the Arctic site of Ny-Ålesund (western Spitzbergen, Svalbard). The measurement campaign was carried out between March and August 2021. Number concentrations for the ultrafine particles increased from March to August, with new-particle-formation processes probably contributing to summertime concentrations, boosted by the enhanced biological emissions in warm months. In contrast, the number concentration for larger-sized particles shows the maximum concentration in March, during the Arctic haze season, and then it decreases to a minimum concentration in June. Turbulent fluxes in the ultrafine particle size range were prevalently downward, especially in summertime when the concentrations are highest. Particle fluxes in the accumulation and quasi-coarse mode were positive (upward) during the colder months, indicating that the site behaves, on average, as a source for particles, even if these fluxes were very small. In the warmer months, particle fluxes were slightly negative, showing that the area under examination acts as a weak sink for particulates in those size ranges. The overall median deposition velocity ($V_{\mathrm{d}}^{+}$) values were 0.90, 0.62, and 4.42 mm s⁻¹, for UFP, ACC, and CRS, respectively. A more pronounced deposition was observed, on average, during the melting and snow-free periods over the whole dimensional range, while the median deposition velocity ($V_{\mathrm{d}}^{+}$) with the snowpack was of 0.73, 0.42, and 3.50 mm s⁻¹. Our results indicate a low sensitivity of the saddle-shape of the deposition velocity (as a function of the geometric mean diameter) to the surface characteristics, showing a fixed minimum in the range of 0.1–0.2 µm and an increase in deposition velocity for particles > 0.3 µm. The deposition velocity growth is probably due to the increase in the roughness length z₀, due to the emerging of obstacles, such as rocks or irregularities and the bare soil, after melting that perturb the flow. The observed deposition velocity depends on the friction velocity; therefore, it is valuable to normalize the deposition velocity with the u_∗ to obtain a better comparison between different measurement sites. The median normalized velocity for the whole period (V_n) was 6.9 × 10⁻³, 5.1 × 10⁻³, and 35.9 × 10⁻³ for UFP, ACC, and CRS, respectively. In general, our observations of deposition velocities are aligned to the previous literature in similar environments (i.e. ice and/or snow), especially for particles in the size range 0.01–1 µm, with a percentage difference lower than 50 % in magnitude. The behaviour of V_d as a function of the particle diameter result was well predicted by the Slinn (1982) parameterization. The S82-derived models considered in this work provide different predictions with similar land covers (snow, tundra, and water ponds) and same environmental conditions, as also reported by Saylor et al. (2019). However, large discrepancies between modelled and observed deposition velocities for the accumulation mode were observed for most of the models used (i.e. Z01, PR11, and E20). The predictions of the available parameterizations generally do not agree with the observations, especially on snow and/or ice surfaces, and the discrepancy in the deposition velocity values can be as high as 2 orders of magnitude. The polynomial parameterization proposed in this study fits our observations and the analysed models, which properly represents the size dependence and magnitude of the deposition velocity. Even if this numerical fit is driven purely by the data, and not by the underlying chemical-physical processes, it could be very useful for implementing new parameterizations for existing chemical transport models considering the role of turbulence on dry deposition, which is a process typically neglected in regional and global models. Overall, our findings indicate that the most common parameterizations for dry depositions in polar areas could contribute significantly to the overall uncertainty in global models simulating short-lived climate forcer (SLCF) climate effects (Whaley et al., 2022). The scarce measurements in the cryosphere and the consequent large gap in deposition velocity modelling demonstrate the need to improve the understanding of deposition processes in polar areas.