The maximum extent and elevation of the Greenland Ice Sheet in southwestern Greenland during the Last Glacial Maximum (LGM, 26–19.5 ka) is poorly constrained. Yet, the size of the Greenland Ice Sheet during the LGM helps to inform estimates of past ice-sheet sensitivity to climate change and provides benchmarks for ice-sheet modeling. Reconstructions of LGM ice extents vary between an inner continental shelf minimum, a mid-shelf position, and a maximum extent at the shelf break. We use three approaches to resolve LGM ice extent in the Sisimiut sector of southwestern Greenland. First, we explore the likelihood of minimum versus maximum Greenland Ice Sheet reconstructions using existing relative sea-level data. We use an empirical relationship between marine limit elevation and distance to LGM terminus established from other Northern Hemisphere Pleistocene ice sheets as context for interpreting marine limit data in southwestern Greenland. Our analysis supports a maximum regional Greenland Ice Sheet extent to the shelf break during the LGM. Second, we apply a simple 1-D crustal rebound model to simulate relative sea-level curves for contrasting ice-sheet sizes and compare these simulated curves with existing relative sea-level data. The only realistic ice-sheet configuration resulting in relative sea-level model-data fit suggests that the Greenland Ice Sheet terminated at the shelf break during the LGM. Lastly, we constrain the LGM ice-sheet thickness using cosmogenic 10Be, 26Al, and 14C exposure dating from two summit areas, one at 381 m above sea level at the coast, and another at 798 m asl 32 km inland. Twenty-four cosmogenic radionuclide measurements, combined with results of our first two approaches, reveal that our targeted summits were likely ice-covered during the LGM and became deglaciated at ca. 11.6 ka. Inventories of in situ 14C in bedrock at one summit point to a small degree of inherited 14C and suggest that the Greenland Ice Sheet advanced to its maximum late Pleistocene extent at 17.1 ± 2.5 ka. Our results point to a configuration where the southwestern part of the Greenland Ice Sheet reached its maximum LGM extent at the continental shelf break.

Reconstructions of past ice-sheet dimensions serve as important prescribed features in transient climatic model simulations (He et al., 2013; Zweck and Huybrechts, 2005) and are also used to independently assess paleoice–sheet model performance (Downs et al., 2020; Lecavalier et al., 2014; Tarasov and Peltier, 2004). Yet, such reconstructions reveal conflicting estimates of the maximum size of the Greenland Ice Sheet during the Last Glacial Maximum (LGM, 26–19.5 ka) around some parts of Greenland. For example, while most reconstructions of Greenland Ice Sheet extent during the LGM depict a terminus located somewhere on the continental shelf, its exact position is uncertain (Fig. 1; Bradley et al., 2018; Fleming and Lambeck, 2004; Funder et al., 2011; Lecavalier et al., 2014). Placing direct constraints on the LGM ice extent requires offshore geophysical surveys, seafloor mapping, and sediment coring; however, only a few areas have yielded detailed knowledge about past ice extent (Dowdeswell et al., 2014; Hogan et al., 2016; Jennings et al., 2017; Newton et al., 2017; Ó Cofaigh et al., 2013; Ó Cofaigh et al., 2004).

Offshore LGM ice extent can also be inferred using terrestrial evidence. Glacial sediments and landforms can be used to constrain former ice extent, but these archives rarely extend back into the LGM on Greenland (Håkansson et al., 2007; Larsen et al., 2018; Young et al., 2020). Relative sea-level curves have been used in combination with glacially weathered trimlines to infer the size and lateral extent of ice sheets, but these estimates are typically used as a minimum ice-sheet estimate (Funder and Hansen, 1996). A growing number of Greenland's coastal regions have been studied using cosmogenic nuclide exposure dating, which provides absolute age constraints on the timing of ice retreat onto land (e.g., Winsor et al., 2015; Roberts et al., 2009; Søndergaard et al., 2020). However, deglaciation chronologies are less straightforward when used in areas where frozen bedded conditions dominated during the LGM, such as in eastern and northern Greenland, or elsewhere at higher elevations (Briner et al., 2013; Corbett et al., 2013; Graham et al., 2019; Kelly, 1980).

