The onset of neoglaciation in Iceland and the 4.2 ka event

Strong similarities in Holocene climate reconstructions derived from multiple proxies (BSi, TOC – total organic carbon, δ¹³C, C∕N, MS – magnetic susceptibility, δ¹⁵N) preserved in sediments from both glacial and non-glacial lakes across Iceland indicate a relatively warm early to mid Holocene from 10 to 6 ka, overprinted with cold excursions presumably related to meltwater impact on North Atlantic circulation until 7.9 ka. Sediment in lakes from glacial catchments indicates their catchments were ice-free during this interval. Statistical treatment of the high-resolution multi-proxy paleoclimate lake records shows that despite great variability in catchment characteristics, the sediment records document more or less synchronous abrupt, cold departures as opposed to the smoothly decreasing trend in Northern Hemisphere summer insolation. Although all lake records document a decline in summer temperature through the Holocene consistent with the regular decline in summer insolation, the onset of significant summer cooling occurs ∼5 ka at high-elevation interior sites but is variably later at sites closer to the coast, suggesting that proximity to the sea may modulate the impact from decreasing summer insolation. The timing of glacier inception during the mid Holocene is determined by the descent of the equilibrium line altitude (ELA), which is dominated by the evolution of summer temperature as summer insolation declined as well as changes in sea surface temperature for coastal glacial systems. The glacial response to the ELA decline is also highly dependent on the local topography. The initial ∼5 ka nucleation of Langjökull in the highlands of Iceland defines the onset of neoglaciation in Iceland. Subsequently, a stepwise expansion of both Langjökull and northeast Vatnajökull occurred between 4.5 and 4.0 ka, with a second abrupt expansion ∼3 ka. Due to its coastal setting and lower topographic threshold, the initial appearance of Drangajökull in the NW of Iceland was delayed until ∼2.3 ka. All lake records reflect abrupt summer temperature and catchment disturbance at ∼4.5 ka, statistically indistinguishable from the global 4.2 ka event, and a second widespread abrupt disturbance at 3.0 ka, similar to the stepwise expansion of Langjökull and northeast Vatnajökull. Both are intervals characterized by large explosive volcanism and tephra distribution in Iceland resulting in intensified local soil erosion. The most widespread increase in glacier advance, landscape instability, and soil erosion occurred shortly after 2 ka, likely due to a complex combination of increased impact from volcanic tephra deposition, cooling climate, and increased sea ice off the coast of Iceland. All lake records indicate a strong decline in temperature ∼1.5 ka, which culminated during the Little Ice Age (1250–1850 CE) when the glaciers reached their maximum Holocene dimensions.

1 Introduction

A growing number of proxy-based climate reconstructions provide essential information for understanding the patterns and mechanisms behind millennial-to-centennial-scale climate variability during the Holocene epoch. Most Northern Hemisphere reconstructions illustrate warmer-than-present conditions during the Holocene thermal maximum (HTM; 11–6 ka) but display spatiotemporal heterogeneity. For example, the HTM in the northern North Atlantic may have been delayed by up to 3000 years compared to the western Arctic (Kaufman et al., 2004). This asynchrony has been attributed to the lingering effects of the residual Laurentide Ice Sheet (LIS) and the Greenland Ice Sheet (GIS), which affected surface energy balances and impacted ocean circulation due to variable meltwater fluxes (Barber et al., 1999; Renssen et al., 2009, 2010; Blaschek and Renssen, 2013; Marcott et al., 2013; Sejrup et al., 2016). Although the primary forcing of the Northern Hemisphere's subsequent cooling was the spatially uniform decline in summer insolation (Berger and Loutre, 1991), orbital forcing also led to nonuniform changes within the climate system (e.g., Wanner et al., 2008). On decadal to multi-centennial timescales, Holocene climate variability is linked to forcing factors such as solar activity and large tropical volcanic eruptions, in addition to coupled modes of internal variability such as the North Atlantic Oscillation (NAO), and changes in the Atlantic Meridional Overturning Circulation (AMOC) and its capacity to transport heat in the North Atlantic. Feedbacks between ocean, atmosphere, sea ice, and vegetation changes also complicate the global response to primary insolation forcing and result in nonlinear climate variability (Trouet et al., 2009; Miller et al., 2012; Lehner et al., 2013; Moreno-Chamarro et al., 2017).

A recent synthesis showed that a general trend of mid Holocene glacier growth in the Northern Hemisphere was a response to declining summer temperatures driven by the orbitally controlled reduction in summer insolation (Solomina et al., 2015). However, spatially divided compilations based on the Arctic Holocene Transitions (AHT) database (Sundquist et al., 2014) indicate heterogeneous responses to the cooling between Arctic regions (Briner et al., 2016; Kaufman et al., 2016; Sejrup et al., 2016). The response of the northern North Atlantic region is considered to be particularly complex. Records from this region show a strongly uniform pattern of sea surface and terrestrial summer temperature responses to changes in summer insolation but also signatures of subregional differences that were apparently related to both the strong influence of the LIS's disintegration and meltwater discharge in the early Holocene and the variability in AMOC strength throughout the remainder of the Holocene (Sejrup et al., 2016).

Figure 1 Location of Iceland in the northern North Atlantic relative to the major ocean surface currents. EGC: East Greenland Current; EIC: East Iceland Current; NIIC: North Icelandic Irminger Current; IC: Irminger Current; NAC: North Atlantic Current.

