Subglacial Freshwater Drainage Increases Simulated Basal Melt of the Totten Ice Shelf

Subglacial freshwater discharge from beneath Antarctic glaciers likely has a strong impact on ice shelf basal melting. However, the difficulty in directly observing subglacial flow highlights the importance of modeling these processes. We use an ocean model of the Totten Ice Shelf cavity into which we inject subglacial discharge derived from a hydrology model applied to Aurora Subglacial Basin. Our results show (a) discharge increases melting in the vicinity of the outflow region, which correlates with features observed in surface elevation maps and satellite‐derived melt maps, with implications for ice shelf stability; (b) the change in melting is driven by the formation of a buoyant plume rather than the addition of heat; and (c) the buoyant plume originating from subglacial discharge‐driven melting is far‐reaching. Basal melting induced by subglacial hydrology is thus important for ice shelf stability, but is absent from almost all ice‐ocean models.


10.1029/2023GL103765
2 of 11 alterations in larger-scale fjord circulation (Slater et al., 2018) demonstrate that the impact of subglacial hydrology can be a significant component of ice-ocean interactions. Motyka et al. (2003) proposed that subglacial discharge drives strong buoyant convection, which forms an ascending plume, entrains warmer waters and drives ice front melting. Jenkins (2011) extended this concept for Antarctica, by unifying the concept of convection-driven melting (where buoyancy generates a plume that subsequently drives melting) with the commonly understood melt driven-convection paradigm (where melting releases a buoyant plume; e.g., the ice pump mechanism, Lewis & Perkin, 1986).
To date, the only model to simulate the input of Antarctic subglacial discharge into a realistic ocean model is Nakayama et al. (2021). SGFW discharge was directed into the ice shelf cavity from topographical routing estimates (as opposed to direct hydrological modeling) and included at two main output locations below Pine Island Ice Shelf in a high resolution ice shelf-ocean model. They showed substantial melt increases in the vicinity of drainage output and suggested that subglacial drainage is necessary to match satellite melt estimates (Nakayama et al., 2021).
In our study, we combine outputs from a subglacial hydrology model into an ice shelf-ocean model to assess not only how the steady-state freshwater flux impacts ice shelf melting, but also why melting is impacted.

Study Area
The Totten Glacier and adjoining Totten Ice Shelf (TIS) are the most rapidly changing region of the East Antarctic (Smith et al., 2020) and currently hold back ∼3.5 m of sea level rise potential from the Aurora Subglacial Basin (Greenbaum et al., 2015). Observations of this region suggest that it is undergoing ice surface lowering, retreat and mass loss (e.g., Flament & Rémy, 2012;Li et al., 2015;Pritchard et al., 2012;Roberts et al., 2018;Velicogna et al., 2020), although complex variability (Adusumilli et al., 2020;Gwyther et al., 2018) limits the ability to detect and attribute recent change. The observed oceanographic regime Rintoul et al., 2016;Silvano et al., 2017Silvano et al., , 2019T. P. Tamura et al., 2022;Williams et al., 2011) and bathymetric connections (e.g., Greenbaum et al., 2015) beneath and adjacent to the TIS indicate the potential for warm, off-shelf waters to access the cavity and drive strong basal melting.
Previous modeling (e.g., Gwyther et al., 2014;Van Achter et al., 2022;Xia et al., 2023) indicated that inflow into the cavity (Figures 1a and 1b) is on the eastern side and can carry relatively warm water (see Figure S1 in Supporting Information S1). Ice shelf meltwater flows out of the cavity adjacent to the western grounding line as a topographically constrained boundary current, before exiting the cavity.
It is unknown how subglacial discharge, which drains into the cavity in relatively high volumes , impacts ice shelf melting and hence stability in this sensitive region.

