Arctic Snow Isotope Hydrology: A Comparative Snow-Water Vapor Study

: The Arctic’s winter water cycle is rapidly changing, with implications for snow moisture sources and transport processes. Stable isotope values ( δ 18 O, δ 2 H, d-excess ) of the Arctic snowpack have potential to provide proxy records of these processes, yet it is unclear how well the isotope values of individual snowfall events are preserved within snow proﬁles. Here, we present water isotope data from multiple taiga and tundra snow proﬁles sampled in Arctic Alaska and Finland, respectively, during winter 2018–2019. We compare the snowpack isotope stratigraphy with meteoric water isotopes (vapor and precipitation) during snowfall days, and combine our measurements with satellite observations and reanalysis data. Our analyses indicate that synoptic-scale atmospheric circulation and regional sea ice coverage are key drivers of the source, amount, and isotopic composition of Arctic snowpacks. We ﬁnd that the western Arctic tundra snowpack proﬁles in Alaska preserved the isotope values for the most recent storm; however, post depositional processes modi-ﬁed the remaining isotope proﬁles. The overall seasonal evolution in the vapor isotope values were better preserved in taiga snow isotope proﬁles in the eastern Arctic, where there is signiﬁcantly less wind-driven redistribution than in the open Alaskan tundra. We demonstrate the potential of the seasonal snowpack to provide a useful proxy for Arctic winter-time moisture sources and propose future analyses.


Introduction
Since the mid-1990s, Arctic surface air temperatures (SATs) have warmed at twice the global mean rate [1]. This amplification of the Arctic is driving pronounced changes in the hydrological cycle, including shifting ocean and atmospheric circulation, altered precipitation and humidity patterns, large-scale permafrost degradation, and precipitous ice-mass loss from glacierized regions throughout the Arctic [2][3][4][5]. In particular, the rapid loss of sea ice [6] has been linked to shifting atmospheric moisture source and transport regimes across the Arctic, with implications for the meridional transport of moisture between mid-and high-latitudes [7][8][9]. These impacts are evident, for example, during atmospheric river events, where plumes of warm water vapor are rapidly transported into the Arctic [10,11], and the advection of moisture from the Arctic to lower latitudes, such as during cold air outbreaks and transient cyclones [12,13]. Consequently, changes in the source, timing, amount, and distribution of snow are evident in western Arctic and are

•
Do the isotope profiles in late winter seasonal snowpack reflect individual precipitation events and their different winter atmospheric moisture sources? • Is post-deposition isotope fractionation over winter distinguishable in snowpack isotope profiles?   (Figure 1a). It is a long-term monitoring site for permafrost hydrology [72][73][74] and tundra snowpack processes [75][76][77], with over 30 years of snow depth and snow water equivalent measurements [14].
The long-term mean annual SAT is −7.7 • C and annual total precipitation is 334 mm   [78]. The snow cover season typically lasts from October to May [14], and the long-term mean end-of-winter snow-water equivalent and snow depth are 124 mm and 50 cm, respectively. Climate data for this study were obtained from the snow telemetry and snow course data network meteorological station (SNOTEL Station ID 968, Figures 1 and 2).
The Imnavait site is in the continuous permafrost region. Soils are clay-rich glacial till with a well-developed organic layer. Treeless tundra vegetation at the site is dominated by tussock sedge and dwarf shrub tundra [79]. The variable microtopography in tussock tundra results in high small-scale variability in snow depths [14]. The regional landscape topography comprises rolling hills and river valleys, where the open, treeless terrain with strong winds enhances the redistribution and sublimation of deposited snow [75]. These processes typically result in areas of relatively deep snow on the slopes west of the creek, and relatively shallow snow on the slopes east of the creek that are highly wind scoured ( Figure 1) [76]. . Key meteorological variables for daily (a) mean air temperature, (b) mean soil temperature, (c) maximum snow depth, (d) cumulative precipitation, (e) mean wind speed, (f) maximum short wave radiation (g) mean relative humidity for winter 2018-2019 at the study sites, with Imnavait in red and Pallas in blue.

Pallas, Finnish Lapland, Eastern Arctic
The Pallas site is located in Finnish Lapland at an elevation range of 270-360 m a.s.l (68.0 • N, 24.2 • E) (Figure 1b). The infrastructure has been utilized as a meteorological monitoring site for over 80 years, with intensive development in atmospheric and ecological research since the 1990s [80], and catchment ecohydrology since 2014.
The annual mean SAT and annual precipitation are 0.4 • C (2002-2019) and 639 mm (2008-2019), respectively. On average, the snow season starts in mid-October and ends in May, with automated snow gauging measurements recording a mean maximum snow depth of 105 cm during the 2007 to 2019 period. Regular snow surveys for snow water equivalent measurements started in 2018 along the same transect as snowpits in this study ( Figure 1). Maximum annual values in 2018, 2019, and 2020 were 279, 153, and 325 mm, respectively. Climate data for this study were recorded at the Finnish Meteorological Institute (FMI) Kenttärova weather station (Station ID 101987, Figures 1 and 2).
Despite its high latitude, the site is not underlain by permafrost and is frozen only seasonally with complete thaw soon after snowmelt. Soils on the hillslopes consist of glacial till, overlain by organic peat soils up to 3 m deep at the valley bottom peatland. Conifer stands of Norway spruce and Scots pine cover most of the site, except for treeless peatlands dominated by mosses and low sedges.
Snow research at the site does not have the long-term historical data as in Imnavait, but some differences in snow conditions between the sites are evident. First, the tree canopy at Pallas largely inhibits wind transport within the forest, though it does occur on the open mire. The tree canopy also shelters the snow from solar incidence and intercepts falling snow. Second, both snow depth and water equivalent are approximately double at Pallas compared to Imnavait.

Snowfall
Imnavait precipitation samples were captured using an open container (25 cm 3 ) and collected the morning following a snowfall event. Precipitation sampling was conducted 13 km west of the snowpit site at the Toolik field station (68.6611 • N, −149.3705 • E), and at an altitude of 843 m a.s.l. that is approximately 40 m lower than the lowest elevation snowpit (Figure 1b). Samples were collected from 12 snowfall events out of the 23 days with precipitation recorded at Imnavait meteorological station from the start of snow accumulation to snowpit sampling. The largest data gap occurred between late December 2018 and mid-January 2019.
Pallas precipitation samples were collected using the same equipment setup as for Imnavait, but more opportunistically due to limited field personnel. Sampling was conducted at the field station (68.0257 • N, 24,1602 • E) 4 km northwest of the snowpit sampling site (Figure 1a), at approximately the same altitude as snowpit 9. We were able to collect 14 samples out of 104 days with recorded snowfall at the Kenttärova meteorological station (Figure 1a) from the start of snow accumulation to snowpit sampling.

