A Holocene history of climate, fire, landscape evolution, and human activity in northeastern Iceland

Paleoclimate reconstructions across Iceland provide a template for past changes in climate across the northern North Atlantic, a crucial region due to its position relative to the global northward heat transport system and its vulnerability to climate change. The roles of orbitally driven summer cooling, volcanism, and human impact as triggers of local environmental changes in the Holocene of Iceland remain debated. While there are indications that human impact may have reduced environmental resilience during late Holocene summer cooling, it is still difficult to resolve to what extent human and natural factors affected Iceland's late Holocene landscape instability. Here, we present a continuous Holocene fire record of northeastern Iceland from proxies archived in Stóra Viðarvatn sediment. We use pyrogenic polycyclic aromatic hydrocarbons (pyroPAHs) to trace shifts in fire regimes, paired with continuous biomarker and bulk geochemical records of soil erosion, lake productivity, and human presence. The molecular composition of pyroPAHs and a wind pattern reconstruction indicate a naturally driven fire signal that is mostly regional. Generally low fire frequency during most of the Holocene significantly increased at 3 ka and again after 1.5 ka BP before known human settlement in Iceland. We propose that shifts in vegetation type caused by cooling summers over the past 3 kyr, in addition to changes in atmospheric circulation, such as shifts in North Atlantic Oscillation (NAO) regime, led to increased aridity and biomass flammability. Our results show no evidence of faecal biomarkers associated with human activity during or after human colonisation in the 9th century CE. Instead, faecal biomarkers follow the pattern described by erosional proxies, pointing toward a negligible human presence and/or a diluted signal in the lake's catchment. However, low post-colonisation levels of pyroPAHs, in contrast to an increasing flux of erosional bulk proxies, suggest that farming and animal husbandry may have suppressed fire frequency by reducing the spread and flammability of fire-prone vegetation (e.g. heathlands).

Overall, our results describe a fire frequency heavily influenced by long-term changes in climate through the Holocene. They also suggest that human colonisation had contrasting effects on the local environment by lowering its resilience to soil erosion while increasing its resilience to fire.

1 Introduction

Iceland is highly sensitive to most mechanisms controlling the evolution of Holocene climate in the North Atlantic, from millennial (e.g. shifts in deep-water formation and ocean current positions) to sub-decadal timescales (e.g. variability of the North Atlantic Oscillation) (Harning et al., 2021; Mjell et al., 2016; Moossen et al., 2015; Petit et al., 2020). Recent lake sedimentary records in Iceland (Alsos et al., 2021; Geirsdóttir et al., 2009a, 2013, 2019, 2020; Harning et al., 2016, 2020; Hiles et al., 2021; Larsen et al., 2011, 2012; Richter et al., 2021) draw a comprehensive picture of Icelandic environments during the Holocene (last 11.7 kyr). These Holocene paleoclimate reconstructions derived from lake sediments in Iceland show first-order millennial trends that reflect orbitally driven changes in Northern Hemisphere summer insolation and millennial to sub-millennial changes that are primarily impacted by northern North Atlantic ocean circulation and in part by local volcanism (e.g. Flowers et al., 2008; Geirsdóttir et al., 2013, 2020; Harning et al., 2018b; Larsen et al., 2012). These Holocene climate reconstructions further indicate a major shift from occasional to increasingly severe landscape instability and soil erosion occurring at least 300 years before the acknowledged settlement of Iceland (ca. 870 CE; The Book of Icelanders “Íslendigabók”, by Ari Thorgilsson, 12th century CE, e.g. Smith, 1995), suggesting that human impact had a secondary role to climate by lowering the resilience of the environment to an already ongoing naturally driven erosion (e.g. Bates et al., 2021; Geirsdóttir et al., 2009b, 2020). The ability to generate high-resolution Holocene terrestrial climate records, along with Iceland's relatively short settlement history, makes Iceland an ideal location to attempt disentangling the impact of natural climate variability and human activities on the changes in the local landscape during the late Holocene.

In this study, we use multiple organic proxies from a Holocene sediment core from the Stóra Viðarvatn lake in northeastern Iceland to investigate the effects of natural and anthropogenic drivers on the local Icelandic environment. First, we focus on tracing the evolution of fire regimes using pyrogenic polycyclic aromatic hydrocarbons (pyroPAHs; Lima et al., 2005). Fires can have a significant impact on ecosystems, affecting vegetation patterns, nutrient cycling, and wildlife habitat (e.g. Goldammer and Furyaev, 1996). The frequency, intensity, and spatial extent of fires can provide insights into past climate and environmental conditions (e.g. Marlon, 2009; Power et al., 2008), and to our knowledge there are no such records for the Holocene in Iceland, while limited data are available for the surrounding regions (Chen et al., 2023; Marlon et al., 2013; Segato et al., 2021; Zennaro et al., 2014). Second, as fire frequency can be influenced by human activities as well (e.g. Marlon et al., 2009, 2013; Zennaro et al., 2015), we also analyse faecal markers of human presence (Vázquez et al., 2021). By analysing these biomarkers from deglaciation to present, we can define their natural, pre-settlement background levels and thus potentially trace anthropogenic impact on the local environment, pinpointing human arrival in the lake catchment.

Finally, by coupling fire and human presence biomarker records with established proxies for environmental change (e.g. soil erosion and primary productivity, e.g. Argiriadis et al., 2018; Geirsdóttir et al., 2013; Gross, 2017), we test what control natural and/or human factors had on the evolution of the Holocene landscape in Iceland.

2 Methods

2.1 Study site

Stóra Viðarvatn (SVID) is a lake (2.6 km² surface area) located in NE Iceland (Fig. 1a and b) at an elevation of 151 m a.s.l. SVID has a maximum depth of 48 m, a catchment area of 17 km² (including the lake surface), and a volume of ca. 3.6×10⁷ m³ (this study, based on data from the National Land Survey of Iceland, Landmælingar Íslands, 2023; Axford et al., 2007). The lake is surrounded by Quaternary age basaltic lavas and glacial hyaloclastites formed by subglacial eruptions, as well as some Holocene soil with several millimetre-to-centimetre thick tephra layers (Hjartarson and Sæmundsson, 2014). The nearby Raufarhöfn station (Icelandic Meteorological Office, 2022) provides weather data for the 1961–1990 CE interval: the mean annual temperature is 2 °C with a maximum in July–August (8 °C), while the lake surface is usually frozen between November and March; the mean annual precipitation is 733 mm a⁻¹, with the lowest values occurring in May (28 mm) and the highest in October (ca. 86 mm), suggesting a lake water residence time between 5 and 9 years.

Figure 1 (a) Study area in NE Iceland. (b) Location of the Raufarhöfn climatological station and Ytra-Áland site (Karlsdóttir et al., 2014), which are 20 km NNW and 13 km ESE from Stóra Viðarvatn (SVID), respectively. (c) Location and catchment area of the Stóra Viðarvatn lake: the 20-SVID-02 core is marked by a pink hexagon, and an older 04-SVID-03 core is marked by a red square (Axford et al., 2007). SVID bathymetry (10 m isolines) is reported by Axford et al. (2007). Watershed catchment and contour lines (10 m) are calculated via ArcGIS (Esri, 2023) based on digital elevation models provided by the National Land Survey of Iceland. The basemap is sourced from Esri.

In February 2020, we recovered a 8.93 m long core 20-SVID-02 (66.236867° N, −15.837837° E; Fig. 1c) from 17.4 m water depth near the centre of the lake (Harning et al., 2023). The sediment was retrieved in seven drives of ∼150 cm each, using lake ice as a coring platform. The core was subsequently stored in a cool room (4 °C) at the Institute of Arctic and Alpine Research, University of Colorado Boulder, where it was subsampled. Previously, two studies have analysed an 8 m long core (04-SVID-03; Fig. 1c) retrieved in February 2004 to trace Holocene temperature (Axford et al., 2007) and δ¹⁸O from chironomid remains, as well as the δD, δ¹³C, δ¹⁵N of total organic matter (Wooller et al., 2008) at a 1–0.2 kyr resolution.

2.2 Tephrochronology

Our sediment core chronology takes advantage of the geochemical fingerprints of visible Icelandic tephra layers and their correlation to marker tephra of known age. A total of 13 tephra layers were sampled along the vertical axis, sieved to isolate glass fragments between 125 and 500 µm, and embedded in epoxy plugs. At the University of Iceland, individual glass shards were analysed on a JEOL JXA-8230 electron microprobe using an acceleration voltage of 15 kV, beam current of 10 nA, and a beam diameter of 10 µm. The international A99 standard was used to monitor for instrumental drift and maintain consistency between measurements. Tephra origin was then assessed using major oxide compositions, following the systematic procedures outlined in Jennings et al. (2014) and Harning et al. (2018a). Briefly, based on SiO₂ wt % vs. total alkali (Na₂O + K₂O) wt %, we determined whether the tephra volcanic source is mafic (tholeiitic or alkalic), intermediate, and/or rhyolitic. From here, we objectively discriminate the source volcanic system through a detailed series of bi-elemental plots produced from available compositional data on Icelandic tephra. Source eruption was then determined using the geochemical fingerprint and relevant stratigraphic information. See the Supplement for complete major oxide compositions and bi-elemental plots.

