Spatial variability of marine-terminating ice sheet retreat in the Puget Lowland

Understanding drivers of marine-terminating ice sheet behavior is important for constraining ice contributions to global sea level rise. In part, the stability of marine-terminating ice is influenced by solid Earth conditions at the grounded-ice margin. While the Cordilleran Ice Sheet (CIS) contributed significantly to global mean sea level during its final post-Last-Glacial-Maximum (LGM) collapse, the drivers and patterns of retreat are not well constrained. Coastal outcrops in the deglaciated Puget Lowland of Washington State – largely below sea level during glacial maxima, then uplifted above sea level via glacial isostatic adjustment (GIA) – record the late Pleistocene history of the CIS. The preservation of LGM glacial and post-LGM deglacial sediments provides a unique opportunity to assess the variability in marine ice sheet behavior of the southernmost CIS. Based on paired stratigraphic and geochronological work, with a newly developed marine reservoir correction for this region, we identify that the late-stage CIS experienced stepwise retreat into a marine environment between 15 000 and 14 000 years before present, consistent with timing of marine incursion into the region reported in earlier works. Standstill of marine-terminating ice for at least 500 years, paired with rapid vertical landscape evolution, was followed by continued retreat of ice in a subaerial environment. These results suggest rapid rates of solid Earth uplift and topographic support (e.g., grounding zone wedges) stabilized the ice margin, supporting final subaerial ice retreat. This work leads to a better understanding of shallow-marine and coastal-ice-sheet retreat and is relevant to sectors of the contemporary Antarctic and Greenland ice sheets and marine-terminating outlet glaciers.

1 Introduction

The terrain and substrate geology beneath ice sheets have the potential to affect the behavior of the overriding ice; they can influence ice flow organization, velocity, and margin positions (Weertman, 1974; Clarke et al., 1977; Clark, 1994; Whillans and van der Veen, 1997; Cuffey and Paterson, 2010; Jamieson et al., 2012; Margold et al., 2015). Coupled ice sheet and solid Earth models indicate that glacial isostatic adjustment (GIA) can stabilize marine-based grounding lines (van der Wal et al., 2015; Whitehouse et al., 2018; Wan et al., 2022), but this relationship has yet to be tested empirically. Due to the difficulty in observing subglacial conditions and solid Earth dynamics beneath modern ice sheets, we turn to the deglacial sediment record of the extinct Cordilleran Ice Sheet (CIS) in the Puget Lowland. Specifically, we consider the marine-based southernmost part of the CIS, the Puget Lobe, which most recently advanced across the Puget Lowland during the Last Glacial Maximum (∼ 20 000 years ago; Mullineaux et al., 1965; Easterbrook et al., 1967; Easterbrook, 1969; Porter and Swanson, 1998). The Puget Lowland records vertical land change due to tectonics and GIA from Puget Lobe advance and retreat in the region, making it an ideal location to study the influence of solid Earth on glacial ice. Margins of Greenland comprise sedimentary basins surrounded by mountains, suggesting the topographically driven deglacial history of the Puget Lobe may be an appropriate analog to study contemporary Greenland Ice Sheet outlet glaciers (Eyles et al., 2018). Additionally, the flexural thickness of the lithosphere and mantle viscosity in Antarctica (Whitehouse et al., 2019) is similar to that of the Puget Lowland (Nield et al., 2014). Contributing to understanding the role of topography and solid Earth conditions on marine-based glacial ice can lead to development of a process-based model on the marine-terminating retreat of modern ice sheets. The findings from this work are relevant to modern glacial systems and have implications for timing of CIS contribution to the global sea level, as well as the routes and timing of human migration into the Americas (Mandryk et al., 2001; Goebel et al., 2011; Lesnek et al., 2018).

1.1 Regional context

The Puget Lowland of Washington state has been glaciated at least six times throughout the Quaternary as a result of CIS advance and retreat in the region. Glaciations occurred during marine isotope stage (MIS) 6 (∼ 97 000 to 150 000 years ago; Easterbrook, 1969), MIS 4 (80 000±20 000 years; Easterbrook et al., 1967; Easterbrook, 1969), and towards the end of MIS 2 in the Last Glacial Maximum (LGM; ∼ 17 500 cal yr BP; Mullineaux et al., 1965; Porter and Swanson, 1998). To provide context for our work, focused on the most recent glaciation of the Puget Lowland (the Fraser glaciation), we share context of previous outcrop and lacustrine work, radiocarbon constraints, and geomorphic analyses from the region.

1.1.1 Pre-Fraser glaciation and Fraser glaciation deposits

A characteristic pre-LGM deposit in the Puget Lowland is the Lawton Clay, formed as the more southern Puget Lowland became a proglacial lake basin from ice advancement into the northern Strait of Juan de Fuca (Mullineaux et al., 1965). Southward-migrating proglacial channels, active 18 000–20 000 years ago, formed extensive outwash plain deposits referred to as the Esperance Sands and mark the oncoming advance of the CIS in the Puget Lowland (Mullineaux et al., 1965; Crandell et al., 1966; Easterbrook, 1969; Clague, 1976; Booth, 1994). The final stage of ice sheet advance during late-stage MIS 2, or the Fraser glaciation, is marked by the deposition of the massive diamicton called the Vashon Till (Willis, 1897; Easterbrook, 1969; Clague, 1981; Domack, 1983; Easterbrook, 1986). Previously radiocarbon-dated wood collected beneath the Vashon Till provides a maximum age for the timing of final ice advance to the latitude of around Seattle (47.608013° N) at ∼ 14 500 ¹⁴C yr BP (∼ 17 500 cal yr BP; Mullineaux et al., 1965; Porter and Swanson, 1998), although timing of maximum ice extent near Olympia, Washington (47.037872° N), is unknown, and the degree of subglacial reworking and erosion of underlying strata is not well understood. The advance of the CIS near British Columbia and southeast Alaska is thought to have occurred after 17 000 cal yr BP (Heaton and Grady, 2003; Ward et al., 2003; Lesnek et al., 2018; Dalton et al., 2020), suggesting the Puget Lobe advanced slightly before or around the same time as the western and northwestern margins of the rest of the CIS. However, rates of terminus advance are hypothesized to be between 80 and 200 m yr⁻¹ (Booth, 1987; Easterbrook, 1992). Additionally, there are disagreements about synchronicity of final advance of the Puget Lobe and the more northern, westward-flowing Juan de Fuca Lobe (Bretz, 1920; Waitt and Thorson, 1983; Easterbrook, 1992).