GlaDS
The Glacier Drainage System (GlaDS) model is a 2D finite element model of subglacial hydrology (Werder et al., 2013). It includes a distributed system that develops on the elements where water depth, pressure and flux are calculated. Channels develop along element edges and form naturally as result of the pressure and flux conditions in the adjacent distributed elements. The application of GlaDS to the Aurora Subglacial Basin feeding into the TIS is described in Dow et al. (2020). Inputs include ice surface and bed topography interpolated from airborne radar lines, along with ice basal velocity and basal melt rates calculated from the Ice Sheet System Model. Flow routing in both the distributed and channelized system, including the location of outflow over the grounding line, is driven by a combination of the topography and the variability of water pressure in the system. The model outputs are from runs that have reached steady state.
Three sensitivity experiments were conducted with a range of flux conditions into the TIS cavity. As the total discharge volume did not vary greatly between experiments, the resultant impact on melt rates is low and we present only the medium conductivity results here (see Figure S2 in Supporting Information S1). The maximum inflow is from a channel at the head of the ice shelf (Figure 1c), with a peak output of 25 m 3 s −1 , contributing ∼78% of the total SGFW inflow.

ROMS
We use the Regional Ocean Modeling System (ROMS; Shchepetkin & McWilliams, 2005) with ice-ocean thermodynamics and mechanical pressure (Dinniman et al., 2007;Galton-Fenzi et al., 2012). Heat and salt exchange across the ice-ocean boundary layer is captured through the three-equation parameterization (Hellmer & Olbers, 1989;Holland & Jenkins, 1999) with velocity-dependent exchange coefficients (McPhee et al., 1987) and a quadratic drag formulation (C D = 0.003). Vertical mixing is parameterized with the K-profile parameterization (Large et al., 1994). The model has 31 sigma-coordinate layers, with approximately 2.5 km horizontal resolution at the ice shelf. The model domain is the same as in Gwyther et al. (2018), and includes the bathymetry estimated from gravity inversion (Greenbaum et al., 2015) and ice draft from airborne radar soundings (Figures 1a and 1b). The experimental conditions include tidal forcing (following Galton-Fenzi et al., 2012), repeat 1992-2017 climatological forcing of the surface and lateral boundaries, with a spin-up period of 20 years. Lateral boundary conditions (velocity, temperature and salinity) are from ECCO2 (Menemenlis et al., 2008); surface boundary conditions include wind forcing (ERA-interim; Dee et al., 2011) and buoyancy fluxes (T. Tamura et al., 2016) to parameterize sea ice formation.
The existing river source module within ROMS (e.g., Zhang et al., 2009) was employed to simulate subglacial discharge. There is an addition to the horizontal momentum term (at the rate described by the river source flux), with associated tracer fluxes diverged by the advection operator, at each cell face identified as a channel outlet. As there is no extra water volume input, it should be thought of as a modification to the existing conditions in the cell. Input was vertically distributed in a linearly increasing manner toward the top of the column. This is to ensure that SGFW entering the model is not buoyantly unstable, which could lead to numerical instabilities. This represents a balance between the realistic entry location (likely at the narrow grounding zone located just above the seabed) and model requirements (where the generation of unstable stratification can lead to numerical instability). Discharge is set to enter at the freshwater pressure melting point with zero salinity.

Satellite Melt Rates
Melt rates are inferred with the mass conservation approach constrained using satellite data following Gourmelen et al. (2017). This method combines satellite observations and climate models balancing mass loss from basal melting, surface lowering through volume divergence, with mass gain from surface processes. Lagrangian surface elevation change is generated from CryoSat-2 swath radar altimetry from 2010 to 2020 (Gourmelen et al., 2018)  (Howat et al., 2019). The background image is the MEaSUREs MODIS Antarctic mosaic . and surface ice velocity is from the ITS_LIVE data set (Gardner et al., 2018), with a full description given in Gourmelen et al. (2017). Surface mass balance is obtained from the Regional Atmospheric Climate Model version 2 (RACMO2; Wessem et al., 2018).

Experiments
To assess the impact of SGFW on melt rates and distribution, a series of experiments are performed with ROMS. All experiments are conducted as sensitivity studies from the end of the spin-up phase, and consist of a single, perturbed year as ocean state evolution beyond this period limits our ability to perform comparative analyses. This study consists of two sets of experiments.
The first set compares a "no-flow" control scenario (addition of no SGFW) to a sensitivity run with the addition of SGFW ("with SGFW"; see GlaDS description and Supporting Information S1).
The second set of experiments compares two further model runs with subglacial discharge to the control experiment (no-flow), but the input SGFW flux is switched between having only the temperature tracer modified ("temperature-only" experiment) or only the salt tracer modified ("salt-only experiment"). The other tracer is unmodified.

