Holocene Hydroclimate Variability in Central Scandinavia Inferred from Flood Layers in Contourite Drift Deposits in Lake Storsjön

Abstract

Despite the societal importance of extreme hydroclimate events, few palaeoenvironmental studies of Scandinavian lake sediments have investigated flood occurrences. Here we present a flood history based on lithological, geochemical and mineral magnetic records of a Holocene sediment sequence collected from contourite drift deposits in Lake Storsjön (63.12° N, 14.37° E). After the last deglaciation, the lake began to form around 9800 cal yr BP, but glacial activity persisted in the catchment for ~250 years. Element concentrations and mineral magnetic properties of the sediments indicate relatively stable sedimentation conditions during the Holocene. However, human impact in the form of expanding agriculture is evident from about 1100 cal yr BP, and intensified in the 20th century. Black layers containing iron sulphide appear irregularly throughout the sequence. The increased influx of organic matter during flood events led to decomposition and oxygen consumption, and eventually to anoxic conditions in the interstitial water preserving these layers. Elevated frequencies of black layer occurrence between 3600 and 1800 cal yr BP reflect vegetation changes in the catchment as well as large-scale climatic change. Soil erosion during snowmelt flood events increased with a tree line descent since the onset of the neoglacial period (~4000 cal yr BP). The peak in black layer occurrence coincides with a prominent solar minimum ~2600 cal yr BP, which may have accentuated the observed pattern due to the prevalence of a negative NAO index, a longer snow accumulation period and consequently stronger snowmelt floods.

1. Introduction

Numerous Holocene palaeoclimate studies have been performed in Scandinavia, many of which are based on lake sediment records [1,2,3,4,5,6,7]. The Scandes Mountains are particularly important for our understanding of climate variations in this region through their wealth of suitable sites along ecological gradients and their proximity to the climatically dynamic North Atlantic region [8,9,10,11,12,13,14]. However, previous studies have almost exclusively investigated sediments deposited in small and relatively shallow lakes, predominantly reflecting local-scale environmental dynamics. In contrast, sediment archives from large lakes with more extensive watersheds integrate across regional hydroclimate regimes, although their responses to environmental changes may be subdued due to longer water residence times and relatively insensitive biota [6,15]. A number of large and deep lakes are located on the lee side of the Scandes Mountains. Close to 20 of these lakes, which were formed by glacial scouring during repeated glaciations, exceed 100 km² in areal extent and many are >100 m deep. However, sediment records from this type of lake have remained largely unexplored as palaeoclimate archives in Scandinavia, although they have played an important role in the development of highly resolved records of Holocene climate change in continental Europe and other parts of the world [16,17,18,19]. In addition, the majority of previously published climate reconstructions from Scandinavia have focused primarily on temperature, while variations in hydroclimate, both gradual changes and extreme events like floods, have been partly overlooked, in spite of their societal importance [20].

Flood occurrences are non-stationary through time and the main driver of their frequency and intensity is the internal and external climate forcing on inter-annual to millennial time-scales [21,22,23,24,25]. However, our knowledge of flood occurrences in Scandinavia is limited due to the general lack of long palaeoclimate records reflecting such extreme events. Reconstructions of past flood occurrences would allow us to better understand their return period as a response to varying climate forcing [26,27,28,29] and may improve the anticipation of future flood risks.

During a flood, fine-grained detrital material is eroded in the catchment and transported in suspension into downstream lakes. When the transport capacity of the inflowing tributary stream diminishes in the water body, the terrigenous material is transported through hypo-, meso-, homo- or hyperpycnal flows and is finally deposited as a distinct detrital flood layer on the lake floor [30,31]. Detrital layers have been explored to reconstruct flood occurrences, but such studies are essentially limited to smaller lakes or sites in the vicinity of tributary rivers [25,26,28,32]. For larger lakes, currents associated with flood events or linked to specific wind regimes can modify the distribution of flood layers [33,34,35,36,37,38,39,40]. Hyperpycnal flows during flood events are strongly controlled by and tend to maintain, the lake floor morphology (i.e., channels). They can export sediment plumes containing river suspended load far away from the delta into deep basins. Lacustrine contourite deposits have been described in lakes with different basin shapes, depths and circulations patterns, e.g., in Lake Superior [41], East African Lakes [42], Lake Baikal [43], Lake Geneva [35] and Patagonian Lakes [38]. These drift deposits are either resulting from the influence of the Coriolis force on the displacement of sediment plumes at the lake floor, and/or from the interaction of hyperpycnal flood events (or turbidity currents) with currents at the lake floor generated by wind-induced internal waves or gyres.

Here we present a 4 m thick Holocene sediment sequence from contourite deposits in Lake Storsjön, a large lake on the lee-side of the Scandes Mountains in central Sweden. A particular focus is placed on black layers irregularly intercalated into the background deposition of grey, clay-rich sediments and their potential to serve as a proxy for past changes in the frequency of strong snowmelt floods. In addition, an ensemble of lithological, chemical and mineral magnetic proxies is used to reconstruct variations in catchment dynamics and in-lake processes, providing new insights into the deglaciation history of the site, millennial-scale climatic and environmental changes in the region during the Holocene and human impact on sediment deposition during recent centuries. Based on this multi-proxy approach, we investigate to what extent proxy variability in the sediments of Lake Storsjön relates to large-scale Holocene climate dynamics in Scandinavia.

2. Materials and Methods

2.1. Study Area

Lake Storsjön is located in central Sweden (63.12° N, 14.37° E; Figure 1a). Average annual temperature at the Frösön meteorological station, located on an island in the lake, is 2.5 °C, with average monthly temperatures ranging from −8.6 °C in January to 13.4 °C in July. The average annual precipitation is 491 mm, with >50% falling between June and September (Figure 1b; all values refer to the 1960–1990 reference period). The surrounding bedrock is mainly composed of greywackes, shales and limestones of Cambro-Silurian age, partly altered by the Caledonian orogeny, covered by a Quaternary clayey diamicton [44]. The catchment area of about 12,000 km² extends from the water divide of the Scandes Mountains at a maximum elevation of 1796 m to the outlet of the lake at 292 m a.s.l. (Figure 1c). It is dominated by coniferous forest (40%) and mires, with some cultivated areas on the clay-rich soils on carbonate bedrock. Lake Storsjön itself covers an area of 456 km², with 13% of the catchment area consisting of lakes.

