Reconstructing subglacial lake activity with an altimetry-based inverse method

Abstract

Subglacial lake water-volume changes produce ice-elevation anomalies that provide clues about water flow beneath glaciers and ice sheets. Significant challenges remain in the quantitative interpretation of these elevation-change anomalies because the surface expression of subglacial lake activity depends on basal conditions, rate of water-volume change, and ice rheology. To address these challenges, we introduce an inverse method that reconstructs subglacial lake activity from altimetry data while accounting for the effects of viscous ice flow. We use a linearized approximation of a Stokes ice-flow model under the assumption that subglacial lake activity only induces small perturbations relative to a reference ice-flow state. We validate this assumption by accurately reconstructing lake activity from synthetic data that are produced with a fully nonlinear model. We then apply the method to estimate the water-volume changes of several active subglacial lakes in Antarctica by inverting data from NASA's Ice, Cloud, and land Elevation Satellite 2 (ICESat-2) laser altimetry mission. The results show that there can be substantial discrepancies (20% or more) between the inversion and traditional estimation methods due to the effects of viscous ice flow. The inverse method will help refine estimates of subglacial water transport and further constrain the role of subglacial hydrology in ice-sheet evolution.

1. Introduction

Ice-sheet surface elevation responds to a variety of time-varying subglacial phenomena, including subglacial-lake volume change, basal-drag variations, and melting or freezing at the ice-water interface. Active subglacial lakes (i.e., those that experience observable volume change in the observational record) in particular have received much attention due to the localized perturbations they produce in ice-sheet surface elevation during volume-change events (e.g., Gray and others, 2005; Wingham and others, 2006; Fricker and others, 2007). NASA's Ice, Cloud, and land Elevation Satellite (ICESat) mission (2003-2009) and the European Space Agency's CryoSat-2 satellite altimetry mission (2010-present) facilitated the detection of over one hundred active subglacial lakes beneath the Antarctic Ice Sheet (e.g., Smith and others, 2009; Wright and Siegert, 2012; Fricker and others, 2016; Livingstone and others, 2022), driving investigations into their possible relation to fast ice flow (e.g., Stearns and others, 2008; Scambos and others, 2011; Siegfried and others, 2016) and into their ability to host microbial ecosystems (e.g., Christner and others, 2014; Achberger and others, 2016; Davis and others, 2023). Fewer subglacial lakes have been discovered beneath the Greenland Ice Sheet based on ice-surface changes, suggesting that there may be significant differences in subglacial hydrological conditions there relative to the Antarctic Ice Sheet (e.g., Bowling and others, 2019; Livingstone and others, 2019, 2022).

High-resolution satellite altimetry data from NASA's ICESat-2 mission (2018-present) presents a valuable opportunity to continue investigating dynamic conditions and constraining water budgets beneath ice sheets (e.g., Markus and others, 2017; Neckel and others, 2021; Siegfried and Fricker, 2021). ICESat-2's temporal resolution (91 day repeat cycle) and spatial resolution (~3.3 km across-track between beam pairs and ~90 m between beams within the pairs) provide unprecedented coverage of active subglacial lakes, which are typically greater than ~10 km in diameter and fill or drain over multiple years (Smith and others, 2009; Siegfried and Fricker, 2021; Stubblefield and others, 2021a). Modeling has shown that accurately estimating subglacial-lake volume change, areal extent, and highstand or lowstand timing from altimetry alone can be complicated by the effects of viscous ice flow (Stubblefield and others, 2021a). Basal vertical velocity anomalies associated with subglacial lake activity can manifest with a wider areal extent and smaller amplitude at the ice-sheet surface when ice flows laterally toward or away from the lake during volume-change events. Ice viscosity, ice thickness, and basal drag exert strong control on ice flow and, therefore, also influence the surface expression of subglacial lake activity (Stubblefield and others, 2021a). Although satellite altimetry data has been incorporated in basal-drag inversions (e.g., Larour and others, 2014; Arthern and others, 2015; Goldberg and others, 2015; Mosbeux and others, 2016), inverse methods that quantify subglacial-lake activity from altimetry and account for ice-flow effects have not yet been developed.

Inversion of time-varying altimetry data necessitates leveraging reduced-order models to alleviate the computational cost associated with repeatedly solving the forward problem. Dimensionality reduction is often achieved using ice-flow models that are based on depth-integrated approximations of the Stokes equations (e.g., Greve and Blatter, 2009). Solving the linearized Stokes equations on simplified domains with spectral methods is an alternative way to achieve computational efficiency when the full stresses in the ice must be resolved (e.g., Budd, 1970; Hutter and others, 1981; Balise and Raymond, 1985; Gudmundsson, 2003; Sergienko, 2012). Previous inversions relying on perturbation methods have not included time-varying data (Gudmundsson and Raymond, 2008; Thorsteinsson and others, 2003). Likewise, a computational method for inverting time-varying elevation data with perturbation-based models would be a valuable step toward quantifying time-varying subglacial lake perturbations. We use this small-perturbation approach as subglacial lake activity typically only induces small perturbations in ice-surface elevation (e.g., a few meters) relative to ice thickness.

Here, we derive, test, and apply an altimetry-based inverse method for quantifying basal vertical velocity perturbations that arise from subglacial lake activity. First, we outline the forward model for the perturbation in ice-surface elevation that is produced by a basal vertical velocity forcing. We then derive the inverse method from a least-squares optimization problem. To verify and validate the method, we present tests with synthetic data from both the linearized and nonlinear models. We then apply the method to a collection of active subglacial lakes in Antarctica (Fig. 1). The results show that ice flow can produce significant discrepancies between the inverse method and a traditional altimetry-based estimation method for calculating changes in subglacial water volume over the current ICESat-2 time period. We conclude by discussing limitations, extensions, and further applications of the method.

