Glacier response to the Little Ice Age during the Neoglacial cooling in Greenland

ice cores, elevation changes of the ice sheet across Greenland and oceanographic re-constructions from marine sediment cores over the past 5,000 years. The data are derived from a comprehensive review of the literature supplemented with unpublished reports. Our review provides a synthesis of the sensitivity of the Greenland ice margins and their variability, which is critical to understanding how Neoglacial glacier activity was interrupted by the current anthropogenic warming. We have reconstructed three distinct periods of glacier expansion from our compilation: two older Neoglacial advances at 2,500 – 1,700 yrs. BP (Before Present = 1950 CE, Common Era) and 1,250 – 950 yrs. BP; followed by a general advance during the younger Neoglacial between 700-50 yrs. BP, which represents the LIA. There is still insufficient data to outline the detailed spatio-temporal relationships between these periods of glacier expansion. Many glaciers advanced early in the Neo-glacial and persisted in close proximity to their present-day position until the end of the LIA. Thus, the LIA response to Northern Hemisphere cooling must be seen within the wider context of the entire Neoglacial period of the past 5,000 years. Ice expansion appears to be closely linked to changes in ice sheet elevation, accumulation, and temperature as well as surface-water cooling in the surrounding oceans. At least for the two youngest Neoglacial advances, volcanic forcing triggering a sea-ice /ocean feedback, could explain their initiation. There are probably several LIA glacier fluctuations since the first culmination close to 1250 CE (Common Era) and available data suggests ice culminations in the 1400s, early to mid-1700s and early to mid-1800s CE. The last LIA maxima lasted until the present deglaciation commenced around 50 yrs. BP (1900 CE). The constraints provided here on the timing and magnitude of LIA glacier fluctuations delivers a more realistic background validation for modelling future ice sheet stability.

In the Northern Hemisphere, an insolation driven Early to Middle Holocene Thermal Maximum was followed by a Neoglacial cooling that culminated during the Little Ice Age (LIA).Here, we review the glacier response to this Neoglacial cooling in Greenland.Changes in the ice margins of outlet glaciers from the Greenland Ice Sheet as well as local glaciers and ice caps are synthesized Greenland-wide.In addition, we compare temperature reconstructions from ice cores, elevation changes of the ice sheet across Greenland and oceanographic reconstructions from marine sediment cores over the past 5,000 years.The data are derived from a comprehensive review of the literature supplemented with unpublished reports.Our review provides a synthesis of the sensitivity of the Greenland ice margins and their variability, which is critical to understanding how Neoglacial glacier activity was interrupted by the current anthropogenic warming.We have reconstructed three distinct periods of glacier expansion from our compilation: two older Neoglacial advances at 2,500 -1,700 yrs.BP (Before Present = 1950 CE, Common Era) and 1,250 -950 yrs.BP; followed by a general advance during the younger Neoglacial between 700-50 yrs.BP, which represents the LIA.There is still insufficient data to outline the detailed spatiotemporal relationships between these periods of glacier expansion.Many glaciers advanced early in the Neoglacial and persisted in close proximity to their present-day position until the end of the LIA.Thus, the LIA response to Northern Hemisphere cooling must be seen within the wider context of the entire Neoglacial period of the past 5,000 years.Ice expansion appears to be closely linked to changes in ice sheet elevation, accumulation, and temperature as well as surface-water cooling in the surrounding oceans.At least for the two youngest Neoglacial advances, volcanic forcing triggering a sea-ice /ocean feedback, could explain their initiation.There are probably several LIA glacier fluctuations since the first culmination close to 1250 CE (Common Era) and available data suggests ice culminations in the 1400s, early to mid-1700s and early to mid-1800s CE.The last LIA maxima lasted until the present deglaciation commenced around 50 yrs.BP (1900 CE).The constraints provided here on the timing and magnitude of LIA glacier fluctuations delivers a more realistic background validation for modelling future ice sheet stability.

