Metamorphism of snow on Arctic sea ice during the melt season: impact on spectral albedo and radiative fluxes through snow

The energy budget of Arctic sea ice is strongly affected by the snow cover. Intensive sampling of snow properties was conducted near Qikiqtarjuaq in Baffin Bay on typical landfast sea ice during two melt seasons in 2015 and 2016. The sampling included stratigraphy, vertical profiles of snow specific surface area (SSA), density and irradiance, and spectral albedo (300–1100 nm). Both years featured four main phases: (I) dry snow cover, (II) surface melting, (III) ripe snowpack, and (IV) melt pond formation. Each phase was characterized by distinctive physical and optical properties. A high SSA value of 49.3 m² kg⁻¹ was measured during phase I on surface wind slabs together with a corresponding broadband albedo (300–3000 nm) of 0.87. Phase II was marked by alternating episodes of surface melting, which dramatically decreased the SSA below 3 m² kg⁻¹, and episodes of snowfall re-establishing pre-melt conditions. Albedo was highly time-variable, with minimum broadband values around 0.70. In phase III, continued melting led to a fully ripe snowpack composed of clustered rounded grains. Albedo began to decrease in the visible as snow thickness decreased but remained steady at longer wavelengths. Moreover, significant spatial variability appeared for the first time following snow depth heterogeneity. Spectral albedo was simulated by radiative transfer using measured SSA and density vertical profiles and estimated impurity contents based on limited measurements. Simulations were most of the time within 1 % of measurements in the visible and within 2 % in the infrared. Simulations allowed the calculations of albedo and of the spectral flux at the snow–ice interface. These showed that photosynthetically active radiation fluxes at the bottom of the snowpack durably exceeded 5 W m⁻² (∼9.2 µmol m⁻² s⁻¹) only when the snowpack thickness started to decrease at the end of phase II.

1 Introduction

Sea ice features and dynamics in the Arctic are undergoing radical changes, including a shift from multi- to first-year ice (Comiso, 2012), a decrease in thickness (Kwok and Rothrock, 2009) and in September areal extent (Comiso et al., 2017), and an earlier breakup and a later freeze-up (Arrigo and van Dijken, 2011). These changes strongly affect air–sea interactions (momentum, heat, gases) with multiple feedback loops (Serreze and Barry, 2011; Stroeve et al., 2012). They also affect marine ecosystems by substantially increasing the amount of sunlight that penetrates into the ocean and supports photosynthesis under sea ice and in the open ocean during the now longer ice-free season (Ardyna and Arrigo, 2020).

The snow cover plays a significant role in Arctic sea ice evolution. Indeed, during winter and early spring, dry snow reflects up to 90 % of incoming solar radiation, which drastically reduces the energy absorbed by the underlying sea ice (Grenfell and Maykut, 1977; Grenfell and Perovich, 2004; Nicolaus et al., 2010). Snow delays the onset of sea ice melt, while its albedo remains high enough and thus directly drives the duration of the melt season, which is itself related to the minimum sea ice extent in September (Perovich et al., 2007). The snow optical properties also control the amount of light reaching the upper ocean column under sea ice, much more so than sea ice itself (Grenfell and Maykut, 1977; Perovich, 1990). It was recently shown that major phytoplankton blooms can take place under sea ice and that snow is the main driver of bloom onset (Ardyna et al., 2020).

The albedo of snow depends mostly on the optical grain size, the density, the snowpack thickness, and the impurity content of the snow (Warren, 1982; Aoki et al., 2003). Snow is a highly scattering medium composed of ice particles that are weakly absorbing in the visible range (Picard et al., 2016a). The albedo increases when the snow particle size decreases, and these changes mostly take place in the infrared. A snow grain size metric relevant to optical studies is the optical diameter d_opt, i.e., the diameter of spheres having the same surface-area-to-volume ratio as the snow (Grenfell and Warren, 1999). At present this is obtained by measuring the snow specific surface area (SSA; m² kg⁻¹). Both are linked with the relationship $\mathrm{SSA}=\mathrm{6}/{\mathit{\rho }}_{\mathrm{ice}}{d}_{\mathrm{opt}}$, with ρ_ice the density of ice. Domine et al. (2006) showed that SSA was well related to snow reflectance in the infrared range, and SSA measurements using reflectance at 1310 nm are now routinely used (Gallet et al., 2009, 2011). Despite the great use of snow SSA, its measurements remain scarce on continental snowpacks, and even more so on sea ice, where they are limited to a few profiles (Domine et al., 2002, 2012), limiting our understanding of the albedo of snow on sea ice.

Once deposited on sea ice, snow grains undergo continuous transformations over time, known as snow metamorphism, which is mostly driven by meteorological conditions. While the snowpack is dry, the main factors responsible for metamorphism are the temperature gradient in the snowpack and wind. The temperature gradient between the warmer sea ice and the colder atmosphere leads to an upward water vapor flux coupled to sublimation–condensation cycles that lead to grain growth and the formation of faceted crystals and ultimately large hollow depth hoar crystals (Colbeck, 1983). The upward vapor flux also leads to mass loss in the basal layers so that the density of depth hoar layers often decreases over the season (Domine et al., 2016b). Wind, in contrast, leads to snow drifting and to the sublimation and fragmentation of grains so that wind processes produce hard, dense wind slabs made of small and mostly rounded grains. Since the temperature gradient is greatest during the beginning of the snow season when the snowpack is thin, and the surface of the thin sea ice is still warm, depth hoar formation is more likely then. Arctic snowpacks thus usually feature basal depth hoar layers topped by wind slabs and occasionally fresh snow before it gets wind-blown (Sturm et al., 2002; Domine et al., 2016b). Overall snow thickness ranges from a few centimeters up to 70 cm, with densities ranging from 70 to 500 kg m⁻³ (Sturm et al., 2002; Domine et al., 2002, 2012; Gallet et al., 2017). The thickest snowpacks are often found on highly deformed sea ice, where features like ridges trap wind-drifted snow (Sturm et al., 2002). When melting starts in spring, snow grains become rounded, and daily freeze–thaw cycles lead to rapid grain growth and to the formation of hard, dense, refrozen layers made of large rounded grains (Colbeck, 1973). In general, snow metamorphism leads to decreases in SSA (Legagneux and Domine, 2005) and consequently in albedo (Domine et al., 2006). Snowpack properties therefore vary over time. Given the large wind-induced spatial heterogeneity (Filhol and Sturm, 2015), the snowpack on sea ice shows large time and space variability, which makes the field study of snow properties and in particular albedo challenging because they require a lot of samplings over a relevant time period.

