Exploring the role of snow metamorphism on the isotopic composition of the surface snow at EastGRIP

Stable water isotopes from polar ice cores are invaluable high-resolution climate proxy records. Recent studies have aimed to improve our understanding of how the climate signal is stored in the stable water isotope record by addressing the influence of post-depositional processes on the isotopic composition of surface snow. In this study, the relationship between surface snow metamorphism and water isotopes during precipitation-free periods is explored using measurements of snow-specific surface area (SSA). Continuous daily SSA measurements from the East Greenland Ice Core Project site (EastGRIP) during the summer seasons of 2017, 2018 and 2019 are used to develop an empirical decay model to describe events of rapid decrease in SSA linked to snow metamorphism. We find that SSA decay during precipitation-free periods at the EastGRIP site is best described by the exponential equation $\mathrm{SSA}\left(t\right)=({\mathrm{SSA}}_{\mathrm{0}}-\mathrm{22})\cdot e^{-\mathit{\alpha }t}+\mathrm{22}$, and has a dependency on wind speed. The relationship between surface snow SSA and snow isotopic composition is primarily explored using empirical orthogonal function analysis. A coherence between SSA and deuterium excess is apparent during 2017 and 2019, suggesting that processes driving change in SSA also influence snow deuterium excess. By contrast, 2018 was characterised by a covariance between SSA and δ¹⁸O highlighting the inter-annual variability in surface regimes. Moreover, we observed changes in isotopic composition consistent with fractionation effects associated with sublimation and vapour diffusion during periods of rapid decrease in SSA. Our findings support recent studies which provide evidence of isotopic fractionation during sublimation, and show that snow deuterium excess is modified during snow metamorphism.

1 Introduction

The traditional interpretation of stable water isotopes in ice cores is based on the linear relationship between local temperature and first-order parameters δ¹⁸O and δD of surface snow on ice sheets (Dansgaard, 1964). Accurate reconstruction requires consideration of precipitation intermittency (Casado et al., 2020; Laepple et al., 2018), past variations in ice-sheet elevation (Vinther et al., 2009), sea ice extent (Faber et al., 2017; Sime et al., 2013), and firn diffusion (Johnsen et al., 2000; Landais et al., 2006; Holme et al., 2018), which influence the water isotopic composition in ice cores. The second-order parameter deuterium excess (d-excess) is defined by the deviation from the near-linear relationship between δ¹⁸O and δD which is driven by non-equilibrium (kinetic) fractionation (d-excess = δD $-\mathrm{8}\cdot {\mathit{\delta }}^{\mathrm{18}}$O). d-excess in ice cores is understood to reflect moisture source conditions (Dansgaard, 1964; Merlivat and Jouzel, 1979; Johnsen et al., 1989) and changes in moisture source region (Masson-Delmotte et al., 2005), and can be modified during snow crystal formation in supersaturated clouds (Ciais and Jouzel, 1994; Sodemann et al., 2008). Recent studies have documented isotopic composition change in the surface snow during precipitation-free periods (Steen-Larsen et al., 2014; Ritter et al., 2016; Casado et al., 2018; Hughes et al., 2021), linked to synoptic variations in atmospheric water vapour composition and subsequent exchange with the surface snow (Steen-Larsen et al., 2014; Ritter et al., 2016; Madsen et al., 2019; Hughes et al., 2021; Wahl et al., 2021; Casado et al., 2021). Post-depositional processes at the surface involve additional kinetic effects adding complexity to the interpretation of d-excess (Hughes et al., 2021; Casado et al., 2021). Here we focus on processes influencing the isotopic composition of the surface snow after deposition while exposed to surface processes, and concentrate on the second-order parameter d-excess due to its sensitivity to kinetic effects.

After deposition, snow grains undergo structural changes known as “snow metamorphism”, which is active at the surface and at greater depths, depending on temperature (gradient) conditions (Colbeck, 1983; Pinzer and Schneebeli, 2009). Surface snow metamorphism is initially driven by a reduction in the snow–air interface to reach thermodynamic stability (Colbeck, 1980; Legagneux and Domine, 2005). The snow–air interface can be described by the widely used parameter snow-specific surface area (SSA), which is dependent on optical grain radius and density of ice (SSA = $\mathrm{3}/{\mathit{\rho }}_{\mathrm{ice}}\cdot r_{\mathrm{opt}}$) (Gallet et al., 2009) and can be used as an indicator for snow metamorphism (Cabanes et al., 2002, 2003; Legagneux et al., 2002). Freshly deposited snow has a high SSA which decreases with time under both isothermal (<10 ^∘C m⁻¹) and temperature gradient (>10 ^∘C m⁻¹) conditions within the snow (Cabanes et al., 2002; Legagneux et al., 2004; Domine et al., 2007; Genthon et al., 2017). A decrease in SSA in dry snow is predominantly the result of water vapour transfer among grains, with smaller grains feeding the growth of larger grains (Legagneux et al., 2004; Flin and Brzoska, 2008; Sokratov and Golubev, 2009; Pinzer et al., 2012). Ventilation by wind can accelerate SSA decrease by enhancing water vapour transfer rates (Picard et al., 2019). Under natural conditions, SSA decrease is driven by a combination of these processes depending on surface conditions (Cabanes et al., 2003; Pinzer and Schneebeli, 2009), each potentially modifying the isotopic composition of the snow (Ebner et al., 2017).

Models can provide a quantitative description of SSA decrease after deposition. Previous studies have proposed SSA decay models using a combination of field measurements and controlled laboratory experiments (Cabanes et al., 2002, 2003; Legagneux et al., 2003, 2004; Flanner and Zender, 2006; Taillandier et al., 2007). Exponential models are documented to produce the best fit to in situ SSA decay data given that they can account for temperature effects under various wind conditions while being simple in their formulation (Cabanes et al., 2003). A subsequent physical-based model was defined by Legagneux et al. (2004) to describe SSA decay based on grain growth theory, which was then further developed by Flanner and Zender (2006), who defined parameters based on surface temperature, temperature gradient and snow density. Existing SSA decay models have rarely been applied to polar ice sheet surface snow (Linow et al., 2012; Carmagnola et al., 2013). Conditions for surface snow on polar ice sheets are not necessarily comparable to other alpine and arctic regions due to negligible melt and the high-latitude radiation budget. Moreover, while continuous surface SSA measurements exist from Antarctica (Gallet et al., 2011, 2014; Picard et al., 2014), those from Greenland focus on the depth evolution of SSA (Linow et al., 2012; Carmagnola et al., 2013). Continuous datasets of daily SSA and corresponding isotopic composition measurements from the accumulation zone of the Greenland Ice Sheet are required for understanding the influence of surface snow metamorphism on surface energy budget (Picard et al., 2012), and for the interpretation of ice core water isotope records (Casado et al., 2021; Wahl et al., 2022). In this study we focus on the latter, which is of particular importance owing to observations of isotopic fractionation during snow metamorphism documented in laboratory studies (Ebner et al., 2017) and field experiments (Casado et al., 2021; Hughes et al., 2021; Wahl et al., 2021). Nonetheless, few studies have focused on the direct relationship between physical snow properties, such as SSA, and post-depositional changes in isotopic composition.

