Reversible ice sheet thinning in the Amundsen Sea Embayment during the Late Holocene
Cosmogenic-nuclide concentrations in subglacial bedrock cores show that the West Antarctic Ice Sheet (WAIS) at a site between Thwaites and Pope glaciers was at least 35 m thinner than present in the past several thousand years and then subsequently thickened. This is important because of concern that present thinning and grounding line retreat at these and nearby glaciers in the Amundsen Sea Embayment may irreversibly lead to deglaciation of significant portions of the WAIS, with decimeter- to meter-scale sea level rise within decades to centuries. A past episode of ice sheet thinning that took place in a similar, although not identical, climate was not irreversible. We propose that the past thinning–thickening cycle was due to a glacioisostatic rebound feedback, similar to that invoked as a possible stabilizing mechanism for current grounding line retreat, in which isostatic uplift caused by Early Holocene thinning led to relative sea level fall favoring grounding line advance.
During the 2019–2020 Antarctic field season, we used a lightweight minerals exploration drill adapted for polar use (Boeckmann et al., 2021) to collect four bedrock cores from 30–40 m beneath the ice sheet surface on the northern side of the Mount Murphy massif, a volcanic edifice located between Thwaites and Pope glaciers in the Amundsen Sea Embayment of West Antarctica (Figs. 1, 2). These and nearby glaciers are experiencing rapid grounding line retreat and thinning, with the result that this sector of Antarctica is the largest contributor to global sea level rise at present (The IMBIE team, 2018; Smith et al., 2020). Retreat and thinning of major Amundsen Sea glaciers is of concern because these glaciers are retreating into overdeepened basins in which ice is grounded well below sea level. Mercer (1978) and, subsequently, many others have proposed that this geometry is inherently unstable because, in general, the rate of ice loss at a marine ice margin increases with the water depth; thus, ice margin retreat into an overdeepening will result in a runaway increase in the retreat rate that continues until the entire overdeepening is fully deglaciated. We refer to this scenario as “irreversible retreat”. As the overdeepening in this case underlies much of the present West Antarctic Ice Sheet (WAIS), this would imply that the present episode of ice thinning and grounding line retreat in this region will inevitably lead to the loss of a large portion of the WAIS, with globally significant sea level impacts. Thus, recent modeling studies (Favier et al., 2014; Joughin et al., 2014; Rosier et al., 2021) have focused on whether or not further retreat inboard of present grounding line positions implies irreversible deglaciation of large segments of West Antarctica. These model simulations and others predict accelerating retreat in the coming decades to centuries, which is consistent with the hypothesis that present grounding line positions are unstable and that irreversible retreat is underway.
The purpose of the drilling program was to apply luminescence dating and cosmogenic-nuclide exposure dating of subglacial bedrock to establish whether or not similar thinning episodes have taken place in the geologically recent past. There is some evidence (Johnson et al., 2022) that portions of the Antarctic ice sheets were thinner than present, with grounding lines inboard of present positions, in the last several thousand years. This implies that past thinning and retreat of grounding lines inboard of present positions in at least some parts of Antarctica were followed by readvance and thickening to the present configuration, rather than continued and irreversible retreat. Thus, determining if and when Amundsen Sea grounding lines have been inboard of their present positions during the geologically recent past is important for understanding the significance of present thinning and retreat. In coastal regions of Antarctica near present grounding lines, ice thickness is closely coupled to grounding line position; therefore, ice thickness changes, which can be investigated by subglacial bedrock exposure dating, are a potential proxy for grounding line position.
Luminescence dating and cosmogenic-nuclide exposure dating of subglacial bedrock are complementary, but not equivalent, methods. Luminescence dating detects past exposure of mineral grains to sunlight (Aitken, 1998). Thus, it can be used to determine whether a rock surface that is ice covered now was ice-free in the past. Detection of past ice sheet thinning using luminescence dating, however, is only possible if the rock surface became completely ice-free (or nearly so, depending on the optical properties of the basal ice). Cosmogenic-nuclide exposure dating, on the other hand, relies on trace nuclides that are produced by cosmic-ray interactions with rocks that are near, but not necessarily exactly at, the surface. Detectable cosmic-ray production extends up to several meters below the ice surface (Dunai, 2010). Thus, it is possible to detect that ice was thinner in the past using cosmogenic-nuclide measurements on subglacial bedrock, even if the rock surface was not completely ice-free.
A much more common use of cosmogenic-nuclide exposure dating is to estimate the deglaciation age of rock surfaces that are ice-free now and were exposed by past deglaciation. This approach has been widely applied to glacial deposits preserved above the present ice surface throughout Antarctica to reconstruct ice sheet change during times in the past when ice sheets were larger than they are now (e.g., Stone et al., 2003). Reconstructing ice sheet change during warm periods when ice sheets were smaller, on the other hand, is necessary for predicting the response to projected climate warming, but it is more difficult because the evidence that an ice sheet was smaller is inaccessible beneath the present ice sheet (Johnson et al., 2022). Here, we show that this obstacle can be overcome by a targeted drilling program.
Exploiting past exposure of subglacial bedrock as a means of detecting ice sheet thinning is only possible at a site where bedrock is above sea level and thinning will progressively increase the extent of exposed rock adjacent to the ice margin (Spector et al., 2018). The Mount Murphy massif is the only such location in the vicinity of Thwaites and Pope glaciers. On the north side of the massif, several ice-free ridges descend towards Pope Glacier and Crosson Ice Shelf (Fig. 1). At one of these ridges, the north ridge of Kay Peak, a subsidiary peak of the Mount Murphy massif (henceforth “Kay Peak Ridge”), we collected a series of bedrock samples on the exposed part of the ridge from the ice margin to 200 m above the ice surface as well as several bedrock cores from beneath 30–40 m of ice on the subglacial extension of the ridge.
Observations in recent decades show that the ice thickness at Kay Peak Ridge is likely coupled to the position of the grounding line of Pope Glacier. Like Thwaites Glacier, Pope Glacier has experienced rapid grounding line retreat in recent decades, from immediately adjacent to our drill site in 1992–1995 to ∼10 km landward of the site in 2014 (Fig. 1; Milillo et al., 2022). Correspondingly, comparison of 1966 aerial photographs with 2018–2019 satellite imagery indicates expansion of the area of exposed rock on the ridge, which requires thinning of adjacent ice by at least 31 m in 52 years (Fig. 3). Lacking surface mass balance data for the 1966–2020 period, it is not possible to conclusively prove that ice thinning at the drill site is solely a dynamic response to nearby grounding line retreat, but the coincidence of thinning and retreat is consistent with a grounding line control on ice thickness in this area.
Exposure-age data from the ice-free portion of the ridge, collected in previous research (Johnson et al., 2020) and supplemented by additional measurements in this work (Methods 2 and Tables S1 and S5 in the Supplement), include measurements of 14C and 10Be produced in situ in quartz (Fig. 4). The important difference between 14C and 10Be exposure ages is that the 5700-year half-life of 14C precludes the survival of any significant 14C concentration from surface exposure prior to the Last Glacial Maximum (LGM) 15 000–25 000 years ago. Although measured 10Be concentrations could be, and often are, the result of surface exposure during both the present interglaciation and one or more previous ones, 14C can only be the result of post-LGM exposure (Nichols et al., 2019). As described by Johnson et al. (2020), 10Be exposure ages from bedrock and erratics on Kay Peak Ridge are spuriously old due to 10Be inheritance, but 14C exposure ages from the ridge agree with a large set of other data from the overall Mount Murphy massif in recording progressive exposure of the ridge by Early Holocene ice sheet thinning from 200 m above the present ice surface at 6–10 ka to near present conditions at 4–7 ka. The 10Be exposure ages from a nearby small outcrop 1 km to the west of Kay Peak Ridge (Fig. 4; Adams et al., 2022) also indicate deglaciation at 6–7 ka. Although the 14C exposure ages from Kay Peak Ridge scatter more than expected from measurement uncertainties, likely because of the complex effect of both regional ice surface lowering and the retreat of fringing ice fields in uncovering the ridge axis, they agree with other exposure-age data sets throughout the Amundsen Sea region that uniformly show rapid Early and Middle Holocene thinning reaching the present ice surface between 4 and 7 ka (Johnson et al., 2014, 2020). Early Holocene thinning was associated with, and presumably caused by, the coeval retreat of glacier grounding lines across the adjacent continental shelf (Hillenbrand et al., 2017).
There are no exposure-age data in the Amundsen Sea region indicating thicker than present ice after 4 ka (Johnson et al., 2022). Thus, existing observations require either (i) zero change in ice thickness in the last several thousand years or (ii) continued thinning below the present ice surface, followed by thickening to the present configuration. Zero change in ice thickness for millennia appears unlikely given dynamic Late Holocene boundary conditions, including relative sea level (RSL) change forced by eustatic and glacioisostatic effects, climate and oceanographic changes (Walker and Holland, 2007; Hillenbrand et al., 2017), and changes in grounding line position elsewhere in Antarctica (Venturelli et al., 2020; King et al., 2022). In addition, apparent 14C exposure ages of 4–7 ka on samples below the 1966 ice surface (Fig. 4) require at least some period in the Holocene during which the ice was thinner than it was in 1966.
Relative sea level in coastal Antarctica rose between the LGM and Early Holocene because meltwater addition to the ocean from disintegrating LGM ice sheets outpaced local isostatic rebound, but this relationship reversed in the Late Holocene as global ice volume stabilized but isostatic equilibration to Early Holocene mass loss was not yet complete (Lambeck et al., 2014). As relative sea level is a primary control on grounding line position at marine ice margins (e.g., Gomez et al., 2010), this sequence of events predicts grounding line advance and associated thickening in the Late Holocene (Johnson et al., 2022), and isostatic response to ice sheet thinning has also been proposed as a potential stabilizing feedback for present grounding line retreat and thinning (e.g., Barletta et al., 2018). The aim of the drilling program was to determine if this process was active in the past by detecting and, if present, quantifying a hypothesized Holocene thinning–thickening cycle.
We located the subglacial extension of Kay Peak Ridge using ice-penetrating radar. Ice within ∼200 m of the lowest exposed rock on the ridge was crevassed and unsafe for drilling, so the survey began at the closest accessible point to the exposed ridge, where ice was ∼25 m thick, and extended away from the ice margin until the ridge crest was below the depth feasible for drilling (∼60 m). This identified an approximately 70 m long section of subglacial ridge where drilling was possible (Figs. 2, 5).
We conducted the survey with a ski-towed radar system consisting of a Geophysical Survey Systems, Inc. (GSSI) SIR 4000 control unit and model 3108 100 MHz monostatic transceiver, coupled to a Garmin GPSMap 78 handheld GPS unit. We also installed a series of snow stakes in a 70 m × 100 m grid with 10 m spacing overlying the subglacial ridge and marked these stakes on radar profiles collected along grid lines (Fig. 5). This provided higher-resolution spatial referencing than possible using the GPS alone. The antenna was oriented perpendicular to the direction of travel and towed at approximately 1 km h−1. Scans were recorded for 2400 ns with 1024 32-bit samples per scan and 30–50 scans per second, translating to ∼90 scans per meter. Data were recorded with range gain and post-processed with (i) high-pass (25 MHz) and low-pass (300 MHz) filters to reduce noise, (ii) distance corrections based on the 10 m spaced grid stakes, and (iii) stacking (5–10 stacks) to improve the signal-to-noise ratio. We then applied a basic Kirchhoff migration of the data using hyperbolae from bedrock and internal reflectors to estimate relative permittivity as well as associated velocities and depths to bedrock. We did not apply more advanced processing techniques in the field because the survey was conducted immediately (1 d) before we began drilling. Hence, we selected drill sites as close as possible to the ridge axis by locating the apparent high point on each ridge crossing in the migrated data.
We collected subglacial bedrock samples using a two-stage drilling process. First, we made an access hole through the ice, extending as close as possible to the bedrock surface, using a Badger-Eclipse cable-suspended ice core drill, which produces a 113 mm borehole. Drilling with this system stopped when the sediment concentration in basal ice became high enough to dull the cutting bit and halt progress. We then occupied the access hole with a “Winkie” pipe-string, normal circulation compact rock coring drill, modified for polar use by the US Ice Drilling Program (Boeckmann et al., 2021). We set casing in the ice at the bottom of the access hole using an inflatable packer to enable fluid circulation and cored through the remaining sediment-rich basal ice into underlying bedrock. The ice–bedrock interface was frozen at all core sites, and we did not observe loss of drilling fluid into the ice–bed interface.
We reached bedrock 36–41 m below the surface at four boreholes (Fig. 2, Table S1). In all boreholes, we observed a firn–ice transition 12–17 m below the surface and then a transition from clean glacier ice to a 40–120 cm layer of sediment-rich basal ice in contact with bedrock. Note that the presence of debris-rich, optically opaque basal ice is important in this study because resetting of luminescence signals could not occur through thin ice; detection of ice thinning using luminescence dating would require complete removal of basal ice and direct sunlight exposure of the bedrock surface. At three holes, we recovered 114–137 cm long, 33 mm diameter cores of the underlying bedrock, whose lithology is a biotite gneiss identical to the bedrock exposed at the ice margin. At the fourth borehole (19-KP-H6) we recovered 7 cm of fragmented rock with the same lithology. However, this section of the core included fragments of both rock and debris-rich basal ice, so it was not possible to determine whether this sample represented (1) part of a densely jointed or rubbly bedrock surface or (2) one or more loose clasts within the basal ice. At two of the boreholes, we were able to extract bedrock core tops from the core collection barrel within a lightproof bag so that they were not exposed to sunlight, permitting infrared stimulated luminescence (IRSL) measurements on the bedrock surface. In the next sections, we describe IRSL measurements on two core tops and measurements of the cosmic-ray-produced radionuclides 10Be, 14C, and 26Al in quartz from three cores.
We made IRSL depth profile measurements (Sohbati et al., 2011, 2012) in the uppermost 4 cm of the two core top samples that we were able to remove from the core barrel without exposure to sunlight (19-KP-H1 and 19-KP-H4). In addition, to establish kinetics of bleaching during sunlight exposure for this lithology, we collected equivalent depth-profile data on a rock surface above the present ice margin that had been newly exposed by collection of a rock sample 4 years earlier (Johnson et al., 2020).
Depth-profile measurements were made by slicing rock samples parallel to the exposed surface using a low-speed, water-cooled lapidary saw with a 0.38 mm fine-kerf blade. Each slice was then approximately broken into quarters, each of which was disaggregated using gentle crushing in a mortar and pestle to yield mineral grains suitable for IRSL measurement, and all were analyzed separately in order to better capture the natural variability in the luminescence response at each depth. We then used a post-infrared infrared-stimulated luminescence (pIRIR) protocol (Buylaert et al., 2009) to determine sensitivity-corrected natural luminescence intensities () as well as the dose response and fading characteristics of each sample. This pIRIR protocol yields measurements of two luminescence signals: the more bleachable IR50 signal, which originates from both Na- and K-rich feldspars and has higher signal-fading rates, and the less bleachable pIRIR225 signal, which is mostly from K-rich feldspars and fades less (Thomsen et al., 2018). Further details of the measurements are given in Methods 1 and Table S2 in the Supplement.
The depth profile in the surface sample shows significant bleaching, manifested as values of for both luminescence signals that approach zero at the surface and are significantly below field saturation values in the uppermost 1 cm, after 4 years of exposure (Fig. 6). If the core tops had been exposed to sunlight during the Holocene, they would also show values of systematically below saturation at shallow depths. They do not. Rather, IRSL signals in core tops are indistinguishable from calculated field saturation at all depths. This means that the bedrock surfaces of the core tops were not exposed to sunlight during the period required to approach saturation. We estimate the period required to approach saturation as 200–280 kyr based on estimated dose rates and measured fading properties (Methods 1 in the Supplement). Therefore, bedrock at these two core sites was not ice-free during the Holocene.
We sectioned bedrock cores into 10 cm segments, separated quartz from each segment, and measured concentrations of cosmic-ray-produced 10Be and 14C in quartz separates from three cores as well as 26Al in quartz from one core. As discussed above in the introduction, the purpose of measuring multiple cosmogenic nuclides in the core samples is to distinguish Holocene from pre-Holocene exposure of bedrock. 10Be and 26Al have relatively long half-lives (1.4 and 0.7 Ma, respectively); thus, if the core sites were not affected by subglacial erosion during the LGM, 26Al and 10Be produced in near-surface rocks in a previous ice-free period prior to the LGM could be preserved to the present. Observing cosmogenic 10Be or 26Al in subglacial bedrock would prove that the surface had been exposed at some time in the past, but it would not prove that the surface had been exposed in the Holocene. In contrast, 14C has a half-life of 5700 years, so 14C produced prior to the LGM would be removed by radioactive decay. Cosmogenic 14C observed in subglacial bedrock unambiguously requires Holocene exposure.
In brief, 26Al and 10Be analysis employs an isotope dilution method in which quartz is dissolved in hydrofluoric acid in the presence of a known quantity of the common isotopes 9Be and 27Al, Be and Al are separated by column chromatography, and 10Be and 26Al concentrations are quantified by accelerator mass spectrometer (AMS) measurement of the and ratios. 14C produced in situ in quartz is measured by fusing quartz with a lithium metaborate flux in a high-purity O2 atmosphere to release 14C as CO2, followed by CO2 purification, barometric measurement of the total amount of CO2, and conversion to graphite for AMS measurement of the ratio. The detection limit for all three nuclides is quantified by measurement of process blanks. The 10Be and 26Al concentrations in core samples are 2000–7000 and 16 000–44 000 atoms g−1, 2–5 times the detection limit inferred from process blanks and within the range routinely measured for exposure-dating applications. The 14C concentrations in quartz in the core samples are in the range of 10 000–50 000 atoms g−1. These concentrations are 2–3 times the typical detection limits for in situ 14C in quartz, but they are significantly lower than routinely measured. Thus, we characterized the distribution of process blanks by analyzing 50 process blanks during the period in which sample analysis took place. In addition, we made replicate analyses of multiple aliquots of quartz for nearly all samples. Analytical methods, blank-correction procedures, and determination of detection limits are described in detail in Methods 2 and Tables S3–S5 in the Supplement as well as in Balco et al. (2022).
10Be and 26Al concentrations decrease with depth in all cores (Fig. 7). Thus, 10Be and 26Al concentrations are the result of near-surface cosmic-ray production. Although 10Be and 26Al concentrations in one core conform to the production ratio and are therefore consistent with recent exposure, the presence of cosmic-ray-produced 10Be and 26Al alone, as discussed above, does not necessarily show that ice at this site was thinner during the Holocene. Observed concentrations could have been produced during an interglacial period prior to the LGM and not subsequently removed by subglacial erosion.
Concentrations of 14C are also above background in core samples and inversely correlated with depth in all cores (Fig. 7), again requiring near-surface cosmic-ray production. The 10Be and 26Al data require cosmic-ray exposure at some point, not necessarily in the Holocene, but the 14C data require Holocene cosmic-ray exposure.
In contrast to the IRSL data, the cosmogenic-nuclide data require that ice at the core sites was thinner sometime during the Holocene than it is now. We show this quantitatively by considering a null hypothesis that (i) any preexisting nuclide inheritance was removed by subglacial erosion during the LGM; (ii) the ice thinned rapidly during the Early to Middle Holocene, as indicated by exposure-age data from ice-free areas in the region; (iii) the ice thickness reached its present value between 4 and 7 ka; and (iv) the ice thickness did not change subsequently. In this hypothesis, ice overlying core sites was never thinner than present, so 10Be, 26Al, and 14C production rates were never higher than ∼0.01, ∼0.1, and ∼0.3 atoms g−1 yr−1, their respective values under ∼40 m of ice. The hypothesis predicts 10Be, 26Al, and 14C concentrations in the range of 100–200, 900–1500, and 1300–2600 atoms g−1, respectively (Fig. 7), and is thus clearly rejected by the observations. Therefore, the ice sheet was thinner during the Holocene than it is now.
The cosmogenic-nuclide data from the cores, with exposure-age data from the ice-free portion of the ridge and observations during drilling, also provide information on whether or not subglacial erosion at the core site is possible. In general, subglacial erosion is an important consideration for subglacial bedrock exposure dating because even if a rock surface was ice-free in the past, if several meters of subglacial erosion took place during subsequent ice cover, the cosmogenic-nuclide inventory produced during the ice-free period could be removed. Thus, failing to find cosmogenic-nuclide concentrations above background in core samples does not prove that the ice sheet was never thinner at the core site. However, the opposite is not true: finding cosmogenic-nuclide concentrations above background in core samples, as we have done here, not only requires that the ice was thinner at the core site in the past but also strictly limits the possible amount of later subglacial erosion. Thus, the presence of 14C above background in core samples requires Holocene thinning and precludes significant subsequent erosion. This is not surprising, as subglacial erosion is effectively glaciologically impossible at this site under Holocene conditions. At present, ice at the core sites is frozen to the bed and, given the mean annual temperature at Kay Peak of near −15 ∘C, basal melting would only be predicted if the ice thickened by hundreds of meters. Glacial-geological observations and exposure-age data from other sites in coastal West Antarctica (Sugden et al., 2005; Balco et al., 2013), where the LGM ice surface was hundreds of meters above present, show this clearly: evidence for basal melting and subglacial erosion during the LGM is only present near the present ice margin where LGM ice cover was thickest, and it is not present at higher elevations, implying that the freezing isotherm at the LGM was near the elevation of the present ice margin. Likewise, at Kay Peak Ridge, the LGM ice surface was hundreds of meters higher than present, but there is no geomorphic evidence of erosion on the ice-free portion of the ridge, and pre-Holocene 10Be exposure ages in erratics and bedrock further preclude subglacial erosion. However, sediment-rich basal ice at the core sites requires entrainment and transport of subglacial sediment at some point, which likely indicates that basal ice was close to the pressure-melting point during LGM conditions. Taking all of this into account, we conclude the following. First, as at other coastal sites in Antarctica, the coincidence of evidence precluding glacial erosion above the present ice margin with debris-rich basal ice at core sites indicates that the basal melting isotherm was near the present ice margin when the ice was hundreds of meters thicker at the LGM. Second, Late Holocene thickness changes that were likely in the range of tens of meters rather than hundreds could not result in basal melting at the core sites. Thus, subglacial erosion did not occur at the core sites after the lowstand required by the cosmogenic-nuclide data.
To quantify Holocene ice thickness change at the core sites, we used a constrained random search algorithm, coupled to a forward model for cosmogenic-nuclide production, to identify ice thickness change histories that fit the observations.
The input to the forward model is a history of ice surface elevation change at a sample site, expressed as a series of (t,H) pairs, where t is the time in years before present and H is the ice surface elevation in meters above sea level. Surface elevation H(t) is converted to an ice and/or firn thickness , where H0 is the current elevation of the sample. This procedure implicitly disregards glacioisostatic elevation change and assumes that the sample elevation is constant. It would be possible to estimate glacioisostatic elevation change by reference to model simulations or to relative sea level data from Pine Island Bay (Braddock et al., 2022) indicating approximately 20 m of relative sea level fall during the Late Holocene. However, elevation change of this magnitude would imply no more than a 2.5 % variation in cosmogenic-nuclide production rates, which is not significant.
The mass thickness of ice and firn z(t) (g cm−2), which is necessary to compute nuclide production rates when a sample is ice covered, is obtained from the following relationship:
where Z(t) has units of centimeters, ρice is the density of glacier ice (0.917 g cm−3), and D (680 g cm−2) and L (1130 cm) are constants obtained from fitting this equation to ice and firn densities measured in the field during drilling (Fig. 8).
Cosmogenic-nuclide production rates P(H,z) (atoms g−1 yr−1) are calculated from surface elevation H(t) and mass thickness z(t) as the sum of near-surface production by high-energy neutron spallation and production by deeply penetrating muons. Spallogenic production is exponential with an attenuation length of 140 g cm−2 (Balco, 2017). Surface production rates due to spallation use the scaling method of Stone (2000) and the Antarctic atmosphere approximation of Radok et al. (1996), as implemented by Balco et al. (2008). Spallogenic 10Be and 26Al production is calibrated using the “primary” calibration data set of Borchers et al. (2016). Spallogenic 14C production is calibrated from measurements at Tulane of the 14C concentration in the CRONUS-A intercomparison standard (Jull et al., 2015), which was collected from a site in the Antarctic Dry Valleys, and the assumption of production–decay equilibrium at that sample site (Balco et al., 2019; Goehring et al., 2019; Nichols et al., 2019). Production rates due to muons use Model 1A of Balco (2017) with cross-sections for 10Be, 26Al, and 14C production as calibrated in that reference. Given the production rate P(H(t),z(t)), the nuclide concentration at the present time is the integral:
where zs is the mass depth of a core sample below the bedrock surface (g cm−2) and λ is the decay constant (year) for the nuclide in question. Samples are assumed to have a zero nuclide concentration at the beginning of the elevation change history (t=tmax). For computational efficiency, P(H,z) is precalculated on a 2 d mesh in H and z and is used as a lookup table during integration. Interpolation in z is log-linear to improve accuracy. We used the default numerical integration routine provided by MATLAB numerical computation software.
The random search algorithm generates piecewise-linear age–elevation histories H(t) that are constrained by geological and geochronological observations, applies the forward model for cosmogenic-nuclide production to compute predicted nuclide concentrations in the core samples corresponding to the age–elevation history, and then compares predicted to observed concentrations to identify H(t) that fit the observations. The randomly generated age–elevation histories have the following properties. First, the ice surface is at 350 m elevation between 35 ka and the Early Holocene (Johnson et al., 2020). An elevation of 350 m is chosen because it is slightly above the elevation of the highest glacial erratics with post-LGM exposure ages found on Kay Peak Ridge (Johnson et al., 2020). The LGM surface may have been higher, but the exact value chosen is not significant in this calculation because cosmogenic-nuclide production rates under >100 m of ice are insignificant in our calculation. Second, the beginning of deglaciation is drawn from a uniform random distribution between 9 and 6 ka. These limits are selected to be consistent with regional exposure-age data (Johnson et al., 2020; Adams et al., 2022). Third, the elevation change history after initial deglaciation consists of three piecewise-linear segments between the beginning of deglaciation and the present time. These are generated by selecting two times from a uniform random distribution within this interval as well as two surface elevations from a uniform random distribution between the elevation of the top of sediment-rich ice at a core site (0.7–1.2 m above the bedrock surface) and 100 m elevation (20–25 m above bedrock at the core sites). The 100 m upper limit was chosen after initial experiments to exclude regions of the parameter space that would be certain to yield poor fits to the data, thus improving the efficiency of the random search algorithm. Constraining the minimum thickness to be the thickness of sediment-rich ice observed at present is based on the assumption that the dirty basal ice is the result of subglacial erosion at the LGM and that if the ice had been thinner following the LGM, this unit would no longer be present. Fourth, the ice at the core site reaches its present (2019–2020) elevation at t=0. Omitting post-1966 thinning to simplify this constraint does not have a significant effect on the results. Finally, following the reasoning in the previous section, subglacial erosion during the Holocene is assumed to be zero. Note that the assumption that nuclide concentrations are zero at t=tmax means that it is not necessary to specify subglacial erosion during LGM conditions. Alternatively, the zero initial condition assumption could be considered an assumption that meter-scale erosion did take place during the LGM.
We quantify the fit of predicted to observed 10Be and 14C concentrations in cores using an error-weighted sum-of-squares statistic M:
where Nm,i is the measured nuclide concentration (atoms g−1) in analysis i, Np,i is the nuclide concentration predicted by the forward model for the sample on which analysis i was made, and σNm,i is the uncertainty in the measured nuclide concentration. This formulation implies a normal uncertainty distribution for Nm,i. As this is not the case for the 14C measurements (Methods 2 in the Supplement), we approximate the uncertainty distribution for 14C measurements as an asymmetrical normal distribution with different values of σNm,i above and below. Because 26Al measurements were only made in one core (19-KP-H1), we use only 10Be and 14C data in computing M, thereby standardizing the calculation for all cores. However, as discussed below, 26Al data are exactly consistent with the 10Be measurements, so including or excluding them in calculating M for core 19-KP-H1 does not affect the results. If we assume that all measurement uncertainties are correctly characterized, the model correctly represents all operative processes, and predicted and observed values differ only due to measurement uncertainties, then the expected value of M is near 1. However, these assumptions are not completely correct. For example, we have disregarded uncertainty in production rates due to muon interactions, which may be 10 % or greater (this was omitted because production rate uncertainties are correlated between core samples, which complicates interpretation of a misfit statistic). Thus, we view minimum values of M near 2 generated by the random search algorithm as acceptable fits. All results of the random search algorithm are shown in Methods 3 in the Supplement.
Ice thickness change histories generated by the random search algorithm that fit the bedrock core data uniformly include a period of several thousand years during which the ice was several meters thick at the core sites (Figs. 9 and 7 as well as Methods 3 in the Supplement). Ice thickness during this lowstand is in the range of 2–5, 3–6, and 4–7 m at 19-KP-H1, 19-KP-H4, and 19-KP-H5, respectively, implying 30–35 m of thinning relative to present. Although we only assessed the model fit against data from the bedrock cores and not against the exposure-age data from the exposed ridge, scenarios that fit the subsurface data follow the older bound of the surface data, which is consistent with the hypothesis that scatter in these data is the result of locally derived snow and ice fields that persisted after ice surface lowering. Scenarios with shorter lowstands, or thicker ice during the lowstands, predict systematically lower cosmogenic-nuclide concentrations than observed (Fig. 7).
Ice thickness change histories that fit the measured nuclide concentrations display a trade-off between lowstand duration and ice thickness during the lowstand. A longer lowstand duration predicts higher nuclide concentrations, and thicker ice during the lowstand predicts lower nuclide concentrations. A longer-duration lowstand with thicker ice therefore yields a similar fit to data as a shorter-duration lowstand with thinner ice during the lowstand. Thus, it is not possible to choose a single best-fitting thickness change history or to accurately date the beginning and end of the lowstand. For example, a 4000-year period with 3 m thick ice fits the data for core H1 as well as a 7000-year period with 5 m thick ice (Fig. 9.). However, all cores require a lowstand period of no less than 3000 years during which ice was no thicker than several meters (Figs. 7, 9).
Finally, the fact that the same Late Holocene lowstand scenarios are inferred from fitting 10Be, 14C, and 26Al measurements separately (Fig. 10) precludes the possibility that any significant fraction of the measured 10Be or 26Al inventory was inherited from pre-Holocene exposure, perhaps in an interglacial period prior to the LGM. If this were the case, corresponding inherited 14C would have been removed by radioactive decay, and a longer period of exposure would be required to fit the 10Be and 26Al data than the 14C data. Subglacial erosion during the last glacial cycle, as also suggested by the presence of debris-rich basal ice and IRSL data precluding sunlight exposure in the last 200–300 ka, likely removed rock surfaces exposed during prior interglaciations.
To summarize, although IRSL data show that bedrock at the core sites was not ice-free in the Holocene, the cosmogenic-nuclide data show that ice at the core sites was 30–35 m thinner than present for at least 3000 years, and possibly as long as 8000 years, during the Middle to Late Holocene. The coincidence of IRSL data that preclude direct exposure of rock surfaces to sunlight and cosmogenic-nuclide data that require significant ice thinning is consistent with the present geometry of ice and rock on the exposed portion of Kay Peak Ridge. Locating the ridge axis using single-channel ice-penetrating radar, as we did, is not expected to be precise at the meter scale because of limits on instrument resolution and the likelihood of off-axis reflections from irregular bedrock topography. The currently exposed portion of the ridge is mantled by ice fields on both sides that extend to within meters of the ridge axis and thicken with distance from the axis (Fig. 2). In addition, lowstand ice thicknesses inferred from the random search algorithm vary between cores over a range of 2–3 m (Fig. 9), but the core top elevations vary over a range of 6 m. These observations suggest that ice with a thickness of several meters present during the lowstand does not reflect a flat ice surface uniformly 35 m below the present surface but instead reflects a condition like what we observe now, in which fringing ice fields of varying thickness covered drill sites that are slightly off the ridge axis. This, in turn, indicates that the estimate of 35 m of thinning is a minimum. Areas farther from the ridge axis nearly certainly thinned more.
Our subglacial bedrock exposure-dating results are direct, unambiguous evidence for ice thinning and subsequent thickening of at least 35 m during the past several thousand years. Ice thinning below present conditions earlier in the Holocene was not followed by further runaway thinning and significant regional ice loss. Instead, it was followed by thickening to the present ice sheet configuration. The significance of this observation with respect to the question of whether or not grounding line retreat and thinning observed in recent decades is or is not irreversible depends on whether or not the ice thickness at Kay Peak is an accurate proxy for grounding line position. If one accepts that simultaneous thinning at Kay Peak and retreat of the Pope Glacier grounding line observed in recent decades are dynamically linked, then our observation that mid-Holocene thinning below the present ice thickness was reversible implies that mid-Holocene grounding line retreat inboard of present positions was also reversible. By analogy with present conditions, thinning and subsequent thickening of 35 m or more may have been associated with grounding line retreat tens of kilometers upstream of present locations, and subsequent readvance.
These observations, with independent evidence for falling relative sea level in the nearby Pine Island Bay region during the past 4000 years (Braddock et al., 2022), most likely indicate that a grounding line retreat–advance cycle driven by glacioisostatic rebound–RSL feedback took place during the Holocene. Observations from elsewhere in Antarctica (Bradley et al., 2015; Kingslake et al., 2018; Venturelli et al., 2020; King et al., 2022; Johnson et al., 2022) indicate that this likely took place at many marine ice margins around the continent. In general terms, therefore, our results provide supporting evidence for the hypothesis that glacioisostatic rebound can provide an important stabilizing feedback on grounding line retreat. However, boundary conditions during a Late Holocene retreat–readvance cycle were not identical to present conditions. Applying evidence for a Holocene retreat–readvance cycle to determine the importance of glacioisostatic rebound feedback in stabilizing present grounding line retreat in the Amundsen Sea region also requires understanding the following: (i) the difference between isostatic responses to large-scale LGM-to-Holocene thinning that took place over thousands of years and thinning of lesser magnitude observed in recent decades; (ii) the relative importance (to the grounding line position) of relative sea level change, ice shelf buttressing, and sub-ice-shelf melting driven by oceanographic conditions, which cannot be deduced from any data presented here; and (iii) past and present bed topography in the region of grounding lines. In addition, the past thinning–thickening cycle that we have detected spanned a duration of at least 3000 years, which is geologically rapid but slow by comparison to the timescale of projected sea level rise impacts of present ice sheet thinning.
All data described in this paper that have not already been published elsewhere are included in the Supplement. The MATLAB code used to carry out all calculations in this paper is also included in the Supplement.
The supplement related to this article is available online at: https://doi.org/10.5194/tc-17-1787-2023-supplement.
The author contributions, following the CRediT authorship guidelines, are as follows: GB, BG, JSJ, DHR, and JW: conceptualization; GB, NB, BG, JSJ, KN, DHR, RAV, and JW: methodology; GB, NB, KN, and RAV: software; GB, NB, BG, JSJ, KN, DHR, RAV, and KW: validation; GB, NB, BG, KN, DHR, RAV, and KW: analysis; JA, GB, NB, SC, BG, JSJ, KN, DHR, RAV, and KW: investigation; GB, SC, BG, JSJ, DHR, and KW: resources; GB, NB, SC, BG, KN, DHR, RAV, and KW: data curation; GB, NB, KN, and RAV: original draft; JA, GB, NB, SB, SC, BG, BH, JSJ, KN, DHR, RAV, JW, and KW: review and editing; GB, NB, and SC: visualization; GB, BG, JSJ, and DHR: supervision; GB, BG, JSJ, KN, DHR, RAV, and JW: administration; GB, SC, BG, BH, JSJ, DHR, JW, and KW: funding acquisition.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work would not have been possible without the efforts of many members of the US Antarctic Program and the British Antarctic Survey, in particular Chris Simmons, Leslie Blank, and Nick Gillett; pilots and crew of Kenn Borek Air, Ltd. and the 109th Airlift Wing of the New York Air National Guard; the US Ice Drilling Program, in particular Grant Boeckmann and Eliot Moravec; the Polar Geospatial Center of the University of Minnesota; and the National Ocean Sciences AMS facility, in particular Mark Roberts. This work is part of the Geological History Constraints project, which is a component of the International Thwaites Glacier Collaboration (ITGC). The US Ice Drilling Program is supported by NSF Cooperative Agreement 1836328, and the Polar Geospatial Center is supported by NSF grant no. OPP-2029685. This paper has LANL publication release no. LA-UR-22-32933 and is ITGC contribution no. 66.
This research has been supported by the US National Science Foundation Office of Polar Programs (grant nos. OPP-1738989 and EAR-1806629), the UK National Environment Research Council (grant nos. NE/S006710/1, NE/S00663X/1, and NE/S006753/1), and the Ann and Gordon Getty Foundation.
This paper was edited by Nicolas Jourdain and reviewed by Jason Briner and Nathaniel A. Lifton.
Adams, J. R., Johnson, J. S., Roberts, S. J., Mason, P. J., Nichols, K. A., Venturelli, R. A., Wilcken, K., Balco, G., Goehring, B., Hall, B., Woodward, J., and Rood, D. H.: New 10Be exposure ages improve Holocene ice sheet thinning history near the grounding line of Pope Glacier, Antarctica, The Cryosphere, 16, 4887–4905, https://doi.org/10.5194/tc-16-4887-2022, 2022. a, b, c
Aitken, M. J.: Introduction to optical dating: the dating of Quaternary sediments by the use of photon-stimulated luminescence, Clarendon Press, ISBN 0-19-85-4092-2, 1998. a
Balco, G., Stone, J. O., Lifton, N. A., and Dunai, T. J.: A complete and easily accessible means of calculating surface exposure ages or erosion rates from 10Be and 26Al measurements, Quat. Geochronol., 3, 174–195, 2008. a
Balco, G., Schaefer, J. M., and the LARISSA group: Exposure-age record of Holocene ice sheet and ice shelf change in the northeast Antarctic Peninsula, Quaternary Sci. Rev., 59, 101–111, https://doi.org/10.1016/j.quascirev.2012.10.022, 2013. a
Balco, G., Todd, C., Goehring, B. M., Moening-Swanson, I., and Nichols, K.: Glacial geology and cosmogenic-nuclide exposure ages from the Tucker Glacier-Whitehall Glacier confluence, northern Victoria Land, Antarctica, Am. J. Sci., 319, 255–286, 2019. a
Balco, G., Brown, N., Nichols, K., Venturelli, R. A., Adams, J., Braddock, S., Campbell, S., Goehring, B., Johnson, J. S., Rood, D. H., Wilcken, K., Hall, B., and Woodward, J.: Author comment on “Reversible ice sheet thinning in the Amundsen Sea Embayment during the Late Holocene”, The Cryosphere Discussions, https://doi.org/10.5194/tc-2022-172-AC2, 2022. a
Barletta, V. R., Bevis, M., Smith, B. E., Wilson, T., Brown, A., Bordoni, A., Willis, M., Khan, S. A., Rovira-Navarro, M., Dalziel, I., Smalley Jr., R., Kendrick, E., Konfal, S., Caccamise II, D. J., Aster, R. C., Nyblade, A., and Wiens, D. A.: Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability, Science, 360, 1335–1339, https://doi.org/10.1126/science.aao1447, 2018. a
Boeckmann, G. V., Gibson, C. J., Kuhl, T. W., Moravec, E., Johnson, J. A., Meulemans, Z., and Slawny, K.: Adaptation of the Winkie Drill for subglacial bedrock sampling, Ann. Glaciol., 62, 109–117, 2021. a, b
Borchers, B., Marrero, S., Balco, G., Caffee, M., Goehring, B., Lifton, N., Nishiizumi, K., Phillips, F., Schaefer, J., and Stone, J.: Geological calibration of spallation production rates in the CRONUS-Earth project, Quat. Geochronol., 31, 188–198, 2016. a
Braddock, S., Hall, B., Johnson, J., Balco, G., Spoth, M., Whitehouse, P., Campbell, S., Goehring, B., Rood, D., and Woodward, J.: Relative sea-level data preclude major late Holocene ice mass change in Pine Island Bay, Nat. Geosci., 15, 568–572, https://doi.org/10.1038/s41561-022-00961-y, 2022. a, b
Bradley, S. L., Hindmarsh, R. C., Whitehouse, P. L., Bentley, M. J., and King, M. A.: Low post-glacial rebound rates in the Weddell Sea due to Late Holocene ice-sheet readvance, Earth Planet. Sc. Lett., 413, 79–89, 2015. a
Buylaert, J.-P., Murray, A. S., Thomsen, K. J., and Jain, M.: Testing the potential of an elevated temperature IRSL signal from K-feldspar, Radiat. Meas., 44, 560–565, 2009. a
Dunai, T. J.: Cosmogenic nuclides: principles, concepts and applications in the Earth surface sciences, Cambridge University Press, ISBN 978-0-521-87380-2, 2010. a
Favier, L., Durand, G., Cornford, S. L., Gudmundsson, G. H., Gagliardini, O., Gillet-Chaulet, F., Zwinger, T., Payne, A., and Le Brocq, A. M.: Retreat of Pine Island Glacier controlled by marine ice-sheet instability, Nat. Clim. Change, 4, 117–121, 2014. a
Goehring, B. M., Balco, G., Todd, C., Moening-Swanson, I., and Nichols, K.: Late-glacial grounding line retreat in the northern Ross Sea, Antarctica, Geology, 47, 291–294, 2019. a
Gomez, N., Mitrovica, J. X., Huybers, P., and Clark, P. U.: Sea level as a stabilizing factor for marine-ice-sheet grounding lines, Nat. Geosci., 3, 850–853, 2010. a
Hillenbrand, C.-D., Smith, J. A., Hodell, D. A., Greaves, M., Poole, C. R., Kender, S., Williams, M., Andersen, T. J., Jernas, P. E., Elderfield, H., Klages, J. P., Roberts, S. J., Gohl, K., Larter, R. D., and Kuhn, G.: West Antarctic Ice Sheet retreat driven by Holocene warm water incursions, Nature, 547, 43–48, https://doi.org/10.1038/nature22995, 2017. a, b
Johnson, J. S., Bentley, M. J., Smith, J. A., Finkel, R., Rood, D., Gohl, K., Balco, G., Larter, R. D., and Schaefer, J.: Rapid thinning of Pine Island Glacier in the early Holocene, Science, 343, 999–1001, 2014. a
Johnson, J. S., Roberts, S. J., Rood, D. H., Pollard, D., Schaefer, J. M., Whitehouse, P. L., Ireland, L. C., Lamp, J. L., Goehring, B. M., Rand, C., and Smith, J. A.: Deglaciation of Pope Glacier implies widespread early Holocene ice sheet thinning in the Amundsen Sea sector of Antarctica, Earth Planet. Sc. Lett., 548, 116501, https://doi.org/10.1016/j.epsl.2020.116501, 2020. a, b, c, d, e, f, g, h
Johnson, J. S., Venturelli, R. A., Balco, G., Allen, C. S., Braddock, S., Campbell, S., Goehring, B. M., Hall, B. L., Neff, P. D., Nichols, K. A., Rood, D. H., Thomas, E. R., and Woodward, J.: Review article: Existing and potential evidence for Holocene grounding line retreat and readvance in Antarctica, The Cryosphere, 16, 1543–1562, https://doi.org/10.5194/tc-16-1543-2022, 2022. a, b, c, d, e
Joughin, I., Smith, B. E., and Medley, B.: Marine ice sheet collapse potentially under way for the Thwaites Glacier Basin, West Antarctica, Science, 344, 735–738, 2014. a
Jull, A. T., Scott, E. M., and Bierman, P.: The CRONUS-Earth inter-comparison for cosmogenic isotope analysis, Quat. Geochronol., 26, 3–10, 2015. a
King, M. A., Watson, C. S., and White, D.: GPS rates of vertical bedrock motion suggest late Holocene ice-sheet readvance in a critical sector of East Antarctica, Geophys. Res. Lett., 49, e2021GL097232, https://doi.org/10.1029/2021GL097232, 2022. a, b
Kingslake, J., Scherer, R., Albrecht, T., Coenen, J., Powell, R., Reese, R., Stansell, N., Tulaczyk, S., Wearing, M., and Whitehouse, P.: Extensive retreat and re-advance of the West Antarctic Ice Sheet during the Holocene, Nature, 558, 430–434, 2018. a
Lambeck, K., Rouby, H., Purcell, A., Sun, Y., and Sambridge, M.: Sea level and global ice volumes from the Last Glacial Maximum to the Holocene, P. Natl. Acad. Sci. USA, 111, 15296–15303, 2014. a
Matsuoka, K., Skoglund, A., Roth, G., de Pomereu, J., Griffiths, H., Headland, R., Herried, B., Katsumata, K., Le Brocq, A., Licht, K., Morgan, F., Neff, P. D., Ritz, C., Scheinert, M., Tamura, T., Van de Putte, A., van den Broeke, M., von Deschwanden, A., Deschamps-Berger, C., Van Liefferinge, B., Tronstad, S., and Melvær, Y.: Quantarctica, an integrated mapping environment for Antarctica, the Southern Ocean, and sub-Antarctic islands, Environ. Model. Softw., 140, 105015, https://doi.org/10.1016/j.envsoft.2021.105015, 2021. a
Mercer, J. H.: West Antarctic ice sheet and CO2 greenhouse effect: a threat of disaster, Nature, 271, 321–325, 1978. a
Milillo, P., Rignot, E., Rizzoli, P., Scheuchl, B., Mouginot, J., Bueso-Bello, J. L., Prats-Iraola, P., and Dini, L.: Rapid glacier retreat rates observed in West Antarctica, Nat. Geosci., 15, 48–53, 2022. a, b
Nichols, K. A., Goehring, B. M., Balco, G., Johnson, J. S., Hein, A. S., and Todd, C.: New Last Glacial Maximum ice thickness constraints for the Weddell Sea Embayment, Antarctica, The Cryosphere, 13, 2935–2951, https://doi.org/10.5194/tc-13-2935-2019, 2019. a, b
Radok, U., Allison, I., and Wendler, G.: Atmospheric surface pressure over the interior of Antarctica, Antarct. Sci., 8, 209–217, 1996. a
Rignot, E., Mouginot, J., and Scheuchl, B.: Antarctic grounding line mapping from differential satellite radar interferometry, Geophysi. Res. Lett., 38, L10504, https://doi.org/10.1029/2011GL047109, 2011. a
Rosier, S. H. R., Reese, R., Donges, J. F., De Rydt, J., Gudmundsson, G. H., and Winkelmann, R.: The tipping points and early warning indicators for Pine Island Glacier, West Antarctica, The Cryosphere, 15, 1501–1516, https://doi.org/10.5194/tc-15-1501-2021, 2021. a
Smith, B., Fricker, H. A., Gardner, A. S., Medley, B., Nilsson, J., Paolo, F. S., Holschuh, N., Adusumilli, S., Brunt, K., Csatho, B., Harbeck, K. Markus, T., Neumann, T., Siegfried, M. R., and Zwally, H. J.: Pervasive ice sheet mass loss reflects competing ocean and atmosphere processes, Science, 368, 1239–1242, https://doi.org/10.1126/science.aaz5845, 2020. a
Sohbati, R., Murray, A. S., Jain, M., Buylaert, J.-P., and Thomsen, K. J.: Investigating the resetting of OSL signals in rock surfaces, Geochronometria, 38, 249–258, 2011. a
Sohbati, R., Murray, A. S., Chapot, M. S., Jain, M., and Pederson, J.: Optically stimulated luminescence (OSL) as a chronometer for surface exposure dating, J. Geophys. Res.-Sol. Ea., 117, B09202, https://doi.org/10.1029/2012JB009383, 2012. a
Spector, P., Stone, J., Pollard, D., Hillebrand, T., Lewis, C., and Gombiner, J.: West Antarctic sites for subglacial drilling to test for past ice-sheet collapse, The Cryosphere, 12, 2741–2757, https://doi.org/10.5194/tc-12-2741-2018, 2018. a
Stone, J. O.: Air pressure and cosmogenic isotope production, J. Geophys. Res.-Sol. Ea., 105, 23753–23759, 2000. a
Stone, J. O., Balco, G. A., Sugden, D. E., Caffee, M. W., Sass III, L. C., Cowdery, S. G., and Siddoway, C.: Holocene deglaciation of Marie Byrd land, west Antarctica, Science, 299, 99–102, 2003. a
Sugden, D. E., Balco, G., Cowdery, S. G., Stone, J. O., and Sass III, L. C.: Selective glacial erosion and weathering zones in the coastal mountains of Marie Byrd Land, Antarctica, Geomorphology, 67, 317–334, 2005. a
Thomsen, K., Kook, M., Murray, A., and Jain, M.: Resolving luminescence in spatial and compositional domains, Radiati. Meas., 120, 260–266, 2018. a
Venturelli, R., Siegfried, M., Roush, K., Li, W., Burnett, J., Zook, R., Fricker, H., Priscu, J., Leventer, A., and Rosenheim, B.: Mid-Holocene grounding line retreat and readvance at Whillans Ice Stream, West Antarctica, Geophys. Res. Lett., 47, e2020GL088476, https://doi.org/10.1029/2020GL088476, 2020. a, b
Walker, R. T. and Holland, D. M.: A two-dimensional coupled model for ice shelf–ocean interaction, Ocean Model., 17, 123–139, 2007. a