Calculated from glacial erratics in the Cascade and Olympic mountains, the maximum thickness of the Puget Lobe is thought to range between 300 and 1200 m, with the glacial lobe thickening northward (Easterbrook, 1963; 1979; Thorson, 1980). Ice surface slope (60 cm km⁻¹) and rate of basal sliding (650 m yr⁻¹), calculated from proposed ice thickness suggest at its LGM extent, suggest the Puget Lobe was comparable to modern temperate glaciers (Booth, 1987; Easterbrook, 1992).

The contact between the outwash and the overlying Vashon Till suggests erosion, at least in part due to extensive glacial scour during till deposition (Bretz, 1913; Easterbrook, 1968, 1992). Till deformation due to high pore pressure beneath the ice (Booth, 1984) is also evidenced by preserved streamlined subglacial bedforms through deformation across the lowlands (Booth, 1984; Goldstein, 1994; McKenzie et al., 2023). Incised channels that served as drainage pathways in the subglacial environment (Booth, 1984) are thought to have co-evolved with the subglacial streamlined bedform field (Goldstein, 1994) during LGM ice occupancy in the Puget Lowland. These geologic and geomorphic records collectively indicate highly dynamic polythermal conditions at the ice–bed interface in association with the Vashon Till.

1.1.2 Deglacial and nonglacial Holocene deposits

Recording the transition between a subglacial and glaciomarine environment is the shell-bearing Everson Glaciomarine Drift deposits, present along southern coastal outcrops in the lowlands (Armstrong et al., 1965; Easterbrook, 1969; Powell, 1980; Thorson, 1980; Pessl et al., 1981; Domack, 1983, 1984; Dethier et al., 1995; Swanson and Caffee, 2001). The Everson Glaciomarine Drift marks intrusion of marine water into the Puget Lobe subglacial environment, primarily grounded below sea level (Thorson, 1980; Dethier et al., 1995; Demet et al., 2019). The oldest marine shells dated from the Everson Glaciomarine Drift suggest the Puget Lowland was partially deglaciated and open to marine influence by 13 470±90 and 13 090±90 ¹⁴C yr BP (∼ 16 275 and 15 750 cal yr BP, respectively; Easterbook, 1992; Dethier et al., 1995; Stuiver et al., 1998; Swanson and Caffee, 2001). The lack of sufficiently documented stratigraphic context for individual ages and a lack of marine reservoir correction (MRC) for this region contribute to uncertainties in this generalized date of deglaciation in the Puget Lowland (e.g., Porter and Swanson, 1998; Swanson and Caffee, 2001). Additionally, conflicting ages from freshwater lacustrine organics on the eastern fringe of the Puget Lowland suggest ice retreat before ∼ 13 600 ¹⁴C yr BP (∼ 16 500 cal yr BP; Rigg and Gould, 1957; Leopold et al., 1982; Anundsen et al., 1994), and cosmogenic nuclide production rates, calculated from wide-ranging radiocarbon ages marking deglaciation, indicate that retreat occurred ∼ 15 500 years ago (Swanson and Caffee, 2001). Disagreements in the timing of ice retreat remain, while emerging research suggests much of the CIS experienced Pleistocene Termination mass loss before significant climate reversals (Menounos et al., 2017; Walczak et al., 2020).

The distribution of the Everson Glaciomarine Drift has been used to suggest marine incursion, inciting a rapid lift-off of grounded Puget Lobe ice (i.e., rapid transition from grounded ice to a floating ice shelf; Thorson, 1980, 1981; Waitt and Thorson, 1983; Booth, 1987; Booth et al., 2003). Synchronous retreat of the Puget Lobe and the largely westward-flowing Juan de Fuca Lobe due to decoupling of the Puget Lobe from its bed due to marine incursion has also been suggested (Easterbrook, 1992). However, major differences in deglacial stratigraphy across the Puget Lowland (Powell, 1980; Pessl et al., 1981; Domack, 1984; Demet et al., 2019) indicate variable patterns of ice-marginal retreat in time and space. The presence of grounding zone wedges (GZWs) across the region, sedimentary deposits found at the margin of marine-terminating ice lobes, also supports the idea of stepwise retreat with periodic stability (Simkins et al., 2017; Demet et al., 2019). The development of these ice-marginal landforms were likely supported by the identified high rates of sedimentation in the region (∼ 2.5 mm yr⁻¹; Simkins et al., 2017, 2018; Demet et al., 2019).

Additionally, the modern elevation of marine limits in the Puget Lowland range from ∼ 125 m above sea level in the northern San Juan Islands to lower than 30 m at the southern end of Whidbey Island (Thorson, 1981, 1989; Dethier et al., 1995; Kovanen and Slaymaker, 2004; Polenz et al., 2005), indicating highly variable rates of GIA across the region. Ice retreat, as marked by the deposition of the Everson Glaciomarine Drift, is thought to have occurred at the same time as rapid isostatic uplift between 13 600 and 11 300 ¹⁴C yr BP (Dethier et al., 1995). The rate of landscape emergence due to GIA in the Puget Lowland may have been as high as 10 cm yr⁻¹ during early deglaciation (Dethier et al., 1995). The high rate of GIA-induced uplift during glacial retreat suggests that relative sea level fall in the Puget Lowland outpaced rapid global sea level rise, leading to the emergence of the landscape during the end of the LGM (Shugar et al., 2014; Yokoyama and Purcell, 2021). Emergence of this landscape from below to above sea level is also distinctly marked in post-glacial stratigraphy by thin subaerial deposits (e.g., fluvial sediments and soil) overlying the glacial and glaciomarine deposits (Domack, 1984; Demet et al., 2019). Both pre-existing topography and GIA uplift could have periodically stabilized the Puget Lobe during retreat, as suggested for contemporary ice sheets (Durand et al., 2011; Favier et al., 2016; Alley et al., 2021; Robel et al., 2022), highlighting the importance of elucidating the influence of both variables on ice sheet behavior.