In this paper, the aim is to explore the behaviour of surface snow metamorphism on polar ice sheets using daily SSA measurements from northeast Greenland during summer, and to compare the change in physical properties with the isotopic composition measurements. The primary focus is to document events where changes in SSA occur rapidly over a number of precipitation-free days, which we use as a proxy for snow metamorphism. Events of rapid SSA decrease (SSA decay events) are used to (1) quantify and model surface snow metamorphism in polar snow and (2) assess isotopic change during surface snow metamorphism in situ. The data presented here have great value for our understanding of the influence of post-depositional processes on physical and isotopic changes in the polar ice sheet surface snow. This allows for a deeper understanding of snow properties at remote regions of polar ice sheets and contributes to the interpretation of stable water isotopes in polar ice cores.

2 Study site and methods

2.1 EastGRIP site overview and meteorological data

All data used in this paper were collected as part of the Surface Program corresponding to the international deep ice core drilling project at the East Greenland Ice Core Project site (EastGRIP, 75.65^∘ N, 35.99^∘ W; 2,700 m a.s.l.) during summer field seasons (May–August) of 2017, 2018 and 2019. The local accumulation rate is approximately 12 ± 2 cm w.e. yr⁻¹ (Schaller et al., 2017).

Figure 1 (a) A map of Greenland with a green dot indicating the EastGRIP site (Wahl et al., 2021). (b) A photograph of the clean snow area at the field site (Credit: Bruce Vaughn), black lines indicate the SSA sampling transect with 10 m spacing shown as dashed lines. (c) A photograph of SSA sampling cups (Credit: Sonja Wahl), and (d) an illustration of the sampling device from Klein (2014).

Meteorological data used for this study are from the Program for Monitoring of the Greenland Ice Sheet (PROMICE) Automatic Weather Station set up by the Geological Survey of Denmark and Greenland (GEUS) at EastGRIP in 2016 (Fausto et al., 2021). The data are 10 min mean values for a number of variables. In addition to the surface variables, snow temperature was measured using a thermistor string at 0.1 m intervals during 2017 and 2018 but was modified to 1 m intervals in 2019. An additional thermistor string was thus installed in May 2019, from which we use the 0.1 m measurements. Relevant weather variables for this study are surface temperature (calculated from longwave radiation down and longwave radiation up with longwave emissivity set at 0.97), air temperature and wind speed (Van As, 2011). Mean weather conditions vary between sampling years, as outlined in Table 1, with prevailing westerly winds throughout the sampling seasons.

An eddy covariance (EC) measurement tower was set up at EastGRIP for every summer observation period (Madsen et al., 2019; Wahl et al., 2021). Here we use the 30 min latent heat flux (LE) measurements which are calculated from the measurement of humidity fluxes between the surface and atmosphere. Positive LE indicates upwards energy flux in the form of sublimation in Table 1. Observations of weather conditions such as ground fog, drifting snow and snowfall were documented each day in the EastGRIP field diary.

Table 1 Weather statistics for summer field campaigns in 2017, 2018 and 2019. Mean and standard deviation for weather variables during the three sampling seasons. Surface temperature, relative humidity with reference to ice and wind speed use PROMICE weather station based on 10 min measurements. Latent heat flux, a 30 min mean upwards flux from the eddy covariance tower.

2.2 Accumulation, snow sampling and subsequent measurements

Each summer season of 2017, 2018 and 2019, snow samples were taken once a day, primarily in the morning, from May to August at 10 sampling sites. Each site was marked by a stick, along a 90 m transect running perpendicular to the prevailing wind direction (285^∘ N (Wahl et al., 2021)) with 10 m spacing upwind of the EastGRIP camp to ensure clean snow (Fig. 1b). The specific dates for each season are given in Table 1. The precise location of each sample was marked by a small stick to ensure the adjacent snow is sampled the next day and to avoid sampling snow from different depths. A 6 cm diameter sampling device (cup) collected the top 2.5 cm of surface snow (Fig. 1c and d). Snow density is determined using the weight of each snow sample with a known volume. Sticks were placed at each sampling site at the start of each season to measure snow height, that is, the distance between the snow surface and top of the stick with a ruler with an uncertainty margin of ±0.5 cm. Accumulation, used here to describe the change in snow height (cm), was calculated using the cumulative sum of the daily difference between measurements of snow height from each site. The resultant datasets consist of daily measurements of four parameters at each of the 10 sampling sites: SSA, density, snow accumulation and stable water isotopes. The field season for 2018 started on 5 May, 9 d earlier than 2017 (14 May), and 22 d earlier than 2019 (27 May). The meteorological data are re-sampled to the SSA sampling time periods.

2.3 SSA measurements

Each snow sample is placed into the Ice Cube sampling container below an infra-red (IR) laser diode (wavelength of 1310 nm), where the SSA is calculated based on IR hemispherical reflectance, as explained in Gallet et al. (2009). More information on the Ice Cube device can be found in Zuanon (2013). Gallet et al. (2009) show that SSA measurements for snow with density of 200 kg m⁻³ and SSA of 35 m² kg⁻¹ mostly reflect the top 1 cm of a 2.5 cm snow sample when using 1310 nm radiation, due to the e-folding depth. The properties of each 2.5 cm snow sample will determine the e-folding depth (i.e. the depth to which the light irradiances within the snowpack is reduced to $\mathrm{1}/e$ (approximately 37 %) of its initial value), with higher SSA and density causing a decreased e-folding depth. Given that the mean snow density from all field seasons is 293 kg m⁻³ (307 ± 40, 278 ± 47, 294 ± 50 kg m⁻³ for 2017, 2018 and 2019 respectively) and the mean SSA is 37.5 m² kg⁻¹, the measurement of the SSA values will be weighted towards the top <1 cm of the 2.5 cm sample. However, recent studies have shown that the SSA values obtained from the instrument used here (Ice Cube) agree well with measurements from computed microtomography on the same samples (Martin and Schneebeli, 2022). Further, our data are in fairly good agreement with remote sensing products (Kokhanovsky et al., 2019; Vandecrux et al., 2022). Thus we consider our SSA values as representative for the upper 2.5 cm surface snow.

