Significance Thwaites Glacier is one of the West Antarctica's most prominent, rapidly evolving, and potentially unstable contributors to global sea level rise. Uncertainty in the amount and spatial pattern of geothermal flux and melting beneath this glacier is a major limitation in predicting its future behavior and sea level contribution. In this paper, a combination of radar sounding and subglacial water routing is used to show that large areas at the base of Thwaites Glacier are actively melting in response to geothermal flux consistent with rift-associated magma migration and volcanism. This supports the hypothesis that heterogeneous geothermal flux and local magmatic processes could be critical factors in determining the future behavior of the West Antarctic Ice Sheet.

Abstract Heterogeneous hydrologic, lithologic, and geologic basal boundary conditions can exert strong control on the evolution, stability, and sea level contribution of marine ice sheets. Geothermal flux is one of the most dynamically critical ice sheet boundary conditions but is extremely difficult to constrain at the scale required to understand and predict the behavior of rapidly changing glaciers. This lack of observational constraint on geothermal flux is particularly problematic for the glacier catchments of the West Antarctic Ice Sheet within the low topography of the West Antarctic Rift System where geothermal fluxes are expected to be high, heterogeneous, and possibly transient. We use airborne radar sounding data with a subglacial water routing model to estimate the distribution of basal melting and geothermal flux beneath Thwaites Glacier, West Antarctica. We show that the Thwaites Glacier catchment has a minimum average geothermal flux of ∼114 ± 10 mW/m2 with areas of high flux exceeding 200 mW/m2 consistent with hypothesized rift-associated magmatic migration and volcanism. These areas of highest geothermal flux include the westernmost tributary of Thwaites Glacier adjacent to the subaerial Mount Takahe volcano and the upper reaches of the central tributary near the West Antarctic Ice Sheet Divide ice core drilling site.

Heterogeneous geothermal flux and subglacial volcanism have the potential to modulate ice sheet behavior and stability by providing a large, variable supply of meltwater to the subglacial water system, lubricating and accelerating the overlying ice (1, 2). However, the magnitude and spatial pattern of geothermal flux are extremely difficult to measure, and the catchment-scale constraints derived from seismic tomography (3) and satellite magnetometry (4) produce contradicting spatial patterns and cannot resolve geothermal features relevant to local ice sheet forcing. Despite strong evidence for magma migration (5) and volcanism (5⇓⇓–8) beneath the West Antarctic Ice Sheet (WAIS), the limitations of these heterogeneous estimates have led modeling studies to assume unrealistic spatially uniform geothermal flux distributions (9, 10). Accurate modeling of ice sheet contributions to sea level, site selection for ice core drilling, and enhanced understanding of ice–mantle interactions all require more accurate higher-resolution estimates of the spatial distribution of geothermal flux across critical glacier catchments than are currently available.

Thwaites Glacier is one of the largest, most rapidly changing glaciers on Earth, and its landward sloping bed reaches into the deep interior of the WAIS, making it a leading component in scenarios for rapid deglaciation (9, 11). In addition, the catchment of Thwaites Glacier (Fig. 1A) also lies within the West Antarctic Rift System, a potentially reactivated intracontinental extension zone of low topography where crustal thinning from distributed Cretaceous and narrow-mode Cenozoic rifting produces elevated geothermal flux (5, 6, 8, 14⇓⇓–17). Given the setting and configuration of its catchment, heterogeneous geothermal flux beneath Thwaites Glacier is likely a significant factor in local, regional, and continental-scale ice sheet stability. Thwaites Glacier has been observed by a catchment-wide airborne radar sounding survey (11). To date, the use of radar sounding data to constrain melt rates has been limited to the interpretation of bed echo strengths (6, 18) to infer basal water or radar layer drawdown to infer melted ice loss at the bed (19). However, the interpretation of layer drawdown relies on the existence, persistence, and interpretability of layers in radar sounding profiles as well as constrained accumulation rates (19). Further, the strength of bed echoes is affected by a combination of the material and geometric properties of the ice sheet and bed which introduce ambiguities in quantitative echo interpretation (18, 20⇓–22). Fortunately, the upstream portion of Thwaites Glacier is known to be underlain by a well-quantified subglacial water system of distributed canals (23). Distributed canals have relatively constant average depths (24), and their reflecting interfaces can be modeled as flat plates (23, 25). Therefore, geometrically corrected (18) relative bed echo strengths in the upstream region of Thwaites Glacier will be proportional to the areal coverage (25) and local flux of subglacial water. This specific knowledge of the subglacial interface can be used to overcome the limitations of radar bed echo interpretation and unambiguously establish meltwater quantities with well-bounded uncertainties.