Despite the large amount of research conducted on the Puget Lobe, many of the radiocarbon ages from this region lack detailed stratigraphic context for their proposed constraints (e.g., Easterbrook, 1992, and references therein). Additionally, the continued absence of a local MRC has left some uncertainties in the marine-based radiocarbon used in the region. Pleistocene variability in marine reservoir ages in this region (Schmuck et al., 2021) make creating a widely applicable MRC difficult, but developing an understanding of recent marine reservoir effects local to the Puget Lowland would increase reliability of marine ages. Recent advancements in technologies for radiocarbon dating, sedimentology, and stratigraphic analyses also beg a reanalysis of Puget Lowland deposits that have not been observed using modern understandings of glaciology. Subsequently, the need to clarify spatiotemporal details of ice retreat patterns and drivers of Puget Lobe retreat persists.

1.2 Relevance to solid Earth dynamics and modern ice sheets and glaciers

Based on modeled evidence of GIA control on ice behavior in analogous Antarctic Peninsula glacial catchments (Nield et al., 2014; Whitehouse et al., 2019), in addition to previously identified geomorphic evidence of ice margin standstill in the Puget Lowland (Simkins et al., 2017; Demet et al., 2019), we hypothesize that the landscape position above and below sea level, in response to loading and unloading of the solid Earth, influenced ice margin positions and drove a punctuated retreat of the CIS during the late Pleistocene. In the central Puget Lowland, Whidbey Island spans nearly 100 km in distance along the north–south direction of glacial ice movement and hosts extensive coastal bluff features (Fig. 1). The outcrops, composed of glacial and interglacial sediments, preserve details of ice advance and retreat across the formerly marine landscape, as well as landscape transitions that took place coeval with deglaciation. Except for localized tectonic deformation of surficial sediments (Sherrod et al., 2008), local LGM and subsequent deglacial deposits appear to have little post-depositional reworking (Booth and Hallet, 1993; Kovanen and Slaymaker, 2004; Eyles et al., 2018; Demet et al., 2019; McKenzie et al., 2023).

Figure 1 (a) Regional map showcasing maximum extent of the Cordilleran Ice Sheet (CIS) at 15, 14, and 13 kyr cal yr BP (data synthesized by Henry Haro from Ehlers et al., 2010). This map has been adapted from Earthstar Geographics satellite imagery. Glaciomarine shell radiocarbon dates from the literature (red dots) were recalibrated using Marine20, and new marine reservoir corrections are listed in calendar years before present (cal yr BP). Only glaciomarine shell ages with available metadata and error standards relevant to work presented here were included in the analysis (e.g., Easterbrook, 1992; Dethier et al., 1995; Swanson and Caffee, 2001). Information about ages and recalibration conducted in this work may be found in Table 1. (b) Whidbey Island inset map with sites labeled south to north. Grounding zone wedges (GZWs) identified and inferred from Demet et al. (2019) are mapped in blue. Streaming of the bed visible around the margins of the inset and on the southwest side of Whidbey Island (outlined in gray). This map has been adapted from data provided by the © Esri, Garmin, NaturalVue, Airbus, USGS, National Geodetic Survey (NGS), NASA, Consortium of International Agricultural Research Centers (CGIAR), National Center for Ecological Analysis and Synthesis (NCEAS), NLS, Ordnance Survey (OS), National Mapping Agency (NMA), Geostastyrelsen, GSA, GSI, and the GIS user interface.

In this work, decimeter-scale stratigraphic and sedimentological assessments are complemented by accelerator mass spectrometry radiocarbon (¹⁴C) and optically stimulated luminescence (OSL) dating. While these two dating methods have been utilized in this region for decades (e.g., Rigg and Gould, 1957; Leopold et al., 1982; Easterbrook, 1992; Anundsen et al., 1994; Dethier et al., 1995; Swanson and Caffee, 2001), our hypothesis of the relationship and timing of landscape emergence in relation to ice retreat and periodic stabilization of ice retreat has not been directly assessed. Additionally, with this work, we develop a modern and regional MRC that is applied to the dates for marine deglaciation. Therefore, the application of advances in geochronology paired with a high-resolution stratigraphic assessment of Whidbey Island is a novel approach to elucidating the ice retreat and land emergence across the region.

2 Materials and methods

2.1 Sedimentology and stratigraphy

Sediment samples were collected from the following Whidbey Island outcrops: (a) Double Bluff, (b) Fort Casey, (c) Penn Cove, (d) West Beach, and (e) Cliffside at 10 cm intervals (Fig. 1b; Table S1 in the Supplement), with additional subsamples collected from units with laminations, lenses, or rip-up clasts. Sites were determined based on accessibility from the beachfront. If some outcrop units were pinched out while others continued, these multiple-location data collection areas were considered one site (indicated by white dots on the stratigraphic columns in Fig. 2). However, if continuity between units was not visible or accessible across a single beachfront, these sites were separated into “1” and “2”, such as in the case of Fort Casey and West Beach. Thin (∼ < 0.5 cm thick) horizontally continuous layers are referred to as laminations, while less continuous layers that pinch out are referred to as lenses (e.g., Fig. S1). Over 300 discrete bulk sediment samples were analyzed at the University of Virginia for grain size and magnetic susceptibility (MS). An additional 15 peat, wood, and marine shell samples were excavated for radiocarbon dating. Grain size analyses were conducted via a Bettersizer S3 Plus particle analyzer on sample matrix material (material ≤ 3 mm in diameter), and MS measurements were collected with a Bartington MS2 magnetic susceptibility meter. MS values provide information about the number and size of magnetic grains in each sample, elucidating the continuity and source of biogenic and lithogenic deposits (Thompson and Oldfield, 1986; Verosub and Roberts, 1995; Rosenbaum, 2005; Hatfield et al., 2017; Reilly et al., 2019). Grain sizes from the Bettersizer S3 Plus particle analyzer are categorized as silt or clay (0.06–63 µm), very fine sand (63–125 µm), or fine sand (125–250 µm) (Wentworth, 1922; Folk and Ward, 1957). Results of the Whidbey Island stratigraphy are presented according to latitudinal location, starting with the southernmost site, Double Bluff, followed by the Fort Casey sites, Penn Cove, West Beach sites, and ending with the northernmost Cliffside and Rocky Point sites.

Figure 2 Outcrop sites from south to north: Double Bluff, Fort Casey 1, Fort Casey 2, Penn Cove, West Beach Site 1, West Beach Site 2, and Cliffside, represented by stratigraphic column with collected radiocarbon and OSL and grain size data. Icons indicate where shells or wood were present in the stratigraphy. Not all occurrences of wood or shells were radiocarbon dated. The white dots on the stratigraphic columns indicate the end of one visible region and start of a new location where visible units were mapped. Colors alongside stratigraphic units indicate grain size graph correlations and are not correlated between sites – each site is independently considered. Background colors on the grain size graphs indicate transitions in grain size. From left to right, gray is clay/silt, light brown is very fine sand, and dark brown is fine sand. Variations in sampling resolution are a function of accessibility to outcrops from the beachfront. Some units were more accessible for sampling than others.

