Holocene precipitation seasonality in northern Svalbard: Influence of sea ice and regional ocean surface conditions

Arctic precipitation is predicted to increase in the coming century, due to a combination of enhanced northward atmospheric moisture transport and local surface evaporation from ice-free seas. However, large model uncertainties, limited long-term observations, and high spatiotemporal variability limit our understanding of these mechanisms, emphasizing the need for paleoclimate records of precipitation changes. Here we use lipid biomarkers in lake sediments to reconstruct precipitation seasonality in northern Spitsbergen, Svalbard. We measured the hydrogen isotopic ratios (d2H) of n-alkanoic acids (C20 eC30) from sedimentary leaf waxes in lake Austre Nevlingen, Spitsbergen. We interpret d 2H values of mid-chain (C22) and long-chain (C28) n-alkanoic acids to represent d2H of lake and soil water, respectively. Austre Nevlingen lake water d2H reflects amount-weighted mean annual precipitation d2H. In contrast, soil water is mostly recharged by summer rainfall, and therefore reflects d2H values of summer precipitation. Austre Nevlingen leaf wax d2H values are 2H-depleted in the Early Holocene, suggesting high winter precipitation amounts. This coincides with high summer insolation, strong Atlantic water advection and reduced spring sea-ice cover in surrounding waters. Winter precipitation continued to dominate until c. 6 cal. kyr BP. After 6 cal. kyr BP, the trend in the biomarker record is not as clear. This could be related to colder conditions causing longer duration of seasonal lake-ice cover, thereby influencing the precipitation seasonality registered by the lake water. The Austre Nevlingen record suggests a close relationship between precipitation seasonality and regional ocean surface conditions, consistent with simulations suggesting that Arctic winter sea-ice loss will lead to increased local evaporation. © 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
By the end of the 21 st century, warming in the Arctic is predicted to exceed the global average by a factor of 2.2e2.4 (Collins et al., 2013). This amplified warming will affect the amount and seasonality of precipitation, through hydrological intensification (Rawlins et al., 2010;Collins et al., 2013;Bintanja and Selten, 2014). Increased high-latitude precipitation can be caused by 1) atmospheric circulation changes and enhanced poleward moisture transport, mechanisms mainly influencing summer precipitation (Dickson et al., 2000;Zhang et al., 2013), and 2) increased local surface evaporation, due to warmer Arctic seas and reduced winter sea-ice cover, mainly influencing fall and winter precipitation (Bintanja and Selten, 2014;Kopec et al., 2016). A warmer atmosphere also causes a larger fraction of Arctic precipitation to fall as rain (Førland et al., 2020). These hydrological changes will impact Arctic ecosystems, glacier mass balance, and infrastructure (Bintanja and Andry, 2017;Adakudlu et al., 2019). In order to better assess future Arctic precipitation changes, we need improved understanding of the mechanisms causing precipitation variation in the past (Sundqvist et al., 2014;Linderholm et al., 2018). In the absence of instrumental data, paleoclimate (proxy) data are necessary to fill this critical knowledge gap.
Robust records of Holocene precipitation amount and seasonality are lacking. Leaf wax hydrogen isotope ratios (d 2 H) provide means to evaluate hydrological changes over long time scales. The potential of the method for high-latitude lacustrine records has been demonstrated in Arctic Canada (Thomas et al., 2012), Greenland (Balascio et al., 2013;Thomas et al., 2016Thomas et al., , 2018, Svalbard  and northeastern Russia (Wilkie, 2012).
Leaf waxes are straight-chain hydrocarbon compounds produced by both terrestrial and aquatic plants. The n-alkanoic acid and n-alkane components of leaf waxes are well-preserved in sedimentary records due to their resistance to degradation (Eglinton and Calvin, 1967). Terrestrial plants predominantly produce long-chain wax compounds, whereas aquatic plants produce mid-chain wax compounds (Ficken et al., 2000;Meyers, 2003;Nichols et al., 2009;Gao et al., 2011). The d 2 H values of leaf waxes reflect the d 2 H values of the plant source water, with a biosynthetic fractionation that is largely constant for specific compounds (Sachse et al., 2012;McFarlin et al., 2019). The source water for terrestrial plants is soil water, which is mostly recharged by summer rainfall (Cooper et al., 1991;Throckmorton et al., 2016). The hydrogen isotope ratios of terrestrial leaf waxes (d 2 H terr ) therefore reflects summer precipitation isotope values, influenced by some evaporative enrichment of soil and leaf water (Kahmen et al., 2013). In contrast, the source water for aquatic plants is lake water, which may reflect summer or mean annual precipitation d 2 H values, depending on the residence time and source of lake water Thomas et al., 2020). Hence, leaf wax d 2 H can be used to reconstruct source water d 2 H and ultimately d 2 H of precipitation and other aspects of climate, including evaporation (Rach et al., 2017). The d 2 H of precipitation is influenced by changes in local and source temperature, moisture source location and transport history (Frankenberg et al., 2009), as well as evaporation for terrestrial waxes (e.g., Sachse et al., 2012;Thomas et al., 2016Thomas et al., , 2018. The hydrogen isotope ratios of aquatic plant leaf waxes (d 2 H aq ) reflect the isotopic composition of the lake water during the growing season. Therefore, the precipitation isotope seasonality recorded by d 2 H aq (annual mean or summer) depends on the local conditions determining the residence time of the lake water . Lakes with short residence times have growing season lake water d 2 H values biased towards summer precipitation d 2 H, whereas water isotopes from lakes with long residence times reflect amount-weighted mean annual precipitation d 2 H (Jonsson et al., 2009;Cluett and Thomas, 2020). Balascio et al. (2018) inferred that d 2 H aq changes in lake Hakluytvatnet ( Fig. 1), which has a short residence time (i.e., the lake water during the growing season reflects summer precipitation d 2 H), reflect air temperature and the relative influence of polar and sub-polar air masses. Here, we apply the same concept for another lake in northern Svalbard. By choosing a lake with a longer residence time (i.e., lake water isotopes representing mean annual precipitation d 2 H), we can extract information about seasonal variation in Arctic precipitation. We use d 2 H of C 22 and C 28 n-alkanoic acids to infer past lake water and leaf water d 2 H values, and reconstruct Holocene precipitation seasonality in northern Spitsbergen. By better constraining precipitation seasonality during earlier Holocene warm periods, we can improve our understanding of the mechanisms causing precipitation change in the past, present and in the future.

Regional setting
The Svalbard archipelago is influenced by strong climate gradients, with air masses transporting heat and moisture from the Atlantic regions meeting cold polar air from the Arctic Basin ( Fig. 1a; Førland et al., 2009Førland et al., , 2011Vikhamar-Schuler et al., 2019). South of Svalbard, the North Atlantic Current splits into the West Spitsbergen Current (WSC) and the North Cape Current (NCaC), transporting warm saline water along the west coast of Svalbard and eastward into the Barents Sea (Fig. 1a). The Barents Sea also receives colder, fresher water masses from the Arctic Ocean transported by the East Spitsbergen Current (ESC). During extreme cold events over the Barents Sea in winter, the strong temperature gradient between cold Arctic air and relatively warm Barents Sea surface water can cause convection and development of lowpressure systems (Rasmussen, 1985). Easterly winds caused by these systems bring moisture to eastern Svalbard, resulting in high precipitation (Hagen et al., 1993).
Since the establishment of the first permanent meteorological stations in Svalbard in the 1910s, temperatures have increased during all seasons, most notably in winter and spring (Førland et al., 2011). The trend in precipitation amount for the same period is not as clear. The robustness of instrumental precipitation data on Svalbard is hindered by 1) a lack of stations -especially in the central and eastern parts of the archipelago, and 2) measurement errors -notably due to snowdrift and undercatch (Førland and Hanssen-Bauer, 2000;Hanssen-Bauer, 2002;Humlum, 2002). However, instrumental records generally display positive trends in autumn and winter precipitation, and negative trends in summer and spring (Adakudlu et al., 2019). On an annual basis, all long precipitation amount time series from Svalbard (i.e., >40 yr) show positive trends, with precipitation increases between 2e8% per decade (Førland et al., 2011). Climate model projections for Svalbard agree that precipitation amount will increase in the future as temperature rises and sea ice retreats. Førland et al. (2011) projected an annual precipitation increase from 1961e1990 to 2071e2100 of less than 10% in southwestern Spitsbergen, but more than 40% in northeastern Svalbard. This regional climate model simulation was based on the MPIB2 scenario (Max Planck Institute global ECHAM4 model with SRES B2 emission scenario). Based on the RCP8.5 scenario ("business as usual"; , Adakudlu et al. (2019) projected a 65% annual precipitation increase in Spitsbergen from 1971e2000 to 2071e2100. In the same study, RCP4.5 (reductions after 2040, "medium emissions";  gives a median annual precipitation increase of 45%. Both model projections indicate largest precipitation increases in northeastern Svalbard, with greatest change in winter (þ54e90%) and least in summer (þ27e48%; Adakudlu et al., 2019). The fact that all models agree that the greatest precipitation increases will occur in the northeast, highlights the precipitation sensitivity of our study area. Model uncertainties and high spatial and temporal variability in precipitation emphasize the need for paleoclimate records to elucidate the mechanisms causing precipitation variability on Svalbard.

Study area
Lake Austre Nevlingen is located 2 km east of the 108-km-long Wijdefjorden, close to the mouth of the fjord (Fig. 1). Wijdefjorden is associated with the Billefjorden Fault Zone, separating Devonian sandstones and clastic sedimentary rocks in the west from Precambrian crystalline bedrock in the east (Dallmann, 2015). The landscape around Austre Nevlingen is characterized by undulating terrain with lakes in depressions ( Fig. 1b and c). The timing of deglaciation of the study area remains poorly constrained. Cosmogenic exposure ages from Tyrkampen (450 m above sea level (a.s.l.), 25 km to the south), suggest that this highland area was ice free by 18.3 ± 1.2 to 17.9 ± 1.2 kyr ago (Hormes et al., 2013).

Lake catchment
Austre Nevlingen (79 47 0 N, 15 47 0 E; Fig. 1) has a surface area of 0.13 km 2 , a maximum depth of 18 m, and sits at an elevation of 41 m a.s.l. This is below the local marine limit, which can be estimated to c. 60 m a.s.l. based on reconstructed isobases of uplift (Forman, 1990). The Austre Nevlingen catchment is c. 0.48 km 2 , sitting between 153 and 282 m a.s.l. summits (Fig. 1b). The catchment is dominated by frost-shattered bedrock and boulders with almost no soil cover, and dry non-vegetated to sparsely vegetated slopes. During the August 2015 and August 2018 field campaigns, we observed terrestrial vegetation dominated by graminoids (e.g., Luzula confusa) and mosses (e.g., Polytrichum). Submerged bryophytes, but no other aquatic plants, were observed growing in the and lake water d 2 H from Austre Nevlingen (2018). (C) Monthly average precipitation calculated for all available data from the meteorological stations Svalbard Airport, Ny-Ålesund and Isfjord Radio, for the period AD 1990e2019 (MET Norway, 2020). For each box plot, the middle line displays the median precipitation, the box represents the 25% to 75% quartile range and whiskers the maximum and minimum values. lake. The growing season is from mid-June to early September (ORNL DAAC, 2018) and the lake is estimated to be ice-covered from October to the beginning of July (Holm et al., 2012).
The lake is fed by surface runoff, but has no distinct inflow stream. There is also no distinct outflow stream on the surface. To estimate the proportion of lake water replaced annually and seasonally, we multiply the precipitation amount at Svalbard Airport during the different seasons (spring melt season ¼ October to June, ice-free season ¼ July to September) with the catchment area, and divide the product by the lake volume (Jonsson et al., 2009;Thomas et al., 2020). Svalbard Airport (Fig. 1d) has the most similar precipitation amount and seasonality to our study site. On average, approximately 11% of the lake volume is replaced each year: 7% during the spring melt season (i.e., by runoff representing precipitation from months with lake-ice cover) and 4% during the ice-free season (i.e., by runoff during the ice-free months). These values indicate a relatively long residence time, and that the lake water isotopes could be influenced by surface evaporation. Yet, isotope values of lake surface water collected in August 2018 indicate minimal evaporative enrichment. This surface water has a d 2 H value of À96.3 ± 0.2‰ and d 18 O of À13.2 ± 0.02‰ (Figs 2b and 3). The deuterium excess (d-excess ¼ d 2 H e 8 x d 18 O; Dansgaard, 1964) is 9.0‰. These values are well within the range of precipitation dexcess at Ny-Ålesund and Isfjord Radio (IAEA/WMO, 2019) and close to the GMWL (d-excess ¼ 10‰). The lack of isotopic evidence for evaporative enrichment suggests that the lake water may be completely flushed each year , perhaps via subsurface runoff into the lake and subsurface drainage out of the lake through the frost-shattered bedrock. Lichen-free boulders up to 0.5 m above the current lake level indicate that the lake may have experienced higher water levels in the past, which could indicate subsurface drainage during periods of thawed active layer.

Bathymetry and sediment core collection
The bathymetry of Austre Nevlingen was surveyed in August 2018, with a Garmin ECHOMAP™ Plus 73SV with a CV52HW-TM transducer and a 5 Hz receiver, using the Quick Draw contour function. Coring in 2015 was guided in real-time with depth sounding data from a Hondex PS-7 Transducer. To minimize disturbance, sediment cores were obtained from the central, deepest part of the basin (Fig. 1b). Coring was performed through a hole in the floor of a small zodiac, which was anchored in a stable position on the lake surface. A 70-cm-long surface core (ANS1) was collected with a Universal surface corer (120 cm long and 68 mm diameter coring tubes). A 129-cm-long piston core (ANP3) was retrieved using a hand-held lightweight piston corer (200 cm long and 60 mm diameter coring tubes). To record modern conditions in the catchment, lake water and leaves from a selection of the most common plants growing in the area were collected in August 2018. Leaves were sampled from multiple plants of the same species to get representative samples, and chosen to represent different plant groups, including rushes (Luzula confusa), dwarf shrubs (Salix polaris), terrestrial mosses (Polytrichum sp.) and aquatic mosses (unidentified).

Lithology and stratigraphy
The lake sediment cores were split open, logged and analyzed in the sediment lab and ITRAX core facility at the University of Copenhagen. ITRAX scanning was performed to record visual and radiographic imagery and variations in the element-geochemical properties throughout the cores. A rhodium (Rh) tube with a 1 mm resolution and 30 s exposure time was used for X-ray fluorescence (XRF; e.g., Kylander et al., 2011) and 4 mm resolution for magnetic susceptibility (MS; e.g., Sandgren and Snowball, 2002). Information on the ITRAX core scanner is given by Croudace et al. (2006). We present the Ti signal and the Ca/Fe ratio to show variation in minerogenic input, as these can be used to detect glacial meltwater input (Kylander et al., 2011) and lake basin isolation (Larsen et al., 2017), respectively. The Ti peak area was normalized against the incoherent (inc) and coherent (coh) Rh scatter (Ti/ (inc þ coh)) to remove instrumental effects, as suggested by Kylander et al. (2011).
After scanning, the lithology and stratigraphy of the cores were visually inspected and logged. In order to determine whether the lowermost minerogenic unit was deposited in a marine or lacustrine environment, we searched for marine microfossils (e.g., foraminifera) using a stereo microscope. The total organic content was determined through loss-on-ignition (LOI; Heiri et al., 2001). For this purpose, samples (1 cm 3 ) were collected every 2 cm, dried at 110 C for 24 h and ignited at 550 C for 4 h.

Lipid biomarker extraction and analysis
Biogeochemical and isotopic analyses were performed in the Organic and Stable Isotope Biogeochemistry Laboratory at the University at Buffalo, NY, USA. Samples (3e4 cm 3 ) for lipid biomarker analysis were collected every 5e8 cm from ANS1 and every 4 cm from the organic part (i.e., the uppermost 97 cm) of ANP3, resulting in thirty-six samples for the composite record. The samples were freeze-dried and homogenized, and free lipids were extracted with an Accelerated Solvent Extractor (ASE) 200 (Dionex) using dichloromethane (DCM):methanol 9:1 (v:v). An internal standard (C 20:1 n-alkanoic acid, Fisher Scientific, 4.2 mg) was added to the total lipid extract (TLE) after ASE extraction, allowing biomarker quantification. Preparation for n-alkanoic acid analysis followed previously published methods (Thomas et al., 2012). Hydrocarbons were separated by flash-column chromatography using NH 2 -functionalized silica gel, DCM:isopropanol 2:1 (v:v) as the neutral eluent and 4% acetic acid in DCM as the acid eluent. The acid fraction was methylated at 60 C overnight using 5% anhydrous HCl in methanol of known isotopic composition, and purified using silica gel columns, removing hydroxyl-carboxylic acid esters with hexane and recovering the fatty acid methyl esters (FAMEs) in DCM.
The FAMEs were quantified using a Thermo Scientific Trace 1310 gas chromatograph (GC) equipped with two AI1310 autosamplers, two split/splitless injectors, and two flame ionization detectors (FIDs) operated in parallel for higher throughput. For all analyses, the inlets were held at 250 C and operated in splitless mode for the first 45 s, after which split flow was turned on at 14 mL/min. Hydrogen carrier gas was used with a constant flow rate of 3.6 mL/ min. The oven program was as follows: initial temperature of 70 C held for 1 min, then a fast ramp of 27 C/min to 230 C, followed by a slower ramp of 6 C/min to 315 C, with a final hold at 315 C for 10 min. Compound masses were calculated via external calibration curves determined separately for each detector using a C 28 FAME standard (Fisher Scientific), and those masses were normalized to the mass of extracted sediment as well as the recovery of the internal standard to determine the concentration of waxes in sediment.
The d 2 H values of the FAMEs were measured using a Thermo Scientific Delta V Plus isotope ratio mass spectrometer (IRMS) coupled to a Trace 1310 GC via an Isolink II and Conflo IV. GC inlet flow settings and oven temperature program were identical to those used in GC-FID analysis except for the carrier gas, which was helium at a constant flow rate of 1.5 mL/min. The HTC reactor was held at 1420 C for all analyses. The H 3 þ factor was monitored regularly at the beginning of every sequence.
where C n is mg/g sediment of each n-alkanoic acid with n carbon atoms.

Chronology and core correlation
The chronology of the sediment cores was established based on accelerator mass spectrometry (AMS) 14 C measurements on four terrestrial and eight aquatic plant macrofossils. The macrofossils were identified and isolated from 0.5 mm sieving residues, and measured at the Ångstr€ om Laboratory, Uppsala University and Lund University Radiocarbon Dating Laboratory, Sweden. All radiocarbon ages were calibrated using the online OxCal software (v. 4.3;Bronk Ramsey, 2009) and the IntCal13 dataset (Reimer et al., 2013). Calibrated ages are presented in calibrated year before present (cal. yr BP; BP ¼ 1950).
For stratigraphic correlation between the cores, we used visual trends in the XRF data together with the radiocarbon age constraints. The cores were aligned in AnalySeries (v. 2.0.8; Paillard et al., 1996) using tie points in the elemental data. The Ca/Fe ratio was the main signal used (Fig. A, Appendix A), but multiple elemental ratios were compared to construct a common stratigraphic depth scale.
To establish an age-depth relationship for the composite core and incorporate age uncertainty into our biomarker records, we entered all proxy and chronology data into Linked PaleoData (LiPD) files (McKay and Emile-Geay, 2016). The data were analyzed in R (v. 3.6.3; R Core Team, 2020) using the GeoChronR package  and the IntCal13 dataset (Reimer et al., 2013). An agedepth model was constructed using a prior mean accumulation rate (acc. mean) of 100 years/cm and upper (d.min) and lower (d.max) depths set to the upper-and lowermost depths for the biomarker data (Fig. 4). The age ensemble was mapped to the paleo data, allowing us to display the biomarker data with age uncertainty (Figs 6 and 8) and plot the d 2 H records for both cores together (Fig. 8b). Fine lines in Figs 6 and 8 show the raw data plotted on the median of each age point, bold lines represent median values of the age model ensembles, and the light and dark shading show the 1s and 2s age uncertainty, respectively.

Chronology and core correlation
Radiocarbon ages, calibrated median ages and 2s age intervals for ANS1 and ANP3 are presented in Table 1 and Figs 4 and 5. The uppermost age in ANP3 (LuS 12222) appeared to be an outlier compared to the rest of the age-depth model, which was otherwise fairly linear, despite different terrestrial and aquatic sources of the samples. We therefore excluded this sample when running the final age-depth model (Fig. 4). Our biomarker record spans 11.5e2.2 cal. kyr BP for ANP3 and 9.2e1.9. cal. kyr BP for ANS1. Figure 5 shows optical and X-ray imagery, sedimentological logs, LOI and selected XRF data for the two cores. The cores contain three lithological units.

Lithology and stratigraphy
Unit L1 is only present in ANP3 (129e97 cm; Fig. 5), and consists of light grey, clayey-silty diamict with interlaminated organic material. The LOI is below 2%, and Ti/(inc þ coh) is higher (0.20e0.60) than in the other units, where values are generally below 0.04. Also, Ca/Fe is higher than in the overlying unit, especially in a coarser grained interval between c. 112e103 cm depth. L1 could be interpreted to represent minerogenic-rich sedimentation, driven by inflow of glacial meltwater, or deposited in a shallow glaciomarine setting. No foraminifera were observed in this unit. The transition from unit L1 to L2 at 97 cm depth, with decreasing minerogenic fines (silt and clay) in the lower 15 cm, indicates either an abrupt termination of glacial meltwater inflow or isolation from the fjord. Both these interpretations are supported by a sharp decrease in Ti/ (inc þ coh), and gradually decreasing Ca/Fe ratio.
Unit L2 (c. 97e34 cm in ANP3, 70e33 cm in ANS1; Fig. 5) consists of light brown to light olive grey stratified gyttja, with sporadic interbedded aquatic bryophytes and high organic content, representing accumulation of organic material with minimal minerogenic input. The organic content is higher in ANP3, increasing gradually from 13% at the lower boundary to 30% near 89 cm depth, and remaining between 30e40% for the rest of the unit. In ANS1, the same unit has LOI values between 17e27%. In both cores, Ti/ (inc þ coh) is generally below 0.03, and Ca/Fe ratio low and stable around 0.015e0.030.
Unit L3 (34e0 cm in ANP3, 33e0 cm in ANS1; Fig. 5) consists of dark brown gyttja, with a gradational lower boundary, and organicrich beds interbedded with abundant aquatic bryophytes. LOI continues to gradually increase to 37% in ANS1 and just below 50% in ANP3, in the uppermost part of the cores. Ca/Fe ratios are higher and more variable (0.025e0.155), whereas Ti/(inc þ coh) ratios show values similar to L2, although decreasing in ANS1.

n-alkanoic acid sources and productivity
Large isotopic differences between mid-chain (C 20 , C 22 and C 24 ) and long-chain (C 26 , C 28 and C 30 ) n-alkanoic acids, but similar values within these groups, suggests different sources of mid-and long-chain n-alkanoic acids (Fig. 6). This supports the hypothesis that Austre Nevlingen records both summer (terrestrial) and mean annual (aquatic) precipitation isotopes, and that hydrogen isotope ratios of C 22 and C 28 leaf waxes can be used as indicators of different plant sources (Gao et al., 2011;Thomas et al., 2016, in press) and to reflect different seasonality.
The analysis of modern plants suggests that different leaf wax sources (i.e., aquatic or terrestrial plants) have different chain length distribution. Luzula confusa contains the highest relative C 28 concentration, followed by Salix polaris (Fig. 7). These two terrestrial plant species also show the lowest concentrations of midchain waxes. The two mosses, Polytrichum sp. and an aquatic bryophyte, record the highest C 20 and C 22 concentrations, and relatively less C 28 . The leaf wax concentrations were highest in the Early and Late Holocene (Fig. 6), suggesting highest productivity during these periods. Long-chain leaf wax concentrations are higher than midchain concentrations throughout most of the record, reflecting dominance of terrestrial vegetation or higher leaf wax productivity for terrestrial plants. Generally, the sediment record is dominated by C 28 , C 26 and C 24 n-alkanoic acids, with relatively higher C 20 and C 22 concentrations in ANP3 than in ANS1 (Fig. 7).
The average chain length follows a similar trend as the leaf wax concentrations, with high values (>26) during periods of high C 28 concentration (Fig. 6). This correlation is weaker for ANS1, where the concentration of mid-chain waxes (especially C 24 ) generally is higher than in ANP3.

d 2 H of sedimentary n-alkanoic acids
Based on the evidence presented above, we interpret d 2 H C28 to reflect terrestrial plant wax d 2 H (d 2 H terr ) and d 2 H C22 to reflect aquatic plant wax d 2 H (d 2 H aq ). d 2 H terr is more stable throughout the Holocene than d 2 H aq (Fig. 8). The difference in d 2 H terr between the cores is generally less than 10‰ for most of the record, whereas d 2 H aq differs as much as 80‰ at times. When interpreting the d 2 H records for the two sediment cores separately, they suggest different stories.
ANP3 covers a longer period and extends further back in time than ANS1 (Fig. 8a) and À145‰ for most of the record, with a slightly decreasing trend.
In ANS1, d 2 H aq is relatively 2 H-enriched prior to c. 7.9 cal. kyr BP. On the contrary, d 2 H aq in ANP3 becomes 2 H-depleted during this period. d 2 H aq varies by around 40‰ (À300 to À260‰) for the rest of the ANS1 record, with values comparable to the 2 H-depleted intervals in ANP3. d 2 H terr is relatively 2 H-enriched prior to c.
6.2 cal. kyr BP, with values between À160 and À150‰, compared to À180 to À165‰ for the later part of the record.

Discussion
The isotopic composition of n-alkanoic acids in the Austre Nevlingen record reveals environmental and climatic changes throughout the Holocene. Even though we are able to correlate ANS1 and ANP3 based on the XRF data (Fig. A, Appendix A), the d 2 H records of each core differ significantly, especially for the mid-chain compounds. Despite the large d 2 H aq variability in the upper portion of ANP3, there are several data points with d 2 H values similar to ANS1, suggesting that this variability is real. The data do not indicate that one record should be favored over the other. Therefore, we discuss the d 2 H records for each core separately (Fig. 8a) and for the composite record (Fig. 8b). The composite record includes both the oldest sediments from ANP3 and the youngest sediments from ANS1, allowing us to discuss our results in full.

Influence of lake system dynamics on d 2 H
In Austre Nevlingen, shifts between relatively 2 H-enriched and 2 H-depleted d 2 H aq values during the Mid and Late Holocene ( Fig. 8a and b) suggest that the lake is close to an isotopic threshold. This means that the lake water isotopes could shift between recording summer and mean annual precipitation isotope d 2 H depending on processes acting in the catchment Thomas et al., 2020).  (Reimer et al., 2013). Details of each radiocarbon age are given in Table 1.
We suggest that one of the main mechanisms causing the high amplitude (>80‰) variability in Austre Nevlingen d 2 H aq could be the duration of ice on the lake. If lake-ice cover persists during the main spring snowmelt, most of the 2 H-depleted water could bypass the lake in these years (MacDonald et al., 2016). Furthermore, late lake-ice melt could delay the peak in primary production of the lake, and the 2 H-depleted winter precipitation isotope signal may not be incorporated into the leaf waxes. On the contrary, years without summer lake-ice cover would cause the aquatic plants to receive and incorporate more 2 H-depleted water. Another possible factor affecting the seasonal signal preserved by the leaf waxes is fluctuating lake level, affecting the residence time of water in the lake. Previously higher lake level (and therefore larger lake volume) is suggested by the light-colored 'bathtub' rim around the lake, implying a previously longer residence time. On the other hand, a lower lake level (and therefore smaller lake  Table 1), LOI (%), Ca/Fe and Ti normalized by the incoherent and coherent signal (Ti/(inc þ coh)). The XRF data are plotted as raw data and with a 25-point running average. volume) could result in a shorter residence time. If the lake was fully flushed by spring melt in some years but not in others (i.e., because of changes in lake volume), lake water, and therefore d 2 H aq , could have shifted between reflecting summer and mean annual precipitation d 2 H values. Such lake level shifts, although at a much larger scale, have been described from Hakluytvatnet (Fig. 1). Balascio et al. (2018) interpret a mid-Holocene (7.5e5 cal. kyr BP) hiatus in their sedimentary record to represent desiccation of the lake as a response to warm and/or dry conditions at the time. Today, Austre Nevlingen lake water is not affected by evaporative enrichment, as the d-excess measured in August 2018 (9.0‰) is within the range of precipitation d-excess. It is possible that evaporation had a significant impact on lake water d 2 H during earlier Holocene warm periods as suggested by Balascio et al. (2018), if the ice-free periods were longer and the summers were warmer. Another possible way to switch the lake water from a mean annual to a summer signal would be by subterranean drainage during periods of deepened active layer. This process could theoretically decrease the lake volume enough to flush the lake completely in summer.
Another explanation for the up to 80‰ difference in d 2 H aq between the Austre Nevlingen cores and between adjacent samples in a single core could be the relative abundance of aquatic bryophytes in the samples, since moss mats occur at different depths in the cores. If some samples mainly contain interbedded bryophytes rather than bulk sediments, d 2 H aq in these samples could be dominated by one source rather than a mix of different sources. However, we could not confirm any correlation between samples dominated by aquatic mosses and significantly 2 H-enriched or 2 Hdepleted d 2 H aq values, when revisiting the residues after lipid extraction. Furthermore, there are no distinct shifts in the relative chain length concentration associated with the shifts in d 2 H aq values that would suggest different plant sources (Fig. 7). Therefore, the presence of aquatic plant macrofossils does not explain the observed isotopic differences.
More variable d 2 H aq than d 2 H terr ( Fig. 8a and b) could reflect that the most prominent changes in precipitation have occurred during winter, since d 2 H aq in Austre Nevlingen is interpreted to represent a mix of summer and winter precipitation, whereas d 2 H terr represents summer rainfall. Additionally, d 2 H aq is relatively more sensitive to local variation within the lake, whereas d 2 H terr represents conditions integrated throughout the catchment.
We compare our isotope values from Austre Nevlingen ( Fig. 8a  and b) to those from Hakluytvatnet ( Fig. 8c; Balascio et al., 2018), and find some significant differences between the two records. The precipitation isotope seasonality recorded by the two lakes differ because of the large difference in residence time between the lakes. Hakluytvatnet is much smaller than its catchment , and therefore has a short residence time (i.e., it is likely fully flushed by spring melt and again by summer precipitation each year). The aquatic leaf wax d 2 H record from Hakluytvatnet should therefore be interpreted in terms of summer precipitation isotopes. This could explain the similar d 2 H patterns for C 25 and C 29 n-alkanes in that record (Fig. 8c). Austre Nevlingen has a longer residence time, and is not completely flushed by summer precipitation each year, so likely reflects the mean annual precipitation d 2 H. Because d 2 H for C 22 and C 28 in Austre Nevlingen represent different seasons, they show different patterns ( Fig. 8a and b). Larger isotopic variation at Austre Nevlingen compared to Hakluytvatnet could also reflect greater seasonality at our site, as it is situated farther from the main moisture source pathway along the west coast of Spitsbergen. Furthermore, the difference in analyzed leaf wax compounds (i.e., n-alkanes for Hakluytvatnet and n-alkanoic acids for Austre Nevlingen) could explain the different trends (Curtin et al., 2019).
The discussed factors highlight the need to interpret lake d 2 H records carefully, depending on the local conditions. For this specific system, we interpret d 2 H terr to reflect summer precipitation and evaporative enrichment and d 2 H aq to reflect mean annual precipitation, sometimes modified by changing lake dynamics. Mean annual precipitation isotopes mainly change as a function of moisture transport (Frankenberg et al., 2009). Even though the main signal preserved can be assumed to reflect climatic changes, the influence of various internal lake dynamics must also be taken into account.
5.2. Early Holocene: 11.7e8.2 cal. kyr BP Decreasing d 2 H aq and increasing d 2 H terr from the base of our d 2 H record at c. 11.5 to 9.5 cal. kyr BP ( Fig. 8a and b) coincide with a June solstice insolation maximum at 80 N ( Fig. 8f; Laskar et al., 2004) and rapidly increasing eastern Fram Strait surface and subsurface temperatures (Fig. 8e;Hald et al., 2007;Werner et al., 2016). Similar subsurface temperature trends are reported from the Barents Sea Margin (Risebrobakken et al., 2011). Hald et al. (2007, and Risebrobakken et al. (2011) interpret this warming as a response to increased influx of warm Atlantic water (AW) to the waters around Svalbard. This led to reduced spring sea-ice cover, evident by strongly decreasing PIP 25 indices in the eastern Fram Strait (Figs 1a and 8d; P B IP 25 and P D IP 25 , based on brassicasterol and dinosterol, respectively; Werner et al., 2016). Additionally, Allaart et al. (2020) suggested an Early Holocene decline in sea-ice cover in the nearby Wijdefjorden (Fig. 1c), based on IP 25 in marine sediments. Increased radiative forcing and greater AW heat transfer has also been suggested to cause warmer-than-present conditions in shallow waters around Svalbard between 11e9 cal. kyr BP (Mangerud and Svendsen, 2018) and peak warmth in three northern Svalbard lakes around 10 cal. kyr BP (van der Bilt et al., 2019).
Increasing Austre Nevlingen d 2 H terr between 11e9.5 cal. kyr BP supports Early Holocene summer warming, suggesting enhanced evaporative enrichment during summer and/or a more proximal source for growing season precipitation. A more proximal moisture source could be explained by reduced spring sea-ice cover and increased local evaporation. Balascio et al. (2018) also interpret increasing d 2 H (for C 25 to C 29 n-alkanes) to reflect the high summer insolation and greater influence of mild air masses from the south.
These processes could also explain the strong decreasing d 2 H aq during the same period. Mean annual temperatures would have to Fig. 7. Relative concentration distribution of even-chained C 20 eC 30 n-alkanoic acids in sediment cores ANS1 and ANP3, as well as selection of modern plants growing in the catchment. Fig. 8. Austre Nevlingen leaf wax hydrogen isotope data compared to selected regional Holocene climate records. For locations, see Fig. 1. Lines, shading and vertical colored bars shown as in Fig. 6. (A) Austre Nevlingen leaf wax d 2 H, C 22 (blue; aquatic) and C 28 (green; terrestrial) n-alkanoic acids for sediment cores ANS1 and ANP3. (B) Composite leaf wax d 2 H record for Austre Nevlingen. (C) Lake Hakluytvatnet leaf wax d 2 H, C 25 (blue) and C 29 (green) n-alkanes . Dashed lines in the mid-Holocene denote a possible hiatus, suggested by Balascio et al. (2018). (D) Eastern Fram Strait sea-ice proxies P B IP 25 (black) and P D IP 25 (grey) . (E) Eastern Fram Strait subsurface temperature based on planktic foraminiferal fauna assemblages (fine line; Werner et al., 2016), including 3-point running means (bold line; Husum and Hald, 2012). (F) June 21 st insolation at 80 N (Laskar et al., 2004). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) be more than 30 C cooler to explain the 100‰ depletion, considering a temperature-d 18 O relationship of 0.34e0.40%/ C (Kotlyakov et al., 2004), translated to d 2 H by multiplying by 8 (Dansgaard, 1964). More likely, the 2 H-depletion reflects a strong winter signal as a response to greater winter ocean evaporation due to the reduced sea-ice cover or more distal moisture transport during this interval.
These interpretations support an early Holocene Climate Optimum, which is in agreement with other recent studies in Svalbard Mangerud and Svendsen, 2018;Røthe et al., 2018;de Wet et al., 2018;van der Bilt et al., 2018van der Bilt et al., , 2019Voldstad et al., 2020), and elsewhere in the Arctic (Lecavalier et al., 2017;McFarlin et al., 2018). Our d 2 H record suggests that the regional Early Holocene warming had a strong effect on precipitation at the time, with local moisture from open seas leading to increased winter snowfall. This could also explain why some glaciers advanced in this otherwise warm period (e.g., Farnsworth et al., 2017Farnsworth et al., , 2018. After c. 9.5 cal. kyr BP, d 2 H terr stabilizes around À160‰, indicating stable summer conditions during this period. d 2 H aq continues to decrease until c. 8.5 cal. kyr BP, reflecting continued sea-ice decrease Allaart et al., 2020) and enhanced local evaporation.

