Reconstructing spatially variable mass balances from past ice extents by inverse modeling

Abstract

With the conclusion of the Last Glacial Maximum (LGM), about 20 000 years ago, ended the most recent long-lasting cold phase in Earth history. This last glacial advance left a strong observable imprint on the landscape, such as moraines, trimlines and other glacial geomorphic features. These features reflect the extent of former glaciers and ice caps, which in turn provides information on past temperature and precipitation conditions. Here we present an inverse approach to reconstruct the equilibrium line altitudes (E) from observed ice extents. The ice-flow model is developed solving the mass conservation equation using the shallow ice approximation and implemented using Graphical Processing Units (GPUs). We present the theoretical basis of the inversion method, which relies on a Tikhonov regularization, and demonstrate its ability to constrain spatial variations in mass balance with idealized and real glaciers.

1. INTRODUCTION

The Quaternary period, spanning the last 2.6 million years of Earth history, was marked by warm and short interludes between long and cold periods. While literature provides tight constraints on how climate varied during glacial periods in marine settings (Zachos and others, 2001; Lisiecki and Raymo, 2005; Herbert and others, 2016), changes in temperature and precipitation on continents remain poorly constrained. Use and availability of proxies such as pollen (Peyron and others, 1998) and speleothems (Luetscher and others, 2015) can help shed more light on past climate conditions. However, past glacial advances have left a strong imprint on the Earth's surface in the form of moraines, trimlines and other glacial geomorphic features, which provide us invaluable information on past climate on continents (Penck, 1905; Plummer and Phillips, 2003; Laabs and Carson, 2005). Furthermore, thanks to geochronological methods (Gosse and Phillips, 2001), such as¹⁴C (Hajdas, 2009), cosmogenic nuclide (Blanckenburg and Willenbring, 2014) and optically stimulated luminescence dating (Rhodes, 2011), the timing of past glacial advances can be dated.

Ice-flow models have been used before to infer past climate conditions on continents (Plummer and Phillips, 2003; Anderson and others, 2006; Kessler and others, 2006; Jouvet and others, 2008; Golledge and others, 2012; Rowan and others, 2014; Eaves and others, 2016; Mey and others, 2016). These studies consisted of running ice-flow models forward and fitting observed moraines by optimizing some of the model parameters. A key component in such models is the mass balance, which defines ice accumulation and ablation, and is thus primarily a function of precipitation and temperature (Cuffey and Paterson, 2010). Several models exist to simulate mass balance, from the simple positive degree-day model (Braithwaite, 1995) to more sophisticated energy mass-balance models (Oerlemans, 1992, 1997; Becker and others, 2016; Pellicciotti and others, 2016). These models have been successfully applied to reconstruct glacial history (Becker and others, 2016; Seguinot and others, 2018), but require knowledge of a large number of uncertain parameters, which might change over a glacial cycle. In-depth analysis of ice-flow models used for paleoclimatic reconstructions can be found in the paper by Kirchner and others (2016).

In practice, past ice thicknesses are poorly constrained. Recent studies on the European Alps have emphasized the systematic mismatch between paleo ice-flow model and reconstruction from trimlines (e.g., Cohen and others, 2017; Seguinot and others, 2018), in terms of ice thickness. Here we introduce an inversion method to infer spatially varying mass-balance information from the former extent of glaciers derived from geomorphological observations, such as moraines. We apply a gradient descent method to minimize a cost function, which includes the misfit between modeled and geomorphological maximal glacial extents and a regularization term to ensure sufficient smoothness of the inverted variables. As a result, the method consists of running a forward model a large number of times with a tuning of the mass balance at each iteration. The ice-flow model is modeled by solving the mass conservation equation using the shallow ice approximation (SIA) (Mahaffy, 1976; Hutter, 1983) and the developed codes are accelerated using Graphic Processing Units (GPUs).

Below, we first introduce the basic equations governing our forward model. Then we present the theoretical basis behind the method and provide a step-by-step recipe for its numerical implementation. We then illustrate the method's efficiency for constraining spatial variations in mass balance over mountain ranges with a series of examples.

2. FORWARD MODEL

2.1. Ice dynamics

At the time scale relevant for ice-sheet dynamics, ice can be described as a viscous and incompressible non-Newtonian fluid (Cuffey and Paterson, 2010). Ice thickness can be modeled by solving the mass conservation equation:

