Governing equations

Ice flow is governed by the Stokes equations, a reduced form of the momentum equations that neglects inertial terms, yielding a balance between the stress divergence and body forces

(1)$${\partial \sigma_{ij}\over \partial x_j} + \rho g\hat{z}_i = 0\comma\; $$

where σ_ij is the Cauchy stress tensor, x _j is the direction in the coordinate system, ρ is the mass density of the ice, g is the gravitational acceleration, $\hat {z}$ is the vertical unit vector, and summation is implied for repeating indices. The Cauchy stress σ_ij is related to the deviatoric stress τ_ij by τ_ij = σ_ij − pδ_ij, where p = (1/3)σ_kk is the mean isotropic pressure. We assume that ice is incompressible, such that ∂u _i/∂x _i = 0, where u _i is the velocity vector.

The momentum balance in Eqn (1) can be vertically integrated to obtain a reduced form of the momentum balance equations, known as SSA. The primary assumptions of SSA are that the vertical shear and bridging stresses (horizontal gradients of the vertical shear stresses) are negligible, ice thickness is much smaller than the length and width of the ice stream (the thin-film approximation), and the normal stresses are cryostatic, meaning they are of the form σ_zz = −ρg(s − z), where s is the height of the ice surface and z is parallel to the gravity vector and positive upward, with z = 0 at the bed. The cryostatic assumption ensures σ_zz is zero at the ice surface to second order. The resulting SSA equations are as follows:

(2)$${\partial\over \partial x}\lsqb H\lpar 2\tau_{xx} + \tau_{yy}\rpar \rsqb + {\partial\over \partial y}\lsqb H\tau_{xy}\rsqb + \tau_{b_x} = \tau_{d_x}\comma\; $$

(3)$${\partial\over \partial y}\lsqb H\lpar 2\tau_{yy} + \tau_{xx}\rpar \rsqb + {\partial\over \partial x}\lsqb H\tau_{xy}\rsqb + \tau_{b_y} = \tau_{d_y}\comma\; $$

where H is the ice thickness, $\tau _{b_i} = \lsqb \tau _{b_x}\comma\; \, \tau _{b_y}\rsqb$ is the basal drag vector and $\tau _{d_i} = \lsqb \tau _{d_x}\comma\; \, \tau _{d_y}\rsqb = -\rho g H\lpar {\partial s}/{\partial x_i}\rpar$ is the driving stress (MacAyeal, Reference MacAyeal1989). The deviatoric stresses are related to strain rates by (Glen, Reference Glen1955)

(4)$$\tau_{ij} = 2 \eta \dot{\epsilon}_{ij}\comma\; $$

(5)$$\eta = {1\over 2}A^{-{1}/{n}} \dot{\epsilon}_e^{\lpar {1-n}\rpar /{n}}\comma\; $$

where $\dot {\epsilon }_{ij} = \lpar {1}/{2}\rpar \lpar \lpar {\partial u_i}/{\partial x_j}\rpar + \lpar {\partial u_j}/{\partial x_i}\rpar \rpar$ is the strain rate tensor, $\dot {\epsilon }_e = \sqrt {\lpar {1}/{2}\rpar \dot {\epsilon }_{ij} \dot {\epsilon }_{ij}}$ is the effective strain rate, A is the flow-rate parameter, η is the dynamic viscosity, and n is the stress exponent commonly taken to be 3, a reasonable assumption for the temperatures and pressures relevant to this study (Jezek and others, Reference Jezek, Alley and Thomas1985; Barnes and others, Reference Barnes, Tabor and Walker1971). Our interest in this study is in the net softening or stiffening of glacier ice, and so we take A, n, and η to be scalars and leave for future work an explicit consideration of anisotropic constitutive relations. Applying the constitutive relation (Eqn (4)), the momentum equations (Eqns (2) and (3)) can be rewritten in terms of dynamic viscosity η (and consequently, the flow-rate parameter A) and velocity gradients:

(6)$$ {\partial\over \partial x}\Big[2 \eta H \Big(2 {\partial u_x\over \partial x} + {\partial u_y\over \partial y} \Big)\Big] + {\partial\over \partial y}\Big[\eta H \Big({\partial u_x\over \partial y} + {\partial u_y\over \partial x} \Big)\Big] + \tau_{b_x} = \tau_{d_x}\comma\; $$

(7)$${\partial\over \partial y}\Big[2 \eta H\Big(2 {\partial u_y\over \partial y} + {\partial u_x\over \partial x}\Big)\Big] + {\partial\over \partial x}\Big[\eta H \Big({\partial u_x\over \partial y} + {\partial u_y\over \partial x}\Big)\Big] + \tau_{b_y} = \tau_{d_y}.$$

Boundary condition

Slip at the bed is characterized by the sliding law, which relates basal velocity $u_{b_i}$ to the basal drag vector $\tau _{b_i}$

(8)$$u_{b_i} = -c \lpar \tau_{b_k} \tau_{b_k}\rpar ^{\lpar {m-1}\rpar /{2}}\tau_{b_i}\comma\; $$

where m is the scalar sliding law exponent and c is the scalar basal slipperiness (representing the lack of resistance the bed provides to the ice slipping over it). The sliding law exponent can vary widely, from negative values to infinity for a perfectly plastic bed (Schoof, Reference Schoof2005, Reference Schoof2010; Minchew and others, Reference Minchew2016). This exponent is commonly assumed to be m = n, which reduces the number of free parameters in the model and represents the viscous flow of ice around obstacles in the absence of Leeward cavitation. In practice, m = 3 describes a sliding law used for hard-bed sliding (Weertman, Reference Weertman1957), currently the most commonly-used in ice-sheet modeling and the value used in this study. Constraining the value of the sliding exponent is an area of active research (Joughin and others, Reference Joughin, Smith and Schoof2019).

Inverse method

The inverse method is implemented in Úa and the setup is summarized here. This method minimizes a cost function J = I + R, where I is the misfit of the surface velocities:

(9)$$ \eqalignb { I \lpar u_x\lpar A\comma\; c\rpar \comma\; u_y\lpar A\comma\; c\rpar \rpar = {1\over 2 \Omega} \displaystyle {\int\int} \Bigg[\Bigg({\vert u_x-u_{x_{{\rm obs}}}\vert \over u_{x_{{\rm err}}}} \Bigg)^{2}+ \Bigg( {\vert u_y-u_{y_{{\rm obs}}}\vert \over u_{y_{{\rm err}}}}\Bigg)^2 \Bigg]{\rm d}x\, {\rm d}y\comma\; }$$

where (u _x,  u _y) are the modeled surface velocities, $\lpar u_{x_{{\rm obs}}}\comma\; \, u_{y_{{\rm obs}}}\rpar$ are the observed surface velocities, $\lpar u_{x_{{\rm err}}}\comma\; \, u_{y_{{\rm err}}}\rpar$ are the formal errors in the observed surface velocities, and Ω is the area of the model domain. These inversions are ill-posed, and thus regularization is needed, particularly in the presence of measurement errors or errors in estimation of ice thickness (Habermann and others, Reference Habermann, Maxwell and Truffer2012). Here, we use a Tikhonov regularization term defined as

(10)$$ \eqalign{ R \lpar A\comma\; c\rpar = {1\over 2 \Omega} \int\int g_{a_A}\lpar A-A_p\rpar ^2 + g_{a_c}\lpar c-c_p\rpar ^2 \cr \quad + g_{s_A} \vert \nabla\lpar A-A_p\rpar \vert ^2 + g_{s_c}\vert \nabla\lpar c-c_p\rpar \vert ^2 {\rm d}x\, {\rm d}y\comma\; }$$

