Enhanced subglacial discharge amplifies Petermann Ice Shelf melting when ocean thermal forcing saturates

Domain: earth_science
Published: 2025-05-06
Editorial outcome: major_revision

Abstract
Increased basal melting of the Petermann Ice Shelf is typically attributed to rising ocean temperatures. While subglacial discharge is known to intensify basal melt, the underlying mechanisms and their evolution in a warming climate remain unresolved. Using a 3-D numerical regional ice shelf–ocean model centered on the Petermann Fjord, we identify a regime shift in heat flux efficiency within the ice shelf cavity when discharge exceeds the current peak summer value. In this regime, thermal driving saturates, and discharge-intensified currents increase melt by enhancing shear-driven turbulent mixing across the ice shelf-ocean boundary layer. Increases in melt are most profound at the crests of basal channels, where vigorous meltwater confluence amplifies friction velocity. Challenging conventional attributions of increased ice shelf basal melting to ocean warming alone, our results demonstrate how atmospheric warming exacerbates ocean-driven melt processes and is likely to play a dominant role in amplifying future basal melt.

Introduction
The Greenland Ice Sheet is now the dominant single source of barystatic sea-level rise1,2. However, observational estimates of contemporary mass loss likely represent conservative baselines3, as process-level uncertainties in ice-ocean interactions and surface melt amplification suggest future rates will exceed current projections4,5. The northern sector of the Greenland Ice Sheet hosts several of the last remaining marine-terminating outlet glaciers featuring floating ice tongues (also called ice shelves) that exert buttressing forces to the inland ice. During 2001–2021, the dominating factors of mass loss from these ice shelves were (in decreasing order) basal melt, calving, and negative surface mass balance. Between 2003 and 2010, collapses were observed at Hagen Bræ, C.H. Ostenfeld glacier, and Zachariæ Isstrøm6. Dynamic thinning, retreat and even collapse combined with perturbations of the outlet glacier’s grounding line can lead to a loss of backstress and accelerated mass loss from the northern Greenland Ice Sheet sectors7 holding an ice volume of ca. 4m global mean sea level equivalent8,9.

Along the Nares Strait coast of the northern Greenland Ice Sheet, Petermann Ice Shelf terminates into the Petermann Fjord. Formerly the longest ice shelf in the northern hemisphere at ca. 80km in length, two recent massive episodic calving events in 2010 and 2012 reduced its length by ca. 30km10. The ice shelf draft thins from ca. 400–600m near the grounding line to ca. 50–100m near the calving front. Evenly spaced across its 15–20km width are 1–2km wide basal melt channels that incise several-100m deep into the draft11. Warming of Atlantic Water in Petermann Fjord by ca. 0.3°C has been observed since 1960, with much (ca. 0.2°C) of the warming occurring in the last two decades12,13. Average basal melt rates at the Petermann Ice Shelf grounding line, at ca. 600m below sea level, are currently reported to exceed 17m/yr, leading studies6,13 to suggest that increasing ice shelf basal melt is following increasing ocean temperatures in Greenland fjords. Numerical modelling study and satellite observations have revealed that the highest basal melt rates occur in the kilometre-size grounding zone of the glacier, driven by tidally modulated intrusions of seawater14,15. Seasonal and climate-warming driven long-term changes in regional sea ice cover are also known to enhance basal melting. According to model results16,17, elevated temperatures and an energetic ocean circulation under the ice shelf, driven by variations in Nares Strait sea ice mobility and thickness, increased basal melting by a factor ranging between 1.2 and 1.7. Also, melting can be enhanced by routing of glacier surface melt to the glacier bed, and further to the grounding line, where it enters the fjord as subglacial discharge. Modelling study suggests that a higher subglacial discharge during summer (June-August) triples basal melt at Petermann Ice Shelf compared to winter18. Furthermore, within the central basal channel, for the region extending from the grounding line to 16km seaward of it, between August and December 2015, an estimated area-averaged melt rate maximum of 170m/yr has been reported19, consistent with episodes of enhanced discharge-driven sub-ice shelf currents. This suggests that subglacial discharge may play a more important role in modulating Petermann Ice Shelf basal melt than sea ice-driven changes in ocean circulation. However, with atmospheric and oceanic warming expected to amplify in the Arctic by the end of the century20,21, a salient question remains unanswered: How would Petermann Ice Shelf basal melt respond to projected increases in discharge under a warmer future climate, and what would be the mechanisms driving this response?

In this study, we extend a high-resolution unstructured grid 3-D ocean–sea ice–ice shelf regional model setup17,22 centred over the Petermann Ice Shelf and Petermann Fjord to include subglacial discharge at the grounding line (Fig. 1). The fjord bathymetry and the ice shelf draft are derived from BedMachine v38. Modifications have been applied to the topographical datasets (discussed further below) to amend an inaccurately represented sub-ice shelf water column thickness, and furthermore, to implement a 540–610m deep inner sill situated ca. 25km from the grounding line, retrieved using aerogravity data from Operation IceBridge23 (Supplementary Fig. E1). Open boundaries in the Lincoln Sea and Baffin Bay are used to provide lateral ocean boundary conditions. These include temperature, salinity, and seawater velocities derived from A4ROMS24 (4km pan-Arctic Regional Ocean Modelling System grid), and merged sea surface elevation derived using A4ROMS and 5km AOTIM25 (Arctic Ocean Tidal Inverse Model). Data from the 5.5km polar (p) version of RACMO2.3p226 (Regional Atmospheric Climate Model) is used to compute the momentum, heat, and freshwater fluxes. Sea ice concentration and thickness, bulk ice salinity, and sea ice velocities are obtained from A4ROMS-CICE27 (Community Ice CodE). Using this setup, we perform a control run with zero subglacial discharge, and three test runs designed to simulate the response of the cavity circulation and ice shelf basal melt rates to present day subglacial discharge, increased discharge under high-emissions Representative Concentration Pathway 8.5 (RCP8.5) climate warming scenario, and an intermediate discharge scenario between the two (Table 1). We show that discharge increases basal melting by (a) accelerating a quasi-estuarine circulation28 which increases the advective heat transport into the cavity, and by (b) enhancing the shear-driven mixing across the ice shelf-ocean boundary layer. The latter exerts a dominant control over the increase in basal melting seen for higher discharge under future warming scenarios while the former (advective heat transport) saturates, suggesting a shift towards a regime where melt increases are primarily driven by increases in friction velocity (i.e., velocity scale representing shear stress in the boundary layer). Our findings challenge the prevailing focus on increased oceanic heat forcing as the primary driver of basal melting at the Petermann Ice Shelf, revealing that heightened subglacial discharge significantly enhances (channelised) basal melt. This suggests that atmospheric warming-induced excitation of oceanic melt processes will critically influence the stability of the ice shelf.

Results
Response of cavity circulation and basal melt to contemporary discharge
In the control run, Petermann Ice Shelf basal melt rates largely range from a few m/yr near the calving front to ca. 40–60m/yr under the deeper (400–600m) regions of the ice shelf base near the grounding line (Fig. 2a). Melt rates decreasing with distance from the grounding line is consistent with results obtained from steady state ice discharge divergence approach11 and satellite remote sensing29. For the period from 2011–2015, remotely sensed estimates reveal melt rates of up to ca. 80m/yr near the grounding line, although they largely range between 40–50m/yr decreasing to ca. 10m/yr under the shallower (<200m) portions of the ice shelf29, and are comparable to our modelled estimates. With the introduction of contemporary summer mean subglacial discharge (Qsg-present) into the ice shelf-fjord system, the modelled mean melt strengthens nearly throughout the domain and appears to be concentrated along the basal channels that run along the length of the ice shelf (Figs. 1d, 2a, b). Spatial variability of our modelled melt on scales similar to the basal channel configuration is consistent with observational findings11,29. Strongest increases of up to ca. 100m/yr are seen within these channels under the deeper regions, which gradually decreases towards the calving front. This increase lies within the in situ conductivity-temperature and radar-derived estimates of discharge-enhanced seasonal melt rate maxima of 40–170m/yr19. The summer mean basal melt rate (mb), friction velocity (u*), and thermal driving (ΔT) are averaged over (a) the entire ice shelf base, and (b) over parts of the ice shelf where draft is deeper than 300m, and are indicated, respectively, by subscripts pgis and 300m appended to the variable in question. Compared to the summer mean modelled melt rates in the control experiment, the contemporary discharge-driven increase in mbpgis and mb300m, respectively, is 41.2% and 68.1%, driven by similar increases in u*pgis (47%) and u*300m (67.4%) and modest increases in ΔTpgis (1.1%) and ΔT300m (4%) (Fig. 1d).

