Comparing Glacial-Geological Evidence and Model Simulations of Ice Sheet Change since the Last Glacial Period in the Amundsen Sea Sector of Antarctica

the Antarctic Ice has undergone extensive changes, resulting in a much smaller present-day configuration. Improving our understanding of basic physical processes that played important roles during that retreat is critical to providing more robust model projections of future retreat and sea-level rise. Here, a limited-area nested ice sheet model was applied to the last deglacial retreat of the West Antarctic Ice Sheet in the Amundsen Sea Embayment (ASE), at 5 km resolution. The ice sheet response to climate and sea-level forcing was examined at two sites along the flowlines of Pine Island Glacier and Pope Glacier, close to the Hudson Mountains and Mount Murphy respectively, and the simulated responses compared with ice sheet thinning histories derived from glacial-geological data. The sensitivity of results to selected model parameters was also assessed. The model simulations predict a broadly similar response to ocean forcing in both the central and eastern ASE, with an initial rapid phase of thinning followed by a slower phase to the modern configuration. Although there is a mismatch of up to 5,000 years between the timing of simulated and observed thinning, the modeling suggests that the upstream geological records of ice surface elevation change reflect a response to retreat near the grounding line. The model-data mismatch could potentially be improved by accounting for regional variations in mantle viscosity, sea-surface heights and basal sliding properties across the continental global sea-level rise from Antarctica. Glacial-geological studies undertaken in two mountain ranges in the region (the Hudson Mountains and Mt Murphy, close to Pine Island Glacier and Pope Glacier) provided records of the magnitude and speed of those changes. Here we compare those records with model simulations that predict response of the ice sheet to different forcing mechanisms (climate and sea-level) over the past 20,000 years. We find that the model predicts a similar response to forcing at both study sites, with initially fast thinning followed by a slower phase as the ice sheet shrinks to its modern shape and size. The similarity in behavior at both sites may be due to similar ocean forcing mechanisms affecting both glaciers. When we compare the model simulations with geological records of ice sheet change, we find a mismatch of up to 5,000 years. Improving how ice-Earth interactions are represented in the model could potentially resolve this mismatch.

