Estimating effective pressures in active subglacial lakes with ICESat-2 satellite altimetry

TL;DR

This study uses satellite data to estimate pressures at the base of glaciers, revealing insights into how water affects ice movement.

Contribution

A new method is developed to estimate effective pressures in subglacial lakes using ICESat-2 altimetry data.

Findings

Effective pressure deviations are typically small, often less than 10 kPa.

Large effective pressures, exceeding 100 kPa, were observed in some subglacial lakes.

The method provides new data to improve subglacial hydrology models.

Abstract

The difference between the ice and water pressures, or the effective pressure, influences water flow and sliding at the ice-bed interface. Effective pressure is typically quantified with subglacial hydrology models because direct measurements of the subglacial environment are sparse. Active subglacial lakes provide an opportunity to constrain effective pressures with altimetry because subglacial water-volume changes manifest at the ice-sheet surface as elevation-change anomalies. Here, we develop a method for estimating effective pressures from altimetry data above active subglacial lakes. We synthesise a previous theory of subglacial lake effective pressure with an altimetry-based inverse method that relates elevation-change data to water-volume changes. We apply the method to elevation-change data from NASA’s ICESat-2 satellite altimetry mission over several active lakes in Antarctica. We find that deviations from flotation (zero effective pressure) are typically a negligible fraction of the overburden (e.g., 10 kPa), although larger deviations can arise when the ice viscosity is large. For example, effective pressures over subglacial lake Byrds10 in East Antarctica locally reached magnitudes on the order of the tensile strength of glacier ice (e.g., over 100 kPa). These effective pressure estimates can constrain subglacial hydrology models in regions with active subglacial lakes and provide new insights into glacier-bed dynamics.

Figures