The lake has a mean depth of 17.3 m, a maximum depth of 74 m and a volume of 8.02 km³. It can be divided into five principal parts (A to E; Figure 1d) according to its bathymetry and exposure to inlets and to the outlet. Due to the large catchment compared to the size of the lake, its water residence time of about 1 year is relatively short, with a more rapid throughflow in the northern parts of the lake, where the in- and outflows are located and a longer residence time in its southern basins [45].

The average ice cover period of Lake Storsjön lasts from mid-December to mid-May. Snowmelt leads to a distinct peak in river discharge in the catchment in May or June [46]. However, river management at >10 hydropower stations upstream of Storsjön’s main inlet masks the natural seasonal discharge pattern at the inlet. The lake has been regulated by a dam for hydropower production since 1938, involving yearly water level variations in the range of 290.5–293.3 m a.s.l. The development of a thermocline at 17–26 m water depth has been recorded in September in several parts (B, C and E) on a number of occasions between 2010 and 2014 [45]. Hydrodynamic simulations revealed significant surface currents associated with prevailing wind patterns during the ice-free season and distinct deeper currents linked to the lake floor morphology [45].

2.2. Seismic Survey

A seismic survey using a portable 4 kHz chirp digital device (Knudsen; Perth, Canada) and a conventional GPS (Garmin, Olathe, KS, USA) installed on-board a small fishing boat was performed during three days in November 2013 for mapping of the bedrock morphology and sediment geometry, as well as for identification of suitable coring sites. In total, 90 km of seismic profiles were acquired based on available bathymetric data [45], essentially in parts C and E of the lake, while strong winds and large waves prevented surveying the central basin.

2.3. Sediment Coring

Sediment cores were retrieved from the ice-covered lake in April 2014 in a 33 m deep basin with no major inlets in its vicinity. The seismic data at this location indicated a high sediment thickness and an undisturbed stratigraphy (Figure 2 and Figure 3). An 86 cm long surface sediment core was obtained using a freeze corer [47]. Seven overlapping 2 m long sediment cores, 8.3 cm in diameter, were obtained from three adjacent drill holes using a UWITEC piston corer (Mondsee, Austria).

2.4. Core Correlation

A continuous composite sediment profile of 7 m thickness (STO14) was constructed from the freeze core and the piston cores (supplementary material, Figure S1). The overlapping parts of the piston cores were aligned using 19 macroscopic lithological marker layers. Two of the cores were aligned at 440 cm based on their coring depths due to absence of clear stratigraphic markers. The freeze core was connected to the uppermost piston core based on six tie points identified in organic matter contents and element concentrations obtained from loss on ignition (LOI) and X-ray fluorescence (XRF) measurements, respectively (see below).

2.5. Dating and Age Modelling

Seven terrestrial macroscopic plant remains were identified under a binocular microscope and used for AMS radiocarbon dating at the LMC14 laboratory in Saclay (France) [48,49,50] and at the Radiocarbon Dating Laboratory at Lund University, Sweden (

Table 1. Radiocarbon dates from the STO14 sediment sequence. Samples were measured at the Laboratoire de Mesure du Carbone 14 (LMC 14), Université Paris-Saclay, France (Sac) and at the Radiocarbon Dating Laboratory, Lund University, Sweden (Lu).

Depth (cm)	Material	14C Age BP ± Error	Cal. Age BP (2σ Range)	Lab. ID
41.6	Wood fragment	1165 ± 30	985–1177	SacA39831
100.4	Twig	2330 ± 30	2213–2431	SacA39833
124.5	Wood fragment	2540 ± 30	2497–2747	SacA39832
180.2	Twig	3410 ± 30	3577–3811	SacA39834
229.8	Wood fragments	4950 ± 35	5603–5740	SacA39836
264.4	Twig	5915 ± 40	6657–6847	SacA39835
322.9	Wood fragment	7785 ± 45	8447–8639	LuS11056

). Measured radiocarbon ages were calibrated using the IntCal13 calibration dataset [51] and are expressed as calibrated years before 1950 CE (cal yr BP).

Marker horizons of pollution lead were identified using Pb concentrations determined by XRF (see below). Enrichment factors of Pb were calculated relative to a conservative lithogenic element (Zr) and a background value of Pb according to [52]. A site-specific Pb background was calculated using the part of the sediment sequence deposited prior to 2500 cal yr BP, i.e., before the earliest Pb pollution signals [53]. The anthropogenic Pb contribution was then calculated following [54]. An age uncertainty of ±50 years was assigned to the inferred Pb pollution dates to account for the sampling resolution of 0.5 to 1 cm.

Activities of ²¹⁰Pb and ¹³⁷Cs were measured on the upper 21 cm of the freeze core with an ORTEC High-Purity Germanium Gamma Detector (AMETEK Inc., Berwin, PA, USA) at the Department of Geology, Lund University. The ²¹⁰Pb activity was measured at its gamma peak at 46.5 KeV, the ²²⁴Ra was determined via its granddaughters ¹¹⁴Pb at 295 KeV and 352 KeV and ¹¹⁴Bi at 609 KeV. Self-absorption corrections were made for ²¹⁰Pb on each measurement [55]. The ¹³⁷Cs activity was measured at the 662 KeV gamma peak.

The age-depth model for the sediment sequence was created with Clam version 2.2 [56] using a smoothing spline with a smooth level of 0.3, including seven ¹⁴C dates, three pollution Pb dates and one ¹³⁷Cs date (²¹⁰Pb dating was not possible due to insufficient sample sizes). The sampling year 2014 (−64 cal yr BP) was assigned to the top of the composite profile.