The snowmelt period leads to major and sudden changes over sea ice. It extends from the first surface melt event to the formation of melt ponds with typical durations ranging from 10 d to 1 month (Sturm et al., 2002; Perovich et al., 2007; Perovich and Polashenski, 2012). It can be triggered by weather conditions such as positive temperatures or rain events (Sturm et al., 2002; Nicolaus et al., 2010). Surface melting results in the formation of a thin surface layer of rounded grains, which tends to thicken with further melting. Once wet metamorphism reaches the bottom of the snowpack, the remaining snow layers melt rapidly (Gallet et al., 2017). As snow grains grow, and SSA decreases, albedo drops remarkably, by about 0.15, almost doubling the radiative energy absorption (Perovich et al., 2002; Nicolaus et al., 2010). This acts as a positive feedback enhancing further melting. The combined effects of surface melting and atmosphere warming enhance the air moisture content, often producing persistent overcast conditions leading to snow precipitation (Maksimovich and Vihma, 2012; Mortin et al., 2016). Perovich et al. (2002, 2017) and Gallet et al. (2017) observed sudden increases in albedo after such fresh snow precipitation, which suddenly increased the snow albedo and interrupted melt progression. Furthermore, sufficient summer snowfalls occasionally allow the snowpack to persist through the entire summer (Warren et al., 1999). Melt onset appears to be a chaotic transient period on the Arctic ocean, as changing meteorology can significantly lengthen the melting period.

Over the last decades, considerable effort has been made to better understand the radiative properties of snow on sea ice and their evolution across seasons. Snow albedo drives processes which control the energy budget of sea ice and the radiative flux to the underlying ocean, and albedo itself depends on snow properties. However, studies which aim to link physical and optical properties of snow never measured all variables required to calculate albedo and relied on the estimation of some physical variables, in particular snow SSA. Today, data are thus lacking to fully quantify and model the global radiative transfer in sea ice because we do not have time series of the snow properties required to understand albedo evolution. Moreover, present datasets do not include systematic combined measurements of snow optical and physical properties at the same spot. This lack of data is particularly detrimental during the melting period, when albedo is highly time-variable following alternation of freezing and melting events and precipitation.

The purpose of this paper is to contribute to filling these gaps. We simultaneously documented the temporal evolution of snow physical properties and albedo during two melting periods (2015 and 2016) on typical Arctic landfast sea ice on the east coast of Baffin Island. One or two snow pits were sampled almost every day. Measurements included, for the first time over sea ice, the time evolution of the SSA vertical profile and the corresponding spectral albedo. The stratigraphy and snow density vertical profile were also documented. In addition, the light absorption caused by light-absorbing particles (LAPs) was determined using vertical irradiance profiles in sufficiently thick snowpacks and measuring absorption coefficients of LAPs isolated from snow samples. We first aim to link surface conditions characterized by snow physical properties (SSA and density) and impurity content with albedo. The second objective is to verify that measurements of vertical profiles of snow properties enable reliable simulations of the albedo of snow-covered sea ice, especially during the chaotic melting period. Finally, we briefly compare simulated radiative fluxes at the snow–ice interface with published underwater measurements (Massicotte et al., 2020) of photosynthetically active radiation and chlorophyll a concentration to evaluate the potential of snow physical modeling in determining phytoplankton blooms under sea ice.

2 Materials and methods

2.1 Study area

Field sampling was conducted close to Qikiqtarjuaq Island in Baffin Bay (Fig. 1) from 12 May to 18 June 2015 and from 17 May to 25 June 2016 as part of the Green Edge project (Massicotte et al., 2020). All measurements were made on typical landfast sea ice a few hundred meters around an ice camp (same location for both years). In 2015 the ice was very smooth, whereas in 2016 the ice surface was deformed with small reliefs and ridges because of an early breakup in December 2015. The sampling periods were selected to capture the melting stage, but complete snow meltout and pond formations were only observed in 2016. A meteorological station was set up about 100 m from the ice camp and provided continuous measurements of 2 m air temperature and of snow thickness (Massicotte et al., 2020).

Figure 1 Location of the measurement site close to Qikiqtarjuaq Island (67^∘33^′29^′′ N, 64^∘01^′29^′′ W), east coast of Baffin Island, Canada.

2.2 Optical properties

2.2.1 Albedo measurements