The light reflected from the snow samples is converted into inter-hemispheric IR reflectance using a calibration curve based on methane absorption methods (Gallet et al., 2009). A radiative-transfer model is used to retrieve SSA from inter-hemispherical IR reflectance. To avoid influence from solar radiation, SSA was measured inside a white tent or in a snow cave kept at temperatures between −5 and −20 ^∘C. We assume an uncertainty of 10 % for SSA measurements between 5 and 130 m² kg⁻¹ (Gallet et al., 2009).

2.4 Surface snow isotopes

Individual SSA samples were put in separate bags and subsequently measured for water isotopic composition. Thus, every day the 10 SSA samples have a corresponding isotopic composition. The resultant isotope value is the average composition over the top 2.5 cm of snow. Each sample was kept frozen during transportation and storage. The samples were then analysed at Alfred Wegener Institute in Bremerhaven using cavity ring-down spectroscopy instruments (Picarro L-2120-i and L-2140-i) following the protocol of Van Geldern and Barth (2012). This technique produces δ¹⁸O and δD measurements with an estimated uncertainty of 0.15 ‰ and 0.8 ‰ respectively. The values calculated for d-excess have an estimated uncertainty of 1 ‰.

2.5 Data analysis

2.5.1 Defining SSA decay events

Freshly deposited snow has a high SSA that slowly decreases through time due to snow metamorphism. Based on this understanding, two terms are defined:

1. SSA increase: Increases in SSA indicate deposition events in the form of precipitation or deposited snowdrift.

2. SSA decrease: Decreases in SSA are due to snow metamorphism and other post-depositional processes such as wind scouring and, in a few cases, surface hoar formation, where the SSA decreases.

SSA decays are primarily identified in the time series to quantify snow metamorphism and corresponding isotopic change during precipitation-free periods. Surface conditions are considered to remove events with surface snow layer perturbations via deposition or erosion by the wind.

A threshold is derived to systematically identify periods of rapid SSA decay – hereafter referred to as “SSA decay events”. SSA decay events captured by this threshold are defined by the peak SSA value (Day-0), through to the next increase in SSA. A set of criteria are applied to the SSA decay events to avoid events with wind-perturbed surfaces. While in Antarctica, drifting of unconsolidated snow has been observed at mean hourly wind speeds as low as 4.5 m s⁻¹ at 2 m height above the surface (Birnbaum et al., 2010), a study from northeast Greenland, with similar conditions to EastGRIP, documented snowdrift starting at 6 m s⁻¹ (Christiansen, 2001), due to warmer temperatures facilitating bonding of the surface snow (Li and Pomeroy, 1997). Additional field-diary observations from EastGRIP document significant snowdrift when wind speeds exceed 7 m s⁻¹. Based on these observations, two wind speed categories are defined.

1. Low-wind events: The first includes events with daily maximum wind speed (computed from 10 min averaged wind speed) consistently below 6 m s⁻¹, hereafter referred to as “low-wind events”, where negligible surface perturbation is ensured.

2. Moderate-wind events: A second category considers events with daily maximum wind speed between 6 and 7 m s⁻¹, hereafter “moderate-wind events”. The inclusion of these events facilitates an assessment of the influence of wind speed on SSA decay.

Subsequent isotopic analysis is first broadly applied to both low- and moderate-wind events over 1 and 2 d periods, followed by a focused assessment of isotopic change, and corresponding temperature flux is applied to low-wind events alone given the assurance of unperturbed snow. All events with wind speed above 7 m s⁻¹ are excluded from analysis.

2.5.2 Modelling surface snow metamorphism

The first empirical SSA decay model, Eq. (1), was proposed by Cabanes et al. (2003), who described a temperature-dependent exponential decay based on snow samples collected from the Alps (Cabanes et al., 2003) and the Canadian Arctic (Cabanes et al., 2002). Legagneux et al. (2003) found Eq. (2) to best describe experimental SSA decay under controlled laboratory conditions.

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{SSA}\left(t\right)={\mathrm{SSA}}_{\mathrm{0}}\cdot e^{-\mathit{\alpha }t}\text{(2)}& {\displaystyle }\mathrm{SSA}\left(t\right)=B-A\cdot \ln (t+\mathrm{\Delta }t)\end{array}$$

In both Eqs. (1) and (2), SSA(t) is the SSA value at a time, t, in days since the initial SSA value, SSA₀, and α is the decay rate. Parameters A and B in Eq. (2) were found to be arbitrarily related to the decay rate and initial SSA of each sample. To improve the physical basis of the model, the theory of Ostwald ripening, describing grain growth driven by thermodynamic instability, was implemented into the model (Legagneux et al., 2004). Equation (3) has two parameters, τ and n; τ is the decay rate and n relates to theoretical grain growth. The physical model was further developed by Flanner and Zender (2006) to incorporate a physical quantification of the parameters including information about temperature, temperature gradient, and density. Based on these three conditions, they created a look-up table for τ and n.

$$\begin{array}{}\text{(3)}& \mathrm{SSA}\left(t\right)={\mathrm{SSA}}_{\mathrm{0}}\cdot {\left({\displaystyle \frac{\mathit{\tau }}{t+\mathit{\tau }}}\right)}^{\mathrm{1}/n}\end{array}$$