2.2 Accelerator mass spectrometry radiocarbon analysis

Assuming a constant cosmically produced ¹⁴C to ¹²C ratio, the variation in this ratio can be used to determine the amount of time since the death of formerly living specimens. Samples were run at the National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) Laboratory at Woods Hole Oceanographic Institute. The unprocessed wood material underwent a series of six to eight acid–base–acid leaches to remove contamination and inorganic carbon prior to combustion. The carbonate shell samples underwent carbonate hydrolysis, and the resulting carbon combustion reacted with an Fe catalyst along vacuum-sealed lines to produce graphite (NOSAMS, 2023). Resulting graphite pellets were pressed into targets and analyzed by accelerator mass spectrometry in addition to standard and processing blanks (Roberts et al., 2019). The AMS measurements determined the ratio of ¹⁴C to ¹²C in each of the pellets, which was then used to calculate the radiocarbon age using the Libby ¹⁴C half-life of 5568 years (Stuiver and Polach, 1977; Stuiver, 1980).

Conversion of radiocarbon years to calendar years before present was conducted using the Int20 curve for terrestrial carbon samples and the Marine20 curve for marine shell samples using the Calib 8.2 interface (Heaton et al., 2020). Marine20 is the baseline marine curve used for Calib 8.2 and is the most up-to-date, internationally agreed marine radiocarbon age calibration curve for non-polar global average marine records (Heaton et al., 2020). A modern marine reservoir correction was calculated in Calib 8.2 and applied to all carbonate shell samples, both newly presented in this work and to previously published radiocarbon ages, using contemporary shells with known pre-1955 (i.e., prior to nuclear bomb testing) collected dates from the Burke Museum in Seattle, Washington.

The modern (pre-1955) shells from the Burke Institute were live-collected between 1911 and 1931. Species of the modern bivalves include Modiolus rectus, Musculus niger, Cardita ventricosa, Macoma carlottensis, Mya arenaria, and Macoma nasuta. The radiocarbon ages calculated from these specimens range from 815±15 to 925±20 ¹⁴C yr (Table S2). Utilizing the marine reservoir correction curve developed by Calib 8.2, an average marine reservoir correction for this region is 264 ¹⁴C yr, which was applied to other marine shells from this work and previously published work (Fig. 1). While there is a narrow range of marine reservoir effects between 211 and 318 ¹⁴C yr, a species-specific effect was not observed. It is also important to note that due to known Pleistocene variability in marine-reservoir ages in this region (Schmuck et al., 2021), this marine reservoir correction is most accurate for marine organisms aged between 1911 and 1931. However, while these limitations exist, this new reservoir age is highly localized to the Puget Lowland and was therefore extended for use in older collected marine shells.

2.3 Optically stimulated luminescence

In depositional environments, minerals are exposed to radiation from in situ uranium (Ur), thorium (Th), potassium (K), and cosmic rays (Rhodes, 2011; Duller, 2015). Incoming radiation excites electrons which are trapped in structure deformities of quartz and feldspar grains (Rhodes, 2011). When exposed to sunlight, electrons are released from the traps. In returning to their original states, they emit luminescence and the mineral is reset. Upon burial, trapped electrons re-accumulate, and the number is proportional to the burial time and the radiation exposure, termed “dose”. The rate of irradiation, the “dose rate”, can be calculated from the cosmic flux as well as the U, Th, and ⁴⁰K concentrations of the surrounding sediments. The OSL signal is proportional to the dose and can be measured by exposing the mineral to light in a controlled setting. An age since burial can be determined by dividing the dose by the dose rate (Tables S3, S4, and S5).

Materials from glacial environments present challenges due to the potential of the OSL signal not being fully reset between transport and deposition (Wallinga and Cunningham, 2015). Additionally, extensive overburden pressure from glacial ice has the potential to partially or completely reset OSL signatures, which could provide large error to the final OSL stage (Bateman et al., 2012; King et al., 2014). Subglacial environments, especially those under ice streams, have a presence of significant meltwater which can saturate sediment pore space and influence quartz and feldspar exposure to radiation at the time of and for an extended period of time after deposition (Wallinga and Cunningham, 2015; Duller, 2013).

While a detailed description of the OSL procedure can be found in the Supplement (Sect. S1), a summary is provided here. In order to avoid premature bleaching of samples, they were collected before sunrise or after sunset, only exposed to low-energy red light, and wrapped in dark black plastic before being transported to East Carolina University (ECU) for preparation and processing. Samples were prepared for OSL analysis under darkroom conditions, using standard procedures to extract 63–212 µm quartz. Due to feldspar contamination, a post-IR (infrared) blue single-aliquot regenerative dose (SAR) procedure was used to measure the quartz equivalent dose (Murray and Wintle, 2000; Wallinga et al., 2002; Wintle and Murray, 2006).

Bulk sediment was collected from outcrops for high-resolution gamma spectrometry measurements and stored for at least 4 weeks prior to measurement. OSL samples were taken at unit boundaries, while dose rate samples were only taken from the same unit as the OSL samples. Therefore, the gamma dose rates reflect the sample unit only and contain no information about adjacent, underlying, or overlying units. Uranium concentrations determined from ²³⁴Th were all significantly higher than concentrations determined from ²¹⁴Pb and ²¹⁴Bi. We assumed that ²³⁴U was leached out of the sample due to in situ water presence.

The sample ages, calculated in calendar years, were calculated by dividing the dose by the dose rate (Tables S3, S4, S5). For samples with feldspar contamination that showed fading, the ages were corrected as suggested by Auclair et al. (2003). While ¹⁴C ages are reported in thousands of calendar years before present (cal BP; 1955), all OSL ages are reported in kyr based on the date of collection (2020). OSL ages in kyr can be directly compared to kyr cal BP by subtracting 72 years from the OSL age. Cases of large overdispersion (> 20 %) indicate the mixing of grains or non-uniform resetting of the samples. Possible reasons could be from bioturbation, mixing with light-exposed surface grains during collection, or incomplete resetting during deposition. In Sect. 3, we cannot attribute much relevance to overdispersed samples.

