2.1 Ice model

In the coupled ice–ocean model, ice velocities are calculated using the Wavelet-based, Adaptive-grid, Vertically Integrated ice sheet model, WAVI (Arthern et al., 2015; Arthern and Williams, 2017). WAVI is a finite-volume ice sheet model including a treatment of both membrane and simplified vertical shear stresses as described by Goldberg (2011). We use a numerical model with a time step of 20 d and a 2 km horizontal resolution covering the whole Amundsen Sea sector domain shown in Fig. 1a. This horizontal resolution has been found to be appropriate in a recent study, currently under review, testing the impact of resolution on grounding-line retreat and SLR contributions, with a nearly identical configuration of the Amundsen Sea sector (Williams et al., 2024). The ice rheology is described using Glen's flow law, with an exponent of n=3. A Weertman sliding law is used (Weertman, 1964), for which basal sliding drag scales with the cube root of sliding velocity, multiplied by a basal sliding drag coefficient C. As described in Arthern et al. (2015), WAVI uses a data assimilation method to match modelled ice velocities and rates of thickness changes with observations of the MEaSUREs 2014/2015 surface velocities (Mouginot et al., 2017a, b) and rate of change of surface elevation from Smith et al. (2020), resulting in an initial state representing conditions in approximately 2015. Partially grounded cells are utilised using a sub-grid parameterisation to better represent the grounding line (e.g. Arthern and Williams, 2017; Pattyn et al., 2006; Cornford et al., 2012; Seroussi et al., 2014), where the grounding fraction is used to proportionally apply the Weertman sliding drag coefficient. An ice model relaxation is then run for a set period of time (4000 years). During this relaxation the grounding line and the thickness of ice shelves remain fixed, but the grounded ice thickness is allowed to change (see Arthern et al., 2015, for full details). This brings the flux divergence into much better agreement with observations of accumulation and rates of ice thickness change but at the cost of the surface elevation and ice velocities agreeing less well with observations (see Appendix A). Ice thickness prior to relaxation and bathymetry fields are from BedMachine V3 (Morlighem et al., 2020; Morlighem, 2022), where a minimum thickness of 50 m is applied to determine the initial ice extent, which includes the current gaps between the east and west of Thwaites Ice Shelf. The ice front and outer catchment boundaries are kept fixed throughout the simulations. Accumulation and englacial temperatures are kept constant in the forward simulations, using data sets from Arthern et al. (2006) and Pattyn (2010), respectively.

2.2 Ocean model

We use the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997) to simulate the ocean circulation in the Amundsen Sea. The model grid uses polar stereographic Cartesian coordinates with horizontal resolution of 2 km and vertical resolution of 20 m, and we use a time step of 200 s. The ocean model domain extent is shown by the black box in Fig. 1a. Note that this does not include the whole continental shelf, so we apply ocean boundary conditions and velocities to the northern and western open boundaries of the domain to impose prescribed Amundsen Sea conditions in idealised “warm” and “cold” experiments. We also examine a third “no melting” case in which no ocean-driven melting is applied to the ice shelves. This study focuses solely on the ice sheet, ice shelf and ocean interactions driven by wider Amundsen Sea conditions, so sea ice and other freshwater sources/sinks are not included in the simulations, and no atmospheric forcing is applied over the model domain.

The boundary conditions in warm and cold experiments are as follows. Following previous studies, an idealised, piecewise linear, vertical profile is applied to replicate the warm CDW layer below cold Winter Water, with a thermocline 400 m thick (De Rydt et al., 2014). The base of the thermocline is placed at 600 m depth in the warm scenario and lowered to 800 m depth in the cold scenario (Fig. 1c, d). These profiles correspond to the warmest and coldest observed Amundsen Sea conditions, in 2009 and 2012, respectively (Dutrieux et al., 2014). On the northern boundary the warm CDW layer has a temperature of 1.2 °C and a salinity of 34.7 (PSU), whereas the western boundary is forced with a more modified CDW layer, with temperature 0.6 °C and salinity 34.6 (PSU) (Fig. 1c, d). In all simulations, at the boundaries, we apply average ocean velocities (1965–2015) from a larger regional model (Holland et al., 2022; Naughten et al., 2022), and we restore sea level to zero.

