Late Glacial deglaciation of the Zackenberg area, NE Greenland

The Greenland Ice Sheet (GrIS) is a key component of the global climate system. However, our current understanding of the spatio-temporal oscillations and landscape transformation of the GrIS margins since the last gla-cialcycleisstillincomplete.TheobjectiveofthisworkistostudythedeglaciationoftheZackenbergValley(74°N, 20°E), NE Greenland, and the origin of the derivedlandforms.Basedon extensive ﬁ eldwork and high-detail geo-morphological mapping we identi ﬁ ed the different types of landforms, from which those of glacial and paraglacial origin were used to understand the processes driving regional environmental evolution. We applied cosmic-ray exposure (CRE) dating to 32 samples taken from erosive and depositional glacial landforms distributed across the valley. Geomorphological evidence shows that >800-m-thick Late Quaternary glacier ﬁ lled the valleys and fjords and covered mountain summits. In subsequent phases, as ice thickness decreased, the glacier was limited to the interior of the valley, leaving several lateral moraines. The deglaciation of the Zackenberg Valleythatstartedby~13.7 – 12.5kaalsoacceleratedslopeparaglacialprocesses.Manyblocksfromlateralmoraines wereremobilizedandfell,reachingthevalley ﬂ oorwheretheycoveredthethinningglaciertongue;transformingit into a debris-covered glacier that subsequently melted gradually. By ca. 10.5 ka, the last remnants of glacial ice dis-appearedfromtheZackenbergValley ﬂ oor,achronologyofdeglaciationthatissimilartothatobservedinothersites across NE Greenland.Theresultsofthisworkmust beconsideredinsimilarstudies,reinforcing theneedtosupport CRE ages of the different geomorphological phases with paleoclimatic data from other sedimentary records. © 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Polar regions are crucial components of the complex global climate system.Changes in these regions are not confined to the high latitudes, but have effects throughout the planet triggered by an intricate set of feedback processes between the atmosphere, ocean, sea ice, ice sheets, and land surfaces (Goosse et al., 2018).The Greenland Ice Sheet (GrIS) is considered a tipping element in Earth's climate (Lenton et al., 2008) because it has played a major role in the stability of the climate system since the last glacial maximum (LGM; 26-19 ka; Cohen and Gibbard, 2019).This mostly landbased ice sheet is the second largest body of ice globally, and the only ice sheet in the Northern Hemisphere, storing an ice mass with a sea level equivalents of ~7.4 m (Bamber et al., 2013).To predict future sea level rise and other consequences of accelerated GrIS melting, an accurate monitoring of modern ice sheet mass balance as well as a better comprehension of GrIS past dynamics are essential.
Over the last several decades, ice cores obtained from the interior of the GrIS, terrestrial climate records from ice-free areas, and marine sediment records collected from the adjacent sea floor have generated an accurate picture of the Late Quaternary climatic evolution and associated environmental changes (Briner et al., 2016).Together with Antarctica, the GrIS is one of the only ice sheets that survived the last deglaciation during Termination-1 (T-1; ~19-11 ka).The GrIS persisted during the Last Interglacial period (~130-116 ka) when global temperatures were significantly higher than present, although it was particularly reduced in its NE sector (Vasskog et al., 2015).The GrIS significantly expanded during the last glacial cycle (115-11.7ka), reaching its maximum volume during the LGM (~26.5-19ka) (Clark et al., 2009;Vasskog et al., 2015).Subsequently, as temperatures increased following the LGM, the GrIS shrank and became mostly land based.However, paleoclimate records reveal abrupt temperature shifts during T-1 (5-15 °C), with strong seasonality (Buizert et al., 2014;Vasskog et al., 2015) that must have driven changes in GrIS volume.Changes in the ice volume stored in Greenland influence freshwater delivery to source areas of North Atlantic deep water formation, and rapid changes may thus affect the stability of the climate system (Broecker, 2018).During T-1, these temperature shifts favoured the reorganisation of the thermohaline circulation in the Southern Ocean that led to rapid global CO 2 rise along with massive deglacial environmental and biotic changes in the polar regions (Denton et al., 2010).A better understanding of GrIS fluctuations within the last glacial cycle and during T-1 can thus provide insights about GrIS response to rapidly changing climate conditions at the present day.However, important knowledge gaps still exist as to how the GrIS and peripheral glaciers respond to climate switches, their sensitivity to climate, and spatio-temporal patterns of past glacial oscillations (e.g. Kelly and Lowell, 2009;Vasskog et al., 2015;Larocca et al., 2020aLarocca et al., , 2020b)).This is particularly true for areas including the NE Greenland coastal region (in which this research focuses), where the chronology of glacial fluctuations and landscape changes during T-1 are still poorly understood.
The modern ice-free land areas beyond the margin of the GrIS include sources of paleoenvironmental information (glacial records, lake sediments, fens, deltas, etc.) that can be used to reconstruct the evolution of the GrIS and the glaciers at its periphery since deglaciation.This is the case of the Zackenberg Valley, NE Greenland, where the glacial and periglacial geomorphology is well-known (Christiansen and Humlum, 1993;Christiansen, 1994Christiansen, , 1998;;Bennike et al., 2008;Cable et al., 2018), but the chronology of glacial fluctuations and associated paraglacial dynamics is yet to be determined.Here, the reconstruction of glacial history and landscape dynamics is hampered by (i) very active slope processes on the hillsides surrounding the main valley floor (solifluction, debris flows, nivation, etc.); (ii) the intensity of paraglacial processes following deglaciation; (iii) the transformation of the debris-free glacier into a debris-covered glacier during the last stages of glacial retreat and the intense paraglacial readjustment; and (iv), the timing of wastage of debris-covered ice, and how it affected the development of the present-day hummocky terrain (Christiansen and Humlum, 1993;Cable et al., 2018).These processes have all affected the stability of glacial landforms since their deposition and thus represent a challenge to the successful application of cosmic-ray exposure (CRE) dating, which has rarely been applied to glacial landforms in steep valleys and debris-covered glaciers in the (sub)polar regions (e.g.Tanarro et al., 2019;Fernández-Fernández et al., 2020;Charton et al., 2020).
In this study, we combined a detailed geomorphological survey with CRE dating of erosive and depositional landforms left by outlet glaciers in the Zackenberg Valley that calved into Young Sund fjord during the last glacial cycle (Christiansen and Humlum, 1993).Our goals were: (i) to examine the limits of CRE dating for establishing time constraints in highly dynamic glacial-paraglacial-periglacial environments; (ii) to reconstruct the spatio-temporal patterns of glacial culminations and retreats in such a geomorphologically active area; and (iii) to interpret the glacial chronology considering the complex geomorphological evolution of the valley.To achieve these goals, we specifically addressed the following questions: -To what extent can CRE dating be successfully applied in this area, where post-glacial erosive and depositional processes are widespread?-What information do the ages of glacial geomorphic features give about the environmental history of the Zackenberg area?Are they representative of the glacial chronology or have they been intensely reworked by paraglacial processes?-When did the GrIS margin start shrinking during the last glacial cycle?-What were the phases of major glacial advances/stillstands and retreats?-Is there a synchronous pattern of glacial advances/stillstands and retreats during T-1 in NE Greenland and other regions of Greenland?2. Regional setting