2.6. Sediment Subsampling

The sediment sequence was subsampled for lithological, chemical and mineral magnetic analyses. A continuous series of sediment samples was extracted from the freeze core at 1 cm resolution from 1.5 to 10.5 cm depth (top sample 0–1.5 cm), at 0.5 cm resolution from 10.5 to 20 cm depth and at 1 cm resolution from 20 to 81 cm depth. The piston cores were sampled using cubic boxes (7 cm³) at 3 cm resolution, with the topmost sample centred at 1.5 cm depth. In addition to the regular interval samples, 14 samples were extracted from individual black layers, with a thickness between 0.4 and 0.7 cm.

2.7. Physical and Chemical Properties

Grey level values were extracted from a digital greyscale photograph (resolution 300 dpi) of the sediment sequence, along the vertical axis and averaged over a 2 cm wide horizontal window. Grey level values are expressed on a scale between 0 (black) and 255 (white).

Grain size distributions were determined for six black layers and six regular sediment (clay gyttja) samples of about 2.5 g dry weight. Samples were sieved with 250, 125 and 63 µm mesh sizes and each fraction was weighed with a precision balance. Three of the black layers and three of the clay gyttja samples were inspected at 15 kV with a Hitachi S-3400N (Tokyo, Japan) scanning electron microscope (SEM) at the Department of Geology, Lund University.

XRF and LOI were measured on the regular interval samples throughout the sequence and on the 14 black layer samples. XRF was measured using a Thermo Scientific Niton XL3t Goldd+ XRF analyser (ThermoFisher Scientific Inc., Waltham, MA, USA) in mining mode (Cu/Zn) for 180 s. All measurements were duplicated and the presented values are averages of the replicate measurements. Water and organic matter contents were determined by LOI using a Nabertherm muffle furnace (Lilienthal, Germany), following [57]. Samples were first heated to 105 °C for two hours and then combusted at 550 °C for four hours. The samples were weighed before and after each step.

2.8. Mineral Magnetic Properties

Mineral magnetic measurements were carried out on the regular interval samples at the Lund Paleomagnetic Laboratory. A first set of samples covering the composite sediment sequence was measured in September 2014 and a second set, comprising the overlapping parts of the piston cores, was measured in August 2016. Volume-specific magnetic susceptibilities (κ) were measured on all samples using a Geofyzica Brno KLY-2 kappabridge (Brno, Czech Republic). Anhysteretic remanent magnetizations (ARMs) were induced with a peak alternating field (AF) of 80mT with a direct current (DC) bias field of 50μT and progressively AF demagnetized at 5, 10, 15, 20, 40, 60 and 80 mT. Volume-specific susceptibility of ARM (κ_ARM) was calculated by normalizing the ARM intensity by the DC bias field and median destructive field of the ARM (MDF_ARM) determined as the AF required to remove 50% of the induced ARM. ARM measurements were not carried out for samples collected from the freeze core. Saturation isothermal remanent magnetizations (SIRMs) were induced using a Redcliffe model BSM-700 pulse magnetizer (Redcliffe Magtronics Ltd., Bristol, UK) with a DC field of 1T and measured using a Molspin Minispin magnetometer. Backfield IRMs were subsequently induced with a DC field of 100mT using a Molspin pulse magnetizer and measured using a Molspin Minispin magnetometer (Bartington Instruments Ltd., Witney, UK). S-ratios were determined by dividing the backfield IRM_-100mT with the SIRM. After completion of the magnetic measurements, the samples from the composite sequence were freeze dried in order to calculate dry density.

3. Results

3.1. Seismic Survey

Seismic profiles allowed mapping of the bedrock morphology corresponding to the acoustic substratum (Figure 2 and Figure 3). The complex morphology of the bedrock is typical for lakes of glacial origin, showing numerous over-deepened sub-basins delimited by sills or islands [58]. The bedrock crops out frequently at the lake floor (Figure 2) and is covered by unevenly distributed sediments constituting a succession of seismic units with contrasting acoustic facies. Seismic Unit 1 (SU1) is characterized by a varying thickness across the studied parts of the lake basin, ranging from 0 to 10 milliseconds two-way travel time (ms TWT) and consisting of transparent to chaotic acoustic facies with local high amplitude and continuous reflections (Figure 2). SU1 exhibits maximum thickness in Part C. Seismic Unit 2 (SU2) is characterized by a continuous and high-amplitude reflection and is generally thin (<1 ms TWT). In contrast, Seismic Unit 3 (SU3) is characterized by transparent acoustic facies (<10 ms TWT thick) and by the development of several low amplitude continuous reflections parallel to the lake floor that are locally developing onlap terminations on SU2, SU1 or the bedrock (Figure 2 and Figure 3). In addition, SU3 displays a channel-levee configuration typical of (patchy) drift or contourite deposits.

3.2. Lithological Description

Three main lithological units (I, II and III) were identified in the sediment sequence (Figure S1). These sedimentary units correspond to the seismic units (SU) identified in 0: Unit I corresponds to SU1, Unit II to SU2 and Unit III to SU3 (Figure 3). Unit I (715–399 cm depth) consists of dark grey, poorly sorted, sandy-silty diamictic sediments with abundant gravel and small stones (Subunits Ia, Ib and Id), interrupted by layers of predominantly clayey-silty sediments at two levels (Subunits Ic and Ie). Unit II (399–354 cm) is sharply distinguished from Unit I and is composed of finely laminated couplets of dark and light grey clay, with thicknesses between 1 and 2 mm. Unit III (354–0 cm) consists of brownish grey clay gyttja with occasional macroscopic terrestrial plant remains. Throughout Unit III, 335 black layers of 0.1 to 0.7 cm thickness appear at intervals between 0.2 and 3.5 cm (Figure 4). After the core opening, their colour changed gradually into reddish-brown. Vivianite concretions around macroscopic plant remains were identified visually (based on the bright blue colour upon exposure to air) within black layers at 153.5, 171, 180.5, 231 and 230 cm depth. Unit III is subdivided into three parts (Figure 5). Subunit IIIa (354–190 cm) is of grey colour and shows a steadily increasing organic matter content. Above 190 cm depth the sediment colour changes gradually from grey to brown, illustrated also by a decrease in grey level value (Figure 5) and the organic matter content remains relatively stable (Subunit IIIb). Subunit IIIc (43–0 cm) is defined by the onset of increasing concentrations in several elements (see below). Only Units II and III were analysed in detail for this study, as Unit I does not represent lacustrine sediments and lacks chronological information (see Section 4.1).

