1 Introduction

Glaciers have been identified as good indicators of climate change and have been closely monitored over the past two decades. Of paramount interest is the amount of ice contained in glaciers over time. Clearly, in order to establish this quantity, it is necessary to know the location of the glacier’s bedrock, the free surface, and the ice-density distribution. The field measurement of these quantities is time-consuming, costly, and potentially dangerous in such extreme glacial environment. A preferred strategy is to use remote airborne radar sensing. Such techniques are best able to provide information at the free surface of the glacier such as its surface elevation or ice surface velocity. Obtaining reliable information at the bedrock is however considerably more difficult which explains in part why a high-resolution map of the Antarctica ice-thickness distribution is still not available at present. In spite of much recent effort, the resolution of the bedrock data set for Antarctica is still too coarse to capture mountain topographies such as of those of the Antarctic Peninsula or the Transantarctic Mountains with its narrow subglacial valleys and high ice thickness variability [1,2,3,4].

Much is known and understood about how ice flows and there exists a range of mathematical models with different levels of approximation able to reproduce the dynamics of glaciers and ice-shelves [5,6,7,8]. Such models generally require as an input the bedrock location, surface mass balance, ice rheology, and bedrock conditions (geothermal heat flux and bed morphology). Clearly, the higher the level of sophistication, the greater the hope of correctly modelling the ice dynamics. Given all the above inputs, the task of solving for the ice flow dynamics and therefore inferring the ice thickness is generally thought of as the direct problem. On the other hand, the task of identifying one of the unknown inputs from field observations is known as the inverse problem. Such problems are known to be more challenging to solve because the existence of a solution is not guaranteed and if it exists, it is not necessarily unique or stable [9].

We explore here the hypothesis that the free surface is an observable signature of the flow conditions at the bedrock and we propose a numerical algorithm which enables one to decipher this signature, in principle. Specifically, we tackle the inverse problem of identifying the bedrock and the free surface location from the free surface velocity field assuming known basal conditions and surface mass balance. This study builds on our recent work that proposed a new methodology to tackle such an inverse problem in the plane flow regime [10,11].

