| 02 Mar 2026
Modelling the evolution of Thwaites Glacier over the 20th century
Abstract. Thwaites Glacier is rapidly evolving and could make large sea-level contributions in the coming centuries, making it essential to understand the drivers of the ongoing ice loss. Sediment-core analysis suggests that Thwaites Glacier was in a relatively steady state for millennia before its western pinning point ungrounded in the 1940s. Here, we include a first analysis of 1947 aerial imagery of Thwaites Ice Shelf, which shows that it was relatively undamaged, contrasting with the highly-damaged present-day. Additionally, the main outflow and shear margin were displaced ~15 km westwards compared to the present day. We use the MITgcm-WAVI coupled ocean-ice sheet model to create example quasi-steady pre-1940s configurations for Thwaites Glacier, including a most plausible pre-1940s state, finding that healing the damaged ice shelf is necessary to achieve this. Next, we trigger ice retreat and highlight key processes as the model evolves into the present-day configuration, including ice damage, pinning-point ungrounding driven by ocean melting, and ice piracy between eastern and western parts of Thwaites Glacier. By conducting reversibility experiments during the retreat, we find that multiple quasi-steady ice-sheet states are possible under the same ocean forcing, demonstrating the potential for tipping points in the Thwaites system. Either ice damage or increased ocean forcing can eliminate these quasi-steady states, prompting retreat resembling that observed today. Taken together, these results demonstrate that the sea-level contribution from Thwaites Glacier is not simply controlled by ocean warming in the Amundsen Sea, and is highly sensitive to ice-damage feedbacks, which must be incorporated into sea-level projections.
General comments:
Bett et al. successfully use coupled ice sheet-ocean modelling and historical observations to investigate the response of Thwaites Glacier to idealised perturbations in ocean forcing and ice-shelf damage. The authors find that an advanced glacier is only achievable with an extended ice front and ‘healed’ ice-shelf damage, allowing the western ice shelf to ground on a shallow bathymetric high. When ocean conditions are altered, the glacier retreats towards its present-day state and beyond, at a rate that increases with added damage. The authors also demonstrate ‘ice piracy’, with eastern Thwaites starved of ice during retreat, which could explain observed thinning despite low melt rates. Reversibility experiments indicate several different quasi-steady states depending on the level of ice-shelf damage and ocean forcing, highlighting the complex feedbacks that may have governed glacier evolution over the 20th century. These results have direct relevance for other modelling studies focussed on the evolution of Thwaites, whether initialised from its present-day damaged state or a historic undamaged configuration.
Specific comments:
Overall, the manuscript is generally well-written, but some sentences need to be reworded, and descriptions of model simulations could be condensed to help readability. I recommend using experiment names for the different simulations to help the reader keep track of what is being discussed and reduce the amount of words being used to describe damage and ocean forcing. The repeated use of 'ocean-only forced case' is also potentially confusing given the concurrent use of 'ice-only simulation'.
Whilst the authors state this is an example 1940s configuration, the title frames the work as modelling Thwaites' evolution over the 20th century. However, given the idealised nature of both the ocean and damage forcing, this characterisation seems difficult to justify. The authors show thathealing the ice-shelf damage is required to advance the grounding line, but the impact of the initialisation procedure, whereby the model is matched to present-day velocities and surface elevation changes, remains unclear. The model setup is also insufficiently described for both the full and reduced physics configurations; this detail could be included in an appendix. Similarly, the damage calculation methodology lacks sufficient detail, despite being central to the study.
The reversibility experiments section requires further justification. It is unclear whether 250-year simulations are sufficient to characterise tipping point behaviour or map the stability landscape meaningfully. Furthermore, buttressing from other sources, notably Crosson Ice Shelf and Bear Island, receives only brief mention in the discussion, despite their potentially significant influence on Thwaites dynamics.
The figures could be improved, particularly regarding distinguishing east versus west Thwaites and the pinning points under discussion.
Technical comments:
L1: The title is not very specific or descriptive for this study, there is no mention of ice-shelf damage, which you show is one of the key drivers of historical retreat
L20: “demonstrating the potential for tipping points in the Thwaites system” - do you mean tipping points in the historical position or in a future state?
L40: “modelled to occur”?
L66-67: Reword this sentence because it sounds like you’re referring to the “sketched outlines” in the “image” in figure 1.
L68: “relatively undamaged for centuries” - are you assuming this from flow speeds of km’s per year and the length of the iceberg in the 1973 image?
L71-74: This sentence seems out of place, I would expect this at the end of your introduction.
L76-77: “Changes in ... melting” compared to when?
L78: modified Circumpolar Deep Water (mCDW)?
L90: Remove ‘the’ before ‘increasing’: “the effect of the increasing ice damage on ice flow speed”
L98: Why have you switched from using “ocean-ice sheet model” in line 15 to “ice/ocean” here and throughout the paper? (e.g., 105, 147, 189, ...). Then in line 101 you use “coupled ice-ocean model”.
L98: Remove the extra “coupled” before “model”
L98: As you only use idealised ocean boundary conditions (and forcing?), is it correct to say the simulations are over the “20th century”?
L102: Does “These simulations” refer to the simulations mentioned in line 98? (i.e., not in Bett et al. 2024?)
L112-113: Put this sentence after the next one. Introduce figure 2 first, then say that the lines are approximate and you can identify key features.
L119: Change “is much more advanced” to “has a longer extent” or “with an ice front that is much more advanced”.
L120-121: Difficult to see where the “main western ice shelf” is in Fig 2 – Can you add a 1973 coastline to aid visualisation?
L125-126: If the main 1940s grounding line was more advanced (as you suggest later in the paper), could the damage in (g) and (f) be due to fracturing at the grounding lines and pinning point, respectively, as we see in present day (e.g., Wang et al. 2025, Wild et al. 2022)? Or is this more indicative of shear margin damage?
L129: What is the uncertainty in 1940s image location? How does it compare with the 15 km between possible shear margin locations?
L130-131: Refer to figure 1 here
L157-160: What is the duration of the coupling period for velocities, if it is the same as the timestep does it vary throughout the simulation or remain constant? How often are the melt rates given to the ice model, is that on the ocean model timestep?
L162: and Kohler, Pope? Do you keep present-day ‘damage’ on Crosson and Dotson ice shelves?
L162: You use “PIG” here but haven’t said what that is beforehand.
L164: Put “e.g.,” inside the citation brackets
L168-169: Does the ice-sheet model include the surface/thickness elevation change (as well as velocities) in the inversion process?
L169: Fix citation formatting
L170: What is this C value representative of, how did you select it? Unless I missed it, I couldn’t find it in Bett et al. (2024).
L172: What is the match with observations after the relaxation step?
L174: Is this temperature field what you later refer to in line 214 in the calculation of the viscosity enhancement factor? Could you show a plot of the field in an appendix/supplement?
L175: Why have you changed the bathymetry for Pine Island Glacier cavity?
L188: You mention different types of “simulations” but then call them different “models” in the next sentence which is a bit confusing.
L188: What is the reduced-physics coupled model, is it still WAVI-MITgcm? You’ve written “see below” but this doesn’t seem to be described anywhere, unless the only difference is the eddy viscosity. Also, is the full-physics coupled model also WAVI-MITgcm?
L190: Remove the extra “initial”
L198: It might be helpful to mention this coupling period in line 157-158 (see earlier comment)
L207: Why a thickness of 50m?
L208: What does this mean: “removing restrictions that would have been placed by attempting to explicitly recreate the convoluted pre-1940s ice front”
L213: How is this “viscosity enhancement factor (E)” calculated?
L215: So, assuming “damage is the primary source of viscosity enhancement”, using E=1-D, more damage yields lower viscosity (and vice versa)?
L216: If you are using Bedmachine (2015) geometry in the initialisation, would it better to say the ‘2015 state’ rather than ‘present-day’, particularly if the 2015 and 2024 damage fields are different. In Fig 1 you show a highly damaged (non-existent) western tongue but presumably that isn’t what your damage field in Fig 3 represents? Can you include a picture of the approximate 2015 state in a supplement? This should then show a damaged shear margin that coincides with your damage field in Fig 3d.
L217: Are the “high damage values” close to 1?
L218: “as might be expected (Figure 1b)” – but the image you show in Fig1b shows a non-existent western ice shelf with virtually no shear margin
L220: Are there tunable parameters in the model initialisation that would impact the resulting damage field?
L225: Change ‘our’ to ‘the’: “During our evolving…”
L225-227: Needs rewording, it is difficult to understand what you are saying here. I presume the “plus” doesn’t mean you are adding something to the present-day damage field, but refers to the other simulation.
L227-228: I’m not sure what you mean here, “1 decimal point rounded average damage…” – can you make it clearer
L228: Are you assuming that the average grounded present-day damage is a maximum damage value for the 1940s?
L236-237: Does the damage change over time (“each coupling period”)?
L238-242: This is a little confusing and doesn’t explain what you’re doing in your later experiments.
L245: Where does “see below” refer to?
L262: Is the damage intensity factor ever set to values other than 0 or 1?
L264: Is the healing of the ice shelves instantaneous?
L269: The present-day configuration is the state after initialisation (inversion and relaxation)?
L271: Why is there a reference to Holland et al. 2023 here, when you are showing your model result?
L272: “Thwaites retreats” - Do you mean different examples of retreat?
L274: As you refer to the “reference pre-1940s" case throughout the paper, can you use experiment names to help the reader keep track of the different simulations? Try to include the thermocline depth and damage value in the name. A summary table would help too.
L296: “unedited ice geometry” - Can you just say 2015 geometry (or ice front)? Unless you are changing other geometric variables?
L298: Same as the above comment, does “edited setup” just mean the extended ice front and healing? Can you make this more explicit/obvious?
L300-302: Although you do not show it, have you tested whether only healing or extending the ice front (both with no melting) causes advance?
L309: Remove “as labelled in Figure 4c”
L335: “levels of damage” - Do you consider damage fields between intensity factor of 0 and 1?
L339: The ramping up of the damage over 5 years is just for the second case? Does it ramp up from the ‘healed’ case to present-day damage?
L353: Remove extra bracket at the end of citation.
L353: Add ‘a’: “We also observe a separate slow ...”
L362: Can you rename this subtitle to be like 4.2.1, such as: “Retreating Thwaites Glacier with ocean melting and damage”
L365-366: Clear up the wording slightly, maybe: “During the damage plus ocean forced simulation, there is faster grounding-line retreat back to the observed 1996 position...”
L366: Why is there a smaller modelled grounding line on the east ridge in the 1940s state compared to observed 1996 and present day?
L384: “initial pinning point” - which one?
L385: Refer to fig 6a here
L388: Use “observed” instead of “real”
L391-393: Please improve the wording of this sentence, particularly: “... have a large retreat focussed on...”
L400-401: Is the model initialisation method (surface elevation change) causing this retreat to happen?
L401: It is difficult to distinguish in Fig 6a,b where the ‘east’ and ‘west’ regions are that you are describing. Is it either side of the inlet?
L409: Move the reference to Fig 4a earlier in the sentence, into line 407 after “Thwaites Ice Shelf”
L434: Is this the increase in flux around 85 years?
L438: Change “present day/future” to “present day and future”
L441: Why is there a reference to Bett et al., (2024) here? If you are referring to the similar result in that paper, you should say that
L451-452: Is it true that you “vary ice damage and ocean forcing through time” or is it just separate instantaneous changes as shown in Figure 8?
L459: Why do you start your reversibility simulations from the 700m simulation and not the 600m?
L466-467: Remove the extra “initial”
L477: This opening sentence was confusing to read at first “...with ice damage retreat...”. Consider rewording this. Maybe: “All simulations starting...ice damage at time zero retreat past...”
L477: All simulations with damage retreat past present-day except 900m?
L480: Does “undamaged” also mean “healed” because there are two different words being used here
L480-481: Remove the extra “state”
L488: CDW “depth” instead of “thickness”?
L495: You say that both coldest cases trend to state 1, but from Fig 8a, it looks like the 800m case readvances to state 2. I can see in panel 8b that you have solid circles (quasi-steady state) at state 1 but you haven’t referred to this in the text. Would it be helpful to refer to Fig 8 b,cearlier?
L499: How do you know that the three states could have been “applicable for the pre-1940s period” if you haven’t run the simulations from those states to present day?
L502: This is the first time you mention the duration of these simulations, did all reversed simulations run for 250 years?
L503: Are the “ocean only forced” simulations the same as those you’ve called “healed” in Fig 8? It’s confusing having different phrases for all these model runs.
L514: What do you mean by “within the area we are considering”?
L514: Fig 8b: How long were the simulations in 8b run for? Because the solid circles for 900 and 800 do not match the final grounded area in 8a for simulations run from 110 years
L521: Add ‘the’: “with the grounding line”
L545: “westward shift” instead of “western shift”
L558: “varying degrees of damage” - unless I have misunderstood the simulations, you only considered two states of damage: present-day “damaged” and a “healed” state. “Varying degrees” sounds like you have modelled more than two.
L562: What “threshold”? Does this refer to the 600m thermocline depth?
L565-566: Does MITgcm produce reasonable melt rates at the grounding-line?
L568: Is this needed: “with another study finding”?
L570-571: Is it only the change in ocean currents driving these higher melt rates, or is there contribution from deeper grounding lines as they retreat?
L578: Add ‘is’: “as is visible”
L579: Capitalise “figure”
L584: Were they “projections” in Bett et al. (2024) if they used idealised ocean forcing?
L589: Remove the extra “:” after the citation.
L606: “evolution” instead of “shape”
L612: Remove the out of place citation
L613: Remove the extra bracket in front of the citation
L610-615: Running simulations for only 250 years doesn’t give you enough information about tipping points or stability. I also don’t think this is “complex behaviour”: you applied warmer ocean temperatures to a weaker ice shelf which caused it to retreat, and when you reduced those temperatures and strengthened the ice it re-advanced.
L626: Remove the dash after second citation
L641: Change “damage/calving” to “damage and calving”
L676: Should “Figure 11 c,d,g” be e,f,g?
L679: What do you mean by “shown here”? Figure 11h?
L679: Do you mean an “eastward” movement?
L684: So, the damage visible in Fig 11b could be rifts caused at the grounding-line and not a shear margin?
L690: This wording is confusing: “1940s to the present day and even before the 1973 satellite image”
L697: I don’t think you can say that you’ve focussed on the “stability landscape” of Thwaites Glacier considering the short model simulations
L715: “quasi-steady states” not “steady states”
L750: This is the first time you mention other sources of buttressing on Thwaites Glacier. Is there any buttressing being provided by Bear Island or Crosson ice shelf during your different simulations, and is this likely to have a bigger or smaller impact compared to the small pinning points?
L755-756: You say “In this study we find...” that only healing the present-day ice shelf and no melt is enough to advance the grounding line, but in lines 300-301 (see above comment) you say both healing and extension are required.
L770-: The present-tense in this conclusion is a little confusing (e.g., “The next step is to run...” on line 782 sounds like future work)
L779: How do you know that it is the western pinning point grounding that allows the grounding line to advance rather than buttressing from Crosson or Bear island?
L786: Where did you show that the “acceleration led to increasing damage”? I don’t recall this.
L798: Add ‘-day’: “present-day grounding line”
Figures:
Fig 1: The yellow 2009 GL is not very visible and it probably isn’t needed anyway. Can you add ice front locations to help show the outlines of floating shelves (e.g., MacGregor et al. 2012). In the caption add a note that the distance arrow shows the width of the iceberg/previous ice shelf and is used again in fig 2.
Fig 2: Make orange flight lines more visible. Does the dashed arrow (64 km) refer to the distance between ‘a’ and ‘e’, rather than anything on the background image? If it is not too messy, add outlines of shelves in the centre panel.
Fig 3: The thick dashed lines showing zoom boxes dominate the subplots and get confused with other model outlines and features. Would the colormap in (d) and (e) be better as a red-white-blue to show the positive and negative values clearer? The caption says the ice damage is saturated to 1 and -1, does this mean there are values outside of this range?
Fig 4: Include arrows in the legend in (a). For (c) and (d) it is not immediately clear what these lines are, can you add ‘GL’ to the legend labels or title.
Fig 5: Mention in the caption that the melt rates are on a log scale. The cyan grounding line is quite hard to see.
Fig 6: In (c) and (d), put the coloured arrows and grounding lines in a legend. In the corners of c,e and d,f respectively add t=455 yrs, and t=85 years. For (c) and (d) consider using a different colormap or not having white for the thickest grounded ice and thickness below 250m. Why do you exclude thickness less than 250m? The lines in (a) and (b) might be more visible if thicker.
Fig 7: Do you need to say “Plotted over approximate area shown in Figure 3b by green dashed box” (also mentioned in other figures).
Fig 8: Remove 950, 850 etc from the x-axis. Add a small legend in 8a for the solid and dashed lines with “Healed” and “Damaged”, or “Ocean only” and “Ocean and damage” (whichever is consistent)
Fig 9: Panel (f) needs the title fixing
Fig 11: Can you put consistent label formatting on these panels