3 Results

We will be moving through results from the southernmost to the northernmost site. Numerical schemes to describe units at each site are independent and do not correlate between sites (Fig. 2). Stratigraphic columns were developed to represent our interpretation of physical data present at several locations across these sites (Fig. 2).

3.1 Double Bluff

The stratigraphically lowermost unit visible at Double Bluff, Unit 4, is a visually well-sorted sand with sparse rounded gravel lenses. Unit 4 is normally graded with clasts ranging from granule to pebbles, with a consistently horizonal long-axis orientation and occasional silt rip-ups from nonvisible underlying units. A gradational boundary leads into the overlying sandy silt and fine clayey silt of Unit 3. This unit contains wavy laminations and woody debris dated to be 46.7+ kyr cal BP (i.e., “radiocarbon dead”; Table 1, NOSAMS receipt no. 171378). Unit 3 generally fines upwards but with variable matrix grain size modes from 10–500 µm (Fig. 2). Unit 2 is composed of massive diamicton with a clay and fine-silt matrix marked by a matrix grain size mode of 8 µm and a mix of angular and rounded granule to cobble-sized clasts without a preferred long-axis orientation. There is a gradational contact between Unit 2 and Unit 1. Unit 1 consists of diamicton with a matrix varying between sandy silt and silty sand with woody debris dated to 48.0+ kyr cal BP in age (i.e., radiocarbon dead; Table 1, NOSAMS receipt no. 176245) and clasts that are predominantly aligned parallel to bedding and evidence of soft-sediment deformation. This uppermost unit has interbedded silt and clay, as well as marine shells in the upper 50 cm of silt, that were inaccessible for sampling. MS values in Unit 3 are distinctly lower than the other units (Fig. S2).

Table 1 Radiocarbon sample descriptions and data. Gray rows indicate previously published radiocarbon data that have been recalibrated using Marine20 and our MRC. Note that n/a stands for not applicable. The abbreviations s.h. and a.s. stand for “shell hash” and “articulated shell”, respectively.

3.2 Fort Casey

The lowermost visible unit, Unit 3, at Fort Casey Site 1 consists of massive diamicton with a fine-silt and clay matrix and randomly oriented pebble to cobble-sized angular and rounded clasts. Interbedded with the massive diamicton are discrete gravel and sand laminations at the base of Unit 3 and silt and clay laminations with rip-ups and woody debris toward the top of Unit 3. Unit 2 consists of fine-sand to pebble-sized clasts in a sandy silt matrix with vertically oriented and reverse-graded angular clasts. Unit 2 has a remarkably consistent matrix grain size throughout the unit. An OSL sample from this unit (Table 2, sample no. 1) could not be dated reliably due to extremely low signals. This unit also contains sand and silt lenses with mud and plant rip-ups (Fig. 2). A gradational boundary leads to Unit 1, which is massive diamicton similar to Unit 3 but with a matrix distinctly lighter in color.

At Fort Casey Site 2, the lower visible unit, Unit 5, contains interbedded clay and sand with reverse grading (Fig. 2). Unit 4, in which no samples were collected, consists of diamicton with concentrated granule to pebble lenses and clay and silt lenses, as well as evidence of soft-sediment deformation. Unit 3 is a massive clay, followed by the Unit 2 layer of silt about 20 cm thick, continuous across an irregular, undulating, and most likely erosional contact. OSL dates at the top of Unit 2 and base of Unit 1 were found to be 40.8±8.2 and 56.6±15.5 kyr (Table 2, sample nos. 3 and 2). The overlying Unit 1 is a diamicton with very fine sand to cobble-sized angular and rounded clasts. Normal grading is present in the matrix of Unit 1 with fractured (i.e., seemingly crushed) granite clasts.

Table 2 OSL age data with overdispersion percentages and total dose rate values. Final sample ages are bolded. To directly compare OSL and ¹⁴C ages, it would be necessary to subtract 72 years from the OSL ages. This correction is considerably smaller than the uncertainty in the ages and can therefore be neglected. Additional dose and dose rate data may be found in Tables S4 and S5.

3.3 Penn Cove

The lowest visible unit at this site, Unit 5, comprises a reverse-graded diamicton with a coarsening upward sand matrix and rounded granules and pebbles (Fig. 2). Following a sharp boundary with Unit 5, Unit 4 consists of silt and sand laminations with cross-bedded sands near the top. Unit 4 deposits were OSL-dated to ages 56.6±4.1 and 44.4±2.8 kyr (Table 2, sample nos. 10 and 11), where the sample with older age showed a slightly large overdispersion (> 20 %) which we attribute to incomplete resetting. The grain size modes for Unit 4 matrix are predominantly between 500–700 µm (Fig. 2). An erosional boundary at the top of Unit 4 leads to the massive clayey silt diamicton of Unit 3 with rounded fine- to cobble-sized clasts and occasional sandy silt and silt lenses. A gradational boundary separates Units 3 and 2, which is a massive clay diamicton with rounded fine sand to cobble grains and articulated shells. Six shells from Unit 2 were radiocarbon dated, with ages spanning 14.8±0.3 to 14.1±0.3 kyr cal BP (Table 1, NOSAMS receipt nos. 176239-176242, 171380, and 171381). Unit 2 also contains sand lenses and wood fragments. Unit 2 has a sharp contact with Unit 1, which consists of normally graded gravel with rounded and angular small to large pebbles with no predominant long-axis orientation. A mode of clay-sized grains is visible in Units 2 and 3 but is not visible in Unit 1 (Fig. 2).

3.4 West Beach

At West Beach Site 1, the lowest unit, Unit 5, consists of matrix-supported diamicton with randomly oriented clasts and two matrix grain size modes at 8 and 20 µm (Fig. 2). This unit has a sandy-silt lamination that interrupts the diamicton. The diamicton above the silty-sand lamination, however, contains highly irregular dips and soft-sediment deformation. Unit 5 has a gradational boundary with Unit 4 – a light clay layer deposited on a laterally irregular surface marked by normal grading or fining upward (Fig. 2). Unit 3 consists of a thick, 0.25 m clast-supported gravel layer with poorly sorted fine-sand to cobble-sized clasts oriented parallel to the depositional bed. A sharp, horizontally regular contact occurs between Unit 3 to the 0.75 m thick, well-sorted sand of Unit 2, with OSL ages of 6.2±0.6 and 4.1±1.8 kyr (Table 2, sample nos. 5 and 4); however, these ages experienced large overdispersion (> 20 %). Unit 2 has a gradational contact with Unit 1, which is a modern soil on top of a basal shell hash dating between 0.62±0.1 and 0.39±0.1 kyr cal BP (Table 1, NOSAMS receipt nos. 176236–176238). MS values are similar throughout Units 5, 4, 2, and 1 but decrease in Unit 3 (Fig. S2).