Figure 1. Map of ICESat-2 ATL15 gridded product (Smith and others, 2022) showing the elevation change of the Antarctic Ice Sheet between October 2018 and April 2022. The map-plane (x,  y) coordinates in the ATL15 dataset correspond to the Antarctic Polar Stereographic Projection (EPSG:3031). Insets show the locations of the subglacial lakes targeted as examples in this study. Subglacial lake boundaries derived from surface altimetry are shown as gray lines (Siegfried and Fricker, 2018). Regional thinning occurs around Thwaites Lake 170 (Thw₁₇₀) and regional thickening occurs around Mercer Subglacial Lake (SLM). Regional elevation-change trends around Slessor Glacier (lake Slessor₂₃), MacAyeal Ice Stream (lake Mac1), and Byrd Glacier (lake Byrd_s10) are less pronounced. We remove regional trends to produce elevation-change anomalies that are used in the inversions.

2. Method

In this section, we derive the forward model and the associated inverse method. First, we outline the general Stokes flow problem to highlight the governing equations and simplifying assumptions. Then, we outline a derivation of the small-perturbation model that is used in the inverse method. Finally, we derive the inverse method with a least-squares optimization approach.

2.1. Stokes flow

We assume that ice deforms as a viscous fluid according to the incompressible Stokes equations, which are given by

(1)$$-\nabla \cdot \pmb{\sigma} = \rho_{\rm i}\pmb{g},\; $$

(2)$$\nabla \cdot \pmb{u} = 0,\; $$

where ρ _i is the ice density, $\pmb {u} = [ u,\; \, v,\; \, w] ^{\rm T}$ is the ice velocity, and $\pmb {g} = g[ 0,\; \, 0,\; \, -1] ^{\rm T}$ denotes gravitational acceleration with magnitude g = 9.81 m s⁻² (Greve and Blatter, 2009; Stubblefield and others, 2021b). We have excluded possible elastic components of ice deformation by assuming a viscous rheology and revisit this choice in the discussion. The stress tensor $\pmb {\sigma }$ is defined via

(3)$$\pmb{\sigma} = -p{\sf I} + \eta\left(\nabla \pmb{u} + \nabla \pmb{u}^{\rm T} \right),\; $$

where p is the pressure, ${\sf I}$ is the identity tensor, and η is the viscosity. At the ice-bed boundary we assume a sliding law of the form

(4)$${\sf T}\pmb{\sigma}\pmb{n} = -\beta {\sf T}\pmb{u},\; $$

where β is the basal drag coefficient, $\pmb {n}$ is an outward-pointing unit normal to the boundary, and ${\sf T} = {\sf I} - \pmb {n}\pmb {n}^{\rm T}$ is a projection tangential to the ice-sheet surface (Stubblefield and others, 2021b). Although the small-perturbation model used in the inversions assumes a Newtonian viscosity and linear sliding law (i.e., constant η and β), we will also consider synthetic data produced by a fully nonlinear model with Glen's law viscosity (Glen, 1955) and a nonlinear Weertman-style sliding law (Weertman, 1957) to test the validity of these simplifications.

The upper surface of the ice-sheet z = h(x,  y,  t) evolves over time according to the kinematic equation

(5)$${\partial h\over \partial t} + u {\partial h\over \partial x} + v{\partial h\over \partial y} = w,\; $$

where the velocity components are evaluated at the surface (z = h). We assume that a stress-free condition,

(6)$$\pmb{\sigma}\pmb{n} = \pmb{0},\; $$

holds at the upper surface of the ice sheet. We approximate the spatial domain as a horizontally unbounded slab because the ice-sheet extent is much greater than areal extent of the subglacial lakes. Away from the lake, we assume that all quantities approach an appropriate far-field reference state that is based on data and available ice-sheet model output.

2.2. Small-perturbation model

Now we will describe the forward model that is used in the inverse method. Although small-perturbation models have been derived previously, we outline a derivation here to highlight the assumptions underlying the inverse method (Balise and Raymond, 1985; Gudmundsson, 2003). Our goal is to find the basal vertical velocity perturbation w _b that produces the surface elevation-change anomaly Δh _a under the assumption that these anomalies arise from subglacial lake activity (Fig. 2). We could also incorporate a basal drag anomaly to represent a slippery spot over a lake in the small-perturbation framework (e.g., Gudmundsson, 2003; Stubblefield, 2022), but the resulting dipolar elevation-change anomaly (Sergienko and others, 2007) is not discernible in any of the active lakes considered herein. We revisit this idea in the discussion.

Figure 2. Sketch of linearized model setup. The horizontal (map-plane) coordinates are (x,  y) with the y direction pointing into the page. The basal vertical velocity anomaly w _b produces an elevation-change anomaly Δh _a. The ice thickness is $\bar {H}$ and the horizontal surface velocity is $\bar {\pmb {u}}$ in the reference flow state. The ice flow is aligned with the x axis here for simplicity but generally also has a component in the y direction. The volume change estimated from the elevation-change anomaly Δh _a can deviate significantly from the subglacial water-volume change (Stubblefield and others, 2021a).

To derive a simplified model for this system, we assume that Δh _a and w _b are small perturbations from a known reference state that is (approximately) characterized by a constant ice thickness $\bar {H}$, horizontal surface velocity $\bar {\pmb {u}} = [ \bar {u},\; \, \bar {v}] ^T$, ice viscosity $\bar {\eta }$, and basal drag coefficient $\bar {\beta }$. We further assume that the basal surface is horizontal in the reference state and that the ice pressure equals the cryostatic pressure. Strictly speaking, an advective component is only present in the free-slip limit ($\bar {\beta } = 0$) under the assumption of a horizontally uniform Stokes flow over a flat bed subject to the stress boundary conditions (4) and (6). However, we retain a background advective velocity in all cases for consistency with the data.

