Glacier response to North Atlantic climate variability during the Holocene

Introduction Conclusions References

only sparse and fragmentary evidence for higher frequency variations in glacier size because in many Northern Hemisphere regions glacier advances of the past few hundred years were the most extensive and destroyed the geomorphic evidence of ice growth and retreat during the past several thousand years. Thus, most glacier records have been of limited use for investigating centennial scale climate forcing and feedback mechanisms. Here we report a continuous record of glacier activity for the last 9.5 ka from southeast Greenland, derived from high-resolution measurements on a proglacial lake sediment sequence. Physical and geochemical parameters show that the glaciers responded to previously documented Northern Hemisphere climatic excursions, including the "8.2 ka" cooling event, the Holocene Thermal Maximum, Neoglacial cooling, and 15 20th Century warming. In addition, the sediments indicate centennial-scale oscillations in glacier size during the late Holocene. Beginning at 4.1 ka, a series of abrupt glacier advances occurred, each lasting ∼ 100 years and followed by a period of retreat, that were superimposed on a gradual trend toward larger glacier size. Thus, while declining summer insolation caused long-term cooling and glacier expansions during the late

Introduction
Glaciers and ice caps represent a small but important portion of the cryosphere (∼ 785 000 km 2 ; Dyurgerov and Meier, 2005). Their mass wasting during the 20-21st century is responsible for 60 % of the sea-level rise unattributable to ocean warming (Meier, 2007) and they continue to retreat at an exceptional rate (Zemp et al., 2012).
Moreover, because small glaciers and ice caps respond rapidly to climate changes and there is a strong relationship between glacier mass balance and summer temperature (Oerlemans, 2005), past glacier extent can inform us about past climate variability.
Holocene glacier activity in the Arctic is reasonably well documented at millennial timescales (Miller et al., 2010). Northern Hemisphere glaciers receded in the early Holocene and were smaller than present during the mid-Holocene. Centennial-scale variations, however, are not well constrained because there are few high-resolution and continuous records, and because in many regions the most extensive glacier advances since the early Holocene took place within the past few hundred years, destroying evidence of intervening glacier positions. This is generally the case in Greenland, 15 where historical (AD 1200(AD -1940 advances of local glaciers were generally the most extensive since at least the early Holocene (Kelly and Lowell, 2009).
Evidence from the Greenland Ice Sheet (Kobashi et al., 2011), marine sediments (Bond et al., 1997;Thornalley et al., 2009;Moffa-Sánchez et al., 2014;Jiang et al., 2015), and terrestrial archives (D'Andrea et al., 2011;Larsen et al., 2012;Olsen et al., 20 2012) indicate that abrupt changes in atmospheric circulation and ocean dynamics, including abrupt cooling events, have punctuated the Holocene. These episodes have alternately been attributed to solar variability, freshwater forcing, volcanic activity, and/or changes in Atlantic Meridional Overturning Circulation (Wanner et al., 2011). How sensitive were glaciers to these abrupt episodes, and did glaciers throughout the North 25 Atlantic respond uniformly? During the period from AD 1250-1900, often referred to as the "Little Ice Age," well-resolved records from the North Atlantic region suggest coherence in ice cap activity that was potentially driven by volcanic activity coupled with sea-ice/ocean feedbacks . However, prior to the last 1000 years there are sparse data for use in investigating the synchrony of glacier response to climate variability in the North Atlantic region. Lakes that receive meltwater from temperate glaciers can be used to develop continuous records of glacier activity. Bedrock erosion at the base of glaciers provides 5 sediment supply for meltwater transport to proglacial lakes. In catchments where other sources of sediment are limited, such as from mass wasting or the release of stored sediment, there is a strong relationship between sediment properties and glacier size (Nesje et al., 2000;Dahl et al., 2003;Jansson et al., 2005), which also follows from the assumption that large glaciers produce more minerogenic material and meltwater than 10 small glaciers. Measurements of physical and geochemical properties of proglacial lake sediments can therefore be used to reconstruct records of past glacier size. Here we report a continuous 9.5 ka record of glacier activity on Kulusuk Island, southeast Greenland developed using sediment cores recovered from Kulusuk Lake (65.56 • N, 37.11 • W; 202 m) (Fig. 1).

Study site
Kulusuk Lake (0.8 km 2 , 69 m maximum depth) is located below a cirque with two small glaciers and is within a low arctic maritime region (MAT −1 • C, MAP 900 mm). Characteristic erosional features indicate that local glaciers have temperate thermal structures (Humlum and Christiansen, 2008). Distinct moraines defined by sharp crests are 20 located in front of both glaciers and the recently glaciated area accounts for ∼ 50 % of the catchment, which is composed of Archaen gneisses (Bridgwater, 1976) (Fig. 1). Kulusuk Lake is ideally situated to capture and preserve a clear sedimentary record of glacier activity because: (i) it only receives runoff from a very small catchment (catchment : lake area ratio of ∼ 2 : 1), minimizing the potential for long-term storage of sed-Introduction and (iii) the small size of the glaciers makes them sensitive to minor climate variations ( Fig. 1). Therefore, bedrock erosion by the glaciers provides the primary source of minerogenic sediment to the lake and changes in glacier size should clearly be reflected in sediment properties.

Sediment core collection and analysis
Sediment cores were recovered from Kulusuk Lake in April 2010 when the lake was ice covered. Bathymetric measurements were made manually through holes drilled in the ice and sediment cores were collected using Uwitec gravity and percussion coring devices from the deepest location, which has a water depth of 69 m. A composite 3.5 m 10 record was compiled by matching the physical stratigraphy and scanning XRF profiles from a 26 cm gravity core (Kul-10D-B) and multiple overlapping percussion cores (Kul-10G-A1, -B1, -A2). The magnetic susceptibility of the cores was measured every 0.5 cm using a Bartington MS2E sensor. The organic-matter content of the cores was measured by loss- 15 on-ignition (LOI) on contiguous 1 cm 3 samples taken at 1 cm intervals. Organic-matter content was calculated as the difference between the weight of dried 1 cm 3 samples and their weight after heating for 4 h at 550 • C (Dean, 1974). Bulk density measurements (g cm −3 ) and the calculated sedimentation rates (cm year −1 ) were used to determine mass accumulation rates (MAR; g cm −2 year −1 ). Grain size measurements were 20 made at 10 cm increments. Samples were pre-treated with a 30 % hydrogen peroxide solution to digest organic material and analyzed using a Beckman Coulter LS200 particle-size analyzer.

Chronology
An age-depth model was established based on 210 Pb analysis of the upper sediments and AMS radiocarbon dates on macrofossils. The 210 Pb activity of samples taken every 1 cm from the upper 10 cm of the record was measured by Flett Research Ltd.
(Winnipeg, Canada) and ages were modeled from these data using a constant rate of 5 supply model. AMS radiocarbon measurements were made on plant/wood fragments and Daphnia ephippia that were wet sieved from core samples. All radiocarbon ages were calibrated to calendar years using CALIB v. 6.0 (Stuiver and Reimer, 1993) with the IntCal09 calibration dataset (Reimer et al., 2009). Ages are presented in calendar years prior to AD 1950 (BP) unless otherwise indicated.

Scanning X-ray fluorescence
To characterize minerogenic changes at higher resolution and with greater sensitivity, an Itrax ™ XRF core scanner was used to produce profiles of relative elemental compositions (Croudace et al., 2006). The Itrax ™ continuously scans the surface of sediment cores at sub-mm resolution with a micro X-ray beam (20 mm × 100 µm) and the relative 15 concentrations of a range of elements are determined based the detection of dispersive energy spectra. Dispersive energy spectra are acquired across each measured interval and peak area integrals are calculated for each element. Peak area integrals are related to elemental concentrations within the sediment, but can also be influenced by characteristics of the sedimentary matrix and therefore only indicate relative changes 20 in elemental composition (Croudace et al., 2006;Rothwell et al., 2006). Our analysis focused on the elements: K, Ca, Ti, Mn, Fe, Zn, Rb, Sr, which have detection limits that range from 150 to 5 ppm (Croudace et al., 2006). All of the cores were scanned at 200 µm intervals with an exposure time of 10 s, voltage of 30 kV, and current of 55 nA.

Sediment stratigraphy and chronology
The 3.5 m composite sediment record from Kulusuk Lake contains distinct lithologic changes, defined by visual stratigraphy, magnetic susceptibility, organic matter content, and elemental data acquired by scanning XRF. The record can be divided into four 5 lithologic units (Fig. 2). Unit I, 3.0-3.5 m, is a gravelly sand. Unit II, from 3.0-1.8 m, is a massive gray clayey silt. There is an abrupt transition to Unit III, a brown, organicrich sediment that extends from 1.8-1.2 m. Unit IV, 1.2-0 m, is a laminated sequence with frequent sandy layers. Laminations consist of fining-upward sequences and impart strong variability in all datasets. 10 Chronologic data shows that there are significant changes in sedimentation rate that correspond to lithostratigraphic units (Fig. 2). An age-depth model was generated assuming that changes in sedimentation rate occurred at the boundaries of these units. In Unit IV, a third-order polynomial was applied to the radiocarbon ages, the core top date that represents when the cores were collected in AD 2010 (−60 cal yr BP), and the date 15 of the base of the 210Pb profile at 10 cm (111 cal yr BP) (Table 1). This relationship was extrapolated to the base of Unit IV. Linear interpolation between the remaining radiocarbon ages was used to generate the age-depth relationship for Units III and II. There is no chronologic control below 215 cm so we did not interpret sedimentation prior to 9.5 cal ka BP. 20 Magnetic susceptibility, organic matter, and mass accumulation rate (MAR) profiles further define these lithologic changes with higher magnetic susceptibility values across intervals with coarser sediment and with lower organic content (Fig. 3) increase and then display more minor fluctuations across the upper 1.2 m. These intervals are clearly defined by MARs, which incorporate sediment density measurements that range from ∼ 0.8-1.8 g cm −3 , but are primarily controlled by the large sedimentation rate changes (Fig. 3).

5
Elemental scans acquired by scanning XRF show a similar response to magnetic susceptibility with higher values across coarser, clastic intervals. However, the XRF data have a greater sensitivity to minerogenic changes and were measured at higher resolution (0.2 mm) (Fig. 3). We focused our analysis on the elements: K, Ca, Ti, Mn, Fe, Zn, Rb, and Sr, which are common in siliciclastic sediments. Changes in the concentrations of these elements reflect changes in the contribution of minerogenic material eroded from catchment bedrock and delivered to the lake. Statistical analysis of the scanning XRF data indicate that all of the elements are highly correlated and that there is a strong primary trend in the data. Correlation coefficients show the strong significant relationships among the majority of the elements (Table 2). Rather than relying on 15 a single element (e.g., Ti), we used principal component analysis (PCA) to define the leading mode of variability (PC1) among the elemental data. PCA allows for a multidimensional examination of the dataset in order to identify the primary signal(s). PCA results indicate that there is one strong primary trend in the elemental data with the first eigenvector (PC1) accounting for 76 % of the total variance. The factor loadings reveal 20 the high correlations between individual element profiles and PC1 ( Table 2). The trends in PC1 are similar to those in the lower resolution magnetic susceptibility and organic matter content records, justifying use of PC1 data to infer past minerogenic changes (Fig. 3). The choice to use PC1 rather than a single representative element (e.g. Ti) to represent changes in sedimentary minerogenic content has no impact on any of our

Sedimentation in proglacial lakes
Sedimentation in proglacial lakes can be the result of a complex set of physical processes associated with the erosion, storage, and transport of sediment within glacial and proglacial systems (Dahl et al., 2003;Jansson et al., 2005). It is important to con-5 sider these complicating factors when selecting sites for glacier reconstructions and when interpreting sedimentary records. Glaciers fundamentally impact the amount and character of minerogenic sediment in proglacial lakes. Glacier size, erosive ability, and meltwater production directly influence the amount of minerogenic sediment delivered to a proglacial lake. However, sedimentation processes in a proglacial lake can also be 10 impacted by mass wasting processes in paraglacial environments (particularly in landscapes with steep unstable slopes), and by the delayed release of sediment stored along the transport pathway between the glacier and the lake (for example, sediment stored in extensive meltwater stream channels). Relative to minerogenic material, organic sedimentation is typically a minor component in proglacial lakes and is related to 15 the input of organic matter from autochthonous and allochthonous primary productivity, and the preservation therof. In proglacial environments in the Arctic, low temperatures restrict vegetation and soil cover, and minerogenic sediment input to lakes from glacial meltwater results in turbidity that impedes autochthonous productivity. Techniques for analyzing sediment from proglacial lakes therefore focus on investi-20 gating changes in the character of minerogenic sediment. The minerogenic content of lake sediments, used as a proxy for glacier size, is commonly evaluated by measuring the magnetic susceptibility and organic matter content of the sediments. Magnetic susceptibility reflects the amount of magnetic minerals eroded and input to a lake, and organic matter content is a function of dilution by minerogenic input, changes in primary productivity, and preservation. At Kulusuk Lake, processes that can complicate the mechanistic link between minerogenic input and glacier size are fundamentally limited. Input of sediment from 2017 Introduction non-glacial processes is restricted due to the small catchment and small catchment to lake-area ratio (∼ 2 : 1), and the proximity of the glaciers to the lake. These factors also limit the potential for sediment storage between the glaciers and the lake. Furthermore, the landscape surrounding the lake is composed of shallow, low-elevation slopes that minimize the likelihood of mass wasting events (which, if present, would have been 5 easily identified by rapid changes in sediment accumulation, or age reversals, based on 14 C dates). Therefore, at Kulusuk Lake, it is reasonable to interpret changes in minerogenic content as a function of glacier size. Variations in magnetic susceptibility, organic matter content, and scanning XRF data (PC1) represent changes in the relative amount and grain size of minerogenic sediment 10 delivered to Kulusuk over the last 9.5 ka (Fig. 3). Magnetic susceptibility and PC1 are directly related to, and organic content is inversely proportional to, minerogenic content. The highest magnetic susceptibility and PC1 values correspond to intervals with coarser grain size. We therefore interpret the sedimentological system in Kulusuk Lake as follows: during periods of increased glacier size, more coarse minerogenic sediment 15 was eroded from the bedrock and delivered to the lake by meltwater; during periods of smaller glacier size, less minerogenic sediment was deposited and a greater relative proportion of organic matter content accumulated.

Holocene glacier fluctuations in southeast Greenland
Dramatic changes in minerogenic input to Kulusuk Lake over the last 9.5 ka reveal 20 that the size of the Kulusuk glaciers has varied significantly throughout the Holocene (Fig. 3). Beginning ca. 8.7 ka, increasing organic matter content and decreasing minerogenic content, inferred from magnetic susceptibility and XRF data, indicates significant retreat of the Kulusuk glaciers, corresponding closely in time with the deglaciation of a nearby inland catchment ca. 8.4 ka (Balascio et al., 2013), following deglacia-25 tion of local coastal areas ca. 11.1-9.5 ka (Long et al., 2008;Roberts et al., 2008). A brief interval of increased minerogenic input shows that this early Holocene retreat was interrupted by an episode of advance at 8.5 ka, coeval with reductions in sea sur- face temperatures and bottom water circulation in the subpolar North Atlantic, as seen in high-resolution marine records (Ellison et al., 2006;Kissel et al., 2013). Another abrupt episode of minerogenic input ca. 8.2 ka signifies another glacier advance and occurred contemporaneously with the largest abrupt Holocene climate cooling event inferred from Greenland ice core records (Thomas et al., 2007). The temporal reso-5 lution and age control of this section of our record cannot provide new constraints on the exact timing of these events, however it clearly demonstrates the sensitivity of the Kulusuk glaciers to rapid, regional climate events. Between 7.8-4.1 ka, the Kulusuk glaciers were at their minimum Holocene extent, inferred from low minerogenic content, low MAR, and high organic matter content in 10 the lake sediments (Fig. 3). We interpret this as an interval with little to no glacier ice in the catchment, primarily based the XRF and magnetic susceptibility data, which are lowest and show reduced variability at this time, relative to the last 4.1 ka. This interval is also marked by extremely high organic matter content that is greater than 12 % (with a maximum of 19 %), suggesting that this period was accompanied by an increase in 15 primary productivity due to a lack of input of glacial flour. Magnetic susceptibility remains close to zero throughout the mid-Holocene section of the core and and appears insensitive to the minor minerogenic changes inferred from the XRF data. There are two excursions in the PC1 record during this interval (ca. 7.2 and 6.2 ka), which we interpret as sediment influxes from paraglacial activity rather than as glacier advances 20 because they are short-lived and do not match the amplitude of variation observed elsewhere. Thus, this record provides well-dated constraints on the Holocene Thermal Maximum in this area and refines previous estimates for its onset and termination (Kaufman et al., 2004).
At 4.1 ka, a sharp increase in XRF-and MS-inferred minerogenic content and de-25 crease in organic matter content, indicate the glaciers once again grew large enough to contribute minerogenic material to the lake. The regrowth of the Kulusuk glaciers represents the lowering of the regional snowline, and the precise timing could be considered unique to this catchment. However, the timing is contemporaneous with hydro-Introduction logic changes at nearby Flower Valley Lake, likely related to an increase in the duration of lake ice cover (Balascio et al., 2013). We propose that this represents significant cooling and the onset of the regional Neoglacial period. The oscillatory and stepwise increase in minerogenic input (decrease in organic matter content) after 4.1 ka suggests that rather than advancing steadily toward their historical extent, the Kulusuk 5 glaciers episodically advanced and retreated at centennial timescales until ca. 1.3 ka. After ca. 1.3 ka they stabilized and advanced only gradually until their rapid 20th century retreat (Fig. 3). Importantly, the major sedimentological transitions in the record are all located near radiocarbon dates, thereby maximizing the certainty of their timing and the calculations of sediment accumulation rates. The timing of glacier size variations 10 between radiocarbon-dated intervals since 4.1 ka are interpolated, and we estimate the accuracy to be better than ±100 years, the average 2-σ uncertainty of the ages.

Evidence for synchronous regional glacier response during the late Holocene
The Kulusuk glacier reconstruction documents centennial scale episodes of glacier ad- 15 vance during the Neoglacial (4.1 to 1.3 ka) coeval with other records of glacier growth in the North Atlantic region. After 4.1 ka, six major advances of the Kulusuk glaciers occurred (4.1, 3.9, 3.2, 2.8, 2.1, and 1.3 ka) and each successive advance resulted in greater glacier extent (Fig. 4). The progressive increase in glacier size is consistent with declining NH summer insolation. However, each episode of glacier advance was 20 followed by a period of retreat (or at least stabilization), possible suggesting that the glaciers repeatedly grew out of equilibrium with external insolation forcing and then retreated back toward an equilibrium state. The episodic advances of the Kulusuk glaciers during the past 4.1 ka are nearly identical in timing to the cooling episodes in the North Atlantic Ocean inferred from ice-rafted hematite-stained grains identified 25 in marine sediment cores (Bond et al., 1997(Bond et al., , 2001 (Fig. 4). Cooling events at these times have also been documented on the East Greenland and Icelandic Shelves and attributed to increased strength of the East Greenland Current (Giraudeau et al., 2000;2020 Jennings et al., 2002;Ran et al., 2008). Moreover, the Langjökull ice cap in Iceland advanced along with the Kulusuk glaciers and the North Atlantic IRD events   (Fig. 4), and advances of the Bregne ice cap in east Greenland at ca. 2.6 and 1.9 ka (Levy et al., 2014), within chronological uncertainty of the Kulusuk glacier advances ca. 2.8 and 2.1 ka, have also been documented. We propose that continuous 5 records of glacier activity around the North Atlantic during the Neoglacial are beginning to reveal synchronous glacier response to abrupt episodes of climate change.
Although the amplitude of variability in the proxy measurements during the past 1.3 ka are lower than earlier in the Holocene, due to the greater size and stability of the Kulusuk glaciers, it is still worthwhile to examine the changes in the sediment properties during this time. The very high sediment accumulation rates during this interval (0.8 mm year −1 ), allow sub-annual XRF measurements and, if interpreted in the same manner as periods with smaller glacier size, can afford a detailed examination of changes in glacier size using the XRF PC1 data (Fig. 5). The overall trend reveals a small and very gradual glacier expansion after 1.3 ka followed by 20th century retreat, 15 which resembles the overall trend in Arctic temperatures over the last 2 ka (Kaufman et al., 2009).
Multi-decadal variations in inferred glacier size during the past 1.3 ka also appear synchronous with those of other glaciers in the region (Fig. 5). Kulusuk glaciers increased in size ca. AD 1250-1300 and again ca. AD 1350 and AD 1450, precisely 20 when ice caps on Baffin Island  and Iceland (Larsen et al., 2011) were expanding (Fig. 4). After AD 1450 Kulusuk glaciers continued to expand, as did Langjökull on Iceland, while evidence from the Baffin ice caps indicates continuous ice cover  consistent with records from around Greenland that suggest the most extensive glacier advances since the early Holocene occurred between AD 1250-1900, and provide evidence for regionally coherent cooling phases during the Little Ice Age (Grove, 2001). We note that this timing contrasts with evidence from east Greenland that suggests LIA advance of the Istorvet Ice Cap approximately 100 years earlier (ca. AD 1150;Lowell et al., 2013), unless the data are reinterpreted as suggested by Miller et al. (2013). Therefore, we propose that there has been regional coherence in glacier activity not only during the past 1.3 ka, as previously suggested , but also during the past 4.1 ka, and glacier growth in response episodic climate change has been a common feature of glaciers in the North Atlantic region throughout, at least, the 10 last 4.1 ka.
Cold events are an important feature of centennial-scale climate of the Holocene (Wanner et al., 2011). Cooling events in the North Atlantic region are likely directly associated with changes in Atlantic Meridional Overturning Circulation (AMOC) (Denton and Broecker, 2008). IRD data suggest that periodic circulation changes of the North 15 Atlantic Ocean resulted in an advection of cold, fresh surface water south and east during ice-rafting events throughout the Holocene (Bond et al., 1997). Ocean circulation changes associated with IRD events and sea-surface temperatures have been attributed to solar forcing (Bond et al., 2001;Moffa-Sánchez et al., 2014;Jiang et al., 2015) and some modeling studies have confirmed that AMOC can switch between dis-20 tinct modes in response to a small external forcing, such as solar variability (Jongma et al., 2007). However, modeling results are inconsistent and it is also possible that cooling events might have simply resulted from internal ocean dynamics (Schulz and Paul, 2002). Regardless of the mechanism, our results demonstrate that glaciers responded quite actively to natural climate variations of the Holocene.

Rates of glacier change during the Holocene
This well-dated, high-resolution record of changes in the size of the Kulusuk glaciers also allows comparison among the rates of past glacier size variations. We present 2022 Introduction relative rates of change inferred from the first derivative of the XRF PC1 data in 105 year binned averages, an interval chosen using the interval with the lowest resolution (Fig. 6). We acknowledge the caveat that they are based on the assumption that the relationship between minerogenic input and glacier size has remained constant. The analysis indicates that the rate of 20th century retreat of the Kulusuk glaciers was 5 greater than during any other century of the past 1.3 ka, including during the Medieval Climate Anomaly. Furthermore, the 20th century retreat rate was two to three times the rate of any other period of retreat during the past 4.1 ka, and almost twice as rapid as the early Holocene retreat that marked the transition into the regional HTM (Fig. 6). This comparison helps to place the rate of 20th century glacier loss in the context of 10 natural episodes of past glacier activity.

Conclusions
The Kulusuk Lake sediment record was used to generate a high-resolution record of changes in the size of the Kulusuk glaciers over the last 9.5 ka. The record shows that the glaciers were sensitive to a number of previously documented regional climate fluc- 15 tuations and extends our understanding of Holocene climate dynamics in this sector of the Arctic. In particular, the record clearly constrains the Holocene Thermal Maximum at this site to between 7.8 and 4.1 ka, when the glaciers likely completed melted away. The regrowth of the Kulusuk glaciers at 4.1 ka corresponds with regional hydrologic changes and reflects the onset of the Neoglacial Period. The last 4.1 ka is marked 20 by a series of abrupt glacier advances as the size of the Kulusuk glaciers increased. These episodes of glacier growth correspond with ice rafting events in the North Atlantic Ocean, as well as regional ice cap expansion, and demonstrate that glaciers in this sector of the Arctic were very active during the late Holocene in response to abrupt cooling events that punctuated millennial-scale insolation-driven cooling.