Snowpits
Snowpits were sampled in the late 2018-2019 snow season when the majority of the anticipated snowfall had occurred. Sampling was conducted [14][15] March 2019 at Imnavait, and 3-4 April 2019 at Pallas. Different sampling strategies were employed at the two sites due to contrasting snowpack conditions. At Imnavait the sampling protocol was based on snowpit stratigraphy, i.e., samples taken from each visually distinct layer. The Pallas sampling used fixed-depth increments of the snow profile, because a taiga snowpack contains mainly depth hoar (i.e., typically 4-10 mm large, striated, cup-shaped crystals formed within the snowpack by grain-to-grain vapor diffusion; see [81]), which limits the ability to identify individual snow layers, i.e., snowfall events. Therefore, incremental sampling at Pallas facilitated the comparison of profiles between pits.
At Imnavait, ten snowpits were analyzed and sampled. The sampling was conducted along a long-term snow survey monitoring transect [14]. The transect ran perpendicular to valley orientation, starting from the hill crest on the west (Pit 1), crossing the valley bottom, and climbing a west-facing slope to finish atop a water divide on the eastern side of the catchment (Pit 10) (Figure 1). The sampling targeted individual stratigraphic layers in the snow profile that were identified by hand hardness test (subjective, horizontal resistance test of each snowpack layer; see [82]), texture, and visual appearance (e.g., the wind-packed layers appear more bright white and opaque than the coarser glass-like depth hoar layers that are typically found in the bottom). To characterize each layer, we measured the density using density cutters of a 1000 cm 3 or 250 cm 3 volume (Pedersen et al. [83]) and identified the grain form (i.e., the shape of the snow grain/crystal) per layer according to the international classification of seasonal snow on the ground [81].
After the profile was characterized, depth-integrated samples were taken from each layer using a plastic coring tube. The sampling progressed from top to bottom in the snowpack. At Pallas, nine snowpits were sampled and analyzed. The snowpits were dug along the snow-survey monitoring transect, starting from the hilltop at the Kenttärova station (Pit 1) and ending at the valley bottom (Pit 9) (Figure 1b). The pits were dug at Atmosphere 2021, 12, 150 6 of 32 roughly 200 m apart spanning the elevation gradient in the study catchment. The sampling protocol consisted of taking incremental 5 cm snow layer samples using a cubic-shaped 5 cm-deep metal snow sampler. Gaps in the Pallas profiles ( Figure 3) are caused by a batch of samples lost in sample processing. Example images of snow pit stratigraphy and sampling strategy for both sites are presented in the Appendix A, Figure A7. boxplots for d-excess values for snow pit and precipitation samples for Imnavait and Pallas. Note how water line plots depict the generally colder conditions in North Alaska and the lower isotope composition compared to those at Pallas. Slopes and standard errors for regression equations are Imnavait pits: 8.14 ± 0.14, Imnavait precip: 8.13 ± 0.28, Pallas pits: 7.46 ± 0.07, Pallas precip: 7.95 ± 0.28.

Snow Stable Isotope Analysis
All snowfall and snowpit samples were stored initially in the field in zip-lock airtight plastic bags. The bags were transported frozen to the laboratory facility where they were left at room temperature until all snow was melted, typically overnight. The liquid water samples with organic particles were filtered, then transferred to plastic 5 or 15 mL vials, and stored refrigerated at +4 • C until analysis.
Stable isotopes of water were measured at the University of Oulu using laser Cavity Ring-Down Spectroscopy (CRDS) analyzers: Picarro L2130-i and Picarro L2140-i. The Picarro system is composed of an autosampler, a high-precision vaporizer, and a laser CRDS analyzer. Picarro L2130-i and Picarro L2140-i analyzers have guaranteed precision of 0.1 and 0.025‰ for δ 18

Atmospheric Water Vapor Analysis
Atmospheric water vapor isotope measurements were conducted at both sites simultaneously during winter 2018-2019. At Imnavait, data were retrieved from the on-site National Ecological Observatory Network (NEON) eddy tower at Toolik station (68.66109 • N, 149.37047 • W, 832 m a.s.l.). The tower sampled ambient air at four heights (10, 20, 30, and 40 m), and here we utilize data from the top collection height for the best estimate of the isotopic composition of the free troposphere prior to (or after) its interaction with the snowpack. At Pallas, measurements were conducted at the FMI Sammaltunturi station (565 m a.s.l.,~280 m above lowest snowpit), approximately 4 km southwest from the snowpit transect (67.9733 • N, 24.1157 • E). Detailed descriptions of the instrument setup, sampling, and data calibration protocols at each site are presented in Appendix A.
Stable isotope ratios (δ 18 O and δ 2 H) were measured at both sites using a Picarro L2130-i Isotope and Gas Concentration Analyzer (herein analyzer) with a Standards Delivery Module (SDM). In summary, ambient air was pumped continuously to the analyzers where isotopic ratios were measured approximately every second (1 Hz). Isotope data were then calibrated using established multi-step protocols [85,86]. First, the humidity-isotope response of each instrument was established in the field to correct for potential analytical bias [86][87][88]. Using an integrated dry air system, the vapor streams of standard waters with a known isotopic composition were measured across a range of controlled humidity injections (Toolik:~600 to 24,500 ppmv, Pallas:~400 to 10,000 ppmv) that encompass the range of ambient humidity values during the campaign. Results of the humidity experiments are shown in Figures A1 and A2. A non-linear regression was used to determine the isotope correction as a function of cavity humidity, and applied to the ambient data to remove the humidity bias. Second, ambient measurements were corrected and standardized to the Vienna Standard Mean Ocean Water (VSMOW2) scale. Approximately every 12 h at Pallas and 24 h at Toolik, water standards with a known isotopic composition were injected by the SDM and analyzed to calibrate and monitor for potential instrument drift. At both sites, the analyzer remained stable during the measurement period, with relatively small drifts between successful standard runs that were substantially smaller than the natural variations examined in this study. All ambient and standard measurements were corrected for the humidity concentration dependence, and at both sites, the standard error was estimated at <0.3‰ for δ 18 O and <1.0‰ for δ 2 H and d-excess.
Extreme outliers were excluded from our analyses, including d-excess values less than −15‰ and greater than 70‰. Data contaminated by the dry air system used for standard runs (i.e., isotopic values clearly outside the natural range) were also excluded. At Toolik, this included two longer periods between 21 December 2018-4 January 2019 and 29 January-1 February 2019. The 1 Hz isotopic measurements were initially aggregated to 1 min values at Toolik and 5 min values at Pallas. Daily mean isotope values were also calculated at each site for comparing with precipitation isotope data. For full details see Appendix A.

