A critical re‐analysis of constraints on the timing and rate of Laurentide Ice Sheet recession in the northeastern United States

We review geochronological data relating to the timing and rate of Laurentide Ice Sheet recession in the northeastern United States and model ice margin movements in a Bayesian framework using compilations of previously published organic 14C (n = 133) and in situ cosmogenic 10Be (n = 95) ages. We compare the resulting method‐specific chronologies with glacial varve records that serve as independent constraints on the pace of ice recession to: (1) construct a synthesis of deglacial chronology throughout the region; and (2) assess the accuracy of each chronometer for constraining the timing of deglaciation. Near the Last Glacial Maximum terminal moraine zone, 10Be and organic 14C ages disagree by thousands of years and limit determination of the initial recession to a date range of 24–20 ka. We infer that 10Be inherited from pre‐glacial exposure adds 2–6 kyr to many exposure ages near the terminal moraines, whereas macrofossil 14C ages are typically 4–8 kyr too young due to a substantial lag between ice recession and sufficient organic material accumulation for dating in some basins. Age discrepancies between these chronometers decrease with distance from the terminal moraine, due to less 10Be inherited from prior exposure and a reduced lag between ice recession and organic material deposition. 14C and 10Be ages generally agree at locations more than 200 km distal from the terminal moraines and suggest a mostly continuous history of ice recession throughout the region from 18 to 13 ka with a variable pace best documented by varves.


Introduction
Reconstructing ice sheet responses to climate change following the Last Glacial Maximum (LGM) provides valuable insight into ice sheet sensitivity (Stokes et al., 2015) and regional landscape evolution (e.g.Corbett et al., 2019;Koester et al., 2020;Halsted et al., 2022).Potential contributions from shrinking ice sheets are the largest source of uncertainty in projections of future sea-level rise (Bamber et al., 2019); hence, understanding the causes and pacing of ice margin recession in the past is key for reducing uncertainty in future ice recession projections (Fox-Kemper et al., 2021).The development and refinement of geochronologic techniques for dating deglaciation has advanced our understanding of post-LGM ice recession, providing deglacial chronologies that tie ice recession to changing climate conditions derived from dated palaeoclimate records (e.g.Dyke, 2004;Stroeven et al., 2016).Using multiple deglacial chronometers in the same location allows a comparison of independently determined ages and, potentially, the detection of ages biased by violated methodological assumptions (Small et al., 2017).
The northeastern United States (US) has a long history of glacial geologic study and provides an opportunity to reconstruct detailed ice recession patterns using multiple geochronologic methods over hundreds of kilometres.These studies extend back almost 200 years (Hitchcock, 1841;Hitchcock et al., 1861) and attempts to date recession using organic 14 C began in the 1950s (Flint and Deevey, 1951;Flint, 1956).Over the past 70 years, the development of new dating techniques and collection of geochronologic data have produced over 200 organic 14 C and cosmogenic 10 Be ages (e.g.Dalton et al., 2020;Halsted et al., 2022, and references therein) and two calibrated glacial varve chronologies (Ridge et al., 2012;Stanford et al., 2020) that together constrain recession of the Laurentide Ice Sheet (LIS) margin in the region following the LGM (Fig. 1).The large amount and diversity of previously published data provides a unique opportunity for deglacial chronology synthesis and highlights discrepancies between methods in some regions that require reconciliation.
In the sections that follow, we review the chronometers used ('Background: dating ice recession') and the history of glacial geologic studies in the northeastern US ('Background: study area and previous work'), provide a compilation ('Methods') and analysis ('Results') of existing constraints on deglacial chronology, present a new interpretation of the most likely chronology of LIS recession based on all available data ('Discussion') and discuss future strategies for further refinement in this region as well as the wider implications of this work for deglacial chronologies elsewhere ('Implications for dating LIS recession').Most previous investigations of LIS deglaciation in the northeastern US used a single geochronologic technique and focused on local ice recession from landscape features (e.g.Balco et al., 2002;Balco and Schaefer, 2006;Peteet et al., 2012;Bromley et al., 2015;Bromley et al., 2019), or were a subset of a whole-LIS reconstruction effort (e.g.Dyke, 2004;Dalton et al., 2020).In contrast, we use multiple deglacial geochronometers in tandem to reconstruct the timing of southeastern LIS margin recession over 600 km, from its maximum extent at the LGM terminal moraine zone to the present-day US-Canada border.

