Did Holocene climate changes drive West Antarctic grounding line retreat and re-advance?

Knowledge of past ice sheet configurations is useful for informing projections of future ice sheet dynamics and for calibrating ice sheet models. The topology of grounding line retreat in the Ross Sea Sector of Antarctica has been much debated, but it has generally been assumed that the modern ice sheet is as small as it has been for more than 100,000 years (Conway et al., 1999; Lee et al., 2017; Lowry et al., 2019; McKay et al., 2016; Scherer et al., 1998). Recent findings suggest 15 that the West Antarctic Ice Sheet (WAIS) grounding line retreated beyond its current location earlier in the Holocene and subsequently re-advanced to reach its modern position (Bradley et al., 2015; Kingslake et al., 2018). Here, we further constrain the post-LGM grounding line retreat and re-advance in the Ross Sea Sector using a two-phase model of radiocarbon input and decay in subglacial sediments from six sub-ice sampling locations. In addition, we reinterpret high basal temperature gradients, measured previously at three sites in this region (Engelhardt, 2004), which we explain as 20 resulting from recent ice shelf re-grounding accompanying grounding line re-advance. At one location – Subglacial Lake Whillans (SLW) – for which a sediment porewater chemistry profile is known, we estimate the grounding line re-advance by simulating ionic diffusion. Collectively, our analyses indicate that the grounding line retreated over SLW ca. 4000 years ago, and over sites on Whillans Ice Stream (WIS), Kamb Ice Stream (KIS), and Bindschadler Ice Stream (BIS) ca. 4500, ca. 2000, and ca. 2000 years ago respectively. The grounding line only recently re-advanced back over those sites ca. 1000, ca. 25 1100, ca. 500, and ca. 500 years ago for SLW, WIS, KIS, and BIS respectively. The timing of grounding line retreat coincided with a warm period in the midto late-Holocene. Conversely, grounding line re-advance is coincident with climate cooling in the last 1000-2000 years. Our estimates for the timing of grounding line retreat and re-advance are also consistent with relatively low carbon-to-nitrogen ratios measured in our subglacial sediment samples (suggesting a marine source of organic matter) and with the lack of grounding-zone wedges in front of modern grounding lines. Based on these 30 results, we propose that the Siple Coast grounding line motions in the midto late-Holocene were driven by relatively https://doi.org/10.5194/tc-2020-308 Preprint. Discussion started: 19 November 2020 c © Author(s) 2020. CC BY 4.0 License.


Radiocarbon Model Methods
We develop a two-phase model of 14 C and 12 C to examine the timing of grounding line retreat and re-advance for our field sites. We calculate the evolution of 14 C and 12 C separately. The first phase of the model represents the time after 15 the grounding line retreated beyond our sites, and the second represents the time after the grounding line had re-advanced.
To model the concentration of 14 C, n, in the first phase of our model (exposure to seawater), we assume that radiocarbon is being added to the sediments at a constant rate, a, while a fraction of this unstable radioisotope decays: Where l is the decay constant and t is time. To simplify our calculations, we substitute:

25
Where t is the mean lifetime of 14 C (8033 years). Eq. (S2) thus simplifies to: We integrate Eq. S3 using the integrating factor method and multiply both sides of the equation by the factor /0# : /0# !" !# − /0# = /0# (S4) Remembering the product rule, we observe that this operation transforms the left-hand side of Eq. (S4) into a derivative of the product of n and the integrating factor, such that: 35 We obtain a definite integral of both sides of Eq. (S5) using a dummy variable of integration, x.

40
∫ 8 /09 ( ); Which results in: We assume that no = 0 because the sediments have been isolated for long enough that there should be no significant radiocarbon present at the moment the grounding line retreated past our sediment sampling sites, exposing them to accumulation of radiocarbon at rate a. Recognizing this, solving for n(t), and substituting for k using Eq. (S2) we obtain the final equation for changes in 14 C concentration in phase 1 of the model: 50 To model the evolution of the most abundant and stable isotope of carbon, 12 C, we assume that it was also being added at a constant rate, A, during phase one when the sediments at the sampling sites were exposed to seawater after 55 grounding line retreat. This time, the initial amount of 12 C is not assumed to be negligible, due to inheritance of radiocarbon dead organic matter in subglacial sediments from the region (Tulaczyk et al., 1998): In the second phase of the model, when the ice sheet has re-advanced over our sediment sampling sites, we assume that the addition of carbon ceases. Hence, the amount of 12 C is no longer changing but 14 C is experiencing decay, following the standard equation: Where n * is the value of n at the time of ice sheet re-grounding over a field site (when the model switches from phase one to phase two).

