1. Introduction
Climate change is forcing the retreat and diminishment of glaciers worldwide (Haeberli and Beniston, Reference Haeberli and Beniston1998; Dyurgerov and Meier, Reference Dyurgerov and Meier2000; Hock and Huss, Reference Hock and Huss2021; Rounce and others, Reference Rounce2023). The associated water runoff, driven by melting glaciers, presents considerable ecological and societal hazards, including flooding, habitat destruction, depletion of freshwater resources and sea-level rise (Meier, Reference Meier1984; Alley and others, Reference Alley, Clark, Huybrechts and Joughin2005; Zhang and others, Reference Zhang, Zhang, Kang and Cong2017). Moreover, glaciers are currently the dominant contributors to current rates of non-steric sea level rise, outpacing the Antarctic and Greenland ice sheets (Zemp and others, Reference Zemp2019; Hugonnet and others, Reference Hugonnet2021). Projections of sea level rise contributions by alpine glaciers are constrained by total glacier volume and are thus sensitive to bed topography (Radić and Hock, Reference Radić and Hock2011; Marzeion and others, Reference Marzeion, Jarosch and Hofer2012). Determining the volume and mass loss of glacier ice is crucial for projecting the precise timing and magnitude of future effects as glaciers continue to shrink worldwide (Jury and others, Reference Jury2019; Rounce and others, Reference Rounce2023; Yang and others, Reference Yang2023). However, current estimates of glacier volume are often poorly constrained due to a sparsity of input data (such as surface mass balance or velocity measurements) (Farinotti and others, Reference Farinotti2019; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022; Hock and others, Reference Hock, Maussion, Marzeion and Nowicki2023) or reliance on equally sparse in-situ ice thickness measurements required for model calibration. Existing methods often disagree on the ice thickness of glaciers worldwide, and this tension has increased with recent studies suggesting that prior methods may be overestimating glacier volume globally, particularly in regions such as High Mountain Asia, and the Greenland and Antarctic peripheries (Farinotti and others, Reference Farinotti2019; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022).
Estimating total glacial surface area or change in terminus position is relatively straightforward using airborne photography or satellite imagery (Paul and others, Reference Paul2013). Glacier ice thickness, by contrast, is notoriously difficult to quantify, with only sparse direct observations of ice thickness available (Welty and others, Reference Welty2020; Hock and others, Reference Hock, Maussion, Marzeion and Nowicki2023). As a result, determining glacial thickness (and volume) remains challenging for most glaciers, especially for glaciers in remote or data-poor regions. Instead, ice thickness is often inferred using surface characteristics of the glacier (Farinotti and others, Reference Farinotti2017, Reference Farinotti2021).
Many existing glacial thickness estimation methods ingest various measurable surface properties, such as surface mass balance (Farinotti and others, Reference Farinotti, Huss, Bauder, Funk and Truffer2009; Huss and Farinotti, Reference Huss and Farinotti2012; Clarke and others, Reference Clarke2013; Michel and others, Reference Michel, Picasso, Farinotti, Bauder, Funk and Blatter2013; Morlighem and others, Reference Morlighem2017; Rabatel and others, Reference Rabatel, Sanchez, Vincent and Six2018) and/or velocity (Gantayat and others, Reference Gantayat, Kulkarni and Srinivasan2014; Rabatel and others, Reference Rabatel, Sanchez, Vincent and Six2018; Sattar and others, Reference Sattar, Goswami, Kulkarni and Das2019) into equations of mass conservation and use this information to invert for ice thickness (Morlighem and others, Reference Morlighem2017). Many of these methods also require quantities such as ice temperature and basal friction (McNabb and others, Reference McNabb2012; Farinotti and others, Reference Farinotti2017; Woodard and others, Reference Woodard, Zoet, Iverson and Helanow2021) that are difficult to measure and are often inferred as part of the inversion. This is especially true of mass balance, where snowfall in mountainous regions can have large variations between adjacent glaciers or even between branches of the same glacier.
An alternative and much simpler approach utilizes the perfect-plastic approximation to estimate glacier thicknesses and underlying bed topography. Recent studies have shown that ice at both hard and deformable beds obeys a nearly plastic sliding law except for when slip velocities are low (e.g. Minchew and Joughin, Reference Minchew and Joughin2020; Zoet and Iverson, Reference Zoet and Iverson2020; Helanow and others, Reference Helanow, Iverson, Woodard and Zoet2021). Moreover, the perfect-plastic approximation has been previously applied to model glacier evolution (Ultee and Bassis, Reference Ultee and Bassis2017; Bassis and Ultee, Reference Bassis and Ultee2019) and has long been used to invert for glacier thickness (Nye, Reference Nye1952; Linsbauer and others, Reference Linsbauer, Hoelzle, Frey and Haeberli2009; Li and others, Reference Li, Ng, Li, Qin and Cheng2012; Frey and others, Reference Frey2014; Chen and others, Reference Chen2022). Although the perfect-plastic approximation does not require mass balance, velocity or temperature as inputs, previous iterations of the perfect-plastic approximation have required a priori specification of the glacier yield strength (Linsbauer and others, Reference Linsbauer, Hoelzle, Frey and Haeberli2009; Clarke and others, Reference Clarke2013; Frey and others, Reference Frey2014; Chen and others, Reference Chen2022). Given that the yield strength may vary significantly between glaciers, this represents an obstacle to using the perfect-plastic approximation in large-scale glacier inversions. To overcome this barrier, we explore a method based on the perfect-plastic approximation that requires no a priori knowledge or parameterizations of the yield strength.
Here, we extend previous methods that use the perfect-plastic approximation and show that, provided the perfect-plastic approximation is valid and glaciers are changing, we only need two easily obtained quantities to determine the bed topography uniquely: two maps of surface elevation at different times or rates of surface elevation change and a reference elevation map. We demonstrate that these two quantities produce self-consistent bed topography and yield strength estimations without requiring the specification of additional quantities. The structure of this paper is as follows: first, we introduce a set of reference glaciers previously used by the Ice Thickness Models Intercomparison eXperiment (hereafter ITMIX) (Farinotti and others, Reference Farinotti2017). Next, we review the perfect-plastic approximation and show how we exploit measurements of surface elevation change to infer both glacier ice thickness and yield strength simultaneously. Finally, applying our method to a set of reference glaciers, we demonstrate that our method produces comparable inversion results to others in the ITMIX experiments (Farinotti and others, Reference Farinotti2017, Reference Farinotti2021), and discuss how our method may offer the potential for regional and global-scale glacier analysis given the increasingly widespread availability of measurements of glacier surface elevation and thickness change (e.g. Hugonnet and others, Reference Hugonnet2021).
2. Studied glaciers
Following ITMIX (Farinotti and others, Reference Farinotti2017), we test our inversion on several reference glaciers: the Svalbard ice cap Austfonna, the Scandinavian glacier Hellstugubreen, the Swiss glacier Unteraargletscher (hereafter Unteraar) and the many glaciers of Mount Elbrus in Russia. We choose this selection of glaciers because Farinotti and others (Reference Farinotti2017) provides the reference surface digital elevation models (DEMs) and change rates in ice thickness (
$\frac{\partial h}{\partial t}$)(Fig. 1) that we require as inputs. We use these surface elevation and thickness change maps, together with point observations of the bed, also provided by ITMIX, to test the fidelity of our inversions. This geographically and dynamically diverse set of glaciers (Table 1) enables us to assess inversion performance by directly comparing results and error statistics with other inversion methods across the same reference glaciers.
Figure 1. Subset of digital elevation models (left, panels a through d) and rates of thickness change (right, panels e through g) for each of the four reference glaciers we use to test our inversion in this study.
Table 1. Reference glaciers and quantities relevant to the presented inversion. All input data are sourced directly from the Ice Thickness Models Intercomparison eXperiment (ITMIX) (Farinotti and others, Reference Farinotti2017).
3. Statistics
We employ Mean Absolute Error (MAE) and Mean Bias Error (MBE) to assess our inversion’s performance against ice thickness measurements obtained from echo-radar ground truth bed elevation measurements obtained from ITMIX (Farinotti and others Reference Farinotti2017, Reference Farinotti2021), defined as,
\begin{equation}
\text{MAE} = \frac{1}{M} \sum^M_{i=1} |b^i - B^i|,
\end{equation}
\begin{equation}
\text{MBE} = \frac{1}{M} \sum^M_{i=1} b^i - B^i,
\end{equation} where 
$M$ is the total number of ground truth bed observations, 
$b^i$ and 
$B^i$ are estimated and observed ice thicknesses, respectively, and superscript 
$i$ denotes individual points of known thickness. MAE quantifies the average variance of the estimated bed from bed depth observations. In contrast, MBE gauges the bias relative to bed elevation measurements, which is useful for quantifying whether the inversion overestimates or underestimates total ice volume.
Normalizing these error statistics by mean glacier thickness enables a more meaningful comparison of errors across glaciers, as it accounts for the relative significance of the errors relative to the size of the glaciers. We introduce the Coefficient of Variance (CV) counterparts by normalizing MAE and MBE by the mean measured glacier thickness 
$\bar{H}$ at the time nearest to when thickness observations were collected,
\begin{equation}
\text{CV(MAE)} = \frac{1}{\bar{H}} \left(\text{MAE}\right),
\end{equation}
\begin{equation}
\text{CV(MBE)} = \frac{1}{\bar{H}}\left( \text{MBE}\right).
\end{equation}Next, we review the perfect-plastic approximation and demonstrate how it can be exploited to estimate glacier ice thickness and bed topography given only repeat elevation data that are becoming increasingly readily available by aerial and satellite imagery.
4. Theory, method, and results
4.1. Review of the perfect-plastic approximation
First described by Nye (Reference Nye1951), Nye’s original treatment of the perfect-plastic approximation assumes that glacier ice can be treated as a rigid body beneath a material yield strength, 
$\tau_y$. Once the shear stress at the bed exceeds the yield strength, the ice deforms rapidly to reduce the stress to the yield strength (e.g. Ultee and Bassis, Reference Ultee and Bassis2016; Bassis and Ultee, Reference Bassis and Ultee2019). However, because in the shallow-ice approximation, the shear stress increases with depth and reaches a maximum at the ice-bed interface, the perfect-plastic approximation is also appropriate when longitudinal stresses in the glacier are small and the sliding law at the ice-bed interface can be approximated as plastic (Ultee and Bassis, Reference Ultee and Bassis2016, Reference Ultee and Bassis2017; Bassis and Ultee, Reference Bassis and Ultee2019; Zoet and Iverson, Reference Zoet and Iverson2020). The perfect-plastic approximation provides a relationship between ice thickness 
$h(x,y,t)$, surface slope 
$\alpha(x,y,t)$ and yield strength 
$\tau_y$ (Nye, Reference Nye1951; Paterson, Reference Paterson1970):
\begin{equation}
{\tau_y} = f\rho g h \sin{\alpha},
\end{equation} where we have suppressed the dependence of variables on position 
$(x,y)$ and time 
$t$ for notational convenience. We define the surface slope, 
$\alpha$, as a positive definite down-slope quantity between 0 and 90 degrees (Appendix A) using the method of finite differences (Figure A1) , 
$\rho$ is the density of ice and 
$g$ is the acceleration due to gravity (Table 2). The shape factor 
$f$ incorporates the effect of lateral geometry (Nye, Reference Nye1965; Paterson, Reference Paterson1970). This has been included in other perfect-plasticity-based inversions (Li and others, Reference Li, Ng, Li, Qin and Cheng2012; Van Wyk de Vries and others, Reference Van Wyk de Vries, Carchipulla-Morales, Wickert and Minaya2022). The effect of the shape factor, 
$f$, however, introduces an additional parameter which, in our study, can ultimately be combined with the yield strength to result in an effective yield strength 
$\tau = \frac{\tau_y}{f}$ to write: (Nye, Reference Nye1951; Paterson, Reference Paterson1970)
\begin{equation}
{\tau} = \rho g h \sin{\alpha}.
\end{equation}Table 2. Relevant parameters used in our inversion. The values of all quantities are considered variables except, 
$\rho=910\;\text{kg}\cdot\text{m}^{-3}$ and 
$g=9.8\;\text{m}\cdot\text{s}^{-2}$.
As first noted by Nye Reference Nye(1952), if we can measure the ice surface and assume a yield strength, we can rearrange Equation (6) to determine the ice thickness 
$h$:
\begin{equation}
h = \frac{{\tau}}{\rho g \sin{\alpha}}.
\end{equation} More usefully, defining the ice surface 
$s(x,y)$ and bottom 
$b(x,y)$ such that 
$h(x,y) = s(x,y) - b(x,y)$, the bed topography 
$b(x,y)$ can be deduced from Equation (7) solely from measurements of surface elevation (e.g. Linsbauer and others, Reference Linsbauer, Hoelzle, Frey and Haeberli2009, Reference Linsbauer, Paul and Haeberli2012; Frey and others, Reference Frey2014):
\begin{equation}
b(x,y) = s(x,y,t) -\frac{{\tau}}{\rho g \sin{\alpha(x,y,t)}}.
\end{equation}The challenge is that, although surface elevation (and slope) can increasingly be measured from airborne and spaceborne measurements, we have to make assumptions about the yield strength to apply Equation (8) to estimate the bed (e.g. Nye, Reference Nye1952; Haeberli and Hoelzle, Reference Haeberli and Hoelzle1995; Linsbauer and others, Reference Linsbauer, Paul and Haeberli2012; Clarke and others, Reference Clarke2013; Frey and others, Reference Frey2014; Chen and others, Reference Chen2022; Van Wyk de Vries and others, Reference Van Wyk de Vries, Carchipulla-Morales, Wickert and Minaya2022). This is problematic because the yield strength can vary between glaciers—or even spatially within a single glacier—in ways that heuristic parameterizations of the yield strength may not capture (Gowan and others, Reference Gowan, Hinck, Niu, Clason and Lohmann2023).
4.2. Evolving glaciers and the perfect-plastic approximation
If both the bed elevation 
$b(x,y)$ and yield strength 
$\tau$ are time-independent, we can interpret Equation (8) as a single equation for both yield strength 
$\tau$ and bed topography 
$b(x,y)$: if we know one of these quantities, we can infer the other. However, for glaciers with measurements of rates of surface elevation or thickness change, we have additional information that we can use to determine both the bed topography and yield strength simultaneously.
In the following, we present a novel method for estimating bed topography and yield strength. The first step in our routine is to approximate the assumption of perfect plasticity as purely local and spatially variable, and to estimate the bed topography and yield strength for every point on the glacier surface where we have two independent elevation measurements. As we shall see, this approach can easily be contaminated by noise in elevation data, yielding spurious negative yield strengths and ice thicknesses. To rectify this, our method has two additional steps: first, we filter out outliers, then impose a spatially uniform yield strength for each glacier, thereby solving an over-determined system using standard least-squares to find a single yield strength for each glacier.
4.3. Step 1: Spatially variable yield strength with no regularization
If we have independent measurements of glacier surface elevations at times 
$t_1$ and 
$t_2$, we can recast Equation (8) in the form of two equations with two unknowns:
\begin{align}
b(x,y) = s(x,y,t_1) -\frac{{\tau(x,y)}}{\rho g \sin{\alpha(x,y,t_1)}},
\end{align}
\begin{align}
b(x,y) = s(x,y,t_2) -\frac{{\tau(x,y)}}{\rho g \sin{\alpha(x,y,t_2)}}.
\end{align} Equations (9) and (10) provide a system of equations for the unknown bed 
$b(x,y)$ and yield strength 
$\tau(x,y)$. This motivates our approach to using the perfect-plastic approximation for a simultaneous inversion of bed topography and yield strength. For glaciers that are evolving (and taking care to avoid regions of the glacier where the glacier surface is nearly flat and surface slope approaches zero), we can use two surface DEMs 
$s_1 \equiv s(x,y,t_1)$ and 
$s_2 \equiv s(x,y,t_2)$ and then solve Equations (9) and (10) to determine bed 
$b(x,y)$ and yield strength 
$\tau$ simultaneously, avoiding the need to specify the yield strength separately. For example, given two DEMs, we can write Equations (9) and (10) in matrix form:
\begin{equation}
\begin{bmatrix}
1 & \frac{1}{\rho g \sin{\alpha_1^i}}\\
1 & \frac{1}{\rho g \sin{\alpha_2^i}}
\end{bmatrix}
\begin{bmatrix}
b^i\\
\tau^i
\end{bmatrix}
=
\begin{bmatrix}
s_1^i\\
s_2^i
\end{bmatrix},
\end{equation} where the superscript 
$i$ denotes co-registered grid points in a pair of DEMs, and subscripts denote the time intervals such that 
$\alpha_1 \equiv \alpha(x,y,t_1)$ and 
$\alpha_2 \equiv \alpha(x,y,t_2)$. Alternatively, defining the difference in glacier surface elevation between times 
$t_1$ and 
$t_2$, 
$\Delta s^i \equiv s^i_2 - s^i_1,$ and then subtracting Equation (9) from (10) and rearranging, we can find a single equation for the yield strength at every grid point:
\begin{equation}
\tau^i = \frac{\rho g \Delta s^i\sin{\alpha^i_1}\sin{\alpha^i_2}}{\sin{\alpha^i_1}- \sin{\alpha^i_2}}.
\end{equation} After solving for the yield strength using Equation (12), we can use Equation (9) or (10) to find the bed elevation at each grid point. Here, we estimate the surface slope 
$\alpha^i$ using central differences with a minimum slope regularization to prevent issues that arise in areas where the slope becomes small (Appendix A).
Figures 2 and 3 show thickness and yield strength results when applied to Hellstugubreen Glacier. The inversion produces nonphysical negative yield strengths in 
$39\%$ of grid points in Hellstugubreen’s interior. This failure implies issues with either some combination of the assumptions made or the data ingested. For instance, our method relies on the perfect-plastic approximation (a relatively crude approximation of ice dynamics) and the shallow-ice approximation, which limits application to large length scales compared to the ice thickness (e.g. Fowler and Larson, Reference Fowler and Larson1978). Similarly, noise in the input data may contaminate slope approximations and inversion results, suggesting a need to regularize the solution to obtain more sensible results. To avoid erroneous yield strengths and ice thicknesses, we filter grid points that produce spurious results in Equation (12) and assume the yield strength is approximately constant for each glacier, therefore constructing an overdetermined system to find a single, glacier-averaged yield strength. We discuss these additional steps in the following.
Figure 2. Negative yield strengths from Equation (12) contaminate ice thickness results with negative ice thicknesses for Hellstugubreen Glacier.
Figure 3. Yield strengths produced via Equation (12). We present grid points that produce negative yield strengths in red, and positive yield strength values falling outside the 15th and 85th percentiles in black. Here, only the grey grid points are ingested into Equation (14).
4.4. Step 2: Filter spurious results
To prevent noise in surface elevation measurements and slope estimates from skewing our inversion, we first use Equation (12) to identify regions resulting in negative yield strengths. We then filter out these areas and all other grid points that produce yield strength values outside the 15th and 85th percentiles from the remaining entries. We visualize this two-step filtering routine in Fig. 3.
4.5. Step 3: Spatially uniform yield strength
The first step in our method solves for the yield strength independently for every grid point. This makes the strong assumption that the perfect-plastic approximation applies locally at the scale of grid spacing of DEMs, and that the yield strength can also vary over this scale. In the second step, we filter out negative and outlier grid points to prevent spurious results. Next, we assume that the plastic approximation is only valid when averaged over a sufficiently large patch of the glacier bed and seek a glacier-averaged yield strength for simplicity.
We drop the superscript denoting grid points from 
${\tau}$ and add an overline to denote a single, glacier-averaged value for yield strength, 
$\bar{\tau}$. We can then construct an objective function for the mean squared error (MSE) mismatch between modeled and observed surface elevation changes,
\begin{equation}
MSE = \frac{1}{N}\sum^N_{i=1}\left[\left( \frac{1}{\rho g \sin{{\alpha^i_2}}} - \frac{1}{\rho g \sin{{\alpha^i_1}}} \right) \bar{\tau} - \Delta {s^i}\right]^2,
\end{equation} where superscript 
$i$ denotes individual grid-points, subscript denotes time, 
$N$ is the total number of geographically co-registered points in a DEM pair. To find the best-fitting glacier-averaged yield strength 
$\bar{\tau}$, we minimize Equation (13) in a typical least-squares fashion (Björck, Reference Björck, Ciarlet and Lions1990). Remarkably, Equation (13) can be minimized analytically. The result is a single equation for our estimate of 
$\bar{\tau}$:
\begin{equation}
\bar{\tau} =
\frac{\sum^N_{i=1} \left( \frac{\Delta s^i(\sin{\alpha^i_1} - \sin{\alpha^i_2)}}{\rho g \sin{\alpha^i_1}\sin{\alpha^i_2}} \right)}{\sum^N_{i=1} \left( \frac{\sin{\alpha^i_1} - \sin{\alpha^i_2}}{\rho g \sin{\alpha^i_1}\sin{\alpha^i_2}} \right)^2}.
\end{equation} We then use Equation (7), (9) or (10) to determine the ice thickness or bed topography from 
$\bar{\tau}$. Because our estimate of 
$\tau$ is now no longer exactly consistent with the assumption of perfect plasticity everywhere, the beds estimated at times 
$t_1$ and 
$t_2$ from Equations (9) and (10) may differ slightly. In this study, we take the average of the two to create a single estimate for the glacier bed.
As a test, we first applied this approach to Hellstugubreen. Figure 3 shows the spatial distribution of negative and extreme yield strength values filtered out and not considered in Equation (14). Figure 4 panel h shows that when using only the filtered yield strengths, we estimate a physically plausible bed, whereas Equation (12) alone outputs negative yield strengths and ice thicknesses.
Figure 4. Bed topography (left, panels a through d) and ice thickness measurements (right, panels e through g) for the reference glaciers we use to test our inversion using Equation (14).
We next apply Equation (14) to four reference glaciers included in ITMIX: Austfonna, Elbrus, Hellstugubreen and Unteraar. Summary statistics for the reference glaciers are shown in Table 3. In Fig. 4, panels a through d display maps of inferred bed topography, and panels e through g show ice thickness for all the reference glaciers. Together, these illustrate how our inversion produces bed and thickness estimates without negative or spuriously large ice thicknesses in the inversion. The mean basal yield strength across all glaciers is 
$\sim$110 kPa, which is broadly consistent with previous results in the literature (Nye, Reference Nye1951; Fischer, Reference Fischer2012; Li and others, Reference Li, Ng, Li, Qin and Cheng2012; Clarke and others, Reference Clarke2013; Zoet and Iverson, Reference Zoet and Iverson2020; Van Wyk de Vries and others, Reference Van Wyk de Vries, Carchipulla-Morales, Wickert and Minaya2022).
Table 3. Estimated yield strength, mean ice thickness, and error statistics for each of the four reference glaciers.
4.6. Validation with bed observations
Next, we validate ice thickness estimates against ice thickness measurements for each glacier from Farinotti and others (Reference Farinotti2017) to quantify the errors associated with our inversion. The number of thickness point measurements 
$N$ used for validation, together with other relevant glacier data, is presented in Table 1.
Figure 5 shows observed versus estimated ice thicknesses for each of the four reference glaciers. Here, we observe that errors in our inferred ice thickness vary glacier by glacier. Although Hellstugubreen, our best-performing glacier, follows the zero-error line relatively well, our inversion tends to overestimate smaller ice thicknesses and underestimate large ice thicknesses for Austfonna. Elbrus, however, shows no clear trends but produces noisy results near the median ice thickness. We see from Fig. 5 and Table 3 that our inversion produces erroneously thin ice thickness estimations for Unteraar, with a large negative bias showing severe underestimation of ice thickness (
$-88\%$ mean bias), and the yield strength would need to be much greater to match observations. This under-estimation problem is not unique to our inversion method; other approaches in the ITMIX (Farinotti and others, Reference Farinotti2017) suite struggle similarly to produce sufficiently thick ice thickness measurements for Unteraar in particular. Next, we compare our results to those of other methods in the ITMIX suite.
Figure 5. Error distribution heatmaps for each reference glacier studied, along with all points of known glacier bed depth.
4.7. Comparison to existing ice thickness inversions
We compare our results to several ice thickness inversions to determine how well our simple inversion procedure compares to other inversions that participated in ITMIX (Farinotti and others, Reference Farinotti2017). We see that—consistent with the statistics summarized in Table 3—our inversion performs comparably to other ice-thickness inversions, despite our method’s simplicity. Our MAE and MBE fall within the IQR or outperform other inversions in three of the four glaciers studied, which we visualize in Fig. 6 using box-and-whisker plots. Figure 7 shows how our estimated bed compares to other inversions for a handful of glacier profiles. We demonstrate that our method produces bed geometry profiles comparable to other beds submitted to ITMIX.
Figure 6. Box-and-whisker plots of errors for all inversion methods submitted to the Ice Thickness Models Intercomparison eXperiment (ITMIX) (Farinotti and others, Reference Farinotti2017), together with our perfect-plastic approximation (PPA) bed errors for each of the reference glaciers considered. The blue box represents the IQR of model errors, and whiskers extend to the minimum and maximum non-outlier values. The number of ITMIX inversions per glacier are: Austfonna = 7, Elbrus = 6, Hellstugubreen = 8, Unteraar = 15 inversions.
Figure 7. Panel a (10x vertical exaggeration): cross-sectional profile for Elbrus transverse to flow. Panel b (3x vertical exaggeration): cross-sectional profile for Hellstugubreen parallel to flow. A and B denote the transects’ starting points, with A’ and B’ being endpoints, respectively. In both transects, the blue line is the ice surface, the brown area is the known bed reconstructed from point observations of the bed, the red is the estimated bed profile created by our inversion and the unlabeled grey lines are other ice thickness inversions submitted to ITMIX (Farinotti and others, Reference Farinotti2017).
5. Discussion
The perfect-plastic approximation-based inversion we introduce provides a simple method to estimate glacier thickness, bed topography and basal yield strength using solely surface elevation measurements. Despite our method’s simplicity, we produce estimates of basal yield strength that agree with the current literature (O’Neel and others, Reference O’Neel, Pfeffer, Krimmel and Meier2005; Cuffey and Paterson, Reference Cuffey and Paterson2010; Li and others, Reference Li, Ng, Li, Qin and Cheng2012, Ultee and Bassis, Reference Ultee and Bassis2016, Reference Ultee and Bassis2017; Zoet and Iverson, Reference Zoet and Iverson2020; Helanow and others, Reference Helanow, Iverson, Woodard and Zoet2021), and bed topographies comparable to those obtained with more sophisticated inversion methods (Farinotti and others, Reference Farinotti2017, Reference Farinotti2021). We expect the perfect-plastic approximation to be most valid for glaciers on either hard or deformable beds where ice at the bed is sliding rapidly enough that sliding is approximately plastic (Zoet and Iverson, Reference Zoet and Iverson2020; Helanow and others, Reference Helanow, Iverson, Woodard and Zoet2021). By contrast, we anticipate that nearly static glaciers, which largely evolve due to surface mass balance with little ice dynamics, and surging glaciers, where the yield strength changes appreciably over the observational period, may challenge the perfect-plastic approximation. However, despite our test glaciers not neatly fitting into these categories, we find that our method performs as well as other methods for most glaciers tested, with near-zero mean-bias error when compared to bed observations for Austfonna, Elbrus and Hellstugubreen. This near-zero mean bias is likely a result of our choice to use a single yield strength for the entire glacier, which consequently means that we underestimate in some regions and overestimate in others, in ways that may cancel each other out in Equation (2). Nevertheless, this may be beneficial for regional-scale analyses, where the average accuracy over entire regions is preferred over minimizing errors for individual glaciers, as shown in other model description papers (e.g. Huss and Hock, Reference Huss and Hock2015; Maussion and others, Reference Maussion2019). Our inversion does have some limitations, however, which we discuss in the following.
5.1. Inversion limitations
Our assumptions of a plastic rheology, unchanging bed topography, and static yield strength dictate particular relationships between observed changes in glacier surface elevation and slope, as explored in Equations (12) and (14). When these relationships are not satisfied, the inversion may fail to estimate yield strength and ice thickness accurately. This is the main drawback to Equation (12), where every point is treated independently. We suspect this is a consequence of noisy input data and/or applying the perfect-plastic approximation over length scales where the perfect-plastic approximation does not apply. Thus, applying the perfect-plastic approximation to a larger suite of test glaciers may better reveal conditions when the perfect-plastic approximation is likely to fail and provide erroneous results.
The steps explored in this paper, whether treating every point independently or considering all points of the glacier simultaneously, are only limiting cases in a spectrum of possible yield strength regularization routines. A potential drawback to our inversion is that a single glacier-averaged yield strength may not be sufficient to resolve the bed geometry of some particularly complex glaciers. We suspect that including a more generalized regularization that allows for spatially variable yield strength may yield more accurate estimates of ice thickness and bed topography for some glaciers. Still, we cannot discount the potential that our inversion method may perform poorly in some instances. A perfectly plastic rheology may not be suitable in some cases, particularly for glaciers that are surging. Additionally, it is possible that glaciers with changing beds (for example, Columbia Glacier, AK (Boldt Love and others, Reference Boldt Love, Hallet, Pratt and O’Neel2016)) may not satisfy the condition of time-independent bed geometry imposed in Equations (9), (10) and (11) in some areas, which may contribute to errors in bed topography for these glaciers. However, our inversion errors for the bed are often on the order of tens of meters; thus, even in regions with extreme sediment deposition rates, the assumption of a stagnant bed may not be entirely unjustified. Furthermore, given increased satellite coverage, it is possible that future inversions could improve upon ours by using multiple snapshots to infer a time-varying bed.
A major assumption of the perfect-plastic approximation is that subglacial topographical features manifest on the surface. Our inversion works in the reverse direction, with every surface characteristic being linked to subsurface geography. Therefore, estimates of bed topography and ice thickness can easily be contaminated by noise and errors in surface elevation measurements. To quantify the effect of surface elevation error on our results, we ran a suite of simple perturbation tests on our input data for the glacier Hellstugubreen.
5.2. Sensitivity tests
To better understand the effect of noisy input data, we perform a simple sensitivity analysis, shown in Table 4, by perturbing all of Hellstugubreen’s input surface data grid cells (10-meter resolution) with uniformly distributed vertical noise. We introduce an error of between 
$\pm$ 1 meter to the reference surface elevation map and 
$\pm$ 0.1 meters per year to rates of thickness change. We find that heavy surface contamination can significantly decrease the estimated yield strength, resulting in thinner ice estimations that, on average, increase inversion error. We suspect noise in surface elevation measurements may be responsible for the spuriously thin ice thickness estimates for Unteraar, and that smoothing the surface maps before ingestion may improve our estimated glacier bed.
Table 4. Effects of adding uniform noise to surface inputs for Hellstugubreen Glacier on yield strength, Mean Absolute Error (MAE) and Mean Bias Error (MBE).
5.3. Data quality and surface smoothing
We find that smoothing the ITMIX surface elevation data before ingesting it into our inversion, or using different surface inputs entirely, recovers thicker ice thicknesses for Unteraar, thereby reducing the error relative to observations. We smoothed ITMIX surface measurements at both times using a Gaussian-weighted filter with a standard deviation and filter radius comparable in size to Unteraar’s mean observed ice thickness (about 230 m). Smoothing surface inputs in this manner results in a significantly improved bed for Unteraar, with error metrics MAE and MBE reducing to 
$-103$ meters (
$44\%$) and 
$-84$ meters (
$-36\%$), respectively. These improved metrics place our bed within the IQR of all Unteraar ITMIX inversions instead of being the worst-performing method. Of course, we do not necessarily know the ice thickness in advance for many glaciers, and hence, developing a more systematic regularization/smoothing approach may be necessary for estimating glacier ice thicknesses.
Upon ingesting alternative surface elevation measurements (ASTER GDEM 3.0) and thickness change products from Hugonnet and others Reference Hugonnet(2021) into our inversion without smoothing, we again find that the MAE and MBE metrics for Unteraar drop significantly, with MAE decreasing to 
$81$ meters (
$35\%$) and MBE to 
$-33$ meters (
$-14\%$). These analyses, combined with the reality that most methods from the ITMIX suite similarly struggle with Unteraar, suggest that our difficulties with Unteraar may not necessarily reflect sensitivity to characteristics specific to Unteraar, but to something inherent in the ITMIX data.
5.4. Suitability for global-scale glacier analysis
We believe our inversion method is well-suited for global studies of glacier ice thickness and potential sea level rise, owing to the relative abundance of readily available surface topography products that can be used as inversion inputs. There exist several DEMs providing worldwide surface elevation coverage (TanDEM-X, ASTER GDEM, NASADEM, SRTM, Copernicus DEM), regional DEMs (ArcticDEM, EuroDEM, USGS The National Map, among many others) and glacier thickness change products (Hugonnet and others, Reference Hugonnet2021), that together provide surface elevation measurements for practically all glaciers on Earth.
Moreover, our inversion can perform regional and global-scale computations covering entire regions in the Randolf Glacier Inventory at significantly less computational expense than many other inversion methods, needing only to solve Equation (14) once per glacier inverted. This computational affordability provides the flexibility to perform experiments that may be difficult for other global-scale inversions, such as sensitivity tests where we perturb inputs, interpolate data to different resolutions and even perform bootstraps to assess uncertainty.
6. Conclusions
We present a novel extension of the perfect-plastic approximation to infer glacier ice thickness, apply the approach to reference glaciers, and validate results against ice thickness observations. Our inversion method offers promising results for reconstructing subglacial bed topography and estimating glacier ice thickness. We demonstrate that with the perfect-plastic approximation, we create estimates of glacier ice thickness and bed elevation comparable to more sophisticated methods using minimal input data—only two measurements of glacier surface elevation separated in time. This makes our inversion particularly useful for poorly studied glaciers in remote areas and for glaciers that are changing. Given that the only data required are surface DEMs (which are becoming increasingly readily available) and georeferenced glacier outlines (available through the Randolph Glacier Inventory (Pfeffer and others, Reference Pfeffer2014)), we believe this method is suitable for large-scale global glacier analysis. We seek to expand upon this research by applying our inversion method globally to catalog ice thickness estimates for glaciers worldwide, combining knowledge of surface, slope, and yield strength features to infer subglacial conditions and their interactions in an ever-evolving cryosphere.
Acknowledgements
I thank my advisor and co-author, Jeremy Bassis, along with everyone in the Vintage Glaciology Group, for their continued mentorship, support, and helpful discussions. I also thank the wonderful physics faculty at Centre College for encouraging me to pursue a doctoral program in the sciences. We thank Lizz Ultee and another unnamed reviewer for their insightful feedback during the preparation and review of this paper. This work is funded by DoE grant C3710, Framework for Antarctic System Science in E3SM, NASA grant 80NSSC22K0378 by the DOMINOS project, a component of the International Thwaites Glacier Collaboration (ITGC). Support from National Science Foundation (NSF: Grant 1738896) and Natural Environment Research Council (NERC: Grant NE/S006605/1).
Appendix A. Appendix - Numerics
This section contains technical details that underlie our model development.
A.1. Surface slope and regularization
For any given surface 
$s(x,y,t)$, we find the surface slope angle 
$\alpha(x,y,t)$ via approximation at each discrete point on the ice surface using central finite differences,
\begin{equation}
\begin{aligned}
\frac{\partial s(x,y,t)}{\partial x} \approx \frac{s(x+\Delta x,y,t) - s(x - \Delta x,y,t)}{2\Delta x},\\
\\
\frac{\partial s(x,y,t)}{\partial y} \approx \frac{s(x,y+\Delta y,t) - s(x,y-\Delta y,t)}{2\Delta y},
\end{aligned}
\end{equation} where we define 
$x$ as the east-west direction, 
$y$ as the north-south direction and 
$\Delta x$ and 
$\Delta y$ are grid cell spatial resolutions in the 
$x$ and 
$y$ directions in a Cartesian coordinate system, respectively.
Figure A1. The setup for approximating surface slope using central finite differences described in Equation (A1) is visualized.
Near the glacier edges, where central difference gradient estimation is not possible, a simpler forward or backward difference is taken, depending on which is necessary. Equation (A1) can be written as a vector,
\begin{equation}
\boldsymbol{\nabla s} = \begin{bmatrix}
\frac{\partial s^i}{\partial x} \\
\frac{\partial s^i}{\partial y}
\end{bmatrix}.
\end{equation} Taking the arctangent of the magnitude of the gradient vector in Equation (A2) finds the slope angle 
$\theta$ that the surface makes relative to the horizontal.
\begin{equation}
\theta^i = \arctan||\boldsymbol{\nabla s}|| = \arctan \sqrt{\left(\frac{\partial s^i}{\partial x}\right)^2 + \left(\frac{\partial s^i}{\partial y}\right)^2}.
\end{equation} Following other studies (e.g. Farinotti and others, Reference Farinotti, Huss, Bauder, Funk and Truffer2009, Li and others, Reference Li, Ng, Li, Qin and Cheng2012), we regularize all surface slope values by imposing a minimum slope limit (which we define as 
$\theta_0 = 3$ degrees) to avoid spuriously large ice thicknesses caused by nearly flat surface slopes. We define regularized slope values as,
\begin{equation}
\alpha^i = \sqrt{({\theta^i})^2 + ({\theta_0})^2},
\end{equation} This method of regularization ensures that only small slopes are adjusted significantly. If 
$\theta^2 \gt \gt {\theta_0}^2$ then 
$\alpha^i \approx \theta^i$ and this regularization routine has little effect. In the methods section of this paper, all slope 
$\alpha$ values used and presented are regularized using this routine. Although our slope limiting constant could theoretically be tuned on a glacier-by-glacier basis, we define a single slope lower limit for simplicity. While individual glaciers may be sensitive to the choice of the lower limit of the slope, we find that it is not a significant factor when considering numerous glaciers in aggregate.
A.2. Surface input
We used all surface input data sourced from the ITMIX intercomparison project (Farinotti and others, Reference Farinotti2017), which provides a DEM and gridded thickness change rates for each of the reference glaciers we present. A second surface elevation map can easily be produced by adding thickness change rates to a preexisting map of surface elevations,
\begin{equation}
s^i_2 = s^i_1 + \frac{\partial h^i}{\partial t}\Delta t,
\end{equation} where 
$s_2$ is surface elevation, 
$\frac{\partial h^i}{\partial t}$ is the rate of thickness change measured over the elapsed time 
$\Delta t = t_2 - t_1$. Assuming a non-changing bed, the thickness change is equal to the surface elevation change,
\begin{equation}
\frac{\partial h^i}{\partial t} = \frac{\partial s^i}{\partial t},
\end{equation} which is necessary to compute Equation (A5). Provided the time 
$t_1$ for surface elevation map 
$s_1$ is contained within the range which 
$\frac{\partial h^i}{\partial t}$ is valid, the choice of elapsed time 
$\Delta t$ used in Equation (A5) is arbitrary and makes no significant impact on inversion results so long as 
$\Delta t$ is valid within the provided surface measurement range.
Open access
