Suppressed basal melting in the eastern Thwaites Glacier grounding zone – Nature

Borehole observations

A Sea-Bird Scientific SBE 49 FastCAT CTD profiler was used to observe the water-column structure between 9 January and 12 January 2020. A total of 15 individual CTD casts were completed, sampling at a rate of 16 Hz. Before each deployment, the CTD was stored in a warm bath (approximately 5 °C) to minimize icing in the conductivity cell during the profiler’s traverse of the approximately 90-m air-filled portion of the borehole. The CTD data were processed using standard routines in the Sea-Bird data-processing software version 7.26.7.129 and each profile was averaged into 0.1-m vertical bins. Absolute salinity and conservative temperature were computed using the Gibbs SeaWater (GSW) Oceanographic Toolbox for TEOS-10 (refs. ^53,54). The temperature and conductivity sensors were manufacturer-calibrated before deployment and the stated accuracy of the sensors are ±0.002 °C and ±0.0003 S m⁻¹, respectively.

A moored turbulence instrument cluster was deployed about 1.5 m beneath the ice base to observe the small-scale turbulent fluctuations in the ice shelf–ocean boundary layer⁴². Consisting of a Nobska Modular Acoustic Velocity Sensor (MAVS) differential acoustic travel-time 3D velocity sensor, an RBRcoda fast-response temperature sensor and an RBRconcerto inductive conductivity sensor, the turbulence instrument cluster was scheduled to operate in burst mode, sampling at 5 Hz for 15 min every 2 h. For each 15-min burst, the average temperature, conductivity and velocity values were received over an Iridium satellite link. In this study we use the mean values from 4,459 individual bursts collected between 23 January 2020 and 1 February 2021. An analysis of the full 5-Hz turbulence data awaits a future study. The stated uncertainty in the velocity components is 3 mm s⁻¹, whereas the stated uncertainty in the temperature and conductivity data are ±0.002 °C and ±0.0003 S m⁻¹, respectively. A small offset in the conductivity data caused by proximity effects associated with the inductive sensor was removed through regression against the CTD conductivity data.

Lateral gradients in temperature and salinity

CTD profiling beneath TEIS was largely carried out in three separate sessions with approximately 19 h between the first and second sessions and 47 h between the second and third sessions. Dividing the mean absolute difference in temperature and salinity between each session pair by the lateral advective distance over the time separating each pair (assuming a mean flow speed of 3 cm s⁻¹ during the CTD deployment period; Extended Data Fig. 4a) gives a mean lateral gradient of 2.7 × 10⁻³ ± 10⁻⁴ °C km⁻¹ for temperature and 1.1 × 10⁻³ ± 10⁻⁵ g kg⁻¹ km⁻¹ for salinity. The effect of tidal flow has been ignored in these calculations, as the tidal flow speeds are an order of magnitude weaker.

ApRES

An ApRES was established within 10 m of the borehole on 23 January 2020 and set to record a burst of 20 measurements once every 2 h. The data were recovered from the instrument on 27 December 2020. Data from four further ApRES deployments are also presented in this study (Fig. 1): a 5-month record from the first half of 2019 from an instrument deployed 360 m downstream of the borehole and data from three sites from 2020, contemporaneous with the borehole dataset. One of the 2020 sites was 1,310 m downstream of the borehole, another 1,340 m upstream of the borehole and, finally, one instrument was deployed on grounded ice 2,600 m upstream of the borehole, which tracked across the grounding line during 2020 and recorded every 3 h.

ApRES uses frequency-modulated continuous-wave modulation, with a chirp that scans from 200 to 400 MHz over a 1-s period. The measurements in each burst were checked for quality and then averaged. Each averaged burst was processed⁵⁵ to generate a radar return that preserves the signal phase.