The presence of subglacial discharge at the grounding line enhances the fjord-scale overturning circulation (Fig. 3i). We observe a more robust Atlantic Water inflow at depth, extending into the inner fjord basin and up to the grounding line (Fig. 3c, f). Hereafter, we refer to this process as the thermal pump, which promotes the entrainment of heat and salt deep into the ice shelf cavity (columns 1 and 2 in Fig. 3). We note an increase in bottom temperature of 0.03°C at location ca. 10km from the grounding line, and a temperature excess (thermal driving) of 0.1°C under the ice shelf base at 300m depth (Figs. 2j, 3g). A higher thermal driving near the grounding line (Fig. 2i, j) drives more melting at depth. A more energetic outflow/strengthened meltwater plume dynamics, driven by the introduction of discharge, exacerbates the shear-driven turbulent mixing and heat transfer across the ice shelf-ocean interface (Fig. 2e, f), thus driving more melt. The spatial variability of friction velocity anomalies is consistent with the corresponding melt rate anomalies (cf. panels b and f in Fig. 2). Regions exhibiting a strong increase in friction velocity are concentrated along the longitudinal basal channels near the grounding line, and gradually decrease downstream of it. However, a higher meltwater production (Fig. 2b) and stronger outflow (Fig. 3f, i) deposits larger volumes of meltwater from the deeper to the shallower regions of the ice shelf base, considerably lowering the thermal driving locally (Figs. 2j, 3d, e). Thus, while the increase in u*pgis and u*300m is significant, the contribution of thermal driving is largely restricted to the deeper regions of the ice shelf draft (ΔT300m) (Fig. 1d). This implies that for subglacial discharge equivalent to present-day summer mean magnitudes, increase in basal melt under the deeper (>300m) drafts is driven by increases in both thermal driving (increased thermal pumping) and friction velocity (increased shear-driven mixing).

Heightened discharge and saturation of the thermal pump
When transitioning from the control experiment to heightened discharge in the RCP 8.5 scenario (Qsg-RCP 8.5), and an intermediate discharge scenario between the two (Qsg-median) (Table 1), a higher modelled discharge, as expected in a future warmer atmosphere30,31, results in considerably higher increases in basal melt as compared to Qsg-present (Figs. 1d, 2a–d). In particular, an increase in discharge from control to RCP 8.5 results in more than twofold (168.4%) and threefold (217.2%) increase in mbpgis and mb300m, respectively. These increases are primarily driven by increases in u*pgis (191.8%) and u*300m (216.9%), whereas ΔTpgis and ΔT300m show relatively nominal increases of 2.9% and 6.5%, respectively. For the median experiment, 2.3% and 5.9% increase in ΔTpgis and ΔT300m, and 132.5% and 163.4% increase in u*pgis and u*300m generate 116.3% and 163.6% increase in mbpgis and mb300m, respectively. As noted for Qsg-present, the strongest increases occur under the deeper regions of the ice shelf base, which gradually decreases downstream (Fig. 2c, d), however, channelised melting is seen to strengthen substantially (discussed further below).

Combining insights from all the experiments provides a novel perspective on how discharge dictates the interplay between Petermann Ice Shelf basal melt and its drivers: Increase in the heat entrained in the ice shelf cavity as a result of the thermal pump is strongest at the inclusion of present-day summer mean estimate of subglacial discharge, with much of the heat delivered by the thermal pump concentrated under the deeper ice shelf base in the inner fjord basin (Figs. 1d, 3d, g). At a summer mean discharge value which is approximately fourfold higher than the present (i.e., Qsg-median vs. Qsg-present), thermal driving increases at a lower rate (Figs. 1d, 3g; cf. Fig. 3d and Supplementary Fig. C1a). We note an increase of 1.9% in Qsg-median ΔT300m as compared to Qsg-present (Fig. 1d). Moreover, a larger meltwater deposition further lowers the thermal driving under the shallower drafts (Figs. 2b, c, j, k, 3i and Supplementary Fig. C1a–c), which results in a relatively weaker increase in ΔTpgis (Fig. 1d). Lastly, in the RCP 8.5 scenario, discharge is nearly doubled as compared to the Qsg-median experiment, however, the efficiency of the thermal pump in increasing basal melting saturates (Fig. 3g and Supplementary Fig. C1). Here, increases in ΔTpgis and ΔT300m are indiscernible as compared to the Qsg-median experiment (Fig. 1d). Instead, we find that discharge under the high-emissions RCP 8.5 scenario enhances melting by further accelerating the cavity circulation, wherein, 24.1% and 20.4% increase in mbpgis and mb300m are seen, driven predominantly by 25.5% and 20.3% increase in u*pgis and u*300m, respectively.

In sum, our experiments reveal a regime change in heat flux efficiency within the ice shelf cavity. In the control experiment, the ice shelf-ocean boundary exchange consumes a considerable amount of heat from the ocean that is available for melting the ice. As such, the latent heat loss and meltwater release cools the ambient ocean inside the cavity. In such a scenario, an increased fjord-scale overturning circulation due to increased discharge (e.g., ≤ Qsg-median experiment) entrains more heat in the cavity (thermal pump), which increases melting through enhanced thermal driving, in addition to the enhanced shear-driven turbulent heat exchange. This saturates at some point, when more discharge (e.g., > Qsg-median experiment, as seen for the RCP 8.5 scenario) accelerates the overturning circulation such that there is always enough heat coming into the cavity to maintain a high thermal driving (i.e., heat supply at an elevated level is maintained). In this case, melt rates become entirely a function of the friction velocity.

Convergence of meltwater and channel stability
We use surface layer salinity as a proxy for tracking the spatial distribution of meltwater under the deeper regions of the ice shelf (Fig. 4 and Supplementary Fig. D1). Freshening occurs nearly throughout the domain in the discharge-driven experiments as compared to the control run. Furthermore, regions of substantial freshening align with four prominent ice shelf basal channels (Figs. 4a, b, 5a, b and Supplementary Fig. B1) along which the draft thins from 300–450m ca. 5km from the grounding line to 150m ca. 25–30km seaward of it, suggesting a topographically steered horizontal pathway for the meltwater outflow that is determined by the ice shelf basal morphology and which in turn influences it (discussed further below). Changes appear to be most profound within these channels, particularly at the crests, where salinity decreases by up to 1.06 psu in the RCP 8.5 experiment with respect to the control run (Figs. 4, 5a, b). These locations of meltwater convergence are concordant with locations of the largest melt rate increases (discussed further below) of up to 181.5m/yr; a sevenfold increase compared to the control run. Thermal driving is not higher at the crests where meltwater converges, but under the deeper drafts, which is consistent with the pressure-dependent melting point of seawater (Fig. 5c). Although thermal driving may not increase the most at the crests, we observe increases ranging from 0.02 to 0.12°C in the RCP 8.5 experiment compared to the control run across the transverse cross-section near the grounding line (Fig. 5c). We find that the friction velocity profile across this section is (nearly) congruent with the basal melt rate profile, and exhibits highest increases at the crests of up to 0.03m/s (Fig. 5d).