3.3. Age Model and Accumulation Rates

A distinct ¹³⁷Cs peak was observed at 5 cm depth (Figure 6a). Peak activities of ¹³⁷Cs in lake sediments originate from nuclear bomb tests in the 1960s and the explosion of the Chernobyl nuclear power plant in 1986 CE [59]. Large parts of the Storsjön catchment were affected by ¹³⁷Cs deposition after the Chernobyl accident [60]. It has been shown that ¹³⁷Cs from the Chernobyl fallout in varved sediment sequences from northern Sweden diffused downward and masked the nuclear bomb test peak [61]. Therefore, the observed ¹³⁷Cs peak in the STO14 sediments likely corresponds to the 1986 CE fallout and its presence indicates a well-preserved sediment surface.

Pb concentration peaks in lake sediments can be linked to well-dated events of anthropogenic Pb emissions and widespread, airborne pollution, e.g., at ca. 750 and 420 cal yr BP due to intensified metallurgy in Europe and around ‒25 cal yr BP due to the use and subsequent ban of leaded petrol [62,63]. The ‒25 cal yr BP peak can be clearly identified in the STO14 sequence. Two other Pb peaks are interpreted to reflect the 750 and 420 cal yr BP increases in anthropogenic Pb, which have also been identified in other Swedish lakes [53,64] (Figure 6b). Although these Pb peaks are not very distinct in the STO14 sediment sequence, the inferred dates are in good agreement with adjacent ¹⁴C and ¹³⁷Cs dates.

The radiocarbon dates are presented in Table 1. The age model indicates no age reversals or apparent hiatuses (Figure 7). As no plant macrofossils for radiocarbon dating were found below 320 cm depth, the age model was extrapolated below this point. Sediment facies like the finely laminated sediments found in Unit II are commonly interpreted as glaciolacustrine varves in Scandinavian lakes, with couplets of light and dark layers defined as one annual sedimentation cycle reflecting the seasonally varying glaciofluvial discharge in ice-proximal settings [65,66] (see also Section 4.1). Repeated counts of the dark-light couplets in Unit II yielded 250 (±5) varves and a mean annual sedimentation of 1.8 mm. Since the radiocarbon-based age model placed the boundary between Units II and III at 9550 ± 200 cal yr BP, the lower boundary of Unit II is tentatively assigned an age of 9800 ± 205 cal yr BP based on counting of the assumed varves.

Sediment accumulation rates are relatively constant from 9550 to 3500 cal yr BP (on average 0.3 mm/year), with elevated values between 3500 and 1500 cal yr BP (maximum 0.89 mm/year around 2600 cal yr BP). The highest accumulation rates (>2 mm/year) are found in the upper 2 cm of the sediment sequence.

3.4. Chemical Sediment Properties

Organic matter contents are relatively stable around 4.5% throughout Unit II (Figure 5). From the Unit II-III boundary the organic matter content increases steadily throughout Subunit IIIa to about 12.5% around 4000 cal yr BP. In Subunits IIIb and IIIc organic contents become more variable, fluctuating between 10% and 14% but there is no long-term trend in these Subunits. The highest organic matter content of 15% is found in the topmost sample (Figure 5).

Concentrations of Ti (Figure 8) and other lithogenic elements (Zr, Rb, K and Al; not shown) display largely similar patterns throughout the sediment sequence, with high values in Unit II (prior to 9500 cal yr BP), a constant decrease throughout Unit IIIa (to about 4000 cal yr BP), low values in Unit IIIb (between 4000 and 1100 cal yr BP) and an increase throughout Unit IIIc (from ca. 1100 cal yr BP to present). The average correlation between the lithogenic elements is 0.90 (p < 0.001) and they are all negatively correlated with organic matter contents (r = −0.82; p < 0.001). Si concentrations are not strongly correlated with the lithogenic elements (r = −25; p < 0.05), indicating that rather than clastic input, Si reflects diatom productivity [67]. Ca concentrations are highly variable in Unit II and reach up to 12,000 ppm. There is relatively little variation in Unit IIIa, with concentrations around 6100 ppm. In Unit IIIb concentrations are even lower (5500 ppm on average) until they increase again throughout Unit IIIc (Figure 8). Ca concentrations are negatively correlated with organic matter content (r = −0.40) and positively correlated with the lithogenic elements (mean r = 0.51). The transition from Unit II to III is also characterized by elevated Zr/K ratios (Figure 8), pointing to a coarser grain size [68], as well as by an increase in the water content of the sediments (Figure 5). Mo/Al ratios are low throughout Unit II, indicating oxic conditions at the sediment-water interface [69,70], then increase abruptly at the Unit II-III transition, followed by several minor peaks. From ca. 8000 cal yr BP onward, Mo/Al ratios increase slowly and steadily until another peak is reached at the top of the sediment sequence. S concentrations show similar patterns and are significantly correlated with Mo/Al ratios (r = 0.59; p < 0.001).

Overall, synchronous changes or breakpoints can be observed in several element concentrations and ratios at the Unit II-III boundary (ca. 9500 cal yr BP), at the Subunit IIIa-IIIb boundary (ca. 4000 cal yr BP) and at the Subunit IIIb-IIIc boundary (ca. 1100 cal yr BP). Concentrations of several elements, including Pb (Figure 6), Cu, S and P (Figure 8), increase rapidly in the topmost sediment samples.