Iceland lies in the northern North Atlantic, a region strongly affected by the northward heat transport of AMOC (Fig. 1), and within the periphery of the strongest contemporary warming (Overland et al., 2016). Paleoclimate research across Iceland has shown that positive summer insolation anomalies of the early Holocene resulted in an “ice-free” Iceland by ∼9 ka (Larsen et al., 2012; Striberger et al., 2012; Harning et al., 2016a). Glacier modeling experiments suggest that the observed disappearance of Langjökull in central Iceland by 9 ka required summer temperature 3 ^∘C above the late 20th-century average (Flowers et al., 2008), an estimate reinforced by simulations for Drangajökull, an ice cap in northwest Iceland (Anderson et al., 2018). Reconstructions from two high-resolution lake sediment records in Iceland show that summertime cooling during the Holocene had commenced by ∼5 ka (Geirsdóttir et al., 2013). This cooling was likely intensified by lower sea surface temperatures around Iceland and increased export of Arctic Ocean sea ice resulting in renewed glacier nucleation in the mid Holocene (Geirsdóttir et al., 2009a; Larsen et al., 2012; Cabedo-Sanz et al., 2016). The high-resolution lacustrine records indicate that, despite the monotonic decline in summer insolation, Iceland's landscape changes and ice cap expansions were nonlinear with abrupt changes occurring ∼5.0, ∼4.5–4.0, ∼3.0, and ∼1.5 ka (Geirsdóttir et al., 2013), with Iceland's ice caps reaching maximum dimensions during the Little Ice Age (LIA, ∼1250–1850 CE; Larsen et al., 2015, Harning et al., 2016b; Anderson et al., 2018). In order to understand the nonlinear pattern and stepped environmental and climate changes in Iceland after the HTM and how regional temperatures evolved in terms of timing, magnitude, and glacier inception, we focus specifically on the climate steps between 6.0 and 3.0 ka. This time interval includes the 4.2 ka aridification and cooling event recognized at many global locations. We place the 4.2 ka event in the context of our Icelandic Holocene climate reconstruction and knowledge of large Icelandic volcanic eruptions as a way of judging if it is indeed a major climate event.

This paper composites six different normalized climate proxies (i.e., $\text{mean}=\mathrm{0}\pm \mathrm{1}$) in high sedimentation-rate cores from seven lakes in Iceland. Here, we add five new sediment records to the original two-lake composite of Geirsdóttir et al. (2013) to more effectively evaluate whether individual site-specific proxy records reflect regional climate changes or catchment-specific processes that are less directly coupled to climate. We employ a statistical comparison of our proxy reconstructions from the seven lakes to quantitatively evaluate their between-lake correlation over time and, by inference, their utility as proxies for regional climate.

2 Study sites

The temperature and precipitation gradients across Iceland reflect the island's position between cold sea surface currents from the Arctic and warmer, saltier sea surface currents from lower latitudes delivered by discrete limbs of the AMOC (Fig. 1). Consequently, climate perturbations that alter the strength, temperature, and/or latitudinal position of these currents should translate into significant changes in the mean terrestrial climate state of Iceland. As an example, sea ice imported and locally formed off the north coast of Iceland during extreme cold events provides a positive feedback on the duration and severity of short-term perturbations from volcanic eruptions (Miller et al., 2012; Sicre et al., 2013; Slawinska and Robock, 2018).

Due to Iceland's location astride the North Atlantic Ridge, frequent volcanism during the Holocene has influenced Iceland's environmental history. The basaltic bedrock is easily erodible, which results in high sedimentation rates within Icelandic lake catchments and the potential for high-resolution climate records. The frequent production of thick tephra layers has also periodically stripped the landscape of vegetation and triggered intensified soil erosion (cf. Arnalds, 2004; Geirsdóttir et al., 2009b; Larsen et al., 2011; Blair et al., 2015; Eddudóttir et al., 2017). Impacts on the climate and regional temperatures from Icelandic volcanism, however, have yet to be documented. Clear evidence exists for widespread soil erosion across Iceland during the past 1 000 years. While the prevailing hypothesis favors human settlement as the trigger for widespread erosion, a combination of natural late Holocene cooling and volcanic tephra deposition also likely played a role (e.g., Geirsdóttir et al., 2009b).