The idea that free surface features reveal underlying conditions is not new and there have been several important contributions over the years aiming at interpreting these free surface disturbances. The work of Gudmundsson and co-workers has shed significant light on the transmission of bedrock features to the free surface [12,13]. The authors showed that the governing equations could be seen as transfer functions relating the bedrock topography or the basal sliding parameter to free surface observable quantities such as elevation or velocity. By inverting these transfer functions, the authors were able to deduce conditions for which basal perturbations will transfer to the free surface. Early on, Rasmussen described in [14] a method to estimate the bed topography and surface mass balance from aerial photography. The bed topography and surface mass balance are inferred from observed changes in the surface topography and by ensuring that a discrete form of the mass conservation principle is satisfied. Inspired by Rasmussen’s study, McNabb et al. present in [15] a method that allows the computation of the ice depth from the knowledge of the surface mass balance, values of the surface elevation change, the surface velocity field, and an estimated distribution of the ice thickness at the boundary of the domain of interest. The method solves for the ice thickness between adjacent flow lines having to account for cross-flow ice flux in the mass balance equation. In [16], Fastook et al. used ice surface elevation and free surface velocity to estimate the ice thickness and the fraction of the ice velocity due to basal sliding. In [17], Warner and Budd propose a numerical method to infer the bedrock depth using the free surface elevation, the ice-accumulation distribution, ice flow properties, and the assumption of mass balance. This method computes the ice fluxes that would maintain an ice sheet of specified surface elevation in equilibrium with a prescribed ice accumulation distribution. Then, using the standard expression for the magnitude of the depth-integrated ice flux due to ice flow by deformation, the authors are able to calculate a good estimate of the ice thickness distribution for the Lambert glacier basin. The work of Farinotti and co-workers also deserves special mention [18,19,4]. In [18], the authors propose and apply a new approach to calculate the volume and thickness distribution of glaciers. The required input data is the glacier outline and a Digital Elevation Model (DEM). The authors estimate the distribution of apparent mass balance. From this apparent mass balance, the authors are able to estimate the ice flux from which the ice thickness can be deduced using the flow law for ice [20]. This technique allows the authors to estimate the amount of ice stored in glaciers around the globe [19] and to generate a complete bedrock topography map of the Antarctic Peninsula north of 70^oS with a resolution of 100m x 100m, [4]. A related strategy was used by Gantayat et al. in [21] where the authors demonstrate the potential of the method to improve the estimate of ice thickness distribution of glaciers in the Himalayas for which mass balance is not available. Recently, Michel et al. have made several important contributions to the field of bedrock identification [22,23,24]. In [22], Michel et al. introduce the steady inverse method, the quasi-steady inverse method, and the transient inverse method to reconstruct the bedrock of glaciers from the free surface elevation profile using the shallow-ice approximation. The steady inverse method effectively involves solving a regularized, diffusive, transient Partial Differential Equation (PDE), the solution of which approaches the actual glacier bedrock for a suitable choice of the regularization parameter. The quasi-steady inverse method involves replacing the transient term by a suitable estimate which can be computed if the free surface elevation profile is available at two different times. The transient inverse method remedies this drawback by using a regularized fixed point method for which the glacier bedrock is found as the limit of a series of iterations. Michel et al. introduce the shape optimization method in [23] to reconstruct the bedrock of a glacier from the knowledge of the free surface evolution over time using optimal control theory. The method is based on the shallow-ice approximation and effectively relies on finding the “optimal” bedrock profile which minimizes the misfit between the computed and observed free surface elevation. All the methods proposed by Michel et al. are compared in [24]. In [25], Morlighem et al. present a method which allows the generation of a highly resolved bedrock topography map from sparse ice-thickness data and dense ice velocity data. The idea of the method is to optimize the apparent mass balance and the depth-averaged velocity to minimize the error between the observed and modelled ice thickness. The work Haberman et al. in [26] focuses on the reconstruction of basal condition (stickiness) from free surface velocity measurements. The authors used the shallow-shelf approximation to calculate the response of the ice to the basal conditions and iterative methods to minimize the objective function which gives a measure of the misfit between the computed and observed free surface velocity field. In [27], Goldberg and Heimbach derived and implemented the adjoint of an ice flow model which is transient and hybrid meaning that it accounts for vertical shear in the stress balance. It allowed the authors to recover basal lubrication parameters or the simultaneous retrieval of basal topography and sliding coefficient by minimizing the misfit between computed and observed free surface velocity. In [28], van Pelt et al. develop a simple iterative scheme to reconstruct the bedrock topography of glaciers and ice caps. The intuitive idea to assume a particular bedrock profile and run an ice-flow model to calculate the corresponding free surface location. The authors use the PISM ice-flow solver (Parallel Ice Sheet Model, see [29,30]) for this purpose. The bedrock location is simply updated by adding a fraction of the mismatch between the actual and computed free surface elevation. Recently, Martin and Monnier have shown in [31] that it is possible to identify the fluid rheology and basal slip conditions simultaneously from free surface velocity measurements. The authors use the full Stokes equation in this study and the inverse problem is solved in the optimal control framework whereby one aims to minimize the mismatch between computed and observed data. Basal condition of two transantarctic mountain glaciers were also indirectly inferred by Golledge et al. in [32] by using a glacier flowline model constrained by observed ice-surface velocity data.

While undoubtedly successful, the above studies suffer from one or more of the following limitations:

· They use a local relationship between the measured free surface altitude or velocity and the inferred ice thickness. Such a local relationship has a limited applicability.

· They make some assumption about the relationship between the mean ice-velocity and the free surface velocity. Such assumptions have limited applicability.

· They require many iterations to converge to a solution and each iteration typically involves the solution of one or more PDEs which is computationally demanding.

· They use the complex arsenal of optimal control theory which make them difficult to access to the non-specialists.

· They rely on an arbitrary regularization coefficient.

We propose here a new approach inspired by our recent work on topography reconstruction in thin layer flows (see [33] and references therein). The idea behind this approach is to construct a PDE which effectively governs the inverse reconstruction problem. The key advantage of this methodology is that it provides the solution of the inverse reconstruction problem in one-shot in contrast to iterative or optimization-based approaches. The premises of the application of this method to glacier flows were laid out in [10,11] for plane flow but we extend the method here to show that the knowledge of the free surface velocity is sufficient to reconstruct both the free surface and bedrock elevation in three-dimensional domains.

The paper is organized as follows: the next section describes the basis of the mathematical model in the framework of the shallow-ice approximation. The description of the inverse reconstruction method follows and numerical validation tests are presented thereafter. The paper closes with a discussion and concluding remarks.

2 Description of the mathematical model and the inversion strategy

The ice flow dynamics is assumed here to be well-described by the shallow-ice approximation first introduced by Hutter in [34]. This approximation is based on a perturbation expansion of the quantities arising in the mass and momentum balance equations of the ice sheet in terms of the shallowness parameter e expressing the ratio of the characteristic thickness of the ice sheet to its characteristic extent. The main advantage of this approximation is that it transforms the three-dimensional free boundary problem of ice flow modelling into a two dimensional problem for which the unknown is the ice-sheet thickness as a function of position in a fixed two-dimensional domain.