Letting [u _h,  v _h,  w _h]^T denote the perturbation in ice-sheet surface velocity, we insert perturbations to the reference states, $h = \bar {H} + \Delta h_{\rm a}$ and $\pmb {u} = [ \bar {u},\; \, \bar {v},\; \, 0] ^T + [ u_h,\; \, v_h ,\; \, w_h] ^T$, into the kinematic equation (5) to obtain

(7)$${\partial \Delta h_{\rm a}\over \partial t} + \bar{u} {\partial \Delta h_{\rm a}\over \partial x} + \bar{v}{\partial \Delta h_{\rm a}\over \partial y} = w_h.$$

We have neglected terms involving products of perturbations in (7) under the assumption of small perturbations. We solve equation (7) by taking Fourier transforms with respect to the horizontal coordinates (x,  y) to obtain

(8)$${\partial \widehat{\Delta h_{\rm a}}\over \partial t} + ( i\pmb{k}\cdot\bar{\pmb{u}}) \widehat{\Delta h_{\rm a}} = \widehat{w_h},\; $$

where $\pmb {k} = [ k_x,\; \, k_y] ^{\rm T}$ is the horizontal wavevector. The vertical surface velocity is assumed to satisfy the Stokes flow problem (1)–(6), subject to the above simplifications, which allows us to derive a closed-form expression of the solution operator (Balise and Raymond, 1985; Gudmundsson, 2003; Stubblefield and others, 2021a).

We algebraically solve the Fourier-transformed Stokes problem to obtain an expression for the transformed vertical surface velocity,

(9)$$\widehat{w_h} = -{\cal R}\widehat{\Delta h_{\rm a}} + {\cal T}\widehat{w_{\rm b}} ,\; $$

in terms of the basal vertical velocity and surface elevation anomalies (e.g., Stubblefield and others, 2021a, Supporting Information). In Eqn. (9), ${\cal R}$ is a relaxation function that controls the decay rate of the elevation anomaly, and ${\cal T}$ is a transfer function that maps the basal vertical velocity anomaly to its surface expression. These functions depend on the scaled wavevector magnitude $k' = \vert \pmb {k}\vert \bar {H}$ and drag coefficient $\gamma = {\bar {\beta } \bar {H}} / ( {2\bar {\eta } k'} )$ through the relations

(10)$${\cal R} = \left({\rho_{\rm i} g\bar{H}\over 2\bar{\eta} k'}\right){ ( 1 + \gamma) e^{4{k'}} -( 2 + 4\gamma{k'} ) e^{2{k'}} + 1-\gamma \over ( 1 + \gamma) e^{4{k'}} + ( 2\gamma + 4{k'} + 4\gamma{k'}^2) e^{2{k'}} -1 + \gamma },\; $$

and

(11)$${\cal T} = {2( 1 + \gamma) ( {k'} + 1) e^{3{k'}} + 2( 1-\gamma) ( {k'}-1) e^{{k'}} \over ( 1 + \gamma) e^{4{k'}} + ( 2\gamma + 4{k'} + 4\gamma{k'}^2) e^{2{k'}} -1 + \gamma }.$$

For a detailed derivation of the expressions (10) and (11) see, for example, Stubblefield and others (2021a, Supporting Information) and Stubblefield (2022, Appendix E).

Substituting the expression (9) into (8), we find that the ice-surface elevation anomaly Δh _a evolves in frequency space via

(12)$${\partial \widehat{\Delta h_{\rm a}}\over \partial t} + ( i\pmb{k}\cdot\bar{\pmb{u}}) \widehat{\Delta h_{\rm a}} = -{\cal R}\widehat{\Delta h_{\rm a}} + {\cal T}\widehat{w_{\rm b}}.$$

The solution to equation (12) is given by

(13)$$\widehat{\Delta h_{\rm a}} = \widehat{\Delta h_0}e^{-( i\pmb{k}\cdot\bar{\pmb{u}} + {\cal R}) t} + \widehat{w_{\rm b}} \ast {\cal K},\; $$

where * denotes convolution over time and Δh ₀ is the elevation perturbation at the initial time t = 0. The kernel ${\cal K}$, defined by

(14)$${\cal K} = {\cal T}e^{-( i\pmb{k}\cdot\bar{\pmb{u}} + {\cal R}) t},\; $$

controls the decay of the elevation-change anomaly and transfer of the basal anomaly to the surface. The characteristic time scale for the decay of surface-elevation anomalies is

(15)$$t_{\rm relax} = {2\bar{\eta}\over \rho_{\rm i}g\bar{H}},\; $$

which controls the magnitude of the relaxation function ${\cal R}$ (cf. Turcotte and Schubert, 2014, Chapter 6). The effects of viscous ice flow influence the surface expression of lake activity when the viscous relaxation time t _relax is comparable to the lake filling or draining timescale (Stubblefield and others, 2021a). We highlight the importance of the viscous relaxation time in the examples below.

2.3. Inverse method

Now we will outline the inverse method. We let ${\sf F}$ denote the (map-plane) Fourier transform operator and define the relative elevation-change anomaly via

(16)$$d = {\Delta h_{\rm a}} -{\sf F}^{-1}\left(e^{-( i\pmb{k}\cdot\bar{\pmb{u}} + {\cal R}) t}{\sf F}( {\Delta h}_0) \right),\; $$