Taillandier et al. (2007) proposed two equations based on Eq. (2) to define the decay rate under isothermal and temperature gradient conditions where they were able to directly incorporate a mean surface temperature parameter (T_m). Here we use the model for temperature gradient conditions (Eq. 9 in Taillandier et al., 2007).

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathrm{SSA}\left(t\right)& =[\mathrm{0.659}\cdot {\mathrm{SSA}}_{\mathrm{0}}-\mathrm{27.2}\cdot (T_{\mathrm{m}}-\mathrm{2.03}\left)\right]\\ & -[\mathrm{0.0961}\cdot {\mathrm{SSA}}_{\mathrm{0}}-\mathrm{3.44}\cdot (T_{\mathrm{m}}+\mathrm{1.90}\left)\right]\\ & \cdot \ln \mathit{\{}t+e^{\left(\frac{-\mathrm{0.341}\cdot {\mathrm{SSA}}_{\mathrm{0}}-\mathrm{27.2}\cdot (T_{\mathrm{m}}-\mathrm{2.03})}{\mathrm{0.0961}\cdot {\mathrm{SSA}}_{\mathrm{0}}-\mathrm{3.44}\cdot (T_{\mathrm{m}}+\mathrm{1.90})}\right)}\mathit{\}}\end{array}\end{array}$$

Building upon these previous studies, we define an empirical SSA decay model using continuous daily SSA measurements from EastGRIP to describe the behaviour of surface snow SSA in polar summer conditions (Cabanes et al., 2002, 2003; Flanner and Zender, 2006; Legagneux et al., 2002, 2003; Taillandier et al., 2007). A recent study at EastGRIP has shown the significant heterogeneity in surface snow due to post-depositional reworking from the wind (Zuhr et al., 2021), and therefore each sample location is treated individually to avoid the smoothing out of localised signals when averaging.

3 Results

3.1 EastGRIP meteorological conditions

Meteorological variables over the three sampling seasons vary substantially, as shown in Fig. 2. Air temperatures were below −30 ^∘C between 5–8 May in 2018. Such low temperatures were not recorded for 2017 and 2019. Moreover, when comparing the period from 27 May with 5 August of each year (duration of 2019 season), 2018 air temperatures (−13.3 ^∘C) were still 0.5 ^∘C lower than 2017 and 3.2 ^∘C lower than 2019 (Table 1).

Figure 2 Time series of the ambient conditions during the sampling periods. Data cover the specific sampling periods each year. The 10 min mean data from the AWS are presented for air temperature, wind speed and relative humidity. Mean values (bold lines) are shown for air temperature and relative humidity, while daily maximum is shown for wind speed. Latent heat flux data are 30 min averages from the eddy covariance tower, with the daily sum shown in bold. Snow accumulation is presented in the lower panel as the daily mean over the 10 sampling sites (see Sect. 2.2). Grey bars indicate the derived SSA decay events.

The 2017 season was characterised by high wind intrusions of >13 m s⁻¹ at approximately 20 d intervals. Considering all three sampling years, the average daily maximum wind speed is 7 m s⁻¹, with 209 out of the total 237 sampling days having maximum wind speed above 5 m s⁻¹. The distributions of daily maximum wind speed compared to 10 min mean values are found in the Appendix Fig. A1. Net accumulation was observed during all seasons, with 4.8, 9.3 and 8.6 cm accumulated for the 2017, 2018 and 2019 field seasons respectively. In addition, the total sum of LE during 2019 was 30 % greater than in 2018 indicating strong sublimation. Eddy covariance LE measurements are supported by LE from AWS observations which indicate the same magnitude of difference between the years.

3.1.1 Spatial and temporal snow surface variability

Spatial and temporal variability – defined as the daily standard deviation over transect and the standard deviation of the daily mean values over the season respectively – is observed in SSA, δ¹⁸O and d-excess throughout the field seasons of 2017, 2018 and 2019 (Fig. 3), with highest spatial variability in isotopic composition. SSA is characterised by peaks, often corresponding to high spatial variability, followed by gradual decreases over a number of days, the SSA decay feature which is most prominent during 2017 and corresponds to negligible or decreasing accumulation. High SSA values at the start of the 2018 season (daily mean 88 m² kg⁻¹) correspond to low and spatially homogeneous δ¹⁸O.

Figure 3 Time series of SSA (blue), δ¹⁸O (orange), d-excess (green) and the first principal components (PC1) of each variable in the lower panel. Bands in the upper three panels show spatial standard deviation between the 10 sampling sites. The secondary y axis in the top panel shows the average snow accumulation. The grey bars indicate the SSA decay events.

Inter-annual difference is observed in δ¹⁸O, with seasonal mean values of −31.6 ± 2.1 ‰, −32.7 ± 1.3 ‰ and −27.3 ± 2.1 ‰ for 2017, 2018 and 2019 respectively (Fig. 3). Throughout the season δ¹⁸O follows a gradual increasing trend from May to August concurrent with increasing temperatures. Cases of abrupt decreases (−10 ‰) are observed in the late summer, for example, on 12 July in 2018 and 25 July in 2019, originating from late-summer snowfall events. Note that the 2019 field season started approximately 15 d later than 2017 and 2018, resulting in a bias towards mid-summer conditions. No clear seasonal trend is observed in d-excess, but rather there are periods of gradual decrease in d-excess during periods with no accumulation in all years. The most apparent is from 15 May to 14 June in 2017 corresponding to 0 cm net-accumulation (Fig. 3). The maximum daily spread in δ¹⁸O and d-excess is approximately 15 ‰, indicating strong surface heterogeneity.

Empirical orthogonal function (EOF) analysis is applied to the data to identify the dominant modes of variance in both the temporal and spatial dimensions for each parameter – SSA, δ¹⁸O and d-excess. Using a confidence interval of 95 % (p<0.05), the relationship between SSA and isotopic composition is tested for all 10 identical samples, covering the entire measurement period. The spatial and temporal principal components returned for each variable by the EOF are presented in Fig. 3. During 2017, 2018 and 2019 all variables have one dominant mode of variance, or principal component (PC1). PC1 of SSA (PC1_SSA) explains 61 %, 77 % and 72 % of the variance for the respective years, PC1 of δ¹⁸O (PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$) explains 69 %, 83 % and 75 % of the total variance respectively, while PC1 of d-excess (PC1_dxs) explains 47 %, 51 % and 60 %.