Synoptic Climatology Analyses
We used the NCEP/NCAR reanalysis dataset [89] to evaluate the broad synoptic-scale atmospheric conditions during the winter 2018-2019 measurement campaign. NCEP/NCAR reanalysis comprises data with 2.5 • horizontal grid spacing on 45 vertical levels and are available at 4 h intervals. Here, we used the daily mean 850 mb geopotential heights (i.e., average of the 0z-18z 4 h data) to construct composite maps during all snowfall events sampled at Imnavait and Pallas between 30 October 2018 and 13 March 2019 (Table A1). Daily mean 850 mb vector winds were also used to evaluate the atmospheric transport patterns associated with each individual snowfall event (Appendix A.1, Figures A3 and A4). Additionally, daily mean Arctic sea-ice concentration grids (25 km 2 resolution), derived from the Special Sensor Microwave Imager/Sounder (SSMIS), were obtained from the National Snow and Ice Data Center (NSIDC) [90] and used to construct composite sea ice concentration map during the snowfall sampling days.

Data Analysis and Visualization
Snowpit isotope data are presented as vertical isotope stratigraphy (Figures 4 and 5), where zero represents the base of the snowpack, i.e., the snow-ground interface. To facilitate comparison of vertical isotope variability in snowpits within sites, the snow layer depths were normalized to the total snow depth in each pit. This results in the relative depth of the snow layer within the snowpack (0-100%).  To facilitate comparison between precipitation and snowpack, the plotting position for precipitation on the y-axis is calculated as the relative contribution of amount of each snowfall event to the total snow amount. The depth of each precipitation event is recorded as mm snow water equivalent, which is used as the unit of summation. For example, the first two sampled storm events were each 5 mm in magnitude and captured by one snowfall sample (Figure 8a). Each event was comprised of 6.5% of the total 77 mm of accumulated precipitation during the observation period.
The whole profile was normalized to cumulative snowfall, which gives the relative depth of each snowfall event (0-100%) with respect to accumulated snowfall until the sampling. Start and end plotting position of each precipitation sample on the y-axis are given by Equations (1) and (2), respectively. where: pp start is the y-axis start plotting position in the beginning of the sample i pp end is the y-axis end plotting position in the beginning of the sample i p is the precipitation amount in sample i (mm) n is the number of precipitation events For vapor data, the plotting position was calculated using the same procedure. The mean vapor value of each day was assigned to the precipitation amount during days with snowfall.