Southwestern Greenland has an extensive terrestrial record of the Greenland Ice Sheet during the Holocene (Young et al., 2020), but like other regions in Greenland, its glacial extent during the LGM is unclear. A commonly cited reconstruction places the southwestern Greenland Ice Sheet LGM terminus a mere 40 km offshore of the modernday coastline (Funder et al., 2011), yet studies on the wide continental shelf off Disko Bugt support an LGM ice limit at the shelf break (Ó Cofaigh et al., 2013). Terrain in southwestern Greenland above ~650 m above sea level (asl) is often heavily weathered and has apparent 10Be exposure ages that far exceed the LGM period (Kelly and Bennike, 1985; Rinterknecht et al., 2009; Roberts et al., 2009; Winsor et al., 2015). These results have been interpreted to indicate that the surface of the LGM may not have been more than 650 m asl, which would mean that some mountain summits survived as nunataks during the LGM.

The exposure history of summits interpreted to be potential nunataks could be complicated by non-erosive LGM ice cover. For example, weathered uplands, even those with pre-LGM cosmogenicnuclide inventories in bedrock, and sometimes even in erratics, have been shown to be at least intermittently covered by ice during the LGM (Briner et al., 2006; Corbett et al., 2013; Young et al., 2018). In south-western Greenland, if coastal uplands were covered by ice during the LGM, it would point to a terminus position far onto the continental shelf and possibly at the shelf break. These issues collectively have led to an uncertain reconstruction of the lateral extent and maximum elevation of the Greenland Ice Sheet during the LGM in some areas around Greenland.

Here, we provide new evidence to support an extensive Greenland Ice Sheet configuration in southwestern Greenland during the LGM. We investigate the relationship between the elevation of marine limits and their locations relative to LGM ice extent. Next, we develop a 1-D crustal rebound model for a model-data comparison of relative sea-level history. Finally, we constrain ice thickness with cosmogenic nuclide measurements from erratic and bedrock surfaces atop coastal summits in south-western Greenland. Combined, our data suggest that the southwestern Greenland Ice Sheet terminus extended to the shelf break during the LGM.

2.1. Marine Limit Data Compilation

Raised marine deposits in deglaciated environments provide valuable information for understanding former ice-sheet loading and post-glacial uplift history (Peltier, 1998). Early work showed that marine limit elevation (U) is generally related to its distance from the LGM ice-sheet margin (Andrews, 1968a, 1968b; Andrews et al., 1970; Andrews, 1975). This relationship is:
where the independent variable x is the distance from the ice center, C is a constant with a value near 1, L is the distance from the ice margin to the center of the ice sheet during the LGM (Andrews, 1968a), τav is the average basal shear stress under an ice sheet with density ρice = 917 kg/m3, g is the gravitational acceleration, and n is an exponent with value near 0.5 (Andrews, 1968a). We fit n in Equation 1 to the marine limit data using a power trendline.

We explore the relationship between marine limit elevation and its distance to ice-sheet margins by compiling marine limit data from around the Laurentide, Fennoscandia, Barents–Kara, and Greenland Ice Sheets. We use LGM ice-sheet outlines from published compilations (Funder et al., 2011; Dalton et al., 2020; Hughes et al., 2016) and include 62 marine limit observations from northeastern Canada (Dyke et al., 2005), 55 from Greenland, 22 from Svalbard (Forman, 1990), and 42 from Ireland, the UK, and Scandinavia (Hogan et al., 2016) (Fig. 2). We measure the distance from each marine limit location to the mapped LGM ice-sheet margin. In places where the mapped LGM extent did not reach the shelf break, we made a second measurement from the marine limit location to the shelf break.

2.2. One-Dimensional Crustal Rebound Model and Simulated Relative Sea-Level Curves

We compare observations of marine limit elevations and relative sea level with results from a simple model run with contrasting Greenland ice-loading scenarios. To do so, we applied a simple 1-D crustal rebound model to incorporate ice-sheet profiles calculating the ice-sheet thickness, H, at a given distance, x, from the margin (van der Veen, 2013):
where τy is the yield stress. Discretization of Equation 2 using the forward Euler method along a roughly flow-parallel transect allows calculation of the ice-sheet surface elevation, hi, at discrete positions, xi, from the margin, where i represents the index of discretization and Δx the spatial step:
following Benn and Hulton (2010) and using bed elevations B from Morlighem et al. (2017). We constrained the yield stress by applying Equation 3 to a transect from the modern ice divide to the modern terminus and tuning the yield stress until the surface elevation profile h(x) best matched the current ice-sheet surface profile from Morlighem et al. (2017). We found a best-fit yield stress of 110 kPa, which we employed at all paleo-model times. This model represents the Greenland Ice Sheet as a land-terminating ice sheet, although in both paleo scenarios, it was likely marine-terminating (Fig. 3).