Table 1 Marker tephra layers of known age identified in 20-SVID-02 and used to develop the age model.

Using the 13 marker tephra layers of known age (Table 1), we generated a Bayesian age model using the R package rbacon 2.5.7 (Blaauw and Christeny, 2011; R Core Team, 2020) and default model functions (Fig. 2). We used the “slumps” function for the thicker tephra layers (e.g. Hekla 3 and Hekla 4) to reflect their instantaneous deposition on geologic timescales.

Figure 2 Stóra Viðarvatn age model generated in Bacon (Blaauw and Christeny, 2011). Horizontal green lines denote the age and uncertainty of marker tephra layers, the red line reflects the mean values of model iterations, the grey lines denote the 95 % confidence envelope, and darker shading reflects more likely ages. Grey vertical bars mark the “slumps” used for the Hekla 3 and Hekla 4 tephra layers.

2.3 Sample preparation and analysis

At the University of Colorado Boulder, we retrieved a total of 196 sediment core samples at an average spacing of 4.5 cm, providing a temporal resolution of decadal to centennial timescales. We freeze-dried samples for 24–48 h and ground and homogenised them (mean weight 1.5 g, range 0.6–6.6 g) using an agate mortar and pestle. Using 13–70 mg of sediment, we measured total carbon (TC), total nitrogen (TN), and δ¹³C (relative to VPDB) on an elemental analyser linked to a Thermo Delta V isotope ratio mass spectrometer (EA-IRMS) in the Earth Systems Stable Isotope Laboratory at the University of Colorado Boulder; samples were analysed against a suite of secondary laboratory standards that are extensively calibrated to international standard reference materials to correct for size, blank-mixing, linearity, and drift effects (Harning et al., 2018b). We did not decalcify samples as the contribution of inorganic carbon to TC in Icelandic lake sediment is considered to be negligible (see Sect. 5.1; Geirsdóttir et al., 2020). We analysed 9–11 mg of sediment for biogenic silica by diffuse reflectance Fourier transform infrared spectrometry (FTIRS) on a Bruker Vertex 70 with a Praying Mantis diffuse reflectivity accessory (Harrick) and report values in FTIRS absorbance units (e.g. Harning et al., 2018b).

We processed 86 selected samples for organic biomarker analyses. We extracted 0.4–2 g of dry sediment with an accelerated solvent extractor (Dionex ASE350). Sample size, matrix, and richness in organic compounds greatly varied along the core and initial tests occasionally showed colouration in the sample extract even after 4–5 cycles, likely due to the persistence of organics. Thus, to maintain consistency in lipid yields throughout the core, we used dichloromethane (DCM) : methanol (MeOH) 9:1 for six cycles of 5 min (static time), 100 °C, and 1500 psi. After extraction, we spiked the total lipid extract (TLE) with 1000 ng of 3-methyl-heneicosane (CAS no. 6418-47-9, Sigma-Aldrich), 20 ng of p-terphenyl (CAS no. 92-94-4, TCI), and 50 ng of pregnanol (5β-Pregnan-3α-ol, CAS no. 4352-07-2, Steraloids) as internal standards for the quantification of n-alkanes, PAHs, and sterols, respectively. We concentrated the TLE under a gentle flow of nitrogen and then mixed it with HCl-activated copper shots to remove elemental sulfur as copper sulfide precipitates. We then filtered the samples through a NaSO₄-packed Pasteur column to remove any residual water and copper sulfide and concentrated them under N₂. We subsequently separated the TLE into six chromatographic fractions using a Pasteur pipette packed with silica gel (60–200 µm − 60 A) and solvents of increasing polarity. We calculated the column's dead volume (DV) with n-hexane, and then eluted samples with 1.5 DV of n-hexane (F1), 2 DV of n-hexane : DCM 4:1 (F2), 1.5 DV of DCM (F3), 2 DV of DCM : acetonitrile (F4), 1.5 DV of acetonitrile (F5), and 3 DV of MeOH (F6). We derivatised fraction F4, containing the sterols and stanols, using TMS-BSTFA (Supelco) and pyridine (50:50) at 70 °C for 15 min, then dried it under N₂ and redissolved it in n-Hexane. We added 1 ng of p-terphenyl D₁₄ (CAS: 1718-51-0, Sigma-Aldrich) and 50 ng of 5α-cholestane (CAS: 481-21-0, Sigma-Aldrich) to fractions F2 (PAHs) and F4 (sterols), respectively, as injection standards to check the recovery and quantification consistency of analyses.

We analysed the n-alkanes, PAHs, and sterols using a Thermo Scientific Trace 1310 gas chromatograph (GC) equipped with a programmable temperature vaporiser (PTV) inlet and a Restek glass liner interphase to a TSQ8000-Evo triple quadrupole mass spectrometer (MS). We used a 60 m DB1 column (DB-1MS, 0.25 mm, 0.25 µm film thickness, Agilent, USA) to separate n-alkanes, a DB-5 column (DB-5MS, 0.25 mm, 0.25 µm film thickness, Agilent, USA) for PAHs and sterols, and He (1.2 mL min⁻¹) as a carrier gas. For n-alkane analysis, we injected samples in splitless mode at 65 °C, and the PTV was ramped to 400 °C at 3 °C s⁻¹ and held for 5 min. The GC oven temperature was programmed from 60 to 220 °C (25 °C min⁻¹) and then to 315 °C (2.5 °C min⁻¹, held 13 min). The n-alkanes were analysed in full scan (50–600 $m/z$) using the following MS conditions: 300 °C EI source at 70 eV electron energy, 50 uA emission current, and 15 V electron lens voltage, with a transfer line at 315 °C. For PAH analysis, all samples were injected in splitless mode at 45 °C, and the PTV was ramped to 400 °C at 11.6 °C s⁻¹ and held for 2 min. The GC oven temperature was programmed from 60 °C (held 1 min), to 150 °C (40 °C min⁻¹), and to 320 °C (3 °C min⁻¹, held 15 min). MS conditions were as follows: 250 °C EI source at 70 eV electron energy, 50 uA emission current, and 15 V electron lens voltage, with a transfer line at 320 °C. For sterol and stanol analysis, all samples were injected in splitless mode at 90 °C, evaporated at 100 °C (0.1 min), and the PTV was ramped to 400 °C at 8 °C s⁻¹ and held for 1 min. The GC oven temperature was programmed from 80 °C (held 1 min), to 200 °C (20 °C min⁻¹), and to 320 °C (5 °C min⁻¹, held 20 min). MS conditions were the same as for n-alkanes. PAHs and sterols and stanols were analysed in selected reaction monitoring (SRM) using the collision energies and mass transitions reported in Tables A1 and A2.

2.4 Analysis of air parcel back-trajectory patterns

To define the potential regional extent of airborne PAHs arriving to SVID's catchment area, we traced the back-trajectory of air parcels using HYSPLIT (hybrid single particle Lagrangian integrated trajectory; Draxler et al., 1998; Stein et al., 2015). Using a modified version of an R script originally developed to trace precipitation patterns (Caves Rugenstein and Chamberlain, 2018), we analyse data from the NOAA Global Data Assimilation System (GDAS; resolution 1° by 1°) at a 6 h frequency tracing back-trajectories for 3 d (72 h) and 2 weeks (336 h) during 2 years characterised by opposite North Atlantic Oscillation (NAO; Hurrell et al., 2003) configuration (2009–2010, NAO−; 2013–2014, NAO+; NOAA, 2023). PAHs deposition, which is enhanced by low temperatures, occurs not only via precipitation but also in dry conditions (Arellano et al., 2018; Feng et al., 2017; Golomb et al., 2001; Halsall et al., 2001). Thus, we present data for air parcel trajectories that did and did not produce precipitation within 6 h from the endpoint (SVID), initialising the trajectories at four different altitudes: 1000, 1500, and 2000 m a.s.l. (water vapour usually advects within an altitude of 2 m) (Bershaw et al., 2012; Lechler and Galewsky, 2013; Wallace and Hobbs, 2006) and 150 m a.s.l. (SVID surface elevation).

3 Background on proxies

3.1 Polycyclic aromatic hydrocarbons (PAHs)