Figure 2(a) Initial maximum subsurface potential temperature over the ice–ocean domain for the cold case. (c) Initial melt rate under PIG Ice Shelf in the cold forcing case over the green box in (b). The label shows the total initial meltwater flux. (e) Cross-section through PIG Ice Shelf taken along the west–east green line in (a) showing initial potential temperature in the cold case. (g) Initial melt rate over Thwaites Ice Shelf in the cold forcing case over the red box in (b). (i) Cross-section through Thwaites Ice Shelf taken along the south–north red line in (a), showing initial potential temperature in the cold case. Panels (b), (d), (f), (h) and (j) as in (a), (c), (e), (g) and (i), respectively, but for the warm forcing case.

Ocean-only simulations, with fixed initial ice geometry, are simulated for 2 years as a spin-up before coupling. The resultant conditions in the ocean model are designed to match the spatial distribution of present-day observed maximum subsurface temperatures over the region (Dutrieux et al., 2014), as shown in Fig. 2a and b. The warm and cold oceanic forcing cases enable different amounts of CDW to reach the base of the present-day Thwaites and Pine Island Glacier (PIG) ice shelves (Fig. 2e, f, i, j). The initial melt rates, after the 2-year ocean-only spin-up, of the Thwaites and PIG ice shelves for the warm and cold simulations are shown in Fig. 2c, d, g and h.

Ice shelf melting is represented by a standard three-equation formulation (Holland and Jenkins, 1999). We tune the dimensionless ice shelf melting drag coefficient in this parameterisation to 0.008. This drag coefficient parameterises the ocean stress on the ice base as a function of the ocean model's mixed layer velocities in order to calculate turbulent ocean heat and salt fluxes for use in the melting calculation (Jenkins et al., 2010). For the PIG and Thwaites ice shelves combined, this value produces the closest average match between the initial MITgcm ice shelf melt rates and those that are implicit in the WAVI initialisation (see Appendix B). This minimises the “coupling shock” – the response of the ice model to any mismatch between these two fields – which occurs when the coupled model simulation commences. Without this calibration, the ice sheet trajectory could be impacted, potentially for many decades (Goldberg and Holland, 2022), by the adjustment of the ice due to the transition from implicit initialised melt rates to arbitrarily different ocean model melt rates at the start of the simulation. Note also that this tuned value is close to a value of 0.01 derived from observations (Jenkins et al., 2010).

BedMachine V3 is used as the basis for the bed and seabed geometry throughout the model. However, we found that without modification this dataset led to rapid grounding and advance of PIG. Closer inspection revealed that the seabed bathymetry may have been underestimated near the grounding line of this glacier in BedMachine V3 (Appendix C). Therefore, throughout the domain we deepen the seabed wherever needed to achieve a minimum water column thickness of 280 m on the staggered Arakawa C-grid's velocity grid points, located on the grid faces, tapering towards grounded ice over 6 km down to a minimum water column thickness of 140 m. This procedure is only applied to cells with no ice basal sliding drag in the initial state and is only done once at the start of the simulation before the WAVI ice sheet model is initialised and relaxed rather than being an ongoing process. This edited bathymetry is used in both MITgcm and WAVI.

2.3 Ice–ocean coupling

This study uses a new, synchronously coupled ice–ocean model, which combines the WAVI ice sheet model with MITgcm. Crucially, this coupling occurs through the MITgcm ice sheet package STREAMICE (Goldberg and Heimbach, 2013), making use of previous synchronous coupling developments (Jordan et al., 2018; Goldberg et al., 2018), which enables the STREAMICE ice thickness to evolve continuously in the ocean model. The basic concept of the new model is that WAVI solves for the ice velocities; MITgcm STREAMICE solves for the ice thickness; and the MITgcm ocean model solves for all ocean properties, including melting and also where the ice thickness evolves. Dividing up computations in this way allows for a fully synchronous coupling because both the ice thickness and ocean properties are solved within MITgcm on the ocean time step. This has two advantages: firstly, the ice thickness and ocean free-surface equations are solved simultaneously, allowing full conservation of mass in the ice and ocean coupled system. Secondly, the melt rate responds instantly to changes in ice thickness.