Fig. 1. Subglacial hydrologic setting of Thwaites Glacier. (A) Bed topography of the West Antarctic Ice Sheet and Amundsen Sea Embayment (12). (B) Subglacial hydrologic potential (13) for a distributed water system in the upstream region of the Thwaites Glacier catchment (black boundary). (C) Collection of subglacial water routing models that best fit the observed radar bed echo strength distribution (Fig. 2A), where the darkness of grayscale cells is the number of models (out of 50) for which these cells drain at least 10 others upstream.

Methods In this analysis, we determine the mean and confidence interval uncertainties for englacial attenuation rates (26) [for both scattering and reflecting spreading geometries (18)] to produce maps of the mean (Fig. 2A) and range (Fig. 2B) of observed relative bed echo strengths. Because distributed water is in pressure equilibrium with the overlying ice, its routing will be determined by the subglacial hydrologic potential, calculated using radar-derived ice thickness and surface slope (12, 13) (Fig. 1B). We generate a collection of water routing models by adding noise (at the scale of gridding uncertainties) to the bed topography and selecting those routes that best fit the relative bed echo strengths using uniform melt (Fig. 1C). We use these routing models to determine the spatial distribution of melt required to reproduce the pattern of relative echo strengths (Fig. 2A). We then scale the relative melt distribution by the spatial average and variance of routed subglacial water (13, 27) using the total melt from an ice sheet model of the Thwaites Glacier catchment (9). This ice sheet model includes frictional heating, horizontal advection, and an assumed uniform geothermal flux (9). Finally, we subtract the net effect of friction and advection to estimate the geothermal flux required to produce the remaining melt (Fig. 3). Details are given in SI Methods. Fig. 2. Radar sounding bed echo strengths. (A) Mean estimate of observed relative radar bed echo strengths for the Thwaites Glacier catchment (black boundary) corrected for geometric spreading losses. (B) Range of estimates of corrected relative bed echo strengths. Minor banding is due to variations in aircraft height above the ice surface combined with the different geometric loss terms. Bed topography (12) contour interval for Antarctica is 180 m. Fig. 3. Minimum geothermal flux and basal melt values required to reproduce the observed relative bed echo strengths (Fig. 2A) with subglacial water routing models (13, 27) (Fig. 1C) using the total melt water from an ice sheet model for the upstream portion of the Thwaites Glacier catchment (9). The minimum average inferred flux is ∼114 ± 10 mW/m2. High-flux areas exceed 200 mW/m2. A indicates the Mount Takahe volcano. B indicates the WAIS Divide ice core drilling site. High-melt areas are indicated by C in the westernmost tributary, D adjacent to the Crary mountains, and E in the upper portion of the central tributaries (8). Triangles show areas where radar-inferred melt anomalies exceed those generated by ice dynamics (friction and advection) (9) and inferred geothermal flux exceeds 150 mW/m2 (dark magenta) and 200 mW/m2 (light magenta). Bed topography (12) contour interval for Antarctica is 180 m.