3.5. Mineral Magnetic Properties

Variations in dry densities and κ are more or less synchronous, showing high values at the bottom (Unit II), which gradually decrease towards the top (Unit III) (Figure 9). This is broadly interpreted as a decrease in the minerogenic content of the sediments, in agreement with Ti concentrations (Figure 8) and the opposite trend in LOI (Figure 5). Magnetic parameters sensitive to the ferrimagnetic content (κ_ARM and SIRM) do not follow the same trend, which suggests changes in the mineral magnetic assemblage, e.g., due to the introduction of a new detrital or biogenic magnetic component. S-ratios, which provide a measure of the relative amounts of high (“hard”) vs. low (“soft”) coercivity of remanence magnetic minerals, show only subtle variations throughout the sediment sequence. In general, values approaching –1 correspond to a softer magnetic assemblage, typical of (titano)magnetites, while higher values suggest that part of the signal is carried by magnetically harder (antiferrimagnetic) minerals such as haematite or goethite.

Three samples near the base of Unit IIIa show unusually high SIRM values (and high S-ratios), which may at least partially be explained by measurement errors. These measurements, conducted on the first set of samples which have since been used for other analyses, were not reproduced in any of the second set samples from overlapping cores. Relatively high SIRM/κ ratios (15–35 kAm⁻¹) of these three samples could suggest the presence of authigenic greigite [71]. Greigite, which is a ferrimagnetic iron sulphide, dissolves when exposed to oxygen, which could explain the failure to reproduce these measurements in the second sample set. Peaks in sulphide concentrations (Figure 8) coinciding with the depths of these samples provide further support for the presence of greigite.

The ratio of κ_ARM to SIRM, often used as a grain-size indicator [72,73], varies between ca. 0.4 and 0.8 mm/A. The lowest values, corresponding to coarser magnetic grain sizes, are observed in Unit II while Unit III is characterized by generally smaller magnetic grain sizes. MDF_ARM, which is a coercivity sensitive parameter, shows a sharp increase at the boundary between Unit II and Unit III and then gradually increases up to 45 mT in Unit IIIb, indicative of the presence of bacterial magnetite [74]. Overall the κ_ARM/SIRM and MDF_ARM values are consistent with a simple system consisting of a mixture two end-members as defined by [74]: (i) a detrital component (κ_ARM/SIRM ≈ 0.5 mm/A, MDF_ARM ≈ 30 mT), which is essentially Unit II; and (ii) a soft biogenic component (κ_ARM/SIRM > 0.8 mm/A, MDF_ARM ≈ 45 mT) that becomes gradually more important in Unit IIIa and Unit IIIb. In Unit IIIc there is a relative increase in the contribution of the detrital component again. High values of κ_ARM/SIRM in the lower part of Unit IIIa are not satisfactorily explained by this simple model and could potentially indicate the presence of a third magnetic component (e.g., greigite, see above).

κ_ARM/SIRM ratios around 2–3 mm/A have been shown to be diagnostic of magnetosome chains produced by cultured magnetotactic bacteria [75], although lower values have also been reported for occurrences in natural sediments [76]. It has been found that κ_ARM/SIRM ratios of biogenic magnetite are particularly low under anoxic conditions [74], which may explain the generally modest values recorded here.

3.6. Characterization and Temporal Distribution of Black Layers

The three black layer samples and the three background gyttja samples examined by SEM each showed similar characteristics. Black layers are composed of dominantly minerogenic particles and aggregates, whereas the clay gyttja contains more autochthonous biogenic material, e.g., diatoms (Figure 10). The Storsjön sediments are generally fine-grained with ca. 95 weight % belonging to the silt and clay fractions. The black layers are characterized by a slightly larger grain size than the clay gyttja (Figure 11). The grain size proxy Zr/K [68] also indicates a larger average grain size of the black layers compared to bulk sediment samples (Figure 11).

The content of organic matter is slightly higher in the black layers (14%) than in the bulk sediment samples (12%) (Figure 11). Ti (Figure 11) and other lithogenic elements (Zr, Rb, K and Al; not shown) display broadly similar average values and ranges in the black layers as compared to the bulk sediment. Si concentrations are higher in the bulk sediment, in agreement with higher diatom abundances. The elements which show the largest relative differences between black layers and bulk sediment are Fe (60% higher in the black layers) and P (109% higher). S concentrations and Mo/Al ratios are comparable in both types of samples, while As shows about 30% higher concentrations in the black layers. Note that element concentrations in the black layer samples are compared to bulk sediment samples including black layers. The differences in element concentrations between black layers and the surrounding clay gyttja could therefore be larger than the presented values.

The frequency of occurrence of the 335 black layers is variable throughout the sediment sequence (Figure 12). There is a slight increase around 6000 cal yr BP and a marked augmentation between 3600 and 1800 cal yr BP, with a peak around 2600 cal yr BP. The second increase is also characterized by a clear increase in sediment accumulation rate (Figure 12). Throughout the sediment sequence the average thickness of the black layers does not vary significantly.

4. Discussion

4.1. Deglaciation and Early Postglacial Lake Development (Lithological Units I and II)

In the Jämtland region, the last deglaciation proceeded well into the early Holocene. According to the age model, lacustrine sedimentation in Storsjön began at 9800 ± 205 cal yr BP (Unit I-II boundary), with a dominant glacial influence on the lake until 9550 ± 200 cal yr BP (Unit II-III boundary). This temporal framework is in general agreement with previously published estimates of the regional ice recession chronology (reviewed e.g., by [77]), which are based on radiocarbon dates from the basal parts of peat deposits, varve-dated shifts in lacustrine sedimentation [78], or cosmogenic nuclide dates [79].