Mid Holocene: 8.2e4.2 cal. kyr BP
In the Mid Holocene, d 2 H terr remains stable around À160‰ until c. 5.6 cal. kyr BP, indicating stable summer conditions. Relatively 2 H-depleted values after 5.6 cal. kyr BP suggest gradual cooling, with less evaporative enrichment and/or more 2 H-depleted summer precipitation during this period. Mangerud and Svendsen (2018) suggest Svalbard ocean water temperatures c. 4 C warmer than present between 8.2e6 cal. kyr BP, based on the presence of the marine bivalve mollusk Mytilus edulis, many of them along the northern coast. East Fram Strait subsurface temperatures stayed warm (up to 6 C) until c. 5 cal. kyr BP, with a slightly decreasing long-term trend ( Fig. 8e; Werner et al., 2016). Both Mangerud and Svendsen (2018) and Werner et al. (2016) interpret this warming to be attributed to intrusion of AW to the Svalbard shelf areas. High Mid Holocene temperatures are also inferred from lake sediment records, with minimal minerogenic accumulation, suggesting greatly reduced glacial activity or ice-free catchments (Svendsen and Mangerud, 1997;Røthe et al., 2015Røthe et al., , 2018. In Austre Nevlingen, 2 H-depleted d 2 H aq values suggest increased winter precipitation starting c. 9.5 cal. kyr BP and persisting into the Mid Holocene until c. 6 cal. kyr BP. An increase in winter precipitation could be explained by continuously low sea-ice extent. Low P D IP 25 and P B IP 25 values until c. 7 cal. kyr BP support this interpretation ( Fig. 8d; Werner et al., 2016). Higher winter precipitation during the Early-Mid-Holocene is also suggested by Røthe et al. (2018). They propose snowmelt dominated runoff to lake Vårfluesjøen between 10.2e7 cal. kyr BP, inferred from high frequency of 'snowmelt layers' in the sediment core. Vårfluesjøen is located on the western side of Wijdefjorden, c. 27 km WSW of Austre Nevlingen. After 7 cal. kyr BP, a lower frequency of snowmelt layers and the occurrence of aeolian sand in Vårfluesjøen, suggest drier Mid Holocene conditions (Røthe et al., 2018). Drier Mid Holocene climate is also proposed by Balascio et al. (2018), inferred from the hiatus (7.5e5 cal. kyr BP) in their sedimentary record. Our Austre Nevlingen record suggests that winter-dominated precipitation sustained until c. 6 cal. kyr BP, without any indications of desiccation of the lake. d 2 H terr express larger amplitude variability than in the Early Holocene and first part of the Mid Holocene ( Fig. 8a and b). We observe disagreement between the isotopic values in the two cores, but cannot explain exactly why. Some isotopic overlap between the cores (e.g., the 2 H-depleted d 2 H aq value in ANP3 c. 3.9 cal. kyr BP) supports that this large-amplitude variability is real. One potential interpretation is that this variability reflects changes in the precipitation seasonality recorded by the lake. As discussed in Section 5.1, Austre Nevlingen might be close to an isotopic threshold, meaning that the lake water incorporates relatively 2 H-depleted winter precipitation isotopes in some years, but not in others. This process could be explained by changing regional conditions and/or variable conditions in the catchment, amplifying the climatic signal. These processes could affect the seasonality of the isotopes registered by the lake in the Late Holocene d 2 H aq record. Hence, it is difficult to decipher the Late Holocene climate signal, except that the high-amplitude d 2 H aq variability suggests greater climate variability.
In our Austre Nevlingen record, slightly 2 H-depleted d 2 H terr values after 5.6 cal. kyr BP may be due to this regional cooling trend, as cooler conditions would cause less evaporative enrichment and more 2 H-depleted summer precipitation. Neoglacial cooling could also explain the more 2 H-enriched parts of the d 2 H aq record after 6 cal. kyr BP. Regional cooling and concomitant increases in sea-ice cover would have caused less winter precipitation, in turn causing the lake water to be biased to 2 H-enriched summer precipitation.

Conclusions
C Leaf wax hydrogen isotopes from Austre Nevlingen suggest large variability in precipitation seasonality on northern Spitsbergen throughout the Holocene. One of the strengths of our reconstruction is that we can extract both the summer and mean annual precipitation signal. This is possible because terrestrial and aquatic plants use different source water (soil water and lake water, respectively). Furthermore, the lake is not completely flushed by spring melt, with the long residence time allowing winter precipitation isotopes to be incorporated into the aquatic leaf waxes. Contrasting isotopic composition of mid-chain and long-chain n-alkanoic acids suggests different leaf wax sources for these two sets of compounds. C Most prominent Holocene precipitation changes occurred in winter, reflected in greater variability of the aquatic (mean annual) than terrestrial (summer) leaf wax hydrogen isotope ratios. C Early Holocene regional warming had a strong effect on moisture availability and precipitation seasonality in northern Svalbard. High summer insolation and strong Atlantic water influx contributed to reduced sea ice, which we suggest favored greater local winter evaporation leading to increasing winter precipitation. C Between 9.5e6 cal. kyr BP, the record is characterized by 2 Henriched terrestrial and 2 H-depleted aquatic hydrogen isotope ratios, which we interpret to reflect enhanced summer evaporation and/or 2 H-enriched summer precipitation, reduced sea-ice cover and local moisture contributing to high winter precipitation.
C After c. 6 cal. kyr BP, the d 2 H records show greater amplitude variability. We interpret these shifts to be influenced by lake system dynamics, such as the duration of ice cover on the lake, that amplify the climate signal. Hence, we cannot clearly distinguish the amplitude of climate variability in this portion of the record, although the lake variability is most likely in phase with climate variability. C Our results suggest that the precipitation seasonality in northern Spitsbergen is strongly linked to regional ocean surface conditions. As a result, the positive trend in winter precipitation observed in Svalbard today, may amplify because of warming ocean surface waters and reduction in sea ice.

Data availability
All chronological and geochemical data presented in this paper are freely available at the NOAA/World Data Service for Paleoclimatology website: https://www.ncdc.noaa.gov/paleo/study/29852.

Author contributions
AS, EKT, LH and SEK developed the idea; SEK, AS, WRF, LA, SB and OI conducted fieldwork; SEK, WRF, AS and LH described and interpreted the lithology; AS and LH identified plant macrofossils for radiocarbon dating; LA processed the bathymetry data; SEK, OCC and SD conducted organic and stable isotope biogeochemical work; SEK, EKT and NPM did the age modelling and leaf wax data analysis and visualization; all authors contributed with interpretation and discussion of results; SEK wrote the paper with contributions from all co-authors.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements
Fieldwork, radiocarbon dates, and laboratory analyses were funded by the Carlsberg Foundation (grant number CF14-0756 to AS), the Nansen Foundation (to AS), and the Svalbard Environmental Protection Fund (grant number 17/101 to AS and 17/114 to LA). The 2015 fieldwork was partly funded by the University Centre in Svalbard (UNIS) Research Fund (grant to OI). Laboratory analyses and technician support were also funded by National Science Foundation Grant (grant number 1652274 to EKT) and a UB Center for Undergraduate Research and Creative Activities Grant (to SD). We thank Sveinn Brynj olfsson and Sara Mollie Cohen for field assistance, Mike Retelle for providing field equipment, Marie-Louise Siggaard-Andersen for assistance with the ITRAX-scanning, and Alexandra Rouillard for macrofossil identification and data discussion. Finally, we thank Willem van der Bilt and one anonymous reviewer for constructive comments that improved this manuscript.

Fig. A.
Correlation between surface core ANS1 and piston core ANP3 from lake Austre Nevlingen, Svalbard. The records were aligned in AnalySeries (v. 2.0.8; Paillard et al., 1996), based on tie points in the Ca/Fe data.