Results from the EOF analysis reveal distinct differences between the sampling years, most prevalent is the opposing regime from 2018 to 2019. During 2018, PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ and PC1_dxs have an inverse correlation in the spatial dimension ($r=-\mathrm{0.6}$), and a significant positive correlation between PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ and PC1_dssa in the temporal dimension (p<0.05, r=0.5) (Figs. A2 and A3). By contrast, data from 2019 are characterised by significant positive correlations between PC1_SSA and PC1_dxs in both the spatial (p<0.05, r=0.75) and temporal dimensions (p>0.05, r=0.3), while no relationship is observed between PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ and PC1_dxs. During 2017, there is a weak positive correlation for the temporal component of PC1_SSA and PC1_dxs (p<0.05, r=0.3), and the temporal PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ and PC1_dxs (p<0.05, r=0.4). There is an apparent shift after 15 July where PC1_dxs transitions from co-varying with PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ to co-varying with PC1_SSA. Figures A2 and A3 in the Appendix illustrate the spatial and temporal components of the EOF results.

3.2 SSA decay events

3.2.1 Observations

In agreement with Eq. (1), we observe a relationship between the initial SSA and the subsequent magnitude of decrease when evaluating the SSA decay events (Fig. 3). From the years 2017, 2018 and 2019 a total of 21 events are identified that fulfil the SSA decay criteria (as defined in Sect. 2.5.1). These events are named E1, E2 etc. (Table A1). Using the thresholds defined in Sect. 2.5.1, 15 out of the 21 events are excluded from the analysis due to risk of surface perturbations. Of the six events, two are in the low-wind category (E10 and E11), and four in the moderate-wind category (E1, E13, E18 and E19). Both E10 and E11 had consistent clear sky conditions. Note that E11 was preceded by ground fog, and not snowfall (Table A1).

The rate of SSA decay is strongly influenced by the initial SSA during the decay period (Fig. 4), showing that the rate of change is proportional to the absolute value, as described by an exponential decay law ($r=-\mathrm{0.71}$ and $r=-\mathrm{0.91}$ for low- and moderate-wind events respectively). The mean air temperature for all SSA decay events was between −17.3 and −7 ^∘C. The first day of each event is characterised by the largest change in SSA, followed by a decrease in magnitude over the subsequent days, with negligible change in SSA below 22 m² kg⁻¹.

Figure 4 Linear regressions for change in SSA against the SSA for the low-wind (light blue) and moderate-wind (dark blue) SSA decay events (a) considering all individual samples. The observed SSA decays are shown for the moderate-wind events (b), and the low-wind events (c), followed by the modelled SSA decays for the respective events in (d) and (e). The legend in (d) and (e) corresponds to the SSA decay event number in Table A1.

3.2.2 Model

The SSA decay rate is quantified by plotting the rate of change in SSA per day against the absolute SSA value for all 10 sampling sites for low- (E10 and E11) and moderate-wind (E1, E13, E18 and E19) events (Fig. 4a). We observe a linear relationship between the change in SSA from one day to the next (ΔSSA) and SSA.

The SSA decay model for EastGRIP is constructed using the differential equation for the linear relationship between ΔSSA and absolute SSA. Solving the differential equation with respect to time (t) produces the SSA decay model defined as Eq. (5), which follows the equation structure of Eq. (1):

$$\begin{array}{}\text{(5)}& \mathrm{SSA}\left(t\right)=({\mathrm{SSA}}_{\mathrm{0}}-\mathrm{22})\cdot e^{-\mathit{\alpha }\cdot t}+\mathrm{22},\end{array}$$

where SSA(t) is the SSA measurement at a given time in days since the first measurement (initial SSA), SSA₀ is the initial SSA value, and α is the decay rate. An offset of 22 m² kg⁻¹ is required to account for the non-zero asymptote, and is defined as the SSA value where the derivative of SSA is equal to 0 m² kg⁻¹. The decay rate, determined by the slope of the linear regressions in Fig. 4, is higher for moderate-wind SSA decay events (−0.66 m² kg⁻¹ d⁻¹) than for low-wind SSA decay events (−0.41 m² kg⁻¹ d⁻¹). Within the temperature range of low- and moderate-wind SSA decay events from this study, there is no observable temperature dependence of the SSA decay rate. However, such a temperature dependence is clear in events which were excluded due to potential surface perturbations, where there is a decrease in decay rate and an increase in the background SSA value with low temperatures (Fig. A4). Our results indicate a slower rate of decay under decreased wind speed conditions.

3.2.3 Model performance comparison

Model performance is tested by comparing daily predicted decrease with the 10 daily observations, and by comparing model predictions from this study with those from Flanner and Zender (2006), hereafter referred to as “FZ06”, and the model defined by Taillandier et al. (2007), hereafter “T07”, as defined in Sect. 2.5.2. Residuals between our model and the observations are normally distributed, suggesting no systematic errors in model predictions. The RMSE between our model predictions and observed SSA is 2.48 and 2.60 m² kg⁻¹ for the low-wind and moderate-wind SSA decay events respectively.

Table 2 RMSE values for model evaluation. Model predictions are compared with observations for both the daily mean SSA values (Mean) and the 10 individual SSA measurements (Individual). For the low-wind events, α in Eq. (5) is −0.41 m² kg⁻¹ d⁻¹. For moderate-wind events α in Eq. (5) is −0.66 m² kg⁻¹ d⁻¹. FZ06 parameters τ and n are defined from the look-up table in Flanner and Zender (2006). T07 uses the mean surface temperature for each event as an input parameter.

Parameter values for τ and n in FZ06 are defined for each event based on mean density, surface temperature and temperature gradient using the extensive look-up table referenced in Flanner and Zender (2006). FZ06 consistently overestimates the SSA decay rate, with residuals increasing throughout the events (Fig. A5). T07 generally underestimates the SSA decay rate of low- and moderate-wind events, with the largest errors associated with E1, the event featuring the lowest temperatures and highest wind speeds. RMSE values presented in Table 2 show that the model from this study predicts the observed SSA decay with the least error.

3.3 Isotopic change during SSA decay events

3.3.1 Low- and moderate-wind event analysis

The following sections assess the snow isotopic change during the low- and moderate-wind events. For both low- and moderate-wind events, the change in d-excess is plotted against the change in SSA (Fig. 5), after the first and second day of each event. Here, we include the analysis of the change in d-excess after the second day (referred to as “a 2 d change”), presented in Table 3 for each low- and moderate-wind event.

Figure 5 Isotopic changes during all of the analysed events are shown, with each point indicating a specific sampling site. The daily change in d-excess (d_xs) and SSA is presented in (a), change in d-excess and SSA over the first 2 d of each event is shown in (b), and the respective changes in d-excess are plotted against the absolute d-excess values in (c). Linear regressions are included for daily change (orange) and 2 d change (grey).