At the base of West Beach Site 2 are cross-bedded and coarse-sand laminations. OSL dates from the lowermost sand in Unit 8 are dated to 31.3±2.7 and 38.1±9.7 kyr (Table 2, sample nos. 7 and 6), with overlying Unit 7 sediments OSL dated between 30.7±2.5 and 29.2±4.6 kyr (Table 2, sample nos. 8 and 9); however, these ages experienced large overdispersion (> 20 %). A gradational contact leads into Unit 7, consisting of silt and clay with radiocarbon dead woody debris (48.0+ kyr cal BP; NOSAMS receipt no. 176243). Unit 6 consists of sand with wavy bedding and silt laminations. No samples were collected from Units 5 and 4, consisting of a peat layer and a unit of sand and silt laminations, respectively. The Unit 3 diamicton matrix coarsens upwards, and this unit has many grain size modes between 5 and 70 µm (Fig. 2). Unit 2 consists of diamicton with a fine sand matrix and clasts as large as pebbles and is not spatially continuous throughout the site. A gradational boundary leads into the 0.5 m thick layer of Unit 1, consisting of predominantly of silt.

3.5 Rocky Point, Cliffside

The lowest visible unit at Cliffside, Unit 6, consists of fine-sand to cobble-sized rounded clasts. This massive diamicton has no preferential orientation for clast long axes. The matrix changes from clay to sand and includes sediment deformation beneath clasts (Fig. 2). Unit 6 gradationally transitions to Unit 5, which is a normally graded fine-sand to cobble-sized clast diamicton. Unit 5 is normally graded gravel lenses containing clasts with consistent horizontal long-axis orientation. Unit 5 gradually transitions into the granule and sand layer of Unit 4, which includes sand and silt lenses within gravel-rich and wavy laminations. Unit 3 intrudes into Unit 4 and consists of a massive diamicton with rounded, cobble-sized clasts. The matrix of Unit 3 has two grain size modes at 5 and 20 µm (Fig. 2). Two of the lower-unit samples for Cliffside Unit 3 were taken from the more southern Rocky Point site, as the identified Unit 3 is continuous throughout both sites. Unit 3 gradually transitions into Unit 2, which is a laterally discontinuous light-clay unit with silt layers. Unit 1 is comprised of mostly rounded, normally graded crushed material with fine to large cobble-sized clasts.

4 Interpretation and discussion

We use the sedimentological units described in Sect. 3 to establish a facies model that encompasses glaciomarine and coastal sedimentary processes and depositional environments (i.e., emergent or submergent landscape). Aided by geochronological constraints, this facies model is applied to the stratigraphic sequences observed at each site to construct a regional history of ice behavior and landscape evolution before, during, and following the LGM (Fig. 3).

Figure 3 (a) Grouping of facies based on depositional time periods across Whidbey Island. Units with asterisks have radiocarbon or OSL dates associated with them and are shown to the right of the unit numbers. (b) Top schematic drawing indicates Puget Lobe ice retreat within a marine environment at time 1. The bottom schematic showcases hypothesized northernmost retreat into a subaerial environment at time 2 after landscape reemergence from time 1. Puget Lobe ice retreat in a marine environment is only seen to have occurred at southernmost sites of Double Bluff and Penn Cove, while ice retreat into a subaerial environment is proposed for the West Beach and Cliffside sites.

4.1 Facies interpretation