Albedo measurements were performed with a custom-built radiometer (Solalb, developed at the Institut des Géosciences de l'Environnement, IGE, following Picard et al., 2016b, and Belke-Brea et al., 2020). Light was collected using a cosine collector and guided through an optical fiber to a spectrometer (Maya 2000 PRO, Ocean Optics). Irradiance was measured over the 400 to 1100 nm wavelength range with 3 nm resolution. More details about the cosine collector are given in Picard et al. (2016b). The cosine was fixed at the end of a 3 m aluminum pole, which rested on a tripod 70 cm above the surface (Fig. 2b). At the other end, the operator manually controlled the arm and triggered the spectrometer. The distance between the operator and the cosine collector (Fig. 2b) ensures minimal disturbance by the shadow. The horizontality was ensured by the operator within less than 0.3^∘ using an electronic two-axis inclinometer mounted near the cosine collector. Albedo determinations required the measurements of upwelling and downwelling irradiance, which were made sequentially using the same cosine collector with the pole being manually rotated by 180^∘. No absolute or relative calibration was needed, but measurements had to be made under steady-incident-light conditions during the 30 s process, which seldom strictly prevailed during the Arctic spring. The setup therefore included a reference photodiode to measure light fluctuations at all times for subsequent correction. A single operator could manage the entire process, including albedo measurements along linear transects.

Figure 2 Pictures of sea ice at different stages of the 2015 and 2016 melting season: (a) dry snow cover (26 May 2016), (b) albedo measurement with SOLALB (19 May 2015), (c) snow melting (spatial variability begins to be observable; drone picture from 16 June 2016), (d) melt pond formation (drone picture from 25 June 2016), (e) typical snow surface in phase II with a thin layer of fresh snow covering a 1.5 cm thick layer of wet grains (4 June 2016), and (f) picture of refrozen melt ponds on 8 June 2016 covered by a thin layer of fresh snow (glove for scale).

2.2.2 Irradiance profiles

Vertical irradiance profiles in the snow were performed in 2015 using the SOLEXS radiometer detailed in Picard et al. (2016a). Snow dunes were selected for these measurements because thicker layers are more suitable to accurately determine irradiance profiles. The optical absorption by LAPs in the snowpack was determined from the exponential rate of irradiance decrease in a homogeneous layer (Belke-Brea et al., 2021; Tuzet et al., 2019) so that thicker layers yield more reliable measurements. Another reason for selecting thick snowpacks is that irradiance data in the top 7 cm of the profile had to be discarded because of near-surface perturbations (Picard et al., 2016a; Belke-Brea et al., 2021). Briefly, an optical fiber was slowly lowered through a hole preformed in the snowpack. The fiber was guided by a vertical rail mounted on a four-legged structure (Picard et al., 2016a). A profile requires at least 1 min to complete, which was a major issue under changing light conditions. These were monitored by the same photodiode as for albedo measurements. Measurements were repeated until a complete profile could be recorded under constant incident light. Two such valid profiles were recorded for each snow pit. They were averaged and then interpolated over 5 mm intervals. All final profiles were normalized to the irradiance at a depth of 7 cm (Picard et al., 2016a). Analysis of the profiles yields the spectral absorbance of LAPs, which was used to determine the type of absorbing particles (essentially black carbon, BC, or mineral dust). LAPs were also analyzed using chemical methods to obtain absorption coefficients, as detailed in the following section.

2.2.3 Light-absorbing particles in snow

Absorption spectra of light-absorbing particles in snow were measured during the 2016 field campaign. Twelve snow samples (500 mL each) were collected at random locations around the ice camp before complete meltout. Samples were brought to the lab, where they were melted within a closed cooler placed at ambient indoor temperature. For particulate absorption, we followed the so-called T–R method developed by Tassan and Ferrari (1995, 2002). Briefly, water samples were filtered under low vacuum on Whatman GF/F glass fiber filters immediately after snow melting. Optical density (dimensionless) of all particles on these filters, referred to as OD_p(λ), was then determined from transmittance and reflectance measurements from 200 to 800 nm with 1 nm resolution using a Cary 100 spectrophotometer (Agilent Technologies). A baseline correction was applied by subtracting the optical density of a fully hydrated blank filter from OD_p(λ). A correction for pathlength amplification was applied as proposed by Tassan and Ferrari (2002). We then converted these OD_p values into absorption coefficients for total particles (a_p(λ); m⁻¹) using the following equation:

$$\begin{array}{}\text{(1)}& a_{\mathrm{p}}\left(\mathit{\lambda }\right)=\mathrm{2.3}{\mathrm{OD}}_{\mathrm{p}}\left(\mathit{\lambda }\right){\displaystyle \frac{\mathrm{FA}}{V}},\end{array}$$

where FA is the filtration area on the glass fiber filter (m²), and V is the volume of the melted snow sample (m³). The 2.3 factor accounts for the conversion from base 10 (i.e., OD_p) to base e (i.e., a_p) logarithm.

2.3 Snow physical properties

Here, snow physical properties refer to temperature, snow grain shape and visual size, SSA, and density in each snow pit. These were performed in 2015 and not in the melting 2016 snowpacks. We first identified the main stratigraphic layers by visual inspection. For each layer, the average visual snow grain size and shape were determined using a hand lens. Snow temperature was measured with a Pt100 thermistor probe (0.1 ^∘C accuracy) at several depths from the bottom of the cover to up to 10 cm beneath the surface. Negative freeboard was reported when the sea level was above the snow–ice interface. The vertical profile of snow density was measured using a 100 cm³, 3 cm high box cutter. The collected snow sample was weighed using an electronic scale with a resolution of 0.1 g. According to Conger and McClung (2009), this method allows snow density measurements with an accuracy of 11 %. The main uncertainties concern the real volume extracted by the cutter, which depends on the type of snow and is larger for soft snows such as depth hoar and fresh snow. Finally, vertical profiles of SSA were determined from the snow IR reflectance using the DUFISSS instrument (Gallet et al., 2009). Briefly, DUFISSS measures the albedo of a cylindrical snow sample 63 mm in diameter and 25 mm thick at 1310 nm with an integrating sphere. The SSA is deduced from the albedo using a polynomial relationship. The correction concerning the determination of SSA of wet snow introduced by Gallet et al. (2014a), which mostly consists of adding 0.5 m² kg⁻¹ to the value of SSA, was not applied in this study because it did not induce significant changes in albedo simulations (less than 1 % albedo increase at 1000 nm in most cases and negligible change in the visible). The uncertainty in SSA determinations is 12 % under most conditions (Gallet et al., 2009). Melting can occur if the sample is not handled fast enough under warm meteorological conditions, which leads to a lowered SSA value. Special care was taken to keep sampling tools as cold as possible, e.g., by cooling instruments in bottom snow layers when the surface was melting.

2.4 Sampling protocol

Data presented in this study were collected either in snow pits or along transects.

• Snow pits. Albedo was measured first since it requires a pristine area. A minimum of three measurements were made depending on sky conditions and light variations. All of them were performed facing the sun to avoid any shadow from the operator and the equipment. Irradiance profiles were then measured when thick snowpacks were studied, with minimal perturbation to the snow surface. When irradiance profiles were performed, the pit was then dug so as to intersect the hole made for the irradiance profile. All stratigraphic measurements were carried out along a 1 m long shaded trench. Our objective was to conduct all samplings at the same place in order to fully characterize physical and optical properties of the snow at each station.

  One or two snow pits (each requiring 3 h of work, more for thick snowpacks) were sampled each sampling day. Fewer snow pits were sampled in 2016 (10 versus 35 in 2015) because the snowpack was already ripe (i.e., isothermal at 0 ^∘C and melting throughout) before sampling operations. Snow pit locations were randomly chosen around the ice camp. Particular efforts were made to sample the widest possible range of snowpack depths in order to address spatial variability.

• Transects. Albedo and snow depth were also measured every 5 m along transects (100 to 150 m long). All the equipment was placed on a sled to facilitate transport. These transects were performed in 2016 only. As the snowpack was already ripe, the study of spatial variability using transect measurements was favored.

2.5 Data processing

Upwelling and downwelling irradiance raw acquisitions require several processing steps before the albedo can be obtained. During the field campaigns, spectra were visually checked at the end of the sampling day. Unrealistic data, based on qualitative criteria (spectra shape, range of values, signal-to-noise ratio), were rejected. The first step of processing was to remove the systematic offset in both acquisitions caused by dark current and stray light effects. This offset was approximated for each acquisition as the mean signal at short wavelengths (between 200 and 260 nm) because there is no incoming photon in this wavelength range. Dark current was assumed to be constant over the entire wavelength range. Then, spectra were divided by their corresponding integration time. Our cosine collectors had been previously characterized on an optical bench in order to assess their exact angular response (Picard et al., 2016b). This angular response was then used to correct the downwelling irradiance measurements depending on the solar zenith angle (SZA) during the acquisition if the sky was clear. No correction was needed for overcast skies. We excluded any acquisition for which the reference photodiode signal varied by more than 2 % between the upwelling and downwelling irradiance measurements. Below 2 %, spectra were rescaled using the reference photodiode signal assuming that changes in incident light were equivalent over the entire wavelength range. After all these steps, albedo was calculated as the ratio of upwelling over downwelling irradiance. Albedo spectra were finally smoothed using a low-pass filter in order to eliminate the noise inherent in the measurements. For each measurement site, it was checked that all spectra correctly overlapped before being averaged. For the 2015 dataset, the average standard deviation of all integrated albedos (over the 400–1000 nm wavelength range) measured at each snow pit is 0.3 %, with a maximum of 1 %. Thus, in most cases, it is reasonable to assume that the precision in albedo measurements is better than 1 %. This uncertainty does not account for external error sources such as varying illumination conditions or unexpected shadows or reflections.

2.6 Radiative transfer modeling

The numerical simulations of albedo and irradiance profiles were performed using the Two-stream Analytical Radiative TransfEr In Snow (TARTES) model (Picard, 2022, with details in Libois et al., 2013). Briefly, TARTES uses the delta Eddington approximation (Jimenez-Aquino and Varela, 2005) in a layered plane-parallel snowpack. Each layer is characterized by an average SSA, density, and impurity concentration. Following the asymptotic radiative transfer theory (ART) (Kokhanovsky and Zege, 2004), light scattering in snow is also determined by the absorption enhancement parameter B, which characterizes multiple reflections within a snow grain, and by the asymmetry factor g, which describes the preferential scattering direction (e.g., forward or backward).

TARTES solves the radiative transfer equation at all depths. The underlying sea ice is not modeled; only its albedo (measured in the field) is specified at the bottom of the snowpacks. The snow albedo depends on SZA and cloud cover, but a fully diffuse radiation is equivalent to a direct radiation with SZA ∼50^∘ (Warren, 1982). In our study, SZAs were between 47 and 57^∘; therefore simulations were performed by approximating the illumination conditions with diffuse radiation (SZA of 53^∘ in TARTES, according to ART). Doing so, the maximal error in albedo is ∼0.01 at 1000 nm.

The use of TARTES allows the calculation of albedo on a wide wavelength range, which makes the assessment of the broadband albedo α and total energy absorbed at the snow surface A (W m⁻²) possible. Both were calculated as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }=\underset{\mathrm{300}}{\overset{\mathrm{3000}}{\int }}{\mathit{\alpha }}_{\mathrm{s}}\left(\mathit{\lambda }\right)I\left(\mathit{\lambda }\right)d\mathit{\lambda }/\underset{\mathrm{300}}{\overset{\mathrm{3000}}{\int }}I\left(\mathit{\lambda }\right)d\mathit{\lambda }\text{(3)}& {\displaystyle }& {\displaystyle }A=\underset{\mathrm{300}}{\overset{\mathrm{3000}}{\int }}\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}\left(\mathit{\lambda }\right)\right)I\left(\mathit{\lambda }\right)d\mathit{\lambda },\end{array}$$