Table 3 Table of isotopic change for decay events. Mean values on Day-0 and Day-2 of each event, and percentage change over this 2 d period are presented for δ¹⁸O, d-excess and SSA. Low-wind events (bold text) and moderate-wind events are presented.

Both δ¹⁸O and d-excess change from the initial value (Day-0) in all events, with the percentage change in d-excess being the same order of magnitude as for SSA, and an order of magnitude higher than that of δ¹⁸O. Three out of six events are characterised by increasing δ¹⁸O and decreasing d-excess after 2 d. E1, E13 and E19 deviate from this pattern. E13 and E19 both correspond to total increase in δ¹⁸O and d-excess, whereas E11 is characterised by a slight decrease in δ¹⁸O and decrease in d-excess (Table 3). The Δd-excess over a 1 d period indicates a slight negative skew around a mean of −0.3 ‰. Four of six events (61 % of all sample sites) experience a decrease in d-excess by Day-2. The mean is shifted to −1.2 ‰ over a 2 d period with decreases in d-excess documented in 45 out of the 60 samples (75 %) during precipitation-free periods with minimal surface perturbation. Initial d-excess is observed to have a significant influence on the magnitude of d-excess decrease over the defined period (Fig. 5c), with high initial d-excess corresponding to the largest decreases in d-excess.

3.3.2 Low-wind event analysis

The following section focuses on latent heat flux (LE) and near-surface temperature gradients (TG) corresponding to isotopic change during low-wind events only. This is to ensure that the surface layer we are analysing is constant throughout the event to avoid inaccurate interpretation of isotopic change. As mentioned in Sect. 3.2, ground fog preceded the SSA peak in E11, concurrent with negligible accumulation recorded. By contrast, the day prior to E10 corresponds to the observation of snowfall (Table A1). Figure 6 shows the relationship between the daily change in isotopic composition for E10 and E11 (Δd-excess and Δδ¹⁸O) and SSA (ΔSSA). During E10, significant inverse correlations are observed between Δδ¹⁸O and Δd-excess, and between ΔSSA and Δd-excess ($r=-\mathrm{0.5}$, $r=-\mathrm{0.8}$), while Δδ¹⁸O and ΔSSA are positively correlated (r=0.6). By contrast, no significant relationship is observed between the Δ-parameters during E11, where there is negligible change in surface snow δ¹⁸O (<0.7 ‰).

Figure 6 Isotopic change analysis for low-wind events, E10 and E11. Panel (a) shows daily change in d-excess against change in δ¹⁸O for E10 and E11, (b) shows change in d-excess against change in SSA, and (c) shows change in δ¹⁸O and change in SSA. The r and p values for each regression are indicated in the corresponding colours. Only significant linear regressions are indicated with a line.

Figure 7 LE and temperature gradients for low-wind events. Latent heat flux (LE) (grey), air-surface temperature gradient (TG) (bright red) and surface-10 cm subsurface TG (dark red) during June 2018. E10 (blue) and E11 (red) are highlighted. Dark grey shading in LE indicates sublimation and light grey shows deposition.

The direction of vapour fluxes is inferred using temperature gradients determined from air, surface and subsurface (10 cm depth) temperature data, and LE, measured as an upwards flux. Net sublimation is observed during both E10 and E11, with a total sum of 33.9 and 55.8 W m⁻² for the respective events. LE is inversely related to the TG between the air and the surface, with strong sublimation (>10 W m⁻²), corresponding to a negative TG of 2.5 ^∘C between the air and surface on 10 June 2018 (E10). A concurrent upwards vapour flux is observed between the subsurface and surface, most apparent on 11 June 2018 (Day-2 of E10). Negative LEs up to −4 W m⁻¹ are documented each night corresponding to the transition from a negative to positive TG between the air and surface. Net-deposition was recorded between sampling on 9 June at 15:30 UTC and 10 June 10:30 UTC 2018 (first day of E10) corresponding to a decrease in δ¹⁸O and d-excess. The subsequent day, characterised by net sublimation, had a large increase in δ¹⁸O and a larger decrease in d-excess. Contrary to E10, E11 is characterised by a large decrease in d-excess and a small decrease in δ¹⁸O, both concurrent with net-sublimation and strong negative surface-subsurface TG. The air-surface TG during E11 has a lower mean and reduced diurnal amplitude than E10 facilitating sublimation for a longer period.

4 Discussion

The new SSA and snow isotopic composition datasets presented in this study have revealed concurrent decreases in d-excess (and δ¹⁸O to a lesser degree) and SSA during precipitation-free periods with minimal snow drift. A simple empirical model describing the SSA decay rate under different wind regimes reveals more rapid SSA decay when wind speeds are higher. The following sections look first at why existing models tend to be inaccurate when predicting in situ SSA decay at the surface, followed by a second section discussing the possible mechanisms driving the relationship between SSA and d-excess during precipitation-free periods.

4.1 SSA decay at EastGRIP

Events of rapid SSA decay at EastGRIP are best described by an exponential decay function, in agreement with observations from Cabanes et al. (2003). No cold events ($<-\mathrm{20}$ ^∘C) are captured within the event criteria, potentially indicating that warmer temperatures are needed to induce such change in SSA over the time period, or that the local synopticity favours precipitation (bringing snow with high SSA) coincident with low temperature and high winds during the summer. Events captured by the threshold with mean temperatures below −20 ^∘C are likely capturing wind erosion that exposes sub-surface snow with low SSA. The narrow temperature range of SSA decay events does not facilitate a conclusive definition of a temperature-dependent decay rate (Cabanes et al., 2003; Legagneux et al., 2003; Flanner and Zender, 2006; Taillandier et al., 2007). We instead assess the influence of wind speed on the SSA decay rate and observe a more rapid SSA decay with increased wind speed, which can be explained by increased ventilation of saturated pore air acting as a catalyst for snow metamorphism (Cabanes et al., 2003; Flanner and Zender, 2006; Neumann and Waddington, 2004). While wind erosion cannot be definitively ruled out even after applying the criteria described in Sect. 2.5.1, high wind speeds are documented to increase SSA via the re-deposition of fragmentation and sublimated snow crystals after suspension (Domine et al., 2009), meaning that SSA changes of this nature would not be captured in our analysis.