where A _p is the prior value of the flow-rate parameter and c _p is the prior value of the slipperiness coefficient, both of which represent the information known about the parameters prior to the inversion. This form of Tikhonov regularization encourages smooth solutions that are close to the prior values. The regularization coefficients $g_{a_A}\comma\; \, g_{a_c}\comma\; \, g_{s_A}\comma\; \, g_{s_c}$ can be chosen to ensure the regularization terms remain of a similar magnitude to each other and to the misfit term. The regularization coefficient $g_{a_i}$, where i = [A,  c], controls the departure from the prior (A _p,  c _p), and the coefficient $g_{s_i}$ controls the gradient. The minimization problem min _c,A J(c,  A) is solved using a quasi-Newton method (Gudmundsson and others, Reference Gudmundsson, Krug, Durand, Favier and Gagliardini2012). The adjoint method is a computationally effective method of computing the gradient of the cost function with respect to A and c. More details on the use of adjoints and Lagrange multipliers in inverse problems are provided in Joughin and others (Reference Joughin, MacAyeal and Tulaczyk2004) and Morlighem and others (Reference Morlighem, Seroussi, Larour and Rignot2013).

In this study, we set regularization coefficients and minimize the cost function. Once completed, we set the resulting inferences as the starting value and run the inversion again with reduced regularization to keep the regularization balanced with the misfit. Functionally, we are incorporating the information from the resulting inferences into an informed inversion run. With each inversion run, the regularization coefficients $g_{s_i}$ are decreased by approximately an order of magnitude and the coefficient $g_{a_c}$ is reduced by approximately half an order of magnitude. This approach of functionally decreasing the coefficients throughout the inversion ensures that we get close to the minimum misfit and minimum regularization. The initial and final regularization coefficients used in this study are presented in Supplement Table S2.

There is currently no consistent and systematic approach to finding the regularization coefficients in glaciological contexts. Many previous studies have examined potential ways of determining the amount of regularization (Habermann and others, Reference Habermann, Maxwell and Truffer2012), though since this problem involves four coefficients to tune, many of those methods prove intractable for this problem. Thus, the current best approach for tuning these coefficients is a trial approach while attempting to maintain the same amount of regularization for both the flow-rate parameter and the basal slipperiness coefficient, which requires scaling based on the magnitudes of these two parameters. The approach used here begins with a set of regularization coefficients that ensure approximately equal amounts of regularization placed on each parameter and then decreases the value of the coefficients throughout the inversion to ensure a minimum is reached.

However, more work needs to be done in the field of inverse methods and their applications to glaciological problems to find a more complete and systematic way to pick regularization coefficients, particularly in problems with multiple coefficients that require tuning. A worthwhile direction for glaciological inverse problems may be to turn toward data assimilation and variational methods such as those used in meteorological contexts that require no explicit tuning of regularization coefficients (Goldberg and Sergienko, Reference Goldberg and Sergienko2011). This may streamline our approach and potentially enable inversions that are time-dependent (Goldberg and Heimbach, Reference Goldberg and Heimbach2013; Bannister, Reference Bannister2017). These are all interesting and useful directions for future work. For this study, the regularization coefficients are presented in Table S2 of the Supplement and the variation in the coefficients over the ice streams is visualized in Figure S2 of the Supplement.

Synthetic ice stream

To test the inverse methods described above, we constructed a synthetic ice stream that resembles Antarctic Ice Streams in terms of the geometry and mechanical properties. The synthetic ice stream is 20 km wide and 200 km long. The ice is 1 km thick upstream and decreases downstream, similar to ice streams in Antarctica (Fretwell and others, Reference Fretwell2013). The boundary conditions on the lateral margins are no-slip, and the upstream and downstream boundary conditions are given by a prescribed input flux and hydrostatic pressure, respectively (Gudmundsson and others, Reference Gudmundsson, Krug, Durand, Favier and Gagliardini2012). The model domain encompasses the grounded ice stream with the domain boundary located immediately downstream of the grounding line. The resolution of the forward model is higher than the variability of the inferred fields.

The true values of basal slipperiness and the flow-rate parameters are set to mirror natural fluctuations in these parameters (first column of Fig. 1). Since relatively rapid rates of deformation occur in shear margins enabling multiple processes that enhance creep (Jacobson and Raymond, Reference Jacobson and Raymond1998; Suckale and others, Reference Suckale, Platt, Perol and Rice2014; Hewitt and Schoof, Reference Hewitt and Schoof2017; Haseloff and others, Reference Haseloff, Hewitt and Katz2019), we expect ice to be softer in the margins. We approximate the structure of the flow-rate parameter by a transversely-varying cosine curve, in which the flow-rate parameter is highest in the margins of an ice stream and reaches a minimum in the middle of an ice stream. Thus, we set the true value of the flow-rate parameter in the synthetic ice stream to be