Temperature Model Methods
We also develop a two-phase model of ice temperature evolution through ice to compare to observed basal 70 temperature gradients. Phase 1 of the model represents the time when the ice is assumed to have been part of a floating ice shelf, and phase 2 represents the time after the ice shelf has grounded. For both phases we model the temperature profile through the ice using the vertical diffusion-advection equation: Where T is temperature in degrees Celcius, t is time, a * is the diffusion coefficient, w is the vertical velocity, and z is depth from the ice surface. We calculate a * for both phases using the equation: * = 33.8 − 0.3514 (S12) 80 During phase 1, we model the temperature profile of a floating ice shelf, assuming a range of constant ice thicknesses between 500m and 1000m. We initialize the models with a linear temperature distribution through the ice, assuming the temperature at the ice surface to be -25 ºC (Engelhardt, 2004) and calculating the temperature of the ice at the base using the equation (modified from Begeman et al. [2018]): 85 = 0.081 − 0.0568 − 6.858 × 10 /Y (S13) Where S is the salinity, assumed to be 34 PSU from CTD profiles taken at WGZ (Begeman et al., 2018). We use an accumulation rate of 0.15 m/yr (Waddington et al., 2005), and assume a matching basal melt rate (mb) to allow ice shelf 90 thickness to remain constant throughout phase 1. Prior to selecting this basal melt rate we tested multiple values of basal melt rate, which we will discuss further below. Accumulation is equal to the sum of w and change in ice thickness (h). We allow the simulated ice shelf temperature profile to reach steady state, and then begin phase 2 of the model.
Using the resulting profile from phase 1 as the initial condition for phase 2 runs, we model the temperature profile through grounded ice using Eq. (S11) (phase 2). We assume accumulation to be the same (0.15 m/yr [Waddington et al., 95 2005]), but change the basal boundary condition to reflect the freezing temperature of freshwater ice in contact with fresh meltwater: 100 We allow the ice to thicken at the rate of 0.05 m/yr, thus reducing the vertical velocity. In addition, we examined multiple combinations of h and w (Fig. S4), discussed below. We allow the model to run for 10,000 years, keeping track of the temperature gradient of the bottom 100 m for each year. We also measure basal accretion using the equation: (S15) 105 Where HB is the annual change in basal ice thickness, kT is the thermal conductivity of ice (assumed to be 2 W/m K here), TB is the temperature gradient of the bottom 100 m of ice, G is the geothermal flux (assumed to be 0.07 W/m 2 ) here, rI is the density of ice (assumed to be 895 kg/m 3 here), and L is the latent heat of fusion of water (3.34x10 5 J/kg). It should be noted that Eq. (S15) assumes no significant contribution of heat from basal shear heating. 110 We performed sensitivity tests to examine how altering our assumptions for w, h, and mb affected the final basal temperature gradient. For phase 1, we first examined the scenario in which the ice shelf thickness remains unchanged.
Thus, in this scenario accumulation equals both mb and w (Fig. S2). For the four values of w (mb) we examined (0 m/yr, 0.05 m/yr, 0.1 m/yr, and 0.15 m/yr), the model was able to reproduce all of the measured basal temperature gradients only when w (mb) = 0.15 m/yr (Fig. S2). Further analysis revealed that the mb had to be a least 0.13 m/yr to reproduce the steep basal 115 temperature gradients. Additionally, we examined multiple values of mb (0 m/yr, 0.05 m/yr, 0.1 m/yr, and 0.15 m/yr), but kept surface accumulation equal to 0.15 m/yr and allowed the ice shelf to thicken in response to lower mb (Fig. S3). Of the four values of mb that we examined (0 m/yr, 0.05 m/yr, 0.1 m/yr, and 0.15 m/yr), the model only reproduced all measured basal temperature gradients when mb = 0.15 m/yr. Thus, we conclude that we are only able to reproduce the steep basal temperature gradients when basal melt rate of the ice shelf is high, i.e., at least comparable to the surface accumulation rate. 120 We also performed sensitivity tests for w and h for the second phase of the model (Fig. S4). We tested multiple combinations of w and h: w = 0 m/yr, h=0.15 m/yr; w = 0.05 m/yr, h=0.1 m/yr; w = 0.1 m/yr, h=0.05 m/yr; w = 0.15 m/yr, h=0 m/yr. Because accumulation rate (which remained constant) is a sum of w and h, w and h are inversely related to each other. In each scenario we were able to reproduce all measured basal temperature gradients. The steepest basal temperature gradients remained largely unchanged by the changes in w and h. However, for the less steep basal temperature gradients, 125 reducing the amount of ice thickening (h) increases the amount of time the ice has been grounded. Because ice thickening has been observed in this area (Joughin and Tulaczyk, 2002), we discount the possibility of a scenario in which no thickening is happening (h=0 m/yr).