We use pyrogenic PAHs (pyroPAHs) as tracers for the frequency and intensity of fire episodes and the PAH perylene as a biogenic PAH related to terrestrial organic matter input. PAHs are semi-volatile compounds that can be of pyrogenic, petrogenic, or biogenic origin (Kozak et al., 2017; Lima et al., 2005). Low-molecular-weight (LMW; see Table A1 for group definition) PAHs in their non-alkylated form (Page et al., 1999; Yunker et al., 2002) constitute the majority of the PAHs produced by the combustion of plant biomass, while the relative amount of high-molecular-weight (HMW) PAHs increases along with higher fire temperatures (McGrath et al., 2003). LMW PAHs tend to be airborne and show high aqueous solubility and higher volatility, whereas HMW PAHs are usually in a solid phase (associated with either soot or char), show lower volatility, and are likely sourced locally (Hoffmann and Wynder, 1977; Junk and Ford, 1980; Karp et al., 2020; Lammel et al., 2009; Lima et al., 2005; Purushothama et al., 1998). Thus, low contributions of HMW PAHs in environmental samples are often considered indicative of either low-temperature fires (e.g. Denis et al., 2012) or a distal source, while high relative amounts generally point toward a more local signal. Finally, perylene is a 5-hexa-ring PAH often detected in aquatic sediments and considered to be mostly of in situ biogenic origin, probably from precursor compounds present in saprophagous and mycorrhizal fungi (e.g. Aizenshtat, 1973; Jiang et al., 2000; Slater et al., 2013; Wang and Huang, 2021), and thus likely linked to higher organic matter content and terrigenous input (Guo and Liao, 2020; Hanke et al., 2019).

3.2 Sterols and stanols as markers of plant sources and animal digestion

Stanols are saturated isomers of sterols (e.g. Patterson, 1971). When the bacterially mediated reduction in sterol double bounds occurs in an open environment (e.g. soil), it leads almost exclusively to the production of 5α stanol isomers. When the reduction in sterols occurs in the animal's digestive tract, their enteric bacterial flora maximises the production of 5β stanols (Hatcher and McGillivary, 1979; Murtaugh and Bunch, 1967). Humans (and partially other omnivores and carnivores) maximise the production of coprostanol (5β-cholestan-3β-ol) through the saturation of animal-derived cholesterol (5-en-cholest-3β-ol). Ruminants such as sheep and cattle, on the other hand, maximise the production of 5β-stigmastanol and 5β-campestanol (Leeming et al., 1994, 1996) from plant-derived sterols like stigmasterol, sitosterol, and campesterol (e.g. Goad, 1977; Goad and Goodwin, 1966; Pancost et al., 2002). Higher and lower ratios of coprostanol and its derived epimer epi-coprostanol (5β-cholestan-3α-ol; (McCalley et al., 1981; Quirk et al., 1980; Wardroper et al., 1978) to 5β-stigmastanol or 5β-campestanol are considered to be a proxy for higher and lower faecal input from human sources relative to ruminant sources and have been widely applied to samples from modern and ancient sewage material and manured soil (e.g. Birk et al., 2012; Bull et al., 2001, 2002; Cordeiro et al., 2008; Evershed et al., 1997; He et al., 2018; Lerch et al., 2022; Simpson et al., 1999; Tyagi et al., 2009).

3.3 n-Alkanes

Plants synthesise n-alkanes and other n-alkyl lipids as part of their waxy coating with a characteristic strong odd-over-even chain length predominance (Eglinton and Hamilton, 1967), which is summarised by their higher carbon preference index (CPI; Bray and Evans, 1961; Marzi et al., 1993). In contrast, lower CPI values are usually indicative of petrogenic, algal, or bacterial sources (Grimalt and Albaigés, 1987; Han and Calvin, 1969). Aquatic sources such as macrophytes and mosses (e.g. Sphagnum) maximise their leaf wax n-alkane production at mid-length homologues (C_21–25), while terrestrial plants (e.g. grasses, sedges, trees, shrubs) are generally skewed toward longer homologues (C_27–31), allowing for use of source-discriminating ratios and indices such as the aquatic plant index (P_aq, Ficken et al., 2000) and average chain length (ACL, Gagosian and Peltzer, 1986).

3.4 Bulk geochemistry proxies

Aquatic and terrestrial catchment productivity, flux of inorganic sediments, and organic matter preservation are the main factors determining the level of total organic carbon content in lacustrine sediments (Meyers and Ishiwatari, 1993). The molar carbon to nitrogen ratio ($\mathrm{C}/\mathrm{N}$) in plant tissue varies between aquatic plants and phytoplankton (<10) and terrestrial plants and bryophytes (>10; Meyers, 1994). Thus, increases in $\mathrm{C}/\mathrm{N}$ are usually interpreted as an increased catchment erosion and input of terrestrial organic matter and/or as a relative decrease in aquatic plant productivity (Fernández-Martínez et al., 2021; Kaushal and Binford, 1999; Meyers, 1997; Meyers and Teranes, 2001; Rieger et al., 1979). Shifts in the abundance of diatom derived biogenic silica (BSi) can trace lake productivity (Colman et al., 1995; Conley, 1988; Conley and Schelske, 2002). The conservation potential of diatom frustules is strongly related to sedimentation rate, with higher rates leading to better preservation. When sedimentation rates are considered relatively constant, shifts in BSi can reflect qualitative changes in spring and summer temperature in high-latitude lakes, such as Iceland (Geirsdóttir et al., 2009a; McKay et al., 2008). The stable isotopic composition of carbon (δ¹³C) can trace shifts in the relative contribution of organic matter sources, with terrestrial plants (but also bryophytes) and associated soils showing more ¹³C-depleted values (ca. −32 ‰ to −25 ‰), aquatic plants exhibiting more ¹³C-enriched values (ca. −20 ‰ to −10 ‰), and freshwater algae and phytoplankton showing a wider isotopic range (Meyers, 1994; Prokopenko et al., 1993; Rundel et al., 1979; Smith and Epstein, 1971; Geirsdóttir et al., 2020 and refs. therein). The physical mixing or stratification of a lake water column can also influence the carbon isotopic signature of aquatic sources (Hernández et al., 2011).

4 Results

4.1 Age model

Based on major oxide composition and stratigraphical information, we identified 13 marker tephra layers of known age (Table 1). Our Bayesian tephra age model shows nearly constant sediment accumulation rates throughout the Holocene (Fig. 2).

Figure 3 Erosional and primary productivity proxies from the pre-G10 ka series tephra interval (pink vertical band) of the 20-SVID-02 core; all concentrations are on grams of dry samples: (a, black) sum of C_19–35 n-alkane concentration aquatic plant index derived from n-alkanes, (b, green) carbon-to-nitrogen ratio, (c, purple) perylene concentration, (d, black and yellow) stable isotopic composition of total carbon, (e, brick red) percentage of total carbon, (f, dark blue) n-alkane carbon preference index (CPI_19–31), (g, red) aquatic plant index derived from n-alkanes, and (h, light blue) biogenic silica.

There is increased uncertainty in age control between the G10 ka tephra series and Kverkfjöll/Hekla 6200 due to fewer marker tephra layers being present. However, the late Holocene, particularly during the historical period of settlement, features numerous tephra layers that result in substantially lower age estimate uncertainty.

To facilitate the interpretation of downcore records, we present and discuss data (1) divided into nine time intervals (I–IX; Table 2) representing the most distinguishable periods of variability with respect to background values and (2) separately for the sections preceding and following the G10 ka tephra series (Óladóttir et al., 2020). We tried analysing some material from the G10 ka series (a very thick tephra unit; Óladóttir et al., 2020), but the yields of organics were too low to be reliable (as expected; its contents are almost 100 % inorganic). Thus, here we only included one sample at the bottom limit of the G10 ka tephra series and one at the upper limit.

Figure 4 Erosional and primary productivity proxies for the 20-SVID-02 core. (a, black): sum of C_19–35 n-alkane concentration aquatic plant index derived from n-alkanes, (b, green) carbon-to-nitrogen ratio, (c, purple) perylene concentration, (d, black and yellow), stable isotopic composition of carbon, (e, brick red) percentage of total carbon, (f, cyan) Northern Hemisphere summer insolation at 65° N (Berger and Loutre, 1999), (g, red) aquatic plant index derived from n-alkanes, (h, light blue) biogenic silica, and (i, dark blue) carbon preference index (CPI_19–31). The dotted red line marks the conventional age (870 CE) of the settlement of Iceland (Landnám). The vertical dashed black lines mark the subdivision of the 20-SVID-02 record into nine intervals (Table 1).

Table 2 Age intervals (approximate) and descriptions of the nine subdivisions of the 20-SVID-02 record used in this study.

4.2 Bulk geochemistry