(11)$$A_{\rm true} = A_r - {A_r\over 2} \cos\Bigg({2 \pi y\over L_y}\Bigg)\comma\; $$

where L _y is the width of the ice stream (such that −L _y/2 ≤ y ≤ L _y/2) and A _r = 1.6729 × 10⁻⁷ a⁻¹ kPa⁻³, a tabulated value for temperate ice (Cuffey and Paterson, Reference Cuffey and Paterson2010).

Fig. 1. Results from two tests of the inverse method conducted on a synthetic ice stream. The left column shows ‘observed’ velocity (synthetic velocity set as measured velocity in the inversion), the prescribed ‘true’ slipperiness distribution (with slippery spots represented as Gaussian spikes), and the prescribed ‘true’ flow-rate parameter distribution (with high values in the lateral shear margins, where ice is softer due to viscous deformation). The middle and right columns present results of tests, with velocity misfit defined as the difference between the observed and inferred surface velocity, scaled by data errors (in this case, u _err = 1 m a⁻¹). The middle column presents results with spatially constant regularization coefficients, in which the inversion captures the spatial variability of both distributions but there is mixing in the center line of the flow-rate parameter distribution. The right column presents results from an inversion where the flow-rate parameter regularization coefficients are scaled by strain rates (Eqn (14)). Employing strain rates in the inversion reduces the mixing in the flow-rate parameter distribution and improves the estimation of the magnitude of slippery spots in the slipperiness distribution. Figure S2 of the Supplement shows the strain rate field.

The magnitude of the slipperiness coefficient varies spatially. Previous inversions have found sticky patches in which basal slipperiness is low and fast-flowing patches in which basal slipperiness is high (MacAyeal, Reference MacAyeal1992; Joughin and others, Reference Joughin, MacAyeal and Tulaczyk2004). Here, we approximate this spatial variation through Gaussian bumps representing ‘slippery patches’, or localized increases in basal slipperiness:

(12)$$c_{\rm true} = 1 + \sum_i a_i \exp \Bigg\{-{\lpar x-x_i\rpar ^2\over d_x} - {\lpar y-y_i\rpar ^2\over d_y}\Bigg\}\comma\; $$

where a _i are a series of coefficients that determine the magnitude of slipperiness in m a⁻¹ kPa⁻³, x _i and y _i are the spatial coordinates of the Gaussian bump, and d _x,  d _y determine the width of the bump in the x- and y-direction. The values of a _i,  x _i,  y _i are given in Table S1 of the Supplement.

Set-up and results of synthetic experiments

To test whether current inversion methods can resolve ice rheology parameters and basal slipperiness in both space and magnitude, we generate synthetic observations by using the true parameters defined in Eqns (12) and (13) in the forward model. The synthetic observed velocity, which we use without added measurement errors, mirrors the surface velocity of many ice streams in Antarctica: a parabolic across-flow profile arises due to no-slip boundary conditions in the lateral margins and the flow velocity increases toward the grounding line. The surface velocity misfit is defined as the difference between the observed and modeled surface velocity scaled by data errors. For the synthetic tests, there are no added measurement errors to the surface velocity, so to define a surface velocity misfit, we set errors equal to 1 m a⁻¹ everywhere. The surface velocity misfit can be used as a judge of how well the inversion is performed. The inversion has minimized the misfit uniformly enough over the spatial domain to get reasonable estimates of the parameters when the difference between observed and inferred velocity falls below data errors (in the synthetic case, the absolute value of the misfit is less than unity) over the full length of the domain and when the misfit lacks a significant structure.