a](#fig1)). The lower surface of the ice is denoted by $[eqn]$ . We assume that ice deforms as a viscous fluid according to the incompressible Stokes equations with a linear sliding law at the base and a stress-free condition at the ice-air interface (Stubblefield and others, [2023a](#ref55)). For consistency with the linearised framework, we assume a uniform ice viscosity $[eqn]$ , basal drag coefficient $[eqn]$ and ice thickness $[eqn]$ in the vicinity of the lake. The assumptions of a linearised framework and uniform material properties are adopted because subglacial lakes are expected to o

a](#fig3),[b](#fig3)).Figure 3.(a,b) Functions $[eqn]$ and $[eqn]$ (Eqns. 15-16) that determine the Fourier-transformed effective pressure $[eqn]$ (Eqn. 18) as functions of the scaled wavevector magnitude $[eqn]$ for different values of the nondimensional parameter $[eqn]$ . (c,d) Functions $[eqn]$ and $[eqn]$ (Eqns. 23-24) that determine the Fourier-transformed elevation-change anomaly $[eqn]$ (Eqn. 22) for different values of $[eqn]$ . All functions are nondimensional.

a](#fig4)). The out-of-phase component $[eqn]$ increases at slower oscillation frequencies ([Fig. 4b](#fig4)). The spectra of the effective pressure and basal vertical velocity can be positively correlated, negatively correlated or uncorrelated, depending on the oscillation frequency and spatial wavenumber ([Fig. 4c](#fig4)).Figure 4.(a) In phase component $[eqn]$ and (b) out-of-phase component $[eqn]$ of the effective pressure spectrum (28) for different oscillation frequencies $[eqn]$ . The nondimensional parameters are set to $[eqn]$ and $[eqn]$ . (c) Effective pressure spectrum versus vert

b](#fig5)) and the basal vertical velocity ([Fig. 5c](#fig5)), but can be more strongly correlated with one over the other, depending on the parameter values. For these parameters, the mean effective pressure is negatively correlated with the elevation-change anomaly over time ([Fig. 5a](#fig5)). Due to the combined influence of the basal velocity and elevation change, positive and negative values of the effective pressure can exist within the lake boundary simultaneously. For example, a ring-shaped region of negative effective pressure forms near the boundary of the lake at the start of the f

b](#fig6)-[c](#fig6)). The effective pressure distribution has a peak around $[eqn]$ kPa during draining stages ( $[eqn]$ ), accompanied by a long tail of negative values ([Fig. 6b](#fig6)). During filling stages ( $[eqn]$ ), the effective pressure distribution has a peak around $[eqn]$ kPa that is accompanied by a long tail of positive values ([Fig. 6c](#fig6)). For these parameters ( $[eqn]$ and $[eqn]$ ), the effective pressure falls within $[eqn]$ kPa, which is a relatively small deviation from flotation ( $[eqn]$ ).Figure 6.(a) Basal vertical velocity $[eqn]$ versus the effective pressure

c](#fig7)). Maps of the elevation-change anomaly, basal vertical velocity and effective pressure during the filling stage highlight the spatial variability in the system ([Fig. 7b](#fig7),[d](#fig7),[f](#fig7)). We find that the effective pressure has both regions of positive and negative effective pressure, which results from the combined influence of the basal vertical velocity and elevation change ([Fig. 7f](#fig7)). In particular, the centre of the lake has a positive effective pressure that is surrounded by a ring of negative effective pressure that forms as the lake refills. The same typ

d](#fig8)). These complex patterns arise from thickness changes as the ice dynamically adjusts after the drainage event (Appendix B). The effective pressure showed a larger spatial extent than the altimetry-derived lake boundary during the draining event. Mac1 showed a second subsidence event with additional oscillations in the effective pressure during 2023–24.Figure 8.Elevation change ( $[eqn]$ , basal vertical velocity inversion ( $[eqn]$ ) and effective pressure ( $[eqn]$ ) for Mac1. (a) Time series of the mean value of the elevation change over the lake. (b) Map–plane contour plot of the

a](#fig9)). The effective pressure shows a similar spatial extent to the elevation-change anomaly, although the spatial footprint appears to have changed significantly since the ICESat-era drainage event (Smith and others, [2009](#ref50); Siegfried and Fricker, [2018](#ref48); Lindzey and others, [2020](#ref27)). We revisit this discrepancy between the ICESat-era outline and the current boundary in the discussion.Figure 9.Elevation change ( $[eqn]$ , basal vertical velocity inversion ( $[eqn]$ ) and effective pressure ( $[eqn]$ ) for David_s1_. (a) Time series of the mean value of the elevatio

a](#fig10)). The draining is associated with an elevated mean effective pressure within the lake. However, there is a large region of negative effective pressure that forms near the centre of the lake ([Fig. 10f](#fig10)). Negative effective pressures arise when tensile viscous stresses dominate Eqn. (18) during draining events ( $[eqn]$ ) with higher ice viscosity. The spatial extent of the effective pressure is much larger than the elevation-change anomaly due to these high tensile viscous stresses, which shows that subglacial pressure perturbations can extend beyond altimetry-derived lake o

a](#fig1)). The lower surface of the ice is denoted by $[eqn]$ . We assume that ice deforms as a viscous fluid according to the incompressible Stokes equations with a linear sliding law at the base and a stress-free condition at the ice-air interface (Stubblefield and others, [2023a](#ref55)). For consistency with the linearised framework, we assume a uniform ice viscosity $[eqn]$ , basal drag coefficient $[eqn]$ and ice thickness $[eqn]$ in the vicinity of the lake. The assumptions of a linearised framework and uniform material properties are adopted because subglacial lakes are expected to o

a](#fig10)). The draining is associated with an elevated mean effective pressure within the lake. However, there is a large region of negative effective pressure that forms near the centre of the lake ([Fig. 10f](#fig10)). Negative effective pressures arise when tensile viscous stresses dominate Eqn. (18) during draining events ( $[eqn]$ ) with higher ice viscosity. The spatial extent of the effective pressure is much larger than the elevation-change anomaly due to these high tensile viscous stresses, which shows that subglacial pressure perturbations can extend beyond altimetry-derived lake o

a](#fig3),[b](#fig3)).Figure 3.(a,b) Functions $[eqn]$ and $[eqn]$ (Eqns. 15-16) that determine the Fourier-transformed effective pressure $[eqn]$ (Eqn. 18) as functions of the scaled wavevector magnitude $[eqn]$ for different values of the nondimensional parameter $[eqn]$ . (c,d) Functions $[eqn]$ and $[eqn]$ (Eqns. 23-24) that determine the Fourier-transformed elevation-change anomaly $[eqn]$ (Eqn. 22) for different values of $[eqn]$ . All functions are nondimensional.

a](#fig4)). The out-of-phase component $[eqn]$ increases at slower oscillation frequencies ([Fig. 4b](#fig4)). The spectra of the effective pressure and basal vertical velocity can be positively correlated, negatively correlated or uncorrelated, depending on the oscillation frequency and spatial wavenumber ([Fig. 4c](#fig4)).Figure 4.(a) In phase component $[eqn]$ and (b) out-of-phase component $[eqn]$ of the effective pressure spectrum (28) for different oscillation frequencies $[eqn]$ . The nondimensional parameters are set to $[eqn]$ and $[eqn]$ . (c) Effective pressure spectrum versus vert

b](#fig5)) and the basal vertical velocity ([Fig. 5c](#fig5)), but can be more strongly correlated with one over the other, depending on the parameter values. For these parameters, the mean effective pressure is negatively correlated with the elevation-change anomaly over time ([Fig. 5a](#fig5)). Due to the combined influence of the basal velocity and elevation change, positive and negative values of the effective pressure can exist within the lake boundary simultaneously. For example, a ring-shaped region of negative effective pressure forms near the boundary of the lake at the start of the f

Full text

Introduction

Subglacial water flow is regulated by the effective pressure, the difference between the ice and water pressures at the glacier base (Röthlisberger, 1972; Shreve, 1972). The effective pressure controls rates of creep closure in subglacial drainage elements such as channels and cavities, and thereby influences the volumetric discharge by restricting water flow (Nye, 1976; Flowers, 2015). Likewise, the magnitude and direction of water flow are determined by hydraulic potential gradients, which depend on gradients in the effective pressure (Hewitt, 2011). The effective pressure also influences the frictional behaviour at the ice–bed interface and modulates ice–flow speeds (Lliboutry, 1968; Bindschadler, 1983; Schoof, 2005; Zoet and Iverson, 2020). While the effective pressure is a fundamental variable that controls subglacial water flow and sliding at the ice–bed interface, direct measurements obtained via borehole drilling are sparse (Iken and others, 1993; Fountain, 1994; Hubbard and others, 1995; Engelhardt and Kamb, 1997; Meierbachtol and others, 2013; Andrews and others, 2014; Rada and Schoof, 2018).

Subglacial lakes that are observed to episodically fill and drain, often called ‘active’ lakes, present an opportunity for constraining effective pressures with altimetry data because subglacial water-volume changes manifest at the ice-sheet surface as elevation-change anomalies (Fricker and others, 2007; 2010; Fricker and Scambos, 2009; Smith and others, 2009). In particular, the coupling between the effective pressure in a subglacial lake with the surrounding drainage system drives the water-volume oscillations that are expressed at the ice-sheet surface (Evatt and others, 2006; Fowler, 2009; Stubblefield and others, 2019). The ICESat (2003–09) and CryoSat-2 (2010 to present) satellite altimetry missions have detected over one hundred active lakes beneath the Antarctic Ice Sheet, while a smaller number have been found in Greenland, Iceland and various mountain glaciers (Smith and others, 2009; Wright and Siegert, 2012; Siegfried and Fricker, 2018; Livingstone and others, 2022). NASA’s ICESat-2 satellite altimetry mission (2018 to present) has allowed for continued detection and monitoring of active subglacial lakes at high spatial and temporal resolution (Neckel and others, 2021; Siegfried and Fricker, 2021; Livingstone and others, 2022; Fan and others, 2023; Freer and others, 2024; Gray and others, 2024a; 2024b).

Previous modelling work quantified subglacial lake effective pressures with a finite element method (Stubblefield and others, 2021b). A limitation of the finite element method is that the mean effective pressure in the lake is determined numerically with a Lagrange multiplier, which only furnishes an indirect relation between lake activity and viscous ice flow. An alternative approach was developed to estimate subglacial water-volume changes with an inverse method that accounts for the effects of viscous ice flow on surface-elevation changes (Stubblefield and others, 2023a). The inverse method assumes that subglacial lake oscillations represent small perturbations to an ice-flow state that is described by auxiliary parameters such as the ice thickness, viscosity, flow speed and basal sliding coefficient (cf. Gudmundsson, 2003).

In this study, we develop a method for estimating effective pressures from elevation-change anomalies above active subglacial lakes by synthesising the previous modelling approaches (Stubblefield and others, 2021b; 2023a). First, we derive general expressions for the effective pressure that depend on the ice–surface elevation and viscous ice flow. Then, we relate the effective pressure to the elevation change and basal vertical velocity with a linearised model (Stubblefield and others, 2023a). We present semi-analytical and synthetic examples to illustrate the basic behaviour of the method. Finally, we apply the method to a collection of subglacial lakes in Antarctica that have shown activity during the ICESat-2 era.

Model derivation

In this section, we first derive general formulas for the effective pressure, which is defined by

$$[eqn]$$

where $[eqn]$ is the ice pressure and $[eqn]$ is the water pressure. We refer to the condition $[eqn]$ —when the ice and water pressures are balanced—as ‘flotation’. Then, we use a linearised (small perturbation) model to relate the effective pressure in a subglacial lake to the basal vertical velocity anomaly $[eqn]$ that produces the observed surface-elevation change $[eqn]$ (Figs. 1,2).Figure 1.(a) Sketch of a subglacial lake in cross-section highlighting the elevation-change anomaly $[eqn]$ , ice-base elevation $[eqn]$ , and effective pressure $[eqn]$ . The ice layer over the lake is characterised by the thickness $[eqn]$ and viscosity $[eqn]$ , while the ice-bed interface is characterised by the basal drag coefficient $[eqn]$ . (b) Map-view sketch showing the lake area $[eqn]$ , lake boundary $[eqn]$ and normal vector $[eqn]$ .Figure 2.Map of ICESat-2 elevation-change data from Antarctica (ATL15 gridded product; Smith and others, 2024) with insets showing anomalies over subglacial lakes Mac1 (MacAyeal Ice Stream), Mercer Subglacial Lake (MSL; Mercer Ice Stream), Byrd_s10_ (Byrd Glacier) and David_s1_ (David Glacier). The map-plane $[eqn]$ coordinates in the ATL15 dataset correspond to the Antarctic Polar Stereographic Projection (EPSG:3031). Lake outlines from Siegfried and Fricker (2018) are shown in silver on the insets.

We represent the ice-surface elevation $[eqn]$ via $[eqn]$ where $[eqn]$ is the background ice thickness and $[eqn]$ is the elevation-change anomaly over the lake (Fig. 1a). The lower surface of the ice is denoted by $[eqn]$ . We assume that ice deforms as a viscous fluid according to the incompressible Stokes equations with a linear sliding law at the base and a stress-free condition at the ice-air interface (Stubblefield and others, 2023a). For consistency with the linearised framework, we assume a uniform ice viscosity $[eqn]$ , basal drag coefficient $[eqn]$ and ice thickness $[eqn]$ in the vicinity of the lake. The assumptions of a linearised framework and uniform material properties are adopted because subglacial lakes are expected to only generate small changes in ice flow relative to a given background state (Stubblefield and others, 2023a). The primary limitation of this approach is that the basal drag coefficient does not depend on the time-varying effective pressure in the subglacial lake, but is instead assumed to be constant. The relevant components of the linearised model are summarised below while precise details on the derivation are outlined in Appendix A.

Effective pressure formulas

2.1.

Hydrostatic approximation

2.1.1.

We derive an expression for the effective pressure in a subglacial lake by assuming that the water pressure follows a hydrostatic gradient,

$$[eqn]$$

where $[eqn]$ kg m $[eqn]$ is the density of water and $[eqn]$ is gravitational acceleration with magnitude $[eqn]$ m s $[eqn]$ . The water pressure is not known a priori because it depends on the stresses in the overlying ice column (Stubblefield and others, 2021b). Therefore, the water pressure is determined up to a constant that depends only on the map-plane coordinates $[eqn]$ , which we resolve by averaging over the surface area of the lake. In this way, we obtain a formula for the water pressure at the ice–water interface ( $[eqn]$ ),

$$[eqn]$$

where bars denote spatial means over the surface area of the lake (Stubblefield and others, 2021b). Taking the difference of the ice pressure $[eqn]$ with (3), we obtain an expression for the effective pressure $[eqn]$ over the subglacial lake,

$$[eqn]$$ $$[eqn]$$

where we have defined the density difference, $[eqn]$ , and the difference between the ice pressure and cryostatic pressure, $[eqn]$ . Two challenges arise in applying the formula (4) to estimate effective pressures in subglacial lakes. First, the mean effective pressure ( $[eqn]$ ) and the deviation of the ice pressure from the cryostatic pressure ( $[eqn]$ ) are undetermined at this stage, requiring incorporation of the ice-flow dynamics (Stubblefield and others, 2021b). Second, the motion of the ice–water interface $[eqn]$ does not necessarily correspond to the motion of the ice surface $[eqn]$ (Stubblefield and others, 2021a; 2023a). We resolve these challenges by directly incorporating the effects of ice flow into the method below.

The hydrostatic assumption implies that subglacial lakes are still bodies of water that tend to capture water from surrounding areas. Gradients in the subglacial lake effective pressure (4) balance the ‘background’ hydraulic gradient of the subglacial drainage system in the limit of a cryostatic ice pressure,

$$[eqn]$$

which implies that the water flux is zero and that the lake corresponds to a minima in the hydraulic potential (e.g., Hewitt, 2011). In subglacial drainage systems, water flows due to deviations from the background hydraulic gradient, so we cannot estimate effective pressures with the formula (4) outside of the lake boundary. Therefore, we will restrict formulas to the interior of the lake below by introducing an indicator function $[eqn]$ , defined by

$$[eqn]$$

where $[eqn]$ denotes the lake area (Fig. 1). While subglacial lake shorelines can migrate during volume-change events, we assume a fixed lake boundary $[eqn]$ for simplicity (Siegfried and Fricker, 2021; Stubblefield and others, 2021a; 2021b). We revisit this assumption in the discussion. Extending the current implementation to allow for time-varying boundaries (evolving $[eqn]$ ) would be straightforward (see Data Availability statement for code repository).

Small-slope approximation

2.1.2.

We derive a direct relation between the effective pressure and ice dynamics by considering the normal stress at the ice-water interface. Assuming that the slope of the basal surface $[eqn]$ is small over the subglacial lake, the normal stress $[eqn]$ at the base $[eqn]$ is given by

$$[eqn]$$

where $[eqn]$ is the vertical component of the ice velocity (Stubblefield and others, 2023b). The normal stress equals the water pressure at the ice-water boundary, which implies that the effective pressure is given by

$$[eqn]$$

at $[eqn]$ . Incompressibility leads to an expression in terms of the horizontal flow divergence,

$$[eqn]$$

where $[eqn]$ is the horizontal velocity. Eqn. (10) shows that the effective pressure is controlled by viscous flow towards the lake, reflecting the concept of creep closure. The ice is floating ( $[eqn]$ ) when the horizontal flow divergence vanishes over the lake as in the case of rigid motion or cryostatic balance. An approximate analysis using (10) and mass continuity shows that the effective pressure is related to dynamic thickness changes through the logarithmic strain rate of the ice column (Appendix B). The effective pressure formulas (9)-(10) are general in the sense that they hold everywhere there is an ice–water interface with small slope. Below, we combine the effective pressure formula (9) that is based on viscous ice stresses with the hydrostatic formula (4) to estimate effective pressures in active subglacial lakes.

We integrate Eqn. (10) over the lake area $[eqn]$ and use the divergence theorem to obtain an expression for the mean effective pressure,

$$[eqn]$$

where $[eqn]$ is the boundary of the lake, $[eqn]$ is the surface area of the lake, and $[eqn]$ denotes an inward-pointing unit normal to the lake boundary (Fig. 1b). In previous work, the mean effective pressure $[eqn]$ was determined with a Lagrange multiplier, which provided an indirect relation between the effective pressure and ice flow (Stubblefield and others, 2021b). Eqn. (11) instead provides a direct relation between the mean effective pressure in a subglacial lake and viscous flow towards the lake.

As a simple example, Eqn. (11) implies that a lake with circular boundary has a mean effective pressure of $[eqn]$ where $[eqn]$ is the radius of the lake and $[eqn]$ is the average inflow speed. This simple expression shows the sensitivity of the mean effective pressure to the ice viscosity. In the limit of a large lake area ( $[eqn]$ ), the mean effective pressure approaches zero analogous to a floating ice shelf (Stubblefield and others, 2023b). On the other hand, the mean effective pressure can become large in the limit of a small lake area ( $[eqn]$ ), although the inflow speed $[eqn]$ may also diminish in this limit. While lakes that are smaller than the ice thickness are difficult to detect with altimetry-based methods (Stubblefield and others, 2021a), many small lakes that could have elevated effective pressures have been detected with ice-penetrating radar (Wright and Siegert, 2012; MacKie and others, 2020).

Perturbation formulas

2.2.

Next, we introduce a linearised (small perturbation) approach for estimating the effective pressure $[eqn]$ from elevation changes (Stubblefield and others, 2023a). We assume a uniform ice thickness $[eqn]$ in the base state where the elevation of the ice-water interface $[eqn]$ corresponds to $[eqn]$ (Fig. 1). The viscosity $[eqn]$ and basal drag $[eqn]$ are also assumed to be uniform in the vicinity of the lake. The limitations of these assumptions have been discussed and tested against a fully nonlinear model in previous work (Stubblefield and others, 2023a). Formally, we consider $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ to be perturbations to a cryostatic base state characterised by the ice pressure $[eqn]$ in addition to vanishing vertical velocity and effective pressure fields. In the base state, Eqn. (8) implies that the mean water pressure is cryostatic ( $[eqn]$ ) while the mean effective pressure is zero ( $[eqn]$ ). Likewise, the parameter $[eqn]$ vanishes in the cryostatic base state. We test the assumption that the background state is at flotation ( $[eqn]$ ) below by considering dynamic thickness changes away from the lake (Table 1; Appendix B).

We linearise the effective pressure around the base-state lower-surface elevation ( $[eqn]$ ) via

$$[eqn]$$

where the second term represents the directional derivative of $[eqn]$ at $[eqn]$ in the direction of the ice–water interface perturbation $[eqn]$ in the cryostatic base state ( $[eqn]$ . We find from Eqn. (4) that the directional derivative is $[eqn]$ . On the other hand, we resolve the first term in (12) with the expression (9) evaluated at $[eqn]$ .

With these considerations, the linearised expression (12) becomes

$$[eqn]$$

The first term in (13) captures the effects of viscous ice flow while the second term represents perturbations to the hydrostatic normal stress due to movement of the ice base (cf. Stubblefield and others, 2023b). The second term (directional derivative) in (13) has been truncated with the indicator function $[eqn]$ because variations in the effective pressure cannot be constrained outside of the lake boundary with the hydrostatic formula. In particular, we recover the expression for the mean effective pressure (11) from Eqn. (13) because the mean of the hydrostatic term over the lake area vanishes.

We solve the linearised problem by deriving a formula for $[eqn]$ , the Fourier transform of the effective pressure (13) with respect to the map-plane coordinates $[eqn]$ . In Appendix A, we show that

$$[eqn]$$

where we have defined the functions

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ being the wavevector magnitude normalised by the background ice thickness $[eqn]$ . In Eqns. (15)-(16), we have introduced a nondimensional parameter that arises from the sliding law,

$$[eqn]$$

which relates all of the fundamental parameters that characterise the base state. The magnitudes of the functions $[eqn]$ and $[eqn]$ decay at larger values of the sliding parameter $[eqn]$ , which leads to diminished viscous stresses at the base (Fig. 3a,b).Figure 3.(a,b) Functions $[eqn]$ and $[eqn]$ (Eqns. 15-16) that determine the Fourier-transformed effective pressure $[eqn]$ (Eqn. 18) as functions of the scaled wavevector magnitude $[eqn]$ for different values of the nondimensional parameter $[eqn]$ . (c,d) Functions $[eqn]$ and $[eqn]$ (Eqns. 23-24) that determine the Fourier-transformed elevation-change anomaly $[eqn]$ (Eqn. 22) for different values of $[eqn]$ . All functions are nondimensional.

We combine Eqns. (13) and (14) to obtain a formula for the Fourier-transformed effective pressure,

$$[eqn]$$

where we have assumed that the basal surface perturbation vanishes outside the lake boundary ( $[eqn]$ ). The potential correlation between the effective pressure and the basal vertical velocity in (18) is analogous to the bending of an elastic beam (Evatt and Fowler, 2007). In particular, upward motion ( $[eqn]$ ) can result in compression ( $[eqn]$ ) at the base and downward motion can result in tension, although the precise level of correlation depends on the magnitude of the viscous term relative to the elevation-dependent terms. In Appendix C, we show that the opposite behaviour arises where the basal vertical velocity is proportional to the negative effective pressure in the limit of steady creep, which is consistent with commonly used creep closure laws for subglacial channels and other drainage elements (Nye, 1976; Hewitt, 2011). A simple, alternative interpretation relates the effective pressure to dynamic thickness changes, which are calculated implicitly with this method (Appendix B).

Following previous work, we obtain the basal vertical velocity $[eqn]$ from elevation-change data $[eqn]$ with an inverse method that follows the linearised approach taken herein (Stubblefield and others, 2023a). The ice-base elevation $[eqn]$ is determined up to an initial condition $[eqn]$ by

$$[eqn]$$

where we have neglected a small advective component from the background flow under the assumption that the ice–water interface motion is caused by the basal vertical velocity from water-volume changes (Stubblefield and others, 2023a). We can calculate the effective pressure from Eqn. (18) after specifying the initial conditions and auxiliary parameters ( $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ ) that the inverse method requires, where $[eqn]$ is the mean ice surface velocity in the base state (Stubblefield and others, 2023a). We provide a simplified analysis of Eqn. (18) in the following section for illustration, while a complete description of the inverse method is provided by Stubblefield and others (2023a).

Analysis

Before applying the estimation method to Antarctic subglacial lakes (Fig. 2), we first analyse the effective pressure $[eqn]$ (Eqn. 18) in relation to the basal vertical velocity $[eqn]$ . We then explore a synthetic example that falls within the parameter regime of the Antarctic subglacial lakes.

Scaling

3.1.

First, we scale Eqn. (18) to facilitate the analysis below. For simplicity, we refrain from renaming the nondimensional variables. We let $[eqn]$ be the elevation anomaly scale (1 m) and $[eqn]$ be the observational time scale (1 yr). We scale the map-plane coordinates $[eqn]$ by $[eqn]$ , the effective pressure by $[eqn]$ , and the vertical velocity by $[eqn]$ . With this scaling, Eqn. (18) becomes

$$[eqn]$$

where we have defined the flotation factor $[eqn]$ . The parameter $[eqn]$ is defined by the ratio

$$[eqn]$$

where $[eqn]$ is the timescale for viscous relaxation (decay) of topography (Turcotte and Schubert, 2014; Stubblefield and others, 2021a; 2023a). The scaling in Eqn. (20) shows that the precise definition of the lake boundary is of lesser importance in determining the effective pressure since it only manifests in the smaller, hydrostatic term that is multiplied by the flotation factor. However, careful consideration of the lake boundary is important for determining where the underlying approximations are valid. Additionally, any suspended sediment in the water column could potentially increase the fluid density, making the hydrostatic term more important.

Closure relations

3.2.

Next, we examine relations between the effective pressure and basal vertical velocity, which are related in subglacial hydrology models through creep closure laws (e.g., Nye, 1976; Evatt, 2015; Meyer and others, 2016). We outline how our formulation reduces to a similar form as the previously derived closure laws under certain simplifying assumptions in Appendix C. For simplicity, we assume here that there is no ice advection in the background state ( $[eqn]$ ) and no initial surface perturbations ( $[eqn]$ at $[eqn]$ ). The general case has been previously covered (Stubblefield and others, 2023a). Under these assumptions, the elevation change is related to the basal vertical velocity by

$$[eqn]$$

where $[eqn]$ denotes convolution over time (Appendix A).

In Eqn. (22), $[eqn]$ describes viscous decay of surface topography (Fig. 3c) while $[eqn]$ is a base-to-surface transfer function (Fig. 3d). These functions depend on the wavenumber $[eqn]$ and nondimensional sliding parameter $[eqn]$ via

$$[eqn]$$ $$[eqn]$$

The functions (23)-(24) are rewritten slightly from previous work due to the scaling adopted herein (Stubblefield and others, 2023a). The elevation-change formula (22) provides the basis for the least-squares inverse method for obtaining the basal vertical velocity (Stubblefield and others, 2023a). Substituting the relations (19) and (22) into Eqn. (20), we obtain

$$[eqn]$$

which provides a direct relation between the effective pressure and basal vertical velocity. In the limit of steady creep, Eqn. (25) reduces to a similar form relative to previous closure laws, where the closure rate is proportional to the negative effective pressure (Appendix C).

Sinusoidal oscillations

3.3.

As a semi-analytical example, we assume that the basal vertical velocity oscillates in time according to

$$[eqn]$$

where $[eqn]$ denotes the spatial pattern of the basal vertical velocity anomaly. We assume that the lake is radially symmetric with radius $[eqn]$ and surface area $[eqn]$ , which implies that the indicator function $[eqn]$ transforms to

$$[eqn]$$

where $[eqn]$ is the order-one Bessel function of the first kind.

We insert (26)-(27) into Eqn. (25), calculate the integrals and neglect an exponential decay term (proportional to $[eqn]$ ) to obtain the long-time behaviour ( $[eqn]$ )

$$[eqn]$$ $$[eqn]$$

We plot the in-phase component $[eqn]$ and out-of-phase component $[eqn]$ for various oscillation frequencies $[eqn]$ in Figure 4. The in-phase component $[eqn]$ decreases as the oscillation frequency $[eqn]$ decreases. For all oscillation frequencies, the in-phase component decays to zero at short wavelengths and in the long-wavelength limit (Fig. 4a). The out-of-phase component $[eqn]$ increases at slower oscillation frequencies (Fig. 4b). The spectra of the effective pressure and basal vertical velocity can be positively correlated, negatively correlated or uncorrelated, depending on the oscillation frequency and spatial wavenumber (Fig. 4c).Figure 4.(a) In phase component $[eqn]$ and (b) out-of-phase component $[eqn]$ of the effective pressure spectrum (28) for different oscillation frequencies $[eqn]$ . The nondimensional parameters are set to $[eqn]$ and $[eqn]$ . (c) Effective pressure spectrum versus vertical velocity spectrum ( $[eqn]$ component), normalised by the spectral amplitude of the vertical velocity $[eqn]$ . For this value of $[eqn]$ , we set the long-wavelength term to $[eqn]$ in Eqn. (28), which corresponds to the Gaussian-shaped anomaly in Figure 5.

Synthetic example

3.4.

To translate the preceding analysis to physical space, we consider a synthetic example with a Gaussian-shaped basal vertical velocity that exhibits sinusoidal oscillations in time (Fig. 5). In dimensional terms, the basal vertical velocity $[eqn]$ has a width of $[eqn]$ 20 km, amplitude of $[eqn]$ m yr $[eqn]$ and oscillation period of 5 yr. While the relation between the effective pressure $[eqn]$ and basal vertical velocity $[eqn]$ depends on the nondimensional parameters, we consider the values $[eqn]$ and $[eqn]$ , which fall within the range of values for the Antarctic subglacial lakes discussed below (Table 1). The effective pressure is influenced by both the elevation-change anomaly (Fig. 5b) and the basal vertical velocity (Fig. 5c), but can be more strongly correlated with one over the other, depending on the parameter values. For these parameters, the mean effective pressure is negatively correlated with the elevation-change anomaly over time (Fig. 5a). Due to the combined influence of the basal velocity and elevation change, positive and negative values of the effective pressure can exist within the lake boundary simultaneously. For example, a ring-shaped region of negative effective pressure forms near the boundary of the lake at the start of the filling stages (Fig. 5d). Similar behaviour with the opposite sign can occur during the draining stages.Figure 5.Synthetic example with nondimensional parameters $[eqn]$ and $[eqn]$ . (a) Time series of the mean elevation change, basal vertical velocity and effective pressure over the lake. (b)–(d) Map-plane plots of elevation change, basal vertical velocity and effective pressure at the time noted by the dashed vertical line in (a). The lake boundary is shown by the black circle.Table 1.Main parameters used in calculating the effective pressures of the Antarctic subglacial lakes (Figure 2). Data sources are provided in the ‘Data availability’ statement. The ‘Data’ section in the main text describes pre-processing of the elevation-change data and estimation of the ice-flow parameters (viscosity and basal drag). The off-lake pressure estimates $[eqn]$ are defined in Eqn. (30). The nondimensional parameters are defined by $[eqn]$ (Eqn. 17) and $[eqn]$ (Eqn. 21) where the observational timescale is $[eqn]$ yr. Parameter values are multiplied by the amount specified in the ‘units’ column.ParameterunitsMSLMac1David_s1_Byrd_s10_Physical $[eqn]$ : ice thicknessm100792419912678 $[eqn]$ : ice viscosityPa s ( $[eqn]$ )4.75.011.660.0 $[eqn]$ : basal drag coefficientPa s m $[eqn]$ ( $[eqn]$ )0.970.817.212.2 $[eqn]$ : off-lake pressure estimatekPa5.2700.003−0.9370.124Nondimensional $[eqn]$ : $[eqn]$ –0.3020.2620.2420.063 $[eqn]$ : $[eqn]$ –0.0100.0080.0610.027

The largest magnitudes of $[eqn]$ and $[eqn]$ occur near the centre of the lake where the deformation is largest, while the mean behaviour over the lake has a smaller magnitude that more closely corresponds to the behaviour near the lake boundary (Figs. 5 and 6a). In particular, the mean values of $[eqn]$ and $[eqn]$ over the lake show a weaker correlation than the pointwise behaviour near the centre of the lake (Fig. 6a).

We compare histograms of the effective pressure for all spatiotemporal points during the draining stage ( $[eqn]$ ) and the filling stage ( $[eqn]$ ) to highlight the bulk behaviour of the system (Fig. 6b-c). The effective pressure distribution has a peak around $[eqn]$ kPa during draining stages ( $[eqn]$ ), accompanied by a long tail of negative values (Fig. 6b). During filling stages ( $[eqn]$ ), the effective pressure distribution has a peak around $[eqn]$ kPa that is accompanied by a long tail of positive values (Fig. 6c). For these parameters ( $[eqn]$ and $[eqn]$ ), the effective pressure falls within $[eqn]$ kPa, which is a relatively small deviation from flotation ( $[eqn]$ ).Figure 6.(a) Basal vertical velocity $[eqn]$ versus the effective pressure $[eqn]$ in the synthetic example (Figure 5) for different values of $[eqn]$ . The nondimensional parameters are set to $[eqn]$ and $[eqn]$ . The colours of the points show the distance from the centre of the lake normalised by the distance to the boundary. The black ellipse corresponds to the spatial mean over the lake at each point in time. (b) Histogram of the effective pressure during the draining stages ( $[eqn]$ ) normalised by the total number of spatiotemporal points. (c) Histogram of the effective pressure during the filling stages ( $[eqn]$ ) normalised by the total number of spatiotemporal points.

Antarctic examples

Data

4.1.

We use Eqn. (18) to calculate the effective pressure $[eqn]$ from the elevation-change anomalies $[eqn]$ and the basal vertical velocity inversions $[eqn]$ over several Antarctic subglacial lakes (Fig. 2). The preprocessing of the elevation-change data and the inversion for the basal vertical velocity are described in detail by Stubblefield and others (2023a). We use the ICESat-2 ATL15 L3B Gridded Antarctic and Arctic Land Ice Height Change (Version 4) data product (Smith and others, 2024) to obtain elevation-change anomalies over the lakes by removing any regional thinning or thickening trend (Stubblefield and others, 2023a). We linearly interpolate the ATL15 data onto the fine spatiotemporal grid of the inverse method, which could produce errors in the pressure estimates if lake activity occurs more rapidly than the temporal resolution of ICESat-2 (91-day repeat cycle).

Removing any regional thickness-change trend is necessary for the inverse method and coincides with the assumption of flotation ( $[eqn]$ ) in the base state. To provide an independent point of comparison for the results, we use the approximation (B.1) that relates effective pressure to dynamic thickness changes (Appendix B) to estimate an ambient effective pressure via

$$[eqn]$$

where $[eqn]$ is the off-lake thickness change over the duration $[eqn]$ (5 years). We compute the off-lake thickness change by taking the spatial average of the ATL15 elevation-change product over all points that are at a distance greater than 80% from the centroid of the lake to the boundary of the computational domain (Stubblefield and others, 2021a). In the estimate (30), we have assumed a slip ratio of one and that all thickness changes arise from ice dynamics rather than snow accumulation (see Appendix B); removing any contributions from snow accumulation or assuming a smaller slip ratio would decrease the magnitude of the off-lake pressure estimates. We find that these estimates are at most only a few kilopascals in magnitude (Table 1), which is consistent with our results and underlying assumptions.

To invert the elevation-change data, we require the ice thickness $[eqn]$ , ice viscosity $[eqn]$ , basal drag coefficient $[eqn]$ and horizontal ice velocity $[eqn]$ that are associated with the background ice-flow state over the lakes (Fig. 1). We obtain horizontal surface velocity from the MEaSUREs Phase-Based Antarctic Ice Velocity Map (Version 1) (Mouginot and others, 2019a; 2019b) and ice thickness from MEaSUREs BedMachine Antarctica (Version 3) (Morlighem and others, 2020; Morlighem, 2022). The basal drag and ice viscosity are estimated with the Ice-sheet and Sea-level System Model (shelfy-stream approximation) by inverting the surface-velocity data (Morlighem and others, 2010; Larour and others, 2012). The ice flow-law coefficient is estimated by providing the empirical relation from Cuffey and Paterson (2010) (with a flow-law exponent of $[eqn]$ ) with an approximate depth-averaged temperature between the melting point and the surface temperature, which is obtained from the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) climate reanalysis (Global Modeling and Assimilation Office (GMAO), 2015; Gelaro and others, 2017). All ‘background’ values are obtained by averaging the data $[eqn]$ and modelled variables $[eqn]$ over a square region (60 km $[eqn]$ 60 km) surrounding the lake (Stubblefield and others, 2023a). Sources for all datasets and the code repository are provided in the ‘Data availability’ statement. Parameter values for each lake are reported in Table 1.

The Antarctic subglacial lakes that we explore here have been the subject of many previous studies: Mercer Subglacial Lake beneath the confluence of the Whillans Ice Stream and Mercer Ice Stream (Fricker and others, 2007; Fricker and Scambos, 2009; Siegfried and others, 2016); Mac1 beneath MacAyeal Ice Stream (Fricker and others, 2010; Siegfried and Fricker, 2021); Byrd_s10_ beneath Byrd Glacier (Smith and others, 2009; Wright and others, 2014); and David_s1_ (sometimes referred to as D2) beneath David Glacier in East Antarctica (Smith and others, 2009; Lindzey and others, 2020; Malczyk and others, 2023). These lakes cover a wide range of physical parameters arising from different ice flow regimes across East Antarctica and West Antarctica (Fig. 2). In particular, the lakes display a variety of combinations of the nondimensional parameters $[eqn]$ and $[eqn]$ (Table 1). We first discuss the effective pressure estimates for each lake and then compare the behaviours between the lakes. Lake boundaries derived from elevation-change anomalies for these lakes are provided by Siegfried and Fricker (2018); we compare the altimetry-derived boundaries to the spatial extent of the effective pressure estimates below. As noted in the description around Eqn. (20), the particular choice of the lake boundary does not substantially influence the calculated effective pressure but is of primary importance in deciding where the calculation is valid.

Mercer Subglacial Lake

4.2.

Mercer Subglacial Lake (MSL) exists beneath an ice thickness of $[eqn]$ km at the confluence of the Whillans Ice Stream and Mercer Ice Stream along the Siple Coast of West Antarctica (Fig. 2). MSL has filled and drained repeatedly since the beginning of the ICESat era (Fricker and others, 2007; Fricker and Scambos, 2009; Siegfried and others, 2016; Siegfried and Fricker, 2018; 2021). The estimated ice viscosity ( $[eqn]$ Pa s) is an order of magnitude smaller than the East Antarctic lakes considered herein. The basal drag coefficient ( $[eqn]$ Pa s m $[eqn]$ ) is slightly larger than Mac1, but the same order of magnitude or smaller than the East Antarctic lakes (Table 1).

We estimate the effective pressure from the elevation-change anomaly and basal vertical velocity (Fig. 7). Since the beginning of the ICESat-2 period, the elevation-change data show that MSL completed a multi-peaked drainage event and has more recently begun to refill. The mean effective pressure in the lake, $[eqn]$ , was small but positive (up to $[eqn]$ 0.5 kPa) during the draining stage and has become negative during the filling stage as the ice column is lifted upwards (Fig. 7c). Maps of the elevation-change anomaly, basal vertical velocity and effective pressure during the filling stage highlight the spatial variability in the system (Fig. 7b,d,f). We find that the effective pressure has both regions of positive and negative effective pressure, which results from the combined influence of the basal vertical velocity and elevation change (Fig. 7f). In particular, the centre of the lake has a positive effective pressure that is surrounded by a ring of negative effective pressure that forms as the lake refills. The same type of ringed structure is found in the synthetic example (Fig. 5). We consider the possible consequences of this behaviour in the discussion. The altimetry-derived lake boundary from Siegfried and Fricker (2018) closely corresponds to the spatial extent of the effective pressure anomaly over MSL.Figure 7.Elevation change ( $[eqn]$ , basal vertical velocity inversion ( $[eqn]$ ) and effective pressure ( $[eqn]$ ) for Mercer Subglacial Lake. (a) Time series of the mean value of the elevation change over the lake. (b) Map–plane contour plot of the elevation change at the time shown by the vertical dashed line in (a). The dashed black line shows the boundary selected for calculating the effective pressure while the solid grey line shows the boundary from Siegfried and Fricker (2018). (c) Time series of the mean basal vertical velocity and (d) map–plane plot at the time shown by the vertical dashed line. (e) Time series of the mean effective pressure (solid), effective pressure within 2 km of the boundary (long dashed) and the reference pressure $[eqn]$ (short dashed). (f) Map–plane plot of the effective pressure. The green hatched region corresponds to the values used to estimate the effective pressure near the lake boundary.

Mac1

4.3.

Subglacial lake Mac1 lies beneath the MacAyeal Ice Stream and has been observed to fill and drain repeatedly since the beginning of the ICESat era (Fricker and others, 2010; Siegfried and Fricker, 2021). The physical setting of Mac1 is similar to MSL with an ice thickness of $[eqn]$ 1 km and estimated viscosity on the order of $[eqn]$ Pa s. However, the basal drag coefficient is estimated to be slightly smaller (Table 1). These physical parameters lead to smaller values of the nondimensional parameters $[eqn]$ and $[eqn]$ relative to MSL. The beginning of ICESat-2 observations coincided with a quiescent period followed by a draining event with a prolonged period of lowstand that lasted approximately two years (Fig. 8). The mean effective pressure was $[eqn]$ 3 kPa in the lake after the draining event. More complex behaviour is seen in the spatial variability during the draining event where the effective pressure has positive and negative regions simultaneously (Fig. 8d). These complex patterns arise from thickness changes as the ice dynamically adjusts after the drainage event (Appendix B). The effective pressure showed a larger spatial extent than the altimetry-derived lake boundary during the draining event. Mac1 showed a second subsidence event with additional oscillations in the effective pressure during 2023–24.Figure 8.Elevation change ( $[eqn]$ , basal vertical velocity inversion ( $[eqn]$ ) and effective pressure ( $[eqn]$ ) for Mac1. (a) Time series of the mean value of the elevation change over the lake. (b) Map–plane contour plot of the elevation change at the time shown by the vertical dashed line in (a). The dashed black line shows the boundary selected for calculating the effective pressure while the solid grey line shows the boundary from Siegfried and Fricker (2018). (c) Time series of the mean basal vertical velocity and (d) map–plane plot at the time shown by the vertical dashed line. (e) Time series of the mean effective pressure (solid), effective pressure within 2 km of the boundary (long dashed) and the reference pressure $[eqn]$ (short dashed). (f) Map–plane plot of the effective pressure. The green hatched region corresponds to the values used to estimate the effective pressure near the lake boundary.

Davids1

4.4.

Subglacial lake David_s1_ lies beneath David Glacier, which feeds the Drygalski Ice Tongue in East Antarctica. David_s1_ showed activity during the ICESat era, along with several other lakes beneath David Glacier (Smith and others, 2009, referred to as lake D2 therein). Subsequent activity has been characterised by an overall upward trend in the elevation change over the lake (Lindzey and others, 2020; Malczyk and others, 2023). The physical setting in East Antarctica is characterised by a thicker ice cover of $[eqn]$ km and higher viscosity ( $[eqn]$ 10 $[eqn]$ Pa s) relative to the examples from West Antarctica (MSL and Mac1). The basal drag coefficient is an order of magnitude larger than the examples from West Antarctica, corresponding to slow ice flow speeds (Table 1). David_s1_ has continued to fill during the ICESat-2 period, causing an increasingly negative effective pressure within the lake (Fig. 9). The effective pressure and the basal vertical velocity show correlated oscillations over time at this higher viscosity value (Fig. 9c,e). These oscillations correspond to fluctuations in the rate of elevation change that manifest as slope changes in the elevation-change timeseries (Fig. 9a). The effective pressure shows a similar spatial extent to the elevation-change anomaly, although the spatial footprint appears to have changed significantly since the ICESat-era drainage event (Smith and others, 2009; Siegfried and Fricker, 2018; Lindzey and others, 2020). We revisit this discrepancy between the ICESat-era outline and the current boundary in the discussion.Figure 9.Elevation change ( $[eqn]$ , basal vertical velocity inversion ( $[eqn]$ ) and effective pressure ( $[eqn]$ ) for David_s1_. (a) Time series of the mean value of the elevation change over the lake. (b) Map–plane contour plot of the elevation change at the time shown by the vertical dashed line in (a). The dashed black line shows the boundary selected for calculating the effective pressure while the solid grey line shows the boundary from Siegfried and Fricker (2018). (c) Time series of the mean basal vertical velocity and (d) map–plane plot at the time shown by the vertical dashed line. (e) Time series of the mean effective pressure (solid), effective pressure within 2 km of the boundary (long dashed) and the reference pressure $[eqn]$ (short dashed). (f) Map–plane plot of the effective pressure. The green hatched region corresponds to the values used to estimate the effective pressure near the lake boundary.

Byrds10

4.5.

Subglacial lake Byrd_s10_ beneath Byrd Glacier in East Antarctica showed activity during the ICESat period (Smith and others, 2009). The physical setting in East Antarctica is characterised by a thicker ice cover of $[eqn]$ km and higher viscosity ( $[eqn]$ 10 $[eqn]$ Pa s) relative to the examples from West Antarctica (MSL and Mac1). However, the estimated basal drag coefficient is an order of magnitude smaller than for David_s1_ (Table 1). Byrd_s10_ was quiescent at the beginning of the ICESat-2 period and has drained over the course of several years (Fig. 10a). The draining is associated with an elevated mean effective pressure within the lake. However, there is a large region of negative effective pressure that forms near the centre of the lake (Fig. 10f). Negative effective pressures arise when tensile viscous stresses dominate Eqn. (18) during draining events ( $[eqn]$ ) with higher ice viscosity. The spatial extent of the effective pressure is much larger than the elevation-change anomaly due to these high tensile viscous stresses, which shows that subglacial pressure perturbations can extend beyond altimetry-derived lake outlines (Fig. 10d). In the discussion, we describe the consequences of a larger effective pressure footprint and how this phenomenon can be illustrated analytically via a Green’s function analysis (Appendix C).Figure 10.Elevation change ( $[eqn]$ , basal vertical velocity inversion ( $[eqn]$ ) and effective pressure ( $[eqn]$ ) for Byrd_s10_. (a) Time series of the mean value of the elevation change over the lake. (b) Map–plane contour plot of the elevation change at the time shown by the vertical dashed line in (a). The dashed black line shows the boundary selected for calculating the effective pressure while the solid grey line shows the boundary from Siegfried and Fricker (2018). (c) Time series of the mean basal vertical velocity and (d) map–plane plot at the time shown by the vertical dashed line. (e) Time series of the mean effective pressure (solid), effective pressure within 2 km of the boundary (long dashed) and the reference pressure $[eqn]$ (short dashed). (f) Map–plane plot of the effective pressure. The green hatched region corresponds to the values used to estimate the effective pressure near the lake boundary.

Effective pressure distributions

4.6.

The preceding examples show that the effective pressure can be more strongly correlated with the basal vertical velocity or the elevation-change anomaly, depending on the physical properties of the ice and bed surrounding the subglacial lake (Figs. 7-10). We examine the relations between $[eqn]$ and $[eqn]$ for each lake by plotting the values for all spatiotemporal points within the lake boundaries (Fig. 11). The West Antarctic lakes (MSL and Mac1) and David_s1_ show only weak correlations between the effective pressure and the basal vertical velocity while Byrd_s10_ shows a stronger correlation (Fig. 11). The stronger correlation is due to the higher estimated ice viscosity at Byrd_s10_, which also results in effective pressures that reach larger magnitudes. Histograms of the effective pressure for all spatiotemporal points within the lake boundaries show that the West Antarctic lakes and David_s1_ are more closely clustered around flotation ( $[eqn]$ ) while Byrd_s10_ has a wider and shorter distribution (Fig. 11). The effective pressures are generally on the order of $[eqn]$ 10 kPa, although they can locally reach values of $[eqn]$ kPa.Figure 11.Basal vertical velocity $[eqn]$ versus the effective pressure $[eqn]$ for the Antarctic subglacial lakes shown in Figures 7-10. Each point within the lake boundary is plotted for each point in time (blue points). Linear regressions are shown by the dashed black lines. Green histograms show effective pressure distributions normalised by the total number of points and the bin width.

Discussion

Here, we have developed a method for estimating effective pressures in subglacial lakes from altimetry data and applied the method to a collection of active subglacial lakes in Antarctica. We have shown that the effective pressure of a subglacial lake can be correlated with either the basal vertical velocity or the elevation-change anomaly. In particular, at higher values of the ice viscosity, we find that the effective pressure is more strongly correlated with the basal vertical velocity. While the subglacial lakes considered herein tend to have effective pressures within $[eqn]$ 10 kPa of flotation, locally elevated effective pressures can arise when the viscosity is large (e.g., Fig. 11). These results highlight the importance of accounting for viscous stresses when estimating the effective pressure and the potential difficulty of interpreting the basal stress state from elevation-change anomalies alone. Moreover, the computational efficiency of the method could facilitate quantifying uncertainties arising from the ice viscosity, which varies over orders of magnitude, as well as incorporating other constraints, like lake depth, in a statistical framework.

While we focused on how viscosity influences effective pressure estimates, the basal drag coefficient can also influence the results to a lesser degree. We showed that the effective pressure can be approximated by the logarithmic strain rate of the ice column multiplied by the ice viscosity and a slip ratio (Appendix B). In slow-flowing regions, the slip ratio could be small enough to drive the effective pressure towards flotation, but this may be counteracted by higher ice viscosity. The results herein are not strongly influenced by basal drag coefficient variability because all examples had a nondimensional sliding parameter of the same magnitude, resulting in approximately equal slip ratios (Table 1; Fig. B1a).

The subglacial lake effective pressures found here are small compared to other drainage elements such as channels and cavities where the effective pressure can reach megapascals in magnitude, which could have important consequences for the drainage system (e.g., Dow and others, 2022; Sommers and others, 2023). For example, negligible effective pressure near the boundary of a subglacial lake could inhibit closure of adjacent drainage elements, which in turn might facilitate shoreline migration via melting or sediment erosion (Evatt and others, 2006). The smallness of the subglacial lake effective pressure relative to neighbouring drainage elements could also generate large gradients in the hydraulic potential that influence the magnitude of water flow. Conversely, elevated effective pressures due to higher viscous stresses may result in a diminished hydraulic gradient, which could hinder water capture, depending on the pressure in the surrounding drainage system. The effective pressure of a subglacial lake being near flotation could therefore be important for influencing flooding events and water capture, although more comprehensive hydrology modelling would be needed to further constrain these dynamics.

The behaviour near the lake boundary is important in relation to potential interactions with the surrounding drainage system. We found a consistent spatial pattern where the effective pressure changes sign near the lake boundary, which arises from the combined influence of the elevation-change anomaly and basal velocity on dynamic thickness changes (Appendix B). For example, the effective pressure became negative near the boundary of Mercer Subglacial Lake when it began to refill, while positive effective pressures existed near the centre of the lake from the lowstand (Fig. 7). Negative effective pressures near the boundary of a lake could cause detachment of ice from the base and shoreline migration, or inhibit creep closure of a subglacial channel that drains the lake, setting the stage for future flooding events (Evatt and others, 2006; Fowler, 2009; Stubblefield and others, 2019, 2021b). On the other hand, we found positive effective pressures developing during draining events near the boundaries of Mac1 and Byrd_s10_, which could be associated with channel closure and the eventual cessation of flooding.

We also found that the areal extent of the effective pressure can extend beyond the spatial extent of the elevation-change anomaly due to viscous stresses within the ice (Fig. 10). To further analyse this behaviour, we derived a Green’s function for the simplified problem of steady creep and found that the effective pressure depends on the basal vertical velocity over zone that is roughly five ice thicknesses wide (Appendix C, Fig. B1). The presence of pressure perturbations that extend beyond the altimetry-based lake outline raises the issue of how subglacial lake boundaries should be identified. For example, it is difficult to determine if water that exists outside of the altimetry-based outline should be considered part of the drainage system or part of the lake without further modelling to assess the hydraulic state.

In a similar vein, the results show that the current lake boundaries do not necessarily coincide with the previous ICESat-era elevation changes, suggesting that subglacial shorelines have potentially migrated relative to previous positions (Siegfried and Fricker, 2021). David_s1_ provides an example of this behaviour where water is ponding in a different area relative to the ICESat-era observations (Fig. 9). A detailed aerogeophysical survey over the David Glacier lake system showed that a minimum in the hydraulic potential existed north of the ICESat-era outline (Lindzey and others, 2020), which is consistent with the recent ICESat-2 elevation changes (Fig. 9). While theory and observations suggest that shorelines can migrate during volume-change events, incorporating these physics into subglacial hydrology models remains an open challenge (Siegfried and Fricker, 2021; Stubblefield and others, 2021a, 2021b).

The effective pressure estimation method developed herein only applies to subglacial lakes, which we have defined as drainage elements that have a minimal hydraulic potential and small slope at the ice base. Other drainage elements like subglacial channels, cavities or sheets do not generally satisfy these conditions (Flowers, 2015). While these other drainage elements do not produce localised elevation-change anomalies, which cause variations in basal drag through sliding relations that indirectly influence elevation-change patterns by modulating ice flow speeds. Connecting elevation changes to subglacial water flow through a broader array of drainage elements is an open challenge that can be addressed in future work with more comprehensive hydrology models (Werder and others, 2013; Sommers and others, 2018).

As the dynamics of subglacial lakes are closely coupled with the surrounding drainage system, the effective pressure estimates derived herein can be used to further constrain or validate subglacial hydrology models in drainage basins hosting subglacial lakes. Likewise, testing subglacial hydrology models against these estimates would help to further constrain transient changes in subglacial sliding. Future work should therefore focus on constraining hydrology models and ice-flow models with these effective pressure estimates.

Conclusions

Here, we have developed a method for estimating the subglacial effective pressure from elevation-change data above active subglacial lakes. We applied the method to a collection of Antarctic subglacial lakes that have shown activity during the ICESat-2 era and found that the effective pressure in these lakes tended to remain within a few tens of kilopascals of flotation. Our analysis suggests that higher viscous stresses are associated with greater deviations from flotation. Similarly, we found that the effective pressure can be more strongly correlated with either the elevation-change anomaly or the basal vertical velocity over time, depending on the magnitude of the viscous stresses at the base. The effective pressure estimates developed herein can be used to help validate subglacial hydrology models with altimetry data and further constrain the dynamics of the glacier bed.