As the ice sheet retreated from the Storsjön area, large bodies of stagnant ice remained. The landscape is characterized by dead-ice topography and numerous small lakes and the occurrence of ice-dammed lakes has been well documented [78,80,81]. Ice-dammed lakes form when ice margins block the natural drainage of a catchment, leaving behind geomorphological evidence such as shorelines, perched deltas, or overflow channels [77]. The shorelines at Frösön Island, 5 km NE of our sampling site, point to the formation of an extensive ice-dammed lake as the ice sheet was breaking up [78]. The general seismic stratigraphy of the lake, with outcropping bedrock on the lake floor and sediments of variable thickness and acoustic facies, is consistent with this interpretation.

Unit I of the STO14 sediment sequence is dominated by diamictic sediments of glacial origin. The intermittent occurrences of fine-grained, partly laminated sediments (Units Ic and Ie, Figure S1) likely represent temporary phases of glaciolacustrine sediment deposition. However, as these stratigraphic intervals lack organic matter, the timing and duration of these temporary lake phases could not be determined through radiocarbon dating. Due to their presumably different deposition rates as compared to the postglacial lacustrine sedimentation, an extrapolation of our age model through Unit I would be inappropriate. The seismic unit SU1, corresponding to Unit I, shows several acoustic facies with highly variable thicknesses and acoustic amplitudes (Figure 3). These are typical for diamictic deposits in large lakes of glacial origin [82,83,84].

At the transition from Unit I to II, a lake basin began to form, although the sedimentation was still strongly influenced by glacial activity in the catchment. Without stabilizing vegetation, strong fluvial, glacial and glaciofluvial erosion led to high accumulation rates following the deglaciation. The sediments of Unit II are dominantly minerogenic, with negligible organic matter content, high concentrations of lithogenic elements (Figure 5 and Figure 8) and a mineral magnetic signature consistent with a predominantly detrital source (Figure 9). The finely laminated clays of Unit II are interpreted to represent glaciolacustrine varves, although microsedimentological analyses and independent chronological information for the bottom and top of Unit II would be needed to confirm the annual nature of the laminae. Dark-light layer couplets (varves) of glaciolacustrine sediments form in proglacial lake environments due to seasonal variations in erosion, transport and deposition processes. Light layers represent periods of enhanced meltwater discharge during the warm season, transporting coarser minerogenic material in suspension from the catchment into the lake. Dark layers are composed of very fine clay particles deposited during winter when the meltwater flux into the lake is clearly reduced and within-lake currents diminish under ice cover [66,85,86]. The varves in Unit II are preserved despite the supposed oxic conditions inferred from low Mo/Al ratios, because the dominant detrital sedimentation and cold temperatures prevented the establishment of burrowing biota. Unit II corresponds to the very high amplitude seismic facies of seismic unit SU2, characteristic for sediments deposited by glacial meltwater supply in an ice-contact lake [82,83].

4.2. Lake and Catchment Processes During the Holocene (Lithological Unit III)

At the transition from Unit II to Unit III (9550 ± 200 cal yr BP) the laminated clay is replaced by a more organic-rich clay gyttja, which is broadly similar to the present-day sediment facies. This transition coincides with major changes in chemical, mineral magnetic and acoustic parameters, which are linked to substantially altered in-lake processes and sedimentation conditions. An increase in grain size is indicated by increasing Zr/K ratios [68], possibly corroborated by an increase in water content (Figure 5 and Figure 8). The increase in grain size may be related to the cessation of the direct influence of glacial activity in the basin, which can be assumed to have led to primarily clay deposition at the distal position represented by our coring site and replacement with mineral matter of predominantly fluvial origin from the wider catchment.

The sharp drop in Ca concentration at the top of Unit II also reflects changes in catchment erosion processes influencing the sediment source. The dominant source of Ca is detrital input rather than biogenic precipitation. Throughout Unit II, glacial erosion and runoff on the newly exposed bedrock contributed clastic carbonate from the calcareous bedrock, which is mainly found S and E of the lake. In contrast, during deposition of Unit III, the lake received more input from the main inflow and the catchment NW of the lake, which has little carbonate bedrock. Furthermore, glacial activity provided fresh mechanically eroded bedrock material in Unit II, whereas Unit III is dominated by fluvial erosion of soils where carbonates were already dissolved.

The Unit II-III transition is also marked by abrupt increases in S concentration and Mo/Al ratios (Figure 8) and by the potential formation of the magnetic iron sulphide greigite (high SIRM/κ; Figure 9). The Mo/Al peak can be explained by lowered oxygen levels in the interstitial water [69,70], caused by the increased primary production (LOI increase, Figure 5) and consequent oxygen consumption by bacterial respiration. This in turn forced primary biota to reduce sulphate as an energy source, a process with high affinity for Mo enrichment [87]. The subsequent decrease in the Mo/Al ratio therefore represents the return of oxygenated interstitial water conditions, perhaps due to increasing circulation in the basin.

Following the Unit II-III transition, the period of ca. 9550 to 4000 cal yr BP (Subunit IIIa) is characterized by a steady decrease in the concentrations of lithogenic elements, which is probably mainly a dilution effect, as organic matter shows a parallel increase (Figure 5). The input of terrestrial organic matter into the lake started to rise with the progressive development of soils and vegetation in the catchment, likely advancing to higher altitudes during this period [11,88]. Furthermore, warming temperatures led to increased organic productivity in the lake, as reflected by the probable introduction of magnetotactic bacteria (increase in κ_ARM/SIRM and MDF_ARM; Figure 9) and to an increase in Si concentrations (Figure 8) through elevated biogenic silica production, predominantly from diatoms, a sediment component identified in the SEM images (Figure 10).

Subunit IIIb (ca. 4000 to 1100 cal yr BP) is characterized by largely constant lithogenic element concentrations, organic matter content and mineral magnetic parameters and by the absence of Mo/Al peaks, which indicates relatively stable environmental conditions.