The $\mathrm{C}/\mathrm{N}$ ratio (Figs. 3 and 4b) ranged from ∼1 to ∼9.5, showing the lowest values at the beginning of the record (I). At ca. 10.63 ka BP, $\mathrm{C}/\mathrm{N}$ increased sharply and reached mid-range values (5–6), remaining relatively stable throughout most of the remaining Holocene (IV–VI), except for two drops, after the G10 ka series tephra and between 8.8 and 7.75 ka BP (V). In the last 2 kyr, $\mathrm{C}/\mathrm{N}$ values steadily increased, leading to the highest values in the most recent portion of the record. The two periods with decreased $\mathrm{C}/\mathrm{N}$ values, as well as the initial increase (I to II), generally paralleled the behaviour of the total carbon (Figs. 3 and 4e), biogenic silica (Figs. 3 and 4h), and δ¹³C (Figs. 3 and 4d) records. TC increased steadily throughout the Holocene from ∼0 % while BSi rapidly increased at the beginning of the record (30 to 90) to then stabilise at 110–120 for more than 5 kyr (V–VI). Both records peaked at ca. 3.5 ka BP (2.8 %, TC; ∼160, BSi max) and temporarily dropped between ca. 3 and 2.25 ka BP. Subsequently, TC increased to its maximum value (3.6 %, modern) while BSi decreased in a stepwise manner, reaching its lowest value of the last 10 kyr (∼50) at ca. 0.2 ka BP. The δ¹³C record showed the most ¹³C-enriched values (−20 ‰) in the oldest interval (I); it then decreased (to −26 ‰, modern) steadily throughout the Holocene, except for two major drops to its most depleted values (ca. −27 ‰) during periods II and V.

4.3 n-Alkanes

We detected n-alkane homologues from C₁₉ to C_33–35 (Fig. A3) in most samples, with a total sum that ranged from 0.3 to 50 µg g⁻¹ (700–4000 µg g⁻¹ TC; Figs. 3 and 4a). The 10.8 to 4.2 ka BP interval showed relatively low and stable values (∼5 µg g⁻¹); concentrations roughly doubled from 4 to 1.5 ka and then again after 0.5 ka BP, reaching its maximum value at the end of the record. The CPI showed a stable odd-over-even predominance (3 to 6.5) through the whole record, except for low values (1 to 3) seen in the interval preceding the G10 ka series tephra (Figs. 3f and 4i). The most abundant homologues were C_23–27 (45 %) in the 10.8 to 3 ka BP interval and C_29–31 (40 %) in the last 3 kyr. This regime change was highlighted by a shift in P_aq from relatively high values (up to 0.8; avg. 0.6) through the early to mid Holocene to lower values (down to 0.2; avg. 0.3) after 3 ka BP (Fig. 4g).

4.4 Faecal sterols and stanols

We detected the three main faecal stanols: the plant-derived 5β-stigmastanol was ∼5–10 times more abundant (up to 3–500 ng g⁻¹; Figs. 5 and 6a) than coprostanol + epi-coprostanol (10–30 ng g⁻¹; Figs. 5 and 6b). The oldest interval (I) showed the lowest concentrations for all three stanols, while the following interval (II) displayed high (highest for coprostanol and epi-coprostanol) concentrations. All three stanols showed low and stable concentrations between ∼9.5 and 7.5 ka BP, gradually increased from ∼7.5 ka BP before reaching a relative maximum around 1.5 ka BP (VII), dropping again during interval VIII, and peaking during the last 200 to 300 years (IX). Parent sterols (β-sitosterol, β-stigmasterol, except cholesterol) and the α-stanol isomers, follow patterns similar to the β-stanols throughout the Holocene, but with 10 to 100 times higher concentrations (Fig. A2).

Figure 5 Sterols and stanols from the pre-G10 ka series tephra interval (vertical pink band) of the 20-SVID-02 core: (a, red) 5β-stigmastanol and (b, dark blue) sum of coprostanol and epi-coprostanol concentration on grams of dry sample (full line) and on grams of TC (dotted line), (c, black) 5α-stigmastanol-to-β-sitosterol ratio, (d, yellow) β-sitosterol concentration, and (e, light blue) 5α-stigmastanol concentration.

Figure 6 Sterols and stanols of the 20-SVID-02 core: (a, red) 5β-stigmastanol and (b, dark blue) sum of coprostanol and epi-coprostanol concentration on grams of dry sample (full line) and on grams of TC (dotted line), (c, black) 5α-stigmastanol-to-β-sitosterol ratio, (d, yellow) β-sitosterol concentration, and (e, light blue) 5α-stigmastanol concentration. The dotted red line marks the conventional age (870 CE) of the settlement of Iceland (Landnám). The vertical dashed black lines mark an arbitrary subdivision of the 20-SVID-02 record into nine intervals (Table 1).

4.5 Polycyclic aromatic hydrocarbons (PAHs)

PAHs were present in all samples and generally in higher concentrations in more recent samples compared to older ones (Fig. A1). Perylene (Figs. 3 and 4c), which accounts for more than 97 % of detected PAHs, maintains low concentrations (0–20 ng g⁻¹) from 10.5 to 4–3 ka BP. A first increase to 0.5–1.5 µg g⁻¹ occurred between 3.5 and 2.8 ka BP, followed by a decrease loosely coeval to the decrease in TC and BSi (VII). After ca. 1.5 ka BP, perylene increases to a maximum value (∼6 µg g⁻¹, ca. 0.15 ka BP), which generally matches the pattern of TC. The second most abundant compound was phenanthrene (0.01–31.1 ng g⁻¹), followed by pyrene (0.1–11.5 ng g⁻¹) and fluoranthene (0.01–17.6 ng g⁻¹). The least abundant PAHs were naphthalene, acenaphthylene, and acenaphthene. However, since the detection of these three low-molecular-weight compounds could have been influenced by evaporation losses during sample preparation, their reported concentrations are likely to be underestimated. Given the overwhelming dominance of perylene and its likely biogenic rather than pyrogenic origin, we removed it from total PAHs abundance calculations to provide a record with features that were not apparent in perylene's trend. In terms of total pyroPAHs abundance, we observe five distinguishable intervals (Fig. 7a). First, a 3-year-long interval (∼80 years average temporal resolution) starting at ca. 10.5 ka BP displays relatively stable low concentrations (<5 ng g⁻¹, SD σ=2.3). The only exception is a point taken within the G10 ka series tephra (up to 10 ng g⁻¹) at ca. 10 ka BP. Second, from ∼7.5 to 2.9 ka BP (∼320 years average temporal resolution), the total PAHs concentration increases from ∼5 to 10 ng g⁻¹, albeit with enhanced variability (σ=3.7). Third, led by the increase in low-molecular-weight PAHs (LMW, Fig. 7b), during the 2.9 to 0.7 ka BP interval, values fluctuate (σ=7.5) between ∼10 and 20 ng g⁻¹, with two major peaks reaching 45 and 35 ng g⁻¹ at ∼2.7 to 2.3 ka BP and 1.5 to 1.3 ka BP, respectively. Fourth, between 0.7 and 0.3 ka BP there is a relatively brief but clear drop (∼10 ng g⁻¹, σ=1.5) led by both LMW and medium-molecular-weight (MMW) PAHs. In the last 250 years we observe a sharp, 10-fold increase in PAHs concentration leading to the highest recorded values (∼200 ng g⁻¹). When normalised for TC (Fig. A1), the absolute values increase 10- to 100-fold, but the patterns do not substantially change.

Figure 7 PAHs record of the 20-SVID-02 core: (a, black) sum of pyrogenic PAHs concentrations; (b, fuchsia) low-, (orange) medium-, and (blue) high-molecular-weight pyrogenic PAHs concentration. The dotted red line marks the conventional age (870 CE) of the settlement of Iceland (Landnám). The vertical pink band marks the G10 ka series tephra. The vertical dashed black lines mark an arbitrary subdivision of 20-SVID-02 record into nine intervals (Table 1).

4.6 HYSPLIT

We calculated a total of 11 392 air trajectories, which we split by year, season, dry (i.e. not associated with precipitation), or precipitation bearing trajectories (Fig. A4). Most trajectories, even for 2-week intervals, show air parcels originating mostly from Iceland and surrounding areas of the North Atlantic, regardless of season or NAO configuration, while the contribution from nearby terrestrial regions (potential PAHs sources) such as Greenland, the British Isles, or Scandinavia is negligible. As the marine environment is not conducive to combustion nor redeposition of particulates, this implies a dominantly Icelandic signal for PAH production.

Since wildfires are concentrated in the relatively dry, snow-free summer season (Mccarty et al., 2021), we focus particularly on the JJA air trajectory data (Fig. 8). These results show that (1) 95 % of back-trajectories originate from Iceland and its nearby waters, (2) 90 % of back-trajectories originate within 100–150 km radius from SVID, (3) dry trajectories more likely originate inland relative to precipitation-carrying trajectories, (4) 95 % of trajectories from a NAO− year tend to be confined to northern Iceland while during a NAO+ year trajectories more commonly originate from inland, and (5) these patterns are consistent even when scaled from 3 d to 2-week intervals (Fig. A4).