which has the contribution from the initial value in equation (13) removed. From equation (13), we define the operator ${\sf G}$ that maps w _b to the relative elevation change d via

(17)$${\sf G}( w_{\rm b}) = {\sf F}^{-1}( {\sf F}( {w}_b) \ast {\cal K}) ,\; $$

where the kernel ${\cal K}$ is defined in equation (14).

We consider a regularized least-squares objective functional,

(18)$${J}( w_{\rm b}) = {1\over 2} \int_{0}^{t_{\rm f}} \int_{-\infty}^{ + \infty}\int_{-\infty}^{ + \infty} \vert {\sf G}( w_{\rm b}) - d \vert ^2 \; {\rm d}x\, {\rm d}y \, {\rm d}t + {\varepsilon\over 2} \int_{0}^{t_{\rm f}} \int_{-\infty}^{ + \infty}\int_{-\infty}^{ + \infty} \vert \nabla w_{\rm b} \vert ^2 \; {\rm d}x\, {\rm d}y \, {\rm d}t,\; $$

where t _f is the final time and the parameter $\varepsilon$ controls the strength of the regularization term. While the regularization in (18) promotes smoothness, other regularizations could be chosen to promote sparsity of the basal forcing, for example (Stadler, 2009). The minimizer of the objective (18) satisfies the normal equation

(19)$${\sf G}^\ast ( {\sf G}( w_{\rm b}) ) -\varepsilon \nabla^2 w_{\rm b} = {\sf G}^\ast ( d) ,\; $$

which can be derived with variational calculus (Vogel, 2002; Hanke, 2017). The adjoint operator ${\sf G}^\ast$ in (19) is given by

(20)$${\sf G}^\ast ( f) = {\sf F}^{-1}( {\sf F}( f) \star {\cal K}) $$

for any function f, where $\star$ denotes cross-correlation over time.

We solve the equation (19) with the conjugate gradient method to obtain the basal vertical velocity w _b. In using the conjugate gradient method to solve this operator equation, we avoid explicitly constructing matrices corresponding to the forward and adjoint operators, and instead simply require the action of these operators on functions (Atkinson and Han, 2009, Section 5.6). We implemented the discretized inverse method in Python with SciPy's fast Fourier transform and convolution algorithms (Virtanen and others, 2014). The code is openly available (https://doi.org/10.5281/zenodo.8371416).

2.4. Estimation of water-volume change

To compare the inversion with previous estimation methods, we will focus on estimating subglacial water-volume changes. Given the basal vertical velocity inversion w _b, the basal water-volume change over a map-plane area B can be computed via

(21)$$\Delta V_{\rm inv} ( t) = \int_{0}^t \left[\iint_{B} w_{\rm b} \; {\rm d}{x}\, {\rm d}{y}\right]{\rm d}t' .$$

Alternatively, the volume-change has often been estimated in previous studies by integrating the elevation change anomaly over the static outline of a lake (Fricker and Scambos, 2009; Smith and others, 2009). Using this approach, the water-volume change is estimated by

(22)$$\Delta V_{\rm alt} ( t) = \iint_{B} \Delta h_{\rm a} - \Delta h_0 \; {\rm d}{x}\, {\rm d}{y},\; $$

where we have integrated over the same map-plane area B. Although an alternative lake boundary could be identified with the inversion, we use the same boundary to calculate both estimates for consistency in comparison. We revisit this problem in the discussion.

In the limits ${\cal R}\to 0$ and ${\cal T}\to 1$, equation (12) implies that these volume changes are equivalent (i.e., ΔV _inv = ΔV _alt). This “rigid-ice” limit is approached when the ice is viscous enough for the relaxation timescale, t _relax (eq. (15)), to greatly exceed the volume-change timescale (Stubblefield and others, 2021a). Although incompressibility causes these volume changes to be equal when integrating over the entire areal extent of a glacier, this approach is impractical for the Antarctic Ice Sheet due to the presence of multiple lakes and regional thickening or thinning trends. We explore the discrepancy between the inversion-derived estimate (21) and surface-derived estimate (22) for a range of parameters in the examples below.

3. Synthetic examples

Before applying the method to the ICESat-2 altimetry data, we first solve two problems with synthetic data to validate the method and illustrate the range of behaviors found in the ICESat-2 examples below. First, we verify the implementation by inverting synthetic data that is produced by prescribing the linearized model with a known basal vertical velocity field and then adding Gaussian white noise to the resulting elevation change. For consistency with the ICESat-2 examples, we remove a small off-lake elevation-change component, Δh _off, from the elevation change as detailed in the next section. For this example, we choose a basal vertical velocity field that is a Gaussian bump undergoing sinusoidal oscillations in time. The inverse method is able to reconstruct the basal vertical velocity and volume-change time series from the synthetic data (Fig. 3). We find that there is little deviation ($\lesssim 5\%$) between the volume-change estimates (21) and (22) on short timescales (i.e., less than ~2.5 years), whereas large deviations occur over decadal timescales. This behavior arises because the viscosity is $\bar {\eta } = 10^{15}$ Pa s for this example, leading to characteristic relaxation timescale of t _relax ≈ 2.8 yr. These results highlight that there will not be significant deviations between the altimetry-based and inversion-based volume-change estimates over the current ICESat-2 time period if the ice viscosity reaches this magnitude. We provide an example of this behavior below. In all examples herein, we set the regularization parameter to $\varepsilon = 1$ in equation (19), which results in accurate reconstructions of the synthetic examples without over-fitting the data.