where α_s is the spectral albedo calculated with TARTES over the 300–3000 nm wavelength range, and I(λ) is the spectral solar irradiance (W m⁻² nm⁻¹). The solar irradiance spectra were calculated with the SBDART model (Ricchiazzi et al., 1998) for 1 June at 12:00 local time (LT) at Qikiqtarjuaq under clear-sky conditions of Arctic spring on snow-covered areas. Only one solar spectrum was used since the study focused on the broadband albedo more than on the energy budget.

Impurities in TARTES are assumed to be external to the ice grains, and the particle size is small compared to the wavelength. Their absorption was taken into account as detailed in Kokhanovsky and Zege (2004). TARTES is therefore able to consider any kind of LAP whose refractive index is known (Tuzet et al., 2019). Concentrations and types can be adjusted for each layer of the snowpack.

Finally, since snow layers near the surface and down to a depth of at least 10 cm were almost always windpacks or clustered rounded grains, we follow the recommendation of Libois et al. (2014) and use B=1.25 and g=0.895.

3 Results

In this section, the 2015 dataset is entirely presented with all the 35 snow pits which comprise vertical profiles of SSA and density, as well as albedo measurements. Only the 2015 data were used for radiative transfer modeling with TARTES. The 2016 dataset is only composed of albedo measurements and qualitative observations on snow and ice. It is specified in the text whenever we refer to the 2016 dataset.

3.1 General evolution and meteorological conditions

Surface conditions changed drastically during both sampling campaigns as depicted in Fig. 2. About 1 month elapsed between the first day of surface melting and meltout. As previously observed by Perovich et al. (2002) and by Nicolaus et al. (2010), as the melting season progressed, the sea ice surface became darker, and spatial variability increased. The time evolution of albedo at 500 and 1000 nm is presented in Fig. 3, and, similarly to Perovich et al. (2002) and Nicolaus et al. (2010), this evolution clearly shows four main stages confirming visual observations in the field. These phases are defined below.

Figure 3 Time evolution of albedo at 500 nm in blue and 1000 nm in red for the 2015 and 2016 field campaigns in darker and lighter colors, respectively. Boxes are used if more than one measurement per day is available; colored dots are used otherwise. The main phases are specified by horizontal gray bars for both years (phases I, II, and III for 2015 and phases II, III, and IV for 2016). Gray bar graphs at the bottom represent the number of albedo measurements per day. In 2015, the campaign ended on 18 June, so no data are available beyond that date.

Figure 4 Continuous measurements of the sea ice meteorological station. Details are given in Massicotte et al. (2020). (a) Time evolution of air temperature and (b) snow thickness measured for Green Edge 2015 (dark gray) and 2016 (light gray). The zero-level snow depth is that measured just after meltout. Horizontal gray bars in (a) denote the sampling periods for both campaigns. Additionally, the main snowfalls in 2015 are specified in (b) with a distinction between conventional snowfalls and light nighttime snowfalls observed at the end of the season.

• Phase I – cold, dry snow (from the first sampling day on 13 to 24 May 24 2015). Sea ice was covered by a dry winter snowpack never subjected to melting. Air temperature increased during this phase but remained below 0 ^∘C (Fig. 4). A significant snowfall event associated with strong winds occurred before the first sampling day in 2015 (8 and 9 May), building a fresh snow layer at least 10 cm thick. Temperature in snow was first colder at the surface, or at least at mid-depth (−6.5 ^∘C), than at the bottom-most layer, where temperatures remained fairly steady between −5 and −4.5 ^∘C in the daytime. The subsequent increase in air temperature led to surface warming and a reversal of the temperature gradient in the snowpack during this first phase.

• Phase II – surface melting (25 May to 11 June 2015, from the first sampling day on 19 May to 9 June 2016). This phase started with the first surface melting event, which coincided with the first positive air temperature in 2015 (Fig. 4). Coarse rounded grains and wet grains appeared, and albedo decreased in the infrared (Fig. 3). Air temperature fluctuated around 0 ^∘C, and several snowfalls were observed both years during that period (Fig. 4; snowfalls specified only for 2015). Moreover, the weather was cloudier than during the previous phase, and thick fog was more common in the early morning. These meteorological conditions persisted in the next phases. Overall, snow temperatures gradually increased until the 0 ^∘C isothermal state was reached on 10 June 2015.

• Phase III – ripe snowpack (12 June to last sampling day on 18 June 2015, 10 to 17 June 2016). At this stage, all the snow layers were at the melting temperature and were entirely comprised of rounded polycrystals. This phase is characterized by a decrease in albedo over the visible range for the first time of the season (Fig. 3). Snowpack thickness decreased very quickly until meltout (6 d in 2015, 7 d in 2016). In 2015, measurements were stopped on 17 June so that the end of this phase was not studied. The first melt ponds appeared on 19 June 2015.

• Phase IV – melt pond formation (in 2016: 18 June to the last sampling day on 26 June). Snowpacks gave way to a mixture of bare ice and melt ponds. The transition from snow cover and bare ice was progressive because the ice surface was granular and looked similar to the large wet grains observed on the snowpack in the ultimate stages of snowmelt. As previously observed, sea ice was first rapidly flooded by extended shallow ponds before they partially drained and got their final shape. During our last sampling day in 2016, 25 June, a cooling event associated with snowfall temporarily froze the ponds (Fig. 2f) and increased albedo (Fig. 3).

3.2 Snow stratigraphy and physical properties

Only physical properties sampled in 2015 are presented here because they cover the main first three phases, unlike in 2016. Figure 5 presents the main stratigraphic layers observed for each phase. Regarding the physical properties, vertical profiles of SSA and density are presented in Fig. 6, with average values in Table 1.

Figure 5 Main stratigraphic layers observed for each phase and for snow dunes in 2015. The lower layers with low SSA are represented in dark gray. The upper layers with high SSA during phase I are in light gray. Layer numbers are those used in Table 1. A vertical scale is not provided as snow depths were highly variable.

Figure 6 Vertical profiles of SSA (a) and density (b) for each snow pit sampled in 2015 with snow elevation (in centimeters) on the y axis and sampling dates on the x axis. Snow dunes are labeled with hatches, and phases I to III are specified. White regions mean no data.