In principle, this approach is no more expensive than other coupling approaches because the two-dimensional ice thickness equations can be solved on the ocean time step with negligible computational expense compared to that incurred in solving three-dimensional ocean equations. However, the approach does require the MITgcm grid to exist wherever ice may go afloat during the simulation, since grounding-line retreat in MITgcm is accomplished naturally. MITgcm solves the ocean free-surface equation every ocean time step in order to conserve mass, which naturally inflates the water column in the MITgcm ocean model wherever the pressure loading of the ice decreases below floatation (Goldberg et al., 2018). In regions of ice which are not floating, a thin subglacial layer is specified in order to enable the expansion of the ocean column during grounding-line retreat (Goldberg et al., 2018; Jordan et al., 2018). We set it to be 4 m thick, but it could have been set to any relatively small thickness compared to the ocean model vertical resolution. This small value has been previously demonstrated to have no impact on the evolution of the coupled system (Goldberg et al., 2018). This layer, which has no effect on basal sliding drag computed by the ice model, is treated as a porous medium, with Darcy flow used to connect the subglacial cells to each other (Goldberg et al., 2018).

Using an MITgcm domain covering the entire catchment of the Amundsen Sea sector would be very inefficient, since it would include large areas of grounded ice that would never go afloat in the ocean/STREAMICE grid. To avoid this, we introduce the concept of “ice-only” and “coupled” domains (Fig. 1a): in the ice-only domain, WAVI is used to solve for both ice sheet velocities and thickness, while in the coupled domain, WAVI solves for ice velocities and MITgcm solves for ice thickness and ocean properties, as described above. The coupled domain only needs to extend far enough inland to accommodate the grounding-line retreat occurring during a projection, which for this study was determined using test simulations.

The coupled model solution procedure is split into coupling periods, chosen to equal the WAVI time step of 20 d, meaning that the WAVI model state is fixed during each coupling period. Over each coupling period (Fig. 1b) the following takes place: (1) WAVI calculates ice velocities over the whole ice domain for one 20 d time step; (2) WAVI steps forward the ice thickness in the ice-only domain for one 20 d time step; (3) MITgcm receives the ice velocity for the coupled domain and the updated ice thickness on its boundaries and sub-cycles the ice thickness, grounding line and ocean properties within the coupled domain using the 200 s ocean time step for the full 20 d coupling period; and (4) MITgcm passes the new ice thickness and therefore the new ice grounding line for the coupled domain back to WAVI, and the next coupling period commences. Therefore, the choice of the length of the coupling period determines the fastest response time for which the ice velocities can respond to changes in buttressing in the coupled domain, where the ice thickness changes with the ocean time step. The WAVI boundary cells outside of the coupled domain are passed to the MITgcm domain and held fixed while MITgcm runs, and then the thickness throughout the coupled domain is passed from MITgcm to WAVI after each coupling step. This procedure, in addition to WAVI being fixed during each coupling period, keeps the two models'/domains' ice thicknesses from diverging and ensures a smooth transition between the two domains.

For consistency, the same advection scheme that is used in WAVI (Arthern et al., 2015) has been coded into STREAMICE. Ice divergences are updated with the changing ice thickness in STREAMICE every ocean time step, while the ice velocity remains fixed over each coupling period. We adopt the principle that the ocean may only melt ice that has no basal sliding drag applied (Arthern and Williams, 2017), so the WAVI ice basal sliding drag field is passed to MITgcm to decide where melting can occur during each coupling period. This means that no melting can occur on partially grounded cells, as recommend by Seroussi and Morlighem (2018). In addition, this means that if a cell becomes fully ungrounded in MITgcm, ice shelf melting only occurs once the coupling period finishes and a new one starts, which updates the ice thickness in WAVI and subsequently passes back a new basal sliding drag field. However, if, during a coupling period, a grid cell becomes grounded in MITgcm, melting on this cell is immediately switched off. In both the WAVI and STREAMICE models a minimum ice thickness of 50 m is applied.