Mean Melt Rates
Melt rates estimated with the satellite-inferred method (Figure 2a) show highest melt through the western ice shelf and grounding zone and weaker to moderate melt over much of the rest of the shelf. This leads to a TIS-wide area-averaged melt of 12 ± 1.4 m yr −1 or a total melt mass flux of 64 ± 7.4 Gt yr −1 . Autonomous Phase-sensitive Radar Echo Sounding (ApRES) observations taken on the southern grounding line region of the TIS in 2017/2018 (Figures 2a and 2b), showed melt rates of 22 ± 2 m yr −1 (Vaňková et al., 2021). In (c), the percentage difference in time-mean melt rate between the experiments with SGFW and without SGFW is shown. In (a) melt rates have been removed in regions where heavy crevassing impedes the melt rate algorithm (e.g., at the calving front). In (b and c), the model calving front and grounding line are delineated with orange lines. Background image is from the MODIS Mosaic Of Antarctica . Locations mentioned in text are marked in panel (a), as well as the north direction. The TI05 Autonomous Phase-sensitive Radar Echo Sounding site (Vaňková et al., 2021) is marked with a circle on the same color scale in (a and b). X and Y coordinates are on a polar stereographic projection (standard latitude of 71°S).
Melt rates for ROMS with the SGFW flux ( Figure 2b) are high (maximum of 50 m yr −1 ) in the deeper grounding line regions and through the eastern channel region, while melting across the northern, shallower portions of the ice shelf is generally less than ∼5 m yr −1 . The TIS-wide area average melt rate is 8.3 m yr −1 (a total melt mass flux of 54 Gt yr −1 ).
The glaciological/satellite method shows higher melting than the ocean model in the deep grounding line region and on the western ice shelf flank; outside of these regions, both products compare relatively well. The simulated melt rate compares well with the observation of melt taken by radio echo sounding at that location (compare 21 m yr −1 in the simulation and the observation of 22 ± 2 m yr −1 ).
Multiple reasons exist for differences between modeled and observed melt rates. First, the model forcing epoch differs from the satellite-product time period. Ocean models are also reliant on accurate bathymetry (e.g., Goldberg et al., 2019), which is still unmeasured beneath most of the TIS. Furthermore, horizontal resolution, while relatively high in this simulation, may not capture small-scale flow features (Dinniman et al., 2016). Like other ice shelf-ocean models, sub-grid-scale parameterisations are used for capturing many unresolved processes-however, these are likely a source of error. Finally, while we have included SGFW discharge, errors in geothermal heating may lead to an underestimation of subglacial fluxes. Sources of error in the satellite-derived product include from the surface accumulation model, firn compaction, and vertical strain profile (e.g., S. Cook et al., 2023;Vaňková & Nicholls, 2022). Nevertheless, the modeled melt rates display a similar pattern to satellite-inferred melt rates, and agree with the Autonomous Phase-sensitive Radar Echo Sounding observations, making them suitable for exploring the impact of subglacial drainage on melt rate.