Prior values represent our prior knowledge of the parameters, both in spatial structure and in magnitude. We set spatially constant prior values: A _p = A _r,  c _p = 1 m a⁻¹ kPa⁻³. For the case of basal slipperiness, this is equivalent to assuming little knowledge of the spatial structure of basal slipperiness and some knowledge of the background field. This is reasonable within Antarctic Ice Streams, since we have borehole measurements to provide some insight into the mechanical properties of the bed, which enables us to estimate the background slipperiness. Most of the uncertainties within many Antarctic Ice Streams lie in the spatial variability and the localized patches of low slipperiness, which our formulation of the prior assumes no knowledge of. For the case of the flow-rate parameter, this is equivalent to having little knowledge of the spatial structure of the flow rate parameter and knowledge of the background field. Again, for the case of Antarctic Ice Streams, this is reasonable since the flow-rate parameter is dependent upon temperature, parameterized through an Arrhenius relation, and thus from ice temperature, we can compute the flow-rate parameter based on tabulated values for ice of a given temperature (Cuffey and Paterson, Reference Cuffey and Paterson2010). We have enough knowledge of ice temperature, particularly at the surface, to enable reasonable estimates of the flow-rate parameter. Therefore, similar to the basal slipperiness parameter, most of the uncertainty lies in the spatial variability.

The results of two synthetic tests are shown in the second and third columns of Figure 1 respectively. The first test is the standard inversion with spatially constant regularization values (shown in Table S2 of the Supplement). The second test considers the use of strain rates in the regularization coefficients as a method of spatially separating the parameters in the inversion (values also shown in Table S2 of the Supplement).

Regularization with strain rates

We now test whether the inclusion of strain rates in the regularization of the inversion might mitigate mixing between the two parameters we are inverting for. We justify this parameterization of the regularization based on our understanding of the dynamics of ice streams. The flow-rate parameter can vary dramatically in the margins of the ice streams due to a variety of mechanisms, including heating and recrystallization, that are enhanced in areas of rapid deformation. Softening of ice in the shear margins can have significant effects on ice stream dynamics. On the other hand, the variability of the basal slipperiness coefficient in the fast-flowing trunk should dominate the contribution of the bed to ice stream dynamics. It is therefore feasible to spatially partition the domain of the ice stream, in effect inverting primarily for basal slipperiness in the center line and inverting primarily for the flow-rate parameter in the shear margins. Since the results of the classical regularization suggested that inverting for both parameters in the same physical location together does not enable separation of the two parameters, this spatial separation should minimize mixing when inverting for both parameters.

In practice, we accomplish this spatial separation by weighting the coefficients for the flow-rate parameter $g_{s_A}\comma\; \, g_{a_A}$ in grounded ice by the inverse of the strain rates. This effectively penalizes changes in the flow-rate parameter in areas with low strain rates, such as near or along the center line, which creates spatial disjointness in the inversion. The inversion puts more weight on the slipperiness parameter in the center line, where most of the changes in the slipperiness parameter are focused, and more weight on the flow-rate parameter in areas of high deformation rates, such as the lateral shear margins. Rheology can be inferred most accurately in rapidly deforming regions, so prioritizing estimates of the flow-rate parameter in small areas of high shear is physically justifiable as well as methodologically sound.

On the ice shelf, there is a negligible basal drag (c → ∞), so the inversion ceases to be a simultaneous inversion and instead becomes a single inversion of only one parameter, the flow-rate parameter. On ice shelves, there is no need for a spatial separation, so we only apply the strain rates in the regularization to grounded ice. We thus adjust the regularization coefficients for the flow-rate parameter by the following:

(13)$$g_{a_A} = {\tilde{g}_{a_A}\over kf + \lpar 1-f\rpar \dot{\epsilon}_{e}}\comma\; $$

(14)$$g_{s_A} = {\tilde{g}_{s_A}\over kf + \lpar 1-f\rpar \dot{\epsilon}_e}\comma\; $$