(1)$$\displaystyle{{\partial H} \over {\partial t}} =- \nabla \cdot {\bf q} + b,$$

where H is the ice thickness, q is the horizontal ice flux and b represents the mass-balance rate. Horizontal ice flux is defined as vertically integrated velocity field ${\bf q} = \int _{B}^{S}{\bf u}\,{\rm d}z$, and we approximate velocities using the SIA:

(2)$$ \eqalign{{\bf u}(z) &= -\underbrace{A (\rho g)^n [H^{n+1} + (S - z)^{n+1}] |\nabla S|^{n-1} \nabla S}_{\it (i)} \cr &\quad + \underbrace{A_s (\rho g)^n H^{n-1} |\nabla S|^{n-1} \nabla S }_{\it (ii)}, $$

where (i) represents the deformation velocity and (ii) represents the sliding velocity, S is the surface altitude, ρ is density of ice, g is the gravitational acceleration, A is the ice-flow parameter, A _s is a sliding parameter and n is the Glen flow law exponent. Values for the flow parameters A and A _s are taken from Oerlemans (1997) (who used values from Budd and others (1979)), and can be found listed in Table 1, along with the other parameter values used in the model.

Table 1. Ice-flow model parameters

The SIA is based on the assumption that large ice bodies are mainly deformed by vertical shear stress, and therefore neglects longitudinal stresses in the mass conservation equation (Eqn 1). We are using this approximation with the assumption that the horizontal scale of the ice extent is much larger than the vertical scale (H ≪ L). Solving of full set of Stokes equations, a full set of momentum equations, would give more realistic results (especially near the margin and ice divide). However, the ice extent is primarily a function of mass balance, topography and ice dynamics. Since we only consider the extent of ice as forward model output in our inversion method, we restrict our analysis to the SIA.

The mass-balance rate, b, is defined as the annual gain, accumulation, or loss of mass, ablation, at the glacier surface (e.g., Oerlemans, 2001; Cuffey and Paterson, 2010). We define the mass-balance rate, b, as a spatially varying variable:

(3)$$ b=\min(\beta (S(x,y)-E(x,y)),c), $$

where β is the balance rate gradient, E is the equilibrium line altitude (i.e., the location where accumulation is equal to ablation) and c is a maximum ice accumulation rate (e.g., Meier, 1962; Kessler and others, 2006; Oerlemans, 2008). The evolution of the ice cap or glacier surface is in part determined by the mass-balance rate and its values depend on the local climate conditions. Mass-balance rate is positive in the accumulation area and negative in the ablation area. The mass-balance gradient, β, tends to be steeper in the ablation than accumulation area (Mayo, 1984). To keep the problem simple for introducing the inversion method, we choose the balance gradient to be constant.

Solving Eqn (1) using the SIA (Mahaffy, 1976; Hutter, 1983) leads to a non-linear diffusion-reaction type equation that can be solved numerically (Hindmarsh and Payne, 1996). In that case, the ice flux q is defined as

(4)$${\bf q}=-D\nabla{S}, \quad S=H+B, $$

where B is the bedrock and D is the diffusivity expressed as:

(5)$$D=(\rho g)^n[AH^{n+2} + A_{\rm s}H^n]\vert \nabla{S} \vert^{(n-1)}. $$

The main advantage of models based on using SIA to solve the mass conservation equation Eqn (1) is their computational efficiency, which we further enhance using GPUs. The details of the numerical implementation we use can be found in the Appendix A. Note, however, that the inversion method approach we propose can be used with other types of approximations for ice flow, such as solving a full set of Stokes equations or a higher order approximation of these equations (e.g., Pattyn, 2003; Egholm and others, 2011), as long as it is computationally possible. For the moment, running 1000–3000 inversion model runs with higher order models without the use of GPUs is not practical.

3. INVERSION METHOD

In this section, we provide the theoretical basis of the inversion method and its numerical implementation.

3.1. Theoretical basis

Our objective is to invert a known ice extent into a spatially variable E using the forward model explained in Section 2. The forward problem is defined as:

(6)$$h_{\rm m}(x,y) = F(E(x,y),\beta), $$

where h _m(x, y) is the data we are trying to fit (note that h(x, y) is used for ice extent and H(x, y) for ice thickness), F is the forward model (the ice-flow model), and E(x, y) and β are the unknown parameters of the mass-balance equation, Eqn (3), we wish to constrain. Below we derive the method to infer E(x, y) and show it can be extended to β.

As with any inversion method, the primary objective is to minimize the misfit between the model predictions and the observations. To do so, we define a misfit function G(E) as a spatially integrated difference between the ice extent calculated for a forward model, h _m(x, y), and the observations, h _o(x, y) as:

(7)$$G(E) = \displaystyle{1 \over 2}\int_0^{L_x} {\int_0^{L_y} {(h_{\rm o}(} } x,y)-h_{\rm m}(x,y))^2{\rm d}x{\rm d}y,$$

where [0, L _x] × [0, L _y] is our computational domain. The chosen computational domain strictly contains the area where thickness is positive, i.e., we have S = B in the neighborhood of the border of the domain. The goal is then to minimize the misfit function G by varying E(x, y) iteratively.

We now look for a descent direction to define a sequence of E _i, which minimizes the misfit function G. For that purpose, we formally compute 〈G′(E), δE〉, the directional derivative of G with respect to E along a given direction δE:

(8)$$\left\langle {{G}^{\prime}(E),\delta E} \right\rangle = -\int_0^{L_x} {\int_0^{L_y} {(h_{\rm o}(} } x,y)-h_{\rm m}(x,y))\displaystyle{{\partial h_{\rm m}(x,y)} \over {\partial E}}\delta E{\rm d}x{\rm d}y.$$

The term ∂h _m/∂E represents the change of ice extent with respect to E, which are both function of x and y.

The goal of the inversion is to minimize G, and thus make 〈G′(E), δE〉 negative. However, the inverse problem is underdetermined (Hansen, 1994; Zhang and Xu, 2011). We solve this by imposing some form of regularization. Here we assume a smooth solution and impose smoothness on E using a Tikhonov regularization (Tikhonov, 1963; Poplavskii and others, 2001). We do this by defining a new misfit function G, now containing the regularization term:

(9)$$\eqalign{G(E) &= \displaystyle{1 \over 2}\int_0^{L_x} {\int_0^{L_y} \bigg[ } \underbrace{{{(h_{\rm o}(x,y)-h_{\rm m}(x,y))}^2}}_{{(i)}} \cr & \quad + \underbrace{{\tau _2 \bigg[{\bigg(\displaystyle{{\partial E} \over {\partial x}}\bigg)}^2 + {\bigg(\displaystyle{{\partial E} \over {\partial y}}\bigg)}^2}}_{{(ii)}}\bigg]\bigg]{\rm d}x{\rm d}y,}$$

(i) in Eqn (9) is the same as in Eqn (7), while (ii) imposes the smoothness on the solution. The parameter τ ₂ is a positive constant which acts as the additional constraint of smoothness. A large constant would force the solution to be very smooth, while a small or zero constant would cancel this constraint. Following the same derivation as above for (ii):

(10)$${ \left\langle {{G}^{\prime}(E),\delta E} \right\rangle = \int_0^{L_x} {\int_0^{L_y} {\tau _2} } \left( {\displaystyle{{\partial E} \over {\partial x}}\displaystyle{\partial \over {\partial x}}\delta E + \displaystyle{{\partial E} \over {\partial y}}\displaystyle{\partial \over {\partial y}}\delta E} \right){\rm d}x{\rm d}y,$$

and using the integration by parts:

(11)$$\eqalign{& \int_0^{L_x} {\int_0^{L_y} {\left[ {\displaystyle{{\partial E} \over {\partial x}}\displaystyle{\partial \over {\partial x}}\delta E + \displaystyle{{\partial E} \over {\partial y}}\displaystyle{\partial \over {\partial y}}\delta E} \right]} } \,\;{\rm d}x{\rm d}y \cr & \quad = \left[ {\left( {\displaystyle{{\partial E} \over {\partial x}} + \displaystyle{{\partial E} \over {\partial y}}} \right)\delta E} \right]_{0,0}^{L_x,L_y} {\rm -}\int_0^{L_x} {\int_0^{L_y} {\left( {\displaystyle{{\partial ^2E} \over {\partial x^2}}\delta E + \displaystyle{{\partial ^2E} \over {\partial y^2}}\delta E} \right)} } \,{\rm d}x{\rm d}y.} $$

Because the first term on the right-hand side becomes negligible, as there is no ice on the boundaries of the model, we can choose E in the neighborhood of the border such that ∂E/∂x = ∂E/∂y = 0, so 〈G′(E), δE〉 becomes:

(12)$$\left\langle {{G}^{\prime}(E),\delta E} \right\rangle = \int_0^{L_x} {\int_0^{L_y} {\left[ {-\tau _2\left( {\displaystyle{{\partial ^2E} \over {\partial x^2}} + \displaystyle{{\partial ^2E} \over {\partial y^2}}} \right)} \right]} } \delta E{\rm d}x{\rm d}y,$$

which combined with Eqn (8) gives:

(13)$$\eqalign{\left\langle {G^{\prime}(E),\delta E} \right \rangle & = \int_0^{L_x} {\int_0^{L_y} {\left[ {\left( {h_{\rm o}(x,y)-h_{\rm m}(x,y)} \right)\left( {-\displaystyle{{\partial h_{\rm m}(x,y)} \over {\partial E}}} \right)} \right.} } \cr & \quad \left. {- \; \tau _2\left( {\displaystyle{{\partial ^2E} \over {\partial x^2}} + \displaystyle{{\partial ^2E} \over {\partial y^2}}} \right)} \right]\delta E{\rm d}x{\rm d}y.} $$

Following the same logic, one can perform the same development for β instead of E:

(14)$$\eqalign{\left\langle {G^{\prime}(\beta ),\delta \beta } \right\rangle & = \int_0^{L_x} {\int_0^{L_y} {\left[ {\left( {h_{\rm o}(x,y)-h_{\rm m}(x,y)} \right)\left( {-\displaystyle{{\partial h_{\rm m}(x,y)} \over {\partial \beta }}} \right)} \right.} } \cr & \quad \left. {- \; \tau _2\left( {\displaystyle{{\partial ^2\beta } \over {\partial x^2}} + \displaystyle{{\partial ^2\beta } \over {\partial y^2}}} \right)} \right]\delta \beta {\rm d}x{\rm d}y.} $$

3.2. Numerical implementation of the inverse algorithm

To this end, we use the gradient descent method (Press, 2007; Fletcher, 2013), in which E(x, y) is updated with the following increment:

(15)$$\delta E = -\left[ {\tau _1\left( {h_{\rm o}-h_{\rm m}} \right)-\tau _2\left( {\displaystyle{{\partial ^2E} \over {\partial x^2}} + \displaystyle{{\partial ^2E} \over {\partial y^2}}} \right)} \right].$$

This choice of δE leads to a decrease of the misfit G, assuming an appropriate choice of the iteration parameters τ ₁ and τ ₂, which are positive constants. We show below that recombining and substituting Eqn (15) into Eqn (13) results in negative increment of the misfit G:

(16)$$\eqalign{\left\langle {G^{\prime}(E),\delta E} \right\rangle & = \int_0^{L_x} {\int_0^{L_y} {\left[ {-\delta E-\left( {h_{\rm o}(x,y)-h_{\rm m}(x,y)} \right)} \right.} } \cr & \quad\times \left. {\left( {\displaystyle{{\partial h_{\rm m}(x,y)} \over {\partial E}} + \tau _1} \right)} \right]\delta E{\rm d}x{\rm d}y.} $$

Sorting Eqn (16) leads to:

(17)$$\eqalign{\left\langle {{G}^{\prime}(E),\delta E} \right\rangle &= -\int_0^{L_x} {\int_0^{L_y} {{\left[ {\delta E} \right]}^2} } {\rm d}x{\rm d}y \cr &\quad -\int_0^{L_x} {\int_0^{L_y} {\left( {h_{\rm o}(x,y)-h_{\rm m}(x,y)} \right)} } \cr & \quad \times\left( {\displaystyle{{\partial h_{\rm m}(x,y)} \over {\partial E}} + \tau _1} \right)\delta E{\rm d}x{\rm d}y.}$$

One can see that Eqn (17) is negative if the term (h _o(x, y) − h _m(x, y)) is sufficiently small, in other words the method converges if our initial guess is not too far from the final solution.

In practice, we follow the steps of the inversion algorithm presented in Table 2. First, we start by assigning an initial guess for E (step I). We arbitrarily choose an equilibrium line altitude field, E, and add some white noise to avoid imposing a spatial pattern to the solution, and therefore avoid the influence of the initial guess on the solution. Adding random noise enables to test whether the method can retrieve a meaningful solution or not, even with a poor initial guess. Next, the mass conservation equation Eqn (1) is solved using the mass-balance rate field, (b), calculated using the input E, Eqn (3), and the model is run until a steady state is reached, i.e., the ice thickness and ice extent do not change (step II). The output of the forward model, i.e., the ice extent, is used to calculate γ (step III). We define γ as the difference between observed ice extent data and ice extent calculated from a forward model run. Note that although the algorithm works using ice thickness (H), we rarely know it. In that case, we calculate ice extent (h) by setting ice thickness equal to 1 where there is ice, and equal to 0 where there is no ice. Step IV consists of updating the equilibrium line altitude locally by increasing (or decreasing) it where the modeled ice extent/ice thickness was found large (or small). If the glacier is too thick or goes beyond the ice extent, E is raised. Conversely, E is lowered if the glacier extent is too small. We then regularize the problem by applying diffusion on E using τ ₂ (step V) for a chosen number of diffusion iterations, n _sm. Steps IV and V correspond to one step of the descent direction given by Eqn (15). Finally, we return to step II, and iterate until γ is minimized (step VI). In the case where we want to invert for β instead of E, the described steps remain the same. One only needs to replace E with β in the equations since the equations are identical.

Table 2. Implementation of the inversion algorithm

The parameters τ ₁ and τ ₂ are positive scalar values, which must be chosen. Both τ ₁ and τ ₂ are problem dependent. The choice is based mainly on the spatial dimensions of the problem and the spatial variability of the solution.

Furthermore, the choice of the inversion update parameter τ ₁ reflects on the precision of the result. Too high a value will make it difficult for the inversion algorithm to converge to the solution because of too large jumps in value of E after each inversion step. On the other hand, a small value, after applying diffusion, might not result in a change in E that is sufficiently large enough to modify the ice extent of an area, thus not enable to fit the data successfully.

Diffusion in step V represents the smoothness one wishes to impose on E. We do this by iteratively solving the diffusion equation using a fully explicit scheme. The degree of smoothness is set by the number of iterations used in the explicit scheme. A large number of diffusion iterations may lead to a very smooth solution, and might not enable to capture local variations in mass-balance rate. In contrast, a small number of diffusion iterations leads to overfitting the data. We must keep in mind that we are fitting the ice extent by varying only one physical parameter, E, and that differences between observations and our model output could come also from the assumptions made in our forward model. The choice of number of diffusion iteration is illustrated later on.

We would like to stress here, the algorithm shows that the method only requires computing ice thickness and the Laplacian of E.

4. RESULTS

In this section, we illustrate the efficiency of the method with a series of synthetic examples. We start with a simple 2-D problem and then present a 2-D synthetic example using real topography (Uinta Mountains). In the simple 2-D example we show, in separate experiments, that the method can invert the ice extent in order to recover both E and β, while in the Uinta Mountains example we show that one can use ice thickness as well as ice extent to reconstruct E. Finally, we show how the method can be applied to a natural example at the Last Glacial Maximum (LGM) in the South Island of New Zealand. In the following experiments, ice extents are presumed to be reached at the same time, steady state solution, and sliding is parametrized to be constant. The balance gradient, β, is also kept constant in all examples except where we use the method to recover β. Since we can only do our inversion in the area of the map where we have ice, the results will also be presented only within the ice extent.

4.1. Two-dimensional inversion

For this 2-D case, we present two experiments. First, we invert the synthetic ice extent and assess whether the method enables us to recover the prescribed E, and in the second experiment we show that the method is capable to recover β as well. The ice thickness, hence ice extent as well, are simulated by prescribing an arbitrary E and β.

A combination of two simple Gaussian shapes is used for the bedrock:

(18)$$B = B_0\left[ {\exp \left( {-\displaystyle{{X^2} \over {{10}^{10}}}-\displaystyle{{Y^2} \over {{10}^9}}} \right) + \exp \left( {-\displaystyle{{X^2} \over {{10}^9}}-\displaystyle{{(Y-(L_y/8))^2} \over {{10}^{10}}}} \right)} \right],$$

where B ₀ is the maximum bedrock elevation, X and Y represent the matrices of distance vectors, and L _y is the length scale in the y-direction. The prescribed E _o we wish to recover is defined as follows:

(19)$$E_{\rm o} = 2150 + 900\arctan \left( {\displaystyle{Y \over {L_y}}} \right).$$

In this experiment, β is constant and set to 0.01. The initial guess is equal to 3000 m with 400 m amplitude white noise.

In the second experiment, we show that the method enables to recover the mass-balance gradient, β. For this, we use the same setup as above, i.e., bedrock is defined as in Eqn (18), and E as in Eqn (19). We prescribe β _o we wish to recover:

(20)$$\beta _{\rm o} = 0.01 + 0.015\arctan \left( {\displaystyle{X \over {L_x}}} \right),$$

where L _x is the length scale in the x-direction.

Figure 1 shows the calculated E and β (Figs 1a and 1c) and the differences between our inversions and the prescribed values (Figs 1b and 1d). In both cases, during the inversion, convergence is reached when the differences between the forward model output ice extent differ from our synthetic ‘observations’ by <10 points. The results show that convergence is reached after 293 iterations for E and 149 iterations for β. More importantly, the two results demonstrate that it is possible to recover E and β using this approach, and capture the spatial variation that was imposed by Eqns (19) and (20) and illustrated in Figs 1a and 1c. The parameters used in this test are summarized in Table 3.

Fig. 1. (a) shows the calculated (modeled) E using our inversion algorithm, (b) difference between synthetic and calculated E where there is ice (within the ice extent), (c) shows the calculated (modeled) β using our inversion algorithm, (d) difference between synthetic and calculated β where there is ice (within the ice extent).

Table 3. Two-dimensional inversion parameters

4.2. Two-dimensional inversion, Uinta Mountains synthetic test

Now the method is applied to a real topography. In this example, we show that the method is capable of restoring spatially variable E by using the ice extent as well as ice thickness as data. The Uinta Mountains, North America, was used as bedrock topography (Fig. 2a). The topography was extracted from the SRTM data (Jarvis and others, 2008) and smoothened using a median filter.

Fig. 2. (a) shows the bedrock map used (Uinta Mountains) with the synthetic ice extent, (b) is the calculated (modeled) E using our inversion algorithm, (c) difference between synthetic and calculated E where there is ice (within the moraines extent).

First, a forward model was run with a prescribed E _o, which includes some spatial variability. The imposed pattern follows a N–S arctangent form. The ice extent is derived from the calculated ice thickness, which then represents the synthetic observations that we wish to invert for mass balance. The ice extent is calculated by setting the ice thickness equal to 1 where there is ice, and to 0 where there is not. Furthermore, high peaks that pierce through the ice are not included in our misfit calculations. The initial guess for E _m field is 3100 m with 400 m amplitude random white noise added to it. To assess the performance of the inversion scheme, the sum of differences between the synthetic E _o field and the calculated E _m are monitored. All parameters used can be found in Table 4.

Table 4. Uinta Mountains inversion experiment parameters

The results of the inversion for the E _m field are presented in Fig. 2. The differences between the prescribed solution E _o and inversion model solution E _m are <30 m, which falls under 5% difference from the prescribed solution (Fig. 2c).

We now repeat the same test with the same forcing, but instead of using the ice extent we randomly chose ice thickness points on our synthetic glacier. Out of 40 317 ice thickness points within the ice extent we chose 433 (Fig. 9a). The only difference in the approach from the previous example is that now we calculate γ as a difference of ice thickness values in the chosen points between calculated and synthetic ice. The resulting figures (Fig. 9 in Appendix B) and tables with used parameters can be found in Appendix B.

Now we investigate the influence of the number of the diffusion iterations. We report two different end-member cases. The first one uses only five diffusion iterations (Fig. 3a) and the second uses 5000 iterations (Fig. 3b). Fig. 3 shows the differences between the E _o for the two end-member tests. It appears that the choice of the number of diffusion iterations strongly influences the precision of the result. A large number of diffusion iterations reduces the spatial variability in mass balance, and is thus filtering out any local gradient, whereas a small number of diffusion iterations reduces the rate at which the residuals are decreased (Fig. 4).

Fig. 3. (a) is the difference between synthetic and calculated E where there is ice ( within the moraines ice extent) with number of diffusion iterations set to 5, (b) is the difference between synthetic and calculated E where there is ice with number diffusion of iterations set to 5000.

Fig. 4. Residuals calculated as a sum of differences between the synthetic E and the calculated E where there is ice. The red dots represent residuals for the case shown in Figure 3a, green for case shown in Figure 3b and blue the case from Figure 2.

4.3. Inversion of LGM ice extent in the South Island of New Zealand

To verify the method for a real situation, we reconstruct the LGM mass balance using the observed ice extent of the South Island of New Zealand (Barrell, 2011). This example offers us the possibility of comparing the results with previous reconstructions (Golledge and others, 2012).

The South Island's climate is heavily influenced by the Westerlies. As the air flows across the Southern Alps, precipitation is orographically enhanced with up to 10 m a⁻¹ of precipitation on the west side (Henderson and Thompson, 1999). The eastern slopes see considerably lower rates of precipitation (about 1 m a⁻¹) as a consequence of the westerly airflow across the mountain range.

In the past, the Southern Alps have experienced extensive glaciations, including full transitions from glacial to fluvial processes (Whitehouse, 1988). Climate has controlled the growth and decay of mountain glaciers in New Zealand since at least 2.4–2.6 Ma (Suggate, 1990), as testified by an exceptional moraine record across the entire South Island (Barrell, 2011). During the LGM, the North Atlantic meridional overturning circulation was weakened compared with present, the Intertropical Convergence Zone shifted southward, which in the southern mid-latitudes modified the location of wind and temperature (Rojas and others, 2009; Denton and others, 2010; McGlone and others, 2010). The LGM was thus marked by changes in regional atmospheric conditions that affected winds and precipitation in the Southern Alps, which may have led to a change in humidity, possibly up to 25% drier than today (Golledge and others, 2012). Temperatures were between 4^°C and 8^°C lower than today (Newnham and others, 1989; Barrows and others, 2007; Golledge and others, 2012; Seltzer and others, 2015). Most importantly, despite these changes, the orographically enhanced precipitation was preserved during the LGM. This contrast of precipitation is also observed in the position of the equilibrium line altitude, E, on either side of the mountain range (Golledge and others, 2012). Therefore, the spatial variation in E provides us with an appropriate test.

In this inversion, we constrain E(x, y) using the LGM ice extent (Barrell, 2011) shown in Fig. 5a. For simplicity, we use the present day bedrock with all elevations below sea level set to zero. Note, we assume that the ice extent obtained from data is reached at the same time over its entire length. All parameter values used for this test can be found in Table 5.

Fig. 5. Application of the inverse algorithm to the South Island of New Zealand. (a) Bedrock map used (gray) with the LGM ice extent obtained from Barrell (2011) (green), (b) calculated (modeled) E field using our inverse algorithm where there is ice (Test A), (c) E field calculated with the second set of parameters using our inverse algorithm where there is ice (Test B).

Table 5. New Zealand inversion experiment parameters

We present two model results, test A and test B, for two different sets of τ ₁, τ ₂ and n _sm. Both test results (Fig. 5) exhibit a clear north–south gradient in the elevation of E, as well as the west–east gradient. The N–S change could be the result of the temperature change because of the latitudinal differences, making the southern areas colder and with larger accumulation zones, and in turn lower E values. The W–E gradient reflects the influence of the westerly winds bringing precipitation to the west side, thus increasing the accumulation and lowering the altitude of E, while the east side is drier with less accumulation and higher E.

In the second test, test B, we reduce the values of the update parameters, τ ₁ and τ ₂, and the number of diffusion steps n _sm. Reducing n _sm leads to more local variations. Smaller τ ₂ (five times smaller than the first test) reduces the influence of smoothing, while a smaller τ ₁ results in smaller changes in E between iterations, all resulting in a more precise fit. Figs 6a and 6b show the differences between the observed (data) and calculated ice extent (model) for each test. Test B enables to better fit the data, although it does it at the expense of creating local depressions in E. This information is lost in the first test because higher smoothing, n _sm, acts as a filter of local variations. It is worth keeping in mind that we are showing the results of the inversion algorithm varying only one parameter, E. This results in overfitting the data and neglecting the influence of other physical processes ignored in the model as well as approximations made in the forward model, as stated above.

Fig. 6. Differences between the observed and modeled ice extents. (a) Test A. (b) Test B. Yellow represents areas where the data have ice but the model output does not, the red is where the model output calculates ice and the data show there was no ice in those areas.

5. DISCUSSION AND CONCLUSIONS

The main goal of this paper is to present a simple inversion method to reconstruct the spatially variable mass balance and equilibrium line altitude (E) fields using ice extent data. This approach offers the possibility of constraining past climates in glaciated mountain ranges. The inverse method is based on a Tikhonov regularization and computationally, besides the forward model, only requires the calculation of a diffusion equation. In order to show the range of applications of the method, we have constructed six experiments, both synthetic and natural, whose results are presented. To significantly reduce the calculation time necessary for our inverse model, we have implemented the codes to run on GPUs.

In the first experiment, the 2-D synthetic case, we have shown the method is capable of reconstructing a spatially variable E as well as a spatially variable β on a simple topography, consisting of a combination of two Gaussian bumps, using ice extent (Fig. 1). The second experiment is the 2-D synthetic case using real topography (Uinta Mountains) where we show that the method can successfully reconstruct spatially variable E, defined as an arctangent function from N–S, using only ice extent as constraint for our inversion. We first ran our forward ice-flow model with a chosen setup and then used the calculated ice extent as synthetic observations to invert E we originally used. The error while recovering E is less than 5% (Fig. 2c). The same was done in the experiment using ice thickness points presented in Appendix B. The difference is that instead of the ice extent we are using randomly chosen ice thickness points (Fig. 9a). Fig. 9 in Appendix B shows that the method can do the inversion successfully. When analyzing errors of the results obtained using ice extent (Fig. 2c) and ice thickness (Fig. 9c), we can see that using the ice extent the largest error is in the middle of the ice extent and the opposite effect using ice thickness. This happens because the method is more precise in places where it is better constrained by data. In the ice extent case, the errors are smaller on the edges, and in the ice thickness case, errors are larger on the edges and smaller inside the ice because of the distribution of our ice thickness points.

The influence of different number of diffusion iteration has been presented (Figs 3a, 3b). It is clear that too many diffusion iterations will result in losing the spatial variability of the result, while too little iterations will not propagate the update from the edges of the ice efficiently to its middle.

In order to verify the method, we applied it to a real case, the LGM ice extent of the South Island of New Zealand (Barrell, 2011). The assumption made in this experiment is that the entire ice extent is reached at the same time, representing something close to the LGM steady state. Our inversion result for the first test is similar in both pattern and values to the reconstruction done by other authors (Golledge and others, 2012). Fig. 5b shows a clear N–S gradient, which corresponds to the latitudinal change and is thus controlled by temperature, and a clear W–E gradient, which corresponds to the influence of the Westerlies and is thus controlled by precipitation. The second result is an example of overfitting the data. It is important to emphasize that the model setup for this experiment is simplified, which is one of the reasons for the mismatches between the observed and calculated ice extent (Fig. 6b). The bedrock topography used was the current one and the values below sea level were set to 0, thus ignoring calving. Basal sliding was parameterized as a constant value, while basal melting was not taken into account. The influence of isostasy was ignored, which is not a valid assumption for such an ice mass, and is most likely the reason E elevations are higher than expected. Despite these limitations, our method is still able to invert the spatial pattern, both the N–S and W–E gradients, and recover values close to other reconstructions.

The experiments presented show several possible applications. In order to obtain a more accurate reconstruction of mass-balance parameters for large mountain ranges, the forward model needs to be improved to include the effects of isostasy. The choice of bedrock topography for a certain time period remains a setback. The use of present day topography for LGM neglects the change in the amount, distribution and thickness of eroded and deposited sediments since the LGM. Sediment-filled valleys today might have been empty or not even eroded. This influences where our ice can flow in forward model calculations, changing the resulting ice extent and therefore our inverse calculations.

The role of β has not been fully addressed in cases where we attempt to recover E. In the presented experiments the balance gradient β was set to be constant value, and chosen arbitrarily. Future work will include an expansion of the method in order to recover both spatially variable E and β fields. Computationally, adding the inversion of another parameter in the algorithm would only require the calculation of an additional diffusion equation for β.

Similar methods to quantify climate parameters (temperature, precipitation, mass balance, equilibrium line altitude) using ice extents have been developed by other authors (Plummer and Phillips, 2003; Kessler and others, 2006; Golledge and others, 2012; Rowan and others, 2014). The approaches differ in defining the climate forcing (energy balance model, positive degree-day model, net mass balance Eqn (3)), and the ice-flow models used vary in complexity. Some authors use current day conditions or present-day mean annual temperature and precipitation measurements to produce scenarios where these parameters are uniformly changed, temperature reduction by −1^°C or precipitation lowered by 25%, and then compare ice-flow model calculations with observations. The main difference with our approach is that our inversion algorithm works as a spatially variable optimization.

Overall, we show that by calculating only an additional diffusion equation, the method is capable to reproduce the general pattern of a spatially variable E over a mountain range using ice extent data.