Figure 2 Air photographs of the seven lakes: (a) Tröllkonuvatn TRK, (b) Skorarvatn (SKR), (c) Torfdalsvatn TORF, (d) Haukadalsvatn HAK, (e) Vestra Gíslholtsvatn VGHV, (f) Arnarvatn Stóra ARN, (g) Hvítárvatn HVT. Air photographs (from 2005) are high resolution (0.5 m) from Loftmyndir (http://www.loftmyndir.is/).

The seven lakes studied here lie on a transect crossing the temperature and precipitation gradient from south to northwest (Table 1, Figs. 1, 2). Furthermore, they represent a range of catchment styles, including lowland/coastal lakes, highland lakes, and non-glacial and glacial lakes, and they span elevations between 30 and 540 m a.s.l. The lakes are of varied sizes, from 0.2 to 28.9 km², and lake depths range from 3 to 83 m (Fig. 2, Table 1). In most years, the lakes are ice covered from November to April, although early winter thaws and/or late season ice cover may occur. All lakes contain a high-resolution Holocene record with sediment thickness ranging from 2.5 m (SKR) to > 20 m (HVT) and occupy glacially scoured bedrock basins. Vestra Gíslholtsvatn (VGHV), Arnarvatn Stóra (ARN), Torfdalsvatn (TORF), Haukadalsvatn (HAK), and Skorarvatn (SKR) have not received glacial meltwater since deglaciation and preserve organic-rich sediments deposited throughout the remainder of the Holocene. In contrast, Hvítárvatn (HVT) and Tröllkonuvatn (TRK) are glacial lakes that have received glacial meltwater when ice margins reoccupied lake catchments following early Holocene deglaciation. Shrub heath characterizes the vegetation around most of the lakes, with low-lying plants and mosses dominating the shallow soils. The one exception is VGHV, the southernmost lake, which today is predominately surrounded by large planted hayfields. However, shrub heath grows near bedrock outcrops on elevated areas around VGHV, whereas fen/mires prefer locations where water can accumulate in lower elevation areas. Three of the lakes – VGHV, HVT, and ARN – lie within the current active volcanic zone of Iceland and therefore contain and preserve more frequent and thicker tephra layers compared to the other four lakes that all lie distal to the main volcanic zones (Fig. 1).

Table 1 Lake settings.

3 Methods

Geirsdóttir et al. (2013) composited six different normalized climate proxies in high sedimentation-rate (≥1 m ka⁻¹) cores from two lakes: a glacial (HVT) and a non-glacial (HAK) lake. Despite the different catchment-specific processes that characterized each lake's catchment, the composite proxy records show remarkably similar stepped changes toward colder states after the mid Holocene. The common signal in very different catchments indicates that the climate proxies in the sediment records reliably reflect past climate change in Iceland. Here, we add five new Icelandic lake sediment records and test whether we can replicate the two-lake composite record in other lakes across Iceland. First, we generate a single multi-proxy composite record for each lake by applying the same methodology as in Geirsdóttir et al. (2013). Before calculating each regional mean record, all records were standardized over the whole record/period, which permits a reliable comparison of the records with each other.

3.1 Generating composite records from physical and organic proxies

The seven multi-proxy lake composites were generated following Geirsdóttir et al. (2013). Individual composite records were calculated from normalized values for each of six environmental proxies by subtracting the physical proxies (TOC – total organic carbon, MS – magnetic susceptibility, C∕N), which reflect erosional activity and cold and windy seasons within the catchment, from the organic matter proxies (BSi, δ¹³C, δ¹⁵N), which reflect biological productivity (mainly diatom growth) and summer warmth. As explained in Geirsdóttir et al. (2009b), TOC and C∕N are grouped with the physical proxies (catchment instability) in the composites due to the flux of soil carbon into the lakes during periods of severe soil erosion occurring during cold, dry, and windy seasons. The TOC in the sediment is primarily a product of both production and transport. The production term is increased during warm periods due to increased plant growth, but transport is minimized as organic material in the catchment accumulates and remains sequestered in soils. During cold periods, even though the production term is minimized, the transport of previously accumulated organic matter from eroding soils contributes a large influx of OC (organic carbon) to the lake sediment. This more than compensates for any decrease in productivity due to shorter growing seasons and leads to a net increase in sediment TOC during cold periods. The soils of Iceland lack cohesion and are susceptible to erosion, both through aeolian processes and overland flow (Arnalds, 2004). Of these processes, wind transport of soils is widespread and significant in Iceland, as displayed by characteristic “rofabarð” features (Arnalds, 2000). A comparison of modern winter wind speed and lake sediment shows good correspondence during the instrumental record in northwest Iceland (Geirsdóttir et al., 2009b). We do not discount that soil erosion happens due to overland flow or glacier erosion but conclude that wind is the dominant driver.

All proxies are given equal weight in each lake's multi-proxy composite (Geirsdóttir et al., 2013). Analyseries software (Paillard et al., 1996) was used to resample each proxy time series to the same 20-year increments. The 20-year interval reflects the minimum resolution within the least resolved dataset used in the analysis (i.e., the BSi record). The C∕N ratio reflects the ratio of terrestrially versus aquatically derived carbon in the lake sediments and hence the degree of erosion in the lake catchments. The individual proxy data used here have been published previously, and corresponding methods, geochronology, and interpretations are described in the relevant original publications (Table 1).

3.2 Correlating lake sediment cores

The abundant tephra layers archived in Icelandic lake sediments generally hold diagnostic geochemical fingerprints, which allow them to serve as key chronostratigraphic markers. These tephra horizons have been dated using historical information and radiocarbon measurements of soil, lacustrine, and marine archives. The tephra-based chronologies for each lake sediment sequence, together with synchronization of paleomagnetic secular variation between four of the lakes (HVT, HAK, ARN, TORF) and a well-dated marine sediment core off the coast of northern Iceland, provide robust chronologic control and minimal age uncertainties (<100 years) over the Holocene (e.g., Jóhannsdóttir, 2007; Stoner et al., 2007; Larsen et al., 2012; Ólafsdóttir et al., 2013; Harning et al., 2018a). These detailed chronologies thus allow for the development of regional syntheses of climate history through the direct comparison of the records from nearby lakes (Geirsdóttir et al., 2013). The age model for each lake was constructed by fitting control points with a smoothed spline using the CLAM code (Blaauw, 2010) before resampling to the same 20-year increment. All ages in the text are reported as calendar years prior to 2000 CE (b2k).

3.3 Statistical analyses

Our primary goal when reconstructing Holocene climate evolution is to test whether changes on decade-to-century scales are regionally coherent (e.g., Geirsdóttir et al., 2013). This effort requires continuous records from different catchments where consideration is given to the large number of highly resolved climate proxies derived from each of the seven lake sediment records. We are primarily interested in the agreement or otherwise of the trends in the data. In order to compare both within-lake and between-lake variables, we normalized the resampled data (n=3478) so that (1) all variables had equal weight ($\text{mean}=\mathrm{0}\pm \mathrm{1}$) and (2) differences between lakes, possibly due to climatic gradients, are clarified. In order to evaluate the linkages between the lake proxies, we applied R-Mode factor analysis to the 6×3478 data matrix (Davis, 1986; Aabel, 2016) to evaluate the proxy records as indices of regional climate change rather than site-specific environmental conditions that may be decoupled from regional climatic parameters. Factor analysis extracts the dominant signal (first principal component, PC), and the explained variance plus each lakes factor scores give a measure of how the different sites are associated with this main signal, hence allowing us to evaluate the importance of an individual proxy as a climate recorder.

Figure 3 (a) Comparison of normalized composite proxy records of the seven lakes. Hekla 4 (H4; 4.2 ka) and Hekla 3 (H3; 3.0 ka) are marked with orange-colored, hatched lines. Blue hatched lines mark steps in the records at 5.5, 4.5, and 1.5 ka. (b) Comparison of normalized composite BSi records of the seven lakes. Hekla 4 (H4; 4.2 ka) and Hekla 3 (H3; 3.0 ka) are marked with orange-colored, hatched lines. Blue hatched lines mark steps in the records at 5.5, 4.5, and 1.5 ka.

4 Results and interpretation

4.1 Proxy measurements and the composite Holocene records

The multi-proxy composites for each lake reduce proxy-specific signals within each lake, while amplifying those signals that are recorded by most or all proxies (Fig. 3a). By isolating each individual lake composite record, the signal of catchment-specific processes within each lake record is preserved, which validates our comparison between different catchments to test for a regional, Iceland-wide signal.

Table 2 Factor analyses.

4.2 Results from factor analyses

We employ the statistical comparison between our proxy-based reconstructions from the lakes in order to quantitatively evaluate their correlation over time and, by inference, their utility as proxies for regional climate. The results (Table 2) indicate that the first factor explains 56.9 % of the variance and has high loadings (associations) with five of the proxies. MS and δ¹⁵N are positively loaded on this axis whereas TC (total carbon), BSi, and δ¹³C have strong negative loadings. Communalities, a measure of shared information (Davis, 1986; Aabel, 2016), is high for all six proxies, with δ¹⁵N having the highest unique (noise) rating. An additional 19.5 % of the variance is explained by the second factor, which is largely a measure of the C∕N ratio. Varimax factor analysis (Table 2) does not result in any significant changes in the rankings or sign of the proxies, except for BSi, which now has the strongest loading on factor 1. The first two factors explain 76.4 % of the data, indicating that most of the variance in the dataset may be explained by these factors and that the seven lake sediment records are showing a similar signal.

Figure 4 Factor 1 and 2 explain 76.4 %. See Table 2 for details.

A plot of the factor scores (Fig. 4), labeled for each of the seven lakes, indicates a clear environmental gradient along the first-factor axis and a less distinct, but still clear, division on the second-factor axis, for example between HVT and HAK (Fig. 4). These distinct clusters reflect the combined scaling of each 3478 measurements per proxy.

Figure 5 (a) Composite records of all seven lakes combined compared to the two-lake composite of Geirsdóttir et al. (2013), (b) BSi composite of all seven lakes compared to composite of all proxies in seven lakes, (c) BSi all lakes, and C∕N all lakes. Dashed lines mark the step changes identified.

The similarities in timing and direction in variations of climate proxies preserved in all seven composites, our robust age models, and the communalities (Table 2) enable us to consider all seven composite lake records and generate a single time series (Fig. 5a). Compared to our previous composite (Geirsdóttir et al., 2013), the seven-lake composite is near-identical indicating that all lakes experienced similar climatic histories with relatively minor superimposed catchment-specific processes (Fig. 5a). Furthermore, the results from the factor analyses suggest that in addition to a seven-lake all-proxy composite, a simple seven-lake composite of just BSi (relative spring/summer warmth; e.g., Geirsdóttir et al., 2009b) and C∕N (cold or erosional activity) can increase our understanding of the evolution of Holocene climate in Iceland and disentangle the catchment responses (C∕N) to climate change and the temperature forcing (BSi). A comparison of such reconstructions with the composite of all proxies/all lakes is shown in Fig. 5b and c.

5 Discussion

5.1 The 4.2 ka event in Iceland and the role of volcanism

The overall (first-order) trend of both the seven-lake all proxies and BSi composites track the orbital cycle of summer insolation across the Northern Hemisphere after peak insolation ∼11 ka (Fig. 5b, c). Superimposed on the first-order decrease in composite BSi-inferred relative temperature are apparent step changes at ∼5.5, ∼4.5, ∼3, and ∼1.5 ka, with the lowest temperatures culminating during the LIA between 0.7 and 0.1 ka. However, the prominent step ∼4.5–4.0 is the first distinct step that can be correlated between all seven lake records (Fig. 3b). This period spans the well-known ∼4.2 ka event identified as a period of abrupt climate change elsewhere around the globe, manifested by pronounced dry conditions resulting in severe societal collapse in the Eastern Mediterranean (e.g., Weiss, 2016, 2017; Weiss et al., 1993; Cullen et al., 2000) and cooler and/or wetter conditions at higher altitudes in the Alps (Zanchetta et al., 2011, 2016). However, depending on the proxy and/or site chosen, the radiometric ages constraining these records from outside of Iceland indicate a much broader interval from 4.5 to 3.5 ka (Gasse, 2000; Wang et al., 2016; LeRoy et al., 2017; Railsback et al., 2018). Discerning whether the period of major climate in Iceland between 4.5 and 4.0 ka relates to the global expressions of the 4.2 ka climate event remains an open question.

Two steps toward cooling in the Icelandic lacustrine records coincide with two of the largest explosive eruptions of the Holocene in Iceland; Hekla 4 (3826±12¹⁴C BP/4197 cal yr BP; Dugmore et al., 1995) and Hekla 3 (2879±34¹⁴C BP/3006 cal yr BP; Dugmore et al., 1995). Both these Hekla eruptions deposited at least 1 cm of tephra over 80 % of the surface of Iceland and are important tephrochronological markers throughout Europe and, importantly, a tie point between most of our lake sediment records. Because tephra fallout from these eruptions likely disrupted local ecosystems (Larsen et al., 2011; Christensen, 2013; Eddudóttir et al., 2017), it is possible that they artificially imply colder conditions in Icelandic lake sediment records from induced terrestrial erosion (higher C∕N), which we would normally attribute to periods of cold and windy winters. However, Hekla 4 and Hekla 3 tephra are not found in either SKR or TRK (two lakes around Drangajökull in northwest Iceland), and thus the proxy records from these two lakes provide unequivocal evidence for cooling at these times unrelated to tephra-induced catchment disturbance or soil erosion (Harning et al., 2018a, b).

In terms of potential regional climate and temperature changes induced by the Hekla 4 eruption, the bulk of the tephra emitted was fine to very fine ash with a residence time in the atmosphere presumably on the order of days to weeks, which together can lead to a significant impact on a local to regional scale on timescales ranging from weeks to ∼3 years. The volume of the tephra is estimated to be 13.3 km³ (1.8 km³ when calculated as dense rock) deposited from a >30 km high eruption column, which makes it a category VEI5 event (Stevenson et al., 2015). Although the aerosol cloud produced by the Hekla 4 eruption had at least a hemispheric coverage, its impact on atmospheric properties and processes would have been negligible as data indicate that the Hekla 4 magma carried only about 20–70 ppm SO₂ (e.g., Portnyagin et al., 2012), compared to, for example, the Laki eruption of 1783–1784 in Iceland, which carried 1600 ppm (Thordarson et al., 1996). On the other hand, Hekla 4 magma would have carried between 500 and 1000 ppm Cl (Laki, 310 ppm) and 1300–1600 ppm F (Laki, 660 ppm). Although the atmospheric venting of sulfur was negligible (<0.5 Tg), the halogen emissions would have been significant or on the order of 10–15 Tg (Sverrisdóttir et al., 2007; Portnyagin et al., 2012; Stevenson et al., 2015). The halogens released by the Hekla 4 eruption into the stratosphere may thus have resulted in hemispheric depletion of ozone and hence stratospheric cooling on the same scale. However, the current unknown is how that cooling translates into surface cooling and associated changes in the tropospheric weather systems.

Although the origin and climatic impact of the 4.2 ka event remains poorly understood, it has previously been linked to both volcanism (e.g., Antoniades et al., 2018) and ocean–atmospheric circulation changes in the North Atlantic (Bond et al., 2001; Bianchi and McCave, 1999; deMenocal, 2001). The similarity in the timing of cool events at higher latitudes in the Northern Hemisphere at this time further supports a temporary climate shift around the North Atlantic, most likely associated with millennial-scale variability in the AMOC and/or the subpolar gyre (e.g., Risebrobakken et al., 2011; Thornalley et al., 2009; Trouet et al., 2009; Orme et al., 2018; Moreno-Chamarro et al., 2017; Zhong et al., 2018). Potential contributions to local and/or hemispheric climate changes resulting from Hekla 4 magma degassing, however, remain a topic of further research.

Figure 6 (a) BSi and C∕N from HVT and Lögurinn. (b) BSi and C∕N from Lögurinn; (c) BSi from TRK. Blue shaded bars indicates glacier inception at each lake catchment.

5.2 The demise and growth of the Icelandic glaciers during the Holocene

Two of our lakes, HVT and TRK, are currently glacial lakes allowing their sediment records to track the activity of Langjökull and Drangajökull during periods when these ice caps attain sufficient dimensions to occupy the respective lake catchments (Figs. 2, 6). Hence, climate proxies in these lake sediments record the growth and demise of an upstream glacier (e.g., Briner et al., 2010; Larsen et al., 2012; Harning et al., 2018b). Both lake sediment records demonstrate that glacier ice disappeared from their catchments prior to 9 ka (Larsen et al., 2012; Harning et al., 2016b, 2018b). Furthermore, detailed sedimentological analyses of HVT indicate that the lake did not receive glacial meltwater between 7.9 and 5.5 ka (Larsen et al., 2012). Glacier modeling experiments show that the demise of Langjökull and Drangajökull during the early Holocene required summer temperatures to rise ∼3 ^∘C above the late 20th-century average at this time (Flowers et al., 2008; Anderson et al., 2018).

At 5.5 ka, sedimentological analyses of HVT (increase in sediment accumulation rate, higher MS values, and diminished biological productivity) show the first abrupt change toward cooler conditions and glacier occupation of the lake catchment (Larsen et al., 2011; Fig. 3a, b). The hydrologically coupled ice-sheet model employed to simulate the evolution of Langjökull through the Holocene also captures the inception of the modern ice cap prior to 5 ka (Flowers et al., 2008), concurrent with the first abrupt change in the BSi, MS, and the sediment accumulation rate records (Fig. 6a). The second distinct change in HVT between 4.5 and 4.0 ka towards cooler climate is not only reflected in a suppressed algal productivity and increased landscape instability but a shift from weakly stratified to finely laminated sediments signaling the onset of a glacial–lacustrine-dominated HVT catchment. The Hekla 4 tephra layer generated from the explosive Hekla eruption dated to ∼4.2 ka marks the culmination of this second change towards cooler climate that started around 4.5 ka. At 3 ka HVT developed the first distinctly varved sediment record, reflecting increased glacier activity in the catchment (Larsen et al., 2011). This again coincides with the eruption of the Hekla volcano at 3 ka and deposition of Hekla 3 tephra. In both Hekla cases, the impact of the tephra in the HVT catchment has previously been linked to increased soil erosion as reflected in the C∕N record (Larsen et al., 2011, 2012), which may obscure the impacts of global cooling related to the 4.2 ka event at the time of the Hekla 4 eruption. However, the potential impacts from Hekla 4 magma degassing may have contributed to the nonlinear cooling step resulting in further growth of Langjökull under a background of decreasing summer insolation. The steepest decline in algal productivity (BSi in HVT) and the largest rate of glacier advance (Larsen et al., 2011) began later, between 1.8 and 1.5 ka and culminated during the LIA (Fig. 6a). This pattern is supported by the glacier numerical simulations, which suggest Langjökull attained its maximum volume during the Little Ice Age; first around 1840 CE and then around 1890 CE when estimated temperature were 1.5 ^∘C below the 1960–1990 average (Flowers et al., 2007) (Fig. 6c).

Unlike HVT, glacier ice did not reach Tröllkonuvatn's catchment until 1 ka, as illustrated by the lack of meltwater-derived clastic sediment in the lake between 9 and 1 ka (Harning et al., 2018b). However, Drangajökull was actively expanding into another threshold lake catchment on its current southeastern margin by at least 2.3 ka (Harning et al., 2016a). Further evidence for late Holocene advances of Drangajökull come from ¹⁴C dates on dead vegetation emerging from the currently receding northern and southern ice margins, which define the timing of persistent ice margin expansion at these locations (Harning et al., 2016a, 2018b). The threshold lake records combined with ¹⁴C-dated emerging dead vegetation define five periods of increasing ice cap dimensions at ∼2.3, 1.8, 1.4, 1, and 0.5 ka, where the final two ice margin advances (1 and 0.5 ka) are interpreted from Tröllkonuvatn's sediment record. These periods of ice growth are matched by decreases in BSi and an increase in C∕N, suggesting that regional climate change is the primary controlling factor (Harning et al., 2018b). The numerical simulations by Anderson et al. (2018) support Drangajökull's late Holocene appearance, subsequent expansion, and maximum dimensions between 0.5 and 0.3 ka, with temperatures likely 1.0 to 1.2 ^∘C below the 1960–1990 reference temperature.

Striberger et al. (2012) studied the sediments of the glacier-fed Lake Lögurinn in eastern Iceland to infer the Holocene meltwater variability of Eyjabakkajökull, a surge type outlet glacier of NE Vatnajökull (Figs. 2, 6b). The results show that Eyjabakkajökull had ceased to deliver glacial meltwater to the lake by 9.0 ka and that the glacier's dimensions were considerably smaller relative to today at this time. The return of glacial meltwater into the lake at ∼ 4.4 ka indicates the regrowth of Eyjabakkajökull into Lake Lögurinn's catchment and an almost 5000-year-long glacier-free period along NE Vatnajökull during the early to mid Holocene. After 4.4 ka, the Lake Lögurinn BSi record reflects a marked and continuous decrease in the aquatic productivity through to the LIA, which we suggest reflects continued late Holocene summer cooling. Although increasing glacially sourced sediment can dilute the BSi signal, either interpretation (decreasing BSi or increasing glacial sediment) or some combination thereof suggests that Eyjabakkajökull attained its maximum Holocene extension during the LIA as well (Striberger et al., 2012).

The results from the factor analyses illustrate a clear gradient along the first-factor axis (BSi) and a less distinct, but still clear, division on the second-factor axis (C∕N). The fact that this trend is especially evident between HVT, the glacial lake in the highlands, and HAK, the non-glacial coastal lake in western Iceland, supports our previous interpretation (e.g., Geirsdóttir et al., 2013) that these two lakes form two end-members of combined proxy responses to the mid to late Holocene forcing toward cooler summers. Great similarities occur between the two highland lakes (HVT and ARN), which both lie within the volcanic zone of Iceland on either side of Langjökull (Fig. 2), although only HVT is currently affected by the glacier (Larsen et al., 2012; Gunnarson, 2017). Importantly, both records show similar timing and direction of the most distinct perturbations throughout the mid to late Holocene (Fig. 3a, b). Similarly, the coastal lakes HAK, TORF, and SKR share common features with more suppressed punctuations than the highland lakes, perhaps indicating greater impact from the surrounding sea surface temperature and its characteristic thermal inertia. This is supported by the similarity of these records to diatom-based sea surface temperatures (SSTs) from the shelf north of Iceland, particularly with regard to how pronounced cooling was delayed until the late Holocene (Fig. 7; Jiang et al., 2015). TRK, although also a glacial lake, shows more similarity with the coastal lakes than glacial lake HVT (Fig. 3a, b). This is likely due to its shared coastal location outside the active volcanic zone with minimal influence from most volcanic eruptions (Harning et al., 2018b). The lake most affected by volcanism in Iceland is VGHV in the southern lowlands (Blair et al., 2015), and thus it seems to form its own group in the factor analyses (Fig. 3a, b and Fig. 4).

Figure 7 Composite record from all seven lakes compared to (a) 65^∘ N summer insolation (Berger and Loutre, 1991). (b) Volcanic forcing W m⁻² (Kobashi et al., 2017), (c) MD99-2275 diatom-based SST (Jiang et al., 2015); (d) Renland–Agassiz ice core record (normalized δ¹⁸O record and its variance) from Greenland (Vinther et al., 2009); (e) quartz (normalized) from MD99-2269 as an indication of sea ice in the northern North Atlantic (Moros et al., 2006); (f) IP25 record of sea ice (Capedo-Sanz et al., 2016); (g) BSi coastal lakes SKR, TRK, TORF, and HAK; (h) BSi highland lakes HVT+ARN; (i) BSi all lakes. Blue hatched lines show the 5.5, 4.5, and 1.5 ka steps in the Icelandic lake records. Orange hatched lines show the timing of Hekla 4 and Hekla 3 volcanic eruptions.

An explanation for the step shifts in the BSi (and other climate proxies) may in part be related to perturbations from aerosol-induced cooling from Icelandic and/or tropical volcanism, which would likely produce brief cold summers resulting in depressed algal productivity and enhanced catchment instability (e.g., HVT, ARN, VGHV; Fig. 4). But with all climate proxies (all composites) exhibiting severe perturbations at the same time (although of variable magnitude), their failure to return to pre-perturbation values suggests a changed relation between the landscape and/or climate following each perturbation. The different response of the highland/glacier record compared to the coastal lake records suggests differences in the response of catchment-specific processes to the orbitally forced summer temperature decline, which in turn depends on geographic location and the impact of the surrounding water (Fig. 6). Although the lake sediment record from the southern lowlands (VGHV) deviates from this, it records brief proxy perturbations following major Icelandic eruptions but shows limited sustained instability (Blair et al., 2015). The behavior of VGHV's record may suggest the influence of intermittent volcanic activity on the proxies against a backdrop of generally consistently “warm” climate along the southern lowlands induced by the Irminger Current (IC; Fig. 1) flowing from east to west south of Iceland. On the other hand, the coastal lakes in northwest Iceland may be affected by variations in SSTs of the cooler and more variable North Iceland Irminger Current (NIIC; Fig. 1). However, noteworthy here is the synchronous and accelerated decline after 1.5 ka found within all lake records and the concomitant abrupt increase in sea ice along the North Iceland Shelf (Fig. 7f).

5.3 The “onset” of neoglaciation

Comparing the first indication of temperature lowering and glacier formation after the HTM in terrestrial regions in and around the northern North Atlantic reveals certain similarities and shows regionally consistent increases in millennial-scale cooling rates. The monotonic Holocene decline in Northern Hemisphere summer insolation was most likely the primary driver for the first expansion of the cryosphere at ∼5 ka identified in Baffin Island, eastern Greenland, in the highlands of Iceland, western Svalbard, and in western Norway (Funder, 1978; Nesje et al., 2001; Jennings et al., 2002; Masson-Delmotte et al., 2005; Bakke et al., 2005a, b, 2010; Vinther et al., 2009; Larsen et al., 2012; Balascio et al., 2015; Solomina et al., 2015; Røthe et al., 2015, 2018; van der Bilt et al., 2015; Gjerde et al., 2016; Briner et al., 2016; Miller et al., 2017). Although the gradual decline in summer insolation progressively lowered the equilibrium line altitude (ELA), the significant stepwise trend apparent in the Icelandic records and other records around the North Atlantic suggests that strong local to regional feedbacks modulated the primary insolation forcing. The rate of cryosphere expansion at 4.5-4.0 ka and particularly after 1.5 ka documents contemporaneous shifts in the northern North Atlantic region. Such episodic ice expansion cannot be explained by the summer insolation forcing alone and requires additional forcing or internal climate variability. Variations in the strength of the thermohaline circulation, weakening of the northward heat transport of the AMOC, and/or increasing influence of the Arctic waters influence all these locations. Changes in the strength of AMOC and/or the subpolar gyre and changes in the Arctic sea ice extent with the associated meridional heat transport into the Arctic have been related to past cooling events, particularly during the last 2 kyr (Trouet et al., 2009, 2012; Lehner et al., 2013; Cabedo-Sanz et al., 2016; Moreno-Chamarro et al., 2017; Zhong et al., 2018).

The time interval known as “neoglaciation” was defined by Porter and Denton (1967) as “the climatic episode characterized by rebirth and/or growth of glaciers following maximum shrinkage during the Hypsithermal [now HTM] interval”. Porter (2000) noted that neoglaciation is a geologic-climate unit based on physical geological evidence of glacier expansion and that palynological evidence of climate change was excluded from the definition. The use of the phrase “onset of neoglaciation” in the literature has been broader than the original definition and is commonly attributed to the first apparent indication of lowered temperature or increased rate of summer cooling rather than simply the renucleation and/or expansion of glaciers.

Following Porter's (2000) definition, the onset of neoglaciation in Iceland based on our lake records occurred before 5.0 ka for Langjökull, despite the initial growth of Drangajökull occurring much later at ∼2.3 ka. This indicates that the spatiotemporal nucleation of glaciers in Iceland was asynchronous and likely reflects the relation between the regional ELA (primarily controlled by summer temperature) and topography. The nature of the topography (i.e., the hypsometry) controls how quickly the glacier will expand after the ELA intersects the topography. Because the topography under Drangajökull is a plateau-like landscape with a large area at its highest elevation, Drangajökull will grow quickly once the ELA intercepts the topography (Anderson et al., 2018).

The current ELA pattern reflects the patterns of temperature and precipitation across Iceland. Temperature differences from the south to the northwest reflect the prevailing wind direction and proximal ocean surface temperatures. Precipitation also impacts ELAs across Iceland but varies primarily due to the interaction of local topography with prevailing winds (Crochet et al., 2007). Because Icelandic glaciers are most sensitive to temperature, we expect the timing of glacier inception to be controlled by the rate at which Holocene temperatures decline following the HTM and the subglacial topographic setting of each of the current Icelandic ice caps.

Changes in the ELA at present are most sensitive to summer temperature, and assuming the same holds for most of the Holocene, we apply a simple approach to define the “onset” of neoglaciation in Iceland by comparing the Holocene evolution of BSi (a measure of variations in summer temperature) and C∕N (a measure of variations in catchment stability, which independently tracks summer temperature, modulated by volcanism) in the coastal lakes (SKR, TRK, TORF, HAK) and the highland lakes (HVT, ARN) (Fig. 7). Although the multi-proxy lake records document complex changes in terrestrial climate and glacier fluctuations in Iceland during the mid to late Holocene, coherent patterns of change are apparent. Based on glacier nucleation, our records show that initiation of neoglacial cooling took place in the highlands of Iceland (Langjökull) around 5.5 ka, where the current ELA is ∼1170 m (Figs. 6a, 7). This cooling is also mirrored by the retreat of woodland from 6000 to 4000 cal yr BP, reflected by decreased Betula pubescens pollen counts in a lake record from the northwest highlands (Eddudóttir et al., 2016). The rate of glacier growth likely increased between 4.5 and 4.0 ka, when NE Vatnajökull (ELA 1320 m) nucleated (Fig. 6b), and continued near ∼2.5 ka, when Drangajökull nucleated (ELA 675 m) (Fig. 6c). Consistent with this scenario is the distinct first-order cooling trend apparent in all seven lakes, including the combined highland (HVT, ARN) and the coastal (SKR, TRK, TORF, HAK) lake composites (Fig. 7). The proxies in the highland lakes, HVT and ARN, show a more abrupt and greater response to the temperature lowering, most likely due to a greater impact from catchment-specific processes (tephra deposition and/or glacial activity), whereas the coastal lakes show a more subdued response, likely reflecting the moderating effect from SSTs (Fig. 7). These lakes show striking similarities not only with the diatom-based SST record from the shelf north of Iceland (Jiang et al., 2015) but also with the IP₂₅-based sea ice reconstruction (Cabedo-Sanz et al., 2016) and the ice-rafted debris (IRD) record based on quartz grain counts (Moros et al., 2006) (Fig. 7), both off the north coast of Iceland. The IP₂₅-based sea ice reconstruction shows a rise at ∼5 ka from a background state beginning earlier at 8 ka and intensifies after 4.5 ka, broadly in line with the decreased abundance of planktic diatoms and lowering of SST (Jiang et al., 2015) (Fig. 7). Further increases in drift ice were evident during the late Holocene after ∼3.3 ka, with maximum sea ice after ∼1.0 ka and during the Little Ice Age (Moros et al., 2006; Cabedo-Sanz et al., 2016). The intensification around 4.5 ka seen in our lake records is in line with increased strength of cold and fresh polar water via the East Greenland Current (EGC; Fig. 1) at 4.5 ka in the northern North Atlantic inferred from the decreased abundance of planktic foraminifera (e.g., Andersen et al., 2004; Jennings et al., 2011; Ólafsdóttir et al., 2010; Kristjánsdóttir et al., 2016; Perner et al., 2015, 2016).

A diatom-based SST reconstruction from a core retrieved from the Iceland Basin south of Iceland shows pronounced SST cooling between 4 and 2 ka, with warmer temperatures prior to 4 ka but also between 2 and 1.5 ka (Orme et al., 2018). Orme et al. (2018) primarily explain the cool interval between 4 and 2 ka as a response to a persistently negative mode of the North Atlantic Oscillation (NAO) that caused strengthening of the northerly winds east of Greenland, which in turn strengthened the East Greenland Current, bringing cool Arctic water as far south as the Icelandic Basin. The abrupt increases in IP₂₅ at ∼1.5 and 0.7 ka are coincident with the increased rate of cooling identified in the Icelandic lacustrine temperature record, suggesting significant coupling between the marine and terrestrial systems in Icelandic waters. From this treatment, we conclude that the primary driver explaining changes in our climate proxies was the decline in summer insolation modulated by changes in ocean circulation and associated SST around the coast of Iceland. Explosive Icelandic volcanism, whether or not it caused the temperature decline, produced important perturbations to lake catchments through tephra deposition, often resulting in the recovery of catchment systems to a different equilibrium state.

6 Conclusions

The results from factor analyses of six climate proxies from seven Icelandic lake sediment records suggest that a simple composite of BSi (relative spring/summer temperature) and C∕N (cold or erosional activity) serve to increase our understanding of the evolution of Holocene climate in Iceland and disentangle the catchment responses to climate change and the temperature forcing.

The Holocene thermal maximum was warm enough to result in a mostly ice-free Iceland by 9 ka. The regular decline in summer insolation after 11 ka, plausibly amplified by responses elsewhere in the North Atlantic region, led to decreases in summer temperature and the destabilization of catchments in non-glacial lakes.

Based on our seven lake records, the onset of neoglaciation is registered in the highlands of Iceland (Langjökull) around 5.5 ka, when its ELA intercepted the local subglacial topography. The delayed nucleation of northeast Vatnajökull (4.4 ka) and Drangajökull (2.3 ka) can be explained by lower topographic thresholds and the necessity of lower summer temperatures for regional ELAs to intercept their respective topography. Our observed glacier nucleations all coincide with the stepped cooling reflected in all seven lake sediment records.

Episodic glacier expansion between 4.5 and 4.0 ka cannot be explained by the summer insolation forcing alone and require additional forcings, likely linked to ocean circulation and/or expansion of Arctic Ocean sea ice or internal climate variability. Magma degassing from local Icelandic eruptions may also have contributed but requires further research to confirm.

The prominent step toward cooling at 4.5–4.0 ka is statistically indistinguishable from the global ∼4.2 ka event and also coincides with Hekla 4, one of the largest explosive eruptions of the Holocene in Iceland. Although the ash generated from the eruption likely disrupted local catchment stability in Iceland, its impact on global atmospheric properties and processes would have been negligible due to the low sulfur content of the eruption.

North Atlantic circulation features (NAO and AMOC) likely influenced the temperature decline in Iceland during the late Holocene, modulating the hemispherically symmetric forcings (insolation, irradiance, and volcanism).