Exploring the Addition of SGFW
The comparative impact of subglacial drainage on ice shelf basal melting is small overall but strong in localized regions (Figure 2c). The percentage difference in melt rates (where a positive value indicates higher melting with SGFW) shows that while the area-average increase in melting is ∼3.3% (an increase in mass flux of 1.38 Gt yr −1 ), there are strong, local variations due to the input of subglacial water. Change in melting is generally strongest beneath the deepest, thickest sections of the ice shelf (adjacent to the grounding line; Figure 2c) where discharge increases melting by 20%-30% (with peak increases of ∼50%). For example, for the region of the ice shelf where the UTM northings coordinate is greater than −1,030 km, melting increases on average by 29%. This results in a total increase in melt mass flux of 1.0 Gt yr −1 over a relatively small 650 km 2 . Enhanced melting also occurs along the northwestern side of the ice shelf (in the location of the outflowing topographically-constrained boundary current) and in the eastern channel ( Figure 2c). This can be compared with a broader region along the southeastern (inflow) edge of the ice shelf where melt rates decrease by 2%-3% compared to the no-flow run ( Figure 2c). The mean change in melt for all cells exhibiting an increase is 4.9%, while the mean percentage decrease in melt for all cells with decreased melt is 1.6%.
The reduced melting on the eastern side of the ice shelf results from cooler waters at the ice/ocean boundary (Figures S2c and S2d in Supporting Information S1). This occurs due to the increased presence of cold meltwater (Figures S3a and S3b in Supporting Information S1) and increased inflow of colder surface waters from outside the cavity. Increased inflow of warmer water occurs at depth ( Figure S4 in Supporting Information S1), but this flows to the grounding line with minimal contact with the intermediate waters above (Figure S3c in Supporting Information S1).
The spatial distribution of subglacial discharge is shown with a passive tracer dye in Figure 3 (also see Movie S1). The SGFW dye enters at the grounding line ( Figure 3a) and by 30 days has reached the calving front (Figure 3d). On leaving the cavity, the SGFW dye enters the Antarctic Coastal Current, where is advected westwards around Law Dome (Figures 3e-3g). Within 5 months (Figure 3h), the SGFW dye has traveled over 500 km.

Single Tracer Experiments
In Figure 4, the change in melting compared to the "no-flow" run is shown for the inflow of SGFW with only a modified temperature tracer or with only a modified salinity tracer. Enhanced melting is shown to be due to freshening from the SGFW, and with little to no impact from the addition of the heat in the SGFW (Figure 4a).
To better understand the dominance of freshening in the enhanced melting signature, we explore the factors driving the melt rate (Figure 4b). The two parameters that control basal melting are the thermal driving and the 6 of 11 friction velocity (Figures 4c and 4d). The thermal driving is defined here as  (Figure 4c), as shown by the correlation between melting and T * . The higher T * and melting toward the grounding line results from the depression in the freezing point under thicker ice. The influence of ocean currents on melting is less obvious, although regions of increased u * often coincide with increased melting, noticeably, for example, near constricted geometry or pinning points (Figure 4d).
The temperature-only scenario has minimal impact on melting (Figure 4e) compared to the control scenario (mean of −0.2%), with only some minor differences near to the ice front emerging in the thermal driving and velocities (Figures 4f and 4g). Given that this region is most proximate to the dynamic open ocean, it is not surprising that minor differences could evolve here between the control and sensitivity scenarios.
The salinity-only scenario, however, has a very strong impact on melting (mean of 3%), essentially entirely explaining the impact of SGFW on melt rates; for example, compare Figures 2c and 4h. Thermal driving is strongly modified by subglacial discharge (Figure 4i), even though no extra heat is being added. Therefore, this increased T * must result from the increased buoyant circulation (see increased u * in Figure 4j) entraining and supplying more ocean heat from deeper ocean waters to the base of the ice shelf. This demonstrates that the addition of fresh and thus buoyant water into the cavity enhances flow past the ice shelf and is the principal cause of the stronger melting.

Subglacial Discharge Has a Strong Local Impact on Melting
The area-averaged increase in melting due to subglacial flux is relatively weak when examined over the whole TIS (around 3.3% or 1.4 Gt yr −1 of extra ice lost). However, some of the increased mass flux is balanced by a reduction in melting in the "oceanic inflow" region of the ice shelf (Figure 2c), which results from altered cavity circulation and inflow of cooler surface waters.
For areas of increased melt there is a mean increase of 5%, with strong localized impacts of 25%-30% in the deep grounding line and eastern channel regions, where the discharge rate from subglacial channels are ∼25 and  1.4 m 3 s −1 , respectively. This result agrees with Nakayama et al. (2021), who demonstrate very strong localized increases in melting near to subglacial outflow locations. The impact of these concentrated discharge plumes and corresponding high shelf melt rates are demonstrated by the presence of ice shelf channels. The shelf channels can be seen in satellite elevation products (e.g., the Reference Elevation Model of Antarctica Howat et al., 2019) and advect from the grounding line at the locations where subglacial discharge is greater than 1 m 3 s −1 toward the ice shelf terminus (Figure 1c). In addition, satellite-derived elevated ice shelf melt rates and corresponding grounding line retreat are also focussed in these regions of high subglacial discharge and high ROMS melt (Adusumilli et al., 2020;Flament & Rémy, 2012). Here our data indicate that the role of subglacial discharge is key for high melt in these sensitive areas. Ice loss at the grounding line (as opposed to in the middle of the ice shelf) will lead to more rapid grounding line retreat and, depending on the bedrock geometry, the potential initiation of an unstable retreat scenario (Marine Ice Sheet Instability; Tsai et al., 2015;Weertman, 1974). Nakayama et al. (2021) show higher melt rates are primarily explained by increased ocean currents, but suggest increased thermal driving could be due to the injection of subglacial discharge, which is warm relative to the seawater freezing point. Our experiments suggest that the thermal component of subglacial water has very little impact on ice shelf melt. Indeed, if subglacial discharge is at the pressure freezing point, then any mixture between it and ambient waters will likely result in reduced thermal driving, except potentially in the case of supercooled ambient waters. This is due to the resultant SGFW-ambient water mixture being closer to the freezing point line (as illustrated with warm and cold cavity conditions on a θ-S diagram with Gade lines (Wåhlin et al., 2010) in Figure S5 in Supporting Information S1).