By using both amplitude and phase, ApRES can monitor the changing distance between the antennas and the ice-shelf base with millimetre-scale precision. This raw Lagrangian ice-shelf thinning includes both the basal melt signal and the ice-column vertical strain that results from ice flow and snow compaction²¹. As well as the range to the ice base, the range to reflecting horizons within the ice column can be monitored and used to estimate the vertical strain within the ice as follows. The motion of internal reflecting horizons in any given depth interval can be found by cross-correlating the complex return for sequential (two-hourly) measurements for the entire time series. The vertical motion of the ice within the layer from one return to the next was derived from the phase of the cross-correlation and the reliability of that estimate was indicated by its amplitude. As all displacements are measured with respect to the antennas, the vertical displacement of any individual layer is the effect of the integrated strain in the ice above. The ice-shelf thinning rate, as measured from the antennas, is obtained by using a depth interval that tracks the return from the ice base.

For tracking the ice base to obtain the total thinning rate, an assumption in ApRES data processing is that, over the period of the time series, the topography of the ice base local to the radar does not change at length scales at or longer than the wavelength of the radar waves in ice, in this case, at length scales greater than about 0.5 m. This requirement was not met at the borehole site nor at various periods of the time series from the upstream sites. However, it was possible to use the first multiple echo, which is the result of the radar signal travelling from the transmit antenna to the ice base, back to the ice surface and then back to the ice base, before finally returning to the receive antenna. The range to the multiple is largely immune to local topographic evolution in the ice base, presumably because of the much larger effective footprint. The multiple is very much weaker than the first basal return, but its phase can be reliably tracked. For the downstream ApRES deployment, and those sections of the other deployments when the first basal return was not changing its form, the melt-rate time series from the first and second returns yielded a satisfactory match. Short-term variations in derived melt rate (3 to 5 days) were much stronger from the multiple, possibly resulting from snow-accumulation events.

Vertical profiles of vertical ice velocities, averaged across the entire time series, were calculated by dividing the return into 4-m layers, cross-correlating as described above and calculating the mean vertical velocity for each layer. From these profiles a depth interval was selected from the lower part of the ice column from which to calculate the integrated non-melt contribution to the thinning rate. The selection of the depth interval was based on the strength of the time-averaged correlations and, for the downstream site, the borehole site and the first upstream site, was 300–472 m, 304–436 m and 300–550 m, respectively. The strain from one measurement to the next was averaged across the depth interval and the time series differentiated and low-pass-filtered to provide a time series of vertical velocity variability at timescales of 5 days and longer. Although the strain rate in the ice column was expected to evolve slowly as the ice moved downstream, we assumed that non-tidal short-term variations would not be present. The vertical profiles of vertical velocity showed an approximately linear gradient through the selected depth intervals and an offset vertical velocity was determined by extrapolating that variation to the depth of the basal reflector. That offset was then added to the vertical velocities to yield the final non-melt contribution to the time series of ice-shelf thinning rate. The 2019 site, and the initially grounded site, were processed slightly differently. The internal reflections from near the ice base were good enough to allow a deep depth interval to be selected and cross-correlated to yield an integrated vertical strain directly. Intervals from 520 to 580 m and from 515 to 550 m were used for the 2019 and initially grounded sites, respectively, with a minor correction to accommodate a small, approximately linearly increasing strain rate near the base.

The final melt-rate time series was calculated by subtracting the non-melt contribution from the thinning rate and then low-pass-filtering at a cutoff of 15 days.

Icefin remotely operated underwater vehicle

Icefin was equipped with a Neil Brown Ocean Sensors conductivity–temperature (CT) sensor and a Valeport ultraP pressure sensor. The stated manufacturer accuracies are ±0.001 S m⁻¹, ±0.005 °C and 0.1 dbar for conductivity, temperature and pressure, respectively, which translate into uncertainties of ±0.008 g kg⁻¹ for S_A and ±0.018 °C for Θ. All sensors were factory calibrated before deployment and then cross-compared with the SBE 49 CTD profiler to remove offsets in conductivity and temperature of 0.0286 S m⁻¹ and 0.0236 °C. The CT sensor recorded at a frequency of 5 Hz, whereas the pressure sensor recorded at 1 Hz. Pressure measurements were interpolated to match the 5-Hz CT data. Hydrographic data were post-processed by removing outliers that exceeded more than two standard deviations from the mean, as well as data points collected when the vehicle speed was lower than 5 cm s⁻¹. A three-point weighted-mean filter was also applied to the conductivity and temperature data.