Three ice shelf basal channel sections are outlined (Fig. 6a), whose lengths are demarcated based on the apex of the 150m ice shelf draft contour (e.g., Fig. 4a and Table 2). The channel averaged annual mean basal melt rates for the western (Sw), central (Sc), and eastern (Se) sections for the RCP 8.5 experiment are calculated as 45m/yr, 51m/yr, and 47m/yr, respectively. Enhanced melting at the crests will promote channel opening, and we can calculate the changes in channel averaged annual mean ice volume above each channel over time to estimate the incision/erosion timescale for that channel under the respective basal melt rates (Eqn. (5), (6)). To bracket the role of ice advection from upstream of the grounding line, three different rates (AR) are considered for each channel which correspond to the contemporary mean10 (1km/yr), and 50% lower (0.5km/yr) and higher (1.5km/yr) than contemporary mean values. While Sw exhibits the lowest mean basal melt rate, it also features the shallowest mean draft value (Table 2). Notably, its incision timescales are comparable to those of Sc, which has the highest mean melt rate but a relatively deeper mean draft compared to Sw. The number of years (N) required for the mean ice volume above the channels to become zero is estimated to be 5.88, 6.52, and 7.31 years for Sw, in contrast to 5.95, 6.75, and 7.8 years for Sc, at AR of 0.5km/yr, 1km/yr, and 1.5km/yr, respectively (Fig. 6b). It has been shown that the viscous ice response to a basal channel is insufficient to significantly delay the breakthrough of basal channels in ice shelves that are thinner than 400m32. With a lower mean melt rate than Sc and the deepest mean draft value, Se demonstrates the greatest resilience. We estimate N to be 8.06, 9.41, and 11.31 years for AR of 0.5km/yr, 1km/yr, and 1.5km/yr, respectively (Fig. 6b). These estimates suggest that ice advection from upstream of the grounding line (alone) may not be sufficient to compensate for the thinning of the ice shelf due to intense channelised basal melting that may be expected under the RCP 8.5 scenario (discussed further below).

Methods
The Finite Volume Community Ocean Model (FVCOM) setup
The unstructured grid, free-surface, 3-D primitive equation Finite-Volume Community Ocean Model (FVCOM)61, augmented by modules for an ice shelf62 and sea ice22, is used. The model domain is centred on Petermann Fjord, covering the region between 75°N–87°N and 29°W–81°W (Fig. 1a). The grid comprises of unstructured triangle elements, featuring a total of 112709 nodes and 219836 cells, including 309 open boundary nodes in the Lincoln Sea and Baffin Bay. The topological flexibility of the unstructured grid allows accurate fitting of the complex irregular coastal geometry along the northern Greenland Ice Sheet margin. In particular, the grid is partitioned into polygons of varying horizontal resolution: A 200m resolution is provided in the Petermann Fjord region to resolve the steep slopes characteristic of the fjord seafloor and ice shelf basal topography, which is relaxed to 2km in the adjacent areas outside the fjord, and further to 4km toward the open ocean boundary regions. For an extensive regional model setup, such a configuration enables us to strike a balance between computational efficiency and the need for accurate representation of topographical features and their variation over small spatial scales.

Petermann Ice Shelf draft derived from the BedMachine v3 dataset8 (150m resolution) represents a pre-2010 calving geometry10, and is interpolated to the 200m FVCOM model grid22. As such, the basal channels are resolved at 200m scale in our model, which is sufficient to represent their 1-2km width11. Bathymetry outside Petermann Fjord comes from the 1km gridded IBCAO v3.0 dataset63, whereas inside it, the BedMachine v3 dataset is used. In the vertical, the domain between the irregular seafloor topography and the ice shelf draft (or the open ocean surface) is discretized into 23 terrain-following σ-layers. The ice shelf module dynamically adjusts horizontal pressure gradient forces to account for the weight of a freely floating ice shelf, with terrain-following σ-layers adapted to subduct beneath the prescribed ice shelf draft62. To alleviate errors in the discretization of horizontal pressure gradients along the sloping σ-layers, smoothing of the bathymetry is applied as a function of the mesh resolution before its implementation in the model22. Without seafloor depth measurements beneath the ice shelf, the BedMachine bathymetry is inferred by extending the known seafloor depth from seaward of the calving front landward towards the grounding line and aligning it with the depth of the ice shelf draft at the grounding line using natural neighbour interpolation. For Petermann Fjord, this leads to implausible abrupt changes in the water column thickness22. Specifically, the BedMachine bathymetry underestimates the water column thickness under the deeper (>200m) regions of the ice shelf, with sub-100m values along the flanks and in the central zone (Supplementary Fig. E1a–c). More critically, this misrepresentation results in the formation of an artificial inner sill, which, if not accounted for, would block the inflow of the warm and dense (27.9 kg/m3) Atlantic Water at depth (Supplementary Fig. E1d) that contacts the grounding line12,17. Therefore, addressing these artefacts is imperative to ensuring that the model is able to reproduce the fundamental aspects of glacial fjord circulation28. Aerogravity data from Operation IceBridge was utilised to model the bathymetry beneath the ice shelf, revealing the presence of an inner sill, ca. 540–610m deep, located about 25km from the grounding line23. In this study, we implement a smoothed bathymetry product, wherein the BedMachine bathymetry has been modified to account for the unrealistic sub-ice shelf water column thickness and further integrated with the aerogravity data22.

Experiment design
The model is run from July 01, 2014 00:00:00 UTC – January 01, 2017 00:00:00 UTC with hourly oceanic (temperature, salinity, velocities, and sea surface elevation), 3-hourly atmospheric (heat and freshwater fluxes, and wind speeds), and daily averaged sea ice (sea ice concentration, thickness, and velocities, and bulk ice salinity) forcings obtained from respective high-resolution regional models over this period17,22, and without discharge, which represents the Control experiment (Table 1). To allow investigations of the response of ice shelf basal melt to variations in subglacial discharge, our FVCOM setup has been augmented to include discharge across the ca. 20km wide grounding line. Three discharge experiments are initialised from the stable model solution of the Control experiment from January 01, 2016 00:00:00 UTC and run until January 01, 2017 00:00:00 UTC (Table 1). Lacking information regarding the precise location(s) of emergence of subglacial meltwater across the grounding line, discharge is injected uniformly across the 90 cells (hereafter denoted by Δ) representing the grounding line (Fig. 1b). At each Δ, discharge Qsg is vertically uniform across the 23 σ-layers. The salinity of Qsg at each Δ is set to 0, and the temperature is set to the (in situ) pressure-dependent freezing point of freshwater. The June-July-August (JJA)-mean discharge values presented in Table 1 represent the daily discharge values over the summer (JJA) months used in our model. Outside of this period, discharge is set to zero (similar to ref. 18). Results for the Control and discharge experiments are presented from the final year of the simulation (January 01, 2016 – January 01, 2017). We note that the pre-2010 ice shelf geometry provided by BedMachine v3 is not consistent with the forcing period(s). Change in ice shelf geometry (e.g., due to calving) has been shown to impact fjord circulation and basal melt rates64. However, the two prior large calving events in 2010 and 201210 removed the softer and thinner (<150m) outer sections of the ice shelf7. As such, we posit that the findings of our study-discharge-driven regime change in heat flux efficiency in the ice shelf cavity, and amplified channelised melting under the deeper drafts - are robust, irrespective of the geometry used.

The present-day JJA-mean discharge (Qsg-present) and its increase in a future warming climate under the RCP 8.5 scenario (Qsg-RCP 8.5) for the Petermann catchment is obtained from simulations using a regional model43 (Table 1). Here, a subset of the Coupled Model Intercomparison Project (CMIP5) Atmosphere and Ocean General Circulation Models (AOGCMs)65 (Fig. 1c) is used to force the Modéle Atmosphérique Régional (MAR) 3.9.666 for the period 1950–2100. The projected MAR meltwater runoff over the Petermann catchment is then bias corrected to ensure that it is consistent with the best present-day (1995–2014) runoff estimate, and to also enable a smooth transition from present to RCP 8.5 forcing. It is then summed up over the Petermann catchment to give the discharge estimates. We note that the modelled runoff outputs were stored as annual means and were converted to the JJA-mean runoff by multiplying the output by a factor of 365/92, where 365 = number of days in a year, and 92 = number of days in JJA. This coarse conversion introduces a minor discrepancy between the converted and true modelled JJA-mean values; however, we maintain that this does not significantly impact the findings of our study. While inter-model differences can be large (Fig. 1c), all five CMIP5 AOGCMs were used to prepare the multi-model mean discharge used in this study. Note that the median discharge magnitude (Qsg-median) is not obtained from the model43, but is constructed using the Qsg-present and Qsg-RCP 8.5 magnitudes (as the median of Qsg-present and Qsg-RCP 8.5). As such, it represents an intermediate future mean discharge magnitude, however, it does not correspond to any known projection scenario or a particular time period (Table 1).

Diagnostics of basal melt rate and its drivers
For the investigation of ice shelf basal melt rate, its expression in the thermodynamic framework62,67 of the primitive equation setting is recalled:

ρfw mb L = -ρsw ΓT u*ΔT, (1)

Here, ρfw and ρsw are the densities of freshwater and ocean water, respectively. L is the latent heat of fusion of ice, cw is the specific heat capacity of ocean water, and ΓT is a non-dimensional heat-transfer coefficient. mb is basal melt rate, controlled primarily by friction velocity (u*) and thermal driving (ΔT), and where

u* = CD (u2w + u2res)1/2, ΔT = Tb - Tw, (2)

Tb = λ1Sb + λ2 + λ3 Pb, (3)

ρfw mb Sb = -ρsw ΓS u*(Sb - Sw), (4)

Above, CD is the drag coefficient, uw is the velocity magnitude some distance away from the ice shelf-ocean boundary, and ures is the velocity magnitude of small sub-grid scale residual currents68. Tb, Sb, and Pb are the potential temperature, salinity, and pressure at the ice shelf-ocean boundary, and Tw and Sw are the potential temperature and salinity some distance away from it. λ1, λ2, and λ3, respectively, are the slope, intercept, and pressure coefficient of the liquidus, and ΓS is a non-dimensional salt-transfer coefficient. Note that Tw, Sw, and uw are obtained from the first (i.e., uppermost) σ-layer. Furthermore, ΓT and ΓS are application-specific, and are here in tuned according to the applied boundary conditions, Petermann Ice Shelf basal topography, and mixing schemes so as to be consistent with the contemporary observational estimates of the Petermann Ice Shelf basal melt (Supplementary Table A1). The drivers of basal melt rate (mb), namely, friction velocity (u*) and thermal driving (ΔT) are used to investigate the impact of discharge on the summer mean basal melt rate.

Changes in ice volume above basal channels over time
We calculate the changes in the annual mean channel-averaged ice volume above basal channels over time (V(N)) to determine the number of years (N) it would take for V(N) to become 0 as

0 = Vi + (N × Va) - (N × Vm), (5)

N = Vi / (Vm - Va) (6)