Increased Melting Results From Freshening and Addition of Buoyancy
Instead, regions of increased melting result from the increased buoyancy-driven plume (Motyka et al., 2003) and are independent of temperature. The role of subglacial water, at least in the vicinity of subglacial channel outlets, is therefore not to melt directly, but to rapidly initiate strong buoyant plumes, entrain warm ocean waters to the ice shelf boundary layer, and increase melting, as simulated by Jenkins (2011). Beyond the convection-driven melt zone, release of ice shelf meltwater (rather than subglacial discharge) is the primary contributor to buoyancy-driven melting.
The version of ROMS applied here makes the hydrostatic assumption and so cannot resolve non-hydrostatic dynamics. Nevertheless, if there is strong stratification and weaker flow speeds, then by the criterion of validity for the hydrostatic approximation (Marshall et al., 1997), the approximation we make in the vicinity of the channel outlet likely still holds. Furthermore, the dynamics of the buoyant circulation initiated by SGFW discharge, beyond the grounding line and at-or-greater than the grid resolution, are largely identical to other sub-ice shelf plumes, which have been simulated by many similar hydrostatic models (e.g., Nakayama et al., 2021).

Small Volumes With Broad Impacts
The volumes of subglacial water necessary to increase entrainment of ocean waters and hence basal melting, are small. Although the greatest increases in melting occur where the greatest freshwater discharge enters the cavity, there is also enhanced melt in regions where subglacial discharges are as low as 3 cm 3 s −1 . Furthermore, the subsequent outflow plume originating from subglacial-driven melting, identified by the region of increased melt on the western flank of the ice shelf, is ∼5 km wide and extends ∼100 km in length, from the grounding line to the ice shelf terminus. This suggests that, given the widescale subglacial drainage around Antarctica (Le Brocq et al., 2013), enhanced melting from this process is extensive, with the strongest impact on melting occurring at the grounding line in areas adjacent to large subglacial outflow channels. As such, it is important that the impacts of subglacial drainage are included in ocean-ice shelf models to project ice shelf evolution in the Antarctic.
It is unclear at what stage future atmospheric conditions and drainage due to surface melt will begin to significantly influence the volumes of subglacial outflow in Antarctica. Furthermore, it is unclear what the impact of rapid and strong pulses of discharge would be, for example, from lake drainage (e.g., Fricker et al., 2007;Warner et al., 2021).

Conclusion
We find increases in TIS melting resulting from SGFW outflow. Increases are locally strong in regions near the grounding line, likely important for maintaining ice shelf stability. The formation of a strong buoyant plume, initiated by SGFW, ascends the sloping ice shelf base and subsequently increases melting. The spatial impact of discharge and increased melting is shown to reach far beyond the immediate grounding line zone.
Changes in ocean heat supply are likely the most important factor for ice shelf stability. However, to understand current and future melt rates, we need to account for the additional impact of convection-driven melting resulting from subglacial discharge. Should surface melting increase above present-day levels and this additional meltwater reaches the bedrock, the volume of subglacial discharge would likely be significantly more -and the impact on basal melting profoundly different.