Table 1 Average values of SSA and density and corresponding standard deviations for each phase (SSA: m² kg⁻¹; density: kg m⁻³). Layer numbers refer to those in Fig. 5. The number of samples for each value is given in parentheses.

• Phase I – cold, dry snow. The observation of 15 snow pits during this phase revealed a dominant stratigraphy composed of three or four main layers. The bottom-most layer (layer I in Fig. 5), in contact with the underlying sea ice, was indurated depth hoar formed from a wind slab (Derksen et al., 2009; Domine et al., 2016b), as evidenced by the presence of depth hoar crystals embedded in a matrix of small rounded grains and confirmed by its high density of 372±51 kg m⁻³. Its SSA was 8.9±4.4 m² kg⁻¹, typical of depth hoar, whether indurated or not (Domine et al., 2016b). Above this basal layer, a layer of indurated faceted grains (layer II) with average SSA of 12.1±1.8 m² kg⁻¹ and average density of 409±40 kg m⁻³ was observed. The upper part of the snowpack was first comprised of a wind slab layer (layer III) made of small rounded grains characterized by significantly higher SSA values, 33.4±2.6 m² kg⁻¹, and lower densities, 276±38 kg m⁻³. Having basal indurated depth hoar denser than upper wind slabs is not rare on sea ice (Sturm et al., 2002). Occasionally a layer of dendritic crystals or decomposing particles could be observed on the surface (layer IVa). The highest SSA values were recorded in this layer, 49.3±5.9 m² kg⁻¹ (see dark-red areas at the surface in Fig. 6a). Moreover, sublimation crystals (Gallet et al., 2014b) sometimes formed at the surface of the snowpacks when an upward flux of water vapor condensed into colder air at the snow surface. Figure 6 also shows a significant contrast in both vertical profiles, with layers I and II characterized by lower SSA and higher density than layers III and IV. Moreover, SSA in layer III gradually decreased over time. Overall, snow depth ranged from 15 to 54 cm. Snow dunes were studied on 19, 22, 23, and 29 May and 4 June. They corresponded to the thicker snowpacks. Their stratigraphy revealed a similar sequence as elsewhere, with depth hoar, faceted crystals, and wind slab. However, the wind slabs had a much lower SSA than elsewhere, suggesting an earlier formation. The highest densities, reaching 500 kg m⁻³, were measured in layer II of dunes (Fig. 6). Smaller features such as sastrugi (Fig. 2a) and barchan dunes were commonly observed before melt onset. Freeboard was always positive during phase I. Since no major modification affected the snowpack during phase I, except near the surface, a lot of the variability observed in the stratigraphies of Fig. 6 during phase I results from spatial rather than time variability.

• Phase II – surface melting. The first melting was observed on 26 May, 1 cm below the surface, and coincided with a low SSA layer at that depth (Fig. 6). Overall, surface melting was characterized by the formation of a layer of rounded polycrystals (layer Va; Fig. 5), of low SSA (10.6±4.1 m² kg⁻¹). Additionally, as melting conditions persisted this layer got thicker, and its SSA kept decreasing to a minimum of 2.6 m² kg⁻¹ on 13 June (phase III). The alternation of negative and positive temperatures during night and day subjected the surface of the snowpack to a diurnal cycle. During daytime, at the surface, bonds between snow grains melted, leading to the observation of wet, clustered, rounded grain, which partially (at least near the surface) froze during the following night, forming again dry rounded polycrystals and often a thin melt–freeze crust at the surface. Several snowfalls (16, 17, 26, 27, and 29 May and 6 June; Fig. 4) deposited fresh snow layers covering layer Va, which then quickly metamorphized (see the second stratigraphic profile of phase II; Fig. 5). Some of these new layers were thick enough to be distinguishable in Fig. 6 because of their high SSA. Fresh snow tended to accumulate in depressions rather than on top of dunes. Melting and subsequent refreezing increased cohesion between snow grains. This totally stopped erosion of snow by wind. As the weather became cloudier, a thin layer of surface hoar or precipitated needle crystals deposited during the night (3 and 10 June; see Fig. 4) were regularly observed at the beginning of the day before they rapidly melted. The underlying snow layers I and II, unaffected by surface melting, remained nearly unchanged (with SSAs of 10.6±4.1 and 13.8±6.9 m² kg⁻¹ and densities of 370±26 and 418±51 kg m⁻³ for layers I and II, respectively). The SSA of layer III (24.7±4.3 m² kg⁻¹) kept on decreasing during phase II (Fig. 5) until it had completely transformed into wet grains (phase III). Solid melt–freeze layers, almost ice-like (layer Vb; Fig. 5), were first observed on 29 May and became more and more common, to the point that they were present everywhere at the end of phase II, and several of them could be found in the same snow pit. Actual ice layers, formed by percolating meltwater that froze at interfaces between layers because of the change in capillary forces, were also observed starting on 6 June (not shown in Fig. 5). After 4 June, the freeboard became frequently negative, possibly due to the increasing mass of snow at the surface but also likely due to melting bottom ice caused by the warming sea.

• Phase III – ripe snowpack. At this stage, the snowpack was only composed of coarse rounded grains with the lowest SSA values recorded (4.6±1.2 m² kg⁻¹). It was isothermal at 0 ^∘C, and its thickness decreased rapidly. In 2016, in contrast to 2015, a layer of liquid water up to 10 cm thick was found nearly everywhere at the base of the snowpack before snow meltout. That layer was likely the result of the imbalance between the rapid input of snow meltwater and the low drain capacity of sea ice rather than the result of a negative freeboard. The water was slightly salty to the taste, but much less than seawater. The salt may have come from the brine that was inevitably present at the ice surface when ice formed in the fall and wicked up the first snow layer by capillary rise (Alvarez-Aviles et al., 2008; Krnavek et al., 2012).

3.3 Spectral albedo

All albedo spectra from the 2015 and 2016 field campaigns are summarized in Fig. 7. They are displayed by phase in order to better illustrate their corresponding specific spectral signatures. Mean albedo values at 500 and 1000 nm are specified in Table 2 for each phase.

Figure 7 Spectral albedo from 400 to 1000 nm for both 2015 and 2016 campaigns, represented with boxplot graphs and sorted by phase: (a) cold, dry snow (13 samples); (b) surface melting (44 samples); (c) ripe snowpack (87 samples); and (d) melt pond formation (106 samples, most of them mixtures of ice and pond). In (d) one albedo spectrum over bare ice and one spectrum over a melt pond are also shown.

Table 2 Mean albedos and corresponding standard deviations at 500 and 1000 nm during each phase. In phase IV, measurements above ice and pond only (one station for each) are specified.

• Phase I – cold, dry snow. The highest albedos, 0.97±0.01 at 500 nm, were measured above cold winter snow (Figs. 3 and 7). Values slightly decreased during this phase almost only in the infrared from 0.80 to 0.70 at 1000 nm. Spatial variability was low, and the lowest albedo, 0.68 at 1000 nm, was recorded only above snow dunes, where fresh snow did not accumulate. The slopes of the dunes where measurements were made was barely perceptible, most likely lower than 1^∘ in all cases. Picard et al. (2020) have evaluated slope effects on albedo. Under diffuse light, over entirely snow-covered surfaces, a slope has no impact on albedo (their Fig. 3). Considering a SZA of 45^∘ and only direct light, the error caused by a 1^∘ slope is slightly over 0.01. Given that the slopes involved here are <1^∘ and that there is always significant diffuse light, we estimate that the slope-caused error is always <0.01, and in most cases much smaller so that we neglect it.