Figure 3. Inversion results for synthetic data produced with the linearized model. (a) Map-plane elevation anomaly and inversion at t = 7 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series. The ice flow direction is shown by the black arrow in (a). The maximum deviation between the surface-derived volume change and the inversion in (b) is 0.83 km³, or 48% of the maximum amplitude of the surface-derived estimate. The inversion accurately recovers the true water-volume change (dashed black line). The parameters for this example are $\bar {H} = 2500$ m, $\bar {\eta } = 10^{15}$ Pa s, $\bar {\beta } = 10^{11}$ Pa s m⁻¹, $\bar {u} = 200$ m yr⁻¹, and $\bar {v} = 0$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 2.82  yr. The pink line marks the time step shown in (a). See Movie S1 for an animation of the inversion over all time steps.

Next, we show an example with synthetic data produced by a fully nonlinear model to test the assumptions underlying the small-perturbation approach (Stubblefield and others, 2021b, 2021a). The nonlinear model assumes a Glen's law viscosity (Glen, 1955; Cuffey and Paterson, 2010), a nonlinear Weertman-style sliding law (Weertman, 1957), fully nonlinear surface kinematic equations, and vanishing basal drag over the lake. For this example, we have assumed radial symmetry with respect to the map-plane coordinates (x,  y) to facilitate numerical solution in three spatial dimensions. We also prescribe a more complex volume-change time series with a duration of three years in accordance with the current ICESat-2 time span and choose a lower viscosity for this example, $\bar {\eta } = 10^{14}$ Pa s (Fig. 4). Despite the simplifications inherent to the inverse method, the inversion accurately recovers the volume change time series that is produced by the nonlinear model (Fig. 4). Most importantly, the inversion is much more accurate than the surface-based volume change for this parameter regime. This example shows that ice viscosities on the order of $\bar {\eta } = 10^{14}$  Pa s can result in significant volume-change deviations over the current ICESat-2 time span. In particular, the examples in Figure 3 and Figure 4 show that the altimetry-based estimate tends to underestimate the magnitude of the true water volume change, regardless of whether the volume change is positive or negative. Next, we describe the data and preprocessing steps before discussing examples of ICESat-2 inversions.

Figure 4. Inversion results for synthetic data produced with a radially-symmetric nonlinear Stokes model (Stubblefield and others, 2021b). (a) Map-plane elevation anomaly and inversion at t = 1.7 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series. The maximum deviation between the surface-derived volume change and inversion in (b) is 0.15 km³, or 56% of the maximum amplitude of the surface-derived estimate. The inversion accurately recovers the true water-volume change (dashed black line). The parameters for this example are $\bar {H} = 1500$ m, $\bar {\eta } = 10^{14}$ Pa s, $\bar {\beta } = 10^{10}$ Pa s m⁻¹, $\bar {u} = 0$ m yr⁻¹, and $\bar {v} = 0$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 0.47 yr. The pink line marks the time step shown in (a). See Movie S2 for a detailed animation of the nonlinear model and Movie S3 for an animation of the inversion over all time steps.

4. Data and preprocessing

We use the ICESat-2 ATL15 L3B Gridded Antarctic and Arctic Land Ice Height Change (Version 2) data product (Smith and others, 2022) to obtain elevation-change anomalies above the Antarctic subglacial lakes shown in Figure 1. The map-plane (x,  y) coordinates in the ATL15 dataset correspond to the Antarctic Polar Stereographic Projection based on the WGS-84 ellipsoid (EPSG:3031). For the examples explored here, we interpolated the ICESat-2 ATL15 data onto a space-time grid with 100 points in each direction (t,  x,  y) to obtain the same resolution as the numerical model. Alternatively, we could restrict the model-data misfit in (18) to the discrete set of data points, but this could require additional temporal regularization that we have not included in this study. We remove any regional thickening or thinning trends by subtracting the spatially averaged off-lake component, Δh _off, as described below. We also have to establish a reference elevation profile to define the elevation-change anomaly. By default, the elevation changes in ATL15 are relative to the ice-surface elevation on January 1, 2020. In general, the elevation anomaly can be defined relative to any of the ATL15 time points by subtracting the elevation surface at a particular reference time t _ref. Therefore, the elevation change anomaly is derived from the ATL15 elevation change product Δh via

(23)$$\Delta h_{\rm a}( x,\; y,\; t) = \Delta h( x,\; y,\; t) - \Delta h_{\rm off}( t) - \left[\Delta h( x,\; y,\; t_{\rm ref}) - \Delta h_{\rm off}( t_{\rm ref}) \right]$$

where Δh _off is the (time-varying) spatial average of Δh away from the lake. Here, the spatial average is taken over all points that are at a distance greater than $80\%$ from the centroid of the lake to the boundary of the computational domain.

Based on previously identified lake activity, an appropriate reference time t _ref to define the anomalies happens to be the initial time in the ATL15 product, October 1, 2018, for all of the lakes considered here except Mercer Subglacial Lake (SLM). SLM reached an apparent highstand near the end of 2017 before beginning a drainage event during the ICESat-2 period (Siegfried and Fricker, 2021), so the initial time in the ICESat-2 data does not correspond to an elevation anomaly of zero. We elaborate on this decision for each lake in more detail below and provide further commentary on preprocessing considerations in the discussion.