We modeled two versions of the ice sheet, one with an LGM terminus at the shelf break (“Max Ice”) and another with the LGM terminus at the location at the Fiskebanke moraines on the inner shelf (“Min Ice” in Fig. 3). We enforced these configurations from 30 ka until 15 ka, which was long enough to relax the underlying crust into equilibrium with the overlying ice load. We applied the same retreat history to both the Max Ice and Min Ice configurations. We initiated deglaciation at 15 ka based on the timing of retreat from the shelf break farther north (e.g., Ó Cofaigh et al., 2013) and then used five mapped Greenland Ice Sheet margin positions for Holocene time steps at 11.6 ka, 10.4 ka, 9.1 ka, 8.1 ka, and 7.3 ka based on data from Young et al. (2020) and Lesnek et al. (2020), with linear retreat rates between each ice margin position. There is evidence that the Greenland Ice Sheet retreated inland of its current position during the middle Holocene. Although the exact distance remains unknown, it likely was no more than a few tens of kilometers from its current position (Young et al., 2021; Lesnek et al., 2020; Ruskeeniemi et al., 2018; Young and Briner, 2015). We represented this by positioning the margin at the modern terminus at 5.5 ka, then 20 km inland of the modern terminus from 4.0 ka to 1.2 ka, then at its modern position at 0 ka, with linear rates of margin migration between these defined times.

For each configuration of the ice sheet at the LGM, we subtracted the local thickness of the modern ice sheet to obtain a thickness change, ΔH, along our transect. We used ΔH to calculate the visco-elastic displacement of the crust during the LGM using the modern topography as a baseline (Morlighem et al., 2017). Our simple model represents the lithosphere as an elastic beam floating in the mantle, following the classic bending beam problem modeled with the biharmonic equation that was outlined and solved by Turcotte and Schubert (2002). At every model timestep, we calculated the elastic displacement of the crust we(x):
where D is the flexural rigidity parameter for the lithosphere, defined below, which is weighted by a spatially varying line load, v(x), that represents the ice sheet:
where we represent the ice sheet as a series of line loads of width Δx = 500 m, spatially varying thickness, H(x), and ice density, ρice, of 917 kg/m3. The horizontal length scale of the forebulge, α, is defined as:
where mantle density, ρm, is 3200 kg/m3. The flexural rigidity parameter D is:
where the Young's modulus of the elastic lithosphere E is 70 GPa (Turcotte and Schubert, 2002), the lithospheric elastic thickness d is 80 km (Steffen et al., 2018), and Poisson's ratio ν for ice is 0.25 (Turcotte and Schubert, 2002). Our crustal model produces a maximum subglacial depression (~−630 m), a slight proglacial depression immediately beyond the terminus (~−130 m), and a forebulge (~+15 m) at a distance of α ~150 km beyond the terminus.
We also model the viscous response of the lithosphere to the loading and unloading by advance and retreat of the ice sheet. The timescale of this response, Tr, is:
where µ is the viscosity of the upper mantle (5 × 1020 Pa • s) (Wake et al., 2016) and L is the length of the LGM ice sheet (Meinesz, 1937). With these physical constants and the differing ice-sheet lengths L in the “Max Ice” and “Min Ice” configurations, we find that Tr for the “Max Ice” model is 4300 yr and Tr for the “Min Ice” model is 6200 yr.
We model the time evolution of the viscous displacement of the crust, w(t), as exponential decay:
where wmax is the depression of the bed during the LGM. This model allows the crust to rebound as soon as the ice sheet begins to thin, even when still covered in ice (i.e., restrained rebound).

To validate this simple model for visco-elastic deformation of the lithosphere, we ran the Benchmark A comparison test described in Martinec et al. (2018) and Spada et al. (2011). Our modeled vertical displacement under the outer 70% of the glaciated area agreed with the Benchmark A results (Martinec et al., 2018) to within 10 m (~10%).

Eustatic sea level increased simultaneously with ice-sheet retreat (Lambeck et al., 2014). We calculate the relative sea level at a site at a given time, R(t), over the entire model duration (pre-15 ka to present) as follows:
where B0 is the modern crust elevation (Morlighem et al., 2017), B(t) is the elevation of the crust at a given time, and S(t) is global sea level at a given time (Lecavalier et al., 2014). This allows us to simulate a relative sea-level curve for a selected site spanning from deglaciation to today (Fleming et al., 1998).