Figure 8 HYSPLIT back-trajectories of air parcels for 2010 (NAO−) and 2013 (NAO+) for annual and summer (JJA) data. Trajectories are calculated on a 2-week (336 h) interval at a 6 h frequency; “precip” indicates trajectories that produced precipitation within 6 h from the SVID endpoint, while “dry” trajectories did not. Contour colours indicate the frequency at which air parcels part of a trajectory travel above a certain area.

5 Discussion

5.1 Primary aquatic production vs. erosion/terrestrial input

In Icelandic lacustrine environments, δ¹³C and $\mathrm{C}/\mathrm{N}$ are generally considered proxies for the relative contribution of terrestrial vs. aquatic organic matter and shifts in primary productivity, as total carbon is virtually solely of organic origin (e.g. Geirsdóttir et al., 2009a). Iceland's bedrock is dominantly comprised of basaltic bedrock, including the catchment of SVID (Hjartarson and Sæmundsson, 2014), meaning there is negligible carbonate available. Although some dissolved inorganic carbon (DIC) has been measured in Icelandic rivers, it is greatly outweighed by organic carbon (Kardjilov et al., 2006), whereas the amount of inorganic carbon measured in soils is negligible (Mankasingh and Gísladóttir, 2019). This is important to consider because SVID has no river inflow and water inflow is dominated by runoff from the catchment through soil. In addition, while there is some evidence for the transport of Saharan dust to Iceland within the last decade (Varga et al., 2021), there is currently no evidence of such transport during the Holocene. Additional pools of inorganic carbon from aquatic invertebrates, such as ostracods, are also considered to be negligible (Alkalaj et al., 2019). Ostracods crystallise their shells in very close equilibrium to the carbon isotopes of DIC (Decrouy, 2012), which for Iceland is notably enriched relative to bulk organic matter carbon isotopes (Sveinbjörnsdóttir et al., 2020). If ostracods were a substantial contributor to the total carbon pool, we would expect lake sediment carbon isotopes to deviate significantly from modern terrestrial and aquatic plant carbon isotope values. As Icelandic lake sediment bulk geochemistry is consistent with the fields of modern plants (see Geirsdóttir et al., 2020), inorganic carbon from aquatic invertebrates is not considered a significant contributor to the total carbon pool.

In addition, in SVID, the similarity of the $\mathrm{C}/\mathrm{N}$ record to the perylene curve reinforces its significance as a proxy for terrigenous input.

5.1.1 Postglacial warming: 11–7.5 ka BP

Deglaciation in the NE of Iceland set in between 15 and 13 ka BP and proceeded in a stepwise fashion, with two main glacier re-advances at ca. 2.7 (Younger Dryas) and 10.9 ka BP (Preboreal; Geirsdóttir et al., 2009b; Norðdahl and Pétursson, 2005). Our record captures sediment below the Askja S tephra layer (10.83±0.57 ka BP, Bronk Ramsey et al., 2015), showing ice-free conditions and the start of organic sedimentation by 10.85 ka BP at SVID's location. Except for the oldest sample, high δ¹³C and low TC values indicate a primarily aquatic source of carbon during the oldest interval (I) (Fig. 3). This suggests an absence of substantial terrestrial vegetation, consistent with a postglacial landscape and possibly a cooler climate associated with the Preboreal period.

TC, BSi, and $\mathrm{C}/\mathrm{N}$ increase suddenly at ca. 10.65 ka BP, maintaining higher values for 2–3 centuries (vice versa for δ¹³C), indicating an enhanced terrestrial input likely resulting from a retreating glacier, development of soil and vascular plants, and generally warming conditions (Fig. 3). After the G10 ka series tephra, all proxy values decrease, likely due to the destructive impact of substantial volcanic ash fallout on both terrestrial and aquatic vegetation and related water chemistry alteration (e.g. δ¹³C dropping due to acidification; Kilian et al., 2006) (Fig. 4). Following the volcanic event, all proxies increase at ca. 9.75 ka BP, whereas terrestrial- relative to aquatic-sourced carbon temporarily increase (ca. 8.7–7.5 ka BP). The observed decrease in δ¹³C (and, partially, $\mathrm{C}/\mathrm{N}$) between 8.7 and 7.5 ka BP is identified in other Icelandic lake sediment records between 8.8 and 7.9 ka BP (e.g. Eddudóttir et al., 2018; Geirsdóttir et al., 2013; Harning et al., 2018b; Larsen et al., 2012) and has been attributed to the likely impact of meltwater pulses into the northern North Atlantic due to the retreating Laurentide ice sheet and/or local effusive volcanic eruptions (Geirsdóttir et al., 2013; Larsen et al., 2012).

The total sum of n-alkanes (C_19–35; Fig. 4a), which is heavily controlled by C₂₉ and C₃₁ (Fig. A3), increases throughout the Holocene similar to the pattern described by $\mathrm{C}/\mathrm{N}$ and perylene as a result of an increased terrigenous input. As inferred by high CPI values, most n-alkanes in SVID originate from plants (terrestrial and possibly also aquatic) through the Holocene record (Fig. 4i). The relatively low CPI values in the oldest interval (I) indicate a negligible contribution from plant sources (relative to phytoplankton) to the carbon pool, which is consistent with a still cold, relatively barren, deglacial environment.

The CPI curve shows a similar but opposite pattern to the δ¹³C record until ca. 8 ka BP, reinforcing δ¹³C as a proxy mostly controlled by terrigenous input (n-alkanes from aquatics show lower CPI values than terrestrial plants; e.g. Bray and Evans, 1961; Duan et al., 2014; Eglinton and Hamilton, 1967; Li et al., 2020). The reason for the change in the relationship between CPI and δ¹³C (which become positively correlated after 8 ka BP) is unclear. While it matches the timing of increasing temperatures (Axford et al., 2007; Fig. 10d) and Betula expansion in the region (Karlsdóttir et al., 2014; Fig. 10e), its interpretation is complicated by the fact that CPI can also be influenced by factors such as changes in mean annual precipitation, seasonality, plant community, and algal productivity (Li et al., 2020).

5.1.2 Mid Holocene plateau and trend inversion: 7.5–4.2 ka BP

Overall, all proxies suggest that the interval between ∼7.5 and 4.2 ka BP was characterised by relatively stable climatic conditions, generally warmer (Axford et al., 2007) and wetter (Moossen et al., 2015) than both the preceding and following periods. These conditions likely led to an enhanced primary productivity within the lake, as suggested by high values of BSi, P_aq, CPI, and pollen-inferred vegetation communities (Eddudóttir et al., 2016; Karlsdóttir et al., 2014). In fact, modelled reconstructions of Icelandic vegetation cover throughout the Holocene show the highest values during this interval (Ólafsdóttir et al., 2001). In particular, partly due to the retreat and disappearance of most glaciers, between 50 % and 60 % of Iceland was likely covered in vegetation, of which a quarter was birch forest, throughout the mid Holocene, with two peaks, one at the Holocene Thermal Maximum (8–7 ka BP) and one at ∼3.5 ka BP (Ólafsdóttir et al., 2001). The warm and moist climate, paired with the expansion of vegetation, may have stabilised the local environment, reducing erosion. Such stability is consistent with the low $\mathrm{C}/\mathrm{N}$, perylene, and n-alkane values throughout this interval.

Though this four thousand year period is broadly categorised by stability, a more detailed view reveals important inflection points in the long term trends of many proxies. For example, while some proxies keep increasing (e.g. TC, and summer temperature, Axford et al., 2007), some stop rising and remain relatively flat throughout this interval (e.g. $\mathrm{C}/\mathrm{N}$, Perylene), while others even invert their trends (δ¹³C, P_aq, CPI). Another trend inversion occurred around 5 ka BP with the inception of neo-glaciation in Iceland when glaciers started to expand again (Geirsdóttir et al., 2019). This could be interpreted as a slow inertial response of the local environment to the decreasing NH summer insolation, likely reducing its resilience to short-term events such as volcanic eruptions and NAO shifts, until some kind of threshold was finally reached around 4 ka BP (Geirsdóttir et al., 2013, 2019).

5.1.3 Increased erosion in a cooling climate: 4.2 ka BP