To invert the elevation-change data, we also must supply the approximate ice thickness $\bar {H}$, ice viscosity $\bar {\eta }$, basal drag coefficient $\bar {\beta }$, and horizontal ice velocity $\bar {\pmb {u}} = [ \bar {u},\; \, \bar {v}] ^{\rm T}$ that describe the reference ice-flow state (Fig. 2). The viscosity and basal drag estimates are derived from the inversions presented in Arthern and others (2015), which relied on the ALBMAP ice thickness (Le Brocq and others, 2010) and the MEaSUREs InSAR-Based Antarctic Ice Velocity Map (Version 1) (Rignot and others, 2011; Mouginot and others, 2012). However, 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) for greater compatibility with the ICESat-2 epoch. All parameter values are obtained by calculating the mean of these data over the extent of the computational domain. The parameter values for each example are reported in Table 1 and the figure captions. To define the boundaries B in the volume estimation equations (21) and (22), we use the latest subglacial boundary inventory (Siegfried and Fricker, 2018), which is a compilation of static active subglacial lake outlines from a variety of sources that used mixed delineation methods.

Table 1. Parameters used in the inversions of the Antarctic subglacial lakes shown in Figure 1

Data sources are described in the “Data and Preprocessing” section. The velocity and centroid are listed in terms of the Antarctic Polar Stereographic Projection coordinates (EPSG:3031).

5. ICESat-2 examples

Next, we will invert ICESat-2 data (ATL15 gridded elevation-change product) for the subglacial lakes shown in Figure 1: Lake Mac1 beneath the MacAyeal Ice Stream (e.g., Fricker and others, 2010; Siegfried and Fricker, 2018, 2021), Mercer Subglacial Lake at the confluence of Mercer Ice Stream and Whillans Ice Stream (e.g., Fricker and others, 2007; Siegfried and Fricker, 2018, 2021; Siegfried and others, 2023), Slessor₂₃ beneath Slessor Glacier (Siegfried and Fricker, 2018; Siegfried and others, 2021), Thw₁₇₀ beneath Thwaites Glacier (Smith and others, 2017; Hoffman and others, 2020) and Byrd_s10 beneath Byrd Glacier in East Antarctica (Smith and others, 2009; Wright and others, 2014). These lakes have been the subject of numerous previous investigations and represent a wide range of filling-draining patterns, physical conditions, and locations across the Antarctic Ice Sheet (Table 1). For these examples, it is important to consider the reference time t _ref used to define the elevation anomaly in equation (23). We base our choices on the lake activity leading up to the ICESat-2 epoch. For example, subglacial lake Mac1 showed little activity since the beginning of the ICESat-2 epoch in 2018 (Siegfried and Fricker, 2021), suggesting that the initial time in the ATL15 data is an appropriate choice of reference time. For Mac1, there is a maximum discrepancy of ~0.12 km³ between the surface-based and inversion-based volume-change estimates, or 24% of the maximum amplitude of the surface-derived estimate (Fig. 5).

Figure 5. Inversion results for subglacial lake Mac1. (a) Map-plane elevation anomaly and inversion at t = 1.5 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series (Siegfried and Fricker, 2018). The ice flow direction is shown by the black arrow in (a). The maximum deviation between the surface-derived volume change and inversion is 0.09 km³, or 24% of the maximum amplitude of the surface-derived estimate. The parameters for this example are $\bar {H} = 926$  m, $\bar {\eta } = 2.3\times 10^{14}$ Pa s, $\bar {\beta } = 7.4\times 10^{10}$ Pa s m⁻¹, $\bar {u} = 334$ m  yr⁻¹, and $\bar {v} = -178$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 1.73 yr. The pink line marks the time step shown in (a). See Movie S4 for an animation of the inversion over all time steps.

We also show inversions of Mercer Subglacial Lake (SLM), which displays multiple oscillations over the ICESat-2 period (Fig. 6). We set the reference time to be t = 1.3 yr after the initial time (i.e. around the second peak in the time series), as this more closely corresponds to the long-term mean of Mercer Subglacial Lake's oscillation pattern (Siegfried and Fricker, 2021). For this example, we find a maximum discrepancy of ~0.05  km³ between the surface-based and inversion-based volume-change estimates, or 19% of the maximum amplitude of the surface-derived estimate.

Figure 6. Inversion results for Mercer Subglacial Lake (SLM in Fig. 1). (a) Map-plane elevation anomaly and inversion at t = 2.5 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series (Siegfried and Fricker, 2018). The ice flow direction is shown by the black arrow in (a). The maximum deviation between the surface-derived volume change and inversion in (b) is 0.05 km³, or 19% of the maximum amplitude of the surface-derived estimate. The parameters for this example are $\bar {H} = 1003$ m, $\bar {\eta } = 2.2\times 10^{14}$ Pa s, $\bar {\beta } = 3.7\times 10^{11}$ Pa s m⁻¹, $\bar {u} = 172$ m yr⁻¹, and $\bar {v} = -65$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 1.56 yr. The pink line marks the time step shown in (a). See Movie S5 for an animation of the inversion over all time steps.

We also invert elevation anomalies from Slessor Glacier (lake Slessor₂₃) and Thwaites Glacier (lake Thw₁₇₀). Slessor₂₃ shows a discrepancy of ~0.52 km³ between the volume-change estimates, which is 62% of the maximum amplitude of the surface-derived estimate (Fig. 7). Thw₁₇₀ also shows a large discrepancy of ~0.21 km³, or $49\%$ of the maximum in the altimetry-based estimate (Fig. 8). For Slessor₂₃, the initial time in the ICESat-2 data appears to be close to the midpoint of a filling stage, so this reference time seems appropriate for defining the elevation anomaly (Siegfried and Fricker, 2018). On the other hand, Thw₁₇₀ appears to be coming out of a quiescent post-drainage period at the beginning of the ICESat-2 period, so choosing the correct reference time is more ambiguous in this case (Hoffman and others, 2020; Malczyk and others, 2020). For example, setting the reference time to t = 1.5 yr instead results in a maximum discrepancy of ~0.075 km³ between the volume-change estimates for the Thw₁₇₀ inversion. This discrepancy arises because the magnitude of the elevation-change anomaly is diminished when choosing the different reference time and less of the signal is attributed to the basal forcing. We quantify the sensitivity to the reference time more thoroughly in Appendix A and highlight the main issues in the discussion.