2.3. Cosmogenic-Nuclide Exposure Dating

Three previous studies with cosmogenic nuclide measurements are relevant to our work because they include 10Be ages in the Sisimiut area. Here, we recalculated published ages using the methods described below. Samples collected from bedrock above 650 m asl yielded some ages in excess of 100 ka (Roberts et al., 2009), and samples from lower elevations have 10Be ages that are >15 ka and thus could be interpreted to indicate ice-free LGM conditions or early deglaciation (Rinterknecht et al., 2009; Roberts et al., 2009). Roberts et al. (2009) interpreted a boundary at ~650 m asl as a possible maximum LGM ice-sheet elevation for the region and supported a limited maximum lateral extent with the LGM terminus placed just offshore (Kelly and Bennike, 1985; Weidick and Bennike, 2007). Winsor et al. (2015) reported ages of ca. 35 ka from boulders at ~600 m asl, which could be interpreted to indicate an even lower LGM ice surface. Winsor et al. (2015) also reported 10Be ages of boulders that average ca. 14.5 ± 1.0 ka from ~150 m asl and suggested that deglaciation of the Sisimiut area took place at this time. The ca. 14.5 ka age of deglaciation in Winsor et al. (2015) contrasts with a single 10Be age of ca. 20 ka from a boulder at ~400 m asl only ~10 km away reported by Rinterknecht et al. (2009). Finally, Young et al. (2020) reported three 10Be ages from ~230 m asl that average 11.5 ± 0.3 ka, which is notably younger than the nearby ages from Winsor et al. (2015). All of these previously dated sites lie to the west of the Taserqat Moraine, which Young et al. (2020) dated to 11.6 ± 0.4 ka (n = 5).

To further explore the unresolved nature of upland glaciation and former Greenland Ice Sheet thickness in the Sisimiut area, we use multiple nuclides, 10Be, 26Al, and in situ 14C, due to their combined utility in elucidating glacier history in areas with complex exposure histories (e.g., Beel et al., 2016; Bierman et al., 1999; Graham et al., 2019; Young et al., 2018). We collected three bedrock and six boulder samples from a summit adjacent to the ocean with sample elevations ranging from 355 m asl to 397 m asl (Figs. 45). The summit at this “western site” displays evidence of past ice cover due to the presence of glacial erratics. Signs of advanced subaerial weathering, such as protruding quartz veins with ~1–2 cm of relief above the surrounding bedrock, are also visible at this site, which indicates a potentially complex exposure history. We also collected three bedrock samples from a summit at 798 m asl located 32 km inland from the western site (Figs. 45). The summit bed-rock surfaces at this “eastern site” exhibit signs of advanced weathering with deep weathering pits; no boulders suitable for sampling were present.

We collected samples using a rock hammer, chisel, and angle grinder from the top ~2 cm of each surface. The boulders were perched in positions that imply stability since deposition. Bedrock samples were collected from local topographic high points to minimize the effect of shielding from seasonal snow cover and surrounding topography.