Figure 3Warm forcing case evolution over 180 years. Panels show snapshots of (a) ice thickness, (b) ice speed and (c) melt rates plotted on a logarithmic scale, taken every 36 years throughout the simulation. Initial (black) and final (red) grounding lines are included.

3.1 Simulated evolution of the Amundsen Sea sector

Following initialisation, the coupled ice–ocean model is run forward for 180 years under both forcing cases, as well as a zero-melting case. The evolution of the warm simulation is shown in Fig. 3. The most prominent feature of the simulation is the acceleration (Fig. 3b) and grounding-line retreat of Thwaites Glacier (Fig. 3a). At the start of the simulation, Thwaites Ice Shelf hosts the observed fast-flowing western shelf and slower eastern ice shelf (Fig. 3b), where a pinning point restrains the flow of ice. However, as the simulation progresses, both sides of the ice shelf accelerate, reaching speeds of up to ∼10 km yr⁻¹, but then the whole ice shelf shows signs of deceleration between 144 and 180 years. In addition, Thwaites Glacier experiences substantial grounding-line retreat during this simulation, approximately 130 km of north–south retreat, leading to a reduction in the total mass of grounded ice and the formation of a new, larger ice shelf cavity, which features high melt rates near the grounding line (Fig. 3c). The fixed ice front leads to a large and thin future Thwaites Ice Shelf, which may lead to artificially elevated melting overall, though only low ice shelf melt rates occur on the thinner ice (Fig. 3a, c). However, such trends in ice speed and grounding-line retreat are restricted to Thwaites: neither PIG nor Smith Glacier experiences the same level of acceleration or retreat, although we do observe an initial acceleration of the PIG (which is later reversed) and some grounding-line retreat in both PIG and Smith.

In order to analyse SLR contributions, melting and grounded ice area loss from different regions and to compare the three forcing cases, we divide the domain into three areas: Thwaites, PIG and Smith (Fig. 1a). These areas are edited from the Zwally et al. (2012) basins, where the Smith area has been separated out. In addition, the boundary between the PIG and Thwaites areas is edited by hand to address the otherwise mis-assignment of SLR contributions, melting and grounding ice loss to PIG as Thwaites retreats (Fig. 3). The model projects notable SLR contributions from all three regions, even in the zero-melting case for the Thwaites and Smith areas, which has important implications. This result demonstrates that the initialised ice sheet state in WAVI is intrinsically unsteady and will continue the present-day ice loss for some period of time irrespective of future climatic forcing. This implies that a tipping point may have been passed at some point for this particular initialised ice sheet, likely during the 20th century (Mouginot et al., 2014), and some amount of SLR (at least ∼40 mm in our simulations) is now committed. Zero melting can never occur in the real world because even seawater at the surface freezing point will drive melting at depth due to the pressure decrease in the freezing temperature, so this is a very conservative test of the presence of committed ice change. In the PIG area there is a negative SLR contribution for the zero-melting case (mass gain), though the contribution from PIG is dependent on the particular choice of bathymetry deepening that is implemented (Appendix C). Notably, we observe the largest zero-melting SLR contribution from the Thwaites area, leading to ∼45 mm in 180 years, with a near-constant SLR rate. For the Thwaites area we also obtain some grounding-line retreat, though the rate of loss reduces during the simulation, suggesting that the SLR contribution would stabilise after a further period of simulation.

Figure 4Time series shown for the three cases: no melting (black line), cold (blue line) and warm (red line). (a) Cumulative global sea level contribution from the Thwaites area. (b) Rate of change of global sea level contribution for the Thwaites area. (c) Total ice shelf melt rate from the Thwaites area with a 2-year running average applied. (d) Cumulative grounded area loss from the Thwaites area. (e–h) Same as (a–d), but for the PIG area. (i–l) Same as (a–d), but for the Smith area.

Download