• Phase II – surface melting. Following the onset of wet snow metamorphism at the surface, albedo declined mostly in the infrared down to 0.45 at 1000 nm (Fig. 7), while it remained high in the visible, with all values > 0.95 at 500 nm. Figure 7 shows large variations in albedo at 1000 nm. In particular, sudden increases are observed after snowfalls (30 May and 4 June 2015, for instance; Fig. 3). These increases brought the albedo back near values observed in phase I. Despite the wider range of albedo in phase II compared to those in phase I (Fig. 7), spatial variability did not increase during this phase because changes in snow SSA were homogeneous over the sea ice surface.

• Phase III – ripe snowpack. Darker patches and lower albedos were observed (Fig. 2c) as snow thickness declined. Albedo decreased in the visible range, from 0.95 to 0.65 at 500 nm, while comparatively it remained steadier in the infrared, 0.43±0.042 at 1000 nm (Fig. 3). In Figs. 3 and 7c, the albedo ranges in the visible are getting wider over time, showing that space variability appeared in this phase. This spatial variability in albedo followed spatial variability in snow thickness. Albedo in the visible was lower above thinner snowpacks.

• Phase IV – melt pond formation. The transition between phase III and IV was not sharp because melt ponds appeared suddenly, while snow patches were still remaining. The albedo decreased over the whole spectral range (Figs. 3 and 7d) down to 0.74 and 0.25 at 500 and 1000 nm, respectively, over bare ice and 0.32 and 0.06 over melt ponds (Table 2). Spatial variability was maximal during this stage. The cooling event that brought less than 1 cm of snow on 25 June 2016 temporarily enhanced albedo (Fig. 3). This increase was largest in the near-infrared.

3.4 Radiative transfer modeling

3.4.1 Snow impurity content

The particulate absorption spectra in melted snow samples collected during the field campaign were all very similar in terms of shape but differed in their magnitudes, which ranged from 0.03 to 1.10 m⁻¹. The average spectrum derived from the 12 samples is shown in Fig. 8 as a thick gray line. It was divided by the density of water (1000 kg m⁻³) and multiplied by the average density of snow (350 kg m⁻³) in order to obtain an average absorption coefficient of the impurities in the snow. The shape of this spectrum is fairly similar to that of mineral dust (MD) (Wagner et al., 2012; Caponi et al., 2017), indicating that MD was the main LAP in snow at our site. The MD likely originated from the land around the ice camp, as granite cliffs are very common near Qikiqtarjuaq. The first simulations showed that simulations could be improved by adding black carbon (BC) as an extra LAP. Refractive indices of the MD and BC, derived respectively from Wagner et al. (2012) and Bond and Bergstrom (2006), were used in TARTES to simulates the LAPs. As shown in Fig. 8, a mixture (dashed red line) of 1700 ng g⁻¹ of MD and 14.4 ng g⁻¹ of BC was found to give the best values to correctly simulate, in terms of shape and amplitude, the average particulate absorption measured within the snow.

Figure 8 Absorption coefficients of snow impurities (m⁻¹), measured in melted snow. The thick gray line denotes the average of all particulate absorption measurements for snow with a density of 350 kg m⁻³. At 500 nm, values ranged from 0.03 to 0.42 m⁻¹. Blues lines denote the coefficients of MD (dashed) and BC (dotted) and the sum of both (dashed red), used to match the measurements, with 1700 ng g⁻¹ of mineral dust and 14.4 ng g⁻¹ of black carbon.

As particulate absorption of snow was only measured during the 2016 campaign, we used our 2015 vertical irradiance profiles obtained in five snow dunes in order to investigate the presence of such impurities within the snow in 2015. LAP concentrations were obtained from the exponential rate of irradiance decrease in a homogeneous layer following Tuzet et al. (2019) and Belke-Brea et al. (2021). Making our calculation on the most linear parts of the five logarithmic profiles, we found an average concentration of 886 ng g⁻¹ of MD and 14.1 ng g⁻¹ of BC in the top layer of the dunes. Those concentrations are within a factor of 2 of (for MD) or similar to (for BC) those derived from the LAP absorption study the following year.

3.4.2 Albedo modeling

Only data from the 2015 campaign are presented since this dataset is more comprehensive and covers the first three phases. As mentioned in the previous section, snow grain size impacts albedo mostly in the near-infrared, while the impact of snow thickness is mainly observed in the visible. Thus, both measurements and simulation results are compared at 500 and 1000 nm. They are presented in Table 3 and Fig. 9, where simulations with and without LAPs are shown.

Figure 9 Albedo measurements (black squares) and simulations (red and blue dots) at 500 nm (a) and 1000 nm (b) for each sampling station in 2015 (different scale on y axis). Blue dots are simulations without impurities. Red dots include impurities (mineral dust and black carbon). At 500 nm, impurities decrease simulated albedo by 1 % to 2 %, whereas at 1000 nm their effect is hardly detectable so that the red curve conceals the blue one. The shaded gray areas indicate the surface melting period (phase II).

Table 3 Relative deviation between albedo simulations and measurements at 500 and 1000 nm (in percentages) and corresponding standard deviation.

Albedo simulations were performed to assess the adequacy of the snow observations collected for radiative transfer modeling and to quantify the sensitivity of surface albedo to snow surface properties and snow depth. Albedo was simulated for each snow pit using vertical profiles of SSA and density (Fig. 6) as inputs in the TARTES model. A single impurity concentration was used for all albedo calculations because we had no way to determine specific concentrations in each snowpack. It is possible to fit each spectral albedo with optimized LAP concentrations, but without measurement in each pit, this is not very meaningful. We used a mixture of 2000 ng g⁻¹ of MD and 20 ng g⁻¹ of BC for all snow pits. These values were chosen to optimize the albedo fit at 500 nm during phase I, where snow physical measurements are the most reliable. Increments of 500 ng g⁻¹ were tested for MD. The 100:1 ratio between MD and BC was chosen for consistency with measurements from irradiance profiles in snow dunes (ratio 118:1) and data from melted snow (ratio 63:1). Impurity concentrations are not necessarily identical in dunes and in the rest of the snowpack for many reasons, one of which being that layers in dunes did not necessarily deposit at the same time as layers with similar grain shapes elsewhere. Furthermore, LAPs and in particular MD particles were probably for the most part deposited during wind events, and they may have been diluted in the thicker dune layers relative to elsewhere. Lastly, LAP concentrations in 2016 may have been slightly different from those in 2015. The values used, although slightly higher than measurements, are thus reasonable.

The largest differences are for the last two dates. These are during phase III, when melting was extensive and rapidly evolving, with visible short-scale variations. Since albedo measurements cover about 1 m², while density and SSA values are from point measurements, an imperfect coincidence between both datasets is expected. Albedo and snow pit measurements are also done in sequence so that variations between both measurements are possible, perhaps even likely. If we remove these last two measurements from the comparison, Fig. 9 shows that at 500 nm, the difference between measurement and model is <0.01 (absolute value) for 25 out of 31 values. At 1000 nm, the difference is <0.04 (absolute value) also for 25 out of 31 values and <0.02 for 16 values. During phase I, average simulations at 500 nm are essentially identical to measurements, while during phase II they are below measurements by 0.008. At 1000 nm, these values are 0.036 and 0.009.