. Thus, improving our understanding of the timing and pace of the last deglacial retreat is critical to providing more robust model projections of future retreat and sea-level rise. Records of changes in surface elevation of the ice sheet through the last deglacial retreat are archived within geological sequences and can be extracted using chronological dating techniques. When obtained from several locations at the present ice sheet margins and inland, such records provide constraints on the varying size and shape of an ice sheet through time, and are thus invaluable for validating and tuning predictive ice sheet models (Albrecht et al., 2020a(Albrecht et al., , 2020bBriggs et al., 2014;Lowry et al., 2020;Pollard et al., 2016). Here, glacial-geological records of ice sheet surface lowering during the past 20,000 years across the Amundsen Sea sector of the West Antarctic Ice Sheet are compared with model simulations of the rate and magnitude of the last deglacial retreat in the same region.

Background
Glacial-geological records of Antarctic ice sheet retreat since the Last Glacial Maximum (LGM) were comprehensively reviewed and synthesized in 2014 (RAISED consortium, 2014). These were deduced using methods such as radiocarbon dating of marine sediments, bathymetric and seismic surveys of submarine glacial bedforms, and cosmogenic nuclide surface exposure dating of terrestrial rock outcrops near modern ice margins. A number of subsequent studies have further improved our knowledge of last deglacial ice sheet surface elevation changes across the continent (e.g., Antarctic Peninsula: Glasser et al., 2014;Jeong et al., 2018; Amundsen-Bellingshausen sector: Johnson et al., 2017;Johnson et al., 2020; Ross Sea sector: Balco et al., 2019;Goehring et al., 2019;Jones et al., 2015;Smellie et al., 2017;Spector et al., 2017Spector et al., , 2019Yokoyama et al., 2016 and others;East Antarctica: Jones et al., 2017;Strub et al., 2015; Weddell Sea sector: Bentley et al., 2017;Balco et al., 2017;Nichols et al., 2019 and references therein). Marine geological data show that grounding line retreat from the outer continental shelf largely occurred prior to the Holocene, between approximately 20 and 10 ka, driven primarily by rising global sea-levels associated with Northern Hemispheric deglaciation. Continent-wide patterns of ice sheet thinning prior to the Holocene are not yet well-known because terrestrial geological data that provide reliable constraints on thinning from 20 to 10 ka are few (see Figure 11, Small et al., 2019). In contrast, several studies have shown that many hundred meters of ice sheet thinning occurred very rapidly at some locations during the early-to mid-Holocene (e.g., Bentley et al., 2011;Johnson, et al., 2019aJohnson, et al., , 2020Johnson, et al., , 2014Jones et al., 2015;Spector et al., 2017). Grounding line retreat to near-modern positions during that time was probably driven by warming of sub-surface ocean water and an associated increase in sub-ice shelf melting (Golledge et al., 2013;Mackintosh et al., 2011;RAISED consortium, 2014).
The Amundsen Sea sector of the West Antarctic Ice Sheet, in particular the Pine Island-Thwaites Glacier catchment, is presently experiencing ongoing rapid thinning and grounding line retreat (Konrad et al., 2018;Milillo et al., 2019;Rignot et al., 2014), which has increased five-fold in the past two decades (Shepherd et al., 2019). Ice sheet models suggest that this retreat may continue and accelerate in coming centuries, driven by Marine Ice Sheet Instability with increasing ice flux across deepening grounding lines, and/or by Marine Ice Cliff Instability with failure of tall ice cliffs (e.g., Arthern & Williams, 2017;Feldmann & Levermann, 2015;Golledge et al., 2019Golledge et al., , 2014. Detailed knowledge of how this sector of the ice sheet varied in both extent and thickness in the last deglacial period is critical for understanding how it can respond to future climate change. A summary of the current state of knowledge is as follows: In the eastern ASE, by the early Holocene (10.6 ka), the grounding line had retreated from a position near the continental shelf edge at the LGM to the inner shelf . Between 12 and 10.6 ka, an ice shelf spanning Pine Island Bay buttressed the grounding line of Pine Island Glacier (Kirshner et al., 2012). High rates of Circumpolar Deep Water (CDW) upwelling onto the continental shelf during the early Holocene may have provided a mechanism for collapse of this ice shelf by 7.5 ka , and may also have been connected with rapid thinning of ice in the Hudson Mountains at ∼8 ka (Johnson et al., 2014). In the central ASE, grounding line retreat occurred with a similar style and timing as in the east. The grounding line had retreated from its maximum LGM extent at the continental shelf edge to near its present location on the inner shelf by 10.1 ka (Smith et al., 2011) and several hundred meters of ice surface lowering occurred shortly after, between 9 and 6 ka (Johnson et al., 2020).

Rationale and Objectives
Although some modeling studies of the last deglacial retreat of the Antarctic ice sheet have examined specific regions, most have been concerned with large-scale behavior (Albrecht et al., 2020a(Albrecht et al., , 2020bBriggs et al., 2013Briggs et al., , 2014Golledge et al., 2014;Gomez et al., 2013Gomez et al., , 2018Lowry et al., 2020;Mackintosh et al., 2011;Maris et al., 2015;Pollard et al., 2016Pollard et al., , 2017Whitehouse et al., 2012). These studies have found broad-scale agreement with the observational pictures described above, but to our knowledge, no last-deglacial studies have focused exclusively on the ASE.
In this study, a nested ice sheet model spanning the ASE is applied to the last deglacial retreat of the Antarctic ice sheet to assess sensitivity to different environmental forcings, to investigate links between grounding line retreat and inland thinning, and to compare the predicted timing and magnitude of retreat in the eastern versus central ASE. The existing body of data described in Section 2 is used as context, but with a focus on comparisons with surface exposure age data from rock outcrops at two sites: the Hudson Mountains adjacent to Pine Island Glacier in the eastern ASE, and Mt Murphy adjacent to Pope Glacier in the central ASE ( Figure 1). Constraints on past ice thicknesses at these sites (Johnson et al., 2014(Johnson et al., , 2020 provide useful limits on the magnitude of inland ice variations.

Compilation of Surface Exposure Age Data
In order to make a comparison between ice sheet thinning histories derived from terrestrial glacial-geological data and model simulations of thinning at the same sites in the ASE over the same last deglacial timescale (20,000 years to present), existing exposure age data were compiled from two key sites in the ASE (Figure 2). Exposure ages are calculated by measuring the abundance of cosmogenic nuclides that are produced JOHNSON ET AL.   (Mouginot et al., 2012;Rignot et al., 2011aRignot et al., , 2017 are overlain on the Landsat Image Mosaic of Antarctica (Bindschadler et al., 2008). Bathymetry is from the International Bathymetric Chart of the Southern Ocean (Arndt et al., 2013). The grounding line (black) is from Rignot et al. (2011b), and the coastline (blue) from the Antarctic Digital Database version 7.1. Dashed yellow lines indicate approximate boundaries between outer (OUT), middle (MID) and inner (IN) shelf regions. by cosmic ray bombardment of an ice-free rock surface (or surface free of other glacial debris or sediments). In simple terms, the higher the measured abundance of the nuclide, the longer the surface has been ice-free. Thus, measurements in rock surfaces along an altitudinal transect act as a "dipstick" for the former height of an ice sheet surface as it thins through time (e.g., Stone et al., 2003). The Hudson Mountains and Mt Murphy exposure age datasets (Johnson et al., 2014(Johnson et al., , 2020 are detailed enough to provide evidence of several hundred meters of ice sheet surface lowering in the ASE during the last deglacial period (Figure 2). A few isolated exposure ages exist from elsewhere in the ASE (Figure 8, Johnson et al., 2020), but are not included in the compilation because they only provide a single deglaciation age at a particular site; several exposure ages would be needed to identify outliers. All exposure ages discussed in the present study were calculated using the online cosmogenic nuclide exposure age calculator, version 3 (Balco et al., 2008; hosted at https:// hess.ess.washington.edu). For 10 Be, the calculations use the default global calibration data set, reference sea-level-high-latitude production rate (3.92 ± 0.31 at.g −1 .yr −1 ), and the "LSDn" scaling scheme (Lifton et al., 2014), which yields smaller systematic biases for Antarctic sites than other commonly used scaling models. For in-situ 14 C age calculations, the production rate is estimated from measurements of the CRO-NUS-A quartz inter-comparison standard (see Johnson et al., 2020), while all other variables are the same as for 10 Be. For details of how the best-fit thinning trajectories were determined, see Supporting Information.
JOHNSON ET AL.

Ice Sheet Modelling
An ice sheet model (PSUICE3D) was applied to the last deglacial retreat of the Antarctic ice sheet. Model simulations are compared with trajectories for ice sheet surface lowering derived from the exposure age dating. The model is described further in the Supporting Information, and detailed formulations are given in Pollard & DeConto (2012a) and Pollard et al. (2015). This section describes aspects specific to this paper.
The model uses a polar stereographic grid centred on the South Pole, over a limited-area nested domain spanning the ASE to allow relatively fine horizontal resolution of 5 km. Model maps below show the entire nested domain. It is run from 30 ka to the present, with lateral boundary conditions obtained from a previous long-term continental-scale simulation at coarser resolution. The model is initialized at 30 ka from the continental-scale run to allow for ∼10,000 years of spin up; results are shown here just for the last 20,000 years of the simulations. The continental-scale run (unpublished) is similar to that in Pollard & De-Conto (2009), except run from 125 to 0 ka at 20 km resolution with the current model version; it simulates LGM expansion and deglaciation to modern reasonably well, as in that paper. The long spin-up time, and the large buffer zones between the nested domain edges and the central region of interest, mean that the details of the continental run have very little impact on the results here.
Past atmospheric forcing is prescribed simply as uniform perturbations to modern climatology, as in previous applications (Pollard & DeConto, 2012a;Pollard et al., 2016Pollard et al., , 2017. Air temperatures are obtained from a modern climatological Antarctic data set (Le Brocq et al., 2010), plus an imposed sinusoidal seasonal cycle, and spatially uniform past perturbations proportional to a deep-sea core δ 18 O record (Lisiecki & Raymo, 2005) and orbital insolation variations. Lapse-rate corrections are applied to account for the difference between model and climate data set surface elevations. The observed annual mean precipitation (LeBrocq et al., 2010) is similarly corrected following temperature changes using a Clausius Clapeyron-like relation. The resulting temperature and precipitation variations over the last deglacial period agree well with those at the WAIS Divide core site (Cuffey et al., 2016;WAIS Divide Project Members, 2013;shown in Johnson et al., 2020, Figure 9d), both in timing and amplitude, and also with climate model output at that location (Lowry et al., 2019a, Figures 2 and 3).
Past oceanic forcing is prescribed from an archived atmospheric-oceanic global climate model (A/OGCM) simulation of the last 22 kyr with oceanic resolution of 3.6° longitude and ∼1.5° latitude in high latitudes (Liu et al., 2009), using their water temperatures at 400 m depth corresponding roughly with CDW, spatially interpolated from their coarse global grid and extrapolated to the central ASE and under floating ice shelves (see Johnson et al., 2020, Figure 9b). A warming pulse occurs in the data set at 14 ka, with ∼0.2°C running-mean amplitude and several hundred years duration. It corresponds in time with an imposed event in the Liu et al. (2009) experiment related to Meltwater Pulse 1A, and may be model dependent; it is removed in our forcing, but has negligible effect on our results. Global mean sea-level variations, driven mainly by Northern Hemispheric ice-sheet melting, are prescribed from the eustatic sea-level component of the ICE-5G data set (Peltier, 2004).
For the present study, a set of four simulations was used to show the basic range of model behavior over wide but reasonable variations of two uncertain model parameters: sub-ice oceanic melt coefficient (O; an overall factor multiplying the oceanic melt rate at the base of floating ice with nominal value 1 corresponding to K in Pollard & DeConto, et al., 2012a;Pollard et al., 2015) and e-folding time (τ; yr) of asthenospheric relaxation in the Elastic Lithosphere Relaxing Asthenosphere (ELRA) bed model (Pollard & DeConto, 2012a), as follows: These two parameters and their ranges were chosen on the basis of previous continental-scale studies with much larger ensembles (Pollard et al., , 2017; running such large ensembles for the ASE would be beyond the scope of this paper. The two parameters provide a preliminary exploration of the envelope of mod-el behavior. Five additional simulations extending this exploration to other model parameters are described primarily in Supporting Information, but are also included in the discussion below: Simulation F: O = 10, τ = 10 Simulation G: O = 10, τ = 100, alternate sea-level forcing (see Supporting Information) Simulation H: O = 10, τ = 100, basal sliding coefficient on continental shelf C shelf = 10 −7 (changed from 10 −5 ) m yr −1 Pa −2 Simulation H': As H, except C shelf = 10 −6 m yr −1 Pa −2 Figure 3 shows maps of ice distribution for the two model simulations with the most and least ice retreat overall (B and C). Figure 3a shows maps of ice distribution at selected times for model simulation B, which yields the most realistic modern ice configuration of the four simulations tested, with D a close second (judged mainly by the differences from modern shown in Figure 4). In this simulation, there is rapid and widespread retreat of grounded ice over the outer and middle continental shelf in the central and eastern ASE between 14 and 13 ka. After 13 ka, grounding lines recede more gradually to their modern positions, with large ice shelves persisting over the inner-and mid-shelf regions until after 10 ka. This behavior and associated feedbacks are discussed in Section 5.2. Simulation C ( Figure 3b) incorporates less ocean melt and faster bed rebound. Grounding line retreat is slower and later than in simulation B, with grounding lines still located in the mid-shelf region of the ASE at ∼12 to 10 ka. This is in better general agreement with marine geological and geophysical records . However, at the end of simulation C, modern grounded and floating ice extents are considerably greater than observed, especially in the eastern ASE. Figure 4 shows the modern state at the end of all four simulations. The most influential parameter for the modern state is the ocean melt coefficient (O); the simulations with larger values (B and D, with O = 10) produce relatively realistic modern ice distributions, whereas those with smaller values (A and C, with O = 1) produce modern ice extents that are unrealistically large. This is because, in all simulations, ocean melting becomes important later in the runs, during the Holocene (Section 5.2). Bed rebound is most important in the earlier pre-Holocene stages, so the influence of the asthenospheric relaxation time scale (τ) on the modern states in Figure 4 is small. Figure 5 shows modeled ice thickness variations at the Hudson Mountains and Mt Murphy for all four simulations. Changes relative to present-day model thickness are shown, since these correspond directly to the changes recorded by the cosmogenic exposure age data. The model predicts a similar response to forcing of both Pine Island Glacier and Pope Glacier in all simulations, with an initially fast rate of ice surface lowering prior to the Holocene (i.e., before ∼12 kyr BP), followed by a slower phase with only minor thinning through the Holocene to the present-day. In all simulations, most of the thinning occurs between 15 and 10 ka. This is much earlier than indicated by the glacial-geological data (Section 5.4).

Rapid Abrupt Episodes of Thinning Prior to the Holocene
The two simulations with large asthenospheric relaxation times (A and B) show rapid and abrupt thinning episodes at both sites simultaneously (∼13-12 ka for A, ∼14.5-13.5 ka for B; Figure 5). These episodes are also evident in other figures, including the map views for simulation B in Figure 3a. The prescribed sea-level rise, predominantly due to Northern Hemispheric ice melt occurring between ∼15 and 10 ka (Peltier, 2004), may trigger these episodes. When sea-level rises sufficiently, it causes grounding lines to retreat rapidly across the relatively flat seafloor topography of the outer ASE continental shelf to submarine bedrock ridges of the mid shelf, which pin and stabilize the grounding line (Figure 7). This behavior is consistent with the widely recognized repeated pinning and unpinning of grounding lines between bedrock ridges during retreat, found in other modeling studies (e.g., Parizek et al., 2013;Seroussi et al., 2017), and is also suggested by the sudden jumps in grounding line positions in Figure 6, below.
The major thinning episodes are more sudden in simulations A and B than in C and D ( Figure 5). This is because the slower bedrock rebound (τ = 3,000 years in simulations A and B, compared with τ = 100 years in C and D) under the retreating ice allows bedrock elevations at grounding lines to remain deeper for longer. This produces larger ice fluxes across the grounding lines (which depend strongly on local bedrock depth; Schoof, 2007; see Supporting Information), and consequently faster drawdown of interior ice. In simulations C and D with faster bedrock rebound, there is still widespread grounding line retreat and thinning, but it occurs more gradually, slowed by shallower grounding line depths and less ice flux from upstream (cf. Gomez et al., 2013Gomez et al., , 2015. Nevertheless, most of the inland thinning at the study sites still occurs before   Fretwell et al., 2013) and ice speeds (Rignot et al., 2011a), with both datasets aggregated to the model's 5 km grid. The grounding line is shown by a thick black line, and the two field study sites are marked by red stars. For modern ice distributions, the contour scales for grounded (blue) versus floating (pink) ice can be distinguished as in Figure 3.

Holocene Thinning
After these major retreat episodes, the much slower retreat and ice thinning in the Holocene (post 10 ka) is likely forced by a combination of ocean warming and increased sub-ice-shelf melting, and the remaining minor amount of sea-level rise. Ocean temperatures in the model are prescribed using the archived A/ OGCM (400 m) results of Liu et al. (2009), in which nearly all of the warming in the proximal ASE ocean occurs between 10 and 5 ka, with very little before 11 ka (see Johnson et al., 2020, Figure 9b). After 10 ka, global mean sea-level rises by only ∼20 m (compared to a ∼95 m rise between 20 and 10 ka; ICE-5G, Peltier, 2004), but this may play a role in helping to unpin grounding lines from the closely spaced ridges in the inner ASE (see Figure 7, below).

Ice Sheet Response Along Flowlines
The model results in Figure 6 show that the timing of inland thinning of grounded ice at both study sites (Hudson Mts and Mt Murphy) coincides with that of grounding line retreat, for both the large and relatively  (b), the x-axis shows distance downstream (km) from the start of the flowline. The y-axis shows elevation (m) relative to modern sea-level. The ocean surface is not shown, but would be horizontal lines extending to the right from the ice edge, at elevations given in the figure legend (m, relative to modern). The approximate positions of closest approach to the exposure age data sites are shown by vertical gray bars. Vertical dashed lines between labels IS, MS, and OS (inner, mid, and outer shelves, respectively) indicate the locations where the transects cross latitudes corresponding to the inner-to mid-shelf transition, mid-to outer-shelf transition, and outer-shelf edge. (c) Transect paths (brown lines, labeled A-B and A′-B′), superimposed on a map of modern observed ice speeds (Rignot et al., 2011a) aggregated to the model's 5 km grid, with the two field sites marked by red stars. abrupt retreats before 10 ka, and for the more gradual later Holocene retreats. Figure 7 provides another view of the retreat along flowlines, showing model transects for selected timeslices through the deglacial period along Pine Island Glacier and Pope Glacier, and extended to the outer continental shelf. Results for simulation D are shown, which along with B produces quite realistic modern configurations (Figure 4), but the later retreat is in less conflict with data ( Figures 5,6 and 8). Each transect passes close to the respective field sites. For both transects, the greatest extent of grounding line retreat, from the outer-to mid-continental shelf, occurs before 10 ka. This is consistent with the retreat histories in Figures 5-6; much of this retreat occurs in a single jump, especially for the Pope Glacier transect. As discusse for Figure 5, we consider this to be forced by sea-level rise causing retreat back to pinning points on the mid-shelf. Subsequent grounding line retreat from 10 ka to present is less extensive and more gradual overall, caused by ocean warming increasing sub-ice melting and thinning of floating ice shelves, and augmented slightly by the remaining sea-level rise (as discussed for Figure 5). This later retreat also tends to occur in steps, reflecting pinning on the succession of bathymetric ridges on the inner shelf of the ASE (Figure 7).
The concurrence in the model simulations of inland thinning of grounded ice and grounding line retreat suggests that there is a correlation between retreat and thinning at the data sites themselves. Similarly, within a given time interval, the amount of upstream thinning at the data sites closely correlates with the JOHNSON ET AL.   (Johnson et al., 2014(Johnson et al., , 2020 are shown as black dots for comparison, and, in contrast to the data in Figure 2, are plotted relative to the present-day ice surface for comparison with the model outputs (i.e., elevations are height above the modern ice surface; ice surface heights used for Mt Murphy sample elevations are those in Figure 5, Johnson et al., 2020). Dashed lines show the best-fit thinning trajectories to the data, as in Figure 2 but for ice thickness changes, instead of altitudes. The labeled horizontal dashed lines show the amount of ice thickening, relative to present, required to submerge the highest currently exposed peak or ice-free location visited. distance over which the grounding line has retreated. This correlation has a theoretical basis in studies of waves of adjustment propagating inland from marginal perturbations, with propagation times of less than ∼500 years over the first ∼500 km for an East Antarctic flowline (Figure 2, Alley & Whillans, 1984). This result is important because exposure age data from upstream/inland sites are frequently utilized to infer changes hundreds of kilometers downstream without necessarily having evidence that upstream ice surface changes reflect changes in grounding line position. However, whilst the relationship between timing and distance of grounding line retreat and ice sheet thinning appears to hold along most transects shown in Figure 7, where topographic highs occur (such as in Figure 7b), it may be interrupted upstream of those peaks, resulting in the ice sheet response to downstream changes not being transmitted as far inland. Thus, while exposure age data from the upstream nunataks provide reliable information on the downstream timing of grounding line retreat in this study, the relationship between retreat and thinning modeled here may not be the same under different bathymetric or other boundary conditions elsewhere.

Model-Data Comparison
When compared with the ice sheet thinning history derived from the surface exposure age data, the model outputs imply a potential mismatch with some of the geological data, especially at Mt Murphy. The most obvious differences between the model simulations and the ASE cosmogenic exposure age data involves the timing of the retreat (see Figure 8, where the exposure age data are plotted together with the modeled thinning at the two data sites).
In the model, ice thinning at both sites mainly occurs prior to the Holocene, between ∼15 and 10 ka. This is in reasonable agreement with the exposure age data from lower elevations in the Hudson Mountains (i.e., below ∼90 m, Figure 8c); at higher elevations at that site, thinning appears to have been faster than in the model, but model-data comparison is limited both by the lack of exposure ages above 370 m asl and any geological constraints on the timing of thinning prior to the early Holocene ( Figure 2). In contrast with the Hudson Mountains, exposure ages from Mt Murphy do not form a well-defined trend, with samples at the same elevation often yielding ages that differ by several thousand years. This scatter is most likely due to exposure of upstream bedrock surfaces (the source of erratic cobbles) for long periods prior to glacial erosion, resulting in the nuclide inventory accumulated during previous ice-free periods remaining at the time of erratic emplacement; such "nuclide inheritance" is common in Antarctica (e.g., Balco, 2011). Thus, the choice of ages to include in the regression analysis is left to expert judgment. Here, the published regression line for a subset of samples from Mt Murphy (Johnson et al., 2020;Supporting Information) is used as the best estimate for the timing and pace of thinning. Using different samples could lower the gradient of the thinning trend shown in Figures 2 and 8b, permitting a better fit to the model output, but this would not be geologically reasonable because the potential for nuclide inheritance at this site means that older ages are unlikely to represent the timing of deglaciation. Model-data comparison using the published regression line implies that most thinning at Mt Murphy occurred up to 5,000 years later than in the model, between ∼10 and 5 ka (Figure 8b).
A number of additional factors might contribute to model-data mismatch, as follows (some are explored further in Supporting Information by additional simulations): (1) The exposure ages may not be correct.
Exposure age calculations rely on knowing how the production rate of cosmogenic nuclides in rock surfaces varies with time, location and elevation. To apply any production rate to a particular sample site, it is necessary to scale it for the altitude and latitude of the sample site. Several production rate scaling-schemes exist for 10 Be and in-situ 14 C (the nuclides from which the exposure age data in this study were determined) (Balco et al., 2008). To mitigate the possibility of model-data mismatch resulting from incorrect exposure age calculations, we used the scaling-scheme that is generally considered the most appropriate for Antarctic sites (LSDn; Lifton et al., 2014). When other published scaling-schemes and production rates are used, the resulting exposure ages differ by no more than ∼1.4 ka, which is within the external errors of the ages. Thus, we do not consider the choice of production rate/scaling scheme to be a plausible explanation for model-data mismatch.
Other factors, such as snow cover or cover by glacial till, may cause apparent exposure ages to be too young. For Mt Murphy and the Hudson Mountains, the impact of snow cover on exposure age calculations was assessed in the original publications (Johnson et al., 2020 and2014, respectively). For samples with uncorrected exposure ages in the 0-16 ka range, cover by 100 cm of snow with density of 0.35 g cm −3 for 8 months of every year (a best estimate based on the authors' experience of snow accumulation at these sites during the austral summer and assuming greater, more persistent, accumulation during the winter) makes the ages no more than 1,900 years older. Alternatively, if a covering of till of 10 cm thickness with density of 1.8 g cm −3 is assumed, the difference is less than 2,000 years. Neither of these alone is sufficient to explain a model-data mismatch of several thousand years. We consider a situation where both occurred simultaneously-which could bring the ages close to the model thinning trajectory-extremely unlikely because no continuous accumulations of till have been observed in the area (the glacial deposits are scattered cobbles perched on scoured bedrock surfaces).
Glacial isostatic adjustment (GIA) can also have an impact on exposure age calculations because it causes the elevation of a sample site, relative to sea-level, to change through time. This in turn affects cosmogenic nuclide production. Using the online tool "iceTEA" (http://ice-tea.org/en/tools/correct-elevation-change/; Jones et al., 2019a), a GIA correction was applied to the exposure age data to assess the potential impact of post-exposure elevation change on the age calculations (see Supporting Information). When corrected for the effect of elevation change predicted by three GIA models (ICE-5G, Peltier, 2004;ICE-6G, Peltier et al., 2015;W12, Whitehouse et al., 2012), the exposure ages in the range 0-10 ka (i.e., when the last deglaciation occurred; Figure 2) were no more than 600 years older, and for 5-16 ka (the range of the majority of exposure ages from this region; Figure 2) they were less than 2,000 years older. This difference is not sufficient to explain a >2,000-year model-data mismatch. The combined effect of snow cover and GIA would make the exposure ages up to 3,900 years older (although without data on depth and persistence of winter snow accumulation at the sample sites, it is difficult to assess the likelihood of such a shift). This could explain a model-data mismatch of more than a couple of thousand years, but it would not be sufficient to explain a mismatch of as much as 5,000 years (e.g., Figure 8b).
(2) Ice thicknesses at the data sites could be influenced by local small-scale topography or other micro-climatic factors.
An example is described in Johnson et al. (2020), where the presence of a prominent windscoop near one of the study sites at Mt Murphy (Turtle Rock) appears to have resulted in deglaciation of that site approximately 4,000 years earlier than the timing of regional deglaciation. However, comparable features are absent from other Mt Murphy sites, and also from the Hudson Mountains. Conversely, the model results at the grid cell closest to each site could be anomalous and not representative of the other nearby cells; however, no such pronounced cell-to-cell variations are seen in the model. Thus small-scale local topographic and micro-climatic variations are not likely explanations for model-data mismatch.
(3) The prescribed climatological forcing and its timing may be fundamentally unrealistic.
Since eustatic sea-level over the last 20,000 years is relatively well known, with most sea-level rise occurring before 10 ka (Peltier, 2004), this component of the model forcing is probably realistic, assuming regional sea-level changes in the ASE are close to the global mean. Small changes to the sea-level forcing do not significantly change the overall results (see simulation G in Figure S2, Supporting Information). The atmospheric forcing (Section 4.2) is also sufficiently well-known that its uncertainties would only have small effects on the ice-sheet changes addressed here (cf. Lowry et al., 2019b). However, the prescribed oceanic temperatures from Liu et al. (2009) are more questionable. To our knowledge, Liu et al. (2009) andLowry et al. (2019a) are the only applicable atmosphere-ocean model integrations through the last deglaciation to date, and these coarse-grid global ocean models may not capture important small-scale variations of CDW upwelling onto the ASE continental shelf. CDW influence is thought to have increased significantly after 10.4 ka until 7.5 ka ; this increase is in rough agreement with the proximal oceanic temperature trends in Liu et al. (2009) (see Figure 9b, Johnson et al., 2020) but not with those further offshore in Lowry et al. (2019a) (see their Figures 7 and 8). In order to improve how small-scale changes in CDW upwelling in the forcing are captured in the present study, a record of Holocene ocean temperature changes in the ASE is needed. This is the subject of ongoing research, but results are not yet available. However, simulation E-with no oceanic melting at all-produces little overall change in the model ice thinning and retreat, which is still dominated by sea-level rise at ∼15 to 10 ka ( Figure S2, Supporting Information). In summary, the data used for ocean temperatures may not adequately capture small-scale variations in CDW influence in the ASE region, but this is probably not the main cause of the model-data mismatch.
In the model, grounding line retreat and upstream thinning between ∼15 and 10 ka are forced mainly by sea-level rise, and there is little subsequent response to the later ocean warming between ∼10 to 5 ka. It is reasonable to question the model sensitivity to these forcings, since the response to sea-level in the model is strongly controlled by an uncertain parameterization of ice flux across the grounding line (Schoof, 2007), and its treatment of sub-ice-shelf oceanic melting is relatively simple. In larger-scale Antarctic simulations, fair agreement with data has been found with this and other ice sheet models (Golledge et al., 2014;Pollard & DeConto, 2009;Pollard et al., 2016), but scale-dependent shortcomings for smaller regions are possible.
For instance, the horizontal grid size in the model, 5 km, may not adequately resolve the bathymetry of the Amundsen Sea shelf. Previous sensitivity modeling studies in this region have found that grounding line retreat may be influenced significantly by fine-scale ridges in basal topography, with intervals of pinning on each ridge alternating with rapid retreat between ridges. These require horizontal resolutions of several hundreds of meters to ∼1 km for accurate simulations (Durand et al., 2011;Parizek et al., 2013;Seroussi et al., 2017), which are not feasible here. However, these studies focused on modern grounding line locations and retreat upstream from them, and their results may not necessarily apply to early grounding line retreat over the mid-to outer-continental shelf of the ASE where the bedrock topography is less rugged, with smooth sedimentary strata and fewer potential pinning points than the predominantly crystalline inner shelf region (Gohl et al., 2013;Graham et al., 2009).
(5) The model does not account for spatial variations in sea-surface height.
In the model used here, sea-level is simply prescribed as a spatially uniform elevation, following the eustatic sea-level component of the ICE-5G data set (Peltier, 2004). This does not account for sea-surface geoid distortions due to ice-ocean-solid Earth gravitational interactions or Earth rotational effects. Coupling an ice sheet model with a global Earth-sea-level model is needed to capture these effects; for instance, as an ice mass retreats, the reduced gravitational attraction causes shoaling of the nearby ocean, shallower grounding line depths, reduced ice fluxes, and slower ice retreat in a negative feedback (Gomez et al., 2013(Gomez et al., , 2015(Gomez et al., , 2018. In ensembles of deglacial Antarctic simulations coupled with a global Earthsea-level model, Pollard et al. (2017) found that results agreed reasonably closely with those from the uncoupled model, suggesting use of the latter may be adequate for continent-scale experiments. On the other hand, it may be inadequate at the finer-scales of the present study: sea-surface heights in the ASE will have departed from the global mean during deglaciation, and may have differed spatially between the eastern and western sides of the basin, with potentially strong effects on local ice retreat that can only be captured by a high-resolution coupled model. Therefore, it is possible that improvements to representation of sea-surface height in the model may help to resolve the model-data mismatch.
(6) The bed model does not include the effects of spatial variations in Earth rheology.
The locally relaxing asthenosphere and elastic lithospheric plate of the ELRA model provides a rough approximation to the more realistic Earth deformation captured by global Earth models, and as mentioned above, yields results in other studies that agree reasonably well with those using global Earth models. However, recent geophysical evidence identifies very weak low-viscosity zones in the mantle under parts of West Antarctica, including the ASE Barletta et al., 2018;Heeszel et al., 2016, O'Donnell et al., 2017. The fast bedrock rebound associated with very weak mantle viscosity can have a strong negative feedback on deglacial ice retreat: as the grounding line retreats, rapid bedrock rebound shallows local ocean depths, causing thinner ice at the grounding line, reduced ice flux from upstream, and slower drawdown of interior ice (Gomez et al., 2013(Gomez et al., , 2015. The short isostatic relaxation timescale in simulations C and D (τ = 100 years; Figure 8), or even the extreme value of τ = 10 years used in simulation F ( Figure S2, Supporting Information), may be representative of the weakest mantle viscosities suggested by the recent data, but this parameter is applied uniformly over the model domain, and does not capture the observed spatial variations. It is possible that regional variations in Earth properties could have larger effects on the timing of local ice retreat within the ASE, especially when combined with regional variations in sea-surface heights (discussed above), which together control ocean depths at grounding lines. Again, capturing these effects will require coupling with regional Earth-sea-level models that include lateral variations in Earth properties (Gomez et al., 2018).
(7) The modeled LGM ice sheet configuration influences subsequent retreat.
A recent study by Albrecht et al. (2020aAlbrecht et al. ( , 2020b compared an ensemble of continental Antarctic simulations over the last two glacial cycles, generated using the PISM ice sheet model, with various geological data in order to better constrain model parameter ranges. In contrast to the results of the present study, Albrecht et al. predict that most grounding line retreat occurred after 10 ka, including in the ASE sector, which concurs with the Holocene timing of thinning implied by the ASE exposure age data described above. This later timing of modeled grounding line retreat may be related to their greater LGM ice thicknesses predicted around coastal areas (Albrecht et al., 2020b). In contrast, large-ensemble studies with the present model (Pollard et al., , 2017 found that relatively slippery bed properties are appropriate on continental shelves. This produces relatively thin ice in coastal areas during the LGM, making those regions sensitive to sea-level forcing and susceptible to early retreat. Better constraints on LGM ice thicknesses would help to resolve model discrepancies and refine model parameters (see Sections 5.5.3 and 6, below).
To examine the potential for thicker LGM ice to resolve the model-data mismatch in timing of ice thinning, two additional simulations H and H′ were performed with reduced sliding coefficients on the continental shelf (see Supporting Information). In these, the main model thinning is significantly delayed at both data sites, and simulation H improves the fit to the exposure age data at Mt Murphy ( Figure S2). However, modeled modern ice thicknesses are unrealistic in these simulations, with far too thick ice at the data sites; also, initial grounding line retreat across the outer ASE occurs much later than indicated by geological data (see Text S2, Figures S3 and S4, Supporting Information). Nevertheless, the effects of the two simulations are different at the two data sites ( Figure S2), which suggests that simulations with spatially varying sliding coefficients could could help to resolve the model-data mismatch.
In summary, we conclude that the most likely contributing factors to the model-data mismatch are that (a) the model does not capture regional variations in sea-surface heights, or strong variations in mantle viscosity within the ASE, and (b) a uniform basal sliding coefficient is prescribed for the whole continental shelf, missing potentially important spatial variability.

Mechanism for Ice Sheet Thinning in the Central ASE
Upwelling of warm CDW onto the inner shelf of Pine Island Bay increased during the early Holocene (10.4-7.5 ka; Hillenbrand et al., 2017) and may have triggered grounding line retreat and collapse of an ice shelf in the eastern ASE, resulting in rapid ice surface lowering in the Hudson Mountains (Johnson et al., 2014). As yet, however, there are no marine geological data available from the inner shelf in front of Pope Glacier (i.e., immediately north of the Crosson Ice Shelf; Figure 1) that could show whether a similar effect was experienced there. Nevertheless, the existing exposure age datasets from Mt Murphy and the Hudson Mountains imply that several hundred meters of ice sheet thinning occurred during the early-to mid-Holocene (∼9-5 ka) in both the eastern and central ASE (Figure 2). If CDW upwelling onto the shelf occurred across a broader area of the ASE, that is, not just in Pine Island Bay, it could have caused grounding line retreat and subsequent thinning at both sites simultaneously. We therefore speculate that the similar model responses to forcing at both sites ( Figure 5) and concurrent timing of thinning determined from exposure age dating ( Figure 2) point to early Holocene loss of an ice shelf in the central ASE as well as in the eastern ASE.

Late Holocene Record of Ice Sheet Thinning
No glacial deposits or bedrock surfaces across the whole ASE have yet yielded reliable deglaciation ages younger than 5 ka (Johnson et al., 2020). The model simulations undertaken here all predict that ice surface lowering and associated grounding line retreat were essentially complete by 5 ka in both the eastern and central ASE (Figure 3 and Figures 5-8). If correct, this implies that the ice sheet surface has been stable throughout the past 5,000 years in the Hudson Mountains and at Mt Murphy, concurring with the lack of geological evidence for recent thinning, but contrasting with dramatic rates of grounding line retreat and flow acceleration that have occurred on both Pine Island and Pope Glaciers during the past few decades (Konrad et al., 2018;Rignot et al., 2014;Shepherd et al., 2019). Alternatively, if the ice sheet has not been stable for the past 5,000 years, but instead experienced at least one episode of thinning below present and re-thickening during that time, the glacial-geological record of this would lie below the modern ice sheet surface. Such cycles have been identified in recent Antarctic-wide modeling studies (cf. Albrecht et al., 2020b;Kingslake et al., 2018;Lowry et al., 2020) and could explain why no glacial-geological record of late Holocene ice surface lowering has yet been found above the modern ice sheet surface anywhere in the ASE (Johnson et al., 2020). Glacial-geological evidence to test this hypothesis can only be obtained by subglacial bedrock recovery drilling (Spector et al., 2018).

Maximum LGM Ice Thickness
Reliable estimates of the maximum Antarctic ice sheet thickness at the LGM are essential for testing the validity of ice sheet models. However, thus far, no samples collected from the ASE have yielded surface exposure ages older than the LGM (none are >25 ka; most are <16 ka). Consequently, nothing is known about the maximum LGM ice sheet thickness in this sector of the Antarctic ice sheet. All model simulations except H and H′ imply that the Hudson Mountains would have been just submerged or close to the ice surface at the LGM (Figure 8 and Figure S2), with the two simulations (B and D) that yield more realistic modern conditions implying submergence by ∼100 m. In contrast, all simulations imply that the highest accessible exposed rock surface at Mt Murphy would have been at least ∼400 m above the ice surface through the LGM. Identification of the precise location of the LGM limit on both peaks, and sampling for surface exposure dating along a transect of all ice-free surfaces above that, would help to resolve current discrepancies in model simulations of the LGM ice sheet configuration (see Section 5.4 (7)).

Conclusions and Further Work
A limited-area nested ice sheet model (Pollard & DeConto, 2012a) was applied to the last deglacial retreat of the West Antarctic Ice Sheet in the ASE, at 5 km resolution. The ice sheet response to ocean forcing was examined at the Hudson Mountains, Mt Murphy, and along the flowlines of the adjacent Pine Island and Pope Glaciers. By comparing model outputs with thinning histories derived from cosmogenic nuclide exposure dating, we draw the following conclusions: • The model simulations suggest that the glacial-geological records of ice surface lowering at both sites reflect the pace and magnitude of downstream grounding line retreat. • The model simulations predict a broadly similar ice sheet response to forcing at both sites, with an initial rapid phase of ice surface lowering followed by a slower phase to the modern configuration. Although the timing of modeled thinning prior to the Holocene contrasts with the early-to mid-Holocene ice surface lowering suggested by the exposure age data, the similar response at both sites is consistent with a hypothesis that both Pine Island and Pope Glaciers were responding to similar (ocean) forcing mechanisms. • The model simulations predict several hundred meters of surface lowering at both sites. This is consistent with the exposure age data, which suggest at least 140 and 560 m of ice surface lowering during the early-to mid-Holocene in the Hudson Mountains and at Mt Murphy, respectively. • There is a mismatch between the timing of ice surface lowering recorded by the exposure age data at both sites and the timing of thinning predicted by the model: depending on how assumptions are applied in the exposure age calculations (relating to snow cover and GIA, for example), the initial rapid phase of thinning occurs up to ∼5,000 years too early in the model. With the model in its present form, this discrepancy is not resolved by changing the role of ocean melting, faster bedrock rebound, or by regional sea-level variations (see Supporting Information).
Further research would be valuable for improving both the ice sheet model and our knowledge of the pace of the last deglacial retreat in the ASE region. Although our experiments suggest that the role of climate forcing is secondary to sea-level forcing (Section 5.4 (3)), the simple prescribed atmospheric and oceanic forcing used here may not capture regional variations over the whole ASE region (such as small-scale changes in CDW influence on grounding line retreat), and should therefore be improved. Obtaining a data set of measured, rather than modeled, ocean temperatures for the Holocene in the ASE would improve this situation. High-resolution regional atmospheric and oceanic climate models would be an alternative, but long-term simulations over the last ∼20,000 years are not available to our knowledge. The model-data mismatch described here could be improved by better representations of Earth-ice interactions which could potentially change the timing of the model response to forcing. For example, global Earth-sea-level models incorporating laterally varying Earth properties (Powell et al., 2020) could be coupled to the ice sheet model to capture regional sea-surface geoid variations as well as the effects of spatially variable weak mantle viscosities below the ASE (Gomez et al., 2018). Information on spatial constraints in basal sliding properties on the continental shelf would also be valuable, and would have the potential to produce distinct model responses in different locations across the ASE. Lastly, improved glacial-geological constraints on LGM ice thicknesses in the ASE coastal region would help to resolve inter-model differences (cf. Albrecht et al., 2020aAlbrecht et al., , 2020b and further constrain basal sliding coefficients. In terms of improving the glacial-geological data set, rock samples from a broader elevation range in the Hudson Mountains would help determine whether the peaks were submerged by ice at the LGM. A subglacial bedrock recovery drilling campaign in the Hudson Mountains forms part of the International Thwaites Glacier Collaboration (2018-2023), and will investigate late Holocene history of grounding line change in the eastern ASE. Finally, exposure dating of glacially eroded surfaces from the highest accessible ridge of Mt Murphy would provide a valuable constraint on maximum LGM thickness of the ice sheet in the ASE.