The cold and warm cases are far more realistic since they are based on observed Amundsen Sea conditions. Both scenarios contribute much more SLR than the zero-melting case (Fig. 4a, e, i). In both cases, the Thwaites area dominates SLR contributions from the sector, providing approximately two-thirds of the total, and responds differently to melting than the PIG and Smith areas. In the PIG and Smith areas, the rates of SLR with melting are approximately constant but are higher than the rates of SLR in the zero-melting case, increasing to ∼0.3 and ∼0.1 mm yr⁻¹ in the warm case, respectively (Fig. 4f, j). However, in the Thwaites area with melting a trend emerges in the SLR rate (Fig. 4b), increasing from ∼0.2 mm yr⁻¹ at the start of the simulation to ∼1.1 mm yr⁻¹ at the end in the warm case. In addition, we observe rapid jumps in the SLR rate for Thwaites that are not present in the other areas. The increase in the SLR rate in the Thwaites results in an approximately quadratic SLR contribution (Fig. 4a) rather than the more linear SLR contributions obtained for PIG and Smith. In all three areas “noise” is present in the SLR rates due to the instant effects of the changing basal sliding drag field, though this is harder to see for the Thwaites area due to the larger y-axis scale. In addition, we observe an increase in the total melting that occurs from the Thwaites Ice Shelf as the simulation progresses: over the 180 years of the simulation, the total melt flux from Thwaites increases by an order of magnitude, from ∼45 (∼25, respectively) Gt yr⁻¹ at the start of the simulation to over ∼440 (∼210) Gt yr⁻¹ at the end in the warm (cold) case (Fig. 4c), which occurs despite constant oceanic boundary forcings during the simulation. This is explored further in Sect. 3.3. In the PIG and Smith areas, the total melt flux is approximately constant (Fig. 4g, k). However, we do observe a strong correlation between the cumulative integral of melting and the SLR contribution individually for each region (Appendix D). In addition, while the effect of the inclusion is small, we do include ice shelf melting that occurs on ice that has reached the minimum thickness, as while this melting is not allowed to thin the ice further, it still has a partial glaciological effect by stopping the ice from thickening. Furthermore, this melting still has an oceanographic effect in the simulation by cooling and freshening the ocean.

For all regions, we obtain consistently higher SLR rates and grounded ice area loss in the warm forcing case compared to the cold forcing case (Fig. 4b, f, j). However, as a fraction of the total SLR in these models, the difference between warm and cold scenarios is remarkably small, only ∼27 % of total SLR from all regions by the end of the simulations. Given that these scenarios bracket the coldest and warmest ocean conditions on record in the Amundsen Sea, this suggests that the future SLR from this region is only weakly influenced by variations within the observed range of present-day ocean conditions: while melting is important to the SLR, its typical climatic variations are less so. Despite this, in the Thwaites area we observe different timings and strengths in the jumps in the SLR rate between the cold and warm cases (Fig. 4b). We also observe an increasing difference between the cold and warm total melt rates for Thwaites Ice Shelf (Fig. 4c), while for the PIG and Smith areas this difference remains approximately constant (Fig. 4g, k). Additionally, there is a larger relative difference between the melting and no melting cases for the PIG area compared to Thwaites. In the PIG area, ice shelf melting prevents the ice shelf from thickening and re-grounding on the bathymetric ridge below it. In the no melting case, the ice shelf re-grounds on this ridge; the buttressing provided by this re-grounding leads to the large differences between the no melting and melting cases for the PIG area (though this is dependent on bathymetry deepening). While there is a smaller relative difference for the Thwaites area, the presence of melting still has a significant impact, changing the SLR contribution from linearly increasing to quadratic.

Overall, we conclude that melting and SLR have a fundamentally different response in the Thwaites area than the other areas, leading to an increasing SLR rate and total melting rate. In addition, the Thwaites area dominates the SLR contribution for the Amundsen Sea sector over the 180 years simulated. Therefore, the remainder of this study focusses upon the processes underlying this behaviour in the Thwaites area.