The transition to Subunit IIIc at ca. 1100 cal yr BP marks the onset of increasing trends in several lithogenic elements (Zr, Rb, Ti, Ca, K, Al) as well as an increase in magnetic grain sizes (detrital source) (Figure 8 and Figure 9), which most likely reflects human impact on the lake related to slowly expanding agricultural activity on fertile soils in the vicinity. Little is known about the settlement and agrarian history of the Storsjön area, but historical data indicate an initial expansion phase in the 12th century CE, followed by widespread farm abandonment in the 14th century CE [89]. Pollen data from a site about 100 km to the south do not show any major agricultural expansion before the 16th century CE [90], in general agreement with enhanced input of lithogenic matter in our sediment sequence at about this stage. A more intense phase of human impact on the lake ecosystem was initiated in the 20th century CE as evident from rising organic matter content and an associated, distinct increase in sediment accumulation rate (Figure 4). These changes were likely related to sewage delivery following urbanization and the use of fertilizers, as indicated by the recent increase in P concentration (Figure 8). This development led to the return of hypoxic conditions in interstitial waters at the coring site for the first time since the early Holocene, as evidenced by rising Mo/Al ratios and S concentration, mounting further evidence of the last decades as a period of anomalous anthropogenic perturbation of previously largely unaffected aquatic ecosystems, even in relatively remote and resilient settings [91,92]. The strong increase in Cu concentration in the recent past (Figure 8) reflects airborne metal pollution related to mining and metallurgy [93] and broadly follows the trend of increasing Pb concentration (Figure 6).

4.3. Palaeoenvironmental Significance of the Black Layers

4.3.1. Sediment Deposition at the Coring Site

The sediments in Unit III at our coring site constitute the levee of a drift deposit located northward of a 200 m wide channel where seismic unit 1 (SU1) is outcropping (Figure 2 and Figure 3), indicating that sedimentation is controlled by deep-water circulation processes [94]. Drifts and moats in lakes, or lacustrine contourite deposits and channels, have previously been documented by seismic profiles and sediment cores in large lakes exposed to strong winds, e.g., Lake Superior in North America [41], Lago Cardiel in Argentina [39] and Lake Geneva in Switzerland [35].

Large tributary flood events can develop long-lasting hyperpycnal flows in both marine and lacustrine environments [31,95]. In lakes, the resulting hyperpycnal flood deposits are layers of variable thickness (millimetres to centimetres), essentially composed of clastic material and rich in organic matter of terrestrial origin [26]. In proximal deltaic environments, the deposits are graded and silty, whereas in deeper distal basins, they are commonly sharp-based and very fine-grained (clayey-silty) [23,31,96]. Hyperpycnal flood deposits in large lakes are frequently documented in deltaic environments and in distal fan lobes from the last deglaciation, or from Holocene cold and wet periods, such as the Little Ice Age, associated with alpine glacier fluctuations and enhanced glaciofluvial regimes [19]. In the case of Storsjön, snowmelt floods could carry fine-grained suspended sediment load from the catchment and develop long-lasting hyperpycnal flows when entering the lake basin. These turbulent flows follow the lake floor morphology and are exposed to the Coriolis force. If they occur frequently, they can control the distribution of clastic sediments in the lake and form channel-levee systems, such as the one we cored at site STO14. Hence, we hypothesize that the observed black layers represent distal snowmelt flood layers deposited in lacustrine contourites.

4.3.2. Mechanisms of Black Layer Formation and Preservation

The formation of black layers in sediments linked to low oxygen concentrations in bottom waters or surface sediments has been documented in a number of contexts. Black bands of iron sulphide have been observed in sediments of the Great Lakes deposited during the Holocene Thermal Optimum at ca. 8000 to 4000 cal yr BP [97]. These layers are thought to represent events of hypoxic or anoxic conditions caused by climate warming that led to intensified stratification and the absence of the seasonal turnover. In sediments along the southern coast of Finland, black layers containing iron sulphides have been identified, corresponding to anoxic and reduced bottom water conditions due to temperature and salinity stratification [98]. Anoxic conditions in surface sediments can also be caused by enhanced sedimentation after discharge events, when rapid burial prevents exchange processes at the sediment-water interface [99].

The black layers in the Storsjön sediments reflect brief and transient events of increased fluvial inflow and sediment transport into the lake, caused by strong snowmelt floods. These meltwater floods generally occur in May/June in the region and can cause soil erosion and increased sediment accumulation and grain size, as shown by studies of varved sediments in Scandinavian lakes [100,101,102]. The larger grain size of the black layers as compared to the surrounding clay gyttja is only minor, probably because of the distal location of the coring site, 30 km from the lake’s main inflow (Figure 1). It could be argued that the black layer itself is not the flood layer but rather that it was present at the sediment surface before the spring snowmelt flood. A rapid burial by flood deposits would reduce the exchange with bottom water, leading to rapid oxygen depletion and preservation of the black colour. This would lead to the higher black layer frequency during times of higher accumulation rates. However, this explanation is contradicted by the larger grain size and the absence of diatoms and other macroscopic organism remains in the black layers (Figure 10), which is also reflected by the low Si concentrations (Figure 11). Together with the higher organic matter contents, this also suggests that the black layer material originates from the catchment rather than from re-distribution of previously deposited sediments by lake currents. We therefore assume that the black layer itself represents a rapid, event-like deposition of allochthonous material.

Snowmelt floods lead to an increased delivery of terrestrial organic detritus from the catchment and perhaps augment the oxygenation of the lake water and consequently the primary production in the lake. The increased accumulation of organic matter is initially followed by decomposition and, due to the oxygen consumption by bacterial respiration, eventually leads to hypoxic or anoxic conditions in the interstitial water. Therefore, not all organic matter is decomposed, leading to the relatively high organic matter content of the black layers (Figure 11).