Climate and Snow Conditions at Imnavait and Pallas
At Imnavait, mean air temperature was −12.1 • C during the study period ( Soil temperatures at 5 cm depth showed marked differences between sites (Figure 2b). At Imnavait, soil temperature decreased linearly from −1.0 • C in October to −6.5 • C in mid-January, after which it varied between −6.5 • C and −4.0 • C, and with a seasonal mean value of −3.9 • C. At Pallas, the mean soil temperature was warmer (−0.5 • C) and remained consistent throughout winter. These differences in soil temperature reflect colder air temperature at Imnavait (Figure 2a), coupled with deeper snow cover at Pallas from December onwards (Figure 2c) that effectively insulated and decoupled the soil from cold air temperatures.
We observed a typical snow depth pattern at Imnavait, which is persistently similar between years as reported in previous studies [14,76]. The deeper snow occurred within depositional areas in the bottom of the catchment, whereas shallower snowpacks covered the wind-scoured hillslopes. The snow gauge recorded a depth of 41 cm during the snowpit sampling in April 2019 (Figure 2c), whereas the average across all 10 pits was 34.6 cm (min 19 cm, max 54 cm). The values are below the long-term snow survey value of 50 cm (1985-2017), which is, however, done later in the snow season allowing more snow accumulation [14]. At Pallas, the gauge-recorded snow depth during sampling was 93 cm, below the long-term mean annual maximum snow depth (105 cm, 2007 to 2019). Mean snow depth in pits was 85.1 cm (min 75 cm, max 93 cm). High-resolution snow depth mapping at the site by Meriö et al. [personal communication] showed that snow depth at the automated gauge is consistently higher than that of the rest of the study area, as recorded in the snowpits.
The data indicate that snowfall and snow accumulation started earlier at Imnavait and were consistently increasing during October and November ( Figure 2c). From February 2019 onward, the snow depth at Imnavait remained at approximately 40 cm, whereas the snow depth at Pallas increased from 50 to 100 cm between early February and the end of March 2019. Accordingly, a similar evolution was observed for snowfall. Total water equivalent of snowfall received during the observation period was 23.5 mm at Pallas and 77 mm at Imnavait.

Isotopic Composition of Snowpack and Precipitation
The snow and precipitation δ 18 O values at Imnavait were lower compared to Pallas ( Figure 3a and Table A1). The differences for both δ 18 O and d-excess were statistically significant between sites (Imnavait and Pallas), but not between snowpit and precipitation samples at either Imnavait or Pallas (Tukey's post hoc test, p = 0.05). The δ 18 O composition of snowpit samples at Imnavait were slightly higher (the mean −27.4‰, standard error of the mean ±0.6‰) compared to precipitation (−29.0, ± 1.3‰), but the range in snowpit and precipitation samples was similar (−35 to −20‰) (Figure 3a). At Pallas, the mean δ 18 O values for snowpit (−20.6, ± 0.2‰) and precipitation (−20.5 ± 1.3‰) were similar, but the range was larger in precipitation (−30 to −12.5‰) than in snowpack (−25 to −15‰). This is despite fewer precipitation samples (n = 14) compared to the number of snow samples (n = 111).
Regression lines (R 2 > 0.98) of Imnavait snowpit and precipitation and Pallas precipitation were close to the global meteoric water line (GMWL, δ 2 H = 8 × δ 18 O) (Figure 3a), with slopes ± 7.95 to 8.14), whereas the Pallas snowpit regression slope of 7.45 was slightly lower than the GMWL. Both sites exhibited a positive correlation between air temperature and δ 18 O values (Appendix A, Figure A6). At Pallas we observed a temperature-δ 18 O slope of 0.60‰/ • C for vapor and 0.68‰/ • C for snowfall. At Imnavait the temperature-δ 18 O slopes were lower, with 0.47‰/ • C for vapor and 0.09‰/ • C for snowfall.
The d-excess values at Imnavait ranged between 0 and 15‰ in both the snowpack and precipitation, with mean values of 7.3 and 8.0‰, respectively. However, the median d-excess was lower in the pit samples compared with precipitation with an offset of 1.7‰ ( Figure 3b). At Pallas, mean snowpack and precipitation d-excess values were 10.4 and 11.2‰, respectively, and collectively ranged between 5 and 17‰. Overall, the Pallas d-excess values median and range between precipitation and pits were similar, with an offset of 0.3‰. Figure 4 illustrates the high variability in δ 18 O between the snowpit profiles at both Imnavait and Pallas. Snow depth at Imnavait was lower and more variable among individual pits compared to at Pallas. The δ 2 H values in snowpit profiles are shown in Appendix A, Figure A5. Because δ 18 O and δ 2 H values are closely correlated ( Figure 3a) and visually difficult to find any differences in the profile plots, we present only the δ 18 O and d-excess values in the subsequent plots.

Isotope Stratigraphy in the Snowpits
To facilitate comparison within and among sites, in Figure 5 we also present combined snowpit profiles. At both sites we found an overall lowering δ 18 O trend with increasing depth ( Figure 5). At Pallas, the typical range in isotope variability between pits at all depths is~5‰, with higher variability observed closer to the top in the snow profiles (above 50 cm snow depth). At Imnavait, the first 10 cm at the base of the profiles has consistently higher δ 18 O values compared to the top of the snowpack. Variability in δ 18 O values increases above the snow depth of~10 cm, without any clear patterns or trend as a function of snow depth.
Isotope stratigraphy of snow pits at Imnavait is highlighted by plotting snowpack isotope samples in each snowpit relative to the maximum snow depth of the pit and visualizing snow density (Figure 6a), and grouping samples based on their snow grain type (Figure 7). The top of the snowpack (~75-100% of total depth) comprised low-density snow (~100 kg m −3 ) with low δ 18 O values, most less than −30‰ (Figure 6a). The middle of the snowpack (~25-75% of total depth) consisted of a high density (~300-400 kg m −3 ) layer where the δ 18 O values varied between −25 and −30‰. The base of the snowpack (~0-25% of total depth) consisted of a moderate dense (~200 kg m −3 ) layer, where the snow δ 18 O values were higher relative to layers above it, with most samples ranging between −20 to −25‰. One of the key factors that appears to explain the δ 18 O values of snow is the grain type ( Figure 7). For instance, at Imnavait rounded-grain snow crystals found in the upper layers of the snowpack had typically low δ 18 O values compared to those in lower layers of the snowpack ("wind packed" or "depth hoar"), differing by up to 20‰. Plots of the snow profile isotopes relative to total depth at Pallas (Figure 6b) did not reveal much more consistency in the layering compared to absolute snow depth (Figure 5b).
At Imnavait the early season precipitation (30% of total snowfall, occurring between October and December 2018) had δ 18 O values that ranged between −26 and −30‰ (Figure 9a). However, the corresponding snowpit samples at the base of the snowpack exhibit higher δ 18 O values (−20 and −25‰). Gaps in precipitation sampling in the middle of the season, corresponding to the wind-packed snow layers in the middle of the snowpack, did not allow reliable comparison between precipitation and snow layers. The two storms which were sampled (~40-50% of total snowfall) were characterized by lower     Due to logistical constraints, Pallas precipitation samples were only collected during mid-winter (events that contributed 40-65% of the total cumulative snowfall) (Figure 8b). The δ 18 O composition of these events was comparable to values observed in the middle of the Pallas snowpack, but with higher variability in the precipitation samples. The few precipitation samples over the season precluded rigorous comparison with the rest of snow isotope profiles. The Pallas vapor δ 18 O composition was consistently lower compared to the precipitation and snowpack, (respective means −28.0, −20.4, and −20.6‰, Figure 9b). Unlike Imnavait, the vapor and snow δ 18 O values were less offset, at some points overlapping with the snow profile. The snow profiles and vapor record showed similar isotope dynamics. Both transitioned from relatively high to low δ 18 O values from early to mid-season (0-60% of normalized snow profiles and vapor composition of snowfall days). In the top profile (60-100%), both the snow and vapor isotope composition remain variable but relatively constant vertically.
The method for comparing the snowpack, atmospheric water vapor, and precipitation d-excess values presented in Figure 10, is defined in Equation (1)  At Pallas, a few precipitation samples exhibited d-excess values similar to those recorded in the snowpack at corresponding depths ( Figure 10b). Moreover, we found that the snowpack d-excess profile from 0 to 75% of snow depth closely followed the vapor isotope composition during snowfall days. Even though the overall trends were similar, the snowpack samples were less variable than the corresponding storm day vapor composition. Vapor d-excess composition during snowfall events that comprised 75-100% of the total snowfall (corresponding to events after the beginning of March) was between 12 and 24‰, which were higher compared to the top 75-100% of the snowpack, with values between 5 and 13‰.

Synoptic-Scale Moisture Transport and Isotope Hydrology
Composite maps depicting the mean synoptic climatology during the sampled snowfall events are presented in Figure 11, and individual event-based maps are shown in Figures A3 and A4. In the western Arctic, daily composited 850 mb geopotential heights indicate a large, elongated region of low pressure centered over the Bering Sea region during the Imnavait snowfall events (Figure 11). The net result was strong westerly and southwesterly winds that typically tracked from the North Pacific and across the Gulf of Alaska or Bering Sea, from where air masses transited across continental Alaska to Imnavait ( Figure A3). Individual event-based wind vector analyses indicate that these circulation patterns were associated with snowfall characterized by a range of isotope values between −18.2 and −35.3‰ for δ 18 O, and 1.8 and 14.2‰ for d-excess (Table A1). The final snowfall event sampled at Imnavait on 13 March 2019 reflected the convergence of strong southerly airflow with northeasterly winds from the Beaufort Sea ( Figure A3). The event was associated with snowfall characterized by δ 18 O and d-excess values of −33.3 and 2.6‰, respectively. Overall, we found that mean sea ice concentration in the western Arctic was high during the 2018-2019 measurement campaign, with mean ice cover in the Beaufort and Chukchi Seas at 93 and 60%, respectively, during the snowfall sampling days at Imnavait (Figure 11).  [89], and the shaded region represents daily mean composited sea ice concentration during the corresponding snowfall sampling days [90]. The dashed blue circle depicts the Arctic Circle at~66 • N. Sampling dates used to construct the composite maps are given in Table A1.
In the eastern Arctic we observed a different synoptic configuration during the snowfall events sampled at Pallas. Daily composited 850 mb geopotential heights indicate a large area of low pressure centered over the Barents-Kara Sea region coupled with relatively high surface pressure in the North Atlantic ( Figure 11). This configuration resulted in prevailing west-southwesterly and southerly airflow that transported air masses to Pallas from the North Atlantic region (~80% of the snowfall events) ( Figure A4). Additionally,~20% of the winter 2018-2019 snowfall events we sampled derived from cyclonic north-northwesterly air masses that flowed across the Norwegian and Barents Seas ( Figure A4). Individual event-based wind vector analyses indicate that Pallas winter snowfall events deriving from Atlantic air masses were characterized by mean δ 18 O and d-excess values of −20.9 and 9.8‰, respectively, whereas the northerly air masses delivered snowfall with relatively higher mean δ 18 O and d-excess values of −17.6 and 20.4‰, respectively (Table A1). Overall, compared to the western Arctic we found that regional sea ice concentration in the eastern Arctic was relatively low during the 2018-2019 measurement campaign, with~22% mean ice coverage in the Barents Sea during the snowfall sampling days (Figure 11).

Isotope Fractionation in the Snowpack
We found that the atmospheric vapor δ 18 O at the sites was low relative to the corresponding daily precipitation and snowpack values. The average differences between vapor and precipitation were 11.4 and 7.6‰, and vapor and snowpack were 13.4 and 7.4‰ for Imnavait and Pallas, respectively (Figure 4), and consistent with those observed in other regions of the Arctic [12,91]. These differences reflect equilibrium fractionation during phase change, which results in isotopic enrichment in heavy isotopes in the solid phase relative to vapor. Modeled equilibrium vapor-ice fractionation at Imnavait and Pallas average air temperatures (−12.1 and −7.7 • C) using fractionation factors 1.0166 and 1.0159 [92], respectively, results in average vapor-ice fractionation of 16.1‰ at Imnavait and 15.4‰ at Pallas. The observed differences are lower than the modeled ones, for Pallas in particular, suggesting non-equilibrium conditions during condensation, post-depositional isotope modification, or both.
There are both similarities and differences in post-depositional isotope modification between the western and the eastern Arctic. At both sites there is a general tendency that the base of the snowpack has higher δ 18 O values compared to the top (Figures 5 and 6). This pattern is likely in part due to the seasonality effect [84]; temperatures of snowfall early in the winter (October and November) occur under warmer temperatures and from moisture sources with warmer water temperatures compared to mid-and late winter (December, January, February and March) when snowfall occurs at lower temperatures ( Figure 2) [5]. The seasonality effect is corroborated by the positive correlation between air temperature and isotope values in both precipitation and vapor samples at both sites ( Figure A6). The seasonality effect is evident at the Pallas site in particular, when comparing the snow isotope profile with atmospheric vapor isotopes during days of snowfall (Figure 9b). Both snowpack and vapor isotope composition gradually shift from higher to lower δ 18 O values as the winter advances. This enforces the observation that snowpack may provide an isotope proxy of seasonality in moisture sources in seasonal snowpacks.
At Imnavait, samples of early season precipitation enabled direct comparison of precipitation with the base of snowpack (Figure 9a). The comparison showed that the snowpack base was characterized by higher δ 18 O values compared to the corresponding precipitation events. This suggests that the seasonal isotope signal was not preserved within the snowpack, and there was post-depositional enrichment of the heavy isotopes over winter. Heavy isotope enrichment of basal snow layers has previously been attributed to soil water diffusion into the snowpack over winter, and molecular diffusion in vapor transport through the snowpack. Friedman [58] and Sturm and Benson [56] showed how preventing soil water from diffusing into the snowpack resulted in lower isotope values at the base of the snowpack compared to a undisturbed reference snowpack. Sinclair and Marshall [47] found that the base of snowpacks on a glacial moraine had higher isotope values compared with snowpacks on the glaciers, and attributed this to soil water availability of the moraine snowpacks. At Imnavait, the shallow tundra snowpack had strong temperature gradients between soil and air ( Figure 2) inducing vapor movement through the snowpack [56], as evidenced by the depth hoar layers (Figures 6 and 7). It is therefore very likely that the higher isotope values at the base of the Imnavait snowpack are not only explained solely by seasonality, but also isotope modification and interaction with soil water over the winter. In the boreal snowpack at Pallas, such comparison was not available due to paucity of data (Figure 9b), but the two sampled early season snowfall samples did not indicate a clear shift like that in Imnavait tundra snow.
At Imnavait, the snowpack is characterized by higher δ 18 O values relative to snowfall at the site, and more so than at Pallas (Figure 3). In addition to the discussed potential soil water-snowpack interaction in tundra snow, the differences between snowfall and snowpack can be indicative of sublimation and evaporation isotope fractionation. In general, conditions become more favorable for sublimation in high wind speeds, low air relative humidity, and higher solar insolation [93]. When comparing the two sites (Figure 2), Imnavait experienced higher wind speeds (mean 3.3 m/s, standard error of the mean ± 0.2 m/s) compared to Pallas (2.7 ± 0.1 m/s), and a more variable wind regime with more days of both low and high winds (Figure 2e). Wind transport is known to intensively redistribute snow at the site [75,76]. Relative humidity was on average higher at Pallas (86.8 ± 0.7%) in comparison to Imnavait (72.4 ± 0.9%) (Figure 2g). The Pallas site is dominated by conifer canopies, which reduces wind transport, provides shading from solar insolation, and increases relative humidity at the snow-atmosphere interface.
The most recent precipitation event prior to snowpit sampling at Imnavait provided key insights to potential snow isotope fractionation due to snow sublimation during wind transport. Surface snow in all pits had δ 18 O values that were consistently higher compared to the precipitation event less than two days prior to the snowpit sampling (Figure 9a). One explanation is that sublimation, enhanced by wind transport, immediately fractionated the fresh snow during the two days between the storm and the pit sampling. At Pallas, isotope values of the most recent snow event of~6 mm (2-3 April 2019, i.e., one day before snowpit sampling) was not clearly identifiable at the top of snow pits. This may be due to different sampling strategies between the sites. At Pallas the constant sampling interval of 5 cm was likely to capture snow also from earlier storms. At Imnavait the sampling was based on stratigraphy, where the most recent storm created a low-density layer that was sampled as whole (see Figure A7).
Snowpack samples at both sites lie close to the GMWL, with a δ 2 H-δ 18 O slope near 8 ( Figure 3). Experimental and field studies have found δ 2 H-δ 18 O slopes of~5-6 in snowpacks undergoing kinetic sublimation fractionation , whereas the slope in our snow samples was close to 8, resembling the GMWL. Only the Pallas snowpit samples had a slope lower than the GMWL, which could indicate kinetic fractionation. The lack of kinetic fractionation signal was particularly surprising at Imnavait tundra snowpack, where modeling by Liston and Sturm (1998) estimated 9-22% of the Imnavait snowfall precipitation input to sublime, with higher sublimation rates on wind-scoured hillslopes. The lack of kinetic fractionation signal in Imnavait snow samples could be explained by whole-grain ice-vapor transition in blowing snow [58], in which the whole snow grain sublimates and does not leave a residual higher δ 18 O in the remaining snowpack. This would result in isotope values similar to equilibrium fractionation, i.e., proportional change in both δ 2 H and δ 18 O isotope, as suggested by the Imnavait surface snow layer and snowfall comparison.

Contrasting d-excess Values in Snowpack and Vapor between Sites
Comparing d-excess values between vapor, precipitation, and snowpack provided an additional tracer to understand the isotope fractionation processes in the seasonal snowpack. Condensation of vapor to rain or snowfall is generally considered an equilibrium fractionation process [84]. Even though equilibrium fractionation results in enrichment in heavy isotopes from vapor to liquid/solid phase, the d-excess parameter should remain largely unchanged from vapor to precipitation and snowpack. This model generally seemed to be the case at Pallas (Figure 10b), but surprisingly not at Imnavait (Figure 10a).
At Imnavait, the vapor d-excess values were substantially higher compared to snowpack and precipitation d-excess values. Precipitation d-excess values had a similar range as the snowpack samples (Figure 10a), but individual snow layers d-excess values did not clearly resemble those in the precipitation. Only the top of the snowpack clearly reflected d-excess values from the recent snowfall event (as for δ 18 O in Figure 9a). One possible reason for the difference between vapor to precipitation and snow is that the strong winds at the Imnavait site caused some fractionation in the precipitated snow before the actual sampling.
At Pallas, vapor d-excess was more variable compared to snowpits or precipitation profiles but, unlike at Imnavait, the trend and central values in vapor d-excess resembled the snowpack d-excess remarkably well between 0 and 75% of the total snow depth (Figure 10b). Interestingly, the d-excess values in vapor during snowfall days were higher compared to the snowpack profile in the late winter (>75% of total snowfall/precipitation). Shortwave radiation at the site increases rapidly at the sites in late winter, and the Pallas sampling was conducted three weeks later in spring 2019 (Figure 3f). Increased radiation load can enhance snow sublimation with kinetic isotope fractionation, leading to lower d-excess values in the remaining snowpack [51,96], which may in part explain the difference in Pallas vapor and snowpack d-excess in the late winter.

Regional Drivers of Snow Amount and Snowpack Isotope Signatures
Collectively, our analyses indicate that synoptic-scale atmospheric circulation is a key driver of the source, amount, and isotopic composition of accumulated snowfall in Arctic snowpacks. Overall, Pallas received approximately twice as much snow compared to Imnavait in winter 2018-2019 ( Figure 2). We propose this phenomenon can be partly attributed to regional differences in sea ice cover and thus the availability of open (ice-free) sources of atmospheric moisture. For example, the Beaufort and Chukchi Seas that surround northern Alaska were ice-covered for the duration of our winter 2018-2019 measurement campaign, thus limiting the availability and transport of atmospheric vapor from the north [12]. Consequently, our atmospheric transport analyses indicate that all snowfall events collected at Imnavait were associated with prevailing southerly air masses from the North Pacific region (Figure 11). To the contrary, in the western Arctic, large areas of the Barents Sea to the north of Scandinavia remained ice-free throughout winter 2018-2019 and, together with the Norwegian Sea, provided a constant source of evaporative moisture that could be transported to Pallas and condense as snowfall [13,97,98]. Such differences are also evident when comparing the winter 2018-2019 relative humidity data between sites, with Imnavait characterized by considerably lower relative humidity, overall, compared to Pallas (Figure 2g). These major differences in the distribution of Arctic sea-ice, and thus the availability of evaporative sources of atmospheric moisture for precipitation, are becoming increasingly evident [7,44,91,99], and reflect complex patterns of atmospheric moisture transport into, within, and out of the Arctic [100].
We propose that the higher variability in δ 18 O values at Pallas (Figure 3a) reflects these differences in air mass transport trajectories, together with a wider variety of oceanic vapor source regions that feed continental precipitation. For example, at Pallas the winter snowfall events during 2018-2019 are sourced from air masses transported across the Barents, Norwegian, Atlantic, and Baltic Seas, and continental areas of Eurasia (Appendix A, Figure A4). Alternatively, in Alaska our analyses indicate that nearly all air masses derived from the North Pacific/Gulf of Alaska region. Thus, moisture transported by the air mass to Imnavait is characterized by less isotopic variability, but there is also significant isotopic fractionation along this storm track.
These differences in air mass transport regimes are reflected in the overall tendency for Imnavait vapor and snowfall to be characterized by moisture with δ 18 O values~5‰ lower than at Pallas. This includes (1) south to north fractionation effects with this snow delivery pathway over the Alaskan landmass (continental effect) [12,35,42,43,101], and (2) moisture transport with further altitudinal fractionation over the Alaska Range (altitudinal effect). To the contrary, the relatively close proximity of the Barents and Norwegian Seas to Pallas, coupled with the overall lower topography of Scandinavia compared to Alaska, means that moist air masses are subject to less rain-out and distillation effects during transit, and thus the Pallas vapor, snowfall, and snowpack are all characterized by relatively higher δ 18 O values compared to those at Imnavait (Figure 9).

Uncertainties and Wider Implications
The ranges of typical snow δ 18 O values between the different snowpits at any given depth were 2-5‰ for Pallas and 5-10‰ for Imnavait (Figures 5 and 6). Most existing studies have used a single or a few snow profiles to characterize the regions' snowpack isotopically [56,65,66,102]. However, our analyses indicate that this limited sampling may be misleading. By sampling multiple snowpits, as in our study, common patterns in the regional snow isotope profiles start to emerge, but there is still pronounced spatial variability among pits ( Figure 6). The 5 cm incremental sampling at Pallas allowed better vertical isotope resolution for comparison between pits, which was lost in layer-based sampling at Imnavait. High-frequency incremental sampling with detailed physical snow layer analysis would be the optimal strategy, however, it is work-intensive to sample multiple pits. For a reliable comparison of snowpit isotope stratigraphy with vapor or precipitation records in Arctic snowpacks, we recommend sampling multiple snowpits. Challenging field conditions and remote locations of our Arctic study sites produced unfortunate gaps in both precipitation time series and snowpit profiles. Complete records of vapor, precipitation, and snow profiles are needed to further validate some of our findings.
Our data showed significant isotope variability in snow profiles at short 1-2 km proximity. In contrast, Sinclair and Marshall [47] showed fair agreement in snow isotope stratigraphy in two snowpits 160 km apart in the Canadian Rocky Mountains. Their high-altitude maritime snowpacks had approximately five times the snow compared to Pallas, making the isotope stratigraphy deposited from each large synoptic snowfall event more distinguishable. Guided by their findings, depositional areas of snow might prove useful for sampling snowpack isotope profiles for proxy moisture sources. In tundra snowpacks experiencing significant wind transport, such areas could be natural (lee slopes) or manmade (snow fences).
Our vapor-precipitation-snowpack comparison is complicated by the implicit assumption in Equation (1) that each deposited snow layer preserved its depth over the winter. In reality, snowpack densification is likely to compact the lower snowpack layers relative to the top [56]. This could be addressed by simulating snow layer evolution with physically based snowpack modeling in which snow layer densification is explicitly represented [103]. Frequent storms depositing the snow layers at Pallas (Figure 2d, Figure 9a) and the marked spatial variability in total snow depth and layering at Imnavait (Figure 4a) would complicate any explicit snow layer simulations and layer identification, but this should be further explored.
Isotope diffusion models offer another potential simulation tool to make snowpit and atmospheric isotope data more comparable [48]. The paleoclimate community has highlighted the importance of post-depositional processes modifying the original snow isotope signature in ice and firn profiles [104][105][106]. These isotope modifications may cascade into uncertainties in ice core isotope profiles used for paleoclimate reconstruction. Isotope diffusion models have been successfully used to reproduce the observed diffusive smearing in the multi-season firn and ice-water isotope profile [47,68]. However, model equations and boundary conditions are designed for deep firn and ice profiles (meters to tens of meters). The applicability of firn isotope diffusion models to tundra or taiga snowpacks with high temperature gradients and additional water source from underlying soil needs to be tested.
Post-depositional modification of water isotopes in natural snowpacks has been simulated from the perspective of ice-water isotope exchange (hydrological scope, see, e.g., [53,65,107]) and ice-vapor isotope exchange (paleoclimate scope, see e.g., [47,48,68]. However, the existing models are not suited to capture key processes causing isotope evolution in Arctic seasonal snow layers: vapor diffusion from soil to snowpack, and isotope fractionation during sublimation. A way forward could be the introduction of isotope simulation routines into physically based snow models already characterizing the water fluxes (vapor flow, sublimation) and snow densification [103,107,108].
As the New Arctic is upon us, with shifting moisture source transport, the use of isotope stratigraphy in snowpacks as seasonal reservoirs of water isotopes can help to document how these shifts are manifested across the north. Robust moisture source trajectory modeling coupled with a pan-Arctic snowpack isotope dataset could help understand the differential influence of important factors, such as sea ice retreat on water movement across the Arctic. Snowpack isotope sampling in the late snow season could allow moisture proxies over larger spatial gradients than possible with fixed liquid and vapor water isotope sample locations.

Conclusions
Our analyses indicate that both the total amount of snow and isotope values in Arctic snowpacks are influenced by synoptic-scale atmospheric circulation. Our atmospheric transport analyses associated western Arctic tundra snow with prevailing southerly air masses from the Gulf of Alaska region, and eastern Arctic taiga snowpack with more variable transport regimes, and thus more variable sources of atmospheric moisture. The differences in source areas were seen in the relatively higher δ 18 O in the eastern Arctic snow.
Our intensive in situ atmospheric vapor isotope monitoring allowed comparison of snowpit isotope stratigraphy and meteoric isotope values. The analysis suggested the variability atmospheric moisture sources over winter were to some extent recorded in our eastern Arctic taiga snow profile. Our western Arctic tundra snowpack profile experienced more post-depositional isotope modification, complicating the use of tundra snow profiles as a proxy for seasonal evolution in moisture sources. We propose that isotope studies in tundra snow could be further advanced by using different sampling strategies and numerical modeling.
Whilst these findings also have critical implications for proxy-based studies that assume isotopic preservation within individual snow/firn layers, our analyses suggest that the isotopic composition of seasonal snowpacks in spring may still provide insight to understanding winter snow processes. Overall, our results present the potential seasonal snowpack isotope records for moisture source attribution in remote regions, where continuous sampling of precipitation or atmospheric vapor is not feasible.
This interdisciplinary study is laying the foundation for a completely new era in the Arctic earth sciences because it couples programs in Atmospheric circulation patterns, sea Ice traits, moisture transport pathways, Snow pack structure, and snow's ability to record these processes in its water Isotope properties (AISI). The degree to which these AISI interactions are stored in snow packs, and can be retrospectively recovered across the Arctic, will require a coordinated, synchronized pan-Arctic program of research that would allow validation and model development. AISI interactions would include sampling in permanent snow fields that records multiple years of winter and samples from glaciers on Svalbard, in Alaska, and on the Greenland Ice Sheet. These comparative studies would benefit from in situ water vapor isotope measurements [55,99,109] and event-based sampling of snow across the entire Arctic [110].
Author Contributions: All authors contributed intellectually to the framework, interpretation of the data and the manuscript writing and editing. Collected and processed the snow pit data in Imnavait Alaska, compiled the snow isotope dataset, led the isotope data analysis and manuscript writing, P. Data Availability Statement: Climate data for the study sites are available in the publicly accessible repository cited in the paper. Isotope data are available on request from the corresponding author.

Acknowledgments:
The authors acknowledge Juha Hatakka and the Finnish Meteorological Institute and staff working at the stations in Pallas-Yllästunturi National Park. Valtteri Hyöky helped maintain the Picarro vapor instrumentation at Pallas and assisted with the humidity-isotope calibrations. Tarja Törmänen and Aino Erkinaro assisted with sample analysis at the University of Oulu. Kaj Lynoe and Jari-Pekka Nousu assisted in snowpit field campaigns. This material is based in part upon work supported by the National Science Foundation through the NEON (National Ecological Observatory Network) Program. A special thanks to Chris Florian, Natchaya Durden, and Hongyan Luo for their support on NEON instrument quality control and data processing.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A.
This Appendix includes a description of the atmospheric stable isotope analyses at Imnavait and Pallas, in addition to Figures A1-A7 and Table A1 related to the supporting text and the main manuscript.
Water vapor isotope composition was measured as a part of the National Ecological Observatory Network (NEON). A Picarro L2130-i has been operated nearly continuously at the Toolik field station since the beginning of September 2017. The analyzer is connected to a network of tubing on the eddy tower that enables ambient air to be sampled at four heights-10, 20, 30, and 40 m. For the analysis presented in this study, we examined only the isotopic measurements from the top collection height of 40 m. By limiting data to this height, we attempted to obtain a best estimate of the isotopic composition of the free troposphere prior to (or after) its interaction with the snow pack. At regular intervals, the L2130-i switches between sampling mode, which samples the ambient air, and field validation mode, which uses a set of known liquid water standards to ensure consistency and accuracy in the measurements. During sampling mode, air is pumped for 9 min. One minute between each measurement is used for solenoid switching and purging the instrument, and for stabilization of the instrument. A Standards Delivery Module (SDM) regularly injects three water standards into the L2130-i to validate isotopic measurements that cover much of the range of natural variability. Standards for validating isotope measurements were: Low (δ 18  . The standards were run approximately every 24 h. A dry air system is connected to the SDM to dilute water vapor concentrations to replicate natural concentrations. Although instrument drift does not have a significant impact on the instrument measurements, it is well known that there is an instrument bias presented by the humidity of the air being measured, particularly at low humidity values [86][87][88]. Because of this effect, isotopic measurements must be corrected based on their humidity. NEON personnel conducted humidity bias experiments on the Toolik Picarro on 5 September 2019 using the Med standard (δ 18 O = −18.38‰ and δ 2 H = −136.52). The standard water was injected for 5 min into dry air of 13 different flow rates to measure standard values across a range from 600 to 24,500 ppmv. These water vapor concentrations were computed from the NEON output variable "rtioMoleWetH2o", which is the equivalent to the Picarro measured water vapor ppmv/1,000,000. The 1Hz measurements of H 2 O concentration, δ 18 O, and δ 2 H were aggregated to 1-min resolution, and the standard deviation of these 1 Hz measurements was also calculated for each 1-min window. Because it takes time for the instrument to stabilize following a change of the input gas, we used these standard deviation values to ensure a proper stabilization was reached. These criteria include the 1-min H 2 O concentration standard deviation being less than 25 ppmv and the change of the log 10 (H 2 O) ppmv from one minute to the next being less than 0.02. Data that meet those criteria are plotted in Figure A1, where we show the relationship between δ 18 O vs. log 10 (H 2 O) ( Figure A1a) and δ 2 H vs. log 10 (H 2 O) ( Figure A1b). Strongly significant relationships are observed for each isotope ratio, and it can be seen that there is significant bias in the low humidity measurements, whereas when the humidity decreases, the δ 18 O and δ 2 H measurements become artificially depleted. The 2nd order polynomial equation fitted in Figure A1 was applied to the ambient air measurements to remove the humidity bias. This experiment was conducted from~600 to 24,500 ppm. Although it is possible that the humidity bias changed over time, this pattern is quite common across Picarro devices, thus we assume that a correction based on this experiment is reasonable for the full dataset and does not introduce new sources of significant error (i.e., these errors are not larger than the humidity bias itself).
Based on examination of the data and the typical low water content of the air at Toolik (mean H 2 O concentration for this measurement period was 2200 ppm), we also exclude the first of the nine minutes of a given measurement height to account for extra necessary stabilization time. The 1 Hz isotopic measurements area initially aggregated to 1-min values. The humidity bias correction determined in Figure A1 was applied to this data, which has the effect of increasing δ 18 O and δ 2 H values at low H 2 O concentrations. Extreme d-excess outliers were excluded from this analysis at this point, including values less than −15‰ and greater than 70‰. Data that were clearly contaminated by the dry air system used for standard runs (i.e., isotopic values clearly outside the natural range) were also excluded from this analysis. This included two longer periods (21 December 2018-4 January 2019, 29 January 2019-1 February 2019) and several short-term intervals (i.e., several minutes/hours). USGS-46; δ O: −29.8 ± 0.03; δ H: −235.8 ± 0.70‰). Standards were measured for 12 min at each humidity level, and a non-linear regression was performed on the δX-mixing ratio relationship, where δX is either δ 18 O or δ 2 H. The nonlinear regression was of the form: where δXcorr is the difference between the observed isotopic value and the actual standard isotopic value, q is the water vapor mixing ratio, and a and b are constants. The best-fit curve was used to determine the isotope correction as a function of cavity humidity ( Figure A2). Second, ambient measurements were standardised to the Vienna-Standard Mean Ocean Water (VSMOW2) scale. Every 12 h, two USGS standards were automatically injected by the SDM at a constant flow rate and analysed for 20 min to calibrate and monitor for potential instrument drift. All ambient and standard measurements were corrected for the humidity concentration dependence, and based on the uncertainty of both corrections, the measurement accuracy was estimated at ± 0.2‰ for δ 18 O and ± 0.6‰ for δ 2 H. Measurement precision, estimated from the standard deviation of calibration measurements at a constant humidity level, was ± 0.3 and ± 0.8‰ for δ 18 O and δ 2 H, respectively, at humidity levels > 2000 ppmv.  After quality assurance and controls, this 1-min data was aggregated to a daily average. Typically, over 300 min of sampling that was distributed evenly across the day in the series of 9-min intervals was compiled in a given day, providing a relatively robust approximation of the daily average. The standard error of daily aggregated data was generally <0.2‰ for δ 18 O and <1‰ for δ 2 H and d-excess. The humidity correction was only calibrated to measurements down to 600 ppm, so isotopic measurements at H 2 O concentrations below that exhibit additional possible uncertainty. Given the strength of the relationship in Figure A1, it is likely that this uncertainty is still rather low compared to the natural variability, but it should be considered possible that there are unknown uncertainties at these low humidity time periods. Due to a number of instrument issues and the challenges of setting up the large number and spatially diverse NEON sites, standard runs did not occur consistently across the measurement period, particularly for the window of this study (October 2018-March 2019). However, when standards were run, it can be observed that the Picarro L2130-i was quite stable over the measurement period, with relatively small drifts between successful standard runs. These drift values were substantially smaller than the natural variations examined in this study.

Appendix A.2. Pallas, Finland
At Pallas, atmospheric vapour isotopes and mixing ratios were measured at the FMI Sammaltunturi station (67.973 • N, 24.116 • E; 565 m asl). The Picarro L2130-i and SDM has been installed inside the station building since mid-December 2017 and is connected to the FMI common sampling line. Ambient air was sampled via 12 m length stainless steel tubing (56 mm) connected to a composite polyethylene-aluminium tube extending the sampling line. Temperature inside the building was maintained at 20 • C throughout the measurement campaign, aiding instrument stability and preventing condensation in the line.
Stable isotope ratios (δ 18 O and δ 2 H) were measured approximately every second by the analyzer and calibrated using standard protocols employed at Toolik. First, we established the humidity-isotope response of the instrument in the field to correct for potential analytical bias [85,86]. Using an integrated dry air cylinder (HiQ zero air 5.0), the vapor stream of two standard waters (United States Geological Survey (USGS) 45 and 46) were measured across a range of controlled humidity injections from 400 to 10,000 ppmv that encompass the range of ambient humidity values during the winter 2018-2019 field campaign. The standard waters approximately span the isotopic range of ambient measurements during this period (USGS-45; δ 18 O: −2.2 ± 0.01‰; δ 2 H: −10.3 ± 0.40‰ and USGS-46; δ 18 O: −29.8 ± 0.03; δ 2 H: −235.8 ± 0.70‰). Standards were measured for 12 min at each humidity level, and a non-linear regression was performed on the δX-mixing ratio relationship, where δX is either δ 18 O or δ 2 H. The nonlinear regression was of the form: where δX corr is the difference between the observed isotopic value and the actual standard isotopic value, q is the water vapor mixing ratio, and a and b are constants. The bestfit curve was used to determine the isotope correction as a function of cavity humidity ( Figure A2). Second, ambient measurements were standardised to the Vienna-Standard Mean Ocean Water (VSMOW2) scale. Every 12 h, two USGS standards were automatically injected by the SDM at a constant flow rate and analysed for 20 min to calibrate and monitor for potential instrument drift. All ambient and standard measurements were corrected for the humidity concentration dependence, and based on the uncertainty of both corrections, the measurement accuracy was estimated at ±0.2‰ for δ 18 O and ±0.6‰ for δ 2 H. Measurement precision, estimated from the standard deviation of calibration measurements at a constant humidity level, was ±0.3 and ±0.8‰ for δ 18 O and δ 2 H, respectively, at humidity levels > 2000 ppmv.     Pallas pit six (see Figure 1). For the Imnavait hand hardness test, snow grain type and snow density data are also presented.  Pallas pit six (see Figure 1). For the Imnavait hand hardness test, snow grain type and snow density data are also presented. Figure A7. Snow pit stratigraphy and sampling intervals with associated isotope values at (a) Imnavait pit three and (b) Pallas pit six (see Figure 1). For the Imnavait hand hardness test, snow grain type and snow density data are also presented.