Carbon and Nitrogen Measurements
Total carbon (TC), total organic carbon (TOC), total inorganic carbon (TIC), carbon-to-nitrogen ratios, and d 13 C 130 were measured at the University of California Santa Cruz Stable Isotope Laboratory. The results from those measurements are shown here. In the two-phase model of radiocarbon we use measurements of total organic carbon (TOC) from sediment cores collected at our field sites to constrain our model results (see section 2.3). Carbon-to-nitrogen ratios and measurements of d 13 C were used to glean information about the sources of organic matter found in the sediments and about microbial communities living in the subglacial environments of our field sites. 135

Patterns of Grounding Line Retreat in Ross Sea Embayment
There have been disagreements about post-LGM grounding line retreat in the central Ross Sea (Bart et al., 2018;Conway et al., 1999;Lee et al., 2017;Lowry et al., 2019;McKay et al., 2016;Spector et al., 2017), due to few reliable age constraints from areas covered by the Ross Ice Shelf and the ice sheet itself (Anderson et al., 2014). Previous theories about grounding line retreat in the Ross Embayment have varied, with some following the "swinging gate" model (Conway et al., 140 1999) whereby the grounding line along the Marie Byrd Land side of the Ross Embayment stayed put ("hinged") near the King Edward VII Peninsula while swinging back along the Transantarctic Mountains in the other side of the Ross Embayment. Others have proposed the "saloon door" model, which conjectures that the grounding line began retreating first in the central Ross Embayment and the sides later caught up (Lee et al., 2017;Lowry et al., 2019;McKay et al., 2016).

Exposure age dating along the Scott Coast of the Transantarctic Mountains indicates that the grounding line reached 145
Beardmore and Shakleton glaciers ca. 8000 years ago (Spector et al., 2017). We find that the grounding line retreated over SLW ca. 4000 years ago, which would suggest that the grounding line retreated faster along the Transantarctic Mountains, as described in the "swinging gate" model. However, our results are also consistent with the "saloon door" model. McKay et al. (2016) dated the grounding line to have been located adjacent to Ross Island ca. 8600 years ago. Although our age constraints from the Siple Coast ice streams provide added information about grounding line positions, they do not provide 150 constraints on the geometry of the grounding line during early stages of retreat. Thus, the grounding line could have retreated earliest in the central Ross Sea, while the rapid retreat along the Transantarctic Mountain front between ca. 8600 and ca. 8000 years ago simply allowed the grounding line along margins of the Ross embayment to catch up. Additionally, differences between KIS/BIS and SLW/WIS -namely that the Fm for KIS and BIS are statistically independent from SLW and WIS -suggest that the grounding line behaved differently in those two regions. Additionally, the fact that the Kamb and 155 Bindschadler Ice Stream tills contain similar U-Pb zircon signatures that are absent from Whillans Ice Stream tills (Licht et al., 2014) lends further evidence for differences between KIS/BIS and SLW/WIS.

205
Accumulation is the sum of vertical velocity (w) and ice thickening (h). sediments collected at our field sites. With the exception of the samples collected at UC (italicized), the samples were collected as sediment cores from below either grounded ice or ice shelf. The UC samples were sediments that were melted out of sedimentladen basal ice.