Our SVID proxy datasets generally agree with previous work using bulk geochemistry proxies in Icelandic lakes, which collectively point toward decreasing primary productivity and increasing landscape instability in response to declining Northern Hemisphere summer insolation (e.g. Geirsdóttir et al., 2013, 2019; Harning et al., 2018b, 2020; Larsen et al., 2012). The decreasing trend in δ¹³C is generally anticorrelated with the TC curve, indicating that TC is increasingly controlled by terrestrial input. BSi, TC, and $\mathrm{C}/\mathrm{N}$ values drop or invert their trend after ca. 4 ka BP, consistent with a general decrease in productivity. This is possibly related to a combination of decreasing moisture and/or summer temperatures (Axford et al., 2007) and the effect of the Hekla 4 volcanic event pushing the local environment beyond a threshold (Eddudóttir et al., 2017). Absolute amounts of n-alkanes also increase starting at 4 ka (Fig. 4a), led by an increase in n-C₂₉ (Fig. A3), typical of terrestrial plants. The ratio between the mid- and long-chain homologues (aquatic plant index or P_aq; Ficken et al., 2000) is often interpreted as a proxy for a wetter or dryer environment. However, this is likely an oversimplification as questions remain about the relationship between n-alkane chain length and vegetation source, particularly for aquatic plants, which are often a minor component of the leaf wax pool in Arctic lakes (Dion-Kirschner et al., 2020; Hollister et al., 2022). Nevertheless, when coupled to the concentration data of n-alkanes (Figs. 4a and A3), P_aq can here be more safely interpreted as indicative of lower or higher terrigenous input. The SVID record (Fig. 4h; Sect. 3.3) shows higher P_aq values (mean 0.6) during the 10.8 to 4–3 ka BP interval, indicating an environment with significant aquatic plant production and likely limited erosion or in-wash. The record then switches abruptly to lower values (mean 0.3) at 3 ka BP, highlighting a shift toward greater in-wash of terrestrially derived material.

The massive Hekla 3 eruption (the most severe Hekla eruption of the Holocene; ∼3010 a BP; Larsen, 1977; Larsen and Eiríksson, 2008) was likely the cause, or at least the trigger, of this abrupt shift at 3 ka BP in most SVID proxy records, particularly the ones related to primary production and erosion (Larsen et al., 2011). In fact, the volcanic fallout likely killed terrestrial plants by burning, root suffocation, and reduced photosynthesis (e.g. Ifkirne et al., 2022; Mack, 1981; De Schutter et al., 2015) and had likely similar effects on aquatic flora as well, also through increased turbidity and acidity of lake waters (e.g. Ayris and Delmelle, 2012). The subsequent reduced coverage of terrestrial plants likely exposed the soil to increased erosion, resulting in more terrestrial in-wash (as reflected by a sudden perylene peak), skewed toward the inorganic components of soil (as reflected by a sudden drop in TC and n-alkanes lasting roughly a century). The increased terrestrial in-wash would have further reduced primary productivity within the lake, as suggested by the drop in BSi. At the same time, the short $\mathrm{C}/\mathrm{N}$ peak and the major drop in P_aq seem to indicate that the productivity and contribution of terrestrial plants remained higher than aquatic sources (Larsen et al., 2011). The post-3 ka trend is temporarily interrupted by what appears to be a partial rebound in primary productivity (BSi increases too) and diminished in-wash of terrestrial material until ca. 1.5 ka BP, whereafter it continues to decline.