Approaches to ice recession dating
In the northeastern US, commonly used approaches to constrain the timing and rate of ice recession are: (1) organic 14 C dating, (2) cosmogenic nuclide exposure dating, and (3) glacial varve chronologies. 14C dating and varve chronologies have been used for many decades (Antevs, 1922(Antevs, , 1928;;Flint and Deevey, 1951), while the use of exposure dating is more recent (Balco et al., 2002;Balco, 2011).Each approach has strengths and limitations for reconstructing ice recession patterns; for example, 14 C ages tend to be more precise than exposure ages, but we can measure exposure ages where organic material is absent.However, each approach has inherent methodological assumptions that must be met for ages to be interpreted as ice recession dates.Ages may not correspond to the timing of ice recession if we violate these assumptions.
Organic 14 C dating 14 C dating of organic matter deposited just above glaciogenic sediment in lake and bog sediment cores provides minimumlimiting ages for deglaciation (e.g.Dyke, 2004;Lowell et al., 2009;Fig. 2A).If older, recycled organic material from before glaciation is not part of the sample, 14 C ages reflect the revegetation of land surrounding a lake or bog after glacial ice has receded (Andree et al., 1986).However, sediment remobilisation in lacustrine environments can introduce older or younger organic matter to basal sediment, complicating interpretations of deglacial vegetation changes (Birks et al., 2014).Additionally, aquatic plants can acquire carbon that is not in equilibrium with the atmosphere, especially in the case of lakes in carbonate bedrock areas, known as the hard water effect (Andree et al., 1986). 14C ages are regarded as minimum limits on deglaciation; if used to approximate the timing of deglaciation, the lag between ice recession and organic matter accumulation is assumed to be negligible (e.g.Peteet et al., 2012).The validity of this assumption has been questioned in northern Maine, about 200 km northeast of our study area, due to large differences between minimum-limiting 14 C ages in upland and lowland ponds (Davis and Davis, 1980;Anderson et al., 1986).
Two types of organic carbon-bearing materials have been used to constrain the timing of deglaciation with 14 C dating: bulk sediment and macrofossils.Bulk-sediment 14 C dating, an approach used for over 70 years in the northeastern US (Flint and Deevey, 1951), provides ages from core-bottom sediment containing disseminated organic matter.Prior to the adoption of accelerator mass spectrometry (AMS) in the 1980s ('Laurentide Ice Sheet recessionhistory of previous constraints', below), 14 C dating in sediment cores relied almost exclusively on bulksediment samples because decay-counting required grams of carbon to make a precise measurement (Bronk Ramsey, 2008).However, the use of bulk samples (e.g.Cotter et al., 1984) has been criticised due to the uncertainty about the sources of carbon in sediment (Grimm et al., 2009).Carbon in a bulksediment sample may originate from several sources other than the organic matter directly related to the time of deglaciation (Strunk et al., 2020), including bicarbonate dissolved in lake water and taken up by aquatic organisms, and wind-transported, inorganic terrestrial carbon (Grimm et al., 2009).The 14 C age of bulk sediment therefore represents the mass-weighted mean age of all carbon sources in the sample and may not closely represent the timing of deglaciation.
Unlike the dating of bulk organic material, macrofossil 14 C dating applied to accumulations of identifiable plant or animal fragments led to less ambiguity about the source of carbon being dated.Developed in the 1970s and applied beginning in the 1980s (Andree et al., 1986) as AMS techniques evolved, this method allowed sample sizes to decrease from grams to milligrams.Plant macrofossils include individual leaves, seeds, needles and twigs from deciduous, boreal and tundra vegetation species (Birks, 2003).Bones from Pleistocene megafauna excavated from post-glacial sediment deposits also have been used to provide minimum-date constraints on deglaciation and palaeovegetation changes (Feranec and Kozlowski, 2012;Dalton et al., 2020).Species such as caribou and mammoth inhabited tundra biomes and likely migrated into deglaciated regions shortly after tundra vegetation colonisation (Feranec and Kozlowski, 2016).Macrofossil 14 C dating offers more precise ages than bulk-sediment samples while using samples that are identifiable.However, like bulksediment samples, macrofossil samples can be remobilised within basins, potentially introducing younger or older samples into sediment cores (Birks et al., 2014).
In situ cosmogenic nuclide ( 10 Be) exposure dating Cosmogenic nuclide exposure dating (here using in situ 10 Be) is typically applied in glacial geological studies by measuring isotope concentrations in the top few centimetres of scoured bedrock or glacially transported boulders (Fig. 2B), yielding exposure ages for these features that are interpreted as the timing of ice recession (Balco, 2011).This method assumes that all 10 Be in sampled rock surfaces was produced during the most recent exposure period, and thus the calculated exposure age corresponds to the timing of ice recession.
Several methodological assumptions are made when interpreting 10 Be concentrations as exposure ages.Exposure ages will be too old if 10 Be that formed in rock surfaces during periods of exposure prior to the LGM was not completely removed by glacial erosion (e.g.Briner et al., 2005;Briner et al., 2016).Such atoms are called 'inherited' nuclides.Whereas metres of erosion can occur in settings where ice is warm-based and covers the landscape for many thousands of years, short ice occupation times and/or non-erosive coldbased ice can allow nuclides from prior periods of exposure to remain in rock surfaces when ice recedes (e.g.Bierman et al., 2015;Corbett et al., 2019;Koester et al., 2020).While inherited 10 Be is present on many mountain tops in the region, suggesting that ice cover was cold-based in these locations (Bierman et al., 2015;Corbett et al., 2019;Koester et al., 2020;Halsted et al., 2022), evidence for cold-based ice cover in LGM marginal areas has not been conclusively demonstrated.
In addition to the spallation production described above, 10 Be is also produced in situ by muons that penetrate tens of metres into bedrock (Nishiizumi et al., 1989;Lal, 1991).During long periods of exposure (tens of thousands of years), 10 Be accumulates in concentrations sufficient to alter apparent exposure ages by thousands of years at depths of many metres due to muon production (Briner et al., 2016;Balco, 2017).Surfaces that experienced prolonged interglacial exposure therefore require deeper glacial erosion to remove nuclides created by muon interactions.
Incorrect characterisation of post-glacial subaerial erosion, cover by sediment or snow, or boulder movement that may expose a previously shielded boulder face (Heyman et al., 2011) will cause exposure ages to be biased young relative to the timing of deglaciation (Ivy-Ochs and Kober, 2008).Such young ages have been observed in several boulder samples collected around the northeastern US (Balco et al., 2002;Drebber et al., 2022).

Varve chronologies
During ice recession, sediment with rhythmically alternating fine-and coarse-grained layers, interpreted as varves (annual layers), is deposited in glacial lakes (Antevs, 1922;1928;Fig. 2C).Varves vary in thickness due to annual changes in glacial meltwater input to a lake and distance from the former ice margin.Such thickness variations allow a correlation of varve sequences from different sampling locations and different lake basins (at least 100 km apart in the northeastern US) to create an integrated varve chronology.'Ice-contact' or 'basal' varves lying directly on bedrock or basal glacial deposits are assumed to have been deposited shortly following ice recession from the location (Ridge et al., 2012).Counting up-core from the basal varve provides a chronology of meltwater input to the lake as the ice margin receded, assuming that varves are annually preserved.By correlating varve sequences from multiple locations through patternmatching of diagrams of varve-thickness changes (Ridge et al., 2012), connected varve chronologies sometimes spanning thousands of years have been constructed (e.g.Cato, 1985;Breckenridge, 2007).However, varve chronologies are 'floating' records of ice recession unless they are anchored to the calendar timescale by calibrating the chronology (typically with 14 C ages from organic material, see Ridge et al., 2012).
In the northeastern US, varves and varve sequences (Antevs, 1922(Antevs, , 1928;;Ridge and Benner, 2012;Ridge et al., 2012) have key features that indicate that the couplets are annual, and that the annual layer count is valid.Glacial lakes in the Hudson, Connecticut, Ashuelot and Merrimack valleys where the North American Varve Chronology (see 'Laurentide Ice Sheet recessionhistory of previous constraints', below) was constructed (Ridge et al., 2012), are long, north-south trending water bodies into which a high volume of glacial runoff and sediment flowed from large areas of receding ice.Sediment was deposited over narrow lake floors with strong melt season bottom currents.This process created relatively thick varves as compared with varves in more expansive lakes.The varves have well-defined, contrasting coarse melt season (fine sand to clayey silt) and fine non-melt season (nearly pure clay) layers that are easily distinguished and measured.Individual couplets are not single events such as storm deposits or turbidity (surge) deposits, which do occur as thin beds within the melt season layer of the annual couplets.
The construction of the glacial part of the varve chronology relevant to this paper is based on multiple overlapping, easily matched varve measurements from single lake basins and between hydrologically separate lake basins in the valleys mentioned above.This matching would not be possible if the mechanisms of couplet deposition were non-annual and determined by local sediment input processes.Annual thickness variation in the glacial varves is the result of annual variations in glacial meltwater sediment input that is a primary signal reflecting a regional, annual variation in glacial melting rate (Ridge et al., 2012).
Calibration of the varve chronology in the central and northern part of our study area is based on macrofossil 14 C ages from specific varves (Ridge and Larsen, 1990;Ridge et al., 2012), while the calibration of the varve chronology near the terminal moraine zone is based on bulk-sediment samples and inorganic concretions (Stanford et al., 2020).The macrofossils used for dating were in excellent condition (mostly Dryas leaves) and showed little evidence of decay or erosion during transport.Calibrated ages for the 14 C samples are consistent with the varve year numbers that represent the calibrated timescale with a well-defined offset (Ridge et al., 2012).The calibration of the varve chronology near the terminal moraine is less certain because both bulksediment samples and inorganic concretions are susceptible to contamination by 'old carbon', thus increasing calculated ages (see 'Organic 14C dating', above; Strunk et al., 2020;Grimm et al., 2009).

Glacial geomorphology of the northeastern United States
Our study area encompasses regions of the northeastern US formerly covered by the LIS, specifically northern New Jersey  (Dalton et al., 2020).Ice recession chronometer samples are the same as in Fig. 1.Main: Summed probability density functions for each dating method at selected ice margin positions (identified in inset map).Yellow triangles show the date of ice recession as indicated by varve chronologies at each ice margin position.[Color figure can be viewed at wileyonlinelibrary.com] and Pennsylvania,southern and eastern New York,and New England (Figs. 1,3).Repeated Pleistocene glaciations stripped much of the pre-Quaternary surficial material; the current surficial geology is composed mostly of thin glacial till over bedrock in the uplands and thicker stratified fluvial or lacustrine deposits overlying till and bedrock in the lowlands with infrequent preservation of pre-LGM saprolite and glacial sediment (Stone and Borns, 1986).
The southernmost locations in the study area, south of ~42°N (Figs. 1,3), are relatively flat to hilly, low-lying coastal regions with periglacial, glaciolacustrine and glaciofluvial landforms (Stone and Borns, 1986;Stone et al., 1998).The maximum extent of the southern LIS is marked by a terminal moraine zone that is prominent across northern New Jersey, Long Island (New York) and the islands of Nantucket and Martha's Vineyard (Massachusetts), where the moraine ridges reach 50-100 m in height and 2-10 km in width (Balco et al., 2002;Corbett et al., 2017;Stanford et al., 2020).A series of smaller recessional moraines indicating brief ice margin readvances or stillstands are located 20-40 km up-ice from the terminal moraine in southeastern Massachusetts (Buzzard's Bay moraine; Balco et al., 2002) and Connecticut (Old Saybrook/Ledyard moraine system; Balco and Schaefer, 2006).Boulders lying on the terminal and recessional moraines have been sampled for cosmogenic nuclide exposure dating because of perceived landform stability (Balco et al., 2002;Balco and Schaefer, 2006).
In the central and northern parts of our study area, encompassing areas north of ~42°N (Figs. 1,3; from Massachusetts to northern Vermont and New Hampshire), the landscape is dominated by north-south-oriented mountain ranges with minimal surficial sediment cover separated by wide valleys with deep sediment fills.Previous analyses of multiple cosmogenic nuclides suggest that glacial erosion differed between the uplands and lowlands in this region, with uplands having experienced minimal erosion while lowlands experienced deep excavation (Bierman et al., 2015;Corbett et al., 2019;Koester et al., 2020).During ice recession, large glacial lakes occupied the lowland valleys, and extensive glaciolacustrine deposits, including varves, are found there today (Stone and Borns, 1986;Ridge et al., 2012).Fewer prominent recessional moraines are found in the central and northern parts of our study area, one exception being the White Mountain moraine system in northern New Hampshire (Fig. 3; Thompson et al., 2017).

Laurentide Ice Sheet advance
Dated glacial deposits suggest that prior to the LGM, an ice sheet last overrode the northeastern US completely during the preceding glacial maximum of Marine Isotope Stage 6 (MIS 6), about 130-140 ka.Pre-LGM till located outboard of the LGM terminal moraine in western New Jersey and beneath LGM till elsewhere in New Jersey, southern New England and Long Island are attributed to MIS 6 on the basis of bracketing nonglacial deposits dated to MIS 5 (~70-130 ka) and MIS 9 (~300-340 ka; Stanford et al., 2016;Stanford et al., 2020).Recent modelling efforts tentatively suggest that the LIS margin may have advanced into the northern part of our study area (Vermont, New Hampshire) during the MIS 4 stadial period (~55-70 ka) but receded north and remained out of our field area during most of MIS 3 (~55-35 ka; Dalton et al., 2022).Dating from southeastern Canada support ice-free conditions in the northeastern US during MIS 3 (Stokes et al., 2012;Pico et al., 2018;Dalton et al., 2019), but lacustrine dated deformation tills imply that the LIS margin was in westcentral New York (to the west of our study area in Fig. 1) by ~40 ka (Young and Burr, 2006).
Constraints on LIS growth into the northeastern US are scarce but suggest uninterrupted advance to the terminal moraine from 35 to 27 ka (Table 1).Luminescence dating of pre-glacial cave sediments (Munroe et al., 2016) and 14 C dating of plant fragments in pre-glacial alluvium (Parent and Dube-Loubert, 2017) indicate that the southeastern LIS margin was still north of the US-Canada border before 35 ka.Pre-LGM conifer macrofossils in a lake sediment core (Anderson et al., 1986) and two palaeosol ages (Dorion, 1997) suggest that the LIS margin entered northern Maine by 33 ka and perhaps as late as 28 ka (Bierman et al., 2015, their supplement).Luminescence dating of glaciolacustrine sediment under LIS till indicates that the ice margin reached the Catskill Mountains of southern New York by 28 ka (Rayburn et al., 2015).Several 14 C ages from sediment underlying till near the terminal moraines indicate that the LIS advanced rapidly to its maximum extent, which it reached by 27-25 ka (Sirkin and Stuckenrath, 1980;Stone and Borns, 1986;Stanford et al., 2016;Stanford et al., 2020).

Laurentide Ice Sheet recessionhistory of previous constraints
Early attempts to constrain the timing of LIS recession from the LGM terminal moraine relied exclusively on bulk-sediment 14 C ages from the lowermost organic material accumulated in lake sediment cores (e.g.Flint and Deevey, 1951;Flint, 1956; see 'Organic 14 C dating', above).The number of core-bottom 14 C ages in the region increased throughout the 1970s and 1980s, culminating in the LIS deglaciation map of Dyke and Prest (1987).This map continued the tradition of drawing LIS margins during successive steps of deglaciation, informed by the constraining ages, mostly uncalibrated 14 C ages from lake core bulk sediment or marine fossils, for ice recession available at the time (e.g.Prest, 1969;Mayewski et al., 1981).Dyke and Prest (1987) assigned a date of 21.7 ka for maximum ice extent along most margins of the LIS (Table 2; all 14 C ages and derived dates reported here are calibrated to years before 1950, see 'Data collection and calibration', below).Ice margin recession in the northeastern US was rapid in the reconstruction of Dyke and Prest (1987) who concluded that the ice margin had receded to central Vermont and New Hampshire (~400 km) by 16.8 ka and was at the US-Canada border (~600 km) by 15.5 ka (Table 2).
The advent of AMS for radiocarbon dating in the mid-1980s sparked a new wave of deglacial chronologies in the northeastern US based on core-bottom macrofossil samples (e.g.McWeeney, 1999;Peteet et al., 2012).A consequence of the rapidly expanding suite of AMS macrofossil 14 C ages was the labelling of many existing bulk-sediment 14 C ages as 'anomalously old', with the assertion that they were contaminated by old carbon (Dyke, 2004).New LIS deglacial chronologies from around the ice sheet favoured macrofossil 14 C ages and featured ice margin recession that was delayed in comparison with previous reconstructions.In these reconstructions, ice recession in the northeastern US started around 20.5 ka (Dyke, 2004; Table 2).
By 2004, varve chronologies (see 'Varve chronologies', above) provided constraints on the timing and rate of deglaciation in the northeastern US (Ridge, 2004).The North American Varve Chronology (NAVC; Ridge, 2004;Ridge et al., 2012) begins in southern Connecticut, about 70 km north of the terminal moraine.After the NAVC was anchored by 14 C ages, it appears to record ice recession from 18.6 to 13.4 ka, with the date of initial deglaciation from the terminal moraine estimated as 24-28 ka (Ridge, 2004; Table 2).The NAVC also makes clear that the rate of ice recession throughout the northeastern US was more variable than depicted in earlier, 14 C-based, chronologies.The speed of ice recession varied between 30 and 100 m a −1 from 18.6 to 14.6 ka, while ice recession after 14.6 ka accelerated to ~300 m a −1 until ice receded north of the modern US-Canada border (Ridge et al., 2012).AMS 14 C ages on macrofossils in varves suggested that the ice margin was north of the US-Canada border by ~13.5 ka (Ridge, 2004; Table 2), almost 2 kyr later than suggested in Dyke and Prest (1987).
The introduction of cosmogenic nuclide exposure dating for constraining ice recession (see 'In situ cosmogenic nuclide ( 10 Be) exposure dating', above) in the northeastern US provided ages for inorganic, glaciogenic deposits and landforms (moraines, erratic boulders and sculpted bedrock) that could not be dated directly with 14 C.As more data accumulated, it became clear that 14 C and cosmogenic nuclide exposure chronologies did not always agree.For example, the majority of exposure ages from boulders on the Martha's Vineyard terminal moraine are older than 25 ka (Balco et al., 2002; Table 2; all exposure ages are recalculated here using originally reported nuclide concentrations, see 'Methods'), much older than the ~20.5 ka deglaciation date suggested by Dyke (2004).Exposure ages from recessional moraines in southern Connecticut, between 20 and 40 km upice from the terminal moraine, are also older than suggested by Dyke (2004;19.3 ka), clustering around 21 ka (Balco and Schaefer, 2006; Table 2).
As more AMS macrofossil 14 C and cosmogenic nuclide exposure age samples were analysed, debate over the timing of initial LIS deglaciation drove further investigations.In 2012, a compilation of new and previously published core-bottom macrofossil 14 C ages (n = 13) located near the terminal moraines supported assertions that deglaciation was delayed until after 17 ka (Peteet et al., 2012; Table 2).In contrast, exposure ages sampled from boulders and bedrock 4 km north of the terminal moraine in New Jersey (n = 16; Corbett et al., 2017) centre around 25.0 ± 2.1 ka (Table 2), agreeing with nearby bulk-sediment 14 C ages that cluster between 22 and 27 ka (Sirkin and Stuckenrath, 1980;Cotter et al., 1986;Stone and Borns, 1986;Stone et al., 1989).
A new LGM and deglacial varve chronology from the Hudson River valley in southern New York and northern New Jersey (Stanford et al., 2020;Fig. 3) supports maximum ice extent between 25 and 24 ka followed by rapid recession (80 m a −1 ) from 24 and 23.5 ka and then slowing to ~12 m a −1 from 23.5 to 22.5 ka.The southern Hudson River valley chronology is not linked directly to the NAVC but is thought to be 3.5-4 kyr older (Stanford et al., 2020).However, this chronology is calibrated with 14 C ages from a combination of concretions and bulk-sediment samples, which some believe have a higher probability of old carbon contamination than macrofossils (Grimm et al., 2009;Peteet et al., 2012).
In summary, 14 C-based deglacial chronologies in the northeastern US are most consistent with initial deglaciation around 18-19 ka (and as late as 17 ka), while varve and cosmogenic exposure nuclide-based chronologies suggest ice recession began around 24-25 ka (and potentially as early as 28 ka).This large chronologic disagreement complicates interpretations of ice recession pacing and presents a challenge for determining the forcings that may have initiated abandonment of the terminal moraine zone.

Deglacial palaeoclimate
Changing palaeoclimate conditions in the Northern Hemisphere, thought to be the primary driver of ice sheet mass loss during the last deglacial period (e.g.Gregoire et al., 2015;Ullman et al., 2015;Stokes, 2017;Halsted et al., 2022), were primarily influenced by variations in incoming solar radiation (insolation), ocean circulation and atmospheric greenhouse gas concentrations (Osman et al., 2021).Despite summer insolation at high northern latitudes increasing after its LGM low around 24 ka (Fig. 5; Laskar et al., 2011), sustained cold glacial conditions persisted in the Northern Hemisphere until around 17 ka.This lag was due to the strong radiative cooling effect of sea ice and continental ice sheets and low atmospheric greenhouse gas concentrations (Tierney et al., 2020), with terrestrial annual temperatures perhaps ≥20°C colder than preindustrial conditions across ice sheet-covered areas (see regional reconstruction in Osman et al., 2021).Annual average temperatures in the Northern Hemisphere began to rise after 17 ka (Osman et al., 2021) as atmospheric greenhouse gas concentrations increased (Shakun et al., 2012) and terrestrial albedo decreased due to shrinking ice sheets (Braconnot and Kageyama, 2015;Baggenstos et al., 2019;Tierney et al., 2020).Model simulations suggest that while average annual temperatures in the North Atlantic remained low until 17 ka, summer temperatures may have begun rising as early as 18 ka (Buizert et al., 2014), perhaps reflecting the ongoing increase in summer insolation at high northern latitudes.
The steady deglacial warming that began around 17 ka in the Northern Hemisphere was followed by an abrupt and substantial warming around 14.7 ka (Bølling warming; Osman et al., 2021) that was especially pronounced in the North Atlantic, where average annual temperatures increased by as much as 14°C in a few decades (Lea et al., 2003;Buizert et al., 2014;Kindler et al., 2014;Ivanovic et al., 2016).A resurgence in Atlantic Meridional Overturning Circulation and thus ocean heat transport into the North Atlantic is typically invoked to explain this abrupt Bølling warming (Liu et al., 2009;Obase and Abe-Ouchi, 2019).The Bølling-Allerød warm period persisted from 14.7 to around 13 ka, interrupted by a brief cold interval around 14 ka (Buizert et al., 2014) and followed by a return to glacial-like wintertime conditions in the North Atlantic region during the Younger Dryas cold period at ~12.9 ka (Osman et al., 2021).The LIS margin was very likely north of our study area by the Younger Dryas (e.g.Thompson et al., 2017;Dalton et al., 2020), and we thus do not consider this period in depth, but we note that the return to glacial-like conditions during this period is associated with the formation of prominent moraines as the ice margin readvanced in southeastern Canada (e.g.Couette et al., 2023).
For organic 14 C ages, we classified core-bottom samples as 'bulk sediment' or 'macrofossil' based on descriptions provided in the original publication or data source.We calibrated all 14 C ages to calendar years using MatCal (Lougheed and Obrochta, 2016), the IntCal20 calibration curve (Reimer et al., 2020) and the originally reported 14 C age and uncertainty.For ice recession modelling and data presentation in figures, we used median calibrated ages and half of the 68.2% probability density interval to approximate a 1σ Gaussian age uncertainty.We provide the full calibrated age range for 14 C samples in Tables S2-S3.
For cosmogenic nuclide exposure ages, we recalculated all 10 Be ages using the exposure age calculator formerly known as CRONUS-Earth, code version 3 (CRONUS-V3; Balco et al., 2008; http://hess.ess.washington.edu/),the global production rate calibration dataset (Balco et al., 2008) and the 'LSDn' scaling scheme (Lifton et al., 2014).We used the sample parameters (location, thickness, shielding) reported in the original publications, 10 Be concentration measurements, standards used for normalisation at the time of analysis, and assumed no erosion or snow cover when calculating exposure ages (Table S1).We report all 10 Be exposure ages in the 'Results' and 'Discussion' sections with external uncertainties (i.e.including both analytical and production rate uncertainties) to better facilitate the comparison of ages across different chronological methods.
We chose not to use a regional cosmogenic nuclide production calibration dataset for northeastern North America (Balco et al., 2009), as it uses the NAVC as an independent age constraint for glacially deposited boulders; thus, ages determined with the northeastern North American production rate are not independent from those determined using glacial varve chronologies.We performed a sensitivity test to determine the impact of using different production rate calibration datasets (Supplementary information) and found that exposure ages are <1% higher with the global dataset compared with the northeastern North America dataset when using LSDn scaling.These age differences are less than analytical uncertainty of exposure ages (typically 3-6%).
During data collection, we evaluated every sample for its suitability as a deglacial age constraint and rejected 33 ages from our reconstruction of ice recession (Tables S1-S3).We followed the methodology of Dyke (2004) to re-evaluate samples formerly considered robust in the light of new data.Bogs and lakes were sometimes revisited by researchers as coring equipment and isotope measuring technology improved and allowed deeper cores and more precise ages.In these cases, we discarded 14 C samples from older publications if a more recent publication contains a deeper and older 14 C age from the same or nearby basins (e.g.Cotter et al., 1984;Cotter et al., 1986) (Tables S1-3).We removed 26 bulksediment 14 C samples, two plant macrofossil samples and two megafauna bone samples following this approach.We also removed five 10 Be exposure ages identified as outliers in the original publications (Balco et al., 2002;Drebber et al., 2022).

Deglacial chronology construction
For our deglacial chronology, we adopted the ice margin positions delineated by Dalton et al. (2020;Fig. 3), which are themselves updated from positions suggested by Dyke (2004).These margin positions reflect the presumed importance of topography for controlling ice recession in this high-relief setting, and deglacial ages along these hypothesised margins are expected to be similar.We used ages in our compilation that are south of each ice margin position, but north of the next older ice margin position, to inform dates of ice recession.For each ice margin, individual ages are represented by probability density functions that are summed together, treating ages from bulk-sediment, macrofossil and cosmogenic samples separately, to show the overall age probability distribution for each geochronometer at ice recession positions (Figs. 3, 4).
We employed a Bayesian approach (Blaauw and Christen, 2011) to model ice margin recession through the northeastern US according to each deglacial chronometer, as in Lowell et al. (2021).Like Lowell et al. (2021), we modelled ice recession separately between moraines, which represent ice sheet stillstands or readvances.The moraines in our model are the Buzzard's Bay moraine system in southeastern Massachusetts (ice margin position 2 in Fig. 3), the Old Saybrook/Ledyard moraine system in southern Connecticut (position 3 in Fig. 3), and the White Mountain moraine system in northern New Hampshire (position 10 in Fig. 3).Thus, using the ice margin positions in our Fig. 3, our model recreates ice recession separately between positions 1 and 2, 2 and 3, 3 and 10, and (for radiocarbon ages), between positions 10 and 11.

Be inheritance modelling
To assess the potential presence and magnitude of 10 Be inheritance throughout our study area, we modelled potential in situ 10 Be inheritance concentrations in exposure age samples.We calculated 10 Be concentrations that accumulated in bedrock depth profiles for a range of pre-LGM exposure durations, and determined how much 10 Be would be left in sampled surfaces given different of glacial erosion during LGM ice cover (following the methodology in Briner et al., 2016).
We determined pre-LGM exposure histories for surfaces in the northeastern US using constraints for glacial occupation over the last glacial/interglacial cycle (see 'Laurentide Ice Sheet advance', above), which provide reasonable estimates for a pre-LGM exposure duration of ~80 to ~110 kyr.We calculated 10 Be production rates from both neutrons and muons (Balco, 2017) down to 10 m depth in bedrock (ρ = 2.65 g cm −3 ), assuming a constant, total surface sea level 10 Be production rate for the northeastern US (~42°N) of 4.3 atoms g −1 a −1 (Balco et al., 2008).We calculated the accumulated concentration of 10 Be in these depth profiles (Fig. 6A) and converted 10 Be concentrations into equivalent surface exposure times by dividing by the surface production rate of 4.3 atoms g −1 a −1 , a value appropriate for our predominately lowelevation sites (Fig. 6B).

Deglaciation dates
The agreement between 10 Be and 14 C ages varies over space.There is a wide range of ages constraining the initial recession of the southeastern LIS margin from its terminal position (positions 1-3 in Fig. 3; Tables S1-3), with substantial disagreement between 10 Be and 14 C age estimates.This disagreement narrows with increased distance from the terminal moraine zone.The difference in 10 Be and 14 C ages at the terminal moraine zone is 8-10 kyr, which decreases to 3-6 kyr at recession positions 2 and 3, after the ice margin had receded 20-40 km (Figs. 3, 5).The 14 C and 10 Be ages are similar, considering data scatter and age uncertainties, at locations more than 200 km north of the terminal moraine zone (starting at position 6 in Fig. 3).Ages from all chronometers remain similar into northern New Hampshire and Vermont (positions 9 and 10 in Fig. 3), where the central tendency of ages is between 13 and 14 ka.
Ages from all chronometers exhibit more variance near the LIS terminal moraine zone than farther north (Table 3).Bulksediment 14 C ages in particular exhibit substantial variance in the southern part of our study area.Considered together, bulk 14 C ages within ~100 km of the LIS terminal moraine zone (margin positions 1-4 in Fig. 3) range from 12.3 to 27.2 ka with a relative standard deviation (RSD) of 20.3% (n = 38).For comparison, macrofossil 14 C ages within ~100 km of the terminal moraine zone (n = 20) range from 12.8 to 18.5 ka with an RSD of 9.5%, and 10 Be exposure ages (n = 49) range from 18.9 to 33.6 ka with an RSD of 13.8%.In the northern parts of our study area (margin positions 7-11), ages from all geochronometers exhibit less variance than near the terminal moraine, with RSD values between 2.4 and 8.5% (Table 3).

Bayesian model results of ice recession
Our modelling of ice margin recession implies markedly different deglacial chronologies based on different dating methods (Fig. 5).The timing and rate of ice recession based solely on 10 Be exposure ages is generally similar to calibrated varve chronologies across much of the central and northern parts of our study area, but recession from the terminal moraine starts earlier in the 10 Be chronology and at a slower rate (Fig. 5).The 10 Be-based model suggests moderately paced ice recession (~50 m a −1 ) throughout much of Connecticut and southern Massachusetts followed by an acceleration (up to ~70 m a −1 ) around 18 ka.Ice recession in this model continues at ~70 m a −1 until around 13.4 ka, at which point the ice margin was in northern Vermont and New Hampshire (position 10 in Fig. 3) and no exposure ages constrain its recession further north.
Ice recession chronologies based on macrofossil and bulksediment 14 C ages differ slightly, but both suggest a far later start of ice recession than either exposure ages or varve chronologies, and at rates that mostly exceeded 100 m a −1 (Fig. 5).Recession from the LGM terminal moraine zone based solely on macrofossil 14 C ages starts at ~15 ka and at a rate of ~120 m a −1 , then accelerating to ~355 m a −1 around 14.5 ka, when the ice margin had receded ~100 km to central Connecticut (position 4 in Fig. 3).Recession based on bulksediment ages starts at ~16 ka and at a rate of ~185 m a −1 , with an acceleration after ~100 km of recession up to ~265 m a −1 around 15.5 ka.These faster modelled recession rates after ~100 km of ice margin recession continued virtually unchanged throughout northern New England, although our macrofossil-based model suggests that recession slowed to ~245 m a −1 after ~13 ka, when the ice margin had receded >500 km to northern New Hampshire and Vermont (position 10 in Fig. 3).

Be inheritance and glacial erosion
Assuming 80-110 kyr of pre-LGM exposure when the LIS was absent near the terminal moraine zone, our modelling suggests that 10 Be was present in measurable quantities down to more than 10 m depth in the pre-LGM bedrock (Fig. 6).Changes in ** Relative standard deviation is the quotient of the standard deviation and average of ages from each ice margin, multiplied by 100 to be expressed as a percentage.~2.5 m, where muons dominate 10 Be production (Fig. 6B).Age-equivalent inheritance decreases from 20 kyr near the surface to 4 kyr at 2.5 m depth, and to about 1 kyr at 11 m depth (Fig. 6A).Our inheritance estimates are minima because we modelled only the last glacial/interglacial cycle and not 10 Be accumulated deep in bedrock during multiple previous interglacial exposures.

Discussion
Our compilation of ages recording LIS deglaciation shows that uncertainties remain regarding the timing of LIS margin recession across the northeastern US, despite decades of sampling and analysis.More than 200 km north of the LIS terminal position, in Vermont and New Hampshire, organic 14 C, varve and cosmogenic nuclide chronometers all give similar results.In contrast, there are multiple, conflicting estimates for the age of the terminal and nearby (20-40 km) recessional LIS moraines in southeastern Massachusetts and southern New York.There are also discrepancies between chronometers during the initial ice margin recession across Connecticut (Figs. 3, 5).These conflicting estimates for the timing, and therefore pace, of ice margin recession indicate that some assumptions intrinsic to the application of each dating method were violated.In the sections that follow, we synthesise existing age data, compare the results of each chronometer, assess the probable causes of these differences, evaluate constraints on the timing of initial ice recession from the LIS maximum moraines, and suggest a way forward.

Calibrated glacial varve chronologies
Glacial varve chronologies provide the highest resolution record of ice recession through the northeastern US and serve as important benchmarks against which to compare the pace of ice recession implied by 14 C and 10 Be ages, but uncertainties remain about the timing of deglaciation according to varve records near the terminal moraine zone.Although the varve chronologies considered here (Ridge et al., 2012;Stanford et al., 2020) are continuous and thought to have minimal counting errors, they, like all varve chronologies, are floating.To use varve chronologies for dating ice margin recession they must be tied to the calendar timescale, typically by 14 C dating of organic material, preferably plant macrofossils identified to species level, within the varve sequence (see 'Varve chronologies', above).
This key requirement in this geochronological approach causes complications for inferences about the date of initial recession.The southern Hudson River valley varve chronology (Stanford et al., 2020), which encompasses the first ~60 km of ice margin recession from the terminal moraine in New Jersey and southern New York, suggests recession beginning at 24 ka.But the 14 C calibration samples anchoring this chronology to the calendar timescale, inorganic concretions and organic material in bulk sediment, are both susceptible to age biases from incorporation of non-atmospherically derived 'old carbon' (see 'Organic 14 C dating', above; Grimm et al., 2009;Strunk et al., 2020).Thus, although the Hudson River valley varve chronology provides key information about the pace of ice margin recession from the terminal moraine (initially 80 m a −1 for the first 500 years and then slowing to ~12 m a −1 for 1 kyr), the inferred timing of deglaciation could be too old due to the nature of the calibration samples.
The calibration of the NAVC to calendar years is more certain (Ridge et al., 2012;Balco et al., 2021) with repeated measurements over its ~5.6 kyr extent.However, the NAVC does not extend to the terminal moraine zone; the oldest varves correlated to the NAVC are at Newburgh, New York, approximately 100 km north of the terminal moraine on Staten Island.The NAVC is calibrated using 14 C ages from macrofossils embedded in 54 separate varves that span over 4800 yr of the NAVC.The scatter of residuals between the best-fit calibration age model and individual macrofossil ages suggests a 0.2 kyr uncertainty in the age of varves (Ridge et al., 2012), far less than the spread of either 14 C or 10 Be ages.The macrofossil 14 C samples used to calibrate the NAVC were deposited in the central parts of our study area later in the deglacial period (after ~18 ka) and thus are likely to closely reflect the timing of deglaciation.A recent, independent  attempt to calibrate the NAVC to the calendar timescale using patterns of meteoric 10 Be fallout in varves was inconclusive for much of the varve record but suggested that the 14 C-based calibration is valid due to overlapping best-fit calibration models for sampled portions of the varve chronology (Balco et al., 2021).

Organic 14 C deglaciation ages
There is a mismatch between the LIS recession chronology suggested by our compilation of 14 C ages and what we know about deglacial climate.Taken at face value, the oldest macrofossil 14 C ages near the terminal moraine zone (within 25 km of maximum ice extent) suggest that the LIS margin remained near its maximum position until after 18 ka, despite increasing summer insolation (Laskar et al., 2011;Fig. 5).
Bayesian modelling using all radiocarbon ages near the terminal moraine zone suggests even more recent most-likely dates for ice recession: 16.1 ± 0.2 ka (bulk-sediment samples) or 15.2 ± 0.4 ka (macrofossil samples), several thousand years after other palaeoclimate and sea-level records indicate that the climate was warming significantly (Buizert et al., 2014;Lambeck et al., 2014;Lora and Ibarra, 2019;Osman et al., 2021).Such long-lasting ice seems unlikely given the warming indicated by independent climate records.Many bulk-sediment ages are slightly older than nearby macrofossil ages (Fig. 5).This pattern could indicate that bulksediment ages better reflect the timing of ice recession because they capture earlier trace organic matter accumulation in basins before sufficient macrofossils for dating are present on the landscape (such as in Briner et al., 2013).However, older bulk ages may result from incorporation of carbon that is not related to the timing of deglaciation (e.g.Grimm et al., 2009;Strunk et al., 2020).Such 'old' carbon could be sourced from lake algae that respired non-atmospheric carbon dioxide in lake water (Meyers and Teranes, 2002) or sourced from older particulate organic material incorporated by the LIS as it overran earlier soils (Strunk et al., 2020).
The apparently young ice recession date (~16-18 ka) that results from compiling organic 14 C ages near the LIS terminal moraine zone likely reflects a delay between deglaciation and significant organic matter deposition in basins; thus, minimum-limiting organic 14 C ages probably post-date the timing of deglaciation by up to several thousand years.Such a delay is plausible in the northeastern US, where year-round stadial conditions persisted until ~18 ka (Buizert et al., 2014;Lora and Ibarra, 2019).These cold, dry conditions, and the scant vegetation and occurrences of permafrost (Stone and Ashley, 1992) that resulted during the early deglacial period, were primarily controlled by changes in the North Atlantic regional climate (Shuman et al., 2004;Buizert et al., 2014;Osman et al., 2021).We surmise that during this relatively cold period, as much as several thousand years elapsed after the ice margin receded and before the landscape became vegetated enough that organic material, particularly macrofossils, were deposited in lake basin sediment in sufficient numbers to be found in sediment cores.
The increasing density of vegetation as the climate warmed and deglaciation continued led to more robust 14 C dating during the latter stages of LIS recession.An increase in North Atlantic summer temperatures around 18 ka (Buizert et al., 2014) is associated with the first changes to predominant LGM vegetation patterns outboard of the LIS margin (Stone and Borns, 1986;Peteet et al., 1994).This summer warming presumably led to widespread revegetation of deglaciated landscapes in the northeastern US and thus the deposition of sufficient organic material in lake sediment to allow reliable, consistent dating.This revegetation, rather than ice margin recession, was the likely source of many macrofossil ages of ~18 ka at and near the LIS terminal moraine.
Overall, 14 C dating of organics appears to be biased toward younger ages due to cold, late-glacial climate and its slowing of post-glacial recolonisation of the landscape by plants.Most organic 14 C analyses used to determine the timing of LIS recession provide only minimum-limiting ages for at least the first few thousand years after initial deglaciation, before the climate warmed and before more extensive vegetation became established.In contrast, after about 18 ka, organic 14 C analyses record the timing of ice recession more faithfully because the climate was sufficiently warm to allow for rapid vegetation of deglaciated land surfaces.Therefore, LIS reconstructions that rely entirely on organic 14 C ages are likely to underestimate the timing of initial LIS retreat from the terminal moraine zone, thus inflating retreat rates calculated from these ages.

Be deglaciation ages
Given other constraints on LIS margin recession, many cosmogenic 10 Be exposure ages on the terminal moraine at Martha's Vineyard appear unreasonably old and therefore unlikely to correspond to moraine abandonment by the ice margin.Even after the exclusion of much older outliers, remaining ages range from 30 to 25 ka (Balco et al., 2002), tying ice recession to a time when there were no known climate forcings that would trigger such early deglaciation and when the global sea level was still falling (Lambeck et al., 2014).These older ages could reflect reworking of early MIS 2 glacial deposits into the terminal moraine as the ice margin fluctuated near its maximum extent, or inheritance of nuclides from previous near-surface, interglacial exposure of sampled clasts.In the latter case, boulders on Martha's Vineyard contained varying, and in many cases significant, concentrations of 10 Be when they were first exposed on the LIS terminal moraine, thus biasing exposure ages too old.
The wide variance of exposure ages on Martha's Vineyard (n = 12, RSD = 42%, Balco et al., 2002) is strong evidence of shallow and spatially variable glacial erosion.Inherited nuclide concentrations decline rapidly within a few metres of exposed bedrock surfaces (Fig. 7); as a result, large (>1 m 3 ), boulders that have been glacially exhumed from shallow depths (less than a few metres) carry substantially different concentrations of 10 Be on all their surfaces (Fig. 7).It is this lack of deep erosion that leads to substantial 10 Be inheritance in sampled boulder surfaces on Martha's Vineyard and thus inflated and inaccurate exposure ages.Even as boulders are reoriented during glacial transport and deposition, 10 Be inheritance in sampled surfaces (and adjacent boulders) can differ by thousands of years.Such nuclide inheritance probably occurs on the Martha's Vineyard moraine, as suggested by the high RSD.At depths of 3-5 m below a bedrock surface exposed during the last interglacial, 10 Be concentrations decline more gradually with depth but are still sufficiently large to represent the equivalent of thousands of years of surface exposure (Briner et al., 2016).
Overall, cosmogenic nuclide ages appear to have an older bias due to limited glacial erosion.This bias is greatest toward the terminal moraine zone, where short ice occupation times minimised subglacial erosion (Melanson et al., 2013).Further from the terminal moraine zone, cosmogenic nuclide exposure ages agree with the NAVC and thus we interpret them to be more accurate (Figs. 3,5).LIS reconstructions that rely entirely on cosmogenic ages are likely to overestimate the date of initial LIS retreat from the terminal moraine zone, and thus underestimate initial retreat rates calculated from these ages.
When did the LIS recede from the terminal moraine?
Several lines of evidence, none direct, suggest that the LIS began receding from the terminal moraine between 24 and 20 ka.Evidence for ice recession starting no earlier than 24 ka comes from ice advance constraints and estimates of ice occupation duration at terminal moraine segments.LIS advance constraints (pre-LGM 14 C and OSL ages; see 'Laurentide Ice Sheet advance', above) indicate that the southeastern LIS margin did not reach its maximum extent before 27-26 ka (Sirkin and Stuckenrath, 1980;Munroe et al., 2016) and perhaps as late as 25 ka (Stanford et al., 2020).Ice remained at or near maximum extent for 1-3 kyr, as estimated by varve records outboard of the terminal moraine in New Jersey (Stanford et al., 2020).The internal morphology of the Martha's Vineyard terminal moraine suggests several kyr of ice occupation (Oldale and O'Hara, 1984;Balco et al., 2002).
Abandonment of the terminal moraine zone around 24 ka implies that either the ice margin was sensitive to subtle climatic forcings or regional palaeoclimatic changes were different from those portrayed in existing large-scale reconstructions (Osman et al., 2021).Few plausible climate forcings existed around 24 ka to explain ice margin recession aside from the initial rise in summer insolation at 40°N that began at this time (Fig. 5; Laskar et al., 2011).Coupled climate-ice sheet modelling supports the notion that small insolation changes are sufficient to trigger margin recession for large ice sheets (Abe-Ouchi et al., 2013), and rising insolation has been invoked as a mechanism for deglaciation elsewhere along the southern LIS margin (Ullman et al., 2015).
Evidence for LIS margin recession having started before 20 ka comes primarily from varve chronologies (Ridge et al., 2012;Stanford et al., 2020).The NAVC indicates that ice recession 100 km north of the terminal moraine in New Jersey occurred at about 18.3 ka (Ridge et al., 2012) and at a pace (~50 m a −1 ) that provides a plausible minimum constraint of ~20 ka for the abandonment of the terminal moraine.Initial deglaciation as indicated by the southern Hudson River valley varve chronology (Stanford et al., 2020), estimated at 24 ka but with some calibration uncertainty ('Calibrated glacial varve chronologies', above), is consistent with tightly clustered 10 Be exposure ages (boulders and bedrock) and both organic and inorganic 14 C ages within ~25 km of the New Jersey terminal moraine that suggest a date for initial deglaciation of 25.2 ± 2.1 ka (Corbett et al., 2017).
Compared with other sites near the terminal moraines, the New Jersey study area of Corbett et al. (2017) is unusual because the 14 C and 10 Be ages agree well.The 10 Be exposure ages (n = 16) are from 10 bedrock outcrops and six boulders.There is no difference between the average age of boulder and bedrock samples, suggesting that inherited nuclide concentrations in all surfaces are minimal, and there are no older or younger outliers in the dataset.In this area, five different 14 C ages bracket most of the 10 Be exposure ages and there is little excess scatter in the exposure ages (RSD < 10%), in sharp contrast to the >40% RSD measured on the Martha's Vineyard section terminal moraine.
The low RSD documented in the New Jersey exposure age samples suggests that the depth of erosion at the sampled sites was significant, stripping metres of rock and/or sediment from glaciated bedrock surfaces.The extensive jointing in the quartzite at this site may have promoted the quarrying of thick blocks of rock.To have similar ages, the sampled boulders must also have been quarried from similarly eroded areas; perhaps ice-streaming characterised this part of New Jersey and more effectively eroded the bed than ice flowing over Martha's Vineyard.Organic 14 C biases may be minimised near this site as well; the oldest bulk-sediment 14 C ages (Cotter et al., 1986) agree with nearby 14 C ages from concretions that formed in ice marginal lakes (Stone et al., 1989).These 14 C ages of inorganic material do not require the presence of vegetation.
Even in New Jersey, where the two chronometers agree well, there remains uncertainty about the veracity of the LIS terminal moraine age estimate.The mean 10 Be exposure age in New Jersey (25.2 ka; Corbett et al., 2017) is older by at least 1 kyr than other inferences suggested for the age of the LIS terminal moraine zone.Measurements in other areas of continental glaciation, such as Norway and Greenland, suggest inheritance of cosmogenic nuclides, primarily those generated by deeply penetrating muons, is sufficient even in deeply eroded terrains to increase calculated exposure ages by hundreds of years to a few thousand (Briner et al., 2016).We recognise that the 14 C ages, all from concretions and bulk sediment, could be too old.Therefore, although agreement in New Jersey between chronometers suggests confidence in both the 10 Be exposure and organic 14 C ages, such agreement does not necessarily guarantee accuracy.

Ice recession rates
Our model suggests that ice recession rates increased over time, beginning slowly as the ice retreated from the terminal position and gradually accelerating, reaching maximum values in northern New England during the Bølling and Allerød warm periods (14.5-14.1,13.9-12.9ka, respectively), which were interrupted by a brief cold period around 14 ka.Because the age of the terminal moraine zone remains uncertain, we begin LIS recession rate calculations about 30 km inland and up-ice.We use 10 Be-based ice recession dates to cover the spatial gap in varve chronologies in the southern part of Connecticut; the good agreement between the 10 Be-based recession model and the NAVC further north lends confidence to these calculations.After 18 ka, ice recession rates agree well between the calibrated varve and 10 Be-based chronologies, if we smooth over the high-resolution details of the varve chronology, and thus we base our recession rate calculations on the combined datasets.We calculate that ice recession began in southern Connecticut (after vacating the Old Saybrook/Ledyard moraine system) between 21 and 20 ka and at a rate of ~50 m a −1 .The calibrated varve chronology in the nearby Hudson River valley of southern New York suggests earlier and slower ice recession (Figs. 3, 5;Stanford et al., 2020); this timing disagreement could reflect an early recession of the ice lobe in the southern Hudson River valley compared with the Connecticut River valley that is not reflected in existing deglacial reconstructions of the region (Dalton et al., 2020).
In northern Connecticut and southern/central Massachusetts, both calibrated varve chronologies and our 10 Be-based age model support ice recession at the same rate (~50 m a −1 ) as in southern Connecticut, implying that the ice margin reached central Massachusetts by ~18-17.4ka.While the most likely dates of ice recession throughout this region from 10 Be ages suggest slightly earlier recession compared with the calibrated varve chronology (~19.2-18ka vs. ~18.6-17.4ka), probably due to the presence of inherited nuclides that bias exposure ages too old, recession dates suggested by the two methods overlap when considering the 95% confidence interval of the 10 Be ages (Fig. 5).
Further north, ice margin recession throughout Massachusetts, Vermont and New Hampshire (~200-600 km from the terminal moraine; positions 7-11 in Fig. 3) was mostly continuous from ~18 to ~13 ka (Fig. 5), with a brief readvance around 14 ka, during which the White Mountain moraine system was formed (Thompson et al., 2017).The 10 Be-based calculations suggest that ice recession accelerated around 18 ka (up to ~70 m a −1 ) and maintained this faster rate through the northern part of our study area.However, the varve chronology (NAVC), which has a far higher resolution than the 10 Be chronology, indicates a more complex series of rate changes with readvance and standstill events (Ridge et al., 2012).
According to the NAVC, ice recession slowed from 17.4 to 16 ka (30-40 m a −1 ) before accelerating from 16 to 14.6 ka (80-90 m a −1 ) after the Chicopee Readvance that may have taken a few centuries, and then accelerating once again from 14.6 to 14 ka (~300 m a −1 ) after a delay during the construction of a few moraines near North Charlestown, New Hampshire (Fig. 5; Ridge et al., 2012).This final acceleration to ~300 m a −1 after 14.6 ka is far faster than the ~70 m a −1 rate suggested by the 10 Be-based recession model but fits within the range of rates suggested by both 14 C-based recession models (bulk sediment: 215-345 m a −1 ; macrofossils: 260-510 m a −1 ).

Implications for dating LIS recession
Inferences we draw from this data compilation have implications for the interpretations of organic 14 C and cosmogenic 10 Be ages linked to LIS recession and suggest two important characteristics of the LIS deglacial chronology in the northeastern US.First, ice recession chronometers (organic 14 C, cosmogenic 10 Be, and calibrated glacial varves) do not always agree, yielding substantial timing differences on the initial, post-LGM recession of the LIS.For example, most minimumlimiting 14 C ages are 3-10 kyr younger than nearby 10 Be exposure ages within tens of kilometres of the LIS terminal moraine zone.Second, this chronologic disagreement narrows with increased ice recession distance; all chronometers give similar results at locations more than 200 km north of the terminal moraine zone, thousands of years after deglaciation began.
The discrepancies we report here among calibrated organic 14 C ages, 10 Be exposure ages, and varve chronologies are not unique to the northeastern US.Similar chronometric disagreements, particularly organic 14 C ages trending younger than 10 Be exposure ages, have been reported at other former ice sheet margins, including elsewhere along the southern LIS (Lowell et al., 2021), the Cordilleran Ice Sheet (Briner and Swanson, 1998) and the Fennoscandian Ice Sheet (Briner et al., 2016;Stroeven et al., 2016).
The reduced variance of each geochronometer within populations of nearby ages in the northern part of our study area (Table 3; 'Deglaciation dates', above) suggests that methodological assumptions for both 10 Be exposure dating and organic 14 C dating become increasingly valid as the distance from the terminal position and the time since the LGM increases.Lower and more homogeneous inheritance of interglacially produced 10 Be is found up-ice from the terminal moraine zone, potentially the result of MIS 4 ice occupation (Dalton et al., 2022), and thus reduced pre-LGM surface exposure, and/or deeper glacial erosion, the result of both longer ice occupation times and thicker ice.Palaeoclimate conditions later in the deglacial period, when the LIS margin was in the northern part of our study area, were more favourable for rapid, abundant and denser revegetation of the deglaciated landscape (Shuman et al., 2004;Lora and Ibarra, 2019).This revegetation not only reduced the overall lag time between deglaciation and organic matter deposition, but also provided more organic matter in lake sediment for dating and thus more consistent core-bottom 14 C ages.
Better constraints on the timing of LIS deglaciation from the terminal moraine zone will require sampling from multiple sites and employment of a variety of chronometers.The data presented by Corbett et al. (2017) from New Jersey suggest that there are sites near the terminal moraine zone where cosmogenic and 14 C age estimates are coincident; such sites provide a promising opportunity for more confidently dating the maximum advance of the LIS in the northeastern US.Understanding potential methodological biases and their spatial occurrences will contribute to a better selection of reliable samples and sites, and ultimately reduce inaccuracies.
Improving 14 C dating of organic material in lake sediment will require finding the sparse organic material deposited before substantial revegetation began.It is possible that many previous hand-driven piston cores in the region did not penetrate basal lake sediments deeply enough to retrieve the oldest post-glacial organic material; electric or hydraulicassisted coring technologies could enable penetration of these stiff basal sediments (e.g.Breckenridge et al., 2012).The identification of threshold lakes, as has been done in Greenland (Briner et al., 2010), provides the opportunity to date ice advances from outside the glacial limit where more contemporaneous organic material may be present.
Accurate cosmogenic exposure dating of initial retreat from the LGM terminal moraine zone will require finding places where the ice sheet accomplished erosion sufficient to remove most cosmogenic nuclides formed during prior interglacial exposure.Using New Jersey as an example, jointed bedrock outcrops that facilitate plucking might be preferred over boulders or more uniform lithologies.Additionally, sample sites located in areas of directed ice flow may have experienced deeper erosion.Large populations of exposure age samples provide a more robust means to reject outliers as does the comparison of boulder and bedrock exposure ages.Alternatively, depth profiles of cosmogenic nuclide concentrations below glacially shaped outcrops can better constrain the concentration of inherited nuclides (induced by muons), providing a means to correct ages derived from outcrop samples.

Figure 1 .
Figure 1.Location map showing deglacial chronology samples, varve chronologies (Ridge et al., 2012; Stanford et al., 2020), place names and study area location relative to the whole Laurentide Ice Sheet (red box).State names are abbreviated for Pennsylvania (PA), New Jersey (NJ), Connecticut (CT) and Rhode Island (RI).[Color figure can be viewed at wileyonlinelibrary.com]

Figure 2 .
Figure 2. Ice recession dating methods.(A) Basal sand, gravel and organic matter that can be 14 C-dated in a lake sediment core from Lake Itasca, Minnesota (Photo credit: A. Noren, University of Minnesota).(B) Glacially deposited boulder for cosmogenic exposure dating located on Mount Mansfield, Vermont, USA (Halsted et al., 2022).Scale card is 15 cm long.(C) Rhythmic glacial lake sediments interpreted as varves, from Waterbury, Vermont, USA (Photo credit: C. Halsted).[Color figure can be viewed at wileyonlinelibrary.com]

Figure 3 .
Figure 3. Inset: Map of study area showing geomorphic features and ice margin positions identified in the text.Ice margin morphologies are from a recent reconstruction of LIS deglaciation (Dalton et al., 2020).Ice recession chronometer samples are the same as in Fig. 1.Main: Summed probability density functions for each dating method at selected ice margin positions (identified in inset map).Yellow triangles show the date of ice recession as indicated by varve chronologies at each ice margin position.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 4 .
Figure 4. Demonstration of summed probability density function (PDF) methodology (informally referred to as a 'camel plot').The dashed black curves are individual sample date PDFs with assumed Gaussian distributions.The larger, irregular filled brown curve is the sum of these individual age PDFs.[Color figure can be viewed at wileyonlinelibrary.com]

Figure 5 .
Figure 5. (A) All chronometer data from each technique (treated individually) available in the study area.Distances on the y-axis are measured along the Connecticut River valley (following the North American Varve Chronology) and samples are plotted at the distance where their assigned ice margin positions intersect the valley (Fig. 3).All 14 C ages are calibrated to calendar years BP and error bars encompass the 68.2% confidence interval on calibrated ages.Error bars on 10 Be ages are 1-sigma external (measurement and production rate estimate) uncertainty.Two outliers with exposure ages >50 ka on the Martha's Vineyard moraine not included in this figure.Dashed lines show the most likely recession chronology for each chronometer from Bayesian modelling (Blaauw and Christen, 2011) and coloured bands around dashed lines show the 95% confidence interval.Italicised numbers show most likely rates of ice recession (m a −1 ) with 95% confidence intervals in parentheses.Shaded bars in the background of the figure indicate the Oldest Dryas (OD), Bølling-Allerød (BA), Younger Dryas (YD) and Holocene (H) intervals.(B) 21 July insolation curve for 40°N (Laskar et al., 2011).(C) Greenland temperature reconstruction from the NGRIP ice core (Kindler et al., 2014).[Color figure can be viewed at wileyonlinelibrary.com]

Figure 7 .
Figure 7. Depth profile of 10 Be concentrations (age equivalent) in a bedrock surface with 110 kyr of exposure.Brown squares show a 1 m 3 boulder plucked from 1.5-2.5 m depth, 11-12 m depth and 19-20 m depth.These brown squares represent boulder exhumation depths expected near the LGM terminal moraine, in southern New England and in northern New England, respectively.Numbers indicate the ageequivalent 10 Be inheritance on the top and bottom of these boulders (see southern New England example).[Color figure can be viewed at wileyonlinelibrary.com]

Table 1 .
Laurentide Ice Sheet advance constraints in or near the northeastern United States.

Table 2 .
Previously published Laurentide Ice Sheet deglacial chronologies in the northeastern US.