The glacier is assumed to occupy the domain enclosed between the bedrock located at an altitude above a reference location and the free surface described by . The coordinate system is such that the plane is normal to the direction of gravity and points upwards. Clearly, the glacier thickness is given by . The spatial variables and the bedrock altitude, the glacier free surface and the glacier thickness are given in the unit meters in the following. Considering Figure 1 the conservation of mass requires that

$\frac{{\partial H}}{{\partial t}} = a - {\nabla _ \bot } \cdot \overrightarrow {{q_ \bot }}$

(1)

where ${\nabla _ \bot } = \left( {\frac{\partial }{{\partial {\rm{x}}}}, \frac{\partial }{{\partial {\rm{y}}}}} \right)$. The function ${\rm{a}}\left( {{\rm{x}}, {\rm{y}}, {\rm{t}}} \right)$ is the accumulation-ablation function which measures the influence of ice flux normal to the glacier surface due to precipitation or melting. The flux $\overrightarrow {{{\rm{q}}_ \bot }}$ is obtained by integrating the velocity field $\overrightarrow {{{\rm{u}}_ \bot }}$ over the glacier thickness

$\overrightarrow {{q_ \bot }} = \mathop \smallint \limits_{{z_b}}^S \overrightarrow {{u_ \bot }} dz.$

(2)

Within the framework of the shallow-ice approximation and its underlying assumptions, the velocity is given by (see [35])

$\overrightarrow {{u_ \bot }} = \frac{1}{2}A{\left( {\rho g} \right)^3}{\left| {{\nabla _ \bot }S} \right|^2}{\nabla _ \bot }S\left[{{{\left( {S - z} \right)}^4} - {H^4}} \right]$

(3)

where ${\left| {{\nabla _ \bot }{\rm{S}}} \right|^2} = {\left( {\frac{{\partial {\rm{S}}}}{{\partial {\rm{x}}}}} \right)^2} + {\left( {\frac{{\partial {\rm{S}}}}{{\partial {\rm{y}}}}} \right)^2}$, ${\rm{\rho }}$ is the ice density, assumed constant here, is the acceleration of gravity, and is Glen’s flow parameter, a rheological constant. It is assumed here that the glacier is frozen at the base such that no sliding occurs there. Inserting eq. (3) into eq. (2) gives

$\overrightarrow {{q_ \bot }} = - \frac{2}{5}A{\left( {\rho g} \right)^3}{H^5}{\left| {{\nabla _ \bot }S} \right|^2}{\nabla _ \bot }S.$

(4)

Finally, inserting eq. (4) into eq. (1) yields the well known shallow-ice approximation in its standard form

$\frac{{\partial H}}{{\partial t}} = a + \frac{2}{5}A{\left( {\rho g} \right)^3}\nabla \cdot \left( {{H^5}{{\left| {{\nabla _ \bot }S} \right|}^2}{\nabla _ \bot }S} \right).$

(5)

This equation is solved in order to generate synthetic data which can be used for the inverse problem. Throughout the paper we use examples of synthetic glacier bedrock data and construct direct and inverse solutions where the bedrock is defined by

${z_b} = {z_0} - \alpha x + b\left( {x, y} \right)$

(6)

where ${\rm{\alpha }}$ is the dimensionless slope of the bedrock and ${{\rm{z}}_0}$ the elevation at ${\rm{x}} = 0$. In addition to the first two terms which describe a flat incline we assume a further contribution with a function ${\rm{b}}\left( {{\rm{x}}, {\rm{y}}} \right)$ to add a localized topographic structure to the bedrock. The topographic structure is given by

$b\left( {x, y} \right) = \gamma \exp \left( { - \frac{{{{\left( {x - 2000} \right)}^2}}}{{{\delta ^2}}} - \frac{{{{\left( {y - 1000} \right)}^2}}}{{{\delta ^2}}}} \right),$

(7)

where ${\rm{\gamma }}$ is an amplitude, and ${\rm{\delta }}$ the steepness of the topography. Equation (7) describes a localized bump with height ${\rm{\gamma }}$. In the following, the parameters are ${{\rm{z}}_0} = 860, {\rm{\delta }} = 200$ and ${\rm{\gamma }} = {\rm{}} \pm 100$, all with the unit meters and the dimensionless slope is ${\rm{\alpha }} = 0.2$. Within the solution domain ${\rm{x}} \in \left[{0, {\rm{}}4000} \right]$ and ${\rm{y}} \in \left[{0, {\rm{}}2000} \right]$, the topography has the shape as depicted in Figure 2.