Given these reasonably good fits, simulations were used to calculate three metrics relevant to radiation and energy of the snow–sea-ice system: the broadband albedo, the broadband radiative energy input into the snow, and the flux of visible radiation at the snow–ice interface. Calculations were performed for two wavelength ranges: the solar spectrum (300–3000 nm, i.e., broadband) and photosynthetically active radiation (PAR; 400–700 nm). This last metric is relevant to understanding phytoplankton blooms under sea ice (Ardyna et al., 2020). Time series of these metrics are plotted in Fig. 10. Usual units for PAR are photons fluxes, and these are shown in Fig. 11. Units of watts per square meter are used in Fig. 10 for consistency and to allow easy comparisons.

Figure 10 Simulated radiative metrics of the snow–sea-ice system. (a) Broadband albedo (300–3000 nm) of the system. (b) Broadband radiative fluxes, based on downwelling radiation calculated on 1 June 2015 using the SBDART model. (c) Photosynthetically active radiation (PAR: 400–700 nm) flux at the snow–ice interface.

Figure 11 Simulated PAR fluxes at the snow–ice interface, compared with PAR fluxes under the sea ice at 1.3 m depth and with chlorophyll a concentrations integrated over the water column. These last two datasets are from Massicotte et al. (2020).

The broadband albedo reached a maximum of 0.87 during phase I (dry snow). During phase II (surface melting), its minimum value was 0.76 near the end of the period, inducing an increase in the solar input from 103 to 185 W m⁻² under an incident clear-sky flux of 784 W m⁻², as detailed in Sect. 2.6. The succession of snowfalls and melting episodes caused significant variations in solar radiation transmitted to the system. The layers of fresh snow reduced the energy input by approximately 20 W m⁻² on 30 May, 16 W m⁻² on 4 June, and 11 W m⁻² on 10 June. On the other hand, the rapid metamorphism of the snow resulted in an increase in solar input of 50 W m⁻² from 30 May to 1 June and 45 W m⁻² from 4 to 6 June. Regarding PAR, the flux at the snow–ice interface remains below 2 W m⁻² most of the time during phases I and II. This is because essentially all visible radiation is reflected. Figure 10c indeed shows that the PAR flux at the snow surface is <10 % of the broadband flux. Only when the snowpack thickness decreases below about 20 to 25 cm at the end of phase II can significant visible radiation reach the ice surface and be transmitted downwards. The PAR flux on 13 June reaches 22 W m⁻².

PAR fluxes are used to evaluate biological productivity, and considering them is especially important in this case to relate their increase to the initiation of the phytoplankton bloom in the water column. Figure 11 plots PAR fluxes at the snow–ice interface simulated in this study together with underwater PAR fluxes at 1.3 m depth and water-column-integrated chlorophyll a amounts, reported by Massicotte et al. (2020).

4 Discussion

The previous section details the evolution of physical properties and albedo of snow. It also tests our ability to retrieve albedo by radiative transfer modeling using measured snow properties. In this section we first propose to reconstruct the main steps of formation of the snowpacks observed during samplings. We then discuss the co-evolution of albedo and of snowpack characteristics. Finally, we discuss possible future improvements in the accuracy of albedo simulations.

4.1 Snowpack formation