Figure 7. Inversion results for subglacial lake Slessor₂₃. (a) Map-plane elevation anomaly and inversion at t = 2.7 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series (Siegfried and Fricker, 2018). The ice flow direction is shown by the black arrow in (a). The maximum deviation between the altimetry-derived volume change and inversion in (b) is 0.52 km³, or 62% of the maximum amplitude of the surface-derived estimate. The parameters for this example are $\bar {H} = 1735$ m, $\bar {\eta } = 2.4\times 10^{14}$ Pa s, $\bar {\beta } = 2.7\times 10^{10}$ Pa s m⁻¹, $\bar {u} = -141$ m yr⁻¹, and $\bar {v} = -146$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 0.97 yr. The pink line marks the time step shown in (a). See Movie S6 for an animation of the inversion over all time steps.

Figure 8. Inversion results for subglacial lake Thw₁₇₀. (a) Map-plane elevation anomaly and inversion at t = 2.8 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series (Smith and others, 2017). The ice flow direction is shown by the black arrow in (a). The maximum deviation between the altimetry-derived volume change and inversion is 0.21  km³, or 49% of the maximum amplitude of the surface-derived estimate. The parameters for this example are $\bar {H} = 2558$ m, $\bar {\eta } = 5.7\times 10^{14}$ Pa s, $\bar {\beta } = 1.3\times 10^{10}$ Pa s m⁻¹, $\bar {u} = -130$ m yr⁻¹, and $\bar {v} = -78$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 1.58 yr. The pink line marks the time step shown in (a). See Movie S7 for an animation of the inversion over all time steps.

The common theme of the preceding examples is that they have ice viscosities on the order of $\bar {\eta } = 10^{14}$ Pa s (Table 1) and volume-change discrepancies that are at least $\sim 20\%$ of the maximum of the altimetry-based estimate (Figs. 5–8). The range of basal drag coefficients and ice thicknesses across these examples (Table 1) suggests that the ice viscosity is the primary parameter controlling the volume-change discrepancies. At higher viscosity values, the volume-change discrepancies diminish over the current ICESat-2 time period because the viscous relaxation time exceeds the oscillation timescale. To illustrate this behavior, we inverted ICESat-2 data over subglacial lake Byrd_s10 and found a much smaller discrepancy ($\sim 4\%$) between the surface-based and inversion-based volume estimates (Fig. 9). This lack of discrepancy arises because the ice over this lake has a viscosity of $\bar {\eta } = 5\times 10^{15}$ Pa s, an order of magnitude higher than the preceding ICESat-2 examples. In this case, the surface and basal motion correspond more closely because the viscous relaxation time, t _relax = 13 yr, is much longer than the current ICESat-2 time span. However, over decadal timescales larger discrepancies are still possible for this parameter regime (e.g., Fig. 3) unless the lake oscillation period is small compared to the relaxation time (Stubblefield and others, 2021a).

Figure 9. Inversion results for subglacial lake Byrd_s10. (a) Map-plane elevation anomaly and inversion at t = 2.5 yr. (b) Time series of the surface-derived volume change (ΔV _alt), the inversion-based volume change (ΔV _inv), and the off-lake component (Δh _off) that is removed prior to inversion. The gray contours in (a) and (b) show the boundaries used to compute the volume-change time series. The ice flow direction is shown by the black arrow in (a). The maximum deviation between the altimetry-derived volume change and inversion is 9 × 10⁻³ km³, or 4% of the maximum amplitude of the surface-derived estimate. The parameters for this example are $\bar {H} = 2676$ m, $\bar {\eta } = 5\times 10^{15}$ Pa s, $\bar {\beta } = 1.4\times 10^{11}$ Pa s m⁻¹, $\bar {u} = -9.4$ m yr⁻¹, and $\bar {v} = -9.8$ m yr⁻¹. The viscous relaxation time associated with these parameters is t _relax = 13 yr. The pink line marks the time step shown in (a). See Movie S8 for an animation of the inversion over all time steps.

6. Discussion

Several practical and technical challenges are worth considering when applying the inverse method. From a practical viewpoint, the primary challenge is deriving the elevation anomaly from the altimetry data. For example, the inversion results may be sensitive to the details of how any regional thickening or thinning trends are separated from the lake-related elevation changes (Fricker and Scambos, 2009; Smith and others, 2009; Siegfried and Fricker, 2018, 2021). The reference elevation profile that is used to define the elevation anomaly from the data can also influence the inversion results, as we discussed in the case of subglacial lake Thw₁₇₀. Likewise, choosing an appropriate reference elevation profile may be difficult when the ice-sheet surface profile is heavily textured or the initial time in the data is during a volume-change event. In the latter case, we have relied on records of lake oscillations from previous satellite altimetry missions to choose appropriate reference times (Siegfried and Fricker, 2018, 2021). In Appendix A, we quantify the sensitivity of inversion results to the choice of reference time for the synthetic data (Fig. 10) and Thw₁₇₀ (Fig. 11). The results highlight the importance of carefully considering the reference time or elevation profile that is used to define the elevation-change anomaly (Appendix A). We leave further exploration of the sensitivity of inversion results to preprocessing steps for future work.

Figure 10. Inversion of synthetic data from Figure 3 after redefining the reference time t _ref in equation (23) to a range of incorrect values. The correct reference time in this example is t _ref = 0. Significant deviations between the inversion and true solution can occur if an incorrect reference time is chosen.