where $\dot {\epsilon }_e$ is the observed effective strain rate, f is a floatation mask, where f = 1 denotes floating ice and f = 0 denotes grounded ice. The constant k is approximately the same order of magnitude as the strain rates (k = 10⁻² s⁻¹ in this study), and $\tilde {g}_{a_A}$ and $\tilde {g}_{s_A}$ are constants chosen to ensure that the terms ${\tilde {g}_{a_A}}/{k}$ and ${\tilde {g}_{s_A}}/{k}$ are approximately the regularization coefficients that would be used without strain rates, to ensure that the regularization balance is maintained across the grounding line.

Supplement Figure S1 presents our approach, similar to an L-curve approach, to determining these coefficients in which we vary the regularization coefficients used and evaluate the resulting misfit to determine how well the inversion performed. We show that some combinations of regularization coefficients fail to decrease the misfit to a reasonable range, while others do significantly better. The coefficients used in this study are presented in Table S2 of the Supplement and are of the same magnitudes as the regularization coefficients applied to the classical regularization case. This ensures that any changes we see in the inferences between the two cases are due to the inclusion of the strain rates and the resulting spatial variability of the coefficients, rather than being because of differences in the magnitude of the regularization coefficients we use.

The formulation presented in Eqns (13) and (14) faces difficulty in areas of very small strain rates, since the regularization coefficients will then get arbitrarily large. In Antarctic Ice Streams, this is unlikely to be a significant problem due to the rapid rates of deformation and therefore this form of regularization can be used effectively in ice streams. However, this method cannot be used reliably outside of ice streams, in areas of low rates of deformation.

The third column of Figure 1 shows the result of the inversion that includes strain rates in the regularization. The misfit decreases to a small value in comparison to the average velocity. The inversion continues to capture the spatial variability of the slipperiness field when including the strain rates. The errors in some of the peaks of the slipperiness distribution are smaller with the strain rates than without, suggesting that the strain rates may have enabled a more accurate resolution of magnitude. This inversion also captures the spatial variability of the flow-rate parameter, with relatively high values in the margins. There is decreased mixing in the center line, suggesting that the inclusion of strain rates in the regularization partitions the domain into regions where ice is deforming relatively rapidly and the model is more sensitive to values of the flow-rate parameter, and regions where there is little deformation in the ice and thus low sensitivity to the values of the flow-rate parameter. However, there is a systematic overestimation of the flow-rate parameter near the grounding line by about a factor of 2. This is likely due to increasing strain rates towards the grounding line (Supplement Fig. S2), and may be mitigated by a more spatially variable regularization coefficient. We expect that for such a small overestimate (less than a factor of 2), the effect on the results for Antarctic Ice Streams is likely to be minimal. We reserve further exploration of this overestimate for future work.

This synthetic test considers the case of high- frequency spatial variability in basal slipperiness, which is often the case for some Antarctic Ice Streams. However, many ice streams have lower frequency spatial variability in basal slipperiness. We conducted tests on a slipperiness field with smaller (in the magnitude of the peak) Gaussian spikes and a long-wavelength background field that increases the slipperiness closer to the grounding line. The regularization coefficients are presented in Table S3 of the Supplement, the regularization coefficients in Figure S4 of the Supplement, and the results in Figure S3 of the Supplement. The inversion resolves spatial distributions in both parameters, with less dramatic mixing even in the case without strain rates. The minimization of mixing may be because the two parameters vary in different directions (basal slipperiness varies along with the flow and the flow-rate parameter varies across flow) and thus there is an inherent spatial disjointness in the inversion. Because of the lower magnitude of the Gaussian spikes, the inversion is slightly better at resolving variations in magnitude. We consider here only the case of ice streams and reserve for future work a thorough investigation of more complex spatial variations in basal slipperiness, the flow-rate parameter, and the geometry of the glacier that may be expected in areas other than ice streams.

While in these synthetic tests, we approximate ‘slippery patches’ in the basal slipperiness field, in which the magnitude of the Gaussian bumps in the slipperiness field is larger than the mean- field, many studies have found evidence of ‘sticky patches’ in Antarctic Ice Streams, a localized decrease in basal slipperiness. Figure S5 in the Supplement presents the results of a similar test with ‘sticky patches’, in which the Gaussian bumps are less than the mean-field. Here, we see similar results to the test presented here. The inversion captures all of the spatial variations in basal slipperiness, except one misplaced sticky spot. The inversion captured the sharp increase in the flow-rate parameter in the shear margins and the flow-rate parameter field shows little mixing between the estimates, likely due to the inclusion of strain rates in the regularization.