Introduction
Ongoing Arctic warming has prompted a huge interest in how the Greenland Ice Sheet (GrIS) responded to climate change during the last millennia, evidence which seems crucial for understanding the ice sheet behavior to long-term changes in external forcing.At the moment, predictions on ice sheet stability are hampered by the lack of knowledge about the potentially delayed impact of the LIA cooling and major ice advance on modern ice mass loss.Here, we attempt to synthesize existing data into a compilation that provides a proper baseline assessment.
The Holocene Thermal Maximum (HTM) was a period with Northern Hemisphere temperatures higher than present (Kaufman et al., 2004;Axford et al., 2021).In Greenland, depending on the location and type of proxies it has been dated to between c. 8,000 and 4,000 yrs.BP (Funder et al., 2011;Weidick et al., 2012;Fréchette and de Vernal, 2009;D'Andrea et al., 2011;Axford et al., 2013;Bennike et al., 2010;Kobashi et al., 2017;Buizert et al., 2018;Andresen et al., 2022).The termination of the HTM was succeeded by a cooler Neoglacial period, which has been conventionally subdivided into older and younger phases (Weidick et al., 2012).We adapt to this, but assign the older Neoglacial phase from 5,000 to 950 yrs.BP and the younger known as the LIA from 750 to 50 yrs.BP; these cold periods are intersected by a regional atmospheric warming during the Medieval Warm Period (MWP) c. 950 -750 yrs.BP.
Climatically, the LIA is a period associated with reduced global temperatures and widespread glacier advances (Grove, 2004).Other Fig. 1.Greenland ice sheet in the glaciated northern hemisphere.Grey dots: deep ice cores in Greenland and Ellesmere Island.Seven zones subdivide the ice sheet into main mass balance drainage areas (modified from Rignot and Mouginot, 2012).N (north), NE (northeast), CE (central east), SE (Southeast), SW (southwest), CW (central west), NW (north west).Arrows indicate conceptual patterns of wind driven transport of aerosols and drift ice associated with the large-scale atmospheric patterns.Asian mineral dust aerosol is transported eastward to Greenland in spring, and drift ice is transported southward during summer and fall (purple arrows) associated with high pressure over Greenland, while sea salt aerosol is transported northward to Greenland from the North Atlantic during winter (grey arrow) associated with atmospheric low pressure around Iceland.During the Holocene cold periods, these wind patterns have been intensified, as indicated by increased loadings of mineral dust and sea salt aerosols in Greenland and by more southerly migrations of drift ice.The intensification is combined with a generally southward shifted westerly circulation.associated historical impacts are the North Atlantic Fish Revolution (c.1500 CE) (Holm et al., 2019), the end of the Greenlandic Norse settlements (c.1350 CE) (Dugmore et al., 2012), and expanded sea-ice cover in the Arctic Ocean (c.1700 -1800 CE) (Macias Fauria et al., 2010).Northern Hemisphere temperature reconstructions estimate a LIA cooling between 0.5 and 1.0 • C relative to the global mean (Mann et al., 2008(Mann et al., , 2009;;Ahmed et al., 2013;Kaufman et al., 2009;Briner et al., 2016;Kobashi et al., 2017;Neukom et al., 2019a).Reconstruction of global mean temperature suggests that the early 20th century (1900)(1901)(1902)(1903)(1904)(1905)(1906)(1907)(1908)(1909) was cooler than >95% of the Holocene (Marcott et al., 2013).Instrumental temperature records and model simulations also show significantly warmer conditions after the 1920s, thus effectively terminating the LIA and thus the Neoglacial period (Box et al., 2013;Neukom et al., 2019b).
In Greenland, glaciers also responded to the Neoglacial cooling, but it has generally been assumed that the time when glaciers reached their maximum extent was somewhat delayed relative to the onset of the LIA temperature cooling with maximum glacier advance in the 1700s CE or the end of the 1800s CE (Dahl-Jensen et al., 1998;Weidick et al., 2012) (Fig. 1).In recent years a number of studies have been conducted in Greenland to understand the ice marginal response of the GrIS to the Neoglacial cooling e.g., Young and Briner (2015); Briner et al. (2016) and Larsen et al. (2018), as well as local glaciers and ice caps (Bradley and Bakke, 2019).Furthermore, several ice-, lake-and marine-proxy records have documented the climate variations leading up to the culmination of the LIA (Briner et al., 2016;Axford et al., 2021).However, until now no Greenland-wide assessment of the climate and glacier response for this specific period has been undertaken (Vasskog et al., 2015).The following synthesis falls into three parts.The first part describes the dynamics of the ice sheet and surrounding ocean variability using records from ice cores as well as lacustrine and marine sediment cores.The second part, compiles ice margin changes from a series of point-specific proxy records.The last attempts to reconcile all the proxies considered during the ice sheet advance, culmination, and retreat phases throughout the Neoglacial.
Data compilation is centered around the lastest advanceculmination-retreat phases of glacier ice in Greenland and the associated terrestrial geomorphological and stratigraphic records.For advance phases, we used 14C dates organic material (peat, gyttja, wood, plant microfossils), bone (reindeer antlers, polar bear skulls) and predominantly shells (molluscs and mollusc fragments).Dates from mixed shell samples were excluded due to age variability.In the stratigraphical record, we aimed for LIA sequences covering the last 2000 yrs in threshold lakes.Subfossil plants were compiled from microfossils that recently melted free from retreating glacier fronts.Ice cores rely on the ice-core chronology established for individual sites.Culmination phases rely on cosmogenic exposure ages (10Be) from end moraines within the Neoglacial landscape. 10Be ages follows the original publication and the production rates used.Absolut dates from the retreat phases comes from minimum limiting 14C ages from organic matter above minerogenic sediment in threshold lakes while other estimates are derived from lichen ages and historical maps and aerials.Common for all compiled dates is a census date of 1 February 2020.
The two-part division of the Neoglacial serves as a framework for comparison between the different temporal datasets and as the review proceeds below should become self-evident.We note that the LIA is not assigned any chronostratigraphic time and therefore considered an informal stratigraphic unit similar to the MWP with undefined boundaries and age ranges.This is also clearly reflected in the onset times suggested in the reviewed literature for the LIA cooling spanning between 800 -550 yrs.BP.In this review, we operate with a range of age scales from the Middle Holocene to the 20th century, which constitute a challenge when reporting age within the same age nomenclature.We consequently report all ages in years before present (yrs.BP), but supplement reported ages in the Common Era (CE) when it makes sense in a historical context.

Climate records in Greenland
In the following we briefly review the ice-core records focusing on Holocene temperature reconstructions and atmospheric circulation changes during the Neoglacial cold periods.

Ice-core temperature reconstruction
The composition of stable water isotopes in glacier ice, δ 18 O and δ 2 H, reflects cloud temperatures during past precipitation (Dansgaard et al., 1973), and has been widely used as a paleo-thermometer for air temperatures above the snow surface.The ice-core stable isotope records from GRIP, NGRIP, Dye 3, Camp Century, Renland, and GISP2 in Greenland (Johnsen et al., 2001 and references therein) and from Agassiz, Ellesmere Island in Canada (Fisher et al., 1995) have provided continuous, high-resolution proxy records of Holocene temperature variations.These records were corrected for effects from altitude changes due to ice flow and ice-sheet elevation changes (Vinther et al., 2009;Lecavalier et al., 2017).However, the paleo-thermometer deduced from stable water isotopes may be biased because it reflects the temperature only at times of precipitation, it is influenced by temperature changes at the moisture source and changes in the transport path-way, and the relationship between isotope composition and air temperature is complicated by local factors such as temperature inversion strength (e.g., Jouzel et al., 1997).Therefore, other temperature proxies from ice cores have been used along with the isotopic records, specifically melt layers observed in the ice cores which indicate that temperatures had been high enough to produce melt at the snow surface.The quantity of melt layers has been combined with δ 18 O records to constrain summer temperatures (Alley and Anandakrishnan, 1995;Fisher et al., 1995;Fisher et al., 2011;Herron et al., 1995;Kameda et al., 1995).Snow accumulation is inversely correlated with temperature and can, in similarity with stable isotope ratios, be used as a temperature proxy (Meese et al., 1994).Meese et al. (1994) used a series of annual markers in the GISP2 ice core to construct a snow accumulation record spanning the Holocene period.Gravity-driven fractionations of Ar and N 2 gas isotopes in the snow column respond to vertical temperature gradients produced by changes in both the temperature at the snow surface and to the accumulation rate.Used in combination, δ 40 Ar and δ 15 N from enclosed air bubbles in glacier ice provide a paleothermometer that is independent of the conditions for precipitation (Severinghaus et al., 1998) and the GISP2 ice core records of δ 40 Ar and δ 15 N have provided a Holocene temperature reconstruction for Greenland (Kobashi et al., 2017).Buizert et al. (2018) reconstructed Holocene temperatures using an inverse modelling of the GISP2 δ 15 N record.The ice temperatures in glaciers are themselves archives of past temperatures because thermal conductivity is poor in ice.Therefore, Monte Carlo inverted borehole temperature records have been read as accurate, if poorly resolved, temperatures from the past (Dahl- Jensen et al., 1998), and have served as calibrations of the better resolved temperature proxies from ice cores (Fig. 2).
The temperature records based on δ 18 O from ice cores show a HTM followed by a gradual cooling trend culminating with the LIA (Fig. 2).The timings of HTM and thus the onset of the Neoglacial cooling, varies between the different δ 18 O ice core records with the North GRIP record (Andersen et al., 2004a(Andersen et al., , 2004b;;Vinther et al., 2006;Johnsen et al., 2001) dating the HTM to 8,000 -4,000 yrs.BP, while the records from the Renland and Agassiz ice caps (Vinther et al., 2009) yield dates on HTM to 10,000 -6,000 yrs.BP (Buizert et al., 2018).A more recent δ 18 O temperature reconstruction from the Agassiz ice cap dates the HTM even earlier, to 11,000 -8,000 yrs.BP with maximum temperatures around 2 • C warmer than present (i.e., 2009 CE), indicating a temperature gradient between the Renland and Agassiz ice core locations that is reducing during the Early Holocene (Lecavalier et al., 2017).Pronounced centennial and millennial scale variations are absent in the δ 18 O records except for a very pronounced event at 8,200 yrs.BP (Alley et al., 1997;von Grafenstein et al., 1998;Barber et al., 1999;Alley and Ágústsdóttir, 2005).The discrepancies between the ice-core records with respect to timing of HTM are not only due to latitudinal differences in the Holocene temperature development.The records are also affected by limitations in the temperature reconstructions from water isotopes.For the same reason, there are only limited millennium-scale variations in the temperature records i.e., no indication of a cooling that could identify the onset of a Neoglacial period after the HTM.However, such a cooling can be seen in other Greenland ice-core temperature reconstructions.
The records of ice melt from GISP2 (Alley and Anandakrishnan, 1995) and Agassiz (Lecavalier et al., 2017) show pronounced melting during the HTM that decreased markedly around 6,000 BP indicating a decrease in summer temperatures at that time.The gas fractionation temperature reconstruction from GISP2 (Kobashi et al., 2017) show, in contrast to the δ 18 O records, large centennial to millennial scale climate variability with a peak warmth around 8,000 yrs.BP.This is followed by a general cooling trend with a major cold period during 6,100 -5,000 yrs.BP, which is again followed by a slight temperature increase until around 4,500 yrs.BP (Mid-Holocene Optimum) (Fig. 2).Afterwards, the temperature shows large variability with a slow underlying cooling trend toward the temperature minimum during the LIA -600 yrs.BP (1350 CE).The MWP is indicated in the interval from 1,350 -600 yrs.BP (600 -1350 CE) with a peak warmth around 1,200 yrs.BP (800 CE).The maximum temperatures around 7,860 yrs.BP were 2.9 ± 1.4 • C warmer than recent decades, while the decades 210 -170 yrs.BP (1740 -1780 CE) were 2.9 ± 0.9 • C colder.
The Monte Carlo inverted borehole temperature record from GRIP (Dahl- Jensen et al., 1998) dates the HTM to 8,000 -5,000 yrs.BP, while the borehole temperature reconstruction from Dye 3 shows a much younger HTM 4,000 -5,000 yrs.BP.These reconstructions exhibit decreased temporal resolution back in time, and millennial scale temperature changes in the early Holocene cannot be detected (Dahl-Jensen et al., 1998) and so they may not be the preferred reconstructions of HTM.However, the two borehole records show a similar temperature history for the last 4,000 yrs.with a temperature minimum around 2,000 yrs BP and a temperature maximum around 950 yrs.BP (1000 years CE) (MWP) 1 • C warmer than present in Greenland.In the records, the LIA is made up of two cold periods, with temperature minima at 400 and 100 yrs.BP (1550 and 1850 CE), and with temperatures respectively 0.5 and 0.7 • C below the present.After the LIA, temperatures reach a maximum around 20 yrs.BP (1930 CE).The findings of two cold periods during the LIA are also recognized in the snow accumulation and meltlayer records from GISP2 (i.e., Meese et al., 1994) where two cold phases can be observed at around 400 and 100 yrs.BP (1550 and 1850 CE).In a composite snow accumulation record retrieved from the seasonal isotope data from GRIP, NGRIP and Dye 3, however, occurrence of low accumulation is found but not consistent with cooling periods around 400 and 100 yrs.BP (1550 and 1850 CE) (Andersen et al., 2006).
The δ 18 O signal in the Greenland ice cores can be resolved into annual cycles for the last 2,000 yrs.(Johnsen et al., 2000) and the seasonal δ 18 O values have been shown to provide a better temperature proxy than annual mean δ 18 O values (Vinther et al., 2010;Zheng et al., 2018).In line with the Monte Carlo inverted borehole temperatures from GRIP and Dye 3 (Dahl-Jensen et al., 1998) the winter δ 18 O signal from GRIP, Crete and Dye 3 all show a MWP from 1,150 -950 yrs.BP (800 -1000 CE) and a cold period around 90 yrs.BP (1860 CE) (Vinther, 2011).However, a cold period around 1500 CE that is clearly seen in the borehole temperature records and, in the Dye 3 winter δ 18 O signal is not seen in the winter δ 18 O signals from GRIP and Crete (Vinther, 2011).

Atmospheric circulation across Greenland in Neoglacial time
Records of soluble ions in the GISP2 ice core provide proxy records of atmospheric circulation strength over the Holocene period (O'Brien et al., 1995).More specifically nssK + (Fig. 3A) that records deposition of windblown terrestrial material from the Asian continent during the spring season is thought to be a proxy for the Siberian high, and ssNa + (Fig. 3B) that records depositions of sea salt aerosols from the North Atlantic during the winter season is a proxy for the Icelandic low (Meeker and Mayewski, 2002).
In contrast to the Greenland δ 18 O records, the time series of aerosol species from GISP2 points to quasi-cyclic ~2600-year interval of colder climate throughout the Holocene among which the 8,200 yrs.BP event and the LIA are found (O'Brien et al., 1995).The cold intervals are especially clear in the nssK + fluxes where the dates given are, 8,800 -  (Vinther et al., 2006), (blue) gas-fractionation temperature from GISP2 (Kobashi et al., 2017), (B) Ice Sheet elevation and upstream corrected temperature from isotope records, red: Renland and Agassiz (Vinther et al., 2009), blue: Agassiz (Lecavalier et al., 2017), (C) Borehole temperatures from GRIP (blue) and Dye 3 (green) (Dahl-Jensen et al., 1998).We subdivide the Neoglacial into an older and younger part similar to Weidick et al. (2012), but assign an age to the older part between 5,000-950 yr BP. and the younger part to 750-50 yr BP.
7,800; 6,100 -5,000; 3,100 -2,400; and 600 -0 yrs.BP that coincide with glacier advances and cooling intervals observed worldwide (O'Brien et al., 1995;Mayewski et al., 2004;Rohling et al., 2002).The transition from the MWP into the LIA around 600 yrs.BP is the most rapid transition observed in the record of the Holocene period (Meeker and Mayewski, 2002).O'Brien et al. (1995), explain increased fluxes of Asian dust and sea salt in Greenland by more arid conditions over continents, a southward displacement of the westerlies resulting in increased atmospheric loadings of Asian dust, and increased meridional atmospheric transport of sea salt aerosol in the North Atlantic.Meeker and Mayewski (2002) compared nssK + and ssNa + from the top of the GISP2 ice core with monthly instrumental sea level pressure data.They found that the LIA was associated with an enhanced high-pressure anomaly over Siberia in spring, and an enhanced low-pressure anomaly over south Greenland and the northern North Atlantic during winter and a high-pressure anomaly in northern Canada in both spring and winter.
The cold intervals found in the GISP2 glacio-chemical series have been linked with the "Bond cycles" that are ~1500-yr intervals of drift ice in subpolar North Atlantic areas seen in records of ice rafted debris (IRD) (Fig. 3C) and pointing to a strong ocean-atmosphere coupling during these intervals (Bond, 1997;Wanner et al., 2011).The millennial-scale occurrence of drift ice in subpolar areas is explained by a southward and eastward displacement of the East Greenland Current (EGC) and Labrador Currents (LC), which carry sea ice-laden cold and fresh water southward along the East coast of Greenland and Canada, respectively (Bond, 1997;Bond et al., 2001).These changes in the North Atlantic Ocean circulation were associated with a basin-wide cooling that may have resulted from reduced North Atlantic Deep-Water formation and was accompanied by enhanced northerly winds in the North Atlantic associated with an intensified high-pressure anomaly over Greenland (Bond, 1997, Bond et al., 2001).A strong ocean-atmosphere coupling in the North Atlantic during the Holocene cold events have also been observed in other atmospheric circulation proxies elsewhere.In monsoon proxy records from the Arabian Sea (Gupta et al., 2003) and Lake Qinghai on the Northeastern Tibetan Plateau (An et al., 2012), as well as in stalagmites from the Dongge Cave in south China (Wang et al., 2005), a link between weakened Asian summer monsoon and enhanced IRD in the North Atlantic has been found.Weakened summer monsoons result from a southward displacement of the intertropical convergence zone and the westerlies during summer.Proxy data from Lake Qinghai further shows that strengthened westerly circulation in winter is linked with enhanced IRD in the North Atlantic (An et al., 2012).Also, Sorrel et al. (2012) and Goslin et al. (2018) identified periods of intensified storminess in Europe throughout the Holocene, indicating a strengthening and southward displacement of westerly winds, which was linked with enhanced IRD in the North Atlantic.
A major feature of atmospheric circulation changes across Greenland during the Holocene cold events could be a generally southwarddisplacement and intensification of westerly circulation.During winter, this pattern produced both stronger meridional flow in the North Atlantic region and an intensified Icelandic low.The spring/summer were comparatively dry over continents due to enhanced zonal atmospheric transport.Northerly winds in the North Atlantic flowed from a high-pressure anomaly over Greenland.Such changes in atmospheric circulation must have had implications for precipitation in Greenland, and changes in snow accumulation associated with the LIA cooling have been observed in some ice cores e.g., from GISP2 (Meese et al., 1994) and from NEEM (Masson-Delmotte et al., 2015).However, on hundredyear averages, large spatial differences in accumulation rates east and west of the ice divide in central North Greenland was observed in airborne radar-detected stratigraphy (Karlsson et al., 2016).Over the last 700 yrs., accumulation rates varied little in the area west of NEEM, while accumulation rates varied up to 20% in the area north of NGRIP.These differences suggest that variations in snow accumulation observed in ice cores were local changes, in response to changed circulation, rather than large scale changes in snow accumulation (Karlsson et al., 2016).

Lake sediment records
Relatively few quantitative reconstructions of climatic variables have been produced from climate-sensitive lakes around the margins of the GrIS.The longest and most sensitive sedimentary records are typically found in headwater lakes located above the postglacial marine limit, and that did not receive meltwater inputs from glaciers following Early Holocene deglaciation i.e., non-glacial fed lakes (see .Traditionally, the Holocene development of environments and temperatures were based on pollen analyses from lake sediment records, and a survey of pollen diagrams in all parts of Greenland showed that a reduction in vegetation in most areas began at 4,500 yrs.BP, and as early as 6,000 yrs.BP in East Greenland (Funder and Fredskild, 1989).More recently, quantitative reconstructions of July or summer air temperature and summer lake temperature have been developed using training sets of chironomids (Axford et al., 2013(Axford et al., , 2017(Axford et al., , 2019;;McFarlin et al., 2018), pollen (e.g., Fréchette and de Vernal, 2009;Gajewski, 2015), and algal alkenones (D'Andrea et al., 2011), respectively.Annual paleo-precipitation records and effective annual Age BP Fig. 3. A) and B) show respectively the nssK + and the ssNa + concentrations from the GISP2 ice core where intervals of high concentrations indicate colder climate with more intense atmospheric circulation (O'Brien et al., 1995).The two curves have been used as proxies for intensities of the spring Siberian high and winter Icelandic low atmospheric pressure systems (Meeker and Mayewski, 2002).C) shows the stacked IRD (percent hematite-stained grains) curve from the subpolar North Atlantic (Bond et al., 2001) where higher contents of IRD indicate a southward advection of drift-ice from EGC and LC.The yellow vertical bars are placed at periods of enhanced storminess in coastal northern Europe (Sorrel et al., 2012) that are indicating a southward displacement and intensification of the westerlies.Dashed lines bounds the older neoglacial and younger neoglacial periods.moisture were reconstructed using δ 2 H from leaf waxes (Balascio et al., 2013;McFarlin et al., 2019;Thomas et al., 2016Thomas et al., , 2020)), δ 18 O from chironomids (Lasher et al., 2017;Lasher and Axford, 2019), biogenic silica (Andresen et al., 2004), organic matter and mineral flux (Anderson et al., 2012) as well as diatoms (Perren et al., 2012).Mineral content has also been used to semi-quantitatively reconstruct summer air temperature (Willemse and Törnqvist, 1999).While these records provide a relevant context for Late Holocene terrestrial climate history, most have focused on the Early to Middle Holocene and are resolved on a millennial to century scale, with a single or only a few observation points covering the last ~1500 yrs.The Holocene climate change from many lake records in Greenland have been synthesized elsewhere (Sundqvist et al., 2014;Briner et al., 2016;Axford et al., 2021).Therefore, only a brief summary of the findings relevant to the Neoglacial and LIA periods is provided here along with an overview of more recent work showing synchronous regional changes.

South and Southwest Greenland
In the South and Southwest Greenland, multiple lines of evidence from biological and geochemical proxy data support a cooling trend from the Middle to Late Holocene by as much as 2 to 5 • C culminating in the last millennium (Briner et al., 2016 and references therein).The coldest summer conditions reconstructed from chironomid assemblages for South Greenland during the LIA were slightly warmer than elsewhere in Greenland i.e., ~2 • C colder than modern values (Lasher and Axford, 2019;Chipman et al., 2018) and occurred relatively late into the 18th and 19th centuries (250 -50 yrs.BP).In contrast, summer temperature reconstructions from pollen suggest warmer temperatures after ~4,500 yrs.BP and potentially spatially heterogeneous conditions across the regions (Briner et al., 2016;Fredskild, 1973Fredskild, , 1983a;;Gajewski, 2015;Kelly and Funder, 1974).More negative and variable values of the North Atlantic Oscillation (NAO) index after 650 yrs.BP were suggested to explain changes in hydro-climatic conditions during the LIA in South and Southwest Greenland (Olsen et al., 2012a;Lasher and Axford, 2019).
However, the influence of the NAO on temperature and precipitation during the Late Holocene in the region and further north remains debated (Aebly and Fritz, 2009;D'Andrea et al., 2011;Olsen et al., 2012a;Lasher and Axford, 2019;Lasher and Axford, 2019).Highly variable, yet overall warm conditions 1,050 -550 yrs.BP, thus including the early cold LIA period inferred from a sub-centennially resolved paleotemperature record from South Greenland could not be associated with a positive NAO phase.Rather, this warmer phase appeared consistent with shifts in subpolar ocean currents (Lasher and Axford, 2019).Shifts in vegetation assemblages around 2,000 yrs.BP suggested not only cooler, but also more humid conditions from the Middle Holocene comparable to current conditions in the Godthåbsfjord area of Southwest Greenland (Fredskild, 1983b).

West and Northwest Greenland
In Western and Northwestern Greenland, the Late Holocene is associated with overall cooler but more variable air temperatures and relatively higher effective precipitation compared to the Middle Holocene (McGowan et al., 2003;Anderson and Leng, 2004;Aebly and Fritz, 2009;D'Andrea et al., 2011;Anderson et al., 2012;Perren et al., 2012;Axford et al., 2013;Briner et al., 2016;Lasher et al., 2017).This apparent inverse relation between sediment organic carbon content (OC) and temperature prompts a more careful consideration of the mechanisms at play when inferring paleotemperatures from OC-based proxies in lacustrine records (Anderson et al., 2018).The alkenonebased, high-resolution records of lake temperature developed at Braya Sø and Lake E in West Greenland provide strong evidence for successive and abrupt changes to cold conditions during the Neoglacial (D'Andrea et al., 2011).The coldest summer lake water conditions in the last 1000 years are inferred between ~850 -630 yrs.BP and are reconstructed to be around 3-4 • C, which would be on the lower end of the measured summer range of 3-24 • C at present (D'Andrea et al., 2011).Farther north, the coldest summer and autumn water δ 18 O inferred temperatures of the Holocene record are ~1.5-3• C cooler than present and are observed after 1,200 yrs.BP (Lasher et al., 2017).Summer air temperature inferred from chironomid assemblages were ~4-5 • C lower than modern temperature between 250 -50 yrs.BP in West Greenland (Axford et al., 2013).This timing is in similarity with cooling in South Greenland, but contrasts with earlier cooling in Northwest Greenland i. e., around ~500 yrs.BP (McFarlin et al., 2018;Axford et al., 2019).The cooling was preceded by approximately 2 • C warmer than present conditions at North Lake West Greenland in the 20th century (50 yrs.BP onwards) (Axford et al., 2013).A recent multi-site analysis of catchment related organic carbon (OC) sources using carbon isotope composition, the C/N ratio and OC accumulation rates, showed that the Late Holocene cooling trend in West Greenland can be associated with an increase in terrestrial carbon burial in lake sediments and an overall decrease in authigenic production.This may be due at least partly to increased lateral transport of carbon to lakes with enhanced catchment instability (Anderson et al., 2018).
Some of the paleo-precipitation records in West Greenland suggest conditions in the Late Holocene were more variable than this broad picture.Analysis of δ 2 H from leaf waxes suggests that a change in moisture sources led to decreased winter snowfall from the Middle Holocene until present (Thomas et al., 2016).This decrease would be closely linked to the regional cooling trend in air and ocean surface waters through a weakening of the West Greenland Current and enhanced early winter sea ice cover in northern Baffin Bay and the Labrador Sea (Thomas et al., 2016).Lake water balance modelling from connected lakes at Kangerlussuaq revealed pluvial intervals at 4,600 and 2,000 yrs.BP overlaying overall relatively low lake levels for the Late Holocene, with a particularly low-stand in the last 700 years (~1,250 yrs.BP onwards), also pointing to changes in precipitations seasonality as a likely mechanism (Aebly and Fritz, 2009).Similarly, diatom records from the same area suggest the relatively wet Neoglacial conditions are followed by dryer conditions in the last 500 yrs.(~1,450 yrs.BP onwards; Perren et al., 2012).

North, East and Southeast Greenland
As summarized by Briner et al. (2016), most of the qualitative and semi-quantitative multi-proxy lake records including macro-and microfossil abundance, geochemical and physical data from North to Southeast Greenland also support a cooling trend after the Middle Holocene, some of which clearly culminate during the LIA, or at least from c. 1,800-900 yrs.BP and onwards (Wagner and Melles, 2002;Bennike and Wagner, 2012;Olsen et al., 2012a;Wagner and Bennike, 2015).This cooling trend has been associated with decreasing summer temperatures (Briner et al., 2016).The accelerated rate of Neoglacial summer cooling at Last Chance Lake (East Greenland) of 3-5 • C between peak HTM (~5,500 yrs.BP) and ~500 yrs.BP is large when compared to other types of paleotemperature records for the region (Vinther et al., 2009;Axford et al., 2017).Yet, the summer temperature values reconstructed are within the 4 -5 • C range also inferred in West and Northwest Greenland, centered earlier around 500 yrs.BP (Axford et al., 2017).Different chironomid training sets from modern assemblages in lakes located between the Canadian High Arctic Archipelago to Maine, eastern USA and Iceland have been applied to reconstruct July air temperatures in Greenland (Francis et al., 2006;Langdon et al., 2008), representing a large part of our current knowledge of past quantitative temperature change from lake records in the region.While these chironomid-based reconstructions are supported by good levels of fit to the respective models applied, high taxonomic resolution datasets are required to reliably detect the potentially high spatial and temporal variability of temperature changes during the LIA in Greenland beyond the ice sheet margins (Chipman et al., 2018;Axford et al., 2017;Medeiros et al., 2020).An Arctic-wide chironomid training set, incorporating data from a range of climatic settings, including from Southwest Greenland is currently in development (Brodersen and Anderson, 2002;Medeiros et al., 2020).This could improve model fitness for Greenland-based quantitative reconstructions.Recent work from Northeast Greenland combining macrofossils and biomarkers supports the Neoglacial cooling trend, but also shows only limited potential for branched glycerol dialkyl glycerol tetraethers (brGDGTs)-based air temperature reconstructions in this high latitude region (Kusch et al., 2019).
Quantitative precipitation reconstructions for East Greenland are sparse (Briner et al., 2016).Evidence from δ 2 H from leaf waxes from Southeast Greenland indicate a Late Holocene increase in effective moisture.The factors driving this increase could be longer periods of ice cover due to lower summer temperature which reduced lake water evaporation and increased runoff from higher winter snow accumulation (Balascio et al., 2013).Though paleo-precipitation records still remain too few to demonstrate such large-scale pattern and account for intra-regional variability, this evidence would suggest a precipitation seesaw across Greenland comparable with the observation of Thomas et al. (2016) of decreased winter snow accumulation during the Late Holocene in West Greenland.
The lake sediment records generally offer evidence for millennialscale climate change in different regions during the Late Holocene.A pattern characterized by cooler summers, wetter and generally more variable conditions than those prevailing during the mid-Holocene culminated between 500-50 yrs.BP, particularly in the South, Southwest, West and Northwest of Greenland where most records are found.
This colder climate, also identified from a compilation of annual temperature proxy records across the Arctic (McKay and Kaufman, 2014), was often followed by warmer conditions during the 20th century (50 yrs.BP onwards).Unsurprisingly, the more highly resolved records reflect higher spatial and temporal heterogeneity on decadal to century scales, suggesting terrestrial climate at the margins of the GrIS during the LIA and Late Holocene was probably less homogenous than this broad picture indicates.On a sub-centennial scale, modelled and proxy-based hydroclimatic reconstructions from ice cores and a limited subset of lake records in Greenland suggested wetter conditions occur during warmer periods over the last millennia (Charpentier Ljungqvist et al., 2016).This pattern, however, is opposite to that of the millennialscale hydroclimatic variability reconstructed from the Greenland lake records reviewed here, where the warm Early to Middle Holocene period has been associated with drier conditions relative to the Late Holocene.Some of this variability may be explained by nonlinear changes in climate seasonality throughout the Holocene.Improved spatial coverage of sub-centennial to decadal resolved lake proxy records sensitive to seasonal variability for the period will help refine our mechanistic understanding of large-scale (hydro)climatic dynamics at play.

Marine sediment records
Oceanographic conditions in East Greenland are highly influenced by the fresh and cold East Greenland Current (EGC) transporting sea ice and polar waters from the Arctic Ocean via the Fram Strait southwards along the East Greenland coast (Fig. 4).In East and Northeast Greenland, the subsurface Atlantic Intermediate Water (AIW) is linked with the North Atlantic drift as it originates from a southward turning branch of the warm West Spitsbergen Current and from southward flow of the Atlantic layer from the Arctic Ocean.Off the coast of Southeast Greenland, the warm and saline Irminger Current (IC) slides underneath the EGC at depths of c. 200 m.The IC and EGC mix off Southwest Greenland to become the West Greenland Current (WGC), which propagates northward.Off the coast of Northwest Greenland, the cold and fresh Baffin Bay Current flows southwards via Nares Strait, carrying sea ice (Buch, 2002).Marine geological studies from the Southeast Greenland shelf (Andersen et al., 2004a(Andersen et al., , 2004b;;Kuijpers et al., 2003), East Greenland (Knudsen et al., 2008a(Knudsen et al., , 2008b)), Denmark Strait (Andresen and Björck, 2005) and the North Atlantic (Koç et al., 1993) indicate a regional cooling following the HTM concurrent with a marked intensification of the EGC from around 4000-3500 yrs.BP and onwards to the LIA.This long-term trend, however, is characterized by centennial to multi-decadal scale fluctuations.Studies of marine sediment obtained in the vicinity of Greenland are summarized below (Fig. 4A-D).The summarized data covers the past c.2000 years to provide an overview of the conditions prior to the onset of the LIA.
The onset of the MWP (950 -1250 CE; Lamb, 1977) was rather asynchronous across Greenland and took place between 1,000 -700 yrs.BP.This period is also referred to as the Medieval Climate Anomaly (MCA), and is characterized by decreased sea-ice occurrence in North East Greenland (Spielhagen et al., 2011;Perner et al., 2015), Central East Greenland (Kolling et al., 2017;Andersen et al., 2004a, 2004b, Jennings and Weiner, 1996) and warming of the AIW (Perner et al., 2015(Perner et al., , 2016)).In contrast, a marine sediment core study offshore Southeast Greenland indicates that the EGC temperature was relatively constant between the DACP and the onset of the LIA with no signs of an intervening MWP warming (Andersen et al., 2004a(Andersen et al., , 2004b) ) and foraminiferal studies indicate glaciomarine conditions near Helheim Glacier between 1,450 and 650 yrs.BP (500 CE and 1300 CE) (Andresen et al., 2013) thus also with no traces of a warm MWP.In South Greenland most sites even record a cooling of the IC during the MWP, likely in response to the concurrent frequently positive NAO index (Dawson et al., 2002), which would strengthen the subpolar gyre and allow less IC water to spill onto the South Greenland shelf (Miettinen et al., 2015;Andresen et al., 2017).Although ice core studies indicate atmospheric warming over Greenland during the MWP (Kobashi et al., 2017) marine sediment core studies from Disko Bugt do not consistently track this warming but rather indicate that sea ice occurrence was variable here (Krawczyk et al., 2013;Ribeiro et al., 2012;Seidenkrantz et al., 2008;Moros et al., 2006) and that the WGC temperature fluctuated (Perner et al., 2011;Andresen et al., 2011;Wangner et al., 2018).Contrastingly, long-lasting episodes of surface water warming are observed within inner Disko Bay (c.1,600 -450 yrs.BP; c. 350 -1500 CE; Lloyd, 2006) and around Holsteinsborg Dyb located just south of Disko Bay (1150 -650 yrs.BP; 800 -1300 CE; Sha et al., 2012).Minor short-lived surface water warming around 1,000 -500yrs BP (1000-1500 CE) has been reported along the coast of Northwest Greenland (Mudie et al., 2005;Hansen et al., 2020).

Younger Neoglacial -the Little Ice Age
The period around 1300 CE signifies the early onset of a major oceanographic transition around Greenland (Fig. 4C), i.e., the onset of the LIA.The EGC intensified markedly around this time (Spielhagen et al., 2011;Werner et al., 2011;Bonnet et al., 2010;Divine and Dick, 2006) and this is evident from records tracking all the way down the east coast of Greenland (Kolling et al., 2017;Andersen et al., 2004a, 2004b, Jennings and Weiner, 1996, Miettinen et al., 2015) towards the west coast (Jensen et al., 2004;Sha et al., 2016;Roncaglia and Kuijpers, 2004).Using a large spatial network of well dated marine proxy sea-ice records Miles et al. (2020) has recently found robust evidence for extreme export of sea ice commencing abruptly around 650 yrs.BP (1300 CE).The onset was characterized by a century-long lasting pulse of ice along the East Greenland Current, which resulted in downstream increases of polar waters and ocean stratification thatremained generally sustained during the subsequent, coldest centuries of the LIA.
At the same time the warm IC intensified episodically during the LIA (Perner et al., 2016;Andresen et al., 2013;Jennings et al., 2011;Miettinen et al., 2015;Lassen et al., 2004;Seidenkrantz et al., 2008;Møller et al., 2006).A high-resolution study has shown that the LIA IC water warming on the shelf may have occurred episodically in response to lowered solar activity (Andresen et al., 2017).Lowered solar activity would increase the frequency of blocking highs and a negative NAO index (Barriopedro et al., 2008).This would in return shrink the subpolar gyre and displace it westwards (Barrier et al., 2014) resulting in a warming of waters in South Greenland (Andresen et al., 2017).The strong water mass stratification during the LIA with strong EGC as well as strong IC inflow towards Greenland resembles the conditions recorded during the DACP in South Greenland.However, the causes behind subsurface water warming episodes on the east, southeast and southwest Greenland shelf may not be the same during the LIA, the DACP and the RWP, since the Atlantic water temperature on the Greenland shelf is influenced by both the upstream North Atlantic Drift as well as the dynamics and strength of the subpolar gyre.
It is noted that the EGC intensification events around DACP and during the LIA were coeval with the so-called Bond events, which are surface cooling episodes in the Denmark Strait and North Atlantic region caused by intensive sea ice drift that has been related to decreased solar activity (Bond et al., 2001).

Sea-level estimates of glacier variability in Greenland
During the Last Glacial Maximum (LGM) eustatic sea level was ~130 -135 m lower than the present day (Clark et al., 2009;Lambeck et al., 2014), with the majority of freshwater being held by the Northern Hemisphere ice sheets including the Laurentide, Innuitian and Greenland Ice Sheets.Model reconstructions of the deglaciation suggest that during the LGM, in addition to its present volume of ~7 m sea-level equivalent (SLE) (Morlighem et al., 2017), the GrIS held an additional 4.65 m ± 0.70 m of excess SLE (Khan et al., 2016;Lecavalier et al., 2014).The ice marginal retreat and mass loss following the LGM, and in some areas a subsequent Neoglacial re-advance, combined to produce an integrated long-term visco-elastic signal, known as glacial isostatic adjustment (GIA), in response to the changing ice mass over time.Thus, the local and regional extent of former ice cover, combined with the local and regional earth structure, lead to a distinctive Late Holocene GIA.This was characterized by millennial-scale uplift and marine regression around the majority of the ice sheet and subsidence in the southwest, where an ongoing forebulge collapse, generated by the earlier ablation of the Laurentide Sheet combined with a Neoglacial ice readvance caused subsidence and a local marine transgression (Khan et al., 2016;Lecavalier et al., 2014).This Late Holocene pattern contrasts sharply with the isostatic uplift measured with GPS throughout Greenland today (Bevis et al., 2012(Bevis et al., , 2019)).

Past 2000 years of Greenland relative sea level
There are around 1200 sea-level index points from around Greenland that record changes in relative sea level (RSL) during the Holocene (Lecavalier et al., 2014).The greatest fall in RSL, on the west coast at Sisimiut, is from a peak of ca.130 m a.s.l. at 10,000 yrs.BP to a minimum of ca. 4 m below m.a.s.l at 2,000 yrs.BP (Bennike et al., 2011;Long et al., 2009) (Fig. 1).This pattern of rapid Early Holocene RSL fall to below the modern level in the Middle to Late Holocene, followed by a rise to the present level in the past few thousand years, is evident around much of the coast of Greenland, but is particularly pronounced in the West and Southwest (e.g.Bennike et al., 2011;Long et al., 1999Long et al., , 2003Long et al., , 2006Long et al., , 2009;;Long and Roberts, 2002;Sparrenbom et al., 2006aSparrenbom et al., , 2006b;;Woodroffe and Long, 2013).In summary, the isostatic rebound triggered by ice mass loss in the early Holocene was mitigated in the south and west during the late Holocene by a combination of Laurentide and Innuitian Ice sheet forbulge collapses and a Neoglacial GrIS re-advance.
As they are sequentially submerged, depressions in the bedrock known as isolation basins record the inflection from falling to rising RSL during the middle to late Holocene.The submerged isolation basin record from Nanortalik in Southern Greenland, shows that RSL reached a minimum of ~10 m below the highest modern tide level between 8,000 -6,000 BP, rising steadily from this low-stand to the present sea level through the Middle and Late Holocene at a rate of approximately 1.2 mm/yrs (Fig. 1) (Sparrenbom et al., 2006a).Establishing the elevation of the sills of these basins is not easy, however, and their sediments can be eroded during the marine ingression (Long et al., 2011).Therefore, RSL chronologies provided by ingression contacts in the sediment sequences collected from submerged isolation basins can provide RSL reconstructions of only millennial, metre-scale resolution at best.
An additional source of information on RSL changes during the last 2000 years comes from salt marshes, which can yield high precision RSL reconstructions (typically +/-0.2 -0.1 m in temperate settings by combining microfossil analyses and a transfer function approach; Gehrels, 2000; Horton and Edwards, 2006).Salt marshes in Greenland are found in low energy protected settings, at the head of tidal ponds or fjords, and the mouths of rivers, as far north as at least 75 • in West Greenland and 65 • in East Greenland.Salt marsh sediments have acted like natural, long-term tide gauges over the last millennium, but deeper, older saltmarsh sediments are not preserved in the nearshore environment, and there is often a gap between the end of the isolation basin record and the start of a nearby saltmarsh sequence (Woodroffe and Long, 2009).Studies of salt marsh sediments in Disko Bay, and at sites near Sisimiut, Paamiut and Nanortalik in the southwest, show that RSL started to rise starting ~700 years ago (Fig. 1).After this period RSL slowed down and stabilized from Sisimiut northwards over the past ~400 years.This process occurred later further south (Woodroffe and Long, 2009;Long et al., 2012: Woodroffe, pers. comm.).This rising RSL trend agrees with isolation basin data, and both indicate bedrock subsidence in South and West Greenland over the majority of the last two millennia.However, changes in RSL stability during and since, the LIA has been consistent with regional ice mass loss.The timing of this mass loss is debated as the saltmarsh data suggests it began during the LIA ~350 yrs.BP (~1600 CE) at Sisimiut, ~150 yrs.BP (~1800 CE) at Nanortalik (Long et al., 2012).The discrepancy between RSL observations and other proxies reconstructing the mass loss chronology underscores our incomplete understanding of the impact of GrIS mass balance changes on short and long-term processes related to RSL around Greenland today, as well as the impact of external contributors to sea level change and solid earth processes.

Ice marginal changes
We have divided the Neoglacial ice marginal fluctuations into advance, culmination and retreat phases.Age constraints on glacial expansion and recession are related to a suite of methods depending on the climatic setting.We use evidence from past changes in ice sheet elevation, reworked organic material, threshold lakes, and in-situ moss chronology to bracket the age of the advance phases, while 10 Be surface exposure dating and dated archaeological sites offer a limiting age for the culmination phases at sites in the terminal zone.Retreat phase dates are provided by lichenometric methods and historical documentation from maps and photographs (Fig. 5).

The response from ice elevation in central Greenland to marginal changes
The Greenland deep ice core records show regional differences in δ 18 O values through time that primarily result from different histories of atmospheric temperature, distance to moisture source, and ice-sheet elevation change (Vinther et al., 2009).Based on the assumption that temperature development on a millennial time scale was similar everywhere in Greenland, Vinther et al. (2009) extracted Holocene millennia scale ice sheet elevation histories at Camp Century, Dye 3, GRIP, and NGRIP from their respective δ 18 O ice core records.Later, Lecavalier et al. (2013) extended this work accounting for isostatic influence on the isotope signal assuming a spatial gradient in the temperatures across Greenland.
The dominant trends in Holocene ice sheet elevation records are largely a response to increased snow accumulation, bedrock uplift and glacier retreat at the margin following the termination of the last glacial period (Cuffey and Clow, 1997;Vinther et al., 2009).However, millennial scale oscillations in the records may result from Holocene fluctuations in snow accumulation and glaciation at the margin (e.g., Alley et al., 2010).Response-times for elevation adjustments at the central part of the GrIS have been estimated on the basis of a kinematic wave model from Nye (1960Nye ( , 1963) ) to be up to 5 millennia for regional snow accumulation change (Alley et al., 2010;van der Veen, 2001) and centuries for changes of marginal extent (Alley et al., 1987;Cuffey and Clow, 1997;van der Veen, 2001).This implies that millennia-scale oscillations in the Holocene ice sheet elevation record are more prone to reflect ice marginal fluctuations than regional snow accumulation changes (Cuffey and Clow, 1997).Marginal changes, on the other hand, may be associated with local changes in snow accumulation rates (e.g., Nye, 1960).The magnitude of ice sheet elevation response to changes in glacier extent has been attempted (e.g., Cuffey and Clow, 1997;Alley and Whillans, 1984), using the relationship between glacier height and width derived from analytical modelling (Vialov, 1958;Cuffey and Paterson, 2010).While ice sheet elevation is a result of long-term climate history, the marginal extent at a given time can be estimated inversely from the rate of elevation changes as pointed out by Nye (1963).Here we use the ice core derived Central Greenland elevation records from Lecavalier et al. (2013) in an extremely simplified inverse approach to obtain the history of changes at the ice margin on a millennial time scale.While the method does not quantify the extent of the marginal ice, it does let us approximate the onset of marginal ice responses.We assume that the shape of the GrIS is always adjusting in order to reach a steady state condition.This transient equilibrium state is never constant, but rather changing with a frequency considerably higher than ice sheet response time.The rate of ice sheet elevation changes (that is the first derivative of the ice sheet elevation as a function of time) will reflect the magnitude of imbalance between elevation and glacier extent, while the acceleration (that is the second derivative of elevation as a function of time) should track changes in the glacier extent.A positive acceleration indicates that the margin is advancing, while a negative acceleration indicates that the margin is retreating.When the sign of the acceleration change (that is when the second derivative of elevation as a function of time is zero), it reflects a change in the direction of glacier movement at the margin downstream from the ice core site, defining starting points for glacier advance or retreat, or equivalently maxima and minima for glacier extent.Although the calibration of stable isotope records to ice sheet elevation has an absolute uncertainty of up to 72 meters (Lecavalier et al., 2013), we assume that the elevation records have recorded times of changes in the elevation response to changes in marginal extent.This assumption is supported by comparing the elevation histories from the GRIP and GISP2 stable isotope records (Vinther et al., 2009, supp.).These were used to validate error estimates in the ice sheet elevation records, showing that although the elevations show some difference within the standard error, the two records show very similar millennia trends.In Fig. 6 we show the second derivative of the ice sheet elevation curves for GRIP, NorthGRIP, Dye 3 and Camp Century with periods of ice marginal advance and retreats indicated on the curves.Accounting for a time-delay in elevation response because glacial fluctuations at the margin is followed by a kinematic wave of adjustment migrating from the margin towards the central part the GrIS in 100-200 years (Cuffey and Clow, 1997;Alley et al., 1987), we have moved the advance/retreat starting points derived from the elevation curve back in time from the points of changing sign of acceleration in order to derive the approximate timing of advance or retreat onsets at the ice margin.We simulated the time development of elevation response to an abrupt marginal change using the response function for the GrIS Summit (Cuffey and Clow, 1997) and found that the second derivative of the elevation curve responds with a pulse starting after about 100 years, peaking after 200-250 years, and fading out within the next few centuries.Considering that elevation response delay-times increase with increasing distance from the margin, we have set the delay-times in the elevation response to onsets of marginal advance and retreat to 200 ± 150 yrs.for the GRIP and NorthGRIP, and 100 ± 100 yrs.for the Dye 3 and Camp Century records.Further, considering that the elevation response to marginal changes last a few centuries we have set the delay-times for the termination of advance and retreat to 400 ± 150 yrs.for the GRIP and NorthGRIP, and to 200 ± 100 yrs.for the Dye 3 and Camp Century records.
We also assume that the ice sheet elevation is most sensitive to changes at margins closest to the coring site.Spatially across GrIS, the Dye 3 elevations respond to marginal changes in East Greenland around 65 • N, GRIP elevation responds to marginal changes in East and West Greenland around 72 • N, NorthGRIP elevations respond to marginal changes in East and West Greenland near 75 • N, and Camp Century elevations respond to marginal changes in Northwest Greenland around 77 • N (Fig. 1).Accordingly, the first pronounced ice advance phase after the HTM (8,000 -5,000 yrs.BP) started at around 4,800 BP in Northwest Greenland (77 • N).Probably a brief advance also occurred in both Central and Southeast Greenland within the uncertainty of the timing.A second advance phase began around 3,600 years BP in Southeast Greenland, and around 3,000 years BP in East and West Greenland 72 -75 • N (probably following a short-term advance at around 3,800 years BP) and lastly in the Northwest Greenland (77 • N) some 2,100 years BP.A third advance phase began in Central East and West Greenland around 900 BP years at 75 • N and 700 years BP at 72 • N.
The ice margin history in Northwest Greenland around 77 • N deduced from the Camp Century elevation record seems to be decoupled from the ice margin history from Central and Southeast Greenland.In Northwest Greenland periods of ice-retreats correlate with the cold periods observed in the GISP2 nssK + record (Fig. 3) suggesting a relationship between ice extent and local precipitation at the ice margin in that part of Greenland.

Reworked organic material
Organic remains retrieved from end moraines or exposed on ice surfaces during downwasting records the timing of ice advances over ice-free land.We assembled 14C dates from 132 sites representing every drainage area of the GrIS, published from the 1970s to present (Table 1, Fig. 7).Marine samples have been corrected for a reservoir effect, ΔR, of 0 yrs.for West Greenland, -150 yrs.for East and 50 yrs.for North Greenland following the standard used for Greenland 14 C dates (e.g., Mörner and Funder, 1990).The reservoir effect is related to the extent and annual duration of sea ice, which inhibits atmospheric CO 2 exchange (Tauber, 1979;Hanslik et al., 2010).Therefore, marine ages from the HTM and LIA could have been respectively slightly higher and lower than shown here (Funder et al., 2011).Conventional dates on shells from studies carried out in 1970-80s often consist of mixed assemblages of different shells in order to secure enough material for a date (e.g., Håkansson, 1972;Weidick, 1975;Kelly, 1980).To better understand the age variability of these mixed assemblages, we collected fifteen intact or fragmented shells of Mya arenaria, Hiatella arctica and Astarte sp. in diamict sediment at the terrestrial part of the southwestern front of Humboldt Glacier, Northwest Greenland in 2017 (Fig. 8).These were all dated individually with AMS and show an age range between 3,600 to 500 yrs.BP, which yields a maximum (oldest) age of 500 yrs.BP (1500 CE) for the LIA advance at the Humboldt Glacier (Søndergaard et al., 2020).On the northern flank, Bennike (2002) dated reworked shells and found an advance about 150 years earlier.Either this attests to a later maximum along the northern margin of Humboldt or more likely, that the youngest shells have not been identified and dated.However, it also raises concerns about the earlier mixed shell ages, as these might be a result of individual shells with different ages.In the dataset we have assembled, there are also examples to the contrary, where ages of mixed samples (Hjort, 1997;Hjort et al., 1998) do match recent AMS results (Larsen et al., 2018).Therefore, as a cautionary measure mixed sample dates have been excluded from the probability and stack plots in Figs. 9 and 10.
Radiocarbon ages on reworked organic materials testify to a time when the ice was behind its present extent (smaller than present) with the youngest age providing a maximum limit for the subsequent ice advance at any particular site.Thus, we only use dates from the LIA moraine/sediments, which provide a maximum age of the LIA advance.In general, the distribution of ages across the GrIS shows that the youngest reworked material from the LIA in Northwest and Southeast Greenland is mostly younger than 2,900 yrs.BP.In West Greenland, it is older than 3,000 yrs.BP and in South and Northeastern Greenland, it is older than 5,000 yrs.BP with many Early Holocene ages.We believe this reflects the availability of coastal land during deglaciation whereby the Northwest and Southeast Greenland are steep mountainous terrain, leaving a narrow area of ice-free land at any time.The maximum ages fall within two periods, 677 -607 yrs.BP (1273 -1343 CE) and 341 -308 yrs.BP (1609 -1642 CE), with a slight tendency for the East and Southeast Greenland to be the youngest.
In a cumulative stack (Fig. 9), the first cluster of ages is pre-LGM, 39,000 -26,000 yrs.BP and except for a single date at Frederikshåb Isblink in Southwest Greenland all were collected at Storstømmen Glacier in Northeast Greenland (Fig. 7).Weidick et al. (1996) relate this to a smaller than present ice extent.Larsen et al. (2018) suggest that the Northeast Greenland Ice Stream was ~20 -70 km behind its present extent.The isolated date from Frederikshåb Isblink is from a mixed shell sample and before the result has been reproduced and confirmed, we will not speculate further.The second cluster of dates range from 10,500 -7,700 yrs.BP.There is a paucity of dated material from 7,400 to 6,000 yrs.BP that might relate to glacier expansion during the older Neoglacial.This may have eroded the HTM sediments and organic material.However, the widespread absence of minerogenic sediment in threshold lakes contradicts this possibility.There is a plateau with a few dates from   2,700 to 1,700 yrs.BP, which may also represent a low probability interval associated with glacier expansion (Fig. 10).After the MWP, we find two periods with little reworked organic material during 1,050 -750 yrs.BP (900 -1200 CE) and 500 -125 yrs.BP (1450 -1825 CE), which matches the younger Neoglacial with two overall periods of glacier expansion during the LIA even though this cannot be resolved further given the restricted quantity of dates.

Threshold lakes
Threshold lakes only receive meltwater and glaciogenic sediments (inorganic, minerogenic) during periods in which glacier ice is within the lake catchment, whereas non-glacial organic-rich sediments (gyttja) are deposited during intervals with no glacial ice extent in the catchment (e.g., Briner et al., 2010).The sedimentary records from lakes with catchments known to be sensitive to glacier proximity can provide information on past ice marginal fluctuations.
In general, the GrIS began to grow again after the HTM until it reached a maximum dimension during the LIA (Funder et al., 2011).This is reflected in threshold lake sediments by a succession of organicrich gyttja below a unit of clastic sediments representing the advancing glacier during the LIA.In addition, some threshold lakes also record a return to organic sedimentation above the clastic sediments as the glacier retreated after the LIA (e.g., Kaplan et al., 2002).However, depending on the distance between the lake and its threshold, the record in different lakes will have different local representations of these phases.
We synthesized most of the published threshold lake records that include a LIA sequence (n=48) (Fig. 11), thereby refining and expanding previous reviews of larger spatial and Holocene-scale patterns of change (McKay and Kaufman, 2014;Sundqvist et al., 2014;Briner et al., 2016).It is worth mentioning that the threshold lake records that we chose not to include here due to their unconstrained LIA signal have been instrumental in constraining local glacial history during this period when supported by other types of records such as 10 Be exposure dating of boulders, and 14 C dating of reworked marine mollusks and subfossil plants (e.g., Bennike and Sparrenbom, 2007;Briner et al., 2013aBriner et al., , 2013b;;Kelley et al., 2012;Larsen et al., 2015;Levy et al., 2017;Larsen et al., 2019;Søndergaard et al., 2019).Moreover a few of the selected threshold records have been chronologically improved for this study.These are marked as unpublished.We selected ages closest to the transition to minerogenic sediment.For the majority of the records, these ages represent a maximum limiting age for the glacial onset, as typically most datable materials are found within the underlying organic layer (Table 2).At Primary Lake, Northwest Greenland, for instance, both bulk sediment and terrestrial macrofossils taken within <2 cm of the gyttja/minerogenic interface from replicate cores, yielded ages with 80 years difference i.e., 690 yrs.BP (1260 ± 30 CE) and 610 yrs.BP (1340 ± 50 CE), respectively (Briner et al., 2013a).The age model by Håkansson et al. (2014), which considers different sedimentation rates based on lithostratigraphic units provides a LIA advance age of 240 yrs.BP contrary to 590 yrs.BP, yielding a ~350 yrs.difference for 2 -4 cm sediment depth distance.This type of systematic dating bias may cause overall increasingly older ages for the advance at higher latitudes as organic material accumulation and the sedimentation rate generally decrease with latitude.Lakes with no age constraint on the minerogenic layers (>15 cm offset from transition) were excluded from the analysis (n=6), except for modelled ages from two sites in East Greenland (Round Lake and Emerald Lake; Lowell et al., 2013;Lusas et al., 2017).These were included as references for the distribution compilation only due to the paucity of suitable sites from that region (Fig. 11).Where the surveyed lakes recorded a LIA advance that was unconstrained within a continuous minerogenic layer following an older Neoglacial advance (<2000 yrs.BP), the age for the earliest advance was selected.
Approximately half of the dated materials at the organicminerogenic transition (22 of 42) were not terrestrial macrofossils but rather humic acids, bulk sediment organic carbon or aquatic macrofossils (e.g., moss, Daphnia ephippia).These materials generally yield coherent ages in their respective contexts (i.e., stratigraphic order) for the lakes included here.These chronologies may integrate an overall older signal due to incorporation of older carbon from reworked organic material, e.g., reservoir effects (e.g., aquatic moss Drepanocladus exannulatus; Olsen et al., 2012b) and methane-oxidizing bacteria consumed by invertebrates incorporated in the sample (Schilder et al., 2015).However, this bias is on a centennial to sub-centennial scale and not significant here given the relatively low temporal resolution of the agedepth models.All radiocarbon age estimates were calibrated to years before present using the IntCal13 calibration dataset (Reimer et al., 2013).Though most of the published literature on threshold lakes in Greenland synthesized here report age estimates using mean values, we found that these yields -50 to +28-year offsets compared to median values.We therefore used the latter as a more accurate presentation of the 14 C probability distributions.Where provided by the authors, we tested the effect of including glacier advance age estimates based on their model rather than that of the closest measurement.
As expected, we found that generally younger ages were modelled by as much as 400 years (e.g., Kaplan et al., 2002;Larsen et al., 2011Larsen et al., , 2015Larsen et al., , 2016Larsen et al., , 2017;;Levy et al., 2017).Using age estimates from context-specific and reliable models built consistently across the dataset would likely improve the accuracy of these estimates.However, the error associated with these models, which is larger where no measurements are available, would not necessarily yield more accurate or precise estimates over the most recent ~500 years than those of the probability distributions of the calibrated radiocarbon age determinations used here.
Our synthesis of the threshold lakes dataset (n=42; 7 -815 m.a.s.l.; 61.0 -73.0 Neoglacial advance of the GrIS or local glaciers and ice caps.This was followed by sustained minerogenic glacial inputs to the lakes until the 20th century, ~50 yrs.BP (Fig. 11).Based on the summed probability densities of the radiocarbon ages closest to the transition, this period could extend to as early as 1150 yrs.BP, while the most synchronous glacial advances for all regions combined occurred between 740 and 480 yrs.BP (Fig. 12).Approximately half (n=22) of the threshold lakes surveyed record the largest Late Holocene extent during the LIA.(Fig. 12) In these records, the first glacier advance occurred from 1,150 yrs.BP and was followed by sustained meltwater inputs to lakes until the 20th century (~50 yrs BP) or present, as observed from aerial photographs and satellite imagery in all regions surveyed.A slightly earlier onset of glacier advance was recorded at 'Ymer lake' in Southeast Greenland (omitted from this re-analysis) at a median age estimate of 1,279 yrs.BP (van der Bilt et al., 2018), similar to other sites in the Southeast (Fig. 13).A number of lakes, however (n=10 i.e., Sports Lake, Lake Lucy, Sikuiui Lake, Frederikshåb, IS5, IS21, Loon Lake, Kulusuk, Saqqap Tasersua and Pauiaivik Lake; Table 2) record more variable conditions than described above over the Late Holocene.At these lakes, multiple advances indicating glacier fluctuations during both the older Neoglacial and LIA occurred within the threshold level.The remaining lakes (n=9) record an older Neoglacial advance with no distinct LIA signal within the minerogenic layer.Only two lakes (Storesø and Rundesø) in Southwest Greenland capture older or potentially older, organic sequences that may indicate a larger Neoglacial compared to the LIA threshold scenario (Larsen et al., 2016).However, poor age constraints at these sites preclude establishing the timing of an older Neoglacial onset.Furthermore, these lakes receive meltwater input during glacier fluctuations that were related to local climate and/or topography independent of regional ice sheet fluctuations (Larsen et al., 2016).

In situ plant remain chronology
Dating of formerly ice-entombed in situ moss is another approach to determine the pattern and timing of the LIA advance phase and ice buildup.Studies from West and East Greenland sampled cold-based sectors of polythermal ice caps or passive sectors of the GrIS (Lowell et al., 2013;Schweinsberg et al., 2017Schweinsberg et al., , 2018)).Ages of subfossil plants at the Istorvet Ice Cap, East Greenland were initially interpreted as a time interval with ice-free conditions and smaller than present ice extent.Accordingly, ice cap extents were inferred to be smaller than present until about 1,845 yrs.BP (1025 CE) (Lowell et al., 2013).This interpretation was subsequently debated by Miller et al. (2013), who revaluated the data filtered through a different sampling protocol (Miller et al., 2012).These moss ages are now interpreted as past episodes of snowline lowering that led to ice cap expansion, with the youngest age representing the final kill date of plants before burial.In the case of Istorvet Ice Cap, this implies that sometime after 1,150 yrs.BP (800 CE) the snowline was lowered and persistent snow or ice cover occurred by 1,050 -950 yrs.BP (900 -1000 CE), with the independent ice cap finally reaching its Neoglacial maximum by 825 yrs.BP (1125 CE) (Fig. 14) (Miller et al., 2013).

10 Be Surface exposure dating
Surface exposure dating ( 10 Be dating) is a well-established method for providing an age of deglaciation at terminal moraines.It is a measure of the exposure of cosmic rays and accumulation of 10 Be isotopes in quartz minerals since the ice cover melted away.There are still relatively few studies using 10 Be dating of LIA moraines.Two studies concern LIA moraines related to the GrIS (Reusche et al., 2018;Levy et al., 2018) whereas four studies concern LIA moraines adjacent to local glaciers and ice caps (Kelly et al., 2008;Winsor et al., 2014;Young et al., 2015;Jomelli et al., 2016).
The 10 Be ages from the two GrIS sites show that the LIA culminated between ~300 to 100 yrs.BP (Reusche et al., 2018;Levy et al., 2018).In contrast, the LIA culmination phase of three local glaciers and ice caps is older and dates between ~900 and 1300 yrs.BP.Additionally, there is sporadic evidence of a Neoglacial phase around ~2800 to 1600 yrs.BP where the GrIS and local glaciers were equal to or larger than the LIA extent (Kelly et al., 2008;Winsor et al., 2014;Young et al., 2015;Jomelli et al., 2016).

Archaeological sites
From almost all parts of Greenland, there are ruins and archaeological sites situated in proximity to glacier front.Greenland was first settled by Paleo-Inuit of the Arctic Small Tool tradition around 4500 BP.Several Palaeo-Inuit migrations such as Greenlandic Dorset (2750-1950 yrs.BP) and Late Dorset followed, and around 700 BP Inuit spread from the Nares Strait region north around Greenland as well as south into West Greenland.All of these hunter-gatherers explored and utilized the varied ecological niches and left behind settlements with stone structures including hearths and tent rings used to anchor the tent skin firmly to the ground and sometimes graves or kayak supports like those shown near Nuulliit Glacier in Fig. 15.Architectural features such as these are unlikely to survive overriding glaciers.The location and age of these remains offer evidence of a maximum possible glacier extent at a given time and can in certain instances supply an upper boundary for the LIA extent where no other data is available.
The list of 21 ancient settlements included here show that glacier fronts never advanced beyond these localities since abandonment  (Table 3).In some cases, there is circumstantial evidence suggesting that the glacier or ice sheet margin is likely to have retreated considerably behind its present location to accommodate settlements at these localities, particularly for some of the Paleo-Inuit localities that were abandoned prior to the Neoglacial ice advances (Fig. 15; Table 3).The distances from the sites to the nearest glacier have been extracted by using the measuring tool integrated in http://nunniffiit.natmus.gl/spatialmap.Several errors may have accumulated in these measurements: i) the location of sites in many places has been extracted from older analogous maps, descriptions, publications and hand drawn sketches; ii) the glacier front depicted on the map underlying the digital register is Greenland 1:250.000.Dating of the localities has been inferred by typological observations most often relating to dwelling types and in some instances from artifacts.The broad date ranges of the general culture historical scheme of Greenland (Gulløv, 2004) have then been adjusted according to well-documented local variations to this scheme (Table 4).The Paleo-Inuit occupation phases of northern Greenland are thus many hundred years shorter than the comparable settlement phases in more resource rich West Greenland (Jensen, 2006;Grønnow and Jensen, 2003), and the Thule migration into and subsequent settlement of Northeast Greenland is also known to have had a different trajectory than in the rest of Greenland (Sørensen and Gulløv, 2012).
The proximity to the 1:250.000map glacier position means that the LIA termination was within 500-1500 m from the archaeological sites and the age of the settlements exclude these glaciers from occupation after 225 yrs.BP (1725 CE).Interestingly, sites in western, northwestern and northern Greenland show that they have not been overridden in the earlier Neoglacial as far back as 4,450 yrs.BP.The significance being that pre-LIA-ice advances were close to the LIA termination or may have been less significant and have been overridden by the LIA advance (Fig. 15).

Lichenometric records
Lichens have been used for over 60 years to date the onset of vegetation in Greenland and thereby infer the post-LIA exposure of bedrock or boulder surfaces (Hansen, 2010).While this technique may seem relatively straightforward, many environmental parameters influence the growth rate of the lichens.It has been estimated that lichenometry yields reliable ages for the past 300 years with an approximately 20% uncertainty in West Greenland (Werner, 1990).For older periods the uncertainty of the lichen ages is larger as they are calculated using an unrealistic linear growth rate (c.f.Forman et al., 2007).The recent literature from Greenland was recently reviewed by Hansen (2010) and since then few new studies have been published based on dating moraines using lichenometry.Accordingly, we only provide a brief summary of the main results, which constrain the retreat phase following the LIA maxima.Most of the studies are from local glaciers and ice-caps in Southwest and Southeast in addition to a few studies from the GrIS margin in Southwest Greenland (Figs. 16,17,Table 5).

Geomorphological imprints by LIA glaciers
The LIA cooling left a prominent geomorphological imprint in the landscape when outlet glaciers from GrIS and local glaciers and ice caps expanded beyond their present margins.Areas with channelized ice flow, along fjords and valleys, or unchannelized ice flow, from more passive ice lobes, both show a current landscape with fresh nonvegetated surfaces covered with scattered boulders.This reflects the former extent of the glacier, where ice abraded and eroded the previously deglaciated landscape free of vegetation, soil, and old sediments (Fig. 18).At higher elevations, the transition to the intact vegetated surfaces marks the former ice contact.When traced downstream at marine-terminating glaciers, this contact often ends at sea level or at land-terminating glaciers, often at end moraines that indicate the former ice marginal position.The so-called trim zone reflects the ice volume lost since the ice retreated from its advanced position (Weidick, 1959;Kjeldsen et al., 2015).
We explored the Greenland-wide spatial distribution of trimlines and end moraines using a glacier compilation of 243 outlet glaciers from the ice sheet wider than 1.5 km in size (Rignot and Mouginot, 2012).As these are mostly reported from ice fronts related to tidewater and laketerminating glaciers, we supplemented this record with an additional outlets including land-terminating glaciers (Fig. 19).The location of the 308 outlets were subsequently mapped against a 25 m gridded DEM and a 2 m black-and-white digital orthophotograph (Korsgaard et al., 2016) for evidence of trimlines and end moraine landforms.We find that of the outlets have a clear trimline or end moraine.The remaining  Holtved (1936) in (Holtved, 1954) and Jens Fog Jensen (2015).Note the heavily diminished glacier tongue.See Tables 3 and 4. (C) Principal sketch of human settlements on the island of Pullassuaq and glacier front of Sermersuaq (Humboldt glacier) in Kane Basin.Pullassuaq holds a large number of Paleo-Inuit settlements assumed to date from Saqqaq to Late Dorset (4,450-500 yrs.BP) as well as a few Thule culture settlements (700-225 yrs.BP).Along the southwestern and southern side of the island there are moraine ridges indicating that at some point of time, probably at the maximum LIA advance, the glacier front reached the island, but did not override it.The dwelling ruins were left undisturbed.When Pullassuaq was settled in times of the Paleo-Inuit, the glacier front most likely had retreated even further than its position until 1990'ies when it was approximately one km from the islet.outlets are surge-type glaciers or large outlets with long and wide calving fronts producing unrecognizable trimlines or trimzones.
It is possible to assess post-LIA retreat from these terminal positions to the mid-1980s using the same orthophotographs (Korsgaard et al., 2016).The maximum extent of land-terminating glaciers is mapped at the outermost moraine crest, while marine-terminating glaciers are mapped at their most distal trimline or moraine (Fig. 19).Many glaciers along the northwest and southeast coasts have retreated significantlymore than 10 km -since their LIA culmination.Some larger outlet systems, such as Jakobshavn, Kangiata Nunata Sermia (KNS) and Kofoed-Hansen Brae, have had a maximum retreat since the LIA of more than 20 km.The pattern of these retreats mirrors the well-known spatial distribution of changes reported during for last two decades (Kjaer et al., 2012;Bjørk et al., 2012: Moon et al., 2012;Vijay et al., 2019;Mouginot et al., 2019) and the pattern of post-LIA thinning (Kjeldsen et al., 2015;Khan et al., 2014;Khan et al., 2020).
The surface change assessed by Kjeldsen et al. (2015) can be used to evaluate the surface lowering associated with ice marginal retreat.In order to assess the post-LIA behavior of each outlet glacier of the GrIS, we use the data to estimate the average elevation change within a specific area of each outlet glacier.As the glaciers vary in size, configuration, and setting, and in their post-LIA response, we extracted data using two approaches.For marine terminating glaciers, we created a 2 km buffer around the 250 m height above ellipsoid (HAE) contour based on ArcticDEM (Polar Geospatial Center from DigitalGlobe, Inc. imagery), and extracted elevation data within this area.For land-and laketerminating glaciers, we created a 2 km buffer around the lowest 10 km of the main flowline, relative to the 1978/87 CE ice front position, and extracted elevation data within this area.
In the period between the LIA maximum (assumed to be c.1900 to be comparable to Kjeldsen et al., 2015) and 1983, the glacier retreat matches the surface lowering rate for Northwest and Southeast Greenland, as well as for the large outlets at Jakobshavn Isbrae, KNS, Frederikshåb Isblink, Qassimiut Lobe, Kangerlussuaq Glacier and Helheim Glacier, with rates of vertical change between 1 and 10 m annually.Interestingly, the north and northeast coast, including the marine outlets at Zachariae Isstrøm and Nioghalvfjerdsfjorden Glacier, show restricted surface lowering although some of them have retreated significantly (Figs. 1,19).This is most likely because of detachment of floating glacier tongues, which does not result in an immediate dynamic thinning or marine drawdown effect on the surface elevation before retreat of the grounding line.
At the interface between the former ice surface, or trimline, and the higher lying terrain, lateral moraines might appear.In Greenland, these are commonly expressed as moraine deposits with oversized boulders that have been pushed, or rolled, by the glacier (Fig. 18).Sorted sediment might be found at the interface, sometimes in ponds or plateau deposits, which reflect outwash activity by meltwater flowing downstream.In cases where the glacier has retreated and the surface has lowered, a series of lateral moraines might be visible.We speculate that these reflect consecutive still stands, or maybe even minor advances during a general retreat phase (Fig. 18).In glacier forefields, this landscape is matched by a series of recessional end moraines expressed as successive lobate shaped raised terrain often consisting of unsorted moraine deposits.
The terrain surface has often been eroded since moraines were formed (Fig. 20).Ice-cored moraines show collapse features such as sinkholes or terraced slopes, however mostly, meltwater erosion has removed large parts of deposited landforms (Fig. 21).(See Fig. 22.)

Historical maps and glacier mapping
As the termination of the LIA falls within the era of historical records in Greenland, there is a possibility that the onset of retreat has been captured in these records.A historical review of glacier front positions in West Greenland has previously been published on the basis of historical maps and written accounts (Weidick, 1959(Weidick, , 1968)).The date for these maps varies from before, during, and after what is presumed to be the termination of the LIA (1800 -1920) (Fig. 21, Table 6).In the historical reviews (Weidick, 1959(Weidick, , 1968) glacier lengths and frontal positions are given for numerous glaciers in Southwest Greenland.Glacier positions were established by comparing the old maps with newer maps and aerial photos, as well as with the written testimonies from the original publications.The main purpose for re-assessing the historical maps in this review is to examine whether the maps can provide accurate information regarding the glacier marginal positions in relation to the mapped LIA extent.Here, we present a suite of historical maps from West and   5.
North Greenland.East Greenland is excluded, due to the lack of suitable maps of non-surging glaciers.
We select historical maps, in which the mapping of the glacier margins is original and not copied from older maps.We also select the best available maps from the period prior to 1920.We exclude maps where the general cartography is obviously wrong, or where the travel route during mapping is known and far from the mapped glacier margin (Table 6).Lastly, we also exclude maps of glaciers with surge behavior, where these are identified based on looped moraines, glaciers with crevasse squeeze ridges in the fore field, glaciers advancing more than 50 m yr− 1 and glaciers with a collapsed accumulation zone (Bjørk et al., 2018a(Bjørk et al., , 2018b)).
We georectify the historical maps on the basis of a 2-meter digital orthophoto (Korsgaard et al., 2016), which in this part of Greenland is mainly from imagery recorded in 1985, supplemented with images from 1981.By placing tie-points near the glacier fronts at easily recognizable landforms such as headlands, islands, and river outlets, the historical maps are then tied down to the precise orthophoto with the most suitable transformation, distorting the original map in either one, two or more dimensions.In cases where multiple glaciers appear on the same map, individual rectifications have been carried out on each glacier front.The end-result is a historical map with as precise coordinates as possible in the vicinity of the glacier front.
Using the georectified historical maps, we measure the distance from the mapped ice front and the LIA maximum position established from the 1981/85 aerial images.To account for positional uncertainties, and to give an estimation of the accuracy of the georectified historical map, we consider the line width with which the map has been drawn and measure the residual between known points or coastlines (Bjørk et al., 2012).As also mentioned in the historical review (Weidick, 1968), there is uncertainty related to the use of historical maps.The quality of historical maps varies, which is reflected in the overall accuracy presented.For a number of glaciers, using the historical maps to pinpoint the frontal position is redundant as the sum of errors far outweighs the distance interval for post-LIA fluctuations.
Evaluating the historical maps, we find that the line widths with which the glacier fronts have been drawn, ranges between 25 and 300 m.The average line width is 100 m for all georectified maps.The positional residual between mapped and true geographical feature also varies greatly, from c. 100 m to more than 2 km, depending on the scale of the map, the quality of the map, and the distance to the nearest tiepoint.This leads us to conclude that great carefulness must be used when incorporating glacier front positions from historical maps.
While we assume that mapped features are correct, there are instances where glacier features have been falsely mapped, and therefore omitted from this study.One example is the Jungersen Glacier mapped by Lauge Koch in 1916 during the 2nd Thule Expedition (Koch, 1932).Traveling on a dog-sledge, Lauge Koch reached a massive impenetrable ice wall and mapped it as the calving margin of Jungersen Glacier.The real Jungersen Glacier calving margin, however, is currently 60 km inland, and the fjord landscape shows no evidence of LIA presence in the fjord, or anywhere near Koch's map, nor does later mapping support this advanced position.After reading the expedition journals and their descriptions of the glacial features, we compare the travel routes to the produced maps and exclude glacier mappings where our confidence in their accuracy is low.
As we are primarily interested in the frontal position of the glacier in relation to what is considered to be the maximum position of the LIA (LIA-max), the distance from the glacier front deduced from the georectified historical map to the presumed LIA-maximum position from aerial photos has been measured.We divide the lengths into five groups (Fig. 21); ±500 m, 500-1500 m (larger or smaller), and >1500 m (larger  or smaller).
There is a trend in the spatial distribution of the glaciers that have historically mapped glacier fronts in front of their LIA-max positions.They are concentrated in the north of the Uummannaq region (north of Disko Bay) and in the Melville Bay region.This suggests that, at the time the historical maps were made, the terminal positions of the glaciers in these areas were at a maximum position, or very close to it.In South Greenland, the trend is different.Here, a majority of the georectified maps show the present glacier fronts at a retreated position in relation to their LIA-max positions.Kangiata Nunâta Sermia stands out, with c. 14 km retreat from the LIA-max position in 1889 (Jensen, 1889), which has been speculated to be the LIA maximum position in the 1700s (Weidick et al., 2012).This region is also the region where the oldest maps have been used, which indicates a marked difference in retreat from the LIAmax along the north-south axis.Despite these regional trends, the historical mapping of the Melville Bay region in the early 1920s also shows that neighboring outlet glaciers may have reacted very differently to the post-LIA warming.

Threshold lakes
Threshold lakes also bear witness of ice retreat from the lake catchment and is thus also a retreat phase proxy.In almost 70% of the cases surveyed, the LIA (or earlier) advance recorded by an onset of minerogenic material in the sediment appears to have continued, or for conditions for minerogenic inputs in the lake's catchment to have been maintained until present, as per documented field reports of water turbidity level or historical imagery.In the remainder ~30% of sequences that do register a return to organic sedimentation following the LIA advance, documentation of the timing is most often reported in relation to latest available (or consulted) historical imagery, either being bracketed by two series of images (Iceboom Lake: 1964-1985 CE;Loon Lake: 1959-1985 CE;Briner et al., 2010) or the most recent image e.g., Sports Lake <1985 CE; Young and Briner, 2015; Table 2).The establishment of chronologies using a combined macrofossil-based AMS 14 C calibrated against both IntCal13 and atmospheric 14 C curves (Bomb 13 NH) and continuous 210 Pb and 137 Cs measurements of these organic sequences would contribute to improve constraints on the most recent glacial retreat at these sites.

Discussion
Our assemblage of proxy records from different archives represents a first attempt to generalize Neoglacial and LIA glacier responses in Greenland using data of varying spatial and temporal resolution, as well as uneven distribution across Greenland.With this in mind, however, we have summarized the possible trends during the different glacier phases across the different records and compare them with different climate proxies.

Reconciling the advance phase proxies
Ice-core records predict the possibility of three distinct ice advance events during the last 3500 years (Fig. 6) with one event occurring at 3,000 -2,200 yrs.BP in Central East and Central West Greenland (72 -75 • N), and at 3,600 -1,800 yrs.BP in Southeast Greenland.Another advance event occurred at 2,100 -1,300 yrs.BP in Northwest Greenland (77 • N) and a third event began 900 yrs.BP in Central West, and East Greenland (72 -75 • N).Higher resolution ice fluctuations within the last 500 years are very difficult to resolve from ice-core elevation data alone (Figs 2,3,6 and 23).These three events of ice advance are all matched in the terrestrial records.Dates from reworked material indicates glacier expansion in the marginal area during two periods in the older Neoglacial; c. 2,700 -1,700 and 1,250 -950 yrs.BP, whereas one overall advance occurred in the younger Neoglacial; 700 -50 yrs.BP (1250- 1900 CE), which represents the LIA advance (Fig. 23).The less prominent period of expansion at 1,250 -950 yrs.BP is not seen in the ice core record.
Threshold lakes show minerogenic material in lake catchments from 4,400 yrs.BP until the 20th Century (Fig. 12).The first appearance of glacier-derived sediment in the threshold lakes intensified from 1,150 yrs.BP (800 CE) until 350 yrs.BP (1600 CE).We interpret these threshold lakes to reflect marginal advance from 1,650 yrs.BP and widespread advance after 1.150 yrs.BP in all drainage zones covered by threshold lakes data.The continuation of minerogenic material across the MWP likely signifies a near-maximum position of the ice margin throughout the younger Neoglacial period.The threshold lake compilation also reveals some spatial structures.Since 1,150 yrs.BP, where onsets of glacier-fed lakes have been recorded for all of the Greenland sites, earlier onsets are recorded at more northern latitudes (latitudinal gradient: West<Southwest<South) and at western longitudes (longitudinal gradient: West<East; Southwest<Southeast) (Fig. 12).However, relatively few ages are available to represent the eastern regions and this latter relation should therefore be considered with caution.Onsets of glacier-fed lakes are latest in South Greenland, beginning only around 630 yrs.BP.This is consistent with lake records in the South experiencing relatively late LIA culmination (~50 yrs.BP) and warmest temperature lows compared to more northern regions.In Southwest Greenland, glaciers began delivering meltwater inputs from 930 yrs.BP and as early as 1,150 yrs.BP in the western drainage zones.Because nearly half of the records were retrieved from West Greenland, glacier fed lakes from this region contributed most to the total sum probability densities (Fig. 12).Finally, the in-situ moss chronology shows snowline lowering around 1,900 yrs.BP interrupted by an ice-free interval indicated by a vegetation free period around the time of the MWP (950-750 yrs BP).Spatially, East Greenland seems to react later than the West coast records, but the data recovery of in-situ moss also remains still too scarce to draw any firm conclusion from the variable data.
Although the marine records around Greenland reveal complex variability in the cold surface and warmer subsurface currents, some major trends can be detected; the majority of records indicate a strong transport of warm Atlantic waters to Greenland at the end of the oldest epoch, which matches the European RWP (Figs. 4,23).Most records also detect a general strengthening of the EGC after c. 1,500 yrs.BP during the European DACP.The DACP was furthermore characterized by cooling of Atlantic waters by East and West Greenland, which contrasts with the synchronous strong IC inflow to Southeast and South Greenland.The oceanographic change mimics the advance phases in some coastal records, namely the threshold lakes and in-situ moss chronology, but not the ice-core prediction of marginal change.Most sea ice-influenced time series from East Greenland record milder conditions during the MWP (c.1,000 yrs.BP) than during the DACP, whereas the records from West Greenland indicate a more fluctuating and spatially heterogeneous climate regime during the MWP.Despite this heterogeneity the timing agrees with all advance phase records and there seems to be a pause in the Neoglacial glacial activity before entering the younger Neoglacial -LIA.
The period around 800 -700 yrs.BP thus signifies the early onset of a major oceanographic transition around Greenland (Fig. 4C), which we interpret as the onset of the LIA.The EGC intensified markedly at the time, which is evident from records tracking all the way down the east coast of Greenland towards the west coast.At the same time, the warm IC occurrence by Greenland intensified episodically during the LIA.The WGC, which had started to cool during the DACP, experienced a further intensification of cooling during the LIA near and by Disko Bugtimplying that the LIA intensification of the EGC was relatively more significant than the concurrent intensification of the IC.In Northwest and West Greenland, the Baffin Bay sea-ice production increased during the LIA with the Baffin Bay sea-ice edge moving into the vicinity of Disko Bay.Marine IRD sediment records from Sermilik Fjord by Helheim Glacier in SE Greenland indicate that this large outlet glacier had a stiff melange and/or a floating ice tongue since at least the DACP and lasting until the onset of the 20th century.It has been hypothesized that the glacier did not lose its tongue or melange until sea ice decreased unprecedentedly at the onset of the 20th century implying that sea ice associated cold conditions buffered the glacier in spite of warm IC inflow during the DACP (Andresen et al., 2017).In contrast, Jakobshavn Isbrae was a calving glacier during the warm RWP, DACP and MWP and only started to grow an ice tongue, according to marine IRD sediment records, at the onset of the LIA (Wangner et al., 2018).
Across all proxy data and bearing in mind the cooling signal detected in the ice core, 'climate' lake records and marine data, we suggest initiation of glacier advances in all of Greenland during the LIA commenced ca.700 yrs.BP (ca.1250 CE).Glaciers would likely already be in an advanced position as seen from threshold lake data and were sensitive to moderate climate forcing.This glacier expansion matches the timing for the abrupt ice growth outlined for Icelandic and Arctic Canada glaciers between 675 -650 yrs.BP (1275 -1300 CE), followed by a substantial intensification 500 -495 yrs.BP (1430 -1455 CE) (Miller et al., 2012).The cold North Atlantic climate anomalies have been attributed to either internal climate variability (Drijfhout et al., 2013;Moreno-Chamarro et al., 2017) and/or external forcing from decadally-paced clusters of volcanic eruptions (e.g., Miller et al., 2012) and/or variability in solar activity (Miettinen et al., 2015;Sejrup et al., 2010;Moffa-Sanchez et al., 2014;Sha et al., 2016).
Northern hemisphere albedo can be enhanced through increased atmospheric sulphate aerosol concentrations, associated with increased volcanism, or increased spatial or temporal extent of snow and sea-ice area (Hartmann et al., 2013).Importantly, enhanced regional albedo can trigger regional, rather than global, cool periods.In the Arctic context, there is further suggestion that initial cooling from volcanic sulphate can produce feedback cooling from increased snow and sea ice extent (Miller et al., 2012).Similarly, radiative decline due to frequent volcanism in the period 500 -1000 CE (Toohey and Sigl, 2017) may explain the DACP advances) simultaneously with increased surface water-cooling and southwards expansion of sea ice.In summary, given our knowledge of northern hemisphere volcanic activity in the past four millennia years, as well as sea ice-albedo feedbacks, we regard an interplay between anomalous volcanism and sea-ice feedbacks as the most likely cause of the Older and Younger Neoglacials described here.
The three Neoglacial glacier expansions in Greenland roughly match the Iron Age Cold Periods, DACP and LIA within the European climate division.
We are aware that Northern hemisphere cold periods can be triggered by other different kinds of climate forcing than volcanism.Besides subtle trends in incoming solar radiation and greenhouse gasses fluctuations (Ruddiman et al., 2016;Flückiger et al., 2002).Changes in turbulent heat fluxes, however, can also produce persistent temperature anomalies.Turbulent heat fluxes describe the sensible and latent heat fluxes that are advected with either air or ocean currents.During the instrumental record period, the 1930s Arctic warm period has been attributed to anomalously high poleward heat transport (Overland and Wang, 2005).This is emerging evidence that the Arctic can be periodically more effectively insulated against poleward heat transport by either persistent strengthening of the polar vortex or persistent weakening of the Atlantic meridional overturning circulation (Gardner and Sharp, 2010;Wanamaker et al., 2012).It is possible that decadal-scale variability in poleward turbulent heat transport to the Arctic may have contributed to the initiation of the Older and Younger Neoglacials described here.Substantial uncertainty remains associated with understanding climatic drivers and responses of Late Holocene cool periods (Bader et al., 2020).

Reconciling the culmination phase proxies
The coastal records generally support the advance phase division with three culminations during the last 3000 yrs.The distribution is partly due to relatively few dates obtained from surface exposure dating and the low resolution provided by the archaeological sites.Oceanographic changes and glacier changes in the Vaigat Strait, north of Disko Island, show that more pronounced WGC cooling episodes occurred around 550, 250 and 50 yrs.BP (1400, 1700 and 1900 CE) and large peaks in IRD were observed concordant with a minor short-lived intra-LIA warming of the WGC (Andresen et al., 2011) at the time of the LIA maximum of the Jakobshavn marine-based glacier (Weidick and Bennike, 2007).This indicates a more dynamical LIA situation than predicted from the ice core records.
Shallow ice cores from the GrIS contain information on spatial variations across the ice sheet.Often a series of ice-core records have been stacked in order to infer local variability from the regional climate signal in these records (Vinther et al., 2010).2016) made a stacked record reaching back one millennium and covering a larger area of North Greenland, including twelve North Greenland Traverse ice core records and the NGRIP ice core record.Here, two separate sub-stacks of the North Greenland records; one stack for cores located in Northeast Greenland, shielded from the westerly winds, and another for cores located on the ice divide.Although these sub-stacks are subject to a 30year running mean, they demonstrate large regional differences and no consistent patterns such as anti-correlation or leads and lags.While the sub-stacks show a near anti-correlation during 950 -700 yrs.BP (1000 -1250 CE) they show a good correlation during 550 -40 yrs.BP (1400 -1990 CE).The North Greenland stack was compared to other Arctic records, including the Dye 3 ice-core record and points to a warm anomaly around 530 yrs.BP (1420 CE) that is seen in all northern records but not in the Dye 3. Regional-scale atmospheric variability can be identified by mapping the leading principal components of the isotope records and correlating individual ice cores records with temperature and pressure data (Vinther et al., 2010).Spatial variability in Greenland ice core records becomes evident when reviewing them individually with respect to two cold spells during the first half of the 19th century, which are both found in the PAGES2k temperature reconstruction (Neukom et al., 2019a) The LIA appears in many marine studies off shore Greenland as well as in the wider North Atlantic realm (Lapointe and Bradley, 2021) to be two-peaked or locally even three-peaked with intervening short-lived warming and often fluctuations within the cooler episodes.This could find support in the duplication of lateral moraines within the trimline zone, which is also recorded in written historical documents that shows glacier advance in the early 350 yrs.BP (18th century), at least in some West Greenland localities (Weidick et al., 2012).However, for the marine records, due to time resolution of the data and age uncertainties, it would be premature to draw conclusions about regional synchronicity (Fig. 4).There is certainly a structure to the ice fluctuations during the LIA epoch in Greenland, which from the available data suggests ice advance at ca. 1,700 yrs.BP (1250 CE) and culmination at about 550 yrs.BP (1400 CE) and again in the 18th century.The LIA glacier response must be seen in context with the older Neoglacial, meaning that the GrIS was in an already advanced position before the younger Neoglacial -LIA cooling.As such, the MWP and perhaps the RWP should be considered as short-term punctuations within a generally cold period, starting ca.4,000 yrs.BP because of complex interaction between atmospheric and oceanographic changes.Thus, our compilation does not support the general view of a delay in glacier response to the onset of the LIA (Dahl- Jensen et al., 1998;Weidick et al., 2012).The prominent fresh-looking landscapes in front of many outlet glaciers in Greenland are likely a result of several millennia of glacier activity and climate induced ice fluctuations and not a single LIA advance.Using additional geochronological approaches to augment 14 C determinations in this period e.g., 210 Pb and rock surface OSL are crucial to further constrain the chronology of LIA glacier oscillations.

Reconciling retreat phase proxies
As the LIA terminal area of most Greenland glaciers have been overridden several times by Neoglacial ice advances it is very difficult to recover evidence of multiple retreat phases except for the final period within the LIA.In summary, we find that generally the onset of ice retreat following the LIA and subsequent lichen growth begins between 540 to 65 yrs.BP (1410 to 1885 CE) in Southwest Greenland and 301 to − 15 yrs.BP (1649 to 1965 CE) in Southeast Greenland.Historical maps testify to 114 -104 yrs.BP (1840s-1850s CE) LIA maxima lasting until the present deglaciation commenced around 50 yrs.BP (1900 CE).If the sequence of lateral moraines seen at many outlet glaciers are interpreted rigorously it should indicate a successive lowering from the highest trimlines/lateral moraines generated by vertically and horizontally expanded glaciers to lower elevated trimlines/lateral moraines that record recessive advances or still-stands either over the entire LIA or from a final ice fluctuation i.e., early 19th century.The well-known presentday retreat pattern and surface lowering of major outlet glaciers from the GrIS is generally mirrored in the spatial distribution of the post LIA retreat after the 50 yrs.BP (1900 CE) advance following the warming in 1920s CE.We speculate that it is largely the same sensitivity pattern that governs the retreat from all three Neoglacial events and the minor retreat phases within the LIA.
Global mean sea level (GMSL) variability during the CE prior to the 20th Century very likely ranged between ~ ± 7 cm and ±11 cm in amplitude (Kopp et al., 2016).However, assessment of Greenland's contribution relies primarily on estimates from deglaciation models (e.g.HUY3 (Lecavalier et al., 2014) as extracting the source and magnitude of sea-level contributions as global proxy sea level indicators is inconclusive.Assessment of the GSML budget for the 20th Century in the IPCC AR5 included an unattributed residual for the period 1901 -1990 CE, however, lack of observation-based estimates meant that contributions from the polar ice sheets and ocean thermal expansion was not included (Church et al., 2013).Subsequent alternative statistical reanalysis of the GMSL from tide gauges suggested a smaller rate of sea-level rise during the 20th century than previously proposed, inferring the budget could be closed without invoking any additional contributions from ocean thermal expansion or ice sheets (Hay et al., 2015).
For the GrIS, however, reassessment of historical expedition maps, historic imagery, and landscape geomorphology (e.g. this review, Bjørk et al., 2012;Csatho et al., 2008;Kjeldsen et al., 2015;Lea et al., 2014;Weidick et al., 1996) clearly shows ice marginal retreat.Quantifying the historical mass balance changes has been derived using different methodologies.Box and Colgan (2013) used a surface mass balance reconstruction from 1840 onwards, based on a spatial and temporal statistical assimilation of meteorological station records, ice cores, and a regional climate model output.They modelled iceberg discharge as a function of runoff, to estimate a net eustatic sea level contribution of 25 ± 10 mm SLE during 1840-2010 CE (Box and Colgan, 2013).This  reconstruction estimated a net positive mass balance equivalent to ~1.2 mm SLE during 1840 -1900 CE.This suggests that up until the beginning of the 20th century, the ice sheet was near balance, which agrees with other independent observations.Kjeldsen et al. (2015) used 1900 CE as a Greenland-wide average for when ice started to retreat from its LIA maximum position and used observations from stereo-photogrammetric imagery to provide spatially resolved, observation-based geodetic mass balance estimates for three periods: LIAmax (1900CE -1983CE) 75.1 ± 29.4 Gt/yrs., 1983-2003CE 73.8 ± 40.5 Gt/yrs., and 2003-2010 CE 186.4 ± 18.9 Gt/yrs., leading to a total contribution of 25.0 ± 9.4 mm SLE, while estimates from an updated semi-empirical model, similar to that of Box and Colgan (2013), yielded 28.0 ± 23.8 mm SLE.These semi-empirical models that provide annual estimates on a Greenland-wide scale, suggest variability during the 20th Century.For example, considerable mass loss during the 1920s and 1930s and near balance conditions during the 1970s and 1980s, through lack of spatial information on where mass loss has occurred.Conversely, the photogrammetric-based estimates lack temporal resolution.However, the release of historic satellite imagery, combined with modern high-resolution spatio-temporal data, has enabled mass balance reconstruction at unprecedented scale from 1972 onwards.It is apparent that the GrIS contributed 13.7 ± 1.1 mm to GMSL rise from 1972 to 2018, with the majority of loss occurring after the 1990s following almost two decades of minor gain in the 1970s and loss in the 1980s (Mouginot et al., 2019).The latter increase in mass loss during the late part of the 1990s and continual trend into the 21st century is supported by numerous other studies (e.g.Khan et al., 2015;Shepherd et al., 2012).

Conclusions
• Across Greenland, we find a consensus in proxy data integrating advance and culmination phases of outlet areas of the GrIS for three distinct periods of glacier expansion, two within the older Neoglacial 2,500 -1,700 yrs.BP and 1,250 -950 yrs.BP and one within the younger Neoglacial 700 -50 yrs.BP, which includes the LIA cooling and glacier expansion at ca. 1250 CE. • We find that the LIA response must be seen within the wider context of the entire Neoglacial period over the past 5,000 yrs.Many glaciers advanced early in the Neoglacial and persisted in a close to presentday position and reacted to moderate climate forcing thereafter.• Ice expansion appears to be closely linked with changes in ice sheet elevation, accumulation and temperature as well as surface water-cooling in the surrounding oceans.At least for the two younger Neoglacial advances, volcanic forcing triggering a sea ice /ocean feedback could explain their initiation.• There is a centennial structure to the LIA ice fluctuations and available data suggests ice advance began around 700 yrs.BP (1250 CE), with culminations in the 15th century, early to mid-18th century and early mid-19th century.This implies that there is no delay in glacier response to the onset of regional cooling across most of Greenland and thus there is no lag in response time to the LIA cooling.However, the exact timing of individual fluctuations remains a problem due to poor chronological resolution, not least for the culmination phase.• Evidence for retreat phases during the LIA is largely restricted to the final retreat phase from the last LIA advance due to the transgressive nature of the outlet glaciers over the same landscape.A clear geomorphological imprint on the landscape shows a gradual lowering of the ice surface generating a consecutive series of lateral moraines associated with several recessive end moraines.The GrIS has a distinct sensitivity pattern of retreat in Northwest and Southeast Greenland as well as for the large outlets of Jakobshavn Isbrae, KNS, Frederikshåb Isblink, Qassimiut Lobe, Kangerlussuaq Glacier and Helheim Glacier.We speculate that a similar pattern was acting during retreat within the older Neoglacial retreat phases.• So far, the lack of spatial resolution of LIA data around Greenland hampers a clear view on regional synchronicity between the different advance and culmination phases.• Increased efforts to employ additional geochronological approaches other than 14 C determinations in these periods (e.g., 210 Pb, OSL surface dating) are crucial to further constrain glacier LIA advances and retreat.• Estimates of the rate of contribution from Greenland to global mean sea-level rise since the end of the LIA requires better understanding of the timing of the final culmination and retreat phase on a regional and local scale around Greenland.This study highlights the lack of data currently able to infer regional synchronicity of the final culmination and retreat phases.Until these data gaps are filled, we will not be able to reconstruct pre-satellite era, post-LIA rate mass loss with sufficient certainty.
interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 4 .
Fig. 4. (A) Locations mentioned in the text.Overview of major oceanographic changes during (B) the pre-cursor (500-1000 CE) cold period and (C) The onset of LIA.EGC= East Greenland Current.IC= Irminger Current.See references in text.

Fig. 5 .Fig. 6 .
Fig. 5. (A) Dating approaches to late Holocene glacier sediments and terminal moraines.Methods are related to i. an advance phase succeeding maximum recession and change to positive net mass balance, ii.A culmination phase mirroring maximum advance and a change to negative mass balance, iii.A retreat phase towards maximum recession.Modified from Porter (1979).Lc (Large catchment), Mc (medium catchment), Sc (small catchment).Grey field around the curve in A marks the different phases.(B) Three possible scenarios for ice fluctuation at the Greenland ice sheet margin since HTM.Scenario 1: No older Neoglacial advance, only a LIA advance in the younger Neoglacial.Ice margin near or inside the present margin succeeding the HTM retreat.Scenario 2: Older Neoglacial advance, but smaller than the LIA advance during the younger Neoglacial.Scenario 3: Older Neoglacial advance larger than LIA advance in younger Neoglacial.

Fig. 15 .
Fig. 15.(A) Locality map with archaeological sites, ages and distance to glacier front.(B) Kayak supports on the Nuulliit settlement in front og the Nuulliit glacier, Northern Greenland photographed by E.Holtved  (1936)  in(Holtved, 1954) and Jens FogJensen (2015).Note the heavily diminished glacier tongue.See Tables3 and 4. (C) Principal sketch of human settlements on the island of Pullassuaq and glacier front of Sermersuaq (Humboldt glacier) in Kane Basin.Pullassuaq holds a large number of Paleo-Inuit settlements assumed to date from Saqqaq to Late Dorset (4,450-500 yrs.BP) as well as a few Thule culture settlements (700-225 yrs.BP).Along the southwestern and southern side of the island there are moraine ridges indicating that at some point of time, probably at the maximum LIA advance, the glacier front reached the island, but did not override it.The dwelling ruins were left undisturbed.When Pullassuaq was settled in times of the Paleo-Inuit, the glacier front most likely had retreated even further than its position until 1990'ies when it was approximately one km from the islet.

Fig. 17 .
Fig. 17.Common lichen species used for lichenometry in Greenland.(A) Physcia caesia (Hoffm.)Fürnr.The more or less circular, pale grey thalli often form composite thalli as seen on the figure, where the biggest thallus measures c. 5 cm diam.Sermilik, Southeast Greenland.(B) Umbilicaria virginis Schaer.The species has up to 8 cm broad, often circular, grey; reticulately ridged thall, which often form colonies.The thalli shown on the figure are smaller, up to c. 1 cm broad.Sermilik, Southeast Greenland.(C).Rusavskia elegans (Link) S. Y. Kondr.& Kärnefelt (Syn.Xanthoria elegans (Link) Th.Fr.).The species has up to 5 cm broad, often circular, reddish orange thalli, which can fuse together with other thalli.Zackenberg, North East Greenland The diameter of the biggest thallus on the figure is c. 5 cm.(D).Pseudephebe minuscula (Nyl.ex Arnold) Brodo & D. Hawksworth.The up to 3 cm broad (biggest thallus on the figure), black-brown thalli are like the well-known green Rhizocarpon geographicum (L.) DC. one of the most often used lichens used in lichenometric measurements in Greenland.The species forms more or less circular thalli.Sermilik, Southeast Greenland.

Fig. 19 .
Fig. 19.(A) Retreat from the LIA maximum relative to 1978/87 CE.South Greenland (B) surface lowering from the LIA maximum relative to 1983 CE. N=308 outlet glaciers.
. Masson-Delmotte et al. (2015) found that these two cold periods, dated to 135 -125 yrs.BP (1815 -1825 CE) and 114 -104 yrs BP (1836 -1846 CE), respectively, have different strengths in the individual ice core δ18O records, and that the two spells exhibit different spatial patterns over Greenland.The NEEM record (Masson-Delmotte et al., 2015) as well as the principal components of winter δ 18 O records (Vinther et al., 2010; Ortega et al., 2014) show an intense cooling during the younger of the two cold periods, while cooling during the first cold period is less pronounced.However, ice cores from West and Northeast Greenland show a more intense cooling between 135 -125 yrs.BP (1815 -1825 CE) and a less pronounced or even insignificant cooling during 114 -104 yrs.BP (1836 -1846 CE).It follows that the 135 -125 yrs.BP (1815-1825 CE) cooling is strongest in the North Central part of Greenland and between 114 -104 yrs.BP (1836-1846 CE) the strongest δ 18 O anomalies are encountered along a Northwest-Southeast central Greenland transect (Masson-Delmotte et al., 2015).

Fig. 21 .
Fig. 21.Moraine degradation.Terminal moraine at Tunu glacier, Sermilik fjord, Southeast Greenland.Image recorded in 1933 CE (landform in middle part of the image).(B) same location in 2011.Modern glacier front is 1400 wide and distance between the 1933 and 2011 moraine is ca.1.2 km.Note the fluvial erosion of the terminal moraine.(C) LIA end-moraine degradation at Kangerlussuaq, associated with ice-covered moraines, some 30 m high.

Fig. 22 .
Fig. 22. Historical documentation.Printing year of geographical maps and ice extent from maps relative to the mapped LIA maximum extent.

Fig. 23 .
Fig. 23.Summary diagram.(A) Monte Carlo inverted borehole temperature records from GRIP and Dye 3 (Dahl-Jensen et al., 1998) (B) Scoop Lake, southern Greenland, maximum temperature anomaly inferred from chironomid δ18O (Lasher and Axford, 2019) (C) Interpretation of ice marginal response to 2nd derivative elevation changes in GRIP, NGRIP, Dye 3 and Camp Century ice cores (D) Advance, culmination and retreat phase proxy data (E) Summarising curve of broad relative changes in ocean surface temperature on the Greenland shelf.(F) Interpretation of the ice marginal position in Greenland.Red: LIA>older Neoglacial; Green: older Neoglacial advances > LIA advances; blue: Older Neoglacial at present-day position until LIA advances.
(continued on next page)

Table 1
(continued ) Fig. 7. Distribution of dated reworked material collected in LIA sediments and landforms.See Table1for detailed sample information.N=133.

Table 2
Synthesis of threshold lakes used in this study (n=48), of which 35 were included in analyses.Only lakes with minerogenic sediment suggesting glacial advance within their catchment in the last ~6000 years are included.

Table 3
Archaeological sites.Periods in bold are the more reliable

Table 4
Cultural division of the archaeological period in Greenland.

Table 6
Historical map data.