Structureless diamicton with randomly oriented clasts of variable size, roundness, lithology, and a range in matrix size are classified as till or sediments deposited directly by glaciers in the subglacial environment (Boulton and Deynoux, 1981; Evans et al., 2006; Sengupta, 2017). Some biological material may be incorporated into till in the form of broken shells or woody fragments. This reworked biogenic material may be incorporated into the ice as it moves across the landscape; therefore, radiocarbon ages of biogenic material will be older than glacial occupation. These characteristics are consistent with glaciomarine tills described offshore of West Antarctica (e.g., Kirshner et al., 2012; Prothro et al., 2018; Smith et al., 2019) and western Greenland (Sheldon et al., 2016; O'Regan et al., 2021), as well as tills deposited by the relict British–Irish Ice Sheet (Evans and Thompson, 2010). Lower boundaries of till units are often characterized by erosional contacts, reflecting glacial advance and erosion of pre-existing substrate, and may contain rip-up clasts from underlying units. Due to similarities in facies formerly identified as tills and some presence of reworked woody fragments, units classified as (local) LGM till (i.e., Vashon Till) in the Puget Lowland include Unit 2 from Double Bluff, Unit 3 at Fort Casey Site 1, Unit 1 from Fort Casey Site 2, Unit 3 from Penn Cove, Unit 5 from West Beach Site 1, and Unit 6 from Cliffside (Figs. 2, 3a). Little post-depositional erosion or reworking of this glacial material is consistent with previous work identifying tills in the region (Booth and Hallet, 1993; Kovanen and Slaymaker, 2004; Eyles et al., 2018; Demet et al., 2019).

Glacial outwash is characterized as diamicton with a range of well-rounded and some angular clasts with a parallel to bedding clast orientation that suggests sediment transport via proglacial meltwater from an upstream source of glacial ice (Boulton and Deynoux, 1981). This facies may indicate deposition in a subaerial or subaqueous environment, but importantly, clast orientation distinguishes proglacial outwash from till (Hagg, 2022). The deposits may also exhibit normal grading and/or sedimentary structures indicative of soft-sediment deformation (e.g., loading structures, flame structures, and sediment deformation beneath clasts; Boulton and Deynoux, 1981). Glacial outwash recorded in British Columbia (Clague, 1975) and the forefield of Mýrdalsjökull ice cap in Iceland (Kjær et al., 2004) feature similar structures seen in several units among our Puget Lowland outcrop sites. Using the defined classification of glacial outwash, Units 1 and 2 from Fort Casey Site 1, Units 4 and 5 from Fort Casey Site 2, Units 1 and 5 from Penn Cove, and Units 1 and 3 from Cliffside are interpreted as glacial outwash deposits (Figs. 2, 3a).

A third diamicton, structurally similar to those interpreted as till yet containing articulated and/or broken marine shells, occasional absence of fines from winnowing of fine-matrix material, and sedimentary structures such as wavy laminations, is interpreted as glaciomarine deposits and composed of both glacial and pelagic sediments that accumulate on the ocean floor seaward of the ice margin. The winnowing of fine-grained material may be a product of tidal currents in a submarine or coastal setting (Smith et al., 2019). Such pelagic sediments have been sampled from a geographically diverse population of sediment cores from deglaciated continental margins (e.g., Anderson et al., 1980; Prothro et al., 2018; Smith et al., 2019), although preservation of shells and other carbonate-based materials is less common in Antarctic glaciomarine sediments. Glaciomarine deposits are also identified in coastal outcrop deposits of northern Svalbard with similar characteristics (Alexanderson et al., 2018). Both Unit 1 from Double Bluff and Unit 2 from Penn Cove are consistent with these classifications and closely resemble the structure and composition of the glaciomarine deposits identified on deglaciated continental margins (Figs. 2, 3a; Anderson et al., 1980; Prothro et al., 2018). At Double Bluff and Penn Cove sites, this facies (also known as Everson Glaciomarine Drift) overlays till, indicating ice marginal retreat into a marine setting with sand-rich deposits recording the removal of fines by bottom currents. Conversely, till that stratigraphically transitions upsection into cross-bedded sands or gravels with parallel-to-bed-oriented clasts, occasional wavy laminations and are barren of marine shells indicate retreat into a subaerial environment, as is observed proximal to the Mýrdalsjökull ice cap in Iceland (Kjær et al., 2004). Unit 3 from West Beach Site 1 and Unit 5 from Cliffside record such evidence of subaerial glacial retreat and both meet these classifications (Figs. 2, 3a).

Facies transitions where grain sizes coarsen upward (i.e., reverse grading) and changes in MS values can be associated with marine to coastal environment transition (e.g., Komar, 1977; McCabe, 1986; Sengupta, 2017) as a result of landscape emergence. Regardless of the process(es) explaining the observed grain coarsening, which may include tectonic activity, glacial isostatic response, or a combination of these factors, we would expect such processes to be marked by facies transitions along the coast. In the Puget Lowland, emergence above sea level has been recorded in the stratigraphy by thin subaerial deposits (e.g., fluvial sediments and soil) overlying the glacial and glaciomarine deposits (Domack, 1984; Demet et al., 2019). Tectonic activity as a cause for coarsening-upward trends in the Puget Lowland stratigraphy can be rejected because the till and sedimentary structures, like cross-bedding, have been preserved in the record. Coarsening-upward grain sizes seen in the transition from finer marine sediments to coastal deposits have been identified in coastal outcrops in Northern Ireland (McCabe, 1986) and Svalbard (Alexanderson et al., 2018) and are interpreted to indicate relative sea level fall. While glacial isostatic rebound is not responsible for the shallowing-upward of Svalbard facies (Alexanderson et al., 2018), the facies and coarsening material identified between Units 3 and 2 at Fort Casey Site 2 transition from Unit 5 laminated silt to Unit 4 cross-bedded sand at Penn Cove, and the coarsening of the grain size with peaks and MS across Units 7 and 6 at West Beach Site 2 could be connected to land emergence events (Figs. 2, 3a).

Facies transitions where grain sizes fine upward (i.e., normal grading) correspond with increases or decreases in MS and are accompanied by the appearance of marine shells that reworked wood or terrestrial carbon are associated with landscape submergence (Komar, 1977; Sengupta, 2017). A specific example from the sedimentological record that marks the transition from a subaerial to a submarine environment is from seismic profiles and regional stratigraphic data in the southwestern Pacific in the South Island, New Zealand (Carter et al., 1986). The fining of material between Unit 4 sand deposits to Unit 3 silts, both with the reworked radiocarbon dead wood at Double Bluff, introduction of shells to the fining material between Units 2 and 1 at West Beach Site 1, and fining of grain size across the Unit 2 and 1 boundary at West Beach Site 2 are all interpreted as a transition to a submarine setting (Figs. 2, 3a).

4.2 Pre-LGM landscape evolution

Prior to Puget Lobe glacial advance across Whidbey Island during the LGM, several submergence and emergence facies transitions record dynamic landscape changes. Landscape emergence above sea level prior to LGM glaciation is recorded by outcrops exposed at Penn Cove and Fort Casey Site 2 (Figs. 2, 3a). Penn Cove OSL ages identify this landscape emergence to occur between 56.6±4.1 and 44.4±2.8 kyr. Similar Fort Casey Site 2 OSL ages constrain this transition to having occurred from 56.6±15.5 to 40.8±8.2 kyr, placing the emergence within the MIS 4 glacial and MIS 3 interglacial stages, which may be connected to a lack of ice coverage and reduced CIS loading of the solid Earth at these times.

A sequence of submergent and emergent facies are also observed in the pre-LGM deposits at West Beach Site 2. OSL dates place a submergence event between 38.1±9.7 and 31.3±2.65 kyr, while OSL dates from overlying facies places subsequent emergence between 30.7±2.5 and 29.2±4.6 kyr (Fig. 3a). Both of these events occurred within the MIS 3 interglacial. This rapid transition between landscape submergence and emergence not only identifies high sedimentation rates at this site but also suggests that the Puget Lowland experienced rapid landscape changes during MIS 3. Clay and sand deposits included as part of the emergence and submergence interpretation may have previously been identified and referred to as the well-known Lawton Clay (Mullineaux et al., 1965) and Esperance sands (Mullineaux et al., 1965; Crandell et al., 1966; Easterbrook, 1969; Clague, 1976; Booth, 1994). The pre-LGM timing of proglacial lake basin development, responsible for the deposition of the Lawton Clay (Mullineaux et al., 1965) and subsequent channel deposition of the Esperance sands, has consistently been radiocarbon dated to 18 000–20 000 years ago (Mullineaux et al., 1965; Crandell et al., 1966; Easterbrook, 1969; Clague, 1976; Booth, 1994). It is highly likely that the uncertainties from our OSL dates would contribute to discrepancy with the well-established radiocarbon dates of the Esperance sands (Sect. S1; Easterbrook, 1969). Despite these uncertainties, our OSL ages provide enhanced understanding of sediment deposition and landscape evolution rates.

4.3 LGM glacial advance

Erosional contacts between till (Vashon Till) and underlying facies mark LGM advance of the Puget Lobe into the region. These erosional contacts beneath till are seen at multiple sites across Whidbey Island, including Double Bluff, Fort Casey Site 2, and Penn Cove (Figs. 2, 3b). OSL ages we collect from below the erosional contact of LGM tills (56.6±4.1 and 44.4±2.8 kyr, within the time frame of MIS 3) cannot be used to find precise maximum age of ice extent. However, previously radiocarbon-dated wood material place final LGM advance into the region after 17 500 cal yr BP (Mullineaux et al., 1965; Porter and Swanson, 1998). This major difference in ages suggests a great deal of erosion at the boundary between underlying sediments and till deposition of the Puget Lobe during ice advance. The erosional surface of outwash deposits underlying LGM till, in some areas up to 75 m (Easterbrook, 1992), is thought in part to be due to glacial scour during ice advance (Bretz, 1913; Easterbrook, 1968, 1992).

4.4 Deglaciation

Glaciomarine sediments (Everson Glaciomarine Drift) in the uppermost 50 cm of Double Bluff Unit 1 record retreat of the Puget Lobe within a marine environment (Fig. 3b; Thorson, 1980; Dethier et al., 1995; Demet et al., 2019). At Penn Cove, the presence of articulated shells in growth position (Fig. 2) within a unit increasingly lacking small grain sizes upcore suggests ice retreat in a marine environment with possible tidal influence (Smith et al., 2019). With the newly developed regional MRC of 264 ¹⁴C yr, the six articulated shells found at Penn Cove were radiocarbon dated to a range of dates between 14.8±0.3 and 14.1±0.3 kyr cal yr BP (Table 1). These ages strongly agree with previous literature marking marine incursion beneath the ice (Fig. 1). Based on the 2σ error in glaciomarine radiocarbon dates (both those presented here and those recalculated from the literature using Marine20 and our MRC; Table 1), glacial ice appears to have been stable at Penn Cove for at least 500 years with concordant high sedimentation rates, accumulating 2.5 m during a near-700-year period. The oldest known presence of Everson Glaciomarine Drift was presented in prior work, dated to ∼ 16 275 and 15 750 cal yr BP (Easterbook, 1992; Dethier et al., 1995; Stuiver et al., 1998; Swanson and Caffee, 2001). However, the specific location, stratigraphic context, and proximity between shells within sedimentological data are missing from these reports. The shells presented in this work were collected from within two units spanning 3 m vertically within the outcrop (Fig. 2) and may represent a younger limit of Everson Glaciomarine Drift deposition before final ice retreat.

Deglacial facies seen at the more northern West Beach Site 1 and Cliffside indicate ice retreat within a subaerial environment (Fig. 3b). The change in ice retreat style seen from the more southern Double Bluff and Penn Cove sites to the northern West Beach and Cliffside sites may be due to the substantial standstill of ice at Penn Cove. The presence of a GZW at Penn Cove (Fig. 1; Simkins et al., 2018; Demet et al., 2019) further supports the idea that the marine-terminating ice was pinned at this location for a substantial amount of time. Additionally, the Rocky Point site features a bedrock high (i.e., a potential pinning point of ice; Hogan et al., 2020), and other mapped GZWs (Fig. 1) suggest several points across Whidbey Island could have periodically stabilized ice due to land rebound during the final ice retreat (Simkins et al., 2018; Demet et al., 2019). The paired stratigraphic-, geomorphic-, and geochronological-based evidence of period ice stability on Whidbey Island provides empirical evidence for sedimentation and GIA as possible mechanisms for ice stabilization during retreat. Within this work, we also identify that the Puget Lobe experienced stepwise retreat, rather than catastrophic loss of ice due to rapid unpinning (see Easterbrook, 1992).

4.5 Nonglacial Holocene landscape evolution

In the sediment record following final glacial influence, the Penn Cove and Cliffside sites contain outwash deposits from proglacial fluvial sources. An OSL age within the submergence facies of Unit 2 at West Beach Site 1 marks the subsequent transition from a post-glacial fluvial environment to a submarine environment between 6.2±0.6 and 4.1±1.8 kyr (Figs. 2, 3a). Radiocarbon-dated shell hash sampled from the uppermost unit at this same West Beach Site 1 suggests a highly energetic aquatic marine or coastal environment was present in this location as early as 0.62±0.1 kyr cal BP through at least 0.39±0.1 kyr cal BP (Figs. 2, 3a). After a slow in initial lithospheric rebound from ice loading or a possible local tectonic event, it is feasible that vertical land movement slowed enough to allow local sea level to resubmerge the region between 600 and 390 years ago (Figs. 2, 3a). In our analysis of nonglacial Holocene sediments in the Puget Lowland, we find that this landscape still evolved rapidly due to ongoing effects of GIA and tectonics in the region following deglaciation.

5 Conclusions

This decimeter-scale physical sedimentological assessment, paired with geochronological assessment of seven sites across the deglaciated Puget Lowland, provides spatiotemporal information on landscape emergence and submergence, as well as final ice advance and retreat of the southernmost CIS. Rates of vertical landscape changes constrained through OSL dating indicate the Puget Lowland was a highly dynamic region where a sequence of landscape emergence and submergence occurred within ∼ 1000 years during MIS 3. This work develops a local MRC applied to new and existing glaciomarine radiocarbon dates, all of which agree that the timing of marine inflow to the grounding line occurred between 15 000 and 14 000 yr BP. Radiocarbon dates paired with sedimentology and existing geomorphology show at least 500 years of ice marginal standstill and substantial grounding zone sedimentation during final ice retreat. We show the Puget Lobe experienced stepwise retreat with GIA and grounding line sedimentation as possible mechanisms for stabilization of the ice margin. While more southern sites (e.g., Double Bluff and Penn Cove) record ice retreat within submarine environments, the northernmost sites (e.g., Cliffside and Rocky Point), appear to record ice retreat into a subaerial environment. The similarities between the rheology in this location and the rheology of the Antarctic Peninsula, as well as the topographic similarities between the Puget Lowland and modern margins of the Greenland Ice Sheet, make these findings highly relevant to increasing process-based understanding of solid Earth influence on ice dynamics in contemporary marine-terminating glacial systems.