The above tests are conducted without added errors in observed velocities. Surface velocity datasets, widely available, are extensive in coverage and have errors much smaller than the magnitude of surface velocity on Antarctica (Gardner and others, Reference Gardner2018). However, any uncertainties in the surface velocity datasets would impact the results of the inversion. Results from a synthetic inversion using measurement errors suggest that, while errors in surface velocity would add noise to the inferred distribution, the large-scale structure and magnitude of basal slipperiness and the flow-rate parameter would still be quite accurate.

Data of ice thickness and basal topography, on the other hand, have more significant errors (Fretwell and others, Reference Fretwell2013). Uncertainties in ice thickness in particular are likely to affect results of inversions, especially estimates of the flow-rate parameters (Eqns (6) and (7)). This is a limitation of inverse methods as applied to any glaciological problem, and more work needs to be done to determine the effect that current uncertainties in data related to glacier geometry may have on the results of glaciological inversions. We leave for future work an in-depth examination of the effect of noise in ice thickness and basal topography on estimates of the flow-rate parameter and the basal slipperiness coefficient.

Results of synthetic tests presented here suggest that inverse methods can provide accurate estimates of basal slipperiness and the flow-rate parameter in Antarctic Ice Streams. Our synthetic tests assume an ice stream geometry, in which the width is much smaller than the length of the glacier, and no-slip lateral margins that are consistent with Antarctic Ice Streams. Furthermore, our synthetic tests assume certain structures of the flow-rate parameter and basal slipperiness that are found in observations and previous inversions on Antarctic Ice Streams (MacAyeal, Reference MacAyeal1992; Echelmeyer and others, Reference Echelmeyer, Harrison, Larsen and Mitchell1994; MacAyeal and others, Reference MacAyeal, Bindschadler and Scambos1995; Joughin and others, Reference Joughin, MacAyeal and Tulaczyk2004; Arthern and others, Reference Arthern, Hindmarsh and Williams2015; Gardner and others, Reference Gardner2018). The inclusion of strain rates into regularization is a strategy that takes advantage of our knowledge of ice stream dynamics and is thus not necessarily applicable to other types of glaciers and ice flow. Thus, this is not necessarily a universal approach.

In particular, the need for further studies is important where basal slipperiness and ice rheology might be expected to vary with the same spatial pattern. Figures S6 and S7 of the Supplement show the results of a synthetic test in which the true flow-rate parameter distribution is constant and the true basal slipperiness distribution varies along-flow and varies in localized Gaussian spikes as shown in Figure 1, respectively. Since the flow-rate parameter does not increase at the margins, and the changes in strain rate are not due to variations in the flow-rate parameter, this approach presented here does not mitigate mixing as it does in Figure 1. Further testing would be required to ascertain whether these inverse methods apply to other outlet glaciers with different glacier geometries or to other glaciers that may have different structures of basal slipperiness and the flow-rate parameter. However, for the case of Antarctic Ice Streams, we have a good physical justification for the flow-rate parameter being higher in the margins thus the results shown here enable us to make use of inverse methods to improve our understanding of Antarctic Ice Streams.

Model and data

Previous studies have applied inverse methods similar to those described above with ‘classical’ regularization (e.g., Joughin and others, Reference Joughin, MacAyeal and Tulaczyk2004), but increased quantity and quality of satellite data in the last few years enable more accurate inversions. Inversions require data of bed topography, surface elevation, and surface velocity. Ice-penetrating radar data are used to derive bed topography, and in Figure 2c, we show the radar coverage for Bedmap2, overlaid with contours of observed surface velocity (magenta line; Gardner and others, Reference Gardner2018). While extensive in the ice streams, particularly in Bindschadler Ice Stream, radar coverage contains gaps evident in the shear margins. Fretwell and others (Reference Fretwell2013) estimated the uncertainty in the bed elevation to be ~60 m in the Siple Coast. We use surface topography data from the Reference Elevation Model of Antarctica (REMA) (Fig. 2b), which computes surface elevation from satellite imagery (Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019). Here, we use the filled elevation data with a spatial resolution of 200 m. The elevation error is estimated to be ~0.15–1.2 m in the Siple Coast (Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019). Ice thickness (Fig. 2f) is computed from the difference of REMA elevation data (Fig. 2b) and basal topography from Bedmap2 (Fig. 2d) (Fretwell and others, Reference Fretwell2013; Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019).