Unsurprisingly, given the parameters are fit to the data, the model defined in this study predicts observed SSA decay for the low- and moderate-wind events with the lowest RMSE. T07 (Eq. 4) underestimated the observed SSA decay rate in all the low- and moderate-wind events, except for E18. The largest error is associated with E1, which also had the highest mean wind-speed (6.9 m s⁻¹) of all analysed events. The tendency for T07 to underestimate the observed SSA decay can be explained by the additional influence of wind speed which accelerates SSA decay (Cabanes et al., 2002) but is not considered in either T07 or FZ06. By contrast, FZ06 consistently overestimates the observed SSA decay rate, most pronounced in E10 and E18. The original parameter values τ and n of FZ06 were tuned to data from alpine regions, potentially explaining the poor fit.

The simple empirical model presented here is limited to conditions at EastGRIP within a narrow temperature range (−18 to −7 ^∘C) and therefore might be unsuitable for sites with different conditions. However, large errors when using the models from the literature indicate that the low-latitude tuning and exclusion of wind effects are not optimal for predicting surface snow SSA decay at EastGRIP. Our study is limited by a relatively small number of unperturbed SSA decay events, hampering statistical evaluation over longer timescales.

4.2 Isotopic change during SSA decay events

4.2.1 Low-wind events

In the absence of snowfall or other surface perturbations, multi-day periods of snow metamorphism correspond to change in snow isotopic composition. We observe decreasing d-excess with SSA through time suggesting that the mechanisms driving snow metamorphism also influence the isotopic composition. However, the mechanisms linking these changes are unclear and are not always consistent in space and time. The following section explores the possible mechanisms driving isotopic change during SSA decay events, by assessing the LE and TG conditions during events with minimal surface perturbations.

Two key mechanisms are expected to drive the rapid SSA decay and concurrent change in snow isotopic composition: (1) snow grain growth via diffusion of interstitial water vapour due to near-surface temperature gradients (Colbeck, 1983; Ebner et al., 2017; Touzeau et al., 2018), observed to cause a decrease in d-excess and slight increase in δ¹⁸O in the defined snow layer (Colbeck, 1983; Ebner et al., 2017; Touzeau et al., 2016); (2) grain rounding via sublimation from convex regions of snow grains (Neumann and Waddington, 2004), causing an increase in δ¹⁸O and a significant decrease in d-excess of the remaining snow (Ritter et al., 2016; Madsen et al., 2019; Casado et al., 2021; Hughes et al., 2021; Wahl et al., 2021). Sublimation is enhanced by ventilation of the saturated pore air, known as “wind-pumping” (Neumann and Waddington, 2004). We note that isothermal metamorphism driven by Ostwald ripening also causes a decrease in SSA (Ebner et al., 2015), but is associated with minimal change in bulk isotopic composition which was not observed in our analysis. However, conditions that favour Ostwald ripening were not observed in our analysis.

Increases in δ¹⁸O and decreases in d-excess during E10 can be attributed to a combination of (1) and (2) based on the observation of net-sublimation and high amplitude diurnal TG variability over the course of the event. Net-deposition was measured during the period between 9 June at 15:30 UTC and 10 June 10:30 UTC 2018, corresponding to an overall decrease in δ¹⁸O, agreeing with previous studies (Stenni et al., 2016; Casado et al., 2021; Feher et al., 2021), and a minimal decrease in d-excess, which is not necessarily expected during deposition. However, disequilibrium between water vapour isotopic composition and snow isotopic composition may explain the deviation from expectation (Wahl et al., 2022). Although ground fog was documented on the day preceding E11, no significant deposition is observed in the LE data in the day preceding E11, indicating the absence of lasting surface hoar formation. The 30 % decrease in d-excess concurrent with no change in δ¹⁸O suggests strong kinetic fractionation during E11. Continuous variations in δ¹⁸O and d-excess throughout June 2018 (Fig. 7) show no clear relationship to total LE or temperature gradients. Field experiments looking at sub-diurnal variability show a stronger dependence of snow isotopic composition on LE (Hughes et al., 2021), potentially explaining the lack of a strong diurnal relationship.

Conclusively identifying these mechanisms requires measurements of water vapour isotopes to model the fractionation effects. In the absence of these data, we infer potential explanations for isotopic change during the low-wind events. Our analysis suggests that SSA of the surface snow is strongly influenced by surface-subsurface TG and wind speed, while the changes in isotopic composition are likely to be influenced by other factors, such as the magnitude of vapour–snow isotopic disequilibrium during sublimation (Wahl et al., 2022). Decoupling the influence of sublimation and interstitial diffusion within the snow requires additional measurements of isotopic composition of atmospheric water vapour to model associated fractionation effects (Wahl et al., 2022). Our results show that while snow isotopic composition does indeed change during SSA decay events, predicting the magnitude, and even the sign, of the isotopic change associated with snow metamorphism is not possible when information about the interstitial vapour isotopic composition is missing.

4.2.2 Inter-annual variability

The recurring difference in snow characteristics and temperature conditions during 2019 compared to 2017 and 2018 could be explained by the shifting phase of the North Atlantic Oscillation (NAO), which is in positive phase during 2017 and 2018 (below-average temperatures) and in negative phase during 2019. A similar pattern is observed in the snow isotopes when considering the period between 27 May and 1 August where the mean δ¹⁸O value during 2019 was −27.3 ‰, which is 4.3 ‰ higher than 2018 (−31.6 ‰) and 3.6 ‰ higher than 2017 (−30.9 ‰). While the difference in mean δ¹⁸O can plausibly be explained by a 3.2 ^∘C mean air temperature difference, the isotope-SSA covariance is not so straightforward.