Results The upstream region of the Thwaites Glacier catchment contains several areas of strong relative bed echoes (Fig. 2A) that exceed the mean bed echo strength by significantly more than the uncertainty in those strengths (Fig. 2B). Because the water system in this portion of the catchment is composed of distributed canals (23), high echo strengths can be interpreted as indicating larger quantities of subglacial water. The relative basal melt distribution required to fit the observed bed echo strengths with subglacial water routing models shows that water routing explains some of the strong reflections (and inferred high water quantities) in the trunk. The distribution of melt and geothermal flux (Fig. 3) includes several regions with high melt that are closely related to rift structure and associated volcanism (7, 8). These include the entire westernmost tributary (Fig. 3, location C) that flanks Mount Takahe (Fig. 3, location A), a subaerial volcano active in the Quaternary (28, 29), and several high-flux areas across the catchment adjacent to topographic features that are hypothesized to be volcanic in origin (7, 8) (e.g., Fig. 3, locations D and E). We also observe high geothermal flux in the upper reaches of the central tributaries that are relatively close to the site of the WAIS Divide ice core (Fig. 3, location B), where unexpectedly high melt and geothermal flux have been estimated.* We estimate a minimum average geothermal flux value of about 114 mW/m2 with a notional uncertainty of about 10 mW/m2 for the Thwaites Glacier catchment with areas exceeding 200 mW/m2 (Fig. 3). These values are likely underestimates due to the low uniform geothermal flux value used in the ice sheet model (9) and the compensating effect of enhanced vertical advection of cold shallow ice in high-melt areas. Note that this latter effect also predicts a subtle gradient of underestimated flux from the interior to the trunk as fast flow and associated frictional melting increases.

Discussion Above, we use radar echo strengths to constrain a subglacial water routing model to estimate the pattern of basal melting and geothermal flux for the Thwaites Glacier catchment within the WAIS. The simplifying assumptions in this analysis are rooted in specific knowledge of the geometry of the subglacial water system in this area (23) and conservative treatment of radar echo strength uncertainties. Our results produce high melt values adjacent to known volcanoes and structures that are morphologically suggestive of volcanic origin (7, 8). We believe that both the magnitude and spatial pattern of geothermal flux we present reflect the geologic and glaciological reality of the Thwaites Glacier bed and that contrary to previous modeling (9), our results show regions of high geothermal flux that are in substantial agreement with levels inferred from the ice core drilling site near the ice divide for the Thwaites catchment.* This new approach provides both higher resolution and more geologically realistic boundary conditions for ice sheet modeling than previous estimates from remote sensing techniques (3, 4). These results also demonstrate an approach that can be applied to a wide variety of radar sounders (because it requires only platform stability and not absolute calibrated echo strengths) in areas known to host distributed subglacial water systems. Our results further suggest that the subglacial water system of Thwaites Glacier may be responding to heterogeneous and temporally variable basal melting driven by the evolution of rift-associated volcanism and support the hypothesis that both heterogeneous geothermal flux (6) and local magmatic processes (5) could be critical factors in determining the future behavior of the WAIS.

Acknowledgments Radar processing assistance was provided by S. Kempf, and radar interpretation was performed by J. DeSantos, A. Jones, E. Powell, and M. Williams. We thank I. Joughin for providing his modeled basal melt water values. We thank two anonymous reviewers for their helpful comments on the manuscript. We thank the National Science Foundation (NSF) (Grants PLR-0636724, PLR-0941678, and PLR-1043761), the National Aeronautics and Space Administration (Grant NNX08AN68G), and the G. Unger Vetlesen Foundation for supporting this work. D.M.S. received support from an NSF Graduate Research Fellowship, a University of Texas Recruiting Fellowship, and the University of Texas Institute for Geophysics (UTIG) Gale White Fellowship. This is UTIG contribution 2711.

Footnotes Author contributions: D.M.S. designed research; D.M.S. performed research; D.M.S. contributed new reagents/analytic tools; D.M.S., D.D.B., D.A.Y., and E.Q. analyzed data; and D.M.S., D.D.B., D.A.Y., and E.Q. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

↵*Clow GD, Cuffey K, Waddington E, American Geophysical Union Fall Meeting, December 3–7, 2012, abstr C31A-0577.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1405184111/-/DCSupplemental.