Fig. 2. Model domain and data used in the inversion: (a) the red outline in the inset shows the location of the model domain, with grounded ice in light gray and floating ice in dark gray. Optical imagery of Bindschadler and MacAyeal Ice Streams from MODIS Mosaic of Antarctica (Scambos and others, Reference Scambos, Haran, Fahnestock, Painter and Bohlander2007), (b) surface elevation from the REMA, (c) radar coverage used to produce bed topography, (d) bed topography from Bedmap2, (e) surface velocity derived from a combination of Landsat 7 and Landsat 8 satellite imagery (Gardner and others, Reference Gardner2018), (f) ice thickness computed as the difference of surface elevation from REMA and bed topography from Bedmap2, (g) effective strain rates ($\dot {\epsilon }_e = \sqrt {\lpar {1}/{2}\rpar \dot {\epsilon }_{ij} \dot {\epsilon }_{ij}}$) computed from surface velocity observations (panel e), (f) driving stress computed from the product of ice thickness (panel f), ice density, and surface slope.

Surface velocities (Fig. 2e) are found from Landsat 7 and 8 satellite data from 2013 to 2015 and the velocity uncertainties are taken on average to be 30 m a⁻¹ (Gardner and others, Reference Gardner2018), which is the error we set for the inversion. The fast flowing regions are moving at ~ 600 m a⁻¹, and the branches of the ice stream at around 300 m a⁻¹. The observed grounding lines are defined from Bedmap2 (Fretwell and others, Reference Fretwell2013) and outlined in green in all panels. The effective strain rates (Fig. 2g), computed from the gradient of the observed velocity fields, demonstrate the high levels of deformation found in the margins. The mean driving stress (Fig. 2h) is ~ 20 kPa and shows a generally decreasing trend with increasing surface velocity, indicating that slip along the bed is the dominant flow regime (Supplement Fig. S8).

The model domain, boundary conditions, and mesh represent the grounded ice streams and proximal regions of the ice shelf. The mesh of the model, constructed in Úa, is refined using observed effective strain rates. We use a relatively fine mesh of ~24000 nodes and 48000 elements, resulting in a spatially varying spatial resolution of a few kilometers (between 1 and 10 km, with finer resolution in areas of high-strain rates). The grounding line in Úa is defined by the floatation condition and closely matches the grounding line of Bedmap2 for the given geometry, though it does not capture some pinning points on the ice shelf because these features are not represented in the Bedmap2 bathymetry (Fig. 2d). The boundary conditions are prescribed to be no-slip at the edges of the grounded model domain, which is chosen to lie only in areas with slow-flowing ice, and are set to be the observed velocities where the model domain is over floating ice.

The inversion is run with spatially constant priors. Based on the tabulated values of the temperature-dependent flow-rate parameter (MacAyeal and others, Reference MacAyeal, Bindschadler and Scambos1995; Cuffey and Paterson, Reference Cuffey and Paterson2010) and an understanding of ice temperature in rapidly-deforming regions of ice sheets (Meyer and Minchew, Reference Meyer and Minchew2018), we estimate the ice temperature to be −10°C, leading to a spatially constant value of the flow-rate parameter A ₀ = 1.15 × 10⁻⁸ a⁻¹ kPa⁻³ (Cuffey and Paterson, Reference Cuffey and Paterson2010). We prescribe the prior of the slipperiness coefficient to be spatially constant at c ₀ = 0.01 m a⁻¹ kPa⁻³. We use strain rates in the regularization coefficients as in Eqn (14) and show the regularization coefficients in the Supplement Figure S9.