Ocean current speeds were measured using a LinkQuest NavQuest 600 Micro Doppler Velocity Log, which doubles as an ADCP. The ADCP provides measurements of the current speed in 2-m bins at a variable distance from the vehicle, controlled by gradients in the pitch, roll, heading and speed of the vehicle. Uncertainty in the current velocity is typically 1% of the vehicle’s velocity in its direction of travel. As Icefin travels at speeds ≤50 cm s⁻¹, the uncertainty in velocity recorded by the ADCP in the direction of travel is ≤5 mm s⁻¹. The uncertainty in velocities perpendicular to the direction of travel is typically much lower. The velocity data were recorded at a rate of 1 Hz and were post-processed by removing data points when the vehicle pitch or roll is greater than 30°. A 30-s running mean filter was applied to all data points and measurements were filtered for gradients greater than one standard deviation from the mean in vehicle speed, pitch, roll and individual bin velocity. Finally velocities were bin-averaged into 1-m depth bins and velocities were excluded if they exceeded one standard deviation of the mean for each bin.

Ship-based CTD profiles

A dual-sensors system based on a Sea-Bird 911 CTD was used for conductivity, temperature and pressure measurements outside the ice-shelf cavity from the RVIB Nathaniel B. Palmer in 2019 and 2020. Standard Sea-Bird software Seasave version 7.26.1.8 was used for data collection and conductivity cell thermal mass correction in 2019 and Seasave version 7.26.7.121 in 2020. Manufacturer-recommended values for cell thermal mass correction were used as follows: thermal anomaly amplitude, α = 0.03 and thermal anomaly time constant 1/β = 7.0. Water samples were taken from the CTD rosette and analysed using a Guildline Portasal salinometer to calibrate the primary and secondary conductivity sensors on the CTD profiler.

Density ratio and Turner angle

Double-diffusive convection occurs as a result of the difference in molecular diffusivities between salt and heat⁵⁶. Under Antarctic ice shelves, the presence of cold and fresh meltwater-laden waters above warm and salty modified Circumpolar Deep Water drives a double-diffusive process known as diffusive convection. Strong diffusive convection can lead to the formation of ‘diffusive staircases’, where well-mixed layers in temperature and salinity are separated by sharp interfaces^57,58. Diffusive convection can still occur without staircase formation however. Diffusive convection can exert a first-order control on the rate of ice-shelf basal melting^59,60.

The susceptibility of a water column to diffusive convection can be characterized through the density ratio

which measures the degree of compensation between temperature and salinity gradients in terms of their effect on density stratification. α is the thermal expansion coefficient, β the haline contraction coefficient, (frac{partial varTheta }{partial {rm{z}}}) the vertical gradient of conservative temperature and (frac{partial {S}_{{rm{A}}}}{partial text{z}}) is the vertical gradient of absolute salinity. A water column is susceptible to diffusive convection when R_ρ is between 0 and 1, with the strength of diffusive convection increasing as R_ρ approaches 1. The Turner angle

is related to the density ratio, in which arctan2 is the four-quadrant inverse tangent (tan⁻¹) and the water column is susceptible to diffusive convection when Tu is between −45° and −90°.

Glacial meltwater and subglacial discharge fractions

Meltwater fractions are calculated using the composite tracer method²⁷ and water-mass endmembers derived from the ship-based CTD profiles. In the absence of glacial meltwater or subglacial discharge, it is assumed that the ambient water column beneath TEIS would be composed exclusively of mCDW and WW that mix along a straight line between their corresponding endmembers⁶¹ (Fig. 2c and Extended Data Fig. 7). For each conservative temperature and absolute salinity observation, a composite tracer can be constructed