Figure 5(a, b) Grounding-line retreat rates over the 180-year simulations in the cold (a) and warm (b) cases indicated by colours, with colour bar shown in (b). The red contours in (a), (b), (d) and (e) and the purple contours in (c) and (f) show the presence of isolated pinning points during the simulation. (c) Percentage change in grounding-line retreat rates between the cold and warm cases (i.e. the percentage difference between those retreat rates shown in (a) and (b), respectively). Bathymetry (d) and Weertman C coefficient (e) shown under initially grounded ice with final extents of the grounding-line retreat in the warm (magenta) and cold (cyan) cases. (f) The difference in grounding-line retreat rates between the cold and warm cases (colours), as well as final grounding-line area (grey area). The area shown in each panel corresponds to that displayed as a black box in Fig. 3.

3.2 Thwaites Glacier retreat and ice shelf pinning points

Figure 5 shows quantities related to Thwaites grounding-line retreat, which is one of the key features that differentiates this area from the PIG and Smith areas. Over the simulation, the warm case has a larger area of grounding-line retreat than the cold case, with a faster rate on average (Fig. 5a, b, c, f). Grounding-line retreat rates are calculated from the discrete migration of the grounding line across grid cells. Specifically, the retreat rate in each grid cell is calculated as ΔxGL/ΔtGL, where Δt_GL is the time between when the subject cell first becomes a grounding-line cell (any of the eight adjacent cells are floating) to when that cell ungrounds, and Δx_GL is the grid cell width or diagonal extent depending on which adjacent cell was first floating. Figure 5d and e show the bed depth and “Weertman C” drag coefficient over the ice area that is initially grounded. A lower Weertman coefficient corresponds to a more slippery bed. In addition, the bathymetry in this area generally deepens inland, which could promote grounding-line retreat (Weertman, 1974; Schoof, 2007). In both warm and cold cases, highly heterogeneous retreat rates are observed (Fig. 5a, b), with areas of fast retreat as high as ∼10 km yr⁻¹ in proximity to areas of much slower retreat. Towards the end of the warm case simulation, the grounding line experiences rapid retreat (Fig. 5b) across a deep and slippery bed section (Fig. 5d–e) before slowing down as it encounters a ridge of shallower and less slippery bed (Fig. 5d–e), where it remains until the end of the simulation. These features explain the large variations in the Thwaites area SLR rate in the last ∼50 years of the simulation (Fig. 4b). In particular, the retreat onto the shallower and less slippery ridge at the end of the warm simulation decreases the ice flux across the grounding line and corresponds with the Thwaites region's decreasing SLR rate during the last 25 years of the simulation.

Comparing the cold and warm scenarios we see the effect of increased melting on grounding-line retreat rates. More specifically, we observe areas of elevated grounding-line retreat rate in the northern part of the ungrounding area, below which is an area with minimal sensitivity to melting scenario, and again increased rates in the south (Fig. 5f). There are very large percentage increases in grounding-line retreat rates in two clear bands on the retreated area (Fig. 5c), where two bands of slow retreat rates in the cold case are not present in the warm case.

In addition, we observe the formation of many pinning points as Thwaites Glacier retreats, shown by red and purple contours in Fig. 5 for the different cases. An ice grid cell is flagged as a “pinning point” if it is grounded but separated from the main grounded ice sheet by floating ice. The red and purple contour is then drawn around all cells that are flagged as pinning points in any one of the outputs. The pinning points are generally located on areas of shallower bathymetry that are downstream of areas of deeper bathymetry (Fig. 5d). The pinning points typically feature slower retreat rates, and inshore of these pinning points we generally observe areas of faster retreat (Fig. 5a, b). These pinning points are crucial for the future evolution of Thwaites Glacier, as they provide drag and buttressing as the ice retreats. The importance of these pinning points is determined by their size, position and duration, so we examine these features and their relationship to SLR rates from the Thwaites area.

Figure 6(a) SLR contribution rates for the Thwaites area for warm (red), cold (blue) and no melting (black) cases. The background colour indicates the time in the simulation according to bands of 25 years with a longer starting (ending) band of 45 (35) years. (b) Colours indicate the band of the final year that pinning point cells exist for the cold melt case, according to the colours shown in (a). Note that only pinning points that exist for over 1 year are shown. Red labelled boxes refer to groups of pinning points discussed in the text. (c) Same as (b), but for the warm case. (d, e) Pinning point duration, taken as the time from isolation to ungrounding, for cold (d) and warm (e) cases. (f) Difference in pinning point duration between cold and warm cases (i.e. the difference of the data shown in d and e) where pinning point locations match. In (b)–(e), grey regions indicate areas which do not become ungrounded during the simulation. The area shown in panels (b)–(f) is shown with a black box in Fig. 3.