Figure 9 Holocene sub-Arctic records of faecal stanols in North Atlantic islands: (a) 5β-stigmastanol (dark red) and sum of coprostanol and epi-coprostanol (dark blue) from core 20-SVID-02, NE Iceland (this study); (b) 5β-stigmastanol (pink) and coprostanol (light blue) from cores LILA09-LILC09 from Lilandsvatnet lake, Lofoten Islands (D'Anjou et al., 2012); and (c) 5β-stigmastanol (orange) and coprostanol + epi-coprostanol (green) from core EI-D-01-15 from Eiðisvatn lake, Faroe Islands (Curtin et al., 2021).

5.2 Is there geochemical evidence for human settlement in the SVID catchment area?

In paleoclimate studies, relative shifts to above natural background levels of β-stanols have been used as a proxy for human settlement, marking the appearance of humans and domesticated animals in specific areas of the world (Shillito et al., 2020; Sistiaga et al., 2014), often in lake catchments (Battistel et al., 2016; Callegaro et al., 2018; Raposeiro et al., 2021; Sear et al., 2020; Vachula et al., 2019, 2020). This method has detected the arrival of Viking settlers in other Nordic regions, such as in the Lofoten islands in northern Norway (D'Anjou et al., 2012, Fig. 9b) and in the Faroe Islands (Curtin et al., 2021, Fig. 9c). However, SVID sterol and stanol records show no evident human signals at or around the time of colonisation (i.e. 9th century CE, ca. 1.1 ka BP). While a relative maximum of sterol and stanol concentrations found at ca. 1.4 ka BP resembles the timing of an earlier-than-colonisation stanol signal found in the Faroes (Curtin et al., 2021), as well as an analogous signal in the Lofoten Islands (D'Anjou et al., 2012), we cannot confidently interpret this peak as indication of human presence as its amplitude is comparable in magnitude to the inherent variability in the record (i.e. low signal-to-noise ratio).

The lack of a clear anthropogenic faecal biomarker signal could be explained by either (1) a scarce or null incidence of human activities in the catchment (unlikely, given the archaeological evidence in nearby areas; Gísladóttir et al., 2012; Lebrun et al., 2023) and/or by (2) dilution of the signal due to the relatively large size of the lake paired with a small catchment. The sterol and stanol records show a general increase throughout the Holocene (Fig. A2) in a pattern that matches the $\mathrm{C}/\mathrm{N}$, n-alkanes, and perylene trends, suggesting that the primary driver of SVID's sterol signal is likely landscape instability and soil erosion rather than human or ruminant presence. Furthermore, ratios of sterols to their derived 5β-5α stanols can trace redox conditions in various environments (e.g. Andersson and Meyers, 2012; Canuel and Martens, 1993; Jaffé et al., 1996; Routh et al., 2014) and thus, potentially, human presence in a lake catchment, as anthropogenic activities tend to mobilise more soil and increase in-wash of organic material, fostering reducing conditions (Argiriadis et al., 2018). In SVID, the stanol values (5β-5α) are consistently lower than their respective sterol precursors, suggesting a generally oxidising environment throughout the Holocene (Fig. 6). The only exception to this trend is in the earliest part of the record (ca. 10.8–10.6 ka BP; Fig. 5), indicative of a more reducing environment, though this is not linked to an increased organic input (low TC values) but more likely to lake stratification with deglacial water sinking at the bottom of the lake (Sugiyama et al., 2021).

5.3 Holocene fire frequency

Pyrogenic PAHs are considered a reliable proxy for fire frequency on a local scale both within and around a catchment (Denis et al., 2012). Although other factors can influence the PAH signal in sedimentary archives (e.g. accumulation rates, degradation; Stogiannidis et al., 2015), we interpret SVID pyroPAHs data as a record of NE Iceland fire history through the Holocene.

The trend in pyroPAHs does not match the erosional signal described by bulk geochemical proxies and n-alkanes, suggesting that soil erosion is not a mechanism for Holocene pyroPAHs variability. We exclude chemical degradation as a source of the signal, as PAHs are relatively stable molecules on long timescales (e.g. Johnsen et al., 2005). In fact, the ratio between low-molecular-weight (more prone to chemical degradation and leeching) and high-molecular-weight PAHs (as defined in Fig. A5) remains above 1 in most samples. Similarly, the pyrene $/$ coronene ratio shows high and stable values throughout the record, indicating good preservation, with no significant degradation or preferential removal of less recalcitrant PAHs such as pyrene (Fig. A5 and refs therein).

The pyroPAH record presents two peaks at ca. 2.8 and 1.5 ka BP, both predating acknowledged human settlement. PyroPAH values subsequently drop (VIII) and increase again in the last 2 centuries, reaching maximum values in the present (Fig. 7a). Analysing PAH data subdivided in molecular weight classes (see Sect. 3.1; Table A1; Fig. 7b) can help explain these two features as well as the general trend. HMW pyroPAHs show low and stable relative contributions (12±7 % of total pyroPAHs) through the whole record and rise to 70 % in the last 200 years. This trend, paired with increased pyroPAH concentrations, is consistent with the burning of coal and oil, the use of internal combustion engines, and increased human presence (Abas and Mohamad, 2011; Kozak et al., 2017). MMW pyroPAHs show relatively stable concentrations through most of the record (until ca. 0.2 ka BP). Their relative abundance ($\sim \mathrm{68}\pm \mathrm{35}$ %) slowly decreases through the Holocene, proportional to the increase in LMW pyroPAHs ($\sim \mathrm{20}\pm \mathrm{16}$ %). The latter, which are predominantly present in the gaseous phase (Karp et al., 2020), peak at 3–2.8 and 1.5 ka BP, and substantially control the shape of the pyroPAH record (Fig. 7a and b). Together, the (1) low and stable values of HMW pyroPAHs, the (2) stable MMW values, and the (3) increasing and peaking values of LMW pyroPAHs are consistent with a general increase in the frequency of low-temperature fires (e.g. peat fires or crawling fires) at a regional level. While the Hekla 3 event (3.01 ka BP, the largest rhyolitic eruption during the Holocene; Larsen, 1977; Larsen and Eiríksson, 2008) occurs just before the first pyroPAH peak, it is unlikely to be its unique or even main cause. However, in an environment such as Iceland, it is legitimate to ask if the frequent eruptions might have influenced the PAHs' natural background and if this is detectable in our records. As discussed below, we are confident that volcanic eruptions have had no significant direct impact on the amount and distribution of PAHs in the SVID archive.

Volcanic activity does produce PAHs (both in the gaseous and in the particulate phases) but generally their long term contribution is considered negligible, particularly in the modern world where the main PAHs source is the burning of fossil fuels (Guiñez et al., 2020; Kozak et al., 2017).

Overall, due to the high temperatures involved, volcanic eruptions tend to produce medium-high molecular weight PAHs in gases and particulates (Guiñez et al., 2020; Ilyinskaya et al., 2017). Volcanic layers can contain pyrogenic (unsubstituted but also alkylated) PAHs with a molecular weight distribution resembling modern fossil fuel combustion, dominated by unsubstituted forms (Murchison and Raymond, 1989), as well as traces of nitro- and oxy-PAHs (Guiñez et al., 2020). We do not see this distribution in any of our samples, not even when they include parts of tephra layers. This suggests that (1) either our sampling method does not capture volcanic layers (maybe due to its resolution) or that (2) there is no such signal in the SVID archive.

Regarding the first hypothesis (sampling method not capturing volcanic PAHs), Kozak et al. (2017) analysed the impact of 2010 and 2011 Icelandic eruptions in Svalbard (Arctic Norway), finding high abundances of 4–5 ring PAHs in volcanic mud, great variations in the contribution of different eruptions to the total PAHs detected in sampled surface water, significant increases in the total abundance of PAHs during eruption years, and also that this increase in PAHs abundance and shift in PAH distributions do not seem to last beyond the eruption years. Considering the temporal resolution of our record (maximum of 10–50 years), it is unlikely that any eruption would have impacted it significantly: from a geological point of view, eruptions tend to be short-lived, and unless they relate to a massive event sustained over a long period of time, they are unlikely to have a strong impact on a sample that represents 10–50 years of sedimentation. No discernible correlation arose between the PAHs curve and the detected tephra, the only exceptions being the Hekla 3 and, to a lesser extent, the Hekla 4 tephra layers. However, these major volcanic events, besides marking the beginning of major shifts in most proxies, correlate to the initial phase of major shifts in PAHs of low (and not high nor medium) molecular weight . This suggests that if a connection between these two eruptions and PAH shifts exists, it must be indirect and, more likely, that the two events acted as a general destabilising factor in an environment already subjected to increasing cooling and erosion.

The only possible example of an increased PAH concentration due to volcanic sources could come from the sample obtained from the lower limit of the G10 ka series tephra (Fig. 7), which has a different composition (mostly inorganic) compared to the rest of the organic-rich samples in our record. In fact, HMW PAHs seem to spike here, even if they still exhibit an overall lower concentration than LMW PAHs; this could also be due to its massive nature (Óladóttir et al., 2020).

Regarding the second hypothesis (missing signal of volcanic PAHs), it is possible that no detectable volcanic PAHs were present in the SVID archive due to (1) its geographical location, relatively far from volcanic sources and formations (Hjartarson and Sæmundsson, 2014), and (2) the nature of Icelandic volcanic eruptions, which are characterised by relatively low-viscosity basaltic lava rather than highly explosive pyroclastic flows (Thordarson and Höskuldsson, 2008), thus reducing the chance of ash production and deposition, particularly in distal locations such as SVID.

The effects of tephra fallout on vegetation and related PAH deposition seem to be quite short lived, with fires events likely coeval to the eruption and vegetation recovering within a few decades (Eddudóttir et al., 2017; Pickarski et al., 2023), while SVID pyroPAH peaks are clearly led by increases in LMW PAHs on a longer timescale.

Notably, this shift in fire regime at ca. 3 ka BP in Iceland falls within a wider pattern of increasing fire frequency emerging from the analysis of several Holocene fire records throughout Europe (Marlon et al., 2013). This is linked to either an increase in cultivated land (fire was used to clear land for agriculture) and/or, particularly in Europe, to increasing aridity (Marlon et al., 2013). We hypothesise that the latter is the most probable explanation for the shifts in NE Icelandic fire regimes as discussed in Sect. 5.4.

5.4 Regional drivers of precipitation and their role in fire frequency

Fuel moisture content and, more generally, environmental moisture are the main variables controlling flammability in vegetational communities typical of temperate and sub-Arctic regions (Marino et al., 2010; Plucinski et al., 2010; Santana and Marrs, 2014). The North Atlantic Oscillation (NAO; Hurrell et al., 2003) modulates the intensity of the westerly storm track and thus the amount and source of precipitation in Iceland. Its positive mode (NAO+) brings intervals of higher precipitation resulting in a wetter (and often warmer) climate than NAO− intervals, which are characterised by weaker westerlies, stronger northerly winds, and drier (and often colder) conditions (Hurrell, 1995; Trouet et al., 2009). Major changes in precipitation regimes usually lead to changes in the hydrogen stable isotopic value of environmental water (δD; Dansgaard, 1964) which translate into shifts in the δD of plant waxes (e.g. n-alkanes; Sachse et al., 2012). This relationship has been calibrated in various environments, including the Arctic (e.g. Berke et al., 2019; Bush et al., 2017; McFarlin et al., 2019; Thomas et al., 2016) and applied for paleo-precipitation reconstructions (e.g. Ardenghi et al., 2019; Niedermeyer et al., 2016; Tierney et al., 2017; Wilkie et al., 2013). A C₂₉ n-alkane δD record from a fjord core in NW Iceland (Fig. 10c; Moossen et al., 2015) describes a relatively stable NAO+ configuration (wetter and more D-depleted) throughout the Holocene, and two major shifts toward NAO− (dryer and less D-depleted) conditions at ca. 3–2.5 and 1.5–1.0 ka BP, matching the timing of first SVID pyroPAH peak and at least partially overlapping with the second one. A similar correlation of fire frequency shifts to NAO− configurations has recently been suggested for other Arctic sites, particularly in Svalbard (Chen et al., 2023).

Figure 10 Regional comparison of climatic, fire, and vegetation records: (a, purple) levoglucosan concentrations from the RECAP ice core, south-eastern Greenland (Segato et al., 2021); (b, black) sum of pyrogenic PAH concentrations, where the red line indicates a three-point running average; (c, blue) stable isotopic composition of hydrogen of sedimentary n-alkanes from marine core MD99-2266 off the coast of NW Iceland (Moossen et al., 2015), where the cyan line indicates a three-point running average; (d, red) chironomid derived temperatures from core 04-SVID-03 (Axford et al., 2007); and (e) pollen percentages of Betula (green), Poaceae (yellow), and Ericales (blue) in a peat section from the Ytra-Áland site, NE Iceland (Karlsdóttir et al., 2014).

Biomass typology (i.e. the kind of vegetation on site) also influences fuel flammability (Chandler et al., 1983; Fernandes and Cruz, 2012; Santana et al., 2011; Scarff and Westoby, 2006). In Iceland, many plant taxa appeared shortly after deglaciation (e.g. Alsos et al., 2021; Harning et al., 2023). Increasing summer temperatures led to the expansion of thermophilic woody plant taxa (e.g. birch) during and after the Holocene Thermal Maximum (e.g. Eddudóttir et al., 2016; Geirsdóttir et al., 2022; Karlsdóttir et al., 2014). Since ∼6 ka BP, the birch woodland in the NE region has evolved into more open heathland and peatland, until the birch population decreased around 3 ka BP (Roy et al., 2018), along with a general temperature decrease (Axford et al., 2007). In this context, at Ytra-Áland (Fig. 1b), two major drops in Betula pollen coeval to two increases in Ericales (heather's order) pollen closely follow the NAO− shifts at 3 and 1.5 ka BP (Fig. 10e; Karlsdóttir et al., 2014). Heathlands, especially in low-moisture conditions, are associated with higher flammability, particularly high sustainability (i.e. how well the combustion proceeds). Thus, heathlands are prone to longer, more stable fires, and with an increased potential for igniting higher canopy elements and underlying peat layers (Plucinski et al., 2010; Rein et al., 2008; Santana and Marrs, 2014). More frequent, more stable, and slow-crawling fires involving dense bushes and peat would increase the amount of pyroPAHs produced and deposited in the region while skewing their distribution toward LMW components (George et al., 2016; Iinuma et al., 2007; Siao et al., 2007), as observed in the SVID pyroPAH signal (Fig. 7). Lastly, the observed shifts in pyroPAHs are unlikely to be the result of a change in their source area. Our back-trajectory analysis reveals that more air parcels originate over Iceland or in the surrounding North Atlantic (Fig. 8). The analysis further reveals that trajectories with terrestrial origins are more likely in NAO+ than NAO− regimes. As such, terrestrial trajectories would be the ones responsible for bringing combustion products to SVID, and we might therefore expect a stronger pyroPAH signal during NAO+ intervals. However, our results display the opposite trend, with the initial increase in pyroPAH abundances at ∼3 ka occurring during an NAO− mode and their subsequent drop occurring during a strong NAO+ interval. Thus, the late Holocene increase in pyroPAHs in SVID is likely to record a substantial increase in local fires (driven by vegetation change and NAO-modulated aridity) rather than a shift in the compounds' sources.

Overall, the SVID pyroPAH signal describes a major shift from a relatively stable (i.e. low fire frequency) early and mid Holocene environment to a dryer late Holocene environment at ca. 3 ka BP that is naturally more prone to long-term persistence of low-temperature wildfires. This is likely the result of the combination of (1) recorded cooling and related shifts in vegetational communities and (2) NAO− shifts and associated dryer conditions.

5.5 Evidence for a human influence on fire frequency?

Unlike other proxies, pyroPAHs return to background levels after reaching high values at 1.5 ka BP, remain low through the Medieval Warm Period (ca. 900–1200 CE) and most of the Little Ice Age (1300–1900 CE), before peaking again in the last 150–200 years. A similar drop in fire markers (anhydrosugars) is also observed in eastern Greenland (Segato et al., 2021; Fig. 10a). Low pyroPAH levels during an interval of known human presence in Iceland suggests that human activities might have curbed regional fire frequency, thus modulating the natural signal, which would have otherwise remained relatively high due to the persistence of colder conditions and more flammable plant communities (regardless of NAO shifts). In fact, while increased pressure from grazing lowered environmental resilience to soil erosion (e.g. Bates et al., 2021; Eddudóttir et al., 2016; McGovern et al., 2007), it likely also decreased fuel flammability (which is dependent on the amount of dead biomass; Davies and Legg, 2011; Santana and Marrs, 2014) in the predominant heathland environment (Lake et al., 2001). Additionally, the creation of farmland and pastureland at the expense of areas with woody vegetation and heathland likely reduced the extent of the biomes naturally prone to fires. This reduction in local wood availability is reflected in changes in foraging habits, as settlers shifted to relying more heavily on more abundant fuel sources such as peat and turf, as well as other marine-derived substances (e.g. seal oil, seaweed; Bold, 2012), while driftwood and imported wood (from Europe or North America, often with ad hoc expeditions) satisfied most of the need for timber (e.g. Bold, 2012; Edvardsson, 2010; Mooney, 2016; Pinta, 2021; Sveinbjarnardóttir et al., 2007). A general mechanism for fire suppression due to the expansion of cultivated land has already been proposed for global data (Marlon et al., 2013) but assumed to be likely asynchronous in different regions and strongly influenced by local climatic, environmental, and social conditions.

We speculate that the drop in fire markers in Iceland from the reduction in wildfire risk due to husbandry and farming exceeded the production of fire markers due to human necessities (e.g. warming), resulting in an overall suppressed fire signal. This would also be consistent with the low population density, which started to increase only in the 1800s CE (Iceland Statistical Service, 2023; Jónsson and Magnússon, 1997), matching the coldest interval of the Little Ice Age and the sharp rise in HMW pyroPAHs. From this perspective, the pressure of human activities would have fostered erosion through decreased environmental resilience while at the same time suppressing natural fire frequency.

6 Conclusions

Our multiproxy analysis of Holocene sediments from Stóra Viðarvatn provides new insight into the coupled vegetation, fire, erosion, and climate regimes of NE Iceland.

• Bulk geochemistry proxies show that the general climatic evolution of NE Iceland is primarily driven by summer insolation: an initial deglacial warming followed by a relatively warm and stable climate until ca. 4–3 ka BP, after which declining summer temperatures result in accelerating catchment erosion.

• Faecal biomarkers, traditionally linked to human activities, do not show an elevated signal at or around colonisation (9th century CE). Instead, faecal biomarkers roughly trace the erosional signal described by bulk geochemical proxies. This may result from a combination of (1) low local anthropogenic pressure (although sparse settlements existed a few km from the study area) and (2) signal dilution due to the large lake size and its relatively small watershed. Therefore, we urge caution when interpreting faecal biomarkers as unequivocal proxies for human presence, particularly when highly sensitive analytical tools like the one used in this study are involved.

• PyroPAHs carry a regional (mostly confined to northern and northeastern Iceland) and predominantly natural signal (i.e. controlled by parameters such as precipitation and moisture availability, vegetation typology and flammability). After generally low fire frequency throughout most of the Holocene, we observe major regime changes at 3 ka and 1.5 ka BP, before known human colonisation in Iceland. During this interval, the distribution of pyroPAHs point toward a regional increase in low-temperature-fire frequency. This can be linked to a change in vegetation typology driven by the cooling of the last 4 to 3 kyr, coupled to major shifts in atmospheric circulation (i.e. NAO regimes) that led to increased aridity and thus flammability. Finally, low levels of pyroPAHs characterise the time following known human colonisation, before rising again (but with a molecular composition more distinctive of fossil fuels) in the last ∼200 years. This suggests that human activities, particularly husbandry and farming, may have suppressed fire frequency by reducing the range and flammability of environments more prone to fire, effectively modulating the natural signal while decreasing the resilience of the local environment to soil erosion.

Appendix A

The $\mathrm{L}/\mathrm{H}$ index (black) is a ratio between low- and high-molecular-weight unsubstituted PAHs, defined as follows (Magi et al., 2002; Stogiannidis et al., 2015, and references therein).

$$\begin{array}{rl}\mathrm{L}/\mathrm{H}& =(\mathrm{Phenanthrene}+\mathrm{Anthracene}+\mathrm{Fluoranthene}\\ & +\mathrm{Pyrene})/(\mathrm{Benzo}\left[\mathrm{a}\right]\mathrm{anthracene}+\mathrm{Chrysene}\\ & +\mathrm{Benzo}\left[\mathrm{k}\right]\mathrm{fluoranthene}+\mathrm{Benzo}\left[\mathrm{a}\right]\mathrm{pyrene}\\ & +\mathrm{Indeno}[\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{c},\mathrm{d}]\mathrm{pyrene}+\mathrm{Dibenzo}[\mathrm{a},\mathrm{h}]\mathrm{anthracene}\\ & +\mathrm{Benzo}[\mathrm{g},\mathrm{h},\mathrm{i}]\mathrm{perylene})\end{array}$$

The pyrene–coronene index (orange) defined as pyrene/(pyrene + coronene) is based on the assumption of an higher preservation potential of the HMW coronene over the lighter, more soluble pyrene (Denis, 2016; Denis et al., 2017; May et al., 1978); higher, more stable values point toward good preservation for both HMW and LMW PAHs.

Table A1 Polycyclic aromatic hydrocarbons (PAHs) analysed in this study. Pyrogenic PAHs are grouped into low, medium, and high molecular weight. Elution order and SRM transitions are reported for each compound.

^a Compound(s) difficult to quantify correctly and thus excluded from final sums. ^b Non-pyrogenic PAH.

Figure A1 Concentration curves of all 20-SVID-02 pyrogenic PAHs recovered in this study. All concentrations are expressed as nanograms per gram of dry sample, except for the last curve (dotted red) which is in nanograms per gram of TC. Note that several vertical axes have been adjusted to minimise the rise in the last 2–3 centuries. Compounds are listed in chromatographic order and grouped by molecular weight through colour shading (LMW in magenta, MMW in orange, HMW in blue).

Figure A2 Concentration curves of all 20-SVID-02 faecal sterol and stanols of interest recovered in this study. All concentrations are expressed as nanograms per gram of dry sample. Note that several vertical axes have been adjusted to minimise the rise over the last 2–3 centuries. Compounds are grouped and colour shaded by structure (cholesterol and stanol derivatives in blue purple, sito-stigmasterols and sito-stigmastanols in red orange, campesterol and campestanol in green).

Figure A3 Concentration curves of all 20-SVID odd-numbered n-alkane homologues from n-C₁₉ to n-C₃₅ recovered in this study. All concentrations are expressed as nanograms per gram of dry sample.

Figure A4 HYSPLIT back-trajectories of air parcels for December 2009 to November 2010 (NAO−) and May 2013 to April 2014 (NAO+, except for October 2014). Trajectories are calculated on 3 d (72 h, a–c) and 2-week (336 h, b–d) intervals at a 6 h frequency. “Precipitation” (top) indicates trajectories that produced precipitation within 6 h from the SVID endpoint, while “dry” (middle) indicates the opposite. The bottom plots are the sum of “precipitation” and “dry” trajectories. Contour colours indicate the frequency with which air parcel trajectories pass over a specific point on the Earth's surface.

Figure A5 PAH indices used as an indication of PAH preservation.

Table A2 Faecal sterols and stanols analysed in this study. Elution order and SRM transitions are reported for each compound.