Figure 11. Inversion of the Thw₁₇₀ data from Figure 8 after redefining the reference time t _ref in equation (23) to a range of alternative values.

The primary technical limitations of the perturbation-based inverse method is that the associated forward models are inherently linear, posed on geometrically simple domains, and cannot deviate significantly from the specified reference state. Although we have tested the validity of the method by inverting synthetic data from a simple radially symmetric nonlinear problem (Fig. 4), more complex problems could require alternative methods. For example, we would caution against applying the method to complex stress regimes like an ice-stream shear margin, which is the case for Engelhardt Subglacial Lake that lies beneath the margin of the Whillans Ice Stream (Fricker and others, 2007; Siegfried and others, 2016). We have also neglected any elastic components of ice deformation by assuming that ice flows as a viscous fluid because the filling-draining events considered herein span multiple years. However, elastic deformation can arise on shorter timescales near grounding lines when lake activity is related to tidal cycles (Milillo and others, 2017). Subglacial lakes have also been observed to drain on sub-weekly timescales during jökulhlaups in Iceland (Björnsson, 2002; Evatt and Fowler, 2007). Moreover, we have assumed that all of the elevation anomaly is driven by ice deformation rather than surface mass balance anomalies, which could arise, for example, from snow infilling the lake depression (Malczyk and others, 2020).

We have also assumed that, to first order, the subglacial lakes do not coincide with reductions in basal drag because the characteristic dipolar elevation anomaly associated with such slippery spots is not discernible in the examples considered herein (e.g., Gudmundsson, 2003; Sergienko and others, 2007). However, some large, inactive Antarctic subglacial lakes are known to coincide with slippery spots where the ice surface is flat over most of the lake except on the upstream side where thinning occurs and the downstream side where thickening occurs (Bell and others, 2006, 2007; Wright and Siegert, 2012). On the other hand, several West Antarctic ice streams also have both subglacial lakes and localized regions of anomalously high basal drag (sticky spots) in close proximity (Winberry and others, 2009; Sergienko and Hulbe, 2011; Winberry and others, 2014; Siegfried and others, 2016). This coupling can arise when enhanced basal traction provides a meltwater source for the lakes or when lake drainage causes channelization that removes water from surrounding regions of the bed (Sergienko and Hulbe, 2011). Joint inversion for basal drag variations and lake activity would help to further refine our understanding of these coupled systems and may be tractable if an additional data source like surface velocity is available.

In this study, we have focused primarily on estimating subglacial water-volume changes. Another application of the inverse method will be estimating subglacial lake shorelines or areal extent. Lake boundaries are currently defined using ice-surface deformation extent to generate static lake boundaries (Siegfried and Fricker, 2018); however, these static boundaries were typically generated using lower spatial resolution altimetry instruments than are available today. This static view of lake boundaries has resulted in a number of lake re-delineation attempts (e.g., Fricker and others, 2014; Siegfried and Fricker, 2018) and more recent suggestions of time-variable lake boundaries (e.g., Neckel and others, 2021; Siegfried and Fricker, 2021). In our study, it is clear that static subglacial lake boundaries do not dependably encompass the ICESat-2 surface height change observations (Figs. 5–9) likely because lake shorelines vary temporally. Additionally, recent numerical modeling shows the surface-derived boundaries can have a larger areal extent than the true lake boundary at the base (Stubblefield and others, 2021a). With our inverse method, we could attempt to reconstruct subglacial shoreline evolution by tracking the areal extent of the basal forcing rather than the surface deformation. Improving the accuracy of subglacial-lake shoreline estimates in this way could be valuable for site selection in future subglacial drilling projects (Tulaczyk and others, 2014; Priscu and others, 2021) and thereby provide stronger constraints on subglacial microbial and geochemical processes (Christner and others, 2014; Achberger and others, 2016; Davis and others, 2023).

The inverse method could be extended to estimate other subglacial hydrological quantities besides water-volume changes. For example, the temporal derivative of the volume change can be related to the relative volumetric water discharge into (or out of) the lake (Evatt and Fowler, 2007). The water discharge naturally appears in models of glacial lakes that are coupled to subglacial channel evolution (Fowler, 1999, 2009; Kingslake, 2015; Carter and others, 2017; Stubblefield and others, 2019; Jenson and others, 2022). Finally, an alternative to prescribing a basal vertical velocity anomaly would have been prescribing a basal pressure anomaly. In particular, converting between vertical velocity and pressure perturbations is straightforward with the spectral method employed here (e.g., Gudmundsson, 2003, equation 54). Pressure perturbations could possibly be related to the subglacial effective pressure, the difference between the cryostatic pressure and water pressure in the lake, since we have assumed that the pressure in the reference state is cryostatic (cf. Evatt and Fowler, 2007; Stubblefield and others, 2021b). Estimating the effective pressure or other hydrological quantities with altimetry would be valuable for further constraining the physics of subglacial drainage systems.

7. Conclusions

We have introduced and applied an inverse method for estimating the basal forcing associated with subglacial lake activity from ice-sheet altimetry. We have provided some validation of the small-perturbation approach by inverting synthetic data from a nonlinear subglacial lake model to obtain a basal vertical velocity field and water-volume change time series that agree with the nonlinear model. We then applied the method to a collection of Antarctic subglacial lakes by inverting satellite altimetry data from NASA's ICESat-2 mission. These results illustrate that there can be significant discrepancies between surface-based estimation methods and the inversion due to the effects of viscous ice flow. In particular, the results show that surface-based estimation methods can underestimate changes in subglacial water volume. The inverse method provides a simple way to refine basal water budget contributions derived from active subglacial lakes and further illuminate the physics of subglacial hydrological systems with satellite altimetry.