Figure 6b–c show the ungrounding time of pinning points in the cold and warm cases, respectively. By examining these times, we can compare how the ungrounding of pinning points relates to jumps in SLR rate that we observe from the Thwaites area. The labelled red boxes in Fig. 6b show the key groups of pinning points in the simulations.

Within the first 60 years, the SLR rates of the warm and cold forcing cases diverge (Fig. 6a). In the warm case, pinning point group “a” completely ungrounds between years 25 and 60 (Fig. 6c), with the final ungrounding of the last pinning point in group “a” coinciding with the ungrounding of group “b”. The removal of the combined associate buttressing leads to a large jump in SLR rate in the warm simulation at ∼60 years. However, in the cold simulation, some of these pinning points in group “a” remain grounded until beyond 95 years (Fig. 6b), leading to a period of relatively steady SLR rate in this simulation (Fig. 6a). In the zero-melting case, part of pinning point group “a” and all of “b” remain grounded throughout the simulation (not shown), and SLR rates remain approximately constant as a result (Fig. 6a). In the warm simulation, group “c” becomes ungrounded by ∼130 years, causing a rapid increase in SLR rate, and the loss of group “d” at ∼145 years leads to another large jump. The cold simulation loses group “c” at ∼170 years, leading to a small jump in SLR rate, while group “d” remains grounded at the end of the simulation.

Figure 6d–e show the duration of the pinning points – the time between separation from the main body of grounded ice to ungrounding. There is a large variation in pinning point durations, with some lasting less than a year, while others persist for decades. In general, the duration of pinning points is lower in the warm case (Fig. 6d–f), with an average pinning point duration of ∼6 years over matching pinning locations, compared to ∼10 years in the cold case, leading to the different ungrounding timings and SLR rates as described above. Reducing the duration of pinning points increases the intensity of periods of rapid ice acceleration and grounding-line retreat, leading to the differences in the spatial map of grounding-line retreat rates between the cold and warm cases shown in Fig. 5c and f.

Figure 7Evolution of conditions as grounding line retreats over pinning points. Snapshots every 9 years starting from the date of group “c” pinning point formation in the warm and cold cases (118 and 144 years, respectively). (a–c) Cross-section through the Thwaites Ice Shelf, taken along the red line in (g) extended to ice front for the warm case. (d–f) Same as (a–c), but for the cold case. (g–i) Melt rates for the warm case, with arrows showing ice shelf boundary layer ocean velocities. (j–l) Same as (g–i), but for the cold case. (m) Ice speed snapshots at start of pinning point formation, and (n) and (o) show differences for the warm case. (p–r) Same as (m–o), but for the cold case. Area shown in (i)–(x) is shown in Fig. 3 as a small green box.

To illustrate the key role of pinning points, we now focus on the ungrounding of group “c”, which causes the rapid jump in SLR rate at year 130 in the warm simulation and the smaller jump at year 170 in the cold simulation. Figure 7 shows the evolution of the two simulations, starting from the individual dates of formation of the group “c” pinning points. Figure 7a–l show the ice geometry and ocean conditions throughout the subsequent evolution. To begin with, the grounding line is located in shallower bathymetry at the top of a retrograde slope (Fig. 7a, d). Ungrounding then occurs laterally around this point, including upstream, encircling the pinning point and leaving it isolated from the rest of the grounded ice (Fig. 7h, k). Note that the pinning point is grounded on the side of the shallower bathymetry rather than on top of the bathymetric feature. With the grounding line now in a deeper bed an acceleration is expected, although the pinning point continues to provide basal sliding drag that resists the flow. The advection of thicker ice from the deeper bed upstream enables the ice to remain grounded on the shallower bed beneath the pinning point. Therefore, only a small acceleration is observed at 9 years, while the ice remains grounded on the pinning point (Fig. 7m–r).