The black colour of the layers, indicating the presence of finely distributed pyrite (FeS₂) formed under anoxic conditions, changed to brown due to oxidation of FeS₂ on opening the cores. Before that, the elevated organic matter content kept the interstitial water anaerobic. The significantly higher As concentration of the black layers is an indicator of pyrite [103,104]. The anoxic conditions in the black layers are further confirmed by concretions of vivianite [(Fe₃(PO₄)₂·8H₂O)] [105,106,107], which are visible around macroscopic plant remains in some thicker black layers. Fe and P, both constituents of vivianite, also exhibit significantly higher concentrations in the black layers. Note that the supposed anoxic conditions occur in the sediment and not in the water. Cold temperatures and wind likely prevent the establishment of stable stratification and related anoxic bottom water conditions in Lake Storsjön.

4.3.3. Relations to Climate and Vegetation Changes

The processes of erosion and transport of material in the drainage basin that lead to the deposition of black layers depend on the amount of snowfall during winter and the rate of temperature increase in spring, as well as, in the longer term, the properties of the land surface [102]. Variations in black layer occurrence therefore depend on land cover (vegetation) changes, flood frequency (climate), or a combination of both.

The type of land cover determines the intensity of catchment erosion during snowmelt floods. A decrease in vegetation cover would lead to destabilized soils and consequently enhanced soil erosion even if the magnitude of meltwater floods does not change. A gradual change to cooler and/or wetter conditions has been observed in Scandinavia at the end of the Holocene Thermal Maximum and the onset of the Neoglacial Period around 4000 cal yr BP, with a thicker snow cover in the Scandes Mountains that melted later in the year [11,14]. This climatic change caused a tree line descent in the Scandes Mountains [3,108,109] including the Storsjön catchment [11,88], as well as an increased formation of peatlands [110]. A lowering of the tree line has likely destabilized soils at high elevations in the Storsjön catchment and increased erosion during snowmelt floods, triggering the mechanisms of black layer formation described above. Consequently, an increase in black layer frequency is observed from about 3500 cal yr BP onward.

Additionally, black layer frequency may increase with a more frequent occurrence of strong snowmelt floods, caused by enhanced snow accumulation in the winter, possibly accompanied by a faster warming and melting of snow in the spring as well as spring rainfall. This was the case e.g., during the 1995 spring flood, the largest in central Sweden in the 20th century: A long winter caused snow accumulation in mountain regions that was 20–50% above average and a rapid increase of temperatures in the spring, accompanied by rainfall, led to the onset of snowmelt almost simultaneously in mountains and lowlands [111]. Such an intense meltwater pulse would cause stronger erosion and trigger the processes described in Section 4.3.2. The factors controlling the amount of snow accumulated during winters and the magnitude of subsequent spring meltwater floods are complex. While a generally warmer climate could increase the total amount of precipitation, it would at the same time shorten the period of the year when precipitation falls as snow, as demonstrated for the River Ångermanälven catchment north of the Storsjön catchment [112]. Lake sediment records from southern Norway reveal a co-variability between winter precipitation and floods on millennial time scales but climate projections indicate a decrease in flood frequency with future warming, as this will mitigate that net effect of a wetter climate regime [113].

The phase of maximum snowmelt flood frequency inferred from the Storsjön sediments (Figure 12) coincides, with a grand solar minimum around 2750–2550 cal yr BP. Shifts in atmospheric circulation associated with the solar forcing caused abrupt regional climate change [114], which has been detected in a number of European palaeoclimate records, e.g., [115,116,117,118,119,120]. Scandinavian climate is dominantly controlled by large-scale atmospheric circulation dynamics in the North Atlantic region [12,121,122]. For example, variations in the North Atlantic Oscillation (NAO) index explain 55% of the variance of streamflow in Norway, with a positive NAO state leading to warm and wet conditions [123]. Model simulations suggest that solar radiation minima lead to the predominance of a negative NAO phase, the driver of cold winters in northern Europe [124]. A reconstruction based on lake sediments from Greenland suggests that negative NAO conditions characterized the period from 2800 to 2600 cal yr BP [125], broadly consistent with similar evidence from North Central Sweden [14] and expanding mountain glaciers in southern Norway [126].

Against the background of increasing availability of erodible soil material brought about by a tree line descent in the Storsjön catchment from 3500 cal yr BP onward, the effects of the solar minimum and associated atmospheric circulation changes may have caused the peak in flood layers around 2600 cal yr BP, with relatively cold winters in northern Europe leading to a longer period of snow accumulation and a higher intensity of subsequent snowmelt floods.

5. Conclusions

We have investigated the distribution as well as the chemical and mineral magnetic properties of sediments in Lake Storsjön and explored their palaeoclimatic significance. The bedrock topography and geometry of the basin infill were mapped for the first time. It is important to understand the highly uneven distribution of contourite deposits in this large lake for future sedimentological studies. The processes of deglaciation and early postglacial lake development inferred from the Storsjön sediment record are coherent with other geological evidence from the region. Lake Storsjön began to form about 9800 cal yr BP and significant glacial influence persisted in the catchment for another 250 years. Throughout most of the Holocene, geochemical and mineral magnetic proxies indicate relatively stable environmental conditions in the lake. Since 1100 cal yr BP, human influence becomes visible in the sediments through unprecedented proxy variation in a Holocene perspective, with potentially significant impact on the lake ecosystem during the 20th century.

Hyperpycnal flows following snowmelt flood events led to the deposition of contourites and the frequent formation of black layers containing iron sulphide through the enhanced input of terrestrial organic detritus. The processes of erosion and transport of material that lead to the deposition of black layers depend on two main factors: the amount of snow accumulated during the winter and the properties of the land surface. The increased occurrence of black layers between 3600 and 1800 cal yr BP likely reflects vegetation changes, as a tree line descent in the catchment destabilized soils and provided increased amounts of erodible organic soil material. The peak in black layer occurrence about 2600 cal yr BP coincides with a grand solar minimum. The associated predominance of the negative phase of the North Atlantic Oscillation caused relatively cold winters in northern Europe, increased snow accumulation and more intense snowmelt floods. The flood history reconstructed from the sediments of this large lake therefore reflects not only local catchment-scale processes but also large-scale climatic change.