We processed samples using a series of physical and chemical treatments. At the University at Buffalo (Buffalo, New York, USA), the rock samples were crushed and then sieved to 250–500 µm. Quartz was isolated using froth flotation in an acetic acid solution followed by hydrochloric and hydrofluoric/nitric acid baths. We dissolved clean quartz with beryllium carrier with a concentration of 1074 ppm of 9Be and added 26Al carrier (1001 ppm) to the bedrock samples (Corbett et al., 2016). Once beryllium and aluminum were isolated and purified, they were oxidized and loaded into cathodes for accelerator mass spectrometry (AMS) measurement at PRIME Lab at Purdue University (West Lafayette, Indiana, USA). 10Be and 26Al ages were calculated with the CRONUS v3 online age calculator (, and 10Be and 26Al ages were calculated using the Baffin Bay production rate of Young et al. (2013) using the Lm scaling scheme (Balco et al., 2008; Lal, 1991; Stone, 2000). We also processed six quartz separates for in situ 14C analysis at Lamont-Doherty Earth Observatory following standard procedures detailed in Lamp et al. (2019). To avoid potential carbon contamination from the froth flotation procedure at the University at Buffalo, we leached samples in concentrated nitric acid for several days before processing them further. The measured 14C concentrations were blank-corrected using a long-term laboratory blank (Lamp et al., 2019). In situ 14C ages were calculated using the western Greenland (14C) production rate calibration data set (Young et al., 2014) and Lm scaling (Balco et al., 2008; Lal, 1991; Stone, 2000).

3.1. Marine Limit Data Complication

Our marine limit analysis shows that marine limit elevation is significantly related to distance from the hypothesized LGM margin (Fig. 6). We fit data for all four ice sheets to Equation 1 and found r2 = 0.73 with C = 4.7 and n = 0.50 in agreement with the theoretical value of n = 0.5 and C of order 1 (Andrews, 1968a). Some of the scatter, we assume, is due to the added influence that ice margin history during deglaciation (variable retreat rate, stillstands, and re-advances) has on the relationship. Of particular interest are the data from southwestern Greenland using the Greenland Ice Sheet ice limit from Funder et al. (2011), which fall outside the trend of data from all other ice sheets (Fig. 6). On the other hand, when using an LGM margin at the shelf break, the best-fit equation resembles that of other ice sheets much more closely.

The LGM reconstruction of Greenland Ice Sheet extent in the Sisimiut region proposed by Dyke (2004) and Funder et al. (2011) results in a marine limit versus distance relationship that differs from that of other Northern Hemisphere ice sheets (Fig. 6) and requires a relationship different from the theoretical parabolic profile (Equation 1). On the other hand, our “Max Ice” scenario, in which the ice-sheet terminus reaches the shelf break, the marine limit versus distance plot becomes similar to the results from other ice sheets (Fig. 6) and fits well with a parabolic profile. These findings support an LGM extent at the shelf break off Sisimiut.

3.2. Crustal Rebound Model

The results from our crustal rebound simulations (see animations in Supplemental Material1) show relative sea-level histories similar to those from independent observations. Our simulated relative sea-level curve using “Max Ice” and a deglaciation age of 12 ka yields a marine limit for Sisimiut of 136 m asl. This is close to the observed marine limit of 140 m asl and generally compatible with field constraints of relative sea level (Bennike et al., 2011; Funder and Hansen, 1996; Long et al., 2011). The simulated marine limit using “Min Ice” is 47 m asl and has a poorer fit with observations (Fig. 7). Thus, the overall better fit between our simulated “Max Ice” relative sea-level history and observations favors a “Max Ice” scenario.

Despite the agreement between geological observations and the results from our “Max Ice” crustal rebound simulations, the ice-sheet history we use in our model is not well constrained prior to 11.6 ka. In our “Max Ice” scenario, we force the ice margin to retreat from the continental shelf break at 15 ka (see Section 2.2) and then prescribe a linear retreat rate over the next 4.4 k.y. until the ice margin reaches the next recorded position ~115 km inland at 11.6 ka. In our “Min Ice” scenario, we also begin retreat at 15 ka at a constant rate until 11.6 ka, but the ice margin begins much farther inland (Fig. 3). If the Greenland Ice Sheet had a complex ice-margin history with cycles of repeated advance and retreat (e.g., retreat from the LGM limit during the Bølling followed by re-advance and final retreat after the Younger Dryas), such an ice-margin history could lead to a different loading history. Similarly, our assumption that the crust achieved equilibrium with the LGM ice load may not be valid if the Greenland Ice Sheet was at its maximum extent for a relatively short period, such as arriving to its LGM limit not long before it receded from that limit, a result that is supported by our cosmogenic nuclide results (see Section 3.3). If the maximum LGM ice extent were reached relatively late and occupied for only a brief period of time, then the Greenland Ice Sheet would need to have been thicker than in our “Max Ice” simulation to yield a marine limit that matches independent observations. In any case, despite these uncertainties, our simplified ice-sheet history for “Max Ice” nevertheless produced a reasonably close fit with observations (Fig. 7).

We also explored how thick the Greenland Ice Sheet would have to be in our “Min Ice” scenario to reproduce the observed 140 m asl marine limit near Sisimiut. This requires an ice thickness of ~1400 m at the modern coastline, which would require an increase in the yield stress from 110 kPa to 310 kPa to maintain an inner-shelf terminus position. In addition to being a factor of three higher than the typically accepted value (van der Veen, 2013), a yield stress of 310 kPa results in an ice-sheet thickness over 5 km at the ice divide. We find it unlikely that a minimum ice configuration could result in a reasonable data-model fit.

3.3. Cosmogenic Nuclide Concentrations

From the western site, three bedrock samples have 10Be ages of 24.0 ± 0.6 ka, 18.6 ± 0.5 ka, and 15.8 ± 0.4 ka (Tables 12). The five samples from erratics at the western site have 10Be ages of 14.7 ± 0.5 ka, 11.8 ± 0.4 ka, 11.7 ± 0.4 ka, 11.0 ± 0.4 ka, and 10.3 ± 0.3 ka. The 26Al ages at the western site are in good agreement with the 10Be ages, and the 26Al/10Be ratios of the samples range from 6.77 ± 0.26–7.78 ± 0.31 (Table 1). We measured in situ 14C in four samples from the western site (one bedrock and three erratics) to further constrain the exposure history of the site (Table 2). In situ 14C ages from two erratics are 11.1 ± 1.1 ka and 11.5 ± 1.4 ka, and the bedrock sample has an age of 11.7 ± 1.2 ka. One of the in situ 14C measurements from an erratic resulted in a concentration of ~210 k atoms/g, which is higher than the site-specific saturation value of ~190 k atoms/g. We suspect this sample was contaminated during the extraction process and disregard it from further consideration.

Two of the erratic samples (18GRO-6 and 18GRO-11) at the western site have 10Be ages that are beyond the 1σ age range of the others (one older, one younger; Fig. 4), and the average age of the remaining three erratics is 11.5 ± 0.4 ka. We interpret this age of 11.5 ± 0.4 ka as the timing of deglaciation at the western site.

The in situ 14C ages for the samples from the western site strengthen the case for the 11.5 ± 0.4 ka age best representing the timing of deglaciation. Sample 18GRO-9, an erratic with a 10Be age of 11.0 ± 0.4 ka, has an in situ 14C age of 11.1 ± 1.1 ka. Both 18GRO-4 (bedrock) and 18GRO-11 (erratic) have 10Be and 26Al ages (ca. 15–17 ka) that are older than our interpreted age of deglaciation (ca. 11.5 ka), which we suspect is due to inheritance. The in situ 14C ages for these same samples (11.7 ka and 11.5 ka, respectively) are within 1σ of our erratic 10Be-based age of deglaciation of 11.5 ± 0.4 ka. The average of all three in situ 14C ages from the western site is 11.5 ± 0.3 ka, the same as our average age of deglaciation based on 10Be ages of erratics (11.5 ± 0.4 ka). Overall, these ages support a deglaciation of the western site around the Younger Dryas/Holocene boundary and imply ice cover prior to that time.

At the eastern site, the three bedrock samples (18GRO-1A, 18GRO-, and 18GRO-2) have 10Be ages of 79.4 ± 1.3 ka, 63.3 ± 1.0 ka, and 69.0 ± 1.7 ka, with corresponding 26Al ages of 66.3 ± 2.9 ka, 56.1 ± 1.1 ka, and 55.7 ± 1.4 ka (Table 1; Fig. 4). The three samples have 26Al/10Be ratios ranging from 5.85 ± 0.20–6.40 ± 0.16 ka (Table 1). We also measured in situ 14C in two of the bedrock samples from the eastern site, which are 17.1 ± 2.5 ka and 19.0 ± 3.2 ka (Table 2; Fig. 4).

The cosmogenic nuclide data from the eastern site exhibit evidence for a complex burial and exposure history. The 26Al and 10Be ages from the three bedrock surfaces are far older than the regional timing of deglaciation (e.g., Young et al., 2020). The samples have 26Al/10Be ratios that are 6.10 ± 0.13, significantly lower than 7.3, which has been proposed as an appropriate constant exposure value for high latitudes (Corbett et al., 2017; Halsted et al., 2021). This indicates that these samples experienced episodes of burial, upwards of 250 k.y., likely by overriding ice (Bierman et al., 1999; Beel et al., 2016; Strunk et al., 2017). These durations increase when using a ratio of 7.3. These 26Al/10Be ratios, which are similar to values that others have found in western Greenland uplands (Roberts et al., 2009; Beel et al., 2016; Corbett et al., 2017), do not allow us to definitively say when this burial occurred, but it is possible that a portion of the burial occurred during the LGM (Beel et al., 2016).

Without further information, whether the eastern site was covered during the LGM remains unknown. However, the apparent in situ 14C ages provide additional insight. Samples 18GRO-1a and 18GRO-2 are much younger than their 10Be and 26Al ages and average 18.1 ± 1.3 ka. There are two ways that these results could be interpreted. First, the in situ 14C ages could represent the timing of deglaciation of the eastern site. Second, the ages could contain some amount of inherited 14C due to only a brief, minimally erosive overriding event during the LGM, in which case the site was deglaciated sometime after ca. 18.1 ka. In the next section, we explore these interpretations in light of the larger data set from both our data from the western site and prior work.

4.1. Support for a More Extensive Greenland Ice Sheet in Southwestern Greenland during the LGM

The cosmogenic nuclide ages at the western site agree with 14C-dated sediments collected from lakes in the region. Our ages are older than the basal lake 14C sediment age of 10.5–11.0 cal ka B.P. (re-calibrated) that provides a minimum limiting deglaciation age for the region (Bennike et al., 2011). The deglaciation of the western site also agrees with the average of three 10Be ages from erratic boulders (11.5 ± 0.3 ka) collected ~10 km south of Sisimiut reported by Young et al. (2020) (Fig. 4). The 26Al/10Be ratios from the western site samples reveal insignificant burial even considering the two pairs from bedrock with low levels of inheritance (Fig. 8). Based on this, we interpret the age of deglaciation for the Sisimiut coast to be ca. 11.5 ka.

We suggest that both the western and eastern sites were deglaciated at around the same time, in which case the in situ 14C concentrations are likely due to (1) low-magnitude erosion (if any) of the bed due to non-erosive ice conditions during the LGM and/or (2) an insufficient duration of ice cover to result in complete 14C decay. The 11.6 ka Taserqat Moraine is draped across the massif of the eastern site (Fig. 4), a mere ~100 m lower in elevation than the summit and only ~1 km to the SE. Given this geometry, we find it unlikely that the eastern site would have remained a nunatak during the LGM. Moreover, considering our exposure ages from the western site and across the general field area, we suggest our eastern site likely deglaciated just prior to Taserqat moraine deposition. Both the western site and the island site of Young et al. (2020) deglaciated at ca. 11.5 ka, the same age as the Taserqat Moraine (11.6 ± 0.4 ka), which implies rapid deglaciation of the region at this time. Thus, the Greenland Ice Sheet likely retreated fast enough to negate any difference in ice thickness at our two study sites, causing them to become ice-free at roughly the same time. Finally, although higher in elevation, the eastern site is farther inland, and generalized ice-surface profiles indicate that both sites would deglaciate at similar times. LGM ice probably would have been sufficiently thick to prevent continued cosmogenic nuclide accumulation through ice during the LGM, a phenomenon that has been observed at thin ice caps elsewhere in the Arctic (Pendleton et al., 2019).

Our exposure ages agree with those of Young et al. (2020) and, combined, suggest that the Sisimiut region deglaciated at ca. 11.5 ka. We note that one of the boulders that we dated (18GRO-9) yields a 10Be exposure age of 11.0 ± 0.4 ka, but a previous study reported a 10Be exposure age of 19.6 ± 5.2 ka from the same boulder (GRE-1; Rinterknecht et al., 2009); however, the two ages overlap at 2 sigma. The cause of the age discrepancy may be due to improving laboratory methods like minimizing background 10B isobaric interference and using improved 10Be ion exchange chromatography for increased sample purity (Corbett et al., 2016). Winsor et al. (2015) collected nine samples from boulders for 10Be exposure dating at two sites close to Sisimiut. Four of five samples collected at ~585 m asl yield an average age of 36.2 ± 2.1 ka with one younger age of 17.9 ± 0.8 ka. Winsor et al. (2015) raised the possibility that the site remained above the LGM ice surface; in light of our results, it seems that their alternative interpretation is more likely, which is that the samples contain inheritance and the site was covered by LGM ice. Winsor et al.'s (2015) four additional 10Be ages, each on a boulder from ~150 m asl, average 14.4 ± 0.9 ka after excluding one young age (10.2 ± 0.6 ka), which they interpreted as an outlier. These ages cannot be easily explained because they are at odds with our results and those of Young et al. (2020), but our results and the relatively younger ages from Young et al. (2020) suggest that small amounts of inherited nuclides may explain these slightly older ages.

Finally, we use likely timing of deglaciation (11.5 ka) to invert the in situ 14C concentration for a plausible duration of LGM ice coverage of the eastern site. Generally, in situ 14C becomes either fully saturated when constantly exposed or has decayed to a non-detectable level when fully buried at ~30 k.y. (Graham et al., 2019). Offshore evidence supports a glaciation in southwestern Greenland during Marine Isotope Stage 4 (MIS 4, 75–58 ka) that was more extensive than current LGM reconstructions suggest (Seidenkrantz et al., 2019). Thereafter, the Greenland Ice Sheet is generally thought to have been about as extensive as its current size throughout MIS 3 (58–27 ka) (England, 1999; Larsen et al., 2018). Thus, a plausible glacial history of the eastern site is that it became deglaciated at ca. 58 ka from an MIS 4 overriding event and subsequently remained ice free until the LGM. Following Graham et al. (2019), we model the in situ 14C concentration in the samples from the eastern site assuming they were saturated with respect to 14C upon arrival of the ice sheet during the LGM. Next, by assuming no erosion during ice cover and re-exposure of the surfaces at 11.5 ka, we calculate the onset of ice-margin advance and surface burial to be between 19.6 ka and 14.6 ka (Fig. 9). Regardless of the exact timing and duration of LGM ice cover, our in situ 14C measurements that yield finite exposure ages below the saturation limit mandate that these bedrock surfaces were not exposed throughout the LGM.

The new results presented here point to a more extensive LGM configuration in southwestern Greenland than is currently depicted in most literature. Our comparison of marine limit values from southwestern Greenland to those from other sectors of the ice sheet suggests that an inner-shelf LGM position is incompatible with the relatively high marine limits observed at coastal sites. Our simplified crustal rebound modeling also reveals that a shelf-break LGM configuration is more compatible with existing relative sea-level data than an ice sheet that terminated on the inner shelf. Finally, our new series of cosmogenic nuclide exposure ages supports ice cover of coastal summits and a final deglaciation age for the coast that occurred no earlier than the late Pleistocene/Holocene boundary. Thus, we favor the view that the Greenland Ice Sheet likely terminated at the continental shelf break during the LGM, with coastal deglaciation occurring just prior to ca. 11.6 ka.

Our conclusion is consistent with recent studies that place the LGM terminus position at the shelf break in the Disko Bugt and Uummannaq Fjord regions (e.g., Ó Cofaigh et al., 2013). Furthermore, the results highlight the need for direct sediment dating of the continental shelves around Greenland. Specifically, if the depositional age of the Helle-fiske or Fiskebanke Moraines were known, then it would be possible to constrain the age and location of the terminus of the Greenland Ice Sheet throughout southwestern Greenland more accurately. Our evidence from southwestern Greenland joins other work suggesting that shelf break termini may have been common throughout Greenland (Fig. 1A; England, 1999; Evans et al., 2009; Larsen et al., 2018; Ó Cofaigh et al., 2013; Roberts et al., 2009; Vasskog et al., 2015). Additional offshore studies and further applications of the modeling approaches taken here could test this hypothesis. If the Greenland Ice Sheet terminated at the shelf break throughout Greenland, this maximum extent would be consistent with the LGM extents of other major coastal ice sheets (Hughes et al., 2016; Mangerud et al., 1998; Margold et al., 2018).

Whereas our findings indicate that the coastal mountain summits in southwestern Greenland near Sisimiut were covered during the last glaciation, simple modeling of in situ 14C inventories suggests the duration of ice overriding was only 5.5 ± 2.5 k.y. starting at 17.1 ± 2.5 ka. Such a late maximum in ice volume is consistent with some ice-sheet modeling studies (e.g., Lecavalier et al., 2014; Buizert et al., 2018) that suggest a late LGM ice phase, perhaps correlative with Heinrich Stadial 1 (ca. 18.5–15 ka). Nevertheless, our cosmogenic nuclide measurements, modeled relative sea-level curves, and modeled ice-sheet profiles suggest that the Greenland Ice Sheet terminated at the shelf break at some point during the LGM or shortly thereafter in south-western Greenland.

We acknowledge Joseph Tulenko for assistance with generating maps of the study location. We thank T. Tuna and E. Bard for careful gas-source measurements at Centre de Recherche et d'Enseignement de Géosciences de l'Environnement (CEREGE). We acknowledge support by National Science Foundation Arctic System Sciences grants ARC-1504267 to J.P. Briner and E.K. Thomas and ARC-1503959 to N.E. Young. We extend thanks to Polar Field Services for field logistical support, and to two anonymous reviewers who took time to critically review this work.

All MATLAB code and input data used to model the viscoelastic relaxation of the lithosphere for comparison to relative sealevel data is available at

1Supplemental Material. Animations of ice-sheet recession and crustal rebound for minimum and maximum LGM ice sheet configurations. Please visit to access the supplemental material, and contact with any questions.
Science Editor: Andrea Hampel
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.