in which Θ_mCDW and ({S}_{{{rm{A}}}_{{rm{mCDW}}}}) are the conservative temperature and absolute salinity of the mCDW endmember, respectively, Θ_WW and ({S}_{{{rm{A}}}_{{rm{WW}}}}) are the conservative temperature and absolute salinity of the WW endmember, respectively, and Θ_ob and (,{S}_{{{rm{A}}}_{{rm{ob}}}}) are observed values of conservative temperature and absolute salinity, respectively. If a data point lies on the ambient mCDW–WW mixing line, ({psi }_{{{rm{o}}{rm{b}}}_{{rm{W}}{rm{W}}}}) is equal to zero. The value of ({psi }_{{{rm{o}}{rm{b}}}_{{rm{W}}{rm{W}}}}) will become non-zero, however, if glacial meltwater (MW) causes a data point to move off the ambient mixing line. The value of the composite tracer in pure MW is

in which Θ_MW and ({S}_{{{rm{A}}}_{{rm{MW}}}}) are the conservative temperature and absolute salinity of the MW endmember, respectively. The fraction of glacial meltwater (x_MW) present in the water column can then be calculated from

To quantify the variability in x_MW caused by uncertainty in the water-mass endmembers, 1,000 independent estimates of the meltwater fraction were made using a set of random endmember properties derived from the normal distribution described by the mean and standard deviation of each endmember property (Extended Data Table 2). The observed meltwater fraction is given by the mean of the 1,000 independent estimates, and the uncertainty is given by the standard error of the mean. In general, the uncertainty is two orders of magnitude smaller than the mean.

The conservative temperature and absolute salinity of the mCDW and WW endmembers were extracted from the ship-based CTD casts collected in front of TEIS (Fig. 2c). The properties of the WW endmember (Extended Data Table 2) are set to those found at the depth of the temperature minimum below the surface layer, whereas the mCDW endmember properties are set to those found at the depth of the temperature maximum. The MW endmember has an effective conservative temperature of −90.8 °C (ref. ²⁷) and an absolute salinity of 0 g kg⁻¹, whereas the conservative temperature of the subglacial discharge endmember is set to the pressure-dependent in situ freezing temperature for freshwater at the depth of the grounding line (−0.36 °C) with an absolute salinity of 0 g kg⁻¹. The endmember values are consistent with those used in previous studies^61,62.

Starting in September 2020, a persistent signal of fresh subglacial discharge (SD) appears in the hydrographic data. In Θ–S_A space, individual data points fall above the mCDW–MW mixing line, indicative of the presence of this fourth water mass (Extended Data Fig. 7). With only Θ and S_A available as tracers, it is not possible to solve for all water-mass fractions simultaneously, as the system is underdetermined. Instead, we have to make a necessary assumption that the influence of WW is negligible and that the water column is composed solely of a mix of mCDW, MW and SD. Although this practical assumption cannot be fully justified, WW is typically found above a depth of 400 m in the Amundsen Sea^61,63,64 and is therefore mostly excluded from the grounding-zone region owing to the depth of the ice base (Fig. 4). As a result, the impact of this assumption on the water-mass fractions is probably small. To determine the SD faction for data points that lie outside the mCDW–WW–MW mixing triangle, a composite tracer is constructed that is equal to zero for data points that lie along the mCDW–MW mixing line:

The value of this composite tracer in pure SD is

in which Θ_SD and ({S}_{{{rm{A}}}_{{rm{SD}}}}) are the conservative temperature and absolute salinity of the SD endmember, respectively, and the SD fraction can be calculated from

Similarly, the MW fraction for data points that lie outside the mCDW–WW–MW mixing triangle can be derived by constructing a composite tracer that is equal to zero for data points that lie along the mCDW–SD mixing line:

Taking the value of this composite tracer in pure MW as

the MW fraction can be calculated as

Ice base Ekman boundary layer

Assuming a momentum balance between a steady, uniform, geostrophic flow beneath a flat ice base and the frictional stress exerted by the ice base against this flow^65,66:

in which f is the Coriolis parameter, ρ₀ is the fluid density, p is pressure, ν is the kinematic eddy viscosity and u and v are the horizontal velocity components, then the vertical structure of the horizontal velocity components through the boundary layer are given by the canonical Ekman solution⁶⁶

in which u_g is the magnitude of the far-field geostrophic flow, z is boundary layer depth and d is the Ekman depth:

Using the root mean square error as a cost function, we fit equations (15) and (16) to the u and v boundary-layer velocity profiles from Icefin to determine the value of the eddy viscosity and the Ekman depth beneath TEIS.

Thermal driving and the three-equation model for basal melting

The rate of ice-shelf basal melting is controlled by the divergence of the sensible heat flux at the phase change interface⁴⁰

in which ρ is density, a_b the basal melt rate, L_i the latent heat of fusion of ice, c the specific heat capacity, κ_i the thermal diffusivity of ice, T the temperature and W is the vertical ocean velocity. The subscripts i, b and w refer to ice, ice–ocean boundary and ocean, respectively. The primes refer to turbulent fluctuations and the angled brackets to the time average. The first term on the right-hand side is the conductive heat flux into the ice, whereas the second term represents the vertical turbulent heat flux through the oceanic boundary layer. In the absence of direct turbulence measurements, or in regional or global models that do not resolve the vertical scales of the ice shelf–ocean boundary layer, the second term is quantified through a simple turbulence closure scheme that models the heat flux, (langle {W}^{{prime} }{T}_{w}^{{prime} }rangle ), as a product of the drag coefficient (C_d), a non-dimensional turbulent transfer coefficient for temperature (Γ_T), the horizontal ocean velocity (U), and the thermal driving (T_d), that is given as the difference between the in situ ocean temperature some distance from the ice–ocean boundary (T_w; in this case, the temperature at the ocean mooring deployed 1.5 m beneath the ice shelf) and the temperature at the ice–ocean boundary (T_b) that is assumed to be at the in situ freezing point at salinity S_b and pressure P_b:

The in situ freezing point at the ice–ocean boundary is given by the liquidus condition

in which λ₁, λ₂ and λ₃ are constants and the boundary salinity (S_b) is given by the salt balance at the phase-change interface:

in which Γ_S is the turbulent transfer coefficient for salt and the salinity of ice (Si) is taken to be zero for ice shelves. The turbulent transfer coefficients, Γ_T and Γ_S, assume that the thermal and saline diffusive sublayers (diffusive regions next to the ice base that are dominated by molecular-scale processes) are controlled exclusively by current shear and thin with increasing current velocity; however, there are known to be many ice shelf–ocean environments where this is not the case⁸. The system of equations described by equations (18)–(21) represents the canonical three-equation model for ice-shelf basal melting^40,67. It can be solved for the basal melt rate given observed ocean temperature, salinity and flow speed and physical constants in Extended Data Table 3 (refs. ^40,42).

Because thermal driving is defined as the difference between the ocean temperature recorded by the mooring and the freezing temperature at the ice–ocean boundary, its magnitude is sensitive to the distance between the ocean mooring and the ice base. Basal melting increases the distance between the ocean mooring and the ice base as a function of time, and this drives an apparent increase in thermal driving without a change in the source water mass as the mooring effectively descends into warmer water. If we assume that the temperature profile through the boundary layer is fixed with respect to time, we can use the observed CTD profiles to estimate that the apparent change in thermal driving owing to ice-base recession is about 0.2 °C, or roughly 57% of the observed change (Extended Data Fig. 2a).

Although the rate of basal melting beneath TEIS predicted by the three-equation melt-rate model is linearly related to the magnitude of the thermal driving, it is relatively insensitive to the apparent change in thermal driving owing to ice-base recession. Indeed, the mean predicted melt rate only falls by 8% when thermal driving is corrected for ice-base recession (Extended Data Fig. 2b) and the predicted melt rate remains an order of magnitude higher than the observed values. The small reduction in predicted melt rate is consistent with our assertion that basal melting beneath TEIS is limited by the strong stratification, weak flow speeds and the lack of shear-driven turbulence, rather than the amount of heat available in the boundary layer. Irrespective of ice-base recession, thermal driving continues to exceed 1.8 °C by the beginning of 2021 (Extended Data Fig. 2a), highlighting the substantial amount of heat that is available in the grounding-zone region to drive basal melting. The apparent increase in thermal driving owing to ice-base recession therefore has no impact on our conclusion that the three-equation melt-rate model is unable to accurately predict the melt rate beneath TEIS.