However, melting now occurs in the newly opened cavity upstream of the pinning point. This thins the ice, enlarging this cavity, enabling greater oceanic connection, and leading to higher ocean velocities and melt rates (Fig. 7h, k). This melting and thinning feedback eventually leads to sufficient thinning upstream of the pinning point for it to unground completely (Fig. 7i). The resulting loss in buttressing leads to a large increase in ice speed (Fig. 7o), further rapid grounding-line retreat and a jump in SLR rate. Therefore, the distribution and strength of localised melt rate patterns strongly determine the duration of these pinning points and thus the overall ice retreat. This is shown clearly in Fig. 7g–l, with the warm case leading to higher melt rates and a much faster ungrounding of the pinning point. The melt rates around these pinning points are highly heterogeneous (Fig. 7g–l), with elevated melt rates typically occurring on their eastern side, where rapid ocean currents drive high melt rates. These rapid currents occur where buoyancy-driven meltwater flow is trapped against the grounding line by Coriolis force (Holland and Feltham, 2006).

Figure 8Cross-sections of potential temperature beneath Thwaites Ice Shelf for cold (a) and warm (b) cases, taken every 36 years, along the same section as shown in Figs. 7 and 2i and j, which is represented in Fig. 2a as the red line.

Download

3.3 Thwaites Ice Shelf geometric melting trends

An increasing trend in the total ice shelf melting occurs over the evolving Thwaites Ice Shelf during both the warm and cold forcing simulations (Fig. 4c). The wider oceanic forcings are fixed, so these trends must be driven by geometric changes in the ice shelf cavity. In Fig. 8 we show how the Thwaites cavity geometry and ocean conditions evolve in the two cases. The warm case has a higher thermocline, so the warm CDW is more easily able to flow over ridges and flood the new cavity areas. In the cold case, the access of the warmest CDW is blocked by seabed highs, with only more heavily modified CDW reaching the ice base. However, sufficiently warm water is still able to drive melting close to the grounding line. In both cases the ice base near the grounding line remains steeply sloped throughout the retreat.

Figure 9(a) Total melt flux from the full Thwaites Ice Shelf area (black) and from grid cells located below 600 m depth only (blue). (b) Total melt flux from grid cells below 600 m and total ice area below 600 m. (c) Average melt rate on grid cells located below 600 m depth and spatial standard deviation of these melt rate values.

Download

Figure 9a shows the evolution of total melt flux from the Thwaites Ice Shelf in the warm case for both the entire ice shelf and for ice below 600 m depth only, which is the thermocline depth in this case. For the majority of the simulation, the total melt flux from the entire ice shelf increases, but it decreases during the last 25 years. Most of the trend in ice shelf melting occurs in the deeper ice, with melting below 600 m peaking at an increase of ∼30 times its initial value, which suggests that as the grounding line retreats, an increase in ice shelf base area below the thermocline controls total melting (Fig. 9b). This confirms that the trend in melting does not result from the increasing ice shelf area associated with an artificially fixed ice front. The increase in deep ice area occurs because the groundling line retreats into deeper bathymetry but also because the slope of the ice shelf base, below 600 m, gets shallower during the simulation, as shown for example in Fig. 8. Without an increasing trend in melt flux, the thicker ice advected across the grounding line as the ice sheet retreats into deeper bathymetry would result in re-grounding on pinning points further downstream. However, towards the end of the simulation the groundling line retreats slowly onto a shallower ridge (Fig. 5d), where the shallowing grounding-line depth decreases the ice shelf area below 600 m and subsequently decreases the total amount of melting below 600 m. The close correspondence in the increase between total melting (10 to 197 Gt yr⁻¹) and ice base area (180 to 3560 km²) implies that the average melt rate (m yr⁻¹) beneath deep ice below 600 m is approximately constant at ∼45 m yr⁻¹ (Fig. 9c). However, we do find some temporal variability in this average melt rate, with it varying up to 50 %, and importantly the local melt rates are highly spatially variable, with a spatial standard deviation close to the mean value. This spatial variability is important for the effect of melting on small-scale pinning points, as described above.