Furthermore, we have to define an accumulation ablation function which we set in accordance to [6,11]:

$a\left( {x, y} \right) = {a_0}\left\{ {\begin{array}{*{20}{c}} {(1 - \left( {300 - x} \right)/100}&{, x \le 300}\\ {\left( {2200 - x} \right)/1900}&{, x > 300} \end{array}} \right.$

(8)

where ${{\rm{a}}_0}$ is a constant with the unit meter water equivalent per year $\left[{{\rm{mw}}.{\rm{e}}.{{\rm{a}}^{ - 1}}} \right]$.

In a first step we will derive a discretization scheme for the direct problem which consists in a given bedrock shape, a known ablation function but an unknown glacier surface. All other parameters are constant, see Table 1. We introduce a finite difference scheme with a spatial discretization ${{\rm{x}}_{\rm{i}}} = {\rm{i\Delta x}}$, ${{\rm{y}}_{\rm{j}}} = {\rm{j\Delta y}}$, where ${\rm{\Delta x}}$ and ${\rm{\Delta y}}$ are the spatial grid spacings. For the time discretization we use an explicit Euler scheme and define ${{\rm{t}}_{\rm{k}}} = {\rm{k\Delta t}}$. Using the spatial and temporal discretizations we define ${\rm{H}}_{{\rm{i}}, {\rm{j}}}^{\rm{k}} = {\rm{H}}\left( {{{\rm{x}}_{\rm{i}}}, {\rm{}}{{\rm{y}}_{\rm{j}}}, {{\rm{t}}_{\rm{k}}}} \right)$. Eq. (5) is now rewritten into a diffusion equation with diffusion coefficient ${\rm{D}} = {{\rm{H}}^5}{\left| {{\nabla _ \bot }{\rm{S}}} \right|^2}$:

$\frac{{\partial H}}{{\partial t}} = a + \frac{2}{5}A{\left( {\rho g} \right)^3}\nabla \cdot \left( {D{\nabla _ \bot }S} \right).$

(9)

The spatial discretization in the ${\rm{x}}$-direction is based on a second-order accurate central finite difference scheme where we use a staggered grid in space to allow the computation of the diffusion coefficient ${\rm{D}}$ at grid points which have an offset of ${\rm{\Delta x}}/2$ compared to the grid of ${\rm{S}}$ and ${\rm{H}}$. The following derivatives are all evaluated at a constant time step and we therefore omit the superscript ${\rm{k}}$:

$\frac{\partial }{{\partial x}}\left( {D\frac{{\partial S}}{{\partial x}}} \right) = \frac{{{{\left( {D\frac{{\partial S}}{{\partial x}}} \right)}_{i + \frac{1}{2}, j}} - {{\left( {D\frac{{\partial S}}{{\partial x}}} \right)}_{i - \frac{1}{2}, j}}}}{{\Delta x}}$

(10)

with the definition of the staggered grid:

${\left( {D\frac{{\partial S}}{{\partial x}}} \right)_{i + \frac{1}{2}, j}} = \frac{{\left( {{D_{i + 1, j}} + {D_{i, j}}} \right)}}{2}\frac{{\left( {{S_{i + 1, j}} - {S_{i, j}}} \right)}}{{\Delta x}};{\left( {D\frac{{\partial S}}{{\partial x}}} \right)_{i - \frac{1}{2}, j}} = \frac{{\left( {{D_{i, j}} + {D_{i - 1, j}}} \right)}}{2}\frac{{\left( {{S_{i, j}} - {S_{i - 1, j}}} \right)}}{{\Delta x}}$

(11)

and

${D_{i, j}} = H_{i, j}^5\left[{{{\left( {\frac{{{S_{i + 1, j}} - {S_{i - 1, j}}}}{{2\Delta x}}} \right)}^2} + {{\left( {\frac{{{S_{i, j + 1}} - {S_{i, j - 1}}}}{{2\Delta y}}} \right)}^2}} \right].$

(12)

The derivative in the -direction is written as

$\frac{\partial }{{\partial y}}\left( {D\frac{{\partial S}}{{\partial y}}} \right) = \frac{{{{\left( {D\frac{{\partial S}}{{\partial y}}} \right)}_{i, j + 1}} - {{\left( {D\frac{{\partial S}}{{\partial y}}} \right)}_{i, j - 1}}}}{{2\Delta y}}.$

(13)

Finally, together with the temporal discretization we obtain the following discretization scheme

$H_{i, j}^{k + 1} = H_{i, j}^k + \Delta t{a_{i, j}} + \Delta t\frac{2}{5}A{\left( {\rho g} \right)^3}\left[{\frac{{{{\left( {D\frac{{\partial S}}{{\partial x}}} \right)}_{i + \frac{1}{2}, j}} - {{\left( {D\frac{{\partial S}}{{\partial x}}} \right)}_{i - \frac{1}{2}, j}}}}{{\Delta x}} + \frac{{{{\left( {D\frac{{\partial S}}{{\partial y}}} \right)}_{i, j + 1}} - {{\left( {D\frac{{\partial S}}{{\partial y}}} \right)}_{i, j - 1}}}}{{2\Delta y}}} \right],$

(14)

where the glacier surface is related to the thickness and the bedrock via ${{\rm{S}}_{{\rm{i}}, {\rm{j}}}} = {{\rm{H}}_{{\rm{i}}, {\rm{j}}}} + {{\rm{z}}_{{\rm{b}}, {\rm{i}}, {\rm{j}}}}$.

We note that eq. (14) has to be complemented with appropriate boundary and initial conditions. As initial condition we assume that there is no glacier at ${\rm{t}} = 0$ which is equivalent with ${\rm{H}} \equiv 0$. The choice of initial condition is not as crucial since we seek steady solutions which represent the equilibrium between the accumulation ablation function and the convection of ice. The solutions of the finite difference scheme show that the steady ice thickness is independent of the initial ice thickness for this example bedrock. As boundary conditions we assume periodic boundary conditions in the -direction, namely at ${\rm{y}} = 0$ and ${\rm{y}} = 2000$. This boundary condition is assumed to simplify the numerical code but further boundary conditions in the -direction can be implemented analogously. For the upstream and downstream boundary at the glacier front and back or sides we use an iterative method to track the contact line between ice and the bedrock: For each iteration where ${\rm{H}}_{{\rm{i}}, {\rm{j}}}^{{\rm{k}} + 1}$ becomes negative, we set ${\rm{H}}_{{\rm{i}}, {\rm{j}}}^{{\rm{k}} + 1} = 0$.

The discretization scheme is implemented in Matlab and evaluated until the maximum norm between ${\rm{H}}_{{\rm{i}}, {\rm{j}}}^{{\rm{k}} + 1}$ and ${\rm{H}}_{{\rm{i}}, {\rm{j}}}^{{\rm{k}}}$ drops below 0.1. A typical solution with a discretization with 150 gridpoints in the -direction and 20 gridpoints in the -direction is shown in Figure 3 where we show the solution for the glacier thickness ${\rm{H}}\left( {{\rm{x}}, {\rm{y}}, {\rm{t}}} \right)$ (left column) and the absolute glacier surface position ${\rm{S}}\left( {{\rm{x}}, {\rm{y}}, {\rm{t}}} \right){\rm{}}$ (right column) for different times.

With the knowledge of the glacier surface we are now able to construct a generic example for the inverse problem. Given the surface position of the glacier from solving eq. (5) it is now possible to determine the velocity field at the glacier surface by evaluating eq. (3) at ${\rm{z}} = {\rm{S}}$. The corresponding steady surface velocity field is shown in Figure 4. We note that at ${\rm{x}} \approx 300$ the velocity in the ${\rm{x}}$-direction is zero and the flow direction of the glacier changes direction. This is due to the fact that the ablation accumulation function also changes sign close to this position.

In order to tackle the inverse problem of reconstructing the glacier bedrock from the knowledge of the free surface velocity field, the first step is to rewrite the mass conservation equation in term of the free surface velocity. Considering eqs. (1), (3) and (4), it is straightforward to show that

$\frac{{\partial H}}{{\partial t}} = a - {\nabla _ \bot } \cdot \left( {\frac{4}{5}H\overrightarrow {{u_s}} } \right),$

(15)

where $\overrightarrow {{{\rm{u}}_{\rm{s}}}}$ is the free surface velocity simply obtained by setting ${\rm{z}} = {\rm{H}}$ in eq. (3). It is clear that with the knowledge of the free surface velocity and the surface mass balance, the hyperbolic PDE eq. (15) can be solved for the glacier thickness ${\rm{H}}$. Hyperbolic PDEs suffer from stability issues and special care is required when discretizing and solving this equation. Although we seek steady solutions of (15) we implement a time-dependent discretization method where only the limit for ${\rm{t}} \to \infty$ is of interest. To be independent of discretization we use a stable Euler scheme with implicit time stepping. With the surface velocity $\overrightarrow {{{\rm{u}}_{\rm{s}}}}$, eq. (15) is convective dominant with the convection having different sign over the solution domain as shown in Figure 4. The convection dominance in the ${\rm{x}}$-direction motivates the use of an upwind scheme in the main flow direction. We rewrite eq. (15) into

$\frac{{\partial H}}{{\partial t}} = a - \frac{4}{5}\left[{\frac{{\partial u}}{{\partial x}}H + \frac{{\partial H}}{{\partial x}}u + \frac{{\partial \left( {Hv} \right)}}{{\partial y}}} \right]$

(16)

where $\overrightarrow {{{\rm{u}}_{\rm{s}}}} = \left( {{\rm{u}}, {\rm{v}}} \right)$ are the components of the surface velocity. We discretize the first and the last term with central finite differences

$\frac{{\partial u}}{{\partial x}}H = \frac{{u_{i + 1, j}^{k + 1} - u_{i - 1, j}^{k + 1}}}{{2\Delta x}}H_{i, j}^{k + 1};\frac{{\partial \left( {Hv} \right)}}{{\partial y}} = \frac{{\left( {Hv} \right)_{i, j + 1}^{k + 1} - \left( {Hv} \right)_{i, j - 1}^{k + 1}}}{{2\Delta y}}$

(17)

and use an upwind discretization for $\frac{{\partial {\rm{H}}}}{{\partial {\rm{x}}}}{\rm{u}}$:

$\frac{{\partial H}}{{\partial x}}u = \left\{ {\begin{array}{*{20}{c}} {\frac{{H_{i, j}^{k + 1} - H_{i - 1, j}^{k + 1}}}{{\Delta x}}u_{i, j}^{k + 1}, }&{u_{i, j}^{k + 1} > 0}\\ {\frac{{H_{i + 1, j}^{k + 1} - H_{i, j}^{k + 1}}}{{\Delta x}}u_{i, j}^{k + 1}, }&{u_{i, j}^{k + 1} < 0} \end{array}} \right.$

(18)

Substituting eqs. (18) and (17) into (16) we obtain a set of linear equations for each time step which are solved iteratively. As a first example we take the data from Figure 4 and solve for the glacier thickness. The result for the original data is shown in Figure 5, where we compare the original data with the numerical solution of eq. (15) at two different cross sections. We find that the original data and the inverse solution agree and minor difference can only be found at the steep front of the glacier.

In the reconstruction process, the free surface velocity field $\left( {{\rm{u}}, {\rm{v}}} \right)$ and the glacier thickness ${\rm{H}}$ are now known but the bedrock is still unknown. In the following, we derive a methodology to also solve for the glacier bedrock ${{\rm{z}}_{\rm{b}}}$ and absolute glacier position ${\rm{S}}$. Componentwise, it is clear that the free surface velocity reads

$u = - \frac{1}{2}A{\left( {\rho g} \right)^3}{H^4}{\left| {{\nabla _ \bot }S} \right|^2}\frac{{\partial S}}{{\partial x}}$

(19)

$v = - \frac{1}{2}A{\left( {\rho g} \right)^3}{H^4}{\left| {{\nabla _ \bot }S} \right|^2}\frac{{\partial S}}{{\partial y}}$

(20)

Therefore, it is also clear that

$\frac{u}{v} = \frac{{\frac{{\partial S}}{{\partial x}}}}{{\frac{{\partial S}}{{\partial y}}}} \Rightarrow \frac{{\partial S}}{{\partial y}} = \frac{v}{u}\frac{{\partial S}}{{\partial x}}$

(21)

Replacing the latter expression in eq. (19) and rearranging leads to

$\frac{{\partial S}}{{\partial x}} = \sqrt[3]{{\frac{{ - u}}{{\frac{1}{2}A{{\left( {\rho g} \right)}^3}{H^4}\left( {1 + {{\left( {\frac{v}{u}} \right)}^2}} \right)}}}}$

(22)

Clearly, for a given free surface velocity field and thickness distribution, the ordinary differential equation (22) can be integrated to yield the free surface altitude distribution. The integration involves an integration constant which can in theory be set by knowing the free surface altitude at one particular location. The integration has to be performed in the downstream and upstream direction. It is worth noting that for downslope unidirectional flow, i.e. for ${\rm{v}} = 0$, Glen flow law is recovered.

Continuing with our example above we integrate eq. (22) numerically with the given surface velocity $\left( {{\rm{u}}, {\rm{v}}} \right)$ and the glacier thickness computed from eq. (15). The solution is shown in Figure 6 where the original bedrock data is compared to the inverse solution. The integration constant can be determined at an arbitrary position which is set to the point ${\rm{x}} = 970{\rm{m}}$ in this example. We find a perfect agreement between the original data and the reconstructed data. Only at the upper end of the glacier where the solution is known anyway since ${\rm{H}} = 0$ the integration shows a difference. We note that in Figure 6 the defective point near ${\rm{x}} = 0$ is a result of the discretization of eq. (22) and possibly low values in the denominator lead to high values in eq. (22).

3 Numerical tests

3.1 Example with natural bedrock topography

The example introduced along with the introduction of the inverse methodology has shown that for a simple bedrock topography the reconstruction of the bedrock and the ice thickness distribution based on the surface velocity is successful. In the following the goal is to demonstrate the method on an example where the topography is more complicated. We use the same shape as in eq. (6) but choose a different shape ${\rm{b}}\left( {{\rm{x}}, {\rm{y}}} \right)$, which has smaller scale structures. To maintain the periodic boundary condition in the ${\rm{y}}$-direction the function ${\rm{b}}\left( {{\rm{x}}, {\rm{y}}} \right)$ is also assumed to be periodic. For a higher impact of the bedrock on the glacier surface we reduce the scaling factor of the ablation-accumulation function to ${{\rm{a}}_0} = 2$. The original bedrock shape is shown in Figure 7(a) whereas the corresponding surface velocity field is shown in Figure 7(b).

In the following, we make use of the inverse algorithm and solve for the local glacier thickness and the bedrock shape based on the data in Figure 7(b). The results of the reconstruction are shown in Figure 8. In (a) we compare the original glacier thickness at different cross sections with the inverse solution and find a perfect agreement. The same holds also for the absolute glacier position in (b) where the magnification validates that both solutions coincide. Again as in the previous example the only deviation can be found at the upper end of the glacier where all variables show steep gradients.

3.2 Sensitivity analysis and robustness

The previous results have demonstrated that the inverse method is capable of reconstructing the bedrock topography and the free surface for synthetic topographies with simple bumps but also for a more complicated topography where the surface is corrugated. Nevertheless, all previous examples have in common that the data on which the reconstruction is based is a solution of the direct problem and hence should match the inverse problem. In real applications, the input data in the form of the free surface velocity and the ablation function which are needed to solve the inverse problem as well as all input parameters like ${\rm{A}}$ and ${\rm{\rho }}$ are variables which have to be measured first. The quality of the input data hence depends on the measurement quality and is also affected by experimental errors and possible noise. Since the knowledge of the ice rheology is crucial for the direct problem, and its sensitivity on the solution has already been studied in Martin and Monnier in [30], we focus in the following only on the variation of the field variables ${\rm{a}}$ and $\overrightarrow {{{\rm{u}}_{\rm{s}}}}$ and study the robustness of the inverse solution with respect to a variation thereof. In general, the field variables ${\rm{a}}$ and $\overrightarrow {{{\rm{u}}_{\rm{s}}}}$ are each superimposed with a perturbation of the form

$\begin{array}{l} \tilde a = a + {R_a}\left( {x, y} \right), \\ \widetilde {{{\vec u}_s}} = \overrightarrow {{u_s}} {\rm{}} + {{\vec R}_u}\left( {x, y} \right), \end{array}$

(23)

where ${{\rm{R}}_{\rm{a}}}\left( {{\rm{x}}, {\rm{y}}} \right) \in \left[{{{\rm{a}}_{{\rm{min}}}}, {{\rm{a}}_{{\rm{max}}}}} \right]$ is a perturbation function for the ablation function and ${{\rm{\vec R}}_{\rm{u}}}\left( {{\rm{x}}, {\rm{y}}} \right)$ is a vector valued perturbation for the surface velocity with the first component in the interval $\left[{{{\rm{u}}_{{\rm{min}}}}, {{\rm{u}}_{{\rm{max}}}}} \right]$ and the second component in the interval $\left[{{{\rm{v}}_{{\rm{min}}}}, {{\rm{v}}_{{\rm{max}}}}} \right]$. Each effect in eq. (23) is studied separately and different forms of perturbations are considered. In a first step, perturbations of a short length scale are superimposed such that in the discretization scheme, each point of the function values at the grid points ${{\rm{R}}_{\rm{a}}}\left( {{{\rm{x}}_{\rm{i}}}, {{\rm{y}}_{\rm{i}}}} \right)$ is point wise distributed equally within the interval. This simulates the influence of noise on the inverse solution which could occur in the presence of experimental measurement errors. In a second step, the perturbation is superimposed in the form of a long-scale variation which indicates the influence of systematic measurement errors. The typical shape of the perturbations are summarized in Figure 9, where the perturbed ablation function and velocity component in the ${\rm{x}}$-direction are shown. We note that the maximum amplitude of the long-scale perturbations for the ablation function are reduced compared to the short-scale perturbations since the simulations have shown that the solution of the inverse problem is more sensitive to perturbations of the ablation function than to the velocity field. The perturbations shown in Figure 9 are of the order of 10% to 40% of the maximum value of the ablation function in (a) and in the order of 15% for the velocity field in (b).

For the first step, we consider only the influence of perturbations of the ablation function and velocity field on the reconstruction of the ice thickness. We therefore take again the synthetic example from eq. (6) and (7), solve the direct problem with the given topography and then consider the inverse problem with the perturbed data as in eq. (23). The results of the reconstruction of the glacier thickness are shown in Figure 10. We find from (a) and (c) that for noisy input data where the noise frequency is of the same magnitude as the numerical resolution we can reasonably reconstruct the glacier thickness for most of the cases. The perturbation leads to a noisy solution but the solution remains mostly bounded, except for the perturbation of the velocity field in (c) at the upper end of the glacier. Nearly the same effect is observed for long-scale perturbations in (b) and (d) where deviations only occur at the lower end of the glacier.

The next step involves the reconstruction of the topography from the glacier thickness computed in Figure 10 by solving eq. (22) with the perturbed input data and finally identifying the bedrock by evaluating ${{\rm{z}}_{\rm{b}}} = {\rm{S}} - {\rm{H}}$. The results are shown in Figure 11. We observe that the integration leads to an accumulation of errors and hence the deviation between the original solution and the inverse solution increases in the downslope direction. The integration constant is determined at ${{\rm{x}}_0} = 970{\rm{m}}$ where we assume the glacier bedrock to be known and as a consequence of the summation of errors, the solution quality decreases far away from this location. A very simple way out to significantly improve the quality of the reconstruction is to use another integration point where the bedrock is known which is far away from the initial point. The integration then leads to a solution which has a much better quality in the downslope region. Both solutions are then combined to form a solution over the whole solution domain. It has to be noted that adding further integration points could be used to further improve the quality. An example for an additional point at ${{\rm{x}}_0} = 3000{\rm{m}}$ is added and the combined solution which is obtained by weighting the components is shown in Figure 12. We find that the only major deviations occur at the upper and lower end of the glacier bedrock. The method with adding further points as integration constants and then combining the weighted solution increases the agreement between the original bedrock and the reconstructed bedrock.

4 Conclusions

This study presents a new numerical technique to infer both the free surface elevation and the location of the bedrock from observation of the ice-surface velocity field and the prior knowledge of the surface mass balance and basal conditions. The basis of this numerical technique is the shallow-ice approximation and the corresponding governing PDE. Consequently, the inherent limitations of this approximation also apply to the inverse reconstruction algorithm. The reconstruction algorithm is a two-step process. In the first step, the governing PDE is rearranged such that it takes as an input the ice-surface velocity and solves for the ice-thickness. The second step takes advantage of the explicit relationships between the ice-surface velocity and the glacier surface elevation. When rearranged these relationships lead to an integral equation which when integrated leads to the reconstruction of the glacier surface height. Compared to other techniques developed to solve similar inverse problems, this technique offers the advantage that it delivers the reconstruction through the solution of a single PDE. Many other techniques require an iterative process whereby each iteration involves the solution of a PDE which may end up being a costly process from the computational standpoint. The numerical technique was tested on synthetic data, i.e. data generated by the direct solution of the shallow-ice approximation. A primitive glacier with a bump or a depression and a more realistic “rugged” glacier were tested. For an ideal dataset, free of noise or error, the reconstruction of both the ice-surface elevation and the bedrock is shown to be excellent, as one would expect since the same PDE is solved. Note however that this is not guaranteed when solving such problems using PDE constrained optimization. When noise is added to either the ice-velocity field or the surface mass balance, the accuracy of the reconstruction naturally degrades, particularly for the ice-surface elevation for which the error accumulates due to the integration operator. However, the reconstruction still proves robust with the error in the reconstructed fields bounded and the trend well-captured for error in the input data as high as 40% for the ablation function and 15% for the ice velocity. It is also found that the accumulation of error resulting from the integration can be mitigated by using more than one reference point to fix the integration constant. The true test of the success of this algorithm is with real field data, which will be our future endeavor.