Bedrock reconstruction from free surface data for unidirectional glacier flow with basal slip
Elizabeth K. McGeorge, Mathieu Sellier, Miguel Moyers-Gonzalez, Phillip L. Wilson
mmanuscript No. (will be inserted by the editor)
Bedrock reconstruction from free surface data forunidirectional glacier flow with basal slip.
Elizabeth K. M c George · Mathieu Sellier · MiguelMoyers-Gonzalez · Phillip L. Wilson
Received: date / Accepted: date
Abstract
Glacier ice flow is shaped and defined by several properties, including the bedrockelevation profile and the basal slip distribution. The effect of these two basal properties canpresent in similar ways in the surface. For bedrock recovery this makes distinguishing be- tween them an interesting and complex problem. The results of this paper show that in somesynthetic test cases it is indeed possible to distinguish and recover both bedrock elevationand basal slip given free surface elevation and free surface velocity. The unidirectional shal-low ice approximation is used to compute steady state surface data for a number of syntheticcases with different bedrock profiles and basal slip distributions. A simple inversion methodbased on Newton’s method is applied to the known surface data to return the bedrock pro-file and basal slip distribution. In each synthetic test case, the inversion was successful inrecovering both the bedrock elevation profile and the basal slip distribution variables. These
E. K. M c GeorgeSchool of Mathematics and StatisticsUniversity of CanterburyChristchurch, New ZealandORCID: 0000-0001-5707-640XM. SellierDepartment of Mechanical EngineeringUniversity of CanterburyChristchurch, New ZealandORCID: 0000-0002-5060-1707M. Moyers-Gonzalez ( (cid:0) )School of Mathematics and StatisticsUniversity of CanterburyChristchurch, New ZealandTel.: +64-64-3-364-2987E-mail: [email protected]: 0000-0003-4817-1506P. L. WilsonSchool of Mathematics and StatisticsUniversity of CanterburyChristchurch, New ZealandORCID: 0000-0002-4563-9399 a r X i v : . [ phy s i c s . g e o - ph ] J un Elizabeth K. M c George et al. results imply that there is a unique bedrock profile and basal slip which give rise to a uniquecombination of free surface velocity and free surface elevation.
Keywords
Glacier · Ice flows · Inverse problems · Shallow Ice Approximation · Basal slip
Understanding cryosphere dynamics is key to modelling climate change. The contributionof land ice to global mean sea level (GLMS) rise for medium emissions scenarios is pro-jected to be at least 0.10 m with some models predicting a contribution of up to 0.27 m[10]. Cazenave et al. (2013) identified one of the main contributors to this rise as the meltingof glaciers. Glaciers are also important socially, with millions of people in the Himalaya,Karakoram and Hindu Kush mountains relying on glacial reserves for their drinking water[36]. Given the potentially large impact of glacier dynamics on human livelihood, compre-hensive glacier models are needed. In particular, accurate methods for calculating the totalice mass of glaciers are required. If the bedrock profile of the glacier is known, the resultingice thickness can be used to calculate the mass of the ice for that particular glacier. Havingexplicit knowledge of glacier mass can be useful and influential in policy and natural re-source planning. However, due to the difficulty of measuring the bedrock profile explicitlyin many cases, it is desirable to use surface measurements and an inversion model to predict the bedrock elevation.Surface elevation and free surface velocity data are already recorded for many ice flowsand glaciers. A number of parties collect and collate data such as the World Glacier Inven-tory (WGI), the United States Snow and Ice Data Center (NSIDC) and the Global LandIce Measurements from Space (GLIMS) initiative. Data is collected in a variety of ways,primarily in-situ or via air- or space- borne craft [15]. Due to the large availability of sur-face data, bedrock recovery methods using these free surface measurements are particularlypopular. Another feature that can be measured or calculated from surface data is the ac-cumulation/ablation distribution for the glacier. This distribution describes how the glaciergrows/diminishes over time dues to snow/ice accumulation/ablation over time. This paperassumes this distribution to be measurable, though it may difficult and time consuming to doso. Field measurements can be costly [27, 23]. Accurately predicting the accumulation ratefrom other measurable surface variables is an area of research in and of itself and many dif-ferent methods employing a vast array of techniques have been proposed [e.g. 13, 29, 4, 37].Glaciers exhibit gravity-driven creep flow which is sustained by the underlying slopedgeography. Glacier ice is categorised as an incompressible, nonlinear, viscous, heat con-ducting fluid [21] which can be described mathematically by the full Stokes flow equationstogether with rheological laws. Many methods of approximating the Stokes flow equationshave been proposed in the last century. One of the most widely used approximations isthe shallow-ice approximation (SIA) [24, 18]. In the SIA model, gravity-driven ice flow issolely balanced by basal drag neglecting longitudinal and transverse stresses, as well as ver-tical stress gradients [1]. Due to the complex nature of glacier ice flow, recovering the glacierice thickness from only surface measurements is a non-trivial inverse problem. Variations inrecovered glacier ice thickness can be as large as the ice thickness recovered for differentmodels. The recovered thickness is also very sensitive to input data [16]. In addition, inver-sion methods can have ill-defined solutions and may impose too many assumptions, such asthe no-slip condition at the base [5, 40, 2, 19, 22].Imposition of a no-slip condition simplifies the inverse problem significantly and al-lows much faster computation. However, basal slip is known to be influential in the flow edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 3 behaviour [26] and so is important to include if possible. Since the primary driving force ofglacier flow is gravity, flow speed is modulated by presence, or lack thereof, of friction at theglacier-bedrock boundary [14]. In a temperate glacier, where high temperatures cause melt,or a thick glacier, where increased regelation causes melt, the glacier-bedrock interface iswet which can cause the ice to slide along the interface easily. Conversely in a glacier whichhas a frozen base, the ice flow is stuck to the ground and does not slide [7]. Other factorssuch as till composition also impact the amount of friction at the base. Increased veloc-ity at the base results in a lower steady state glacier surface due to a process called dynamicthinning; the loss of ice due to accelerated ice flow into the ablation zone [6, 17, 35, 38]. Dy-namic thinning can also be caused by a steeper bedrock profile simply due to the increasedcontribution of gravity on the glacial flow.Since the free surface of an ice flow is affected by both basal slip and bedrock topogra-phy, separating the effects of these two factors in the recovery is difficult [32, 33]. This paperseeks to accurately recover the bedrock topography of a synthetic glacier with non-constantbasal slip given known free surface elevation and velocity. The method proposed builds onthe work of [19] by modifying their method to include basal slip.An overview of the governing ice flow model used is given in Sect. 2. Section 3 con-structs the synthetic glacier surface for a number of different synthetic cases. The resultsaffirm that both basal slip distribution and bedrock profile have a significant effect on the re-sultant steady state surface elevation and free surface velocity. Section 4 gives the derivation of the recovery method proposed and the results of implementing this are given in Sect. 5.A brief sensitivity analysis of the method to noisy surface data is given in Sect. 6. Finally,the results are discussed in Sect. 7 and final conclusions are drawn in Sect. 8.
This paper assumes the glacier flow dynamics are well described by the SIA. The SIA sim-plifies the full Stokes equations by performing a scaling analysis to obtain dimensionlessfield equations for the glacier flow. The small parameter used assumes the glacier extent ismuch larger than its thickness. Some properties of the SIA model are; (1) longitudinal andtransverse stresses, as well as vertical stress gradients vanish,(2) the horizontal componentof the velocity points in the direction of steepest descent of the free surface and does notchange with depth, and (3) domes or troughs have no horizontal velocity. Blatter et al. [8]advise caution when applying the SIA to processes on smaller scales where the assump-tions may no longer be valid, for example, anisotropic basal sliding or locally steep basaltopography. Despite potential drawbacks, the SIA is used widely in ice flow modelling asit reduces a three dimensional flow problem into a two dimensional problem. This makescomputationally simple in comparison to higher order models such as (ADD: more higherorder models) the full Stokes where a full force balance has to be calculated at each step.The SIA is typically set up with x -direction along the flow, the y -direction as the transversedirection, and the z -direction as the upward direction normal to the gravitational field. Tosimplify the testing of the new inversion method, the SIA is restricted to the unidirectionalcase which omits the transverse flow. Elizabeth K. M c George et al.
Table 1:
Notation.
Symbol Meaning z Vertical axis, represents height above a reference elevation x Horizontal axis, distance along glacier from an upstream reference t Time H Glacier height S Glacier surface z b Bedrock profile elevation a Accumulation ablation profile for the glacier u Velocity profile of the glacier u b Basal velocity u s Free surface velocity τ Stress τ b Basal shear stress β Basal slip distribution H , is related to the surface S and the bedrock elevation z b via H = S − z b (1)at any time, t . See Fig. 1 for a pictorial description of this relationship.By considering the momentum balance, volume flux, and mass conservation of theglacier, the expression for height evolution in the unidirectional case is ∂ H ∂ t = a ( x ) − ∂∂ x q x , (2)where a ( x ) is the accumulation/ablation function of the glacier in meters of water equivalentper year, and q x = (cid:90) Sz b u x dz (3)describes the ice flux by integrating the velocity of the ice along the x -direction, u x , from thebedrock to the free surface. Following Gessese et al. [19] and adapting to include basal slipvelocity u b , the velocity profile is given by u x ( z ) = A ( ρ g ) (cid:18) ∂ S ∂ x (cid:19) (cid:2) ( S − z ) − H (cid:3) + u b ( x ) , (4)where ρ is the ice density, g is the acceleration due to gravity and A is the creep or flowparameter given in Table 2. Values for these constants are given in Table 2. The value for ρ is taken as the midpoint of the range for glacier ice as recommended by Cuffey and Patterson[14, Table 2.1].The no-slip condition classically imposed [5, 40, 2, 19, 22] forces u b = edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 5 xz Bedrock, z b Surface, S Height, H Accumulation, a Surface speed, u s Basal slip, β Fig. 1:
Glacier flowing downstream with surface S and bedrock z b and thickness or height H . Surface speed, u s is indicated at the glacier surface and basal slip, β is indicated at theglacier base. Accumulation, a , is represented as falling snow.system has only one unknown to recover. However, as discussed in the introduction, basalslip can have significant effect on glacier height which reduces the practical applications ifit is neglected. Here, no such condition is imposed and the glacier is allowed to have variedbasal slip along the base of the glacier.Weertman [39] first proposed a power-type law for basal shear on a hard bed and bothFowler [18] and Lliboutry [31] proposed a more general form of the law for a flow withcavity formation. Budd et al. [9] found this generalised form to be empirically true for iceflow with basal shear described by τ b = A s u b , (5)where τ b is the shear stress, A s the sliding constant given in Table 2, and u b the basal velocity.The value for A s is taken from Gessese et al. [19].Considering also Glen’s empirical law for the shear rate [20] and Nye’s adaption of thelaw to ice flow [34], at the ice flow base gives ∂ u ∂ z (cid:12)(cid:12)(cid:12)(cid:12) z = z b = A τ b . (6)The value assigned to A depends strongly on temperature. The value given in Table 2 is foran ice sheet at − C and was recommended by Cuffey and Patterson [14, Table 3.4]. The Elizabeth K. M c George et al.
Table 2:
Typical values of constants used throughout.
Symbol Name Value A s Sliding coefficient 5 × − m N − yr − A Glen’s law parameter 4.16 × − Pa − yr − ρ Ice density 880 kg m − g Gravitational acceleration 9.81 m s − derivation of the SIA gives τ b = − ρ gH ∂ S ∂ x , (7)which combines with (5)and (6) to give the following expression for basal velocity u b ( x ) = − β ( x ) A s ( ρ g ) H ( x ) (cid:18) ∂ S ( x ) ∂ x (cid:19) , (8)where β ( x ) is the basal slip distribution which regulates the amount of basal slip at theglacier base. Basal slip is restricted such that β ( x ) ∈ [ , ] for all x in the glacier domain.Physically, β ( x ) = β ( x ) = β ( x ) to be constant along the glacier length. Combining (4) and (8) gives a full expression for the velocity profile. This velocityprofile is substituted into (3) to give the ice flux. Finally, substituting this ice flux into themass balance gives a non-linear diffusion equation ∂ H ∂ t = a + ( ρ g ) ∂∂ x (cid:18) D ∂ S ∂ x (cid:19) (9)with non-linear effective diffusion coefficient D given by D = (cid:12)(cid:12)(cid:12)(cid:12) ∂ S ∂ x (cid:12)(cid:12)(cid:12)(cid:12) H (cid:20) AH + β A s (cid:21) . (10)Note that the full velocity profile easily gives an expression for the surface velocity by setting z = S : u s = − ( ρ g ) (cid:18) ∂ S ∂ x (cid:19) H (cid:18) AH + β A s (cid:19) . (11) Investigating the inverse problem requires synthetic surface data for basic ice flows. Thisdata is produced by solving (9) for the steady state using a slight modification of the finitedifference scheme as laid out by Gessese et al. [19]. The scheme discretizes Eq. (9) in timeusing an Euler explicit scheme and spatially using a second-order accurate central finitedifference scheme. The results is H n + i = H ni + ∆ ta i + ∆ t ∆ x ( ρ g ) (cid:26) ( D ni + + D ni ) (cid:18) S ni + − S ni ∆ x (cid:19) − ( D ni + D ni − ) (cid:18) S ni − S ni − ∆ x (cid:19)(cid:27) , (12) edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 7 where D i = H i (cid:12)(cid:12)(cid:12)(cid:12) S i + − S i − ∆ x (cid:12)(cid:12)(cid:12)(cid:12) (cid:18) AH i + β A s (cid:19) . (13)The scheme is implemented in forwards time with a mesh size of ∆ x = ∆ t = . S i = ( z b ) i + H i . The scheme runsuntil a steady state is achieved as defined bymax i | H n + i − H ni | < .
001 m . (14)This steady state represents an equilibrium between the accumulation and ablation due toa steady flow of ice. Note that the scheme was tested for stability with three other gridresolutions to confirm stability; ∆ x =
50 m, ∆ t = . ∆ x =
20 m, ∆ t = .
016 yrs, and ∆ x = ∆ t = .
001 yrs. This steady state surface is computed for a number of benchmarkcases as outlined in Subsect. 3.2.In each case the glacier starts with no height, H ( x , ) =
0, and has Dirichlet boundaryconditions, H ( , t ) = H ( L , t ) =
0. These conditions represent that there is no height at boththe top and bottom of the glacier.3.1 Benchmark casesCombinations of different accumulation/ablation rate, basal slip distribution, and bedrock elevation profile are used as benchmark cases for testing the methods. ( x ) For each benchmark case the accumulation/ablation function is defined as a ( x ) = (cid:40) a (cid:0) − − x (cid:1) if x ≤ a (cid:0) − x (cid:1) if x ≥
300 (15)where a is the maximum value of the accumulation/ablation function and set to 0.5 for allfuture calculations. Adjusting this maximum values simply raises or lowers the steady statesurface [30]. This function gives the most accumulation at the top end of the glacier andthen linearly decreases along it’s length until at the bottom end which has net ablation. β ( x ) Three different types of basal slip distributions are tested. The different distributions againhelp to test the robustness of the method against more realistic scenarios. The three type arelabelled and mathematically defined as follows. Examples of their shapes are given in Fig.2.1. β ( x ) constant.This type of slipping gives constant slip along the glacier as described by β ( x ) = K , (16)where K is a constant in [ , ] . K =
0, for example, could represent a cold-based glacierin which the ice is frozen to the bed which inhibits the flow mechanisms of sedimentdeformation, ice deformation and basal sliding [28, 3, 25, 11].
Elizabeth K. M c George et al. β ( x ) bump.This type of slip gives negligible slip at the top and bottom of the glacier with narrowregion with lots of slip in the middle as described by β ( x ) = γ e − ( x − M ) δ , (17)where M is the midpoint of the bump, γ is the height of the bump and δ is the extent.3. β ( x ) step.This type of slipping gives a transition from negligible basal slip to some slip as de-scribed by β ( x ) = R + e k ( x − M ) , (18)where lim x → ∞ β ( x ) = R and lim x →− ∞ β ( x ) = k is the steepness and M is the midpointof the transition between asymptotes. B a s a l s li p d i s t r bu t i on , ( x ) (x) = 0(x) = 0.5(x) bump(x) step Fig. 2:
Examples of the three forms of basal slip distribution, β ( x ) . For β ( x ) bump, γ = , M = , δ = β ( x ) step, L = , R = , M = , k = . b ( x ) Two different bedrock elevation profiles are defined in the benchmark cases. The two formsare given by: edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 9
1. An inclined flat bed, z b = f z b ( x ) = f ( x ) = z − α x (19)where z is the elevation at the point x = α is the slope.2. An inclined bumpy bed, z b = bz b ( x ) = b ( x ) = z − α x + A sin ( λ x ) ; (20)where z and α are as before, A is the amplitude of the bumps and λ wavelength. For allsimulations, A =
50 and λ = . This bedrock profile is designed to combine two testcases as used by Gessese et al. [19] to test their no-slip recovery method.In each case z =
900 m and α = . β ( x ) =
0, is plotted along side the constant slip case and is consistent with the results of Gessese et al. [19]. These two figures illustrate thatincreased basal slip results in: (1) a glacier with less height and, (2) an elongated glacierdomain. These two effects are consistent with literature [35, 6] and occur for the other typesof basal slip also but are less visibly obvious. Figure 4 shows that the glacier surface profilefollows that of the bedrock.In Fig. 5, a visible dip occurs in the glacier surface at the location of increased basalslip for the bump case. Similarly, the glacier surface is observed to lower once the transitionoccurs to more basal slip in the step case. In Fig. 6 the same properties are exhibited due tochanges in basal slip but are much harder to see due to the surface undulations produced bythe bedrock.These modelled surfaces show that the shape of the glacier is affected by both thebedrock and the basal slip. Without accounting for basal slip, the dips observed in the sur-face will appear to be the result of related dips in the bedrock. These results for the directcase reinforce the importance of including basal slip in the recovery method for bedrockelevation.
The inverse problem seeks to recover the bedrock elevation profile, z b ( x ) and the basal slipdistribution, β ( x ) , for a steady state glacier from two known free surface quantities; (1) thesurface elevation, S , and (2) the free surface velocity, u s . Previous authors [19, 22] haverecovered bedrock data from one free surface input with the assumption of a sticky, no-slipbase where β ( x ) =
0. Using two input variables allows for the recovery of bedrock elevationprofile and basal slip distribution simultaneously. As with the direct problem, the flow isassumed well described by the unidirectional SIA defined in Eq. (9) and the surface velocityapproximated by Eq. (11).Glaciers and ice flows have been studied extensively throughout the past century by avariety of different groups such as NSIDC. Data collection has been pushed particularly due c George et al. S u r f a c e e l e v a t i on ( m ) S(x) for (x) = 0S(x) for (x) = 0.5z b (x) Fig. 3:
Computed steady state glacier surfaces for constant basal slip with flat inclinedbedrock, z b = f .to ice melt acting as a major contributor to sea level rise [12] (Fig. 13.10, 13.13). As suchthe body of data for glaciers and ice flows is ever increasing. It is reasonable to assumethere is, or can be measured, sufficient data for the two free surface variables as well as theaccumulation ablation function.4.1 Method for the inverse problemGiven two observable variables, u s and S , the following will show that it is possible toaccurately recover to unknown variables, H ( x ) and β ( x ) . Recovery of H ( x ) immediatelygives the desired bedrock due to Eq. (1). To solve for these two variables, first consider theequations which define them. The steady state surface is given by Eq. (9) where ∂ H / ∂ t isset to 0, 0 = a + ( ρ g ) ∂∂ x (cid:32)(cid:12)(cid:12)(cid:12)(cid:12) ∂ S ∂ x (cid:12)(cid:12)(cid:12)(cid:12) H (cid:20) AH + β A s (cid:21) ∂ S ∂ x (cid:33) (21) edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 11 S u r f a c e e l e v a t i on ( m ) S(x) for (x) = 0S(x) for (x) = 0.5z b (x) Fig. 4:
Steady state glacier surfaces for constant basal slip with bumpy inclined bedrock, z b = b .which has two unknowns H and β . Since the free surface velocity is also given, rearrangingthe equation for β gives β = − A s u s ( ρ g ) (cid:12)(cid:12)(cid:12) ∂ S ∂ x (cid:12)(cid:12)(cid:12) ∂ S ∂ x H + AH (22)Substituting this expression for β into the steady state surface equation above and integratingresults in 0 = (cid:90) x adx − ( ρ g ) (cid:12)(cid:12)(cid:12)(cid:12) ∂ S ∂ x (cid:12)(cid:12)(cid:12)(cid:12) ∂ S ∂ x AH − u s H + C (23)which is a polynomial equation in only H with constant of integration C = − (cid:90) x adx + ( ρ g ) (cid:34)(cid:12)(cid:12)(cid:12)(cid:12) ∂ S ∂ x (cid:12)(cid:12)(cid:12)(cid:12) ∂ S ∂ x (cid:35)(cid:12)(cid:12)(cid:12) (cid:12) (cid:12) x = x AH + ( u s ) H (24)where x is some point inside the domain of the glacier where the height is known. It isreasonable to assume height can be known at one location from practical measurements.There are numerous methods which could be employed to solve the polynomial for H .Newton’s method is chosen for it’s simplicity and controllability. Hence solving (23) for c George et al. S u r f a c e e l e v a t i on ( m ) S(x) for (x) bumpS(x) for (x) stepz b (x) Fig. 5:
Steady state glacier surfaces for basal slip with a bump and basal slip with a stepwhere z b = f .each H i = H ( x i ) using Newton’s method H n + i = H ni − F ( H ni ) F (cid:48) ( H ni ) (25)with the following functions F ( H i ) = (cid:90) x i adx − ( ρ g ) (cid:18) ∂ S ∂ x (cid:19) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x i AH i − ( u s ) i H i + C (26) F (cid:48) ( H i ) = − ( ρ g ) (cid:18) ∂ S ∂ x (cid:19) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x = x i AH i − ( u s ) i . (27)For Newton’s method to find the correct root of F it is important to start with a nearby guess.Therefore, the method will move away from x to the left and right using H i = H for x i = x (28) H i = H i − for x i > x (29) H i = H i + for x i < x (30)For each glacier, x is chosen to be in the middle of the glacier domain. edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 13 S u r f a c e e l e v a t i on ( m ) S(x) for (x) bumpS(x) for (x) stepz b (x) Fig. 6:
Steady state glacier surfaces for basal slip with a bump and basal slip with a stepwhere z b = b . Figures 7 to 10 show the inversion results for each combination of underlying bedrock andbasal slip distribution. The recovered bedrock profile elevation, z b , rec , is shown in each sub-fig. (a). The recovered basal slip distribution, β rec , is shown in in each subfig. (b). Therecovered variables in each case are overlaid on the true variables for easy comparison.Table 3 gives relative errors for the recovered bedrock and basal slip distribution re-spectively. The relative error between the recovered variables, x rec , and the true value, x , iscalculated by E ( x rec , x ) = || x rec − x || || x || = (cid:112) ∑ i (( x rec ) i − x i ) (cid:113) ∑ i x i , (31)where i runs along the glacier domain. Note that (31) is not defined for x =
0. In this case,the error is defined as E ( x rec , x ) = || x rec || = (cid:114) ∑ i (( x rec ) i ) . (32)For each benchmark case, the bedrock reconstruction is in very good agreement withthe true profile. The relative errors in bedrock recovery, as given in table 3, are of magnitude c George et al. − for each case. This indicates that the has high accuracy for recovering bedrock eleva-tion profiles. This is illustrated in each sub-fig. (a) which show close alignment between therecovered bedrock and the true bedrock.For each benchmark case, the recovered basal slip distribution is in agreement with itstrue distribution. The relative errors, as given in table 3 are of magnitude 10 − or smallerfor each case. While this is a much larger error than for bedrock recovery, it is still of smallmagnitude. In the (b) sub-figs., the recovered variable closely aligns with the true values. Itis visible that areas or largest error are at the top and tail ends of the glacier. E l e v a t i on ( m ) z b,rec (x) for (x) = 0z b,rec (x) for (x) = 0.5z b (x) true (a) B a s a l s li p d i s t r bu t i on , ( x ) rec (x) for (x) = 0 rec (x) for (x) = 0.5(x) true (b) Fig. 7:
Recovered bedrock (a) for non-constant basal slip where z b = f and correspondingrecovered basal slip (b). E l e v a t i on ( m ) z b,rec (x) for (x) = 0z b,rec (x) for (x) = 0.5z b (x) true (a) B a s a l s li p d i s t r bu t i on , ( x ) rec (x) for (x) = 0 rec (x) for (x) = 0.5(x) true (b) Fig. 8:
Recovered bedrock (a) for constant basal slip where z b = f and corresponding re-covered basal slip (b). edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 15 E l e v a t i on ( m ) z b,rec (x) for (x) bumpz b,rec (x) for (x) stepz b (x) true (a) B a s a l s li p d i s t r bu t i on , ( x ) rec (x) for (x) bump rec (x) for (x) step(x) true (b) Fig. 9:
Recovered bedrock (a) for non-constant basal slip where z b = f and correspondingrecovered basal slip (b). E l e v a t i on ( m ) z b,rec (x) for (x) bumpz b,rec (x) for (x) stepz b (x) true (a) B a s a l s li p d i s t r bu t i on , ( x ) rec (x) for (x) bump rec (x) for (x) step(x) true (b) Fig. 10:
Recovered bedrock (a) for non-constant basal slip where z b = b and correspondingrecovered basal slip (b). Glacier surface data in reality always has some noise. Therefore, for practical applicationsthe method should be capable of handling noisy data. To evaluate the effect of noise on theinversion method, noise is added to each of the measured variables; S , u s . Noise is added tosynthetic surface data in the following way: S noise = S true + ε n max ( H true ) and, (33) u s , noise = u s , true + ε n ( max ( u s , true − min u s , true ) , (34)where n ∈ [ − , ] is randomly distributed. The amount of noise to be easily adjusted bychoosing ε ∈ [ , ] where the larger the choice, the more noise. Once noise is added to sur-face data, the result is smoothed with a local regression using weighted linear least squaresand a second degree polynomial model which assigns less weight to outliers in the regres-sion. The local span for the regression is 20 % of the data points. Data outside six mean c George et al.
Table 3:
Associated errors for recovered bedrock profile elevation and basal slip distributionas defined by Eqs. (31) and (32).
Case Relative errors β ( x ) z b E ( z b , rec , z b ) E ( β rec , β ) f b f b f b f b absolute deviations is given zero weight. Smoothing the data is important as ∂ S ∂ x appears reg-ularly in the equations and needs be well defined. Examples of noise added to surface dataas well as their smoothed counterparts can be found in the supplementary material.This process is applied to all benchmark cases with both ε = . ε = .
2. An ex-ample of the results for 100 samples of noisy surface data using constant slip and a bumpy bedrock with ε = . (a) (b) Fig. 11:
Recovered bedrock elevation profile (a) and basal slip distribution (b) for 100 sam-ples of noisy surface data where underlying bedrock is z b = b and β = .
5. The solutionenvelope is given in red, the median solution as blue circles and the true as a black line forfor reference. In this case, ε = . edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 17 For each combination of bedrock elevation profile and basal slip distribution, the recoveryof each input variable was good. Errors in the reconstructed bedrock elevation for all scenar-ios was negligibly small. Similarly for the basal slip distribution, though there were somerelatively larger errors in this recovery. The largest errors in all cases arose at either the topor bottom end of the glacier.At the top end of the glacier, there is a dome where the gradient of the free surfaceelevation is 0, in other words ∂ S ∂ x =
0. Clearly, given Eq. (11) for surface velocity, this resultsin a stagnation point in the free surface. Due to this, no information about β ( x ) or H ( x ) canbe recovered from Eq. (11). In Newton’s method, this stagnation point presents as F (cid:48) ( H ) =
0, which cannot be solved. Additionally, when F (cid:48) ( H ) becomes very small Newton’s methodbecomes unstable due to division by F (cid:48) ( H ) . To combat this, when | ∂ S ∂ x | ≈
0, the previoussolution for H is taken and the method skips to the next horizontal point.At the bottom end of the glacier, | ∂ S ∂ x | → ∞ . Because a finite number of points is used toapproximate the derivative in places where it changes rapidly, such as at the bottom end, theapproximation is worse. These approximation errors transfer across to the inverse solution.In addition, at the very top and very bottom where the ice ends, ∂ S ∂ x is discontinuous. Thisdiscontinuity gives rise to error also. The results show that it is possible to accurately recover both the bedrock elevation profileand basal slip distribution of a glacier for given surface elevation and velocity in certainrealistic synthetic cases. The method, although simple, is robust regardless of the underly-ing bedrock profile and basal slip distribution. This is a key result when considering theapplicability of the method to ‘real world’ problems in which bedrock and basal slip areunlikely to be uniform. A logical next step in developing the method is testing performancein a three-dimensional flow.Previous authors have focused on bedrock recovery in no-slip cases. That simplificationgives rise to many interesting methods but neglects basal slip which plays a vital role inglacial evolution. Indeed, since basal slip can drastically change the height of the glacierdue to dynamic thinning, the bedrock elevation recovery can have large error if this is notconsidered. The method presented here still returns the correct bedrock elevation profile andbasal slip distribution with the more general inclusion of basal slip for certain broadly realis-tic synthetic cases. Further, a previously unknown implication is that a unique combinationof bedrock elevation profile and basal slip distribution gives rise to unique surface elevationand velocity. Further studies into the generalised case are required to prove this. c George et al.
References
1. Adhikari S, J Marshall S (2012) Parameterization of lateral drag in flowline mod-els of glacier dynamics. Journal of Glaciology 58(212):1119–1132, DOI DOI:10.3189/2012JoG12J018, URL
2. Adhikari S, Marshall SJ (2011) Improvements to shear-deformational models of glacier dynamicsthrough a longitudinal stress factor. Journal of Glaciology 57(206), URL
3. Alley RB, Blankenship DD, Bentley CR, Rooney ST (1986) Deformation of till beneath ice stream B,West Antarctica. Nature 322(6074):57–59, DOI 10.1038/322057a0, URL { \ O } strem G (1975) Erts Data in GlaciologyAn Effort to Monitor Glacier Mass Balancefrom Satellite Imagery. Journal of Glaciology 15(73):403–415, DOI 10.3189/S0022143000034511,URL
5. Barcilon V, MacAyeal DR (1993) Steady flow of a viscous ice stream across a no-slip/free-sliptransition at the bed. Journal of Glaciology 39(131):167–185, DOI 10.3189/S0022143000015811,URL
6. Bevan SL, Luckman A, Khan SA, Murray T (2015) Seasonal dynamic thinning at Helheim Glacier.Earth and Planetary Science Letters 415:47–53, DOI 10.1016/J.EPSL.2015.01.031, URL
7. Bierman PR, Montgomery DR (2014) Key Concepts in Geomorphology. W.H.Freeman & Co Ltd8. Blatter H, Greve R, Abe-Ouchi A (2011) Present State and Prospects of Ice Sheet and GlacierModelling. Surveys in Geophysics 32(4-5):555–583, DOI 10.1007/s10712-011-9128-0, URL http://link.springer.com/10.1007/s10712-011-9128-0
9. Budd WF, Keage PL, Blundy NA (1979) Empirical Studies of Ice Sliding. Journal of Glaciology23(89):157–170, DOI 10.3189/S0022143000029804, URL
10. Cazenave A, Gregory JM, Jevrejeva S, Levermann A, Merrifield MA, Milne GA, Nerem RS, Nunn PD,Payne AJ, Pfeffer WT, Stammer D, Unnikrishnan AS (2013) Sea level change. In: Stocker T, Qin D,Plattner GK, Tignor M, Allen S, Boschung J, Nauels A, Xia Y, Bex V, Midgley P (eds) Climate Change2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Reportof the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UnitedKingdom and New York, NY, USA., chap 13, pp 1137–121611. Christoffersen P, Tulaczyk S (2003) Response of subglacial sediments to basal freeze-on 1. Theory andcomparison to observations from beneath the West Antarctic Ice Sheet. Journal of Geophysical Re-search: Solid Earth 108(B4), DOI 10.1029/2002JB001935, URL http://doi.wiley.com/10.1029/2002JB001935
12. Church JA, Clark PU, Cazenave A, Gregory JM, Jevrejeva S, Levermann A, Merrifield MA, Milne GA,Nerem S, Nunn PD, Payne AJ, Pfeffer WT, Stammer D, Unnikrishnan AS (2013) Sea Level Change. In:Stocker TF, Qin G, Plattner K, Tignor M, Allen S, Boschung J, Nauels A, Xia Y, Bex V, Midgley PM(eds) Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the FifthAssessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press,Cambridge, United Kingdom and New York, NY, USA., chap 13, pp 1137 – 121613. Cr¨uger T, Fischer H, von Storch H (2004) What do accumulation records of single ice cores inGreenland represent? Journal of Geophysical Research: Atmospheres 109(D21):n/a–n/a, DOI 10.1029/2004JD005014, URL http://doi.wiley.com/10.1029/2004JD005014
14. Cuffey KM, Paterson WSB (2010) The Physics of Glaciers. Elsevier Science & Technology Books15. Environmeantal Protection Agency (2017) Collecting Snow and Ice Data. URL
16. Farinotti D, Brinkerhoff DJ, Clarke GKC, F¨urst JJ, Frey H, Gantayat P, Gillet-Chaulet F, Girard C, HussM, Leclercq PW, Linsbauer A, Machguth H, Martin C, Maussion F, Morlighem M, Mosbeux C, Pandit A,Portmann A, Rabatel A, Ramsankaran R, Reerink TJ, Sanchez O, Stentoft PA, Singh Kumari S, van PeltWJJ, Anderson B, Benham T, Binder D, Dowdeswell JA, Fischer A, Helfricht K, Kutuzov S, LavrentievI, McNabb R, Gudmundsson GH, Li H, Andreassen LM (2017) How accurate are estimates of glacier icethickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment. The Cryosphere11(2):949–970, DOI 10.5194/tc-11-949-2017, URL edrock reconstruction from free surface data for unidirectional glacier flow with basal slip. 19
17. Flament T, R´emy F (2012) Dynamic thinning of Antarctic glaciers from along-track repeat radar al-timetry. Journal of Glaciology 58(211):830–840, DOI 10.3189/2012JoG11J118, URL
18. Fowler AC (1987) Sliding with Cavity Formation. Journal of Glaciology 33(115):255–267, DOI10.3189/S0022143000008820, URL
19. Gessese A, Heining C, Sellier M, Mc Nish R, Rack W (2015) Direct reconstruction of glacier bedrockfrom known free surface data using the one-dimensional shallow ice approximation. Geomorphology228:356–371, DOI 10.1016/J.GEOMORPH.2014.09.015, URL
20. Glen JW (1952) Experiments on the Deformation of Ice. Journal of Glaciology 2(12):111–114, DOI10.3189/S0022143000034067, URL
21. Greve R, Blatter H (2009) Dynamics of Ice Sheets and Glaciers. Advances in Geophysical andEnvironmental Mechanics and Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, DOI10.1007/978-3-642-03415-2, URL http://link.springer.com/10.1007/978-3-642-03415-2
22. Heining C, Sellier M (2016) Direct Reconstruction of Three-dimensional Glacier Bedrock and SurfaceElevation from Free Surface Velocity. AIMS Geosciences 2(1):45–63, DOI 10.3934/geosciences.2016.1.63, URL
23. Hubbard B, Glasser N (2005) Field Techniques in Glaciology and Glacial Geomorphology24. Hutter K (1981) The Effect of Longitudinal Strain on the Shear Stress of an Ice Sheet: In Defence ofUsing Stretched Coordinates. Journal of Glaciology 27(95):39–56, DOI 10.3189/S0022143000011217,URL
25. Iverson NR, Hanson B, Hooke RL, Jansson P (1995) Flow mechanism of glaciers on soft beds. Science(New York, NY) 267(5194):80–1, DOI 10.1126/science.267.5194.80, URL
26. Jiskoot H (2011) Dynamics of Glaciers. In: Singh VP, Singh P, Haritashya UK (eds) Encyclo-pedia of Snow, Ice and Glaciers, Springer Netherlands, Dordrecht, pp 245–256, DOI 10.1007/978-90-481-2642-2 \ https://doi.org/10.1007/978-90-481-2642-2_127
27. Kaser G, Fountain A, Jansson P (2003) A manual for monitoring the mass balance of mountain glaciers.Unesco, URL https://globalcryospherewatch.org/bestpractices/docs/UNESCO_manual_glaciers_2003.pdf
28. Kleman J, Glasser NF (2007) The subglacial thermal organisation (STO) of ice sheets. QuaternaryScience Reviews 26(5-6):585–597, DOI 10.1016/J.QUASCIREV.2006.12.010, URL
29. Lal D, Nishiizumi K, Arnold JR (1987) In situ cosmogenic 3 H, 14 C, and 10 Be for determining thenet accumulation and ablation rates of ice sheets. Journal of Geophysical Research 92(B6):4947, DOI10.1029/JB092iB06p04947, URL http://doi.wiley.com/10.1029/JB092iB06p04947
30. Le Meur E, Gagliardini O, Zwinger T, Ruokolainen J (2004) Glacier flow modelling: a comparison ofthe Shallow Ice Approximation and the full-Stokes solution. Comptes Rendus Physique 5(7):709–722,DOI 10.1016/J.CRHY.2004.10.001, URL
31. Lliboutry L (1968) General Theory of Subglacial Cavitation and Sliding of Temperate Glaciers. Jour-nal of Glaciology 7(49):21–58, DOI 10.3189/S0022143000020396, URL
32. Martin N, Monnier J (2015) Inverse rheometry and basal properties inference for pseudoplastic geophys-ical flows. European Journal of Mechanics - B/Fluids 50:110–126, DOI 10.1016/j.euromechflu.2014.11.011, URL https://linkinghub.elsevier.com/retrieve/pii/S0997754614001733
33. Monnier J, des Boscs PEE (2017) Inference of the Bottom Properties in Shallow Ice Approxima-tion Models. Inverse Problems 33(11), DOI 10.1088/1361-6420/aa7b92, URL https://doi.org/10.1088/1361-6420/aa7b92
34. Nye J (1953) The flow law of ice from measurements in glacier tunnels, laboratory experiments andthe Jungfraufirn borehole experiment. Proceedings of the Royal Society of London Series A Mathe-matical and Physical Sciences 219(1139):477–489, DOI 10.1098/rspa.1953.0161, URL
35. Pritchard HD, Arthern RJ, Vaughan DG, Edwards LA (2009) Extensive dynamic thinning on the marginsof the Greenland and Antarctic ice sheets. Nature 461(7266):971–975, DOI 10.1038/nature08471, URL c George et al.36. Rowan AV, Quincey DJ, Gibson MJ, Glasser NF, Westoby MJ, Irvine-Fynn TD, Porter PR, Hambrey MJ(2018) The sustainability of water resources in High Mountain Asia in the context of recent and futureglacier change. Geological Society Special Publication 462(1):189–204, DOI 10.1144/SP462.1237. Schwikowski M, Schl¨appi M, Santiba˜nez P, Rivera A, Casassa G (2013) Net accumulation rates derivedfrom ice core stable isotope records of P´ıo XI glacier, Southern Patagonia Icefield. The Cryosphere7:1635–1644, DOI 10.5194/tc-7-1635-2013, URL
38. Shuman CA, Berthier E, Scambos TA (2011) 20012009 elevation and mass losses in the LarsenA and B embayments, Antarctic Peninsula. Journal of Glaciology 57(204):737–754, DOI 10.3189/002214311797409811, URL
39. Weertman J (1957) On the Sliding of Glaciers. Journal of Glaciology 3(21):33–38, DOI10.3189/S0022143000024709, URL
40. Wilchinsky A, Chugunov V (2001) Modelling ice flow in various Glacier zones. Journal of AppliedMathematics and Mechanics 65(3):479–493, DOI 10.1016/S0021-8928(01)00054-5, URL40. Wilchinsky A, Chugunov V (2001) Modelling ice flow in various Glacier zones. Journal of AppliedMathematics and Mechanics 65(3):479–493, DOI 10.1016/S0021-8928(01)00054-5, URL