Arctic snowpacks, whether on land or sea ice, almost always show the same type of structure, with basal layers of depth hoar and sometimes faceted crystals and top wind slabs with occasional fresh or recent snow that usually gets wind-blown to add more sub-layers to the top wind slab (Sturm et al., 2002; Derksen et al., 2009; Domine et al., 2012, 2016b). Depth hoar forms early in the season, when the elevated temperature gradient between the warm soil or young, thin sea ice and the cold air is large (up to 200 K m⁻¹) and generates upward water vapor fluxes that allow the growth of large depth hoar crystals. Depth hoar can be soft or indurated, depending on the initial density of the snow when it started forming (Domine et al., 2016b). On the North American–Greenland side of the Arctic, there is usually more precipitation and snow accumulation in fall and spring, whereas winter is rather dry (Warren et al., 1999; Sturm et al., 2002; Domine et al., 2021). The depth hoar layers (layer I; Fig. 5) thus most likely formed in fall. They were indurated and of high density, indicating that they formed from dense wind slabs. The layers of faceted crystals (layer II) formed subsequently also in the fall, when the temperature gradient had decreased because the snowpack was thicker, and the surface of the sea ice was colder because it was thicker. The fact that they formed in dense wind slabs probably explains why they did not reach the depth hoar stage, as crystal growth under temperature gradient metamorphism is slower at higher snow densities (Marbouty, 1980; Flanner and Zender, 2006). The density of the lower layers indicates windy meteorology in the fall, with winds frequently exceeding 10 m s⁻¹ (https://weatherspark.com/h/y/147443/2015/Historical-Weather-during-2015-at-Qikiqtarjuaq-Airport, last access: 7 March 2022). The windy character of the area in the fall is confirmed by the presence of whaleback dunes (Filhol and Sturm, 2015), which probably formed mostly in the fall, given the predominance of depth hoar and faceted crystals in their stratigraphy.

Wind slabs (layer III) formed when the temperature gradient had decreased too much to allow the formation of depth hoar or even faceted crystals. Wind slabs may have formed on whalebacks as early as late fall, since the snow thickness reduced the temperature gradient. However, Fig. 4 shows significant accumulation in spring, and this is probably when most of the wind slabs formed. This is confirmed by their fairly high SSA, often >30 m² kg⁻¹. This is higher than for typical wind slabs (Domine et al., 2016b, 2012), suggesting the SSA did not have time to decrease after formation.

Depth hoar almost always has a SSA close to 10 m² kg⁻¹, while the SSA of faceted crystals is not much higher, often <15 m² kg⁻¹ (Domine et al., 2002, 2016a; Taillandier et al., 2006) as observed here (Fig. 6). Wind slabs have a higher SSA, usually in the range of 25–35 m² kg⁻¹ in most places, and this enhanced albedo in the infrared (Wiscombe and Warren, 1980; Domine et al., 2006). However, we note with interest that the lowest SSA values for wind slabs are found on whalebacks, with values often <25 m² kg⁻¹. Since SSA of rounded grains decreases over time (Taillandier et al., 2007), this tends to confirm that the wind slabs of whalebacks formed earlier than elsewhere, and probably in late fall as suggested above.

In summary, most of the whaleback snow height accumulated in fall under strong winds and initially elevated temperature gradients, producing depth hoar and faceted crystals. A significant part of the general snowpack also formed in fall, and the strong temperature gradients also produced depth hoar and faceted crystals there. Wind slabs probably started forming in fall, but the high SSA in most spots indicate recent formation, in spring, just before or during our campaigns (Fig. 4). Based on general western Arctic climatology and on the snow depth data of Fig. 4, we propose that snow accumulation in winter was minimal.

4.2 Albedo and surface evolution

Our data illustrate the co-evolution of albedo and of the state of the snowpack. Albedo was initially high as long as the snowpack remained dry. It gradually decreased after melt onset as surface snow grains coarsened, and this decrease continued until the snowpack completely vanished. During our campaigns, albedo was affected by synoptic weather, such as snowfalls. Snowfalls were common during the melting period. They suddenly and temporarily increased the albedo, substantially delaying meltout.

The data shown in Figs. 7 and 9 show a different evolution of spectral albedo during each phase. Phase I sees little change in albedo, either in the visible or infrared ranges, because the snowpack remains dry. Its thickness does not change, and the small drop in snow SSA (Fig. 6) just leads to a modest drop in the infrared albedo. Phases II and III show strikingly different evolutions. In phase II, the albedo decrease is mostly in the infrared. In this wavelength range, albedo is mostly sensitive to grain size in surface layers (Wiscombe and Warren, 1980; Domine et al., 2006), and the melting of these layers results in the decrease in SSA (i.e., increase in grain size), as shown in Fig. 6, which reduces the infrared albedo. In the visible, albedo is much less sensitive to grain size. The physical variable most likely to affect its value is the snowpack thickness. Wiscombe and Warren (1980) (their Fig. 13) show, using calculations, that for large snow grains (1 mm, fairly similar to melt grains we observed) the albedo at 500 nm decreases in a detectable manner when the snowpack water equivalent decreases to 100 mm. This corresponds to a depth of about 28 cm at our site, a value commonly observed (Fig. 6). Thus, snow depth during phase II mostly remains above the threshold where visible albedo would be affected.

In phase III on the other hand, the infrared albedo varies little, as visible in Fig. 3 for both 2015 and 2016 and in Fig. 7c. This is because SSA has essentially reached its terminal value of a few square meters per kilogram. In the visible however, the albedo is affected by the decrease in snowpack thickness. Radiation reaches the sea ice, where it is more absorbed than on snow, and the visible albedo decreases as the snowpack thins down. In phase IV, the albedo decreases in both the infrared and the visible because snow is replaced by ice and melt ponds. The trends in phase II and III were perturbed by snowfalls, but overall the dynamics of the spectral albedo appears to be a good indication of the phase of the snowpack evolution on sea ice.

The results concerning the overall evolution of the broadband albedo (Fig. 10) are in line with previous observations (Perovich et al., 2002; Grenfell and Perovich, 2004; Nicolaus et al., 2010; Gallet et al., 2017). However, our detailed measurements of snow physical properties allow a more detailed interpretation of spectral albedo evolution, in particular effects related to the decrease in SSA.

These results suggest that impurities can be important on the fast ice to predict melt dynamics and transmission through the column. Because LAPs tend to concentrate during melt and because of the enhanced impact of LAPs for low SSA values (Kokhanovsky, 2013), shortwave absorption can be greatly enhanced, further amplifying melt. LAPs also have a peculiar impact on the PAR available under the sea ice, as a very small number of impurities are sufficient to reduce by orders of magnitude the transmittance in the blue–green region (Tuzet et al., 2019). While this effect may concern a small fraction of the Arctic sea ice cover, it seems important enough to be further investigated.

4.3 Albedo modeling, current limitations

Figure 9 shows an overall very good ability of our model to simulate measured spectral albedo. Snow density and SSA used in modeling were measured values used without adjustment. For simplicity, single values for impurity concentrations (one for MD and one for BC) were used for the whole field campaign. Although these values were optimized, they reflect values measured in whalebacks in 2015 and in melted snow samples in 2016 within a factor of 2. This is acceptable given the uncertainties related to LAPs in general (Tuzet et al., 2019). Further improvements in modeling could be obtained by optimizing B and g. We used B and g values for spheres (Libois et al., 2013) because most of the time surface snow crystals were rounded by melting. Today, however, we have little understanding of how grain shapes affect B and g (Libois et al., 2014). Ideally, different values would be used for each snow layer. This clearly has the potential for improving the agreement between simulations and measurements, especially in the infrared. Attempts using different B and g values using TARTES indeed show that improvements are possible. However, since we currently would have little justification for selecting such values, we feel that this purely fitting exercise would be pointless. At present, our simulations are most of the time within 1 % of measurements in the visible and within 2 % in the infrared. Improving this agreement could be done by performing chemical measurements of impurities in each snow layer and by improving our understanding of the relationship between grain shape and the physical variables B and g. This latter suggestion could be realized by either performing field studies similar to Libois et al. (2014) or ray-tracing simulations on microtomography images of real snow samples (Haussener et al., 2012; Letcher et al., 2021).

4.4 Radiative fluxes and phytoplankton blooms

Figure 11 shows that the PAR flux increase at the snow–ice interface around 9 June coincides fairly well with the chlorophyll a increase. Massicotte et al. (2020) report time series of average daily PAR values in the seawater at 1.3 m depth (also shown in Fig. 11) and indicate that the initiation of the phytoplankton bloom was around 8 June. The PAR under-ice increase curve is somewhat different from our modeled curve, but this can probably be explained by different snowmelt phenology at our snow measurement sites, which are different for each day, and at the PAR measurement sites. The coincidence between the increase in snow transmittance and the increase in chlorophyll a is nevertheless pretty good.

5 Conclusion

Snow over sea ice was intensively studied during two melt seasons in Baffin Bay. These studies include spectral albedo measurements and vertical profiles of physical properties of snow that are relevant to radiative transfer modeling: density and SSA. The entire transition from cold and dry snow covers in early spring to ponded sea ice in late spring was recorded.

In both years, albedo evolved following four main phases related to the conditions of the snow cover. During these phases, broadband albedo was first high, up to 0.87, over a dry snowpack (phase I) composed of basal layers of indurated depth hoar and faceted crystals topped by one or several wind slabs with sometimes a layer of fresh snow. Albedo gradually decreased in the near-infrared as snow grains coarsened because of surface melting (phase II). At some point, the snowpack became ripe and isothermal, and its thickness decreased faster, leading to a decrease in albedo in the visible range for the first time (phase III). This drop in albedo was due to the influence of the underlying darker sea ice as the light penetration depth in snow increased, and snow depth decreased. Spatial variability appeared and was directly linked to the snow thickness and the optical properties of the underlying media. Melt ponds formed with snow meltout during phase IV.

Snow physical properties and impurity contents were used as inputs to a radiative transfer model in order to simulate the albedo. Simulations reproduced the measured time series of albedo very well and demonstrated the interest of measuring time series of snow SSA vertical profiles to understand albedo evolution. Other important factors identified here are the impurity content and snow depth, which have a strong impact in the visible range, and the geometric factors B and g, which are today poorly constrained for real snow.

Lastly, this work allowed the calculation of the spectral radiative fluxes at the snow–ice interface. Calculations showed that PAR fluxes at this interface rose durably above 9 µmol m⁻² s⁻¹ only at the end of phase II (11 June), when the snowpack thickness started decreasing. The timing of this PAR increase coincided fairly well with the initiation of the phytoplankton bloom, as seen in Fig. 11, and with the date mentioned by Massicotte et al. (2020), also 11 June in their Fig. 5. This suggests that modeling radiative transfer through snow, based on accurate data on snow physical and chemical properties, may allow the reliable prediction of the timing of phytoplankton blooms.