where Vi is the initial (contemporary) ice volume above a channel expressed as Lc (channel length [m]) × Wc (mean channel width [m]) × Hc (mean channel depth [m]). Va is the annual mean volume of ice advected downstream of the grounding line per year, expressed as AR (mean advection rate [m/yr]) × Wc × Hc. Vm is the annual mean volume of ice melted per year, expressed as mbc (annual mean channel-averaged basal melt rate [m/yr]) × Lc × Wc. Dimensions associated with each channel section are detailed in Table 2. Note that Wc drops out of the calculation, and as such, is not detailed in Table 2.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study the authors investigate the impact of changes in subglacial discharge on submarine melting of the Peterman Glacier ice shelf. Submarine melting is assumed to scale with thermal forcing and water velocity. Increases in subglacial discharge cause the thermal forcing to increase by drawing warm water into the sub-ice shelf cavity, but only to a point. At high values of subglacial discharge the thermal forcing saturates and melting is primarily affected by water velocity. This is an interesting result that has important implications for destabilization of the marine-terminating glaciers and ice sheets. The oceanographic modeling and results are solid, but I am concerned about the implications for ice shelf stability. The ice shelf is not modeled (unless I am missing something), which means that the channels will only grow with time. Channels advect downstream and are eventually calved, and additionally I would expect changes in stress associated with channel enlargement to cause ice to flow toward the channels and at least partially offset melting. Neither of these things are included in the model, which causes me to question to what extent future runoff would destabilize the glacier. I like this paper, although I do also think that it could benefit from some revisions to improve readability. I understand that Nature Communications wants/requires the methods to come at the end of the paper, but it felt to me that the methods was initially written as section 2 (before the results) and only later moved to the end of the paper, with the consequence that the results were difficult to interpret. For example, it took me awhile to understand what was meant by the various model scenarios. I think the introduction would benefit from a few sentences, or perhaps another paragraph, that would more fully describe the model set-up and simulations in layman’s terms. Finally, I found the use of acronyms and variables within text a bit cumbersome. I encourage the authors to eliminate them to the extent possible in order to make the paper more accessible. “basal melting”, “subglacial discharge”, “friction velocity”, and “grounding line” are much easier to read than m_b, Q_sg, u*, and GL. I don’t see why the glacier and fjord names need to be abbreviated, and some acronyms are not used very often (such as GrIS and RG). There are many places where the text can be tightened up if space is an issue. Specific comments: Title: I feel that the title is a little misleading, although I agree with the sentiment, as this study is not actually modeling glacier flow or fracture propagation. L20: How do “channels carve incisions”? Aren’t the channels themselves products of melting? L21: Are the thickness values stated for the areas immediately above the channels, or is this a statement about the average thickness of the ice shelf? Also, I assume that at least some of the thickness gradient is due to longitudinal extension, but this reads as though it is entirely due to submarine melting of channels. L54: Please provide a definition of friction velocity here. L113: What do you mean by “peaks”? Fig. 3 shows that the temperature continues to increase as you increase the subglacial discharge, just at a lower rate. L117-L118: Similar to previous comment. Contribution to what? And do you just mean that the thermal forcing increases at a lower rate? Fig. 3: What does “27.9” refer to? L213-216: I’m not sure that I follow this. Are you essentially saying that the addition of cold freshwater will limit the thermal pump? L370-373: How sensitive is the model to variability in subglacial discharge? Your simulated discharge follows a boxcar function, while we know that subglacial discharge follows more of a gaussian-like pattern during summer and is highly variable. I found myself wondering how important the variability is. If you have a few days or weeks of exceptionally high or low subglacial discharge, can that push the model into a different regime? L419-430 (and Fig. 6): I found this description confusing. It’s not clear to me what exactly you are calculating. Are you calculating the volume of ice above a channel? Maybe I am just hung up on the term “channel ice volume”, when to me the word channel implies a place where there is no ice. Table 1: Is Q_sg-median just the average the present and RCP 8.5 scenarios? Or is that just coincidence?", "reviewer_2": "The impact of subglacial water on ice shelf melt rates has been an understudied topic, but is now gaining prominence in the field with increasing recognition of its impact on ice sheet dynamics. This paper presents interesting insights into the mechanisms by which enhanced subglacial water flow might affect melt rates (i.e. shear-driven turbulence vs. thermal driving) and will therefore be of interest to the research community. The study is well designed, using a very neat set of experiments to draw interesting conclusions, and the paper is well written: clear, and easy to follow. The figures are also well chosen, and very convincingly underpin the conclusions in the results and discussion sections, although they are not as polished and easy to read as I would like. I make some further suggestions for changes to the text below, but they are only minor. Overall, I this this is a very high quality contribution to the field. Minor comments Lines 3-4: “contemporary mass loss rates are underestimated”. Underestimated by whom? Observational studies, or ice sheet models? Line 14: “ice volume of ca. 4m global mean sea level equivalent” I can’t understand from this sentence which sectors of the ice sheet you are referring to. Please check the reference too – is Morlighem et al. (2017) not more appropriate here (BedMachine, see reference list below)? Line 16: “Erstwhile” is not wrong here, but it’s a rarely used word. You could consider replacing with “Formerly” or “Previously” to aid non-native speakers. Line 21: “thinning it from ca. 600 m near the GL to ca. 100 m near the calving front.” Phrased this way, the sentence can be misinterpreted… do you mean “with ice thickness over the channels of 600 m near the GL to ca. 100 m near the calving front”? Figure 1b: I suggest labelling the calving front with a date. Line 34-36: “estimated area-averaged melt rate maximum of up to 170 m/yr” the phrase “area-average” initially implies that this was estimated over a large sector of the shelf. I suggest rephrasing slightly to make clear that this estimate was calculated only within the basal channel. Line 38: “with warming expected to amplify” I suggest clarifying whether you mean atmospheric warming, ocean warming, or both. Final paragraph of introduction (lines 43-59): It’s unusual to see the results summarised so thoroughly at the end of the introduction. It might be more helpful to give a very brief summary of your experiment design instead e.g. from line 49 onwards something along the lines of “We performed a control run with zero Qsg, and three test simulations designed to simulate present day subglacial conditions, increased Qsg under high-end scenario climate warming, and an intermediate Qsg scenario between the two.” Paragraph beginning line 62: This paragraph jumps into discussion of modelled melt rates without describing the results clearly first, which means the reader has to work hard to catch up. I suggest an opening sentence or two describing the melt pattern in the control and preset day runs, before going on to compare them to observations. Figure 2: You could make the model results larger and easier to see by removing unnecessarily duplicated color scales and axis labels (see also general comments on figures below). Figure 3: The isopycnal label (27.9) is large and distracting (see also general comments on figures below). Figure 5: The lines here can be a little hard to see – consider changing the color scheme to improve readability and increasing line thickness (see also general comments on figures below). Line 151 & Figure 6: The two sets of location names are a bit confusing here (Sw, Sc, and Se vs. Lw, Lc, and Le). Are they both needed? Lw, Lc, and Le are barely mentioned in the text, and could perhaps be removed? Figure 6: The pink dots in panel 6a are not described in the caption, and the symbols marking different lines in panel 6b are far too small to be readable (see also general comments on figures below) Line 222: Should read “buoyant meltwater rises” Discussion line 219-240 & Figure 6b: Given the very straight erosion trajectory in figure 6b, I initially thought you must have applied the summer value for Qsg continuously over the whole year. Looking at the methods it seems that Qsg does go to zero outside of the summer, so why don’t we see a seasonal cycle in the channel volume? Presumably over winter the melt substantially decreases but the ice advection should change much less? I think the discussion section could do with a bit more context on seasonality, including how seasonally varying ocean temperatures affect erosion rates. Line 266-267: There has been further recent work on this subject in Antarctica beyond your references, so I think the statement “insights into this process are limited to tacit reasoning” is a little unfair. I suggest reading and citing: Goldberg et al. 2023, Gwyther et al. 2023, and Pelle et al. 2024. All figures: Please take a look at the text within your figures and make sure it meets minimum font size guidelines, but is not excessively large and distracting (as with some subpanel labels). Some axis labels and legends are very small indeed. In general, using consistent font across all figure will also improve perception of your figure quality. You could also improve clarity of figures 2, 3 & 4 by indicating on colour scales where a variable is absolute or relative (e.g. “Salinity” in panel 4a and “Change in salinity” in panel 4b, use a delta symbol if space is tight).", "reviewer_3": "Petermann glacier in NW Greenland is rapidly losing mass and contributing strongly to modern global sea-level rise. Ice flow from Petermann is currently slowed by the presence of a ~30km floating ice tongue which buttresses flow. With ocean warming across the last ~two decades, that ice shelf has thinned, and future breakup of the ice shelf is possible. This paper uses numerical modeling to estimate the role that future atmospheric warming may have on ice shelf break-up. Where atmospheric temperatures warm, ice surface melt increases. This surface melt is able to reach the subglacial environment at moulins, and then flows into fjords, where it influences ocean circulation and melting. Using a coupled ice-shelf/ocean/sea ice model, this work finds that subglacial discharge causes enhanced subglacial melt, particularly within subglacial channels. Enhanced subglacial discharge acts to accelerate melting by Generally speaking, this work is timely, impactful, and worthy of publication in Nature Communications with revisions. The majority of my suggested improvements below are easy to implement or require simple clarification. I have three more major concerns, one of which is easily addressed (point two below) by additional discussion or re-framing of the work. The other two points may require additional modeling work. My primary issue with the work is that it frames itself as an improvement on previous modeling work because it uses more realistic domain geometry. However, the domain geometry used in these model runs is not sufficiently explained. Maybe I missed it, but I cannot find where the authors got ice shelf draft geometry from. I assume it is a radar flight line? From what year and what instrument? More importantly, the ‘improved’ seafloor geometry is a modified form of BedMachine, but the improved seafloor geometry is never shown. What is the cavity shape? Why were the modifications made to BedMachine? It’s also worth noting that — to the best of my knowledge — BedMachine has no data for the ice shelf cavity. It simply interpolates/kreigs between radar flightlines (and the mass-conservation estimate between them) for grounded ice and the downstream seabed bathymetry from ships. So really very little is known about the actual shape of the seafloor for the ocean cavity that is modeled here. To what degree does uncertainty in this cavity geometry impact the modeling results presented here? I recognize that a full suite of modeling runs across a range of seafloor geometries is likely beyond the scope of this work, but even a single run that presents a modified geometry would be useful for determining the sensitivity of the model results to the unknown cavity geometry. The ice geometry used in the model runs was a pre-2010 ice geometry (pre large calving event) but the atmospheric and ocean forcing was for later years (I believe). Given the large change at the glacier in 2010/12 with the calving, what uncertainty is baked in to modeling results as a product of not using the modern ice configuration? Why was this choice made? Some discussion is warranted here. A third concern I have is a lack of discussion of the role of tidal pumping within the grounding zone environment and its potential impact on grounding zone melt rates. Recent work shows that the high levels of subglacial melt at the Peterman grounding zone might be due to seawater movement with the tides (Gadi et al (2023, GRL) and Ciraci et al (2023, PNAS)). The sub-ice shelf community has recently been very focused on the role of tidal seawater pumping and potential melting inland of the grounding zone by this seawater. There is a rich literature describing both the theory (Walker et al., 2013), modeling seawater intrusions (Robel et al., 2022), and observations of intrusions (Rignot et al., 2024). For a nice overview and intro to this subject, Parizek 2024 (GRL) summarizes the current state of thinking of tidal pumping and its role at Peterman. I think it is a significant omission that this work is not considered in this paper, as tidal pumping across the grounding zone may be a large source of melt (and water flow) there. Again, I recognize that a full suite of model runs that incorporate tidal pumping across the grounding zone are beyond the scope of the work presented here, but from my prospective, not even mentioning the role of tidal pumping across the grounding zone is a significant omission of this work: if nothing else, the potential role of this process should be discussed in the introduction and discussion. Line-by-line and more minor comments The abstract leaves me confused about geometry: a mechanism is described that exacerbates melt at the crest of channels (which I assume is the highest elevation point where the ice is the thinnest) but then the next sentence says that this melt will complete erode channels, which sounds like a flattening of channels. Does ‘erode channels’ here mean a complete incision of the shelf? If so, I would change language to ‘We estimate that sustained intensified summer melting over a decade may undermining the stability of PGIS by allowing channels to cut entirely through the ice shelf’ I found the reliance on terms (e.g. Qsg for subglacial discharge) to be occasionally difficult — for the casual readier, this paper presents a significant number of terms that must be held in memory to derive meaning. I might encourage the authors to be more selective in what they make terms for (is it a word-count problem to just say ‘subglacial discharge?’), but I recognize that this is a style choice. 8-11 the two clauses of this sentence appear to be jarringly unrelated, and generally I have a hard time following the logic structure of 8-14 13 Does the Petermann shelf do much in the way of buttressing? I thought that there was minimal grounded ice velocity response to the 2010 and 2012 calving events, which would suggest perhaps not. There may be modeling work that has been done that quantifies this. If so, worth citing and quickly mentioning here (given that the ice shelf community has recently been surprised by the lack of buttressing at other ice shelves like eastern Thwaites). 21 thinning the ice shelf from.. 21 is the 600m/100m geometry presented for *within* the basal channels or regions outside of the channels? 35-36 slightly confused by geometry here, is this an integrated melt rate for the region from GL to 16km downstream? What time of year were those observations made (i.e. sentence before says that melt rates triple in summer). 37 My sense is that this is a controversial statement, so I would soften indicates to suggests 45-46 reword to improved sub-ice shelf bathymetry using BedMachine3 (Morlighem et al., 2017). As written it reads as if Morlighem is a prior ice-ocean model that these results will be improving upon. 53 does estuarine circulation start (i.e. at present time) NOT saturated and then becomes saturated at a point in time (suggesting regime change at that point?) Clarify here. I also don’t have a great sense of what saturated circulation means. 58-59 at risk of repeating myself, I still wish for a different term than ‘complete erosion of the channels’ — this makes it sound like channels are going away. What about ‘complete erosion of the ice shelf within the channels’? 62-64 The numbers for each of these values (your melt, remotely sensed Melt, in-situ measured melt) really need to be shown here… ‘they agree with remote sensing/field data’ is not sufficient - give a detailed comparison. 95 One more sentence of summary/overview would be nice here - a less technical ‘this is what this means’ 97-129 Missing from this section is a physical description of the mechanism responsible for the saturation of the thermal pump. Is it increased (cold) meltwater pooling at locations in the cavity? 141-142 This is not intuitive for me: isn’t meltwater quite fresh and cold (ie minimal thermal capacity, minimal melt)? Why is there highest melt in areas of meltwater trapping, I would expect melt to be *low* there because of cold temps and increased stratification. Explain. 147-148 I would expect that meltwater would pool at the base of channels and melt would be slowest at the channel apex. Why is this not the case? I believe this is what ROV submersible Icefin observed under Thwaites (e.g. Schmidt et al., 2023 and the subsequent work of Peter Washam) - high melt rates on channel sides and slow melt at channel apex. 164 Given that old (pre-2010) ice draft geometry is used here with modern (2024) melt rates, are these numbers realistic? How would they change if modern ice draft geometry was used? 181 Reword: The modeling work of Hill 188 I would suggest adding one more sentence here about WHY this has a significant bearing on PG stability (‘this is significant’ feels like a bit of a punt without more context — however I frequently do it in my writing as well!). 189-196 THIS! Cool. I would consider adding these ideas to the abstract and introduction. 288-292 This is a key insight (no edits, just applauding) 299-302 I’d push back on this — I see how improving our knowledge of grounding zone geometries is being prioritized by funding agencies and research groups. Your study highlights the need for such work, but saying that this work isn’t being undertaken is an overstatement. Many good groups are doing targeted geophysical work to characterized the grounding zone environment, there are just major logistical and financial barriers to doing this on a continent-wide scale. Maybe re-cast this argument as a push for even more targeted airborne and ground-based geophysical campaigns to describe the grounding zone environment where possible. 304 worth throwing in a citation for BedMap3 here too 309 The argument that you’re making — that meter-scale features are essential to capture in content-wide datasets — is admirable, but we’re just so far from that reality in this moment. More useful, I wonder, might be to suggest ways to parameterize the meter-scale effects to make our continent-scale datasets more usable at this time? Realistically we are not going to send Icefin, Autosub, or Ran under every ice shelf to map every feature. 312 Worth noting though that modeling work has shown that distributed Qsg produces more subglacial melt than channelized Qsg (Slater et al., 2015 GRL) 337 is this seafloor geometry of ice base geometry? 340 include a note as to whether or not this geometry is reasonable, with a citation… I.e. how wide are the channels at various points across the ice shelf— is a 200m grid sufficient to capture their geometry? Figure 1 To improve comprehension without needing to read long legend, I would suggest adding text to 1d that signifies that dark regions are deep water column and lighter regions are integrated Figure 5 What is the source of data for the ice draft?", "reviewer_1_attached": "In this study the authors investigate the impact of changes in subglacial discharge on submarine melting of the Peterman Glacier ice shelf. Submarine melting is assumed to scale with thermal forcing and water velocity. Increases in subglacial discharge cause the thermal forcing to increase by drawing warm water into the sub-ice shelf cavity, but only to a point. At high values of subglacial discharge the thermal forcing saturates and melting is primarily affected by water velocity. This is an interesting result that has important implications for destabilization of the marine-terminating glaciers and ice sheets."}, "author_response": "We appreciate your time and effort in reviewing our manuscript. We are glad to hear that you find our work to be important and of value to the scientific community. We also appreciate the constructive feedback and have revised our manuscript following your suggestions. Please note that line numbers provided below refer to the track changes document (appended to this report for convenience), where new text appears in blue and removed text is marked with red strikethrough. Firstly, we would like to take this opportunity to inform you that during the revision process, we noticed a (minor) mistake in our calculation of the drivers of melt, namely, thermal driving and friction velocity. While the updated numbers are slightly different, we note that this does not impact any of the findings of our study, and the revised figures remain qualitatively identical. We humbly request to proceed with these corrected values, and sincerely apologize for the oversight and for any inconvenience this may cause. Please see the revised percentages in Lines 133 – 134, 167 – 170, 184, 195, and 227, and Figures 1(d), 2(e – l), and 5(c,d). Please find our detailed (point-by-point) responses below. The oceanographic modeling and results are solid, but I am concerned about the implications for ice shelf stability. The ice shelf is not modeled (unless I am missing something), which means that the channels will only grow with time. Channels advect downstream and are eventually calved, and additionally I would expect changes in stress associated with channel enlargement to cause ice to flow toward the channels and at least partially offset melting. Neither of these things are included in the model, which causes me to question to what extent future runoff would destabilize the glacier. Thank you for raising this concern. You are right, we have a static ice shelf geometry (now clarified in the revised text). Indeed, ideally there would be downstream advection and calving, and secondary ice-flow which would offset the channel growth. Mean advection rate of ca. 1 km/yr has been reported for Petermann with a seasonal modulation of ca. 0.1 km/yr (Muenchow et al., 2012). However, to our knowledge, its evolution in a future warming climate has not yet been investigated. Therefore, in our study, we envelop the contemporary observational estimates from Meunchow et al. (2012) with ± 0.5 km/yr. Ice-flow modelling study from Wearing et al. (2021) has shown that secondary flow is insufficient to significantly delay the breakthrough of basal channels in ice shelves that are thinner than 400 m. We acknowledge that the Petermann ice shelf draft is deeper than 400 m in certain sectors (e.g., southeast sector near the grounding line). However, with the mean-state analysis (mean advection and basal melt rates) that was feasible to implement here, and without a high- resolution (at the scale of the topographical data) ice-flow model to capture the nuanced dynamical variability, it is difficult to ascertain how localized small-scale channel closures/openings would impact the wider ice dynamics, and as such, the long-term stability of the glacier. Coupled ice-flow - ocean models (frameworks of which have been recently developed, e.g., Zhou et al. (2024)) or offline coupling of an ice-flow model with an ocean model which comprehensively resolves basal melt rates in response to ocean dynamics is indeed a subject of future work in this field. We have revised our discussion to elucidate the context and limitations of our calculations and measures that need to be taken to overcome these limitations in the future. Please see lines 366 – 380. I like this paper, although I do also think that it could benefit from some revisions to improve readability. I understand that Nature Communications wants/requires the methods to come at the end of the paper, but it felt to me that the methods was initially written as section 2 (before the results) and only later moved to the end of the paper, with the consequence that the results were difficult to interpret. For example, it took me awhile to understand what was meant by the various model scenarios. I think the introduction would benefit from a few sentences, or perhaps another paragraph, that would more fully describe the model set-up and simulations in layman’s terms. Thank you for bringing this to our notice. We have addressed it by amending the Introduction chapter to include more information regarding both the model and experiment setup. Please see lines 59 – 100. Finally, I found the use of acronyms and variables within text a bit cumbersome. I encourage the authors to eliminate them to the extent possible in order to make the paper more accessible. “basal melting”, “subglacial discharge”, “friction velocity”, and “grounding line” are much easier to read than m_b, Q_sg, u*, and GL. I don’t see why the glacier and fjord names need to be abbreviated, and some acronyms are not used very often (such as GrIS and RG). There are many places where the text can be tightened up if space is an issue. Thank you for the suggestion. To our knowledge, we are within the word limit guidelines provided by the journal. Therefore, we have removed the abbreviations for “basal melting”, “subglacial discharge”, “thermal driving”, “friction velocity”, “grounding line”, “Greenland Ice Sheet”, and “Ryder Glacier”, as well as all other glacier and fjord names from the text for clarity. Specific comments: Title: I feel that the title is a little misleading, although I agree with the sentiment, as this study is not actually modeling glacier flow or fracture propagation. This is true, we do not model the glacier flow or fracture propagation. Hence, we preferred the use of “may destabilize”, as we posit based on our ice volume change calculations that the high basal rates modelled in the channels (ca. 200 m/yr) in the RCP 8.5 scenario could be unsustainable for the ice shelf to withstand for longer periods of time. L20: How do “channels carve incisions”? Aren’t the channels themselves products of melting? Thank you for identifying this mistake. We have fixed it (see below) together with the comment on L21. L21: Are the thickness values stated for the areas immediately above the channels, or is this a statement about the average thickness of the ice shelf? Also, I assume that at least some of the thickness gradient is due to longitudinal extension, but this reads as though it is entirely due to submarine melting of channels. The statement was regarding the (average) thickness of the ice shelf – it ranges from 400 – 600 m deep near the GL to 50-100 m near the calving front. You are right, by putting it together with the channel description, it appears to come across as being entirely driven by submarine melting. We have fixed this in the revised manuscript. Please see lines 22 – 28. L54: Please provide a definition of friction velocity here. We have added a definition of friction velocity. Please see lines 95 – 96. L113: What do you mean by “peaks”? Fig. 3 shows that the temperature continues to increase as you increase the subglacial discharge, just at a lower rate. You are correct, thank you for identifying this mistake. “Peaks” is not the right word. We have changed this in the revised text. Please see lines 177 – 179. L117-L118: Similar to previous comment. Contribution to what? And do you just mean that the thermal forcing increases at a lower rate? Indeed, this is a more appropriate formulation. We have adapted it in the revised text. Please see lines 181 – 183. Fig. 3: What does “27.9” refer to? Thank you for noticing this. We forgot to describe it in the caption. It has now been added. L213-216: I’m not sure that I follow this. Are you essentially saying that the addition of cold freshwater will limit the thermal pump? No, that is not what we intend to say, we apologize for the lack of clarity. We wanted to conclude the paragraph by presenting a broad discussion regarding the relative contributions of rising ocean temperatures (Atlantic Water reaching Petermann by the end of this century) vs increased subglacial discharge in driving increased melting in a warmer climate. We have made revisions to clarify our argument. Please see lines 310 – 327. L370-373: How sensitive is the model to variability in subglacial discharge? Your simulated discharge follows a boxcar function, while we know that subglacial discharge follows more of a gaussian-like pattern during summer and is highly variable. I found myself wondering how important the variability is. If you have a few days or weeks of exceptionally high or low subglacial discharge, can that push the model into a different regime? Thank you for raising this point, as we have given it due consideration as well. We are currently in the process of setting up the Glacier Drainage System Model (GlaDS) for the Petermann catchment to determine both the magnitude and location(s) of emergence of subglacial discharge. The gaussian-like profiles of subglacial discharge retrieved from GlaDS reveal considerable seasonal and inter-annual variability. However, this work is still in its nascent stages, and we are yet to analyze the corresponding ocean model simulations. It is difficult to quantify the sensitivity without running appropriate experiments using such a setup, however, we have now included a discussion around it, and caveated its speculative nature, suggesting the need for further studies. Please see lines 267 – 275. L419-430 (and Fig. 6): I found this description confusing. It’s not clear to me what exactly you are calculating. Are you calculating the volume of ice above a channel? Maybe I am just hung up on the term “channel ice volume”, when to me the word channel implies a place where there is no ice. This is correct, we are calculating the volume of ice above a channel. Thank you for bringing this to our attention. We have replaced all occurrences of the phrase “channel ice volume” with “ice volume above a channel”. The section is renamed to “Changes in ice volume above basal channels over time”. Please see, for e.g., lines 611 – 624. Table 1: Is Q_sg-median just the average the present and RCP 8.5 scenarios? Or is that just coincidence? It is just the average of the two scenarios. Subglacial discharge from Slater et al. (2020) was only available for present and the RCP 8.5 scenario. We generated Qsg-median manually by taking the median of the present and RCP 8.5 values. We have made adjustments to the text to reflect this better. Please see lines 580 – 585. Thank you for taking the time to read our manuscript. We are pleased to hear that you found it interesting and well-written. We agree with all of the suggestions that you have provided and have revised our manuscript accordingly. Please note that line numbers provided below refer to the track changes document (appended to this report for convenience), where new text appears in blue and removed text is marked with red strikethrough. Firstly, we would like to take this opportunity to inform you that during the revision process, we noticed a (minor) mistake in our calculation of the drivers of melt, namely, thermal driving and friction velocity. While the updated numbers are slightly different, we note that this does not impact any of the findings of our study, and the revised figures remain qualitatively identical. We humbly request to proceed with these corrected values, and sincerely apologize for the oversight and for any inconvenience this may cause. Please see the revised percentages in Lines 133 – 134, 167 – 170, 184, 195, and 227, and Figures 1(d), 2(e – l), and 5(c,d). Please find our detailed (point-by-point) responses below. Minor comments Lines 3-4: “contemporary mass loss rates are underestimated”. Underestimated by whom? Observational studies, or ice sheet models? Thank you for pointing this out. Greene et al., 2024 suggest that observational estimates of GrIS mass balance have underestimated recent (1992 – 2020) mass loss by ca. 20%. We have fixed this in the revised manuscript. Please see lines 3 – 4. Line 14: “ice volume of ca. 4m global mean sea level equivalent” I can’t understand from this sentence which sectors of the ice sheet you are referring to. Please check the reference too – is Morlighem et al. (2017) not more appropriate here (BedMachine, see reference list below)? We agree with the concerns raised above. In Line 5, we begin the thread about the northern sector of the ice sheet, however, the passage is quite long. Furthermore, the use of “adjacent” is confusing, as it implies the exclusion of the sector Petermann corresponds to. Upon reviewing, we do agree that Morlighem et al. (2017) is appropriate. However, we maintain that Mouginot et al. (2019) is also relevant, as they quantified the sectorial mass balance of the GrIS in unprecedented detail. Herein, the integrated contribution from the northern sector amounts to ca. 4 m. We have fixed both the phrasing and citation in the revised manuscript. Please see lines 14 – 17. Line 16: “Erstwhile” is not wrong here, but it’s a rarely used word. You could consider replacing with “Formerly” or “Previously” to aid non-native speakers. We agree, and have replaced it with Formerly. Please see line 20. Line 21: “thinning it from ca. 600 m near the GL to ca. 100 m near the calving front.” Phrased this way, the sentence can be misinterpreted… do you mean “with ice thickness over the channels of 600 m near the GL to ca. 100 m near the calving front”? Reviewer 1 has also raised concerns regarding the phrasing of this sentence, and we apologize for the confusion it has caused. In Lines 18-21, we aimed at detailing the average thickness of the ice shelf, extending from the grounding line (400 – 600 m) to the calving front (50 – 100 m). However, we also do not want to attribute the thinning completely to basal melting (Reviewer 1). We have rephrased this sentence in the revised manuscript, bearing in mind both the comments. We first document the average thickness gradient from the grounding line to the calving front (driven by all factors). Thereafter, we describe the basal melt channel geometries. Please see lines 22 – 28. Figure 1b: I suggest labelling the calving front with a date. Thank you for the suggestion. We have now included the date (2009) for improved clarity. Please see Figure 1(b). Line 34-36: “estimated area-averaged melt rate maximum of up to 170 m/yr” the phrase “area-average” initially implies that this was estimated over a large sector of the shelf. I suggest rephrasing slightly to make clear that this estimate was calculated only within the basal channel. This is true, thank you for noticing this. We have rephrased it in the revised text. Please see lines 44 – 51. Line 38: “with warming expected to amplify” I suggest clarifying whether you mean atmospheric warming, ocean warming, or both. Thank you for noticing this. We meant both, and have clarified this in this revised version. Please see line 54. Final paragraph of introduction (lines 43-59): It’s unusual to see the results summarised so thoroughly at the end of the introduction. It might be more helpful to give a very brief summary of your experiment design instead e.g. from line 49 onwards something along the lines of “We performed a control run with zero Qsg, and three test simulations designed to simulate present day subglacial conditions, increased Qsg under high-end scenario climate warming, and an intermediate Qsg scenario between the two.” Thank you, the results are indeed described in detail here, and without proper context (i.e., experiment setup), its purpose is not really served. We have introduced a summary of our experiment design, and have condensed the description of the results in the revised manuscript. Please see lines 59 – 100, which also combines suggestions from both Reviewers 1 and 3. Paragraph beginning line 62: This paragraph jumps into discussion of modelled melt rates without describing the results clearly first, which means the reader has to work hard to catch up. I suggest an opening sentence or two describing the melt pattern in the control and preset day runs, before going on to compare them to observations. Thank you for providing this perspective. We agree with providing more context first, and to also provide a more detailed comparison between data and our modelled results (see Reviewer 3’s comment). We have fixed this paragraph in the revised manuscript. Please see lines 103 – 127. Figure 2: You could make the model results larger and easier to see by removing unnecessarily duplicated color scales and axis labels (see also general comments on figures below). Thank you for the recommendations. We have removed the duplicated axis labels in favour of common X (Longitude) and Y (Latitude) axis labels. The duplicated color scales have also been replaced with one common color scale for each set of panels representing changes in melt rates, friction velocity, and thermal driving. A delta symbol has also been added to the c-axis labels to denote changes. Please see Figure 2. Figure 3: The isopycnal label (27.9) is large and distracting (see also general comments on figures below). Thank you, we have replaced it with a smaller font size that is consistent with the rest of the labels. Please see Figure 3. Figure 5: The lines here can be a little hard to see – consider changing the color scheme to improve readability and increasing line thickness (see also general comments on figures below). Thank you for the feedback. We have changed the color scheme and increased the line thickness. Please see Figure 5. Line 151 & Figure 6: The two sets of location names are a bit confusing here (Sw, Sc, and Se vs. Lw, Lc, and Le). Are they both needed? Lw, Lc, and Le are barely mentioned in the text, and could perhaps be removed? We agree, having two sets of location names is confusing and not required. We have removed all occurrences of Lw, Lc, and Le. Labels in panels a – d, Figure 5 have also been removed. We have kept the green stripes in Figure 5(a) to highlight the locations of the prominent basal channels. Figure 6: The pink dots in panel 6a are not described in the caption, and the symbols marking different lines in panel 6b are far too small to be readable (see also general comments on figures below) Thank you for noticing this, we have now added the description to the caption. We noticed that the cumbersome legend of Figure 6(b) can be conveniently added to the figure description, and as such has been removed from the panel. Line 222: Should read “buoyant meltwater rises” Thank you for noticing this. Following comments from Reviewer 3, we have rephrased this sentence. As such, this fix is no longer applicable. Discussion line 219-240 & Figure 6b: Given the very straight erosion trajectory in figure 6b, I initially thought you must have applied the summer value for Qsg continuously over the whole year. Looking at the methods it seems that Qsg does go to zero outside of the summer, so why don’t we see a seasonal cycle in the channel volume? Presumably over winter the melt substantially decreases but the ice advection should change much less? I think the discussion section could do with a bit more context on seasonality, including how seasonally varying ocean temperatures affect erosion rates. Thank you for raising this concern. We do not see a seasonal cycle since we use mean annual basal melt rates and advection rates, and calculate the yearly change in ice volume above the channels until it reaches 0. The annual mean basal melt rates in these channels are still considerably high because of the discharge-driven strong summer signal. While it does not explicitly concern the issue raised in this comment, we note that we are also limited by the lack of systematic high (temporal) resolution data for present and future warming scenarios that model the seasonal variability of subglacial discharge (gaussian-like profiles) and ice advection rates. Indeed, during winter, basal melt rates are lower, and ice advection rate does not vary significantly. Muenchow et al. (2012) have reported mean advection rate of ca. 1 km/yr with a modest seasonal modulation of ca. 0.1 km/yr. Our model is consistent with data from the fjord which suggests that seasonal ocean temperatures in Petermann Fjord vary by ca. 0.1 ºC, which does not have as significant an impact on basal melt as subglacial discharge. We have now expanded our discussion to quantify the impact of seasonal changes in the fjord (ocean and subglacial discharge) on erosion rates in the revised text. Please see lines 328 – 339. Line 266-267: There has been further recent work on this subject in Antarctica beyond your references, so I think the statement “insights into this process are limited to tacit reasoning” is a little unfair. I suggest reading and citing: Goldberg et al. 2023, Gwyther et al. 2023, and Pelle et al. 2024. Thank you for providing us with this information. Upon reviewing the recent literature, we have revised the text and added these as references. Please see lines 405 – 415. All figures: Please take a look at the text within your figures and make sure it meets minimum font size guidelines, but is not excessively large and distracting (as with some subpanel labels). Some axis labels and legends are very small indeed. In general, using consistent font across all figure will also improve perception of your figure quality. You could also improve clarity of figures 2, 3 & 4 by indicating on colour scales where a variable is absolute or relative (e.g. “Salinity” in panel 4a and “Change in salinity” in panel 4b, use a delta symbol if space is tight). Thank you for providing detailed suggestions. We have now implemented a consistent font size across all our figures which meets the font size guidelines. Absolute and relative values are now clearly differentiated by the use of a delta symbol. Reviewer #2 (Remarks on code availability): I have checked that the code is accessible, but having never used the FVCOM model I can't say whether it is sufficient to reproduce the results Thank you for raising this concern. We can assure that the model code and input files, together with the detailed instructions that accompany them, are sufficient to reproduce the results."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have sufficiently addressed all of my prior concerns. I find this to be an interesting paper that advances our understanding of glacier-ocean interactions. However, I have three related editorial comments regarding the paper that are in the spirit of trying to increase its impact. 1. To me, the real impact of this paper is showing how subglacial discharge (and changes in subglacial discharge) affects submarine melting of ice shelves and not simply that increases in subglacial discharge will affect the stability of Petermann Ice Shelf. Previous studies, such as the work by Slater and Straneo, as well as that of Motyka et al. and probably many others, have shown that changes in subglacial discharge can lead to large changes in melt rates. What I take away from the paper is that the effect of subglacial discharge (i) at low fluxes is to draw in warm water and breakdown the boundary layer, whereas (ii) at high fluxes the temperature effect becomes saturated and further increases in discharge essentially affect the boundary layer. Much of the paper focuses on this. The uncertainty in the model topography, location and number of discharge outlets, and lack of ice dynamics model causes me to question the assertion that future runoff may destabilize Petermann Ice Shelf. (Not that I disagree, I just don’t think that’s really what this paper shows.) Personally, I think the deeper understanding of how runoff affects submarine melting to be more interesting and robust than whether or not one particular glacier will destabilize, but maybe that is just my personal bias coming through. 2. The number of panels in figures 2 and 3, and the long figure captions, make them somewhat difficult to digest. I don’t think the figures need to be changed all that much, but I do think it would be helpful if they were labeled somehow so as to make it clear which panels correspond to the control run and each of the runoff scenarios (as was done for figure 4) 3. I think the writing can be tightened up and made more succinct in many places. As just an example, the paragraphs in the discussion section are quite long (almost a page long in draft mode). This can make it difficult to determine (or remember!) the main point of each paragraph. There are quite a few clauses and phrases throughout the paper that are either unnecessary or that could be shortened to improve readability.", "reviewer_2": "The authors have been very thorough in their response to reviewers comments, and have made substantial changes to the text and figures. I feel that all of my comments have been very well addressed.", "reviewer_3": "The authors have sufficiently met my concerns in their revised manuscript and I have no further substantive suggestions. I would commend the authors on their complete and thorough addressing of all (not just my) reviewer suggestion. The updated manuscript has improved readability and the key points are more clearly made. The dropping of terms in favor of words (e.g. Qsg —> subglacial discharge) improves readability. The added ‘in sum’ at the end of the results section is helpful. I appreciate the new concluding sentence of the manuscript and agree — the interplay between surface melt and seawater intrusions at the grounding zone are a real research frontier right now. Minor notes for the authors: 2 First sentence scans somewhat awkwardly, suggest rephrasing to clarify 37 ‘studies’ 50 sea-ice ice-shelf (I think?). Is there a coupled ice flow model here? In looking to your methods (438) it looks like there is a coupled ice shelf model, but I would love one more sentence here describing it. For those not familiar with the Zhao ice shelf model, what is it? What ice rheology is used? Is it a higher-order model? 88 grounding line are 135-148 the continued reliance on Qsg, mbpgis, U*300m and other variable-based terminology really does make reading this difficult for someone unfamiliar with the work. I know one round of this has already been done, but I would suggest one more round of ‘decreasing variable terminology’ to improve readability. 228 estimates (somewhat awkward flow here) 234 are currently being experienced 284-89 maybe I just need an afternoon coffee, but this sentence is far too long to be digestible without multiple readings. 305 Idle curiosity: does your model allow for migration of channels, as has been observed on Thwaites and elsewhere? Do you observe migration of channel location at all? E.g. Chartrand et al., 2020. I wonder if this could delay Petermann collapse."}, "author_response": ""}]

Decision letter