The positive mode of PC1_SSA (Fig. 3) can be interpreted as increases in SSA from depositional events, such as precipitation, and wind-fragmented snowdrift (Domine et al., 2009), while the negative mode is associated with snow metamorphism or wind erosion (Cabanes et al., 2002, 2003; Legagneux et al., 2003, 2004; Taillandier et al., 2007; Flanner and Zender, 2006). Based on this interpretation, a covariance between PC1_SSA and PC1_dxs or PC1${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ indicates that the mechanisms controlling SSA variability also influence the isotopic composition. We observe a decoupling of the temporal variance in d-excess from that of δ¹⁸O (Fig. 3) in 2019. This difference may be related to the recorded transition from winter to summer in 2017 and 2018, but not in 2019.

Ongoing work to disentangle the processes driving change in isotopic composition – sublimation from surface or interstitial vapour diffusion between layers in the pore space – is vital for precise climate reconstruction in ice cores (Touzeau et al., 2018; Casado et al., 2021; Hughes et al., 2021; Wahl et al., 2021). Future studies would benefit from obtaining direct measurements of the isotopic composition of precipitation and surface hoar, to determine the fraction of such deposits in the SSA samples. Furthermore, a quantitative representation of vapour fluxes in the surface snow rather than the temperature-gradient-based approximation used in this study would provide a basis from which to quantify the relative influence of fractionation during sublimation and interstitial diffusion.

4.3 Implications and perspectives

Documented changes in snow isotopic composition during surface snow metamorphism have potential implications for interpretation of stable water isotope records from ice cores, given that the current interpretation assumes the precipitation signal is preserved (Dansgaard, 1964). Variations of d-excess in ice core records are attributed to changes in source region, source region conditions and kinetic fractionation occurring during snow crystal formation in supersaturated clouds (Stenni et al., 2010; Jouzel and Merlivat, 1984). However, we show that there are substantial changes in d-excess during precipitation-free periods supporting recent work showing that post-depositional factors add complexity to the traditional interpretation of ice core water isotope records. Greenland ice core records show abrupt decreases in d-excess corresponding to warming transitions, which is attributed to change in moisture source regions (Steffensen et al., 2008). Results presented in our study document decreases in snow d-excess during surface snow metamorphism associated with temperature gradients and sublimation, potentially contributing to the low d-excess values during warm periods.

The findings of this exploratory study reiterate the importance of quantifying the isotopic fractionation effects associated with processes driving snow metamorphism during precipitation-free periods. Moreover, the inter-annual variability observed at EastGRIP between 2018 and 2019 suggests that precipitation intermittency, temperature (gradients) and wind regimes play a role in isotopic change, which is not readily identified in the surface snow SSA data.

5 Conclusions

This study addresses the rapid SSA decay driven by surface snow metamorphism. In particular, the study aims to explore how rapid SSA decay relates to changes in isotopic composition of the surface snow in the dry accumulation zone of the Greenland Ice Sheet. Ten individual snow samples were collected daily over a 90 m transect at EastGRIP in the period between May and August of 2017, 2018 and 2019. SSA and isotopic composition were measured for each sample. Periods of snow metamorphism after deposition events were defined using SSA measurements to extract periods of rapid decreases in SSA. An exponential SSA decay model ($\mathrm{SSA}\left(t\right)=({\mathrm{SSA}}_{\mathrm{0}}-\mathrm{22})\cdot e^{-\mathit{\alpha }\cdot t}+\mathrm{22}$) was constructed to describe empirically surface snow metamorphism in summer conditions for polar snow, with surface temperatures above -20 ^∘C and with minimal surface perturbation. Two event categories were defined based on wind speed, with an upper threshold of 7 m s⁻¹ to minimise chance of snowdrift. Higher wind speeds increase the SSA decay rate due to enhanced snow metamorphism with increased ventilation of the pore space.

Changes in isotopic composition corresponding to post-depositional processes driving SSA decay are observed in all low- and moderate-wind events. A decrease in d-excess from Day-0 to Day-2 is observed in four out of the six events with no precipitation. Further analysis of low-wind SSA decay events indicates that the combined effects of vapour diffusion and diurnal LE variability causes isotopic fractionation of the surface snow in the absence of precipitation. The differing fractionation effects are expected to be the result of vapour–snow isotopic disequilibrium. A strong correlation observed between SSA and d-excess found in 2019 was not present for 2018. We suggest that this is due to strong sublimation corresponding to high temperatures during 2019.

In summary, our results support documentation of fractionation during sublimation and deposition between the snow surface and atmosphere, indicating that the precipitation isotopic composition signal is not always preserved during the processes driving surface snow metamorphism. Observations of post-depositional decrease in d-excess during rapid SSA decay hints at local processes influencing the d-excess signal, and therefore an interpretation as source region signal alone is insufficient.

Appendix A

Table A1 Description of conditions for SSA decays captured by decreased threshold. Duration and conditions for all 21 events defined by the threshold. Observations were made of snowfall, snowdrift and ground fog. Presented here are the dates, event number, mean surface temperature (ST), maximum wind speed (WS), accumulation (A), and surface conditions including those recorded during the day (24 h) preceding each event (Day-1), cloud cover, observations of surface perturbations during the events based on field observations (Day-0–n), and finally the wind speed category for SSA decay event analysis (Event category). In the “Surface conditions” columns, S, SD and F stand for snowfall, snowdrift and ground fog respectively. C in “Cloud cover” indicates clear sky and O indicates overcast.

Figure A1 Histograms showing (a) the daily maximum values and (b) the 10 min mean values for all sampling days of 2017, 2018 and 2019. The black line indicates the mean.

Figure A2 Results from the EOF analysis showing the relationship between SSA, d-excess and δ¹⁸O in the spatial dimension. Linear regressions are included in plots with significant correlations (p<0.05) between variables.

Figure A3 Results from the EOF analysis showing the relationship between SSA, d-excess and δ¹⁸O in the temporal dimension. Linear regressions are included in plots with significant correlations (p<0.05) between variables.

Figure A4 The change in SSA after 1 d (Day-1 to Day-0) plotted against the absolute value from Day-0 for all 21 SSA decay events (E1–E21 in Table A1) captured by the decrease threshold described in Sect. 2.5.1. The markers are coloured by mean air temperature between samplings (i.e. the mean air temperature from the time between sampling on Day-0 and Day-1).

Figure A5 SSA decay model evaluation and a comparison with observations. Included models are the decay model from this study and the existing decay models from Flanner and Zender (2006), FZ06, and Taillandier et al. (2007), T07. The 10 min averaged wind speed is shown on the secondary y axis, with the 6 m² kg⁻¹ thresholds indicated. Low-wind events E10 and E11 are shown in the upper panel, and moderate-wind events are shown in the lower panel (E1, E13, E18 and E19).

Figure A6 Daily mean 2 m air temperature (red), surface temperature (yellow) and 10 cm subsurface temperature (grey). The data are presented for the 2017, 2018 and 2019 measurement campaigns.