Study area
This study focuses on the ice-free Zackenberg Valley, situated in the Wollaston Foreland peninsula in the SE corner of the Northeast Greenland National Park, the world largest national park.Our study area encompasses the Zackenberg Valley floor from Young Sund fjord (south) to the mouth of the Store Sødal valley, the Dombjerg mountain (1492 m above sea level, hereafter referred to as a.s.l.) and Lindemansdal valley (north), as well as the surrounding peaks and slopes descending from Aucellabjerg (985 m a.s.l, east) and Zackenberg (1338 m a.s.l, west) mountains (Fig. 1).
The highest summits surrounding the Zackenberg Valley form relatively horizontal surfaces that are affected by intense periglacial processes.From these high plateaus, steep hillsides descend towards the valley floor on the northern and eastern slope of the Zackenberg mountain and the south-eastern slope of the Dombjerg mountain, with more gentle slopes to the west of the Aucellabjerg mountain.The Zackenberg River drains the 2-3 km wide U-shaped Zackenberg valley fed by snow and glacier meltwater, and forms a large delta where it reaches Young Sund.This flat area has also been affected by glacio-isostatic processes, which have given rise to a sequence of marine terraces largely covered by periglacial slope sediments.The Quaternary marine limit was established at ~40 m, with an intense postglacial crustal rebound during the Early Holocene, and a stabilization of the relative sea level ~3 ka (Christiansen et al., 2002;Pedersen et al., 2011).
The region has a polar tundra climate (Kottek et al., 2006).Between 1996 and 2015, the mean annual air temperature at Zackenberg Research Station was −9.0 °C and average annual precipitation of 367 mm (Højlund Pedersen, 2017), of which only 10% falls as rain during the summer months from June to September (Hasholt et al., 2008).The short summer season is crucial for the development of the valley's scarce vegetation cover.The large lowland areas include a moist to dry tundra dominated by shrubs <15 cm tall with grasslands, fens and interspersed snow patches, while the variety and size of the plants decrease at higher elevations (CAVM Team, 2003).The whole area is underlain by continuous permafrost (200 to 400 m thick) and a spatially variable active layer (45 to 80 cm thick) (Christiansen et al., 2008(Christiansen et al., , 2010)).
The area contains two different bedrock types separated by the Zackenberg and Lindermansdal valleys that extend along a large N-S fault zone (Escher and Watt, 1976).The western side consists of Caledonian crystalline complexes (Early Proterozoic) with abundant orthogneiss resistant to weathering processes, which provides higher and steeper slopes and rather coarse-grained deposits in the valley bottom.By contrast, the eastern fringe across the slopes of Aucellabjerg is composed of Cretaceous to Jurassic sedimentary rocks (mudstones, sandstones and conglomerates) (Henriksen et al., 2009), which results in a gentler relief with less inclined, long hillslopes with a sediment cover.Exposed bedrock is mostly found at the top of the mountains and upper slopes of both valley sides, whereas the lower parts and the valley floors are covered by Late Quaternary sediments (Gilbert et al., 2017;Cable et al., 2018).
The bedrock geology conditions the geomorphological processes prevailing today in the ice-free Zackenberg Valley, with very different landforms and processes on each side of the valley.Erosive and sedimentary glacial landforms are widespread across the valley and have been intensely transformed by paraglacial and periglacial dynamics.The eastern slopes descending from Zackenberg include a wide range of coarse debris landforms, such as talus cones, protalus lobes, stonebanked solifluction lobes and blockfields (Cable et al., 2018).The degree of preservation of these landforms is lower on the Aucellabjerg slopes along the eastern flank of the valley, where periglacial dynamics reshape the landscape more intensely through mass-wasting processes, such as solifluction or mudflows.Snow patches here are more abundant due to the local snow drifting, which results in copious meltwater availability during the summer season.The combination of high moisture supply, thin active layer and abundant fine-grained sediments favours widespread downslope mass movements (Cable et al., 2018).In some cases, debris flows and nivation activity also mobilize large amounts of sediment downslope and generate large alluvial fans extending down to the main valley floor.

The glacial geomorphology and previous chronological knowledge in the Zackenberg area
Glacial evidence in the Zackenberg area is distributed from mountain plateaus to valley floors and has been described in detail in previous studies (Christiansen and Humlum, 1993;Cable et al., 2018).However, despite the widespread, highly detailed geomorphological mapping of glacial features (moraines, erratic boulders and glacial polished surfaces) (Cable et al., 2018), the chronology of the different glacial phases is not yet constrained.
The existence of erratic boulders on the highest parts of Aucellabjerg and the Zackenberg mountains is indicative of a larger GrIS extent in the past, and of the minimum ice thickness during Quaternary glacial periods, when ice extended far onto the shelf (Bennike et al., 2008).In addition, the Zackenberg Valley includes unsorted glacial sediments (till) and several discrete landforms in the form of moraine complexes on the mountain slopes, and in the valley bottom.These landforms indicate phases with larger glaciers than today and other periods of glacial stabilization within the long-term deglaciation trend (Christiansen and Humlum, 1993).In addition, ice-moulded bedrock surfaces between the moraine ridges on the valley bottom also show clear traces of glacial abrasion, including striae in some cases, thus indicating warm/wetbased ice at some point in the past.
The likely oldest glacial evidence in the area are glacio-lacustrine deposits distributed on the lower western slope of Aucellabjerg (Fig. 3), dated by optically stimulated luminescence (OSL) between 84 ± 8 and 114 ± 11 ka, suggesting that during that time part of the Zackenberg Valley was ice-free and probably filled with an icedammed lake (Christiansen et al., 2002).Ice-free conditions were also reported from a pronival basin at 600 m in the Favorit Valley at Zackenberg Mountain, where a thermoluminescence date yielded an age of 66.2 ± 7 ka, thus confirming the absence of glaciers in the upper part of the valley during part of the last glacial cycle (Christiansen, 1994).OSL-dated glaciofluvial deposits indicate that the deglaciation of the lower Zackenberg Valley occurred after 22 ± 3 ka (Christiansen et al., 2002).Based on the distribution of glacial and periglacial landforms, Christiansen and Humlum (1993) proposed a tentative deglaciation chronology in several stages: (i) large outlet valley glaciers from the GrIS covering the Zackenberg area until 10-9.5 ka; (ii) inland retreat of glaciers and formation of a terminal moraine at ca. 9 ka; (iii) a period of readvance at ca. 8 ka as indicated by moraine ridges 1.5 km north of this moraine system.More recently, based on the fjordvalley fill, including sedimentary deltaic sequences and OSL ages, Gilbert et al. (2017) suggested that the Zackenberg lowlands may have been ice free as early as 13-11 ka, with the initial formation of the Zackenberg Delta in the lowest part of the valley taking place at ~10 ka after the major deglaciation of the valley.Christiansen et al. (2002) also reported 14 C dates suggesting that the development of the delta continued actively through the Early-Mid Holocene until 6.3 ka.These dates result from radiocarbon dating of organic fragments reflecting a minimum age after first establishment of soils or vegetation (including an unknown lag time behind ice-margin retreat) or fluvioglacial deposits located downstream from the glacier fronts, which must therefore be younger than the OSL dates.

Methodology
We assembled a geomorphological and geochronological approach in order to reconstruct the chronology of glacial oscillations of Zackenberg Valley and their interaction with paraglacial activity during the last stages of the valley's deglaciation.Fieldwork was carried out during late July and early August 2018, when the snow-free landscape allowed the identification of geomorphological features.

Geomorphological survey and mapping
Before fieldwork, we created a preliminary geomorphological map based on studies of satellite imagery, and adapting the existing cartography conducted by Cable et al. ( 2018) that identified the main glacial, periglacial, fluvial, and alluvial landforms in the Zackenberg area.Our main target in the geomorphological mapping was the spatial distribution of major glacial and paraglacial landforms, to understand the coupling of glacial, periglacial and paraglacial processes in newly exposed terrain and to optimize the sampling strategy for CRE dating.Considering this, we adapted the previous work with the focus on those landforms, mapping more precisely the location of the moraine ridges (hummocky and lateral), polished surfaces, and main erratic boulders (Fig. 3).Other landforms (as alluvial fans and fluvial deposits), which had less significance for our study, were simplified and broadly represented.Once in the field, we validated this map based on in situ observations that also allowed us to trace an overall picture of the deglaciation and the importance of paraglacial dynamics shaping the current landscape of the Zackenberg area.The final geomorphological map was drawn by digitizing landforms in the ArcMap 10.4.1 work environment over orthorectified panchromatic satellite WorldView-3 (0.3 m resolution) imagery from 2019.Our observations were also supported on the visual inspection of the shaded relief derived from the Digital Elevation Model (8 m resolution) provided by the GEUS (Geological Survey of Denmark and Greenland) and applied with transparency in our final map (Fig. 3).

CRE dating and field strategy
Based on the geomorphological map and field evidence, we collected 39 samples for CRE dating using 10 Be from sites judge to hold the best potential for glacial and paraglacial reconstruction, considering the complex coupling between deglaciation and paraglacial dynamics.Samples were obtained from depositional and erosional glacial landforms across the Zackenberg area (Fig. 2).They were acquired from three types of glacial surfaces: boulders from well-defined moraine ridges on mountain slopes and the valley floor (29 samples), scattered erratic boulders distributed across the landscape and well-inserted in the ground (6 samples), and exposed bedrock surfaces polished by glacial ice (4 samples) (Fig. 2).Field data and sample characteristics are listed in Table 1.Approximately 1 kg of bulk rock was extracted from each sample surface using a hammer and chisel.Samples were taken from a shallow section of the uppermost rock surface (<4 cm), which was mostly flattopped.We avoided sampling corners, edges, and steep faces (>20°) to ensure optimal reception of the cosmic ray flux.Moraine boulders were selected for sampling only if they were deeply anchored in the ground to ensure stability and minimize the risk of boulder collapse.Boulders were situated on the moraine crests and protruding such that they were unlikely to have been covered by sediments since their deposition or subject to unusually thick snow cover during the winter season, to minimize any related effects on nuclide production.The position (latitude, longitude) and elevation of samples were determined using a hand-held GPS (vertical accuracy ±10 m).Topographic shielding from the surrounding horizon was measured in situ with compass and clinometer.

Laboratory analytical procedures and exposure age calculation
Sample crushing and sieving were carried out in the "Laboratory of Physical Geography" at the Universidad Complutense de Madrid, Spain, to isolate the 0.25-1 mm fraction.Further physical and chemical sample preparation was carried out at the "Laboratoire National des Nucleides Cosmogéniques" (LN 2 C) of the Centre Européen de Recherche et d'Enseignement des Géosciences de l'Environnement (CEREGE, Aix-en-Provence, France).First, magnetic grains were discarded using a Frantz LB-2 magnetic separator.Quartz was concentrated in the non-magnetic fraction through repeated acid leaching with a 1:2 mixture of concentrated hydrochloric (HCl) and hexafluorosilicic (H 2 SiF 6 ) acids.Potential contamination by atmospheric 10 Be was removed from the separated quartz by performing 4 rounds of partial dissolution with concentrated hydrofluoric acid (HF).This process also removed the remaining impurities not dissolved in the previous steps.
After adding ~150 μL of a 9 Be carrier solution to each sample ("spike", concentration: 3025 ± 9 μg g −1 ; Merchel et al., 2008), the quartz samples were completely dissolved in 48% HF solution (3.6 mL g −1 quartz + 30 mL in excess).The solution was evaporated until dryness on hot plates, and the dry residues were recovered with hydrochloric acid.Following that, beryllium was precipitated to beryllium hydroxide (Be(OH) 2 ) with ammonia (NH 3 ) until the supernatant solution pH was 8. Afterward, samples underwent successive separations through anion and cation exchange resin columns to remove first iron and other metals and then boron (Dowex 1 × 8 and 50WX8, respectively) (Merchel and Herpers, 1999).The resulting eluted Be was precipitated to Be(OH) 2 by adding several drops of ammonia, and after being washed and dissolved again in concentrated HNO 3 , the solution was evaporated in a quartz crucible, and the residue was oxidized to BeO at 700 °C.Finally, the Be targets were prepared as a 1:1 mixture of niobium (Nb) powder and BeO for Accelerator Mass Spectrometer (AMS) measurements.
The measurements of the 10 Be/ 9 Be ratios on the BeO targets were performed at the French 5 MV "Accélérateur pour les Sciences de la Terre, Environnement et Risques" (ASTER) national facility at CEREGE using ion source 2, and 10 Be concentration in the samples was inferred from them.The 10 Be measurements were calibrated against the in-house standard "STD-11", with an assigned 10 Be/ 9 Be ratio of (1.191 ± 0.013) × 10 −11 (Braucher et al., 2015).A 10 Be half-life of (1.387 ± 0.0012) × 10 6 years was used (Chmeleff et al., 2010;Korschinek et al., 2010).The analytical 1σ uncertainties in the 10 Be/ 9 Be ratios include those in AMS counting statistics, the standard 10 Be/ 9 Be ratio and an external AMS error of 0.5% (Arnold et al., 2010).The analytical 1σ uncertainties in the inferred 10 Be concentrations include the propagation of the chemical blank correction.Some samples yielded limited 10 Be counting statistics and thus low precisions (~6-8% or more) due to unstable and/or low 9 Be currents (<1 μA), which can arise either from impurities in the BeO targets or small sizes of the BeO targets.Most samples with low-current AMS 10 Be/ 9 Be measurements yielded relatively high uncertainties, and were therefore considered chemical outliers (i.e., .Such samples were discarded for further exposure age calculations and discussion.For ZAC-C, while bad current values were recorded during the AMS measurement, the sample still provided a good accuracy measurement, and it was therefore retained for the geochronological discussion.All analytical data is shown in Table 2. Exposure ages were calculated with version 3.0 of the CRONUS-Earth online calculator (Balco et al., 2008; https://hess.ess.washington.edu/).After a bibliographic review of 10 Be-based glacial chronologies in the nearby areas, we applied the same production rates and parameters used in recent publications in order to unify criteria and make our exposure age calculations comparable.Therefore, we used the Arcticwide sea-level/high-latitude 10 Be production rate (3.96 ± 0.15 atoms g −1 a −1 ) (Young et al., 2013) and the "Lm" (Lal/Stone) timedependent scaling model (Lal, 1991;Stone, 2000) (Table 2).For all samples, a 2.7 g cm −3 density was assumed, and no corrections of erosion or snow shielding were applied, and thus, exposure ages are reported for zero-erosion and snow-free scenarios.The partial shielding effect of the surrounding topography was corrected for all sampling sites using the "Topographic Shielding Calculator v.2" (http://stoneage.ice-d.org/math/skyline/skyline_in.html).Post-glacial glacio-isostatic rebound (40 m in the last 9.5 kyr; (Christiansen et al., 2002)) was not considered in exposure age calculations, as the resulting age offset represents only ~7% (older) in Greenland (Jones et al., 2019), which does not affect our conclusions.
Exposure ages are given in Table 2 with their internal (only analytical) and external 1σ uncertainties (including production rate uncertainty).In the text and figures, ages are given with their internal uncertainties unless otherwise stated.A chi-squared test (following Ward and Wilson, 1978) considering the analytical uncertainties was applied to different sample populations (e.g.within single moraine ridges, debris-covered glacier deposits, etc.) to identify potential outliers based on statistical criteria.Age outliers (i.e.too old and too young ages; Heyman et al., 2011) may arise from nuclide concentrations inherited from previous exposure periods or post-depositional processes (erosion, exhumation, rock falls etc.).When age outliers were detected, such ages were excluded from the mean age calculations for the corresponding geomorphological units (e.g., moraine ridges).Although these outliers cannot be used for paleoclimate or paleoglacial inferences, they are indicative of the complexity of using CRE dating in these highly dynamic glacial-paraglacial-periglacial environments.The mean ages were calculated arithmetically, and their uncertainties include the standard deviations of the single ages (not excluded) and the squared production rate uncertainties.

Results
The distribution of geomorphological landforms across the landscape is indicative of the different phases that occurred during the deglaciation in the Zackenberg Valley (Fig. 3), which are chronologically constrained by 10 Be CRE ages inferred from the glacial record (Table 2).

Geomorphological setting and sampling strategy
Based on the spatial distribution of glacial, paraglacial and periglacial landforms, we identified three main sectors in the study area:

Mountain plateaus
The bedrock is mostly exposed above 400-600 m a.s.l. on the mountain plateaus.Glacio-nival cirques, nivation hollows and perennial snow patches are also found near the mountain tops and upper slopes.The summit surfaces (e.g., Aucellabjerg and Zackenberg mountain plateaus) include geomorphic evidence of past glaciations in the form of large erratic granite boulders distributed across the currently ice-free periglacial landscape.Two of these, well-anchored in the relatively flat area composed of sedimentary rocks of the highest slopes of Aucellabjerg, were sampled for CRE dating (ZAC-A and ZAC-C) in order to infer the onset of glacial thinning.The relatively rounded edges of these boulders (Fig. 4) provide evidence that they were transported for long distances, which diminish the potential of nuclide inheritance.

Mountain slopes
Discontinuous lateral moraine ridges were found on the slopes surrounding the Zackenberg Valley at varying elevations, recording the thickness and geometry of the glacier during past glacial phases (Fig. 3).The limit of the highest remnants of moraine deposits ranges ~200 and 550 m a.s.l. on the western slope of Aucellabjerg (Christiansen and Humlum, 1993), ~200 and 800 m a.s.l on the southeast slope of Dombjerg, and between ~400 and 200 m a.s.l. on the eastern slope of Zackenberg on the western side of the main valley.
Nine roughly preserved lateral moraines were identified on the western slopes of Aucellabjerg, counting up to twelve discontinuous ridges distributed at elevations between ca. 60 and 400 m a.s.l.Very active periglacial slope processes have dismantled most of the moraine ridges on this slope through widespread nivation activity with debris flows, mudflows and solifluction processes (Cable et al., 2018); remnants of these moraines are found only on mid-low slopes in the form of discontinuous crests.Whereas the moraines at higher elevations are mostly made of fine sediments, the lower, more recent ridges have a greater abundance of large conglomerate and sandstone boulders of these materials.We collected seven samples from orthogneiss moraine boulders distributed across the prevailing sandstone bedrock, from the upper (403 m a.s.l.) to the lowest sector (64 m a.s.l.) of the Aucellabjerg slope (ZAC-11, ZAC-12, ZAC-13, ZAC-14, ZAC-15, ZAC-16, and ZAC-17) to reconstruct past phases representative of stillstand positions of the glacier within its long-term thinning and retreat.
On the eastern slope of Zackenberg, intense frost shattering has favoured rockfall activity, with the development of protalus lobes and talus cones that have dismantled the moraines located at the foot of the rock wall and partially covered the exposed glacially polished bedrock surfaces.Up to four lateral moraines were identified in this area, which descend gradually towards the bottom of Zackenberg Valley (Fig. 3).Three samples were taken from intermediate moraine ridges in this unit, distributed between 80 and 120 m a.s.l.(ZAC-18, ZAC-19, and ZAC-20) (Fig. 3).
In the northern part of the Zackenberg Valley, several orthogneiss boulders are distributed on a heavily weathered rocky sandstone threshold that also constitutes the watershed divide between the Zackenberg and Lindemansdal valleys.Two of these erratic boulders (ZAC-05 and ZAC-06) and the bedrock surface (ZAC-07) were sampled (Fig. 3).
Moraines are specially well-preserved in the northern part of the Zackenberg Valley on the southeast slope of Dombjerg mountain, where eleven fragments of moraine ridges ranging from ~150 to 400 m a.s.l.gradually curve and slope towards the east, particularly the middle and highest ridges (Fig. 3).The intermediate and lower moraine ridges are discontinuous and highly eroded by alluvial processes, with dispersed till across the valley bottom.We collected nine samples from the upper (330 m a.s.l.) to the intermediate moraines (260 m a.s.l.) (ZAC-21 to ZAC-29) (Fig. 3).

Valley floor
The base of Zackenberg Valley contains two very different geomorphic settings according to the glacial evidence.The upper part of the valley includes widespread glacial landforms including a well-defined frontal moraine system with internal hummocky terrain derived from thermokarst processes and is the site of several lakes.We collected six samples from the frontal moraine (ZAC-39, ZAC-40, B1, B2, ZAC-01b, and ZAC-04) (Fig. 3) and four more from an internal moraine ridge (ZAC-30 to ZAC-33) that dams two of these lakes.The hummocky terrain only defines a clear frontal ridge occupying the valley floor, but there are no clear recessional moraines composed of transverse moraine ridges, as observed in the field.This area is composed of randomly distributed hummocks, mounds with lakes in between, and is dissected by inactive river channels.On the surface of the hummocky terrain, there are abundant angular to subangular boulders, with the presence of few slightly rounded blocks.Some exposed polished bedrock surfaces are located between the discontinuous lateral moraine remnants in the western side of the valley that have been removed by rock falls from the steep rock walls.These polished surfaces were sampled (ZAC-01 and ZAC-36), as were several erratic boulders distributed on them (ZAC-03, ZAC-37 and ZAC-38) (Fig. 3).The lowest section of the valley extends ca. 3 km inland from the present-day coastline and comprises relatively flat ground moraine areas shaped by glacio-isostatic processes, delta sedimentation and nivation activity (Cable et al., 2018;Christiansen, 1998).Here, several metre-size erratic boulders are scattered across the landscape, partially covered by aeolian, delta and glaciomarine sediments.

Exposure ages
The samples collected from the three principal geomorphological settings yielded ages ranging from 78.8 ± 2.8 (ZAC-A) to 9.4 ± 0.5 ka (ZAC-06).Here we present the results ordered from the mountain plateaus down to the valley bottoms: 1. Mountain plateaus Two large erratic boulders of granite found on the highest surfaces of Aucellabjerg (above 800 m a.s.l.) yielded 10 Be exposure ages of 78.8 ± 2.8 (ZAC-A) and 50.3 ± 2.6 ka (ZAC-C), respectively (Fig. 4, Table 2).
On the northern fringe of Zackenberg Valley, in the pass dividing the northern and southern sides of the Wollaston Foreland peninsula, the three samples (Fig. 6C), all around 200 m a.s.l, include two erratic boulders that yielded ages of 11.1 ± 0.5 (ZAC-05) and 9.4 ± 0.5 ka (ZAC-06) with a mean age of 10.3 ± 1.3 ka (n = 2), and another from polished bedrock that was dated at 10.2 ± 0.5 ka (ZAC-07).

Discussion
We complemented the previous detailed mapping of the glacial and periglacial landforms distributed across the Zackenberg Valley (Cable et al., 2018) with 10 Be CRE ages of samples collected from glacial features of the slopes and valley floor.This approach enabled us to generate a detailed space-time reconstruction of deglaciation in the Zackenberg region.

General considerations on the chronological process of deglaciation
The 32 10 Be CRE ages from moraines, erratic boulders and exposed polished bedrock surfaces provide a record of glacial oscillations from ca. 80 ka to the Early Holocene (11.7-8.2 ka).Analysis of the exposure ages and their geomorphological setting also reveals that some of them do not apparently follow a logical sequence.
The western slopes of Aucellabjerg (include abundant fine-grained sediments (mostly silts and sands) with frequent gravels), which enhances the remobilization of surface sediments by solifluction dynamics, affecting the stability of boulders and potentially causing them to roll or move (Balco, 2020).Despite sampling boulders that were wellanchored in the ground and formed part of clearly distinguishable moraine ridges (Fig. 5), it is possible that the very active post-depositional slope processes on Aucellabjerg, or even incomplete exposure due to exhumation, may have resulted in underestimated ages in some cases (Figs. 3 and 9).By contrast, older exposure ages than expected may result from nuclide inheritance or downslope mobilization from (previously exposed) upper parts of a slope.This is the case of samples ZAC-14 and ZAC-15, obtained from a moraine ridge in the middle slope of Aucellabjerg, that do not follow a consistent geochronological pattern, as they are much older than the rest of the highest moraines of the valley, and were thus discarded as outliers.However, once these obvious outliers were removed, there is a consistent geochronological sequence that logically matches the geomorphological observations, and from which the overall deglaciation history of this valley can be reconstructed.

Coupling CRE ages and geomorphological evidence
The highest samples, collected from the summit plateaus, returned the oldest ages.Granitic erratic boulders situated near the summit plateaus demonstrate that the minimum ice thickness in the Zackenberg Valley during their deposition was >800 m.As local plateaus are not composed of granites, we infer that the ice must have flowed from another source, the nearest of which is adjacent to the current margins of the GrIS, ca.60 km west of the Zackenberg Valley.According to the ages of the highest samples ZAC-A and ZAC-C, this phase may have occurred sometime between ca.80-50 ka (MIS 4-3).The occurrence of nuclide inheritance is a well-known issue in deglaciated areas throughout the Arctic (i.e.Greenland: Goehring et al., 2010;Søndergaard et al., 2019;Ceperley et al., 2020;the Canadian Arctic: Bierman et al., 1999;Davis et al., 1999;Briner et al., 2003;Kaplan and Miller, 2003), and should not be ruled out, especially in regions of high relief with gentle topography such as those of the study area, which in turn may determine low ice mobility, low erosion rates and consequently the preservation of an inventory of 10 Be accumulated in previous ice-free periods (see Fernández-Fernández et al., 2021).Future studies should combine multiple nuclides (e.g. 26 Al, in-situ 14 C) to explore complex exposure histories and shed light on processes such as exhumation and prior exposures.
However, in NE Greenland, erratic boulders transported tens of km from their sourcesuch as the granitic boulders distributed across Aucellabjerg and the Zackenberg mountainsare not generally affected by nuclide inheritance and have provided CRE ages associated with glacial retreat following the last glacial cycle (Håkansson et al., 2007a;Larsen et al., 2018).This is also valid for moraine boulders that do not usually retain an inheritance signal from past deglacial periods (Levy et al., 2016;Biette et al., 2020aBiette et al., , 2020b)), with few exceptions (Håkansson et al., 2009).A thermoluminescence measurement from a pronival basin at 600 m a.s.l. on the western slope of Aucellabjerg returned an age of 66.2 ± 7 ka, which may also indicate that the areas at higher elevations were already ice-free by that time (Christiansen, 1994).It is therefore reasonable to consider a scenario of pre-LGM glaciation in the valley, with most of the Zackenberg Valley ice-filled and only some nunataks protruding above the ice field.
The rest of the CRE ages confirm that the deglaciation of the slopes and valley floor occurred during T-1 (Fig. 7).The lateral moraines on the western slopes of Aucellabjerg suggest an age of deglaciation of 12.5 ± 1.5 ka, which is very similar to the exposure ages obtained from the eastern slope of Zackenberg Mountain (12.3 ± 0.8 ka).These ages overlap with the slightly older average age reported from highest ridges on the south-eastern slope of Dombjerg (13.7 ± 1.1 ka; Fig. 7), suggesting that the formation of the highest moraines in the Zackenberg Valley occurred between 13.7 and 12.5 ka.
Slightly younger ages were found on the northern side of the study area in the pass dividing the northern and southern sides of the Wollaston Forland peninsula.The boulder from which we obtained the older sample ZAC-05 (11.1 ± 0.5 ka) was more weathered than ZAC-06 (9.4 ± 0.5 ka).The average exposure age of both boulders (10.3 ± 1.3 ka) matched that of nearby polished bedrock (ZAC-07; 10.2 ± 0.5 ka).This indicates the end of the glacial transfluence from Store Sødal towards the north of the peninsula at that time, and thus the separation of the Zackenberg and Lindemansdal as independent glaciers within their respective valley floors, leaving an ice-free pass in between.
The moraine systems on the southeast slope of Dombjerg provided a complete sequence of glacial thinning from 13.7 ± 1.1 to 11.3 ± 1.1 ka (~330-~260 m a.s.l.).The higher and middle moraines sloping towards the east suggest that the ice was flowing from Store Sødal, whereas the thinner glacier that generated the lowest moraines was still connected with the Lindemansdal glacier until 10.2 ka.Indeed, the age of formation of these lower ridges fits well with an erratic boulder collected on a bedrock surface located just above the hummocky terrain (11.7 ± 1.0 ka) at the west bank of the river, providing further evidence for glacial retreat at the onset of the Holocene.
The frontal ridge revealed the culmination of the last glacial advance, and included rocks deposited between 11.9 ± 0.7 (ZAC-04) and 10.2 ± 0.5 ka (B1).Average exposure ages indicated that the frontal moraine system was deposited around 11.2 ± 1.1 ka, similar to the lowest moraine system on the southeast slope of Dombjerg (11.3 ± 1.1 ka).This suggests that units formed simultaneously, when the glacier generated a well-defined frontal moraine system extending across the valley floor and encircling the southeast slope of Dombjerg.
As the glacier thinned and retreated, the slopes surrounding the valley floor became ice-free, triggering paraglacial processes (Ballantyne, 2008;Oliva et al., 2019).We hypothesise that the debuttressing of the steep valley walls delivered large deposits, from lateral moraine boulders in most cases, onto the glacier surface.The terminal section of the glacier, already very thin, gradually transformed into a debris-covered glacier (Fig. 8) following a pattern also observed in other mountain systems (Janke et al., 2015;Anderson et al., 2018;Jones et al., 2019;Kenner, 2019;Mayr and Hagg, 2019).This transformation, from debris-free to debris-covered glacier, reduces ablation and even accelerates glacier flow (Hambrey et al., 2008;Deline et al., 2015;Anderson and Anderson, 2016;Anderson et al., 2018;Mayr and Hagg, 2019).Once the glacier stagnated and began to retreat, ice disintegration favoured the existence of dead ice patches, which can persist for millennia particularly in permafrost environments (Fernández-Fernández et al., 2017).This process can affect the stability of glacial sediments and trigger the readjustment of some boulders, resulting in younger exposure ages (Fernández-Fernández et al., 2020).Other boulders, which would have remained stable for millennia on the glacial surface, were continuously exposed to cosmogenic radiation that resulted in ages pre-dating the stabilization of the deposit (Bibby et al., 2016;Mackay and Marchant, 2016;Amschwand et al., 2020).This may explain the average age of the hummocky terrain samples (12.6 ± 2.2 ka; Fig. 3).
The study of the uncertainties associated with CRE dating in debriscovered glaciers is still incipient with no conclusive evidence, and these uncertainties are highly dependent on local geomorphological settings, as has been observed in mountain environments in the Himalaya (Scherler and Egholm, 2020), central Alaska Range (Dortch et al., 2010) and Iberian Range in Spain (Fernández-Fernández et al., 2017), as well as (sub)polar regions such as in the Kerguelen Archipelago (Charton et al., 2020), Iceland (Tanarro et al., 2021) or in the McMurdo Dry Valleys in Antarctica (Mackay and Marchant, 2016).CRE ages obtained from fossil debris-covered glaciers may not represent the exact time of deglaciation at a specific site, but they may provide wide timescales ranging between the fall of boulders on the glacier surface and the stabilization after complete melting of the buried ice (Fernández-Fernández et al., 2017).Indeed, ice disapperance may take hundreds or even thousands of years in (almost) static glaciers due to the insulating protection of a thick debris cover (Mackay and Marchant, 2016), and does not necessarily follow a specific spatial pattern, which must be taken into account when interpreting CRE ages from debris-covered glaciers.
Our CRE dates cannot provide further chronological control for climatic changes that led to the formation of the moraines, the withdrawal of the glacier and the stabilization of the hummocky terrain, although the overlapping of the CRE ages within uncertainty ranges suggests that it must have been a rapid process.The application of CRE dating to collapsed debris-covered glaciers must therefore take into account a wide range of possibilities (Charton et al., 2020): (i) fallen boulders from the lateral moraines that may have retained nuclide inheritance (Çiner et al., 2017;Dede et al., 2017) (ii) boulders that overturned during massive melting (Fernández-Fernández et al., 2017); and (iii) boulders that have remained static on the glacier surface and have been continuously exposed to cosmic radiation since their emplacement (Bibby et al., 2016;Mackay and Marchant, 2016;Winkler and Lambiel, 2018;Tanarro et al., 2019;Amschwand et al., 2020;Charton et al., 2020;Fernández-Fernández et al., 2020).Given this diversity of possibilities, it is therefore difficult for a debris-covered glacier to yield rather homogeneous CRE ages after the ice melts.The range of boulder ages within the same geomorphological units must therefore be interpreted in light of the competing processes that control the disintegration of a debris-covered glacier, and its transformation into hummocky terrain.The absence of continuous and clearly aligned ridge crests within the hummocky terrain, together with numerous lakes occupying inter-hummock depressions, supports the interpretation that this hilly terrain formed as a result of stagnation and passive melt-out and not as a consequence of active ice margin retreat (Benn, 1992;Palacios et al., 2021;Rodríguez-Mena et al., 2021;Rodríguez-Mena et al., 2021).
Lastly, the three samples collected from two erratic boulders and a polished surface above the hummocky terrain provided similar results of 10.4 ± 0.7 and 10.7 ± 0.5 ka, respectively.These ages suggest that ice persisted at the foot of the northern slope of Zackenberg almost two more millennia after the valley floor became ice-free (Fig. 3).The most plausible hypothesis is that a slope glacier, fed by the snowdrift and avalanches, persisted at the foot of Zackenberg's shaded northeastern slope at the onset of the Holocene, capable of bedrock erosion on this weathered bedrock.Boulders detached from the rock wall probably slid above the ice mass and accumulated at the base of the slope on the dated exposed polished bedrock when ice had shrunk upslope.

Deglaciation chronology of the Zackenberg area
Our results provide evidence that the deglaciation of the Zackenberg Valley occurred during the Late Glacial, with the following phases in the environmental evolution of the area: -Phase 1: Maximum expansion of the glaciers, covering the summits.
Our highest samples indicate that Zackenberg Valley was almost fully filled, likely prior to the LGM of the last glacial cycle at ca. 80-50 ka, and only the highest peaks may have been ice-free (as no erratic boulders have been found across the plateaus) as nunataks that were affected by intense periglacial conditions (Christiansen et al., 2002; Fig. 8A).At that time, the Zackenberg glacier was one of the outlets of the GrIS, and merged with the Tyrolerfjord glacier, another outlet of the GrIS, flowing east towards Young Sund (Christiansen and Humlum, 1993).However, previous research has provided evidence that the GrIS also advanced to the shelf edge during the LGM and probably covered some of the high plateaus in NE Greenland (Håkansson et al., 2007b;Skov et al., 2020).We cannot fully discard that our dated samples may retain some nuclide inheritance, masking a last glacial cycle glacial advance covering most of the Zackenberg Valley.
-Phase 2: Onset of deglaciation of the debris-free glacier.
In NE Greenland, a rapid deglaciation process began after the deglaciation after the global LGM (Clark et al., 2009).The loss of glacier thickness in the Zackenberg Valley occurred parallel to the retreat of the Tyrolerfjord glacier (Fig. 1).As a result, a sequence of recessional moraines formed at different elevations.No robust evidence of LGM or early T-1 moraine deposition has been found, although it is possible that the highest ridges that may have formed during the LGM advance, and the remnants of some moraines were subsequently dismantled due to the very intense slope processes.The long-term recession of the GrIS and surrounding glaciers was interrupted by periods of glacial stillstand/advance of the Zackenberg valley glacier that favoured the formation of the moraine ridges on the slopes surrounding the main valley floor.The highest moraine ridges, dated at 13.7-12.5ka, may be associated with phases of glacial stabilization during the overall retreat between the Bølling-Allerød (B-A) Interstadial (14.6-12.9ka; GI-1 Greenland ice cores; Rasmussen et al., 2014) and the early Younger Dryas (YD) stadial (12.9-11.7 ka; GS-1; Rasmussen et al., 2014).The lowest ridges (11.3 ka), deposited concurrently with the outer moraine ridge enclosing the dated hummocky terrain (11.2 ka), formed at the onset of the Holocene (Fig. 9).
The glacier stagnated and thinned at this time, triggering the intensification of paraglacial dynamics on the recently deglaciated slopes.The abundant sediment delivery from lateral moraines and recently ice-free rock walls onto the glacier surface favoured its transformation into a debris-covered glacier during the Early Holocene at ca. 11-10 ka.
-Phase 4: Degradation and final collapse of the debris-covered glacier.
As temperatures continued to rise during the Early Holocene, the debris-covered glacier underwent an irregular collapse, and its deposits became a hummocky terrain.The glacier disappeared from the Zackenberg Valley by 10.3 ± 1.3 ka, exposing polished bedrock outcrops and abandoning erratic boulders as the glacier terminus receded up the Store Sødal valley.The valley has not had glacial ice during the last 10 ka as no younger glacial landforms have been identified.As a result of the deglaciation, glacio-isostatic processes exposed a sequence of marine terraces and favoured the incision of the Zackenberg River.A rapid delta progradation began ca.13-11 ka and continued until 6.3 ka, accompanied by permafrost aggradation (Gilbert et al., 2017;Christiansen et al., 2002).

Glacial oscillations in NE Greenland since the last glacial cycle
To date, the available glacial reconstructions have focused mostly on the Northeast Greenland Ice Stream (NEGIS), which drains ~12% of the interior GrIS via three marine-terminating outlet glaciers (Larsen et al., 2018), as well as coastal areas of central East Greenland surrounding Scoresby Sund, ~400 km south of Zackenberg (Kelly et al., 2008;Håkansson et al., 2009Håkansson et al., , 2011;;Lowell et al., 2013;Levy et al., 2014Levy et al., , 2016) ) (Fig. 10).Our 10 Be CRE age dataset, focusing on glacial landforms in the Zackenberg area, expands and strengthens the previous glaciation history based on radiocarbon and luminescence dates (Christiansen and Humlum, 1993;Christiansen, 1994Christiansen, , 1998;;Christiansen et al., 2002;Bennike et al., 2008;Gilbert et al., 2017), providing new evidence of the major role that the Late Glacial deglaciation exerted in the environmental evolution of NE Greenland.
Glacial evidence around Greenland shows asynchronous glacial fluctuations of the GrIS and peripheral mountain glaciers during the last glacial cycle and subsequent deglaciation (Kelly and Lowell, 2009;Vasskog et al., 2015).There is a gap of terrestrial empirical evidence of glacial fluctuations prior to the LGM, although available offshore marine records indicate that all GrIS regions underwent significant growth and retreat in millennial-scale phases (Funder and Hansen, 1996;Lecavalier et al., 2014;Vasskog et al., 2015).
To date, it is unclear whether the interfjord plateaus remained icefree throughout some periods within the last glacial cycle when the fjords were occupied by outlet glaciers (Funder et al., 1994), or on the contrary, were not ice-free at any time until the post-LGM ice thinning (Håkansson et al., 2011).In West Greenland, while some CRE ages of erratic boulders indicated that coastal mountain tops (above ca.800 m a.s.l.) have been ice free since ca.130 ka, (and thus were not glaciated during the last glacial cycle; (Roberts et al., 2009;Lane et al., 2014;Roberts et al., 2013), more recent studies suggest that these high-elevation surfaces were ice-free almost ca.90 ka before the LGM and ice-covered during this period (Strunk et al., 2017).Assuming that the highest dated erratic boulders in the Zackenberg Valley do not retain nuclide inheritance, the maximum glacial expansion of the last glacial cycle would have occurred at ca. 80-50 ka, when most of the valley was inundated by ice but the mountain tops remained ice-free.At this time, the glacier front was probably located ca. 30 km from Zackenberg on the outer shelf (Christiansen and Humlum, 1993), as reported in several other areas across North and central East Greenland (Funder et al., 2011;Lecavalier et al., 2014).An age of 79.1 ± 3.1 ka was reported by Håkansson et al. (2007b) from a moraine boulder in the Store Koldewey island, 250 km north of Zackenberg, although the authors cautioned that it might have been subject to nuclide inheritance.The longest glacial chronology of the last glacial cycle in NE Greenland is from the NEGIS, which was smaller than present between ~41-26 ka, with glacier front at least 20-70 km behind the present ice margin (Larsen et al., 2018).
However, if we assume that nuclide inheritance occurred in the glacial erratic boulders from the Aucellabjerg plateau it is reasonable to suggest that the glacial advance which occurred during the LGM covered most of the Zackenberg Valley.According to the GrIS reconstructions, the ice sheet occupied an area 65% larger than present during the LGM, with margins generally reaching at least the continental shelf (Kelly and Antony, 2009;Funder et al., 2011;Vasskog et al., 2015).This is confirmed in the northeast sector, where NEGIS significantly expanded by 26 ka (Larsen et al., 2018).Evidence for the LGM glacial advance also exists further south near the Scoresby Sund mouth, where 10 Be ages of erratic boulders placed the ice margin into the outer shelf until 17.3 ka, with a thickness of ca.250 m (Håkansson et al., 2007a).Despite the lack of geomorphic evidence of the LGM glacial advance in Zackenberg Valley, geophysical modeling carried out in the Wollaston Forland region suggested ice between 500 and 1000 mthick (Fleming and Lambeck, 2004), which was supported by glacial trimlines on the mountain sides indicating that the Zackenberg area was occupied by valley glaciers (Bennike et al., 2008).The highest moraine remnants in the Zackenberg Valley, located at elevations between 500 and 800 m, may thus have formed during this glacial advance (Fig. 10).
Following the LGM, temperatures increased at the onset of T-1 and a more pronounced warming took place ca.17 ka (Kobashi et al., 2017), resulting in a significant reduction of GrIS volume and retreat of its margins as modelled by palaeoglaciological studies (Funder et al., 2011;Vasskog et al., 2015).It is still unclear whether GrIS deglaciation persisted without interruption well into the Holocene, or if this longterm retreat was interspersed with periods of glacial advance/stabilization that led to phases of moraine formation (Vasskog et al., 2015).
Zackenberg data, however, show that the long-term GrIS recession was interrupted by short periods of stillstand or glacial advance.Recent studies have confirmed the occurrence of several glacial oscillations within T-1 across Greenland that led to the development of moraines (Young et al., 2020), confirming previous evidence inferred from marine sediments (Ó Cofaigh et al., 2013).In NE Greenland, few studies have provided CRE ages of deglaciation within T-1.Biette et al. (2020aBiette et al. ( , 2020b) dated a culmination of moraine formation at 16.2 ka in the Clavering Island, ca. 30 km south from our study area.This age is significantly older than the sequence of Zackenberg-area moraine ridges that was dated at 13.7-11.2ka, suggesting a rapid ice thinning of >300 m from the high moraine locations.
Our CRE results highlight the rapid rate of deglaciation during the Late Glacial, although the strong age clustering around 13-11 ka and the associated uncertainty ranges hinders the distinction of whether the deposition of the lowest ridges, and the concurrent outer moraine  2014).The proposed deglaciation age of the summits may be associated to the warm period MIS 5a (A).The ages of the lateral moraines are related to the deglaciation of the valley, when the glacier was still debris-free.The ages, partially scattered by paraglacial dynamics, range from 16 to 10 ka (B), and show the most intense deglaciation (peak of the temperature curve) just at the beginning of the Holocene.(C) The ages of the debris-covered glacier show the final phase of the glacier front between 14 and 10 ka.Boulders from the lateral moraines affected by paraglacial processes fell onto the debris-covered glacier, and thus samples retain cosmogenic inheritance, although the peak of the temperature curve may point to the age of its final stabilization at ~10 ka.(D) Ages from the valley bottom confirm an accelerated glacial retreat by ~10 ka, which is also confirmed by the ages from the Lindemansdal valley, where glacial dynamics ended by ~10 ka (E).Fig. 10.Map of NE Greenland including the locations discussed in the text, and exposure ages from published glacial 10 Be CRE chronologies (Biette et al., 2020a(Biette et al., , 2020b;;Håkansson et al., 2007aHåkansson et al., , 2007bHåkansson et al., , 2007cHåkansson et al., , 2009Håkansson et al., , 2011;;Kelly et al., 2008;Larsen et al., 2018;Levy et al., 2014Levy et al., , 2016;;Lowell et al., 2013;Skov et al., 2020).The geographical division (North-East vs. Central-East Greenland) follows that of the Geological Survey of Denmark and Greenland (Dawes and Glendal, 2007).
ridge of the valley floor, occurred during the B-A, YD or Early Holocene (Fig. 9).Therefore, we cannot confirm if these moraines correspond to a glacial readvance occurred during the YD, which has rarely been documented in Greenland (Funder et al., 2021).These results coincide with the initial deposition of the Zackenberg Delta at 13-11 ka (Gilbert et al., 2017), confirming that glacial shrinking was parallel to glacioisostatic uplift and abundant proglacial sedimentation in the lower Zackenberg Valley.
The series of glacial advances/stillstands and retreats that occurred in the Zackenberg area between 13.7 and 11.2 ka is broadly synchronous with glacial records from other sites across NE Greenland indicating that climate conditions favoured glacial shrinking with ephemeral phases of minor glacial advance during the Late Glacial (Fig. 10).A rapid glacial retreat of ~65 km was recorded in Dove Bugt (Fig. 10) between ~13 ka and ~10 ka at a rate of ~22 m yr −1 (Skov et al., 2020).Other 10 Be moraine ages suggest periods of moraine formation during late glacial time: on Store Koldewey island, Håkansson et al. (2007b) reported a post-LGM advance of the GrIS between 14.6 and 11.5 ka, whereas Kelly et al. (2008) described two advances in the Scoresby Sund region at 13.0-11.6and 11.7-10.6ka driven by the prevailing cold conditions associated with the YD.On the other hand, Levy et al. (2016) dated glacial advances of the GrIS and of the peripheral Milne Land Ice Cap at 11.4 ka near Scoresby Sund (Fig. 10).Similarly, on neighbouring Clavering Island, local mountain glaciers advanced during the onset of the Holocene at 11.3 and 10.8 ka, respectively (Biette et al., 2020a(Biette et al., , 2020b)).These heterogeneous ages provide evidence of very active glacial dynamics at the T-1 to Early Holocene transition, with the formation of moraine systems over different timescales.
Deglaciation was accelerated by warmer temperatures during the Early Holocene (Clark et al., 2012;Buizert et al., 2014) .CRE dates confirm that the retreat of the NEGIS started at 11.7 ka and accelerated until 9.3 ka (Larsen et al., 2018).In central East Greenland, the Bregne Ice Cap shrunk within its Late Holocene extent by 10.7 ka (Levy et al., 2014).Based on radiocarbon ages, GrIS retreat within the modern coastline occurred between ~11 and 10 ka (Bennike and Björck, 2002;Vasskog et al., 2015), also including NE Greenland (Wagner et al., 2010;Lecavalier et al., 2014).In the Zackenberg area, our data indicates that the valley was mostly deglaciated by ~10.5 ka, although dead ice patches may have survived for millennia under the debris cover, as already observed in mid-latitude mountains (Fernández-Fernández et al., 2017) and subpolar environments (Charton et al., 2020;Fernández-Fernández et al., 2020).Regional glacio-isostatic rebound and rapid relative sea level fall (Christiansen et al., 2002) led to the formation of the current Zackenberg landscape, bringing the rapid incision of the Zackenberg River, high sedimentation rates, rapid delta progradation and permafrost formation (Gilbert et al., 2017).The end of the YD also favoured the readjustment of slopes to the new ice-free setting, with active paraglacial processes that generated a debriscovered glacier and finally resulted in an extensive hummocky terrain.This process is similar to that observed in the Héðinsdalur valley in the Tröllaskagi Peninsula (north Iceland) where a collapsed debriscovered glacier occupied the valley bottom, showing the typical ridgeand-furrow landscape interspersed with thermokarst collapse depressions of various typologies (Fernández-Fernández et al., 2020;Rodríguez-Mena et al., 2021;Palacios et al., 2021).

Conclusions
This study provides insights into past changes of the GrIS, introducing new geomorphological and chronological evidence of the deglaciation dynamics at the end of T-1when most of the currently ice-free areas in the southern sector of NE Greenland became deglaciatedincluding the progressive transformation of one of its outlets from a debris-free to a debris-covered glacier during the final stages of deglaciation.Greenland still retains 12% of the planet's debris-covered glaciers, although scientific knowledge of their evolution is still limited.In fact, this is the first study focused on collapsed debris-covered glaciers in Greenland that has shown the close relationship of their origin with the paraglacial processes that accompany deglaciation.
With a dataset of 32 CRE dates from erratic and moraine boulders, as well as from polished bedrock surfaces, we established the spatiotemporal pattern of deglaciation in the Zackenberg Valley and the impact that paraglacial activity had on the slopes and the valley floor.Glaciers occupied the fjords and surrounding valleys during most part of the last glacial cycle.Ice reached >800 m a.s.l.above the current valley floor during the maximum advance of the last glacial cycle, when only the highest peaks protruded the ice sheet as nunataks.This may have occurred at ca. 80-50 ka, or later, depending on assumptions regarding nuclide inheritance in the highest erratic boulders.In this study, no geomorphic evidence of glacial activity during the LGM was found, although some of the highest, most heavily eroded remnants of moraines existing in the area may correspond to that phase.The complete sequence of moraine ridges distributed across the slopes surrounding the Zackenberg Valley floor revealed ice thinning between ca.13.7 and 11.2 ka.By 11.3-11.2ka, the glacier formed the lowest moraines distributed on the slopes and pushed up the terminal moraines crossing the central part of the Zackenberg Valley floor.Concurrently, glacial thinning exposed the rock slopes and the moraine ridges, favouring rapid paraglacial slope readjustment.High rates of debris supply onto the glacier still existing in the valley resulted in the formation of a debris-covered glacier.However, as the glacier retreated and disconnected from the external moraine, a hummocky terrain formed as a consequence of areal deglaciation and the subsequent melting of dead ice masses.Overlapping CRE ages from the outermost moraine ridges of the valley floor and the internal hummocky terrain (ca.12-11 ka) indicate that this was a very rapid process.Since the Early Holocene (ca.10.5 ka), no glaciers have been present in the Zackenberg Valley, although ice disappearance has brought the formation of the current landscape due to a wide range of primarily periglacial postglacial processes, in addition to glacio-isostatic uplift and marine terrace and delta formation, permafrost aggradation, and nivation processes largely forming the slopes and the lowland.
Our results demonstrate that CRE dating: (i) needs to be complemented with highly detailed geomorphological mapping; (ii) is a valid method for tracking deglaciation sequences, and; (iii) enables reconstruction of environmental transformations in very active geomorphological settings, including areas that have shifted from debris-free to debris-covered glaciers.However, our results also show that (iv) CRE dating in these areas requires the dating of a large number of samples and landforms; and (v) the uncertainties, considering the problems typically associated with paraglacial processes, remain a challenge for establishing accurate chronologies, as only approximate time ranges can be established for the occurrence of certain events.As such, results need to be supported by other local and regional paleoenvironmental and paleoclimate proxies.
The main phases inferred from our CRE dataset agree with the few existing glacial chronologies from other sites across NE Greenland.However, our results also open new uncertainties that should be addressed in future studies, including the impact of LGM glacial advance on coastal regions and the extent of Holocene glacial fluctuations.A better understanding of the natural pattern of glacial oscillations in NE Greenland during warm and cold phases over the last several millennia would help to better frame the magnitude of current glacial shrinking trends associated with the current warming scenario.
CRediT authorship contribution statement JG, MO, DP, DA, JRF, VJ, VR, TL and KA conducted field work in Greenland.JG, JMF, LL, and TL prepared the samples at CEREGE.IS led the data interpretation.AT performed AMS measurements.HC and OH contributed to the interpretation of glacial evolution.JG, MO, DP, JMF

Fig. 1 .
Fig. 1. (A) Location of the study area within Greenland; and (B) detail of the study region.

Fig. 2 .
Fig. 2. Different types of glacial landforms sampled in this study: exposed polished bedrock surfaces, erratic boulders and moraine boulders.

Fig. 3 .
Fig. 3. Geomorphological map (modified from Cable et al., 2018), including the main glacial and periglacial landforms relevant to our study together with the CRE results shown in Tables1 and 2as well the relevant ages from previous works.In the purple areas termed 'moraine', the solid line indicates the crest of individual moraine ridges; the solifluction symbol indicates that this area consists of till affected by postglacial solifluction.

JFig. 5 .
Fig. 5. Dated boulders from different moraine systems surrounding the Zackenberg Valley floor.See Fig. 3 for the location of the samples shown in the photographs.

Fig. 6 .
Fig. 6.Dated polished bedrock surfaces indicative of subaerial exposure following glacial retreat.See Fig. 3 for the location of the samples shown in the photographs.

Fig. 7 .
Fig. 7. Probability distribution functions of 10 Be CRE ages (with their internal uncertainties) for the samples from the mountain plateaus, slopes and valley bottom.Black vertical lines indicate the mean ages with full errors.

Fig. 8 .
Fig. 8. Idealized reconstruction model for the glacial evolution of Zackenberg Valley (with the base map derived from the current Digital Elevation Model) based on different phases: (1) Near maximum ice extent of the last glacial cycle; (2) Deglaciation of the slopes: 2a-Progressive ice thinning post-LGM; 2b-Formation of moraine systems during the YD; (3) Development of the debris-covered glacier on the Zackenberg Valley floor; and (4) Formation of hummocky terrain.

Fig. 9 .
Fig.9.Probability distribution functions of exposure ages for each group of samples vs. temperature evolution since the last glacial cycle in the interior of Greenland based on the δ 18 O record of the NGRIP ice core (GICC05modelext), 5-point running mean.Periods are labelled according toRasmussen et al. (2014).The proposed deglaciation age of the summits may be associated to the warm period MIS 5a (A).The ages of the lateral moraines are related to the deglaciation of the valley, when the glacier was still debris-free.The ages, partially scattered by paraglacial dynamics, range from 16 to 10 ka (B), and show the most intense deglaciation (peak of the temperature curve) just at the beginning of the Holocene.(C) The ages of the debris-covered glacier show the final phase of the glacier front between 14 and 10 ka.Boulders from the lateral moraines affected by paraglacial processes fell onto the debris-covered glacier, and thus samples retain cosmogenic inheritance, although the peak of the temperature curve may point to the age of its final stabilization at ~10 ka.(D) Ages from the valley bottom confirm an accelerated glacial retreat by ~10 ka, which is also confirmed by the ages from the Lindemansdal valley, where glacial dynamics ended by ~10 ka (E).

Table 1
Sample locations, topographic shielding factor and sample thickness.

Table 2
AMS analytical data and calculated exposure ages. 10Be/ 9 Be ratios were inferred from measurements at the ASTER AMS facility.No correction of erosion and snow cover have been made.
The outliers are highlighted in italics: samples rejected by the Chi 2 test or geochronological inconsistent*, and those with low precision AMS measurements**.See methods and results text for more explanation.
a 10