Calculating balance velocities with a membrane stress correction

Abstract

Inverse methods, where surface data are ‘inverted’ in order to quantify basal properties of ice sheets, play a major role in initializations. The balance-velocity method is a unique linear initialization, in which accumulation, surface elevation and thickness data are used to calculate the velocities and basal conditions required to maintain the observed ice-sheet altimetric signal, resulting in an estimate of the basal sliding viscosity that is guaranteed to be non-negative. We examine the observation that balance velocities based on the shallow-ice approximation (SIA) are extremely dependent on grid size, showing that Antarctic balance velocities on a 1 km grid are excessively over-channelized. Incorporating the membrane-stress approximation into balance-velocity calculations and comparing them with a simplified analytical solution shows that numerical error monotonically decreases with grid resolution and over-channelization is eliminated for Newtonian and non-Newtonian rheology. In contrast, for the SIA reducing grid size below the membrane-stress coupling length fails to improve accuracy. However, since this approach is nonlinear, a unique viscosity solution is not guaranteed, and in practice ‘sliding viscosity’ estimates are noisy. This raises problems of the sensitivity of these estimates to data and model errors, which may mean using inverse or smoothing techniques in association with balance-velocity methods in many, if not all, practical applications.

List of Symbols

List of Symbols

1. Introduction

Accurate data describing the physical environment at the base of an ice sheet or glacier are fundamental in making predictions of sea-level changes arising from ice-sheet disequilibration. The prediction of ice-sheet behaviour is becoming increasingly important, since observed changes in continental ice masses may be indicators of the operation of instability mechanisms that have the potential to have a significant impact on global sea levels over the next few centuries (e.g. Schoof, 2007). Predictive models need prescriptions of the viscosity and sliding conditions at the base of the ice sheet that are, in practice, largely unobservable.

There are various approaches that can use surface or elevation data to invert for the basal properties given sufficient data, such as the use of control methods for parameter tuning (e.g. MacAyeal, 1992, 1993), the use of a Bayesian inference approach (Gudmundsson and Raymond, 2008; Raymond and Gudmundsson, 2009) and the inverse Robin formulation of the initialization problem (Arthern and Gudmundsson, 2010; Gillet-Chaulet and others, 2012), as well as spin-up approaches (Pollard and DeConto, 2012) and more direct approaches associated with the use of linearized equations (Arthern and Hindmarsh, 2003; Barrand and others, 2013). The decision of which approach to use is partially dependent on the amount of data and computational resources available. The balance-velocity approach, which represents depth-averaged velocities of a steady-state ice sheet that is in balance, is a suitable method in many cases, since it is computationally inexpensive and requires data that are fairly abundant and accurate. In particular, the balance-velocities approach is appropriate for initializing with altimetric data. Such initializations are attractive because they control the rate of change of mass, which is of primary concern for ice-sheet modellers; many current methods give a substantial nonzero elevation rate. To calculate the balance velocities at every point on an ice sheet three datasets are needed: a map of surface elevation, an accumulation rate (or surface mass balance) and the ice thickness at every point.

Hitherto, balance velocities have been calculated using the shallow-ice approximation (SIA), where flow is always parallel to the direction of the slope. Budd and Warner, (1996) presented a computer scheme for rapid calculation of balance-flux distributions and showed how it could be used for the whole of the grounded ice sheet in Antarctica, and the method was further developed by Fricker and others, (2000). This method has been used to calculate balance velocities for Antarctica (e.g. Bamber and others, 2000a; Huybrechts and others, 2000; Wu and Jezek, 2004; Bingham and others, 2007) and Greenland (e.g. Joughin and others, 1997; Bamber and others, 2000b). This approach is based on following flowlines, and an alternative grid-based approach is presented by Hindmarsh, (1997), Arthern and Hindmarsh, (2003) and Leysinger Vieli and others, (2011).

The SIA method is usually applied to coarse resolutions (grid resolution > 5 km) and the surface elevation is usually smoothed over a number of ice thicknesses, in order to remove topographic sinks and account heuristically for the effects of longitudinal stresses which are not included in the approximation and would otherwise be expected to smooth out localized variations in the ice surface. Using a linearized approach, Kamb and Echelmeyer, (1986) show that the effect of slope and thickness perturbations drops off exponentially, and suggest that this exponential decay length, termed the ‘longitudinal coupling length’, should be between four and ten ice thicknesses. An alternative scale estimate of the membrane coupling length (MCL) is given by Hindmarsh (2006a), which represents the distance over which velocity and stress fields are smoothed due to variations in, for example, the bed properties. This quantity is conceptually related to the ‘longitudinal coupling length’ of Kamb and Echelmeyer, (1986), but is quantitatively different, arguably easier to compute and fully reflects the nonlinear nature of the coupling. The coupling length estimate can be written in forms that do not include the ice thickness, and it tends to infinity for ice shelves, as one would expect. Up to now, the choice of coupling length is, in practice, determined a posteriori by inspecting the resulting flux distribution (Le Brocq and others, 2006), and the smoothing is often done over length scales of up to 20 ice thicknesses (e.g. Bamber and others, 2000a; Testut and others, 2003). Consequently, the SIA cannot reliably capture the fine detail of ice flow and this method has been shown to be unsuitable for use with high-resolution datasets (grid size < 5 km) (Le Brocq and others, 2006). Moreover, the SIA method is apparently ill-posed, as results are severely dependent on grid size (e.g. Testut and others, 2003; Le Brocq and others, 2006; Van den Berg and others, 2006) and the assumption of downhill flow used in the SIA is not always appropriate (e.g. Hindmarsh, 2004; Le Meur and others, 2004). Here we use the MCL to explore the spatial distribution of the stress-gradient coupling in the Antarctic ice sheet, in order to assess the effects of the choice of smoothing length scale when using the SIA. Membrane stresses also affect how far inland perturbations in the forcing at the front of an ice shelf can propagate. A recent result is that the propagation distance is rate-dependent, with even very rapid high-frequency forcings being propagated tens of kilometres inland directly via membrane stresses (Williams and others, 2012). In these ways, membrane stresses may act over distances greater than the smoothing length scales used in SIA balance-velocity calculations.

In this study horizontal membrane (or longitudinal) stresses, which are caused by extension and compression directed along a thin body, are explicitly introduced into calculations for balance velocities, in order to address these issues (see Hindmarsh, 2006a, for further details of membrane stresses). Incorporating such stresses could decrease or remove the need for surface smoothing and allow accurate representation of the high-resolution dynamics of ice flow, which is not possible using the SIA. The addition of these stresses also removes the assumption that the flow has to be downhill and is a step towards more realistic, full-Stokes modelling.

We firstly present a SIA balance-velocity approach based on the work of Hindmarsh, (1997) and use it to calculate balance velocities for the whole of Antarctica on a smoothed 1 km grid. On comparing the resulting velocities with the measured velocities of Rignot and others, (2011), over-channelization of flow is seen to be a clear problem. We then calculate the MCL for the whole Antarctic ice sheet, and find that the smoothing required to represent the membrane stresses is not uniform in space. We conclude that over-channelization induces structure at length scales less than the MCL, and that, in consequence, membrane stresses must be included in the mechanical description.

A balance-velocity approach that incorporates membrane stresses (the membrane-stress approximation (MSA)) is then presented and is validated and applied to two simplified cases. Firstly, two analytical solutions for an idealized geometry with linear rheology and sliding are presented, in order to verify the numerical method and quantify the difference between balance velocities calculated with the SIA and MSA. One of these is a smooth solution, which produces good agreement between the numerical method and the analytical solution (Section 5). The second solution is non-smooth, and induces convergence difficulties for the nonlinear equations and non-smooth solutions (Section 8). Secondly, the MSA and SIA balance-velocity approaches are applied to an example of ice flowing with a nonlinear rheology with idealized ice-stream geometry.

Attempts to apply the approach to the full extent of the Antarctic ice sheet were made, but were not successful. We discuss this, and conclude that the method is too sensitive to data errors to allow solutions or convergence of solutions in many practical cases. Consequently, at best it can estimate smoothed sliding viscosities, likely in an inverse framework. The paper ends with a discussion of the MSA balance-velocity approach and its potential in future initialization procedures.

2. Balance Velocities with the Shallow-Ice Approximation

First we present the balance-velocity approach using the SIA, following the method proposed by Hindmarsh, (1997), and in the following section we show how the approach can be altered to include membrane stresses. We write the force balance equation in the form

(1)

where U* (X*, Y*, t*) is the flow velocity vector, H* (X*, Y*, t*) is the ice thickness, s* (x*, y*, t*) is the ice surface, ρ ^* is the density of ice and g ^* is gravity. A full list of symbols is provided in the List of symbols above; asterisk denotes dimensional quantities and no asterisk denotes non-dimensional quantities throughout. β* ² = 0, the is a coefficient of basal sliding: when β*² = 0, the equations describe an ice shelf and when β*² is large the equations represent a strongly grounded ice sheet. M* is a membrane-stress operator, defined in detail in the next section. The SIA is obtained when M* ≡ 0, thus in this case

(2)

The model can be scaled using

(3)

where l ^* is the domain length and u ^* and v ^* are the two horizontal components of the velocity vector U ^* . The scaled equations (with dimensionless quantities shown without an asterisk) then become

(4)

In the vertically integrated approximation, the flux, q , is simply the product of velocity, U , and ice thickness, H, thus

and Eqn (4) can be written as

(5)

Rewriting Eqn (5) gives

(6)

where

(7)

is a dimensionless diffusivity that measures aspects of basal properties, specifically the basal sliding viscosity. Coupling the momentum balance equation with the continuity equation for steady state,

(8)

where a is the dimensionless accumulation rate, and discretizing using a finite-difference scheme (see Appendix) gives the matrix system

where I is the identity matrix, q = ( q_X , q_Y), and θ = ▽sThe matrices D _x , D _y indicate discretized first-differential operators defined at grid centres and we defineθ_x and θ_y as operator matrices for the x and y components of the slope, respectively, such that the slopes are computed downwind at the same locations as for the diffusivity vector, D . This linear system can be solved using standard methods, and gives a unique solution for D , provided there are no sinks in the surface. This is a slightly disguised version of the grid-based method for estimating balance fluxes given by Hindmarsh, (1997). Smoothing is often used to eliminate sinks. Further comments on the numerical implementation are provided in the Appendix.

The recent publication of the 1 km grid size dataset for Antarctica, Bedmap2 (Fretwell and others, 2013), allows high-resolution SIA balance velocities to be calculated, as shown in Figure 1b. In order to construct this figure, we solve the matrix system described above for the case of sliding to compute balance velocities for all the grounded ice of Antarctica. This method has previously been used to calculate balance velocities of Antarctica at 5 km resolution (Leysinger Vieli and others, 2011). The accumulation map of Arthern and others, (2006) interpolated onto the 1 km grid together with the ice surface and bed from Bedmap2 were then used to calculate the 1 km SIA balance velocities for the Antarctic ice sheet, as shown in Figure 1b. The ice shelves were masked out. In this case the surface slopes were smoothed to eliminate the numerous small sinks in the ice surface, whose existence contradicts one assumption of the SIA balance-velocity approach, but unsmoothed thicknesses are used in the balance-velocity calculation (as Leysinger Vieli and others, 2011). The surface smoothing was done by fitting a polynomial of order n _p to a square surface patch of n _p data cells on either side of each point, where n _p was chosen to be sufficiently large to ensure enough smoothing to eliminate most of the sinks. In this case we use n _p ¼ 20, implying a 20th-order polynomial fit to a 40 km 40 km square around each point. This left 82 sinks, but since these were spread out over the whole ice sheet and only accounted for 3: 24 10^-4% of the points, these were left as sinks, causing singularities to remain at these points. It is clear from Figure 1 that this did not have a significant effect on the overall results.

Fig. 1. (a) Antarctic ice velocity derived from satellite radar interferometry measurements of Rignot and others, (2011). (b) Balance velocities calculated using a SIA with 1 km resolution and a 40 km smoothing of surface slopes. The insert in both cases shows velocities for Pine Island Glacier, West Antarctica, the location of which is indicated by the black box in the main figure.

Figure 1a shows Antarctic velocities compiled from satellite radar interferometry datasets (Rignot and others, 2011), and the calculated SIA balance velocities are shown in Figure 1b; both clearly show the widespread complex flow found by Bamber and others (2000a). However, although both velocity maps have the same general features and patterns of ice streams, the ice streams calculated by balance velocities are evidently more concentrated and channelized than the real data suggest. This over-channelization is an expected feature with the SIA (e.g. Le Brocq and others, 2006), and indicates that the SIA cannot properly resolve the dynamics of fast flow at this spatial resolution. Additionally, this over-channelization effect occurs even with a relatively high amount of smoothing, indicating that the smoothing length scale did not have a significant effect on our results.

3. The Membrane Coupling Length for Antarctica

Hindmarsh (2006a,b) proposed a length scale over which the membrane-stress terms in the force balance equation are expected to act (i.e. when the membrane operator term, M^* in Eqn (1), is of comparable size to the driving stress and the drag terms):

(9)

where B ^* is the ice stiffness (or rigidity) given by B* = A*^-1/n and A ^* is the flow-law parameter, which is a function of ice temperature and fabric; " is the surface slope. We refer to this as the MCL, MCL. The MCL can be calculated at every gridpoint on the 1 km Bedmap2 grid of Antarctica (Fretwell and others, 2013), using velocities from Rignot and others, (2011), slopes calculated from the Bedmap2 surface and using surface temperatures interpolated from the 5 km grid provided by Le Brocq and others, (2010), originally from Comiso, (2000), to calculate A ^* (see Cuffey and Paterson, 2010, for details). Results of this calculation are shown in Figure 2. From this it is clear that the MCL varies significantly over the continent, with high values in the central part of East Antarctica (due to low surface slopes) and generally lower MCL around the edge of the continent. Furthermore, high MCL of many tens of kilometres also appears along some of the main ice streams in West Antarctica, such as Pine Island Glacier and the ice streams of the Siple Coast, indicating that for many fast-flowing streams the membrane stresses act over distances of at least a similar magnitude to the smoothing length scales used to account for these stresses in many SIA calculations (e.g. Fricker and others, (2000), use a horizontal smoothing length of 35 km). Typically, the grid sizes and data smoothing used in earlier SIA calculations have heuristically replicated the smoothing effect of the membrane stresses. Consequently, spatially uniform smoothing with the SIA does not properly account for the non-uniform effects of the membrane stresses in all regions and this, together with the fact that the SIA does not resolve the dynamics at high resolution and leads to over-channelization of streams, suggests that to account for this spatial variability, membrane stresses need to be explicitly modelled in balance-velocity calculations. In the following section, we explicitly include these stresses, by introducing the membrane stress operator, M ^*, into the calculation of balance velocities.

Fig. 2. The membrane coupling length, MCL^* (from Hindmarsh, 2006a), calculated on the 1 km grid of Antarctica from Bedmap2 (Fretwell and others, 2013).

4. The Balance-Velocity Approach Incorporating Membrane Stresses

To incorporate membrane stresses into calculations for balance velocities we follow the much-used model formulation of Muszynski and Birchfield, (1987) and MacAyeal, (1989) and consider the vertically integrated approximation to give the membrane-stress form of Eqn (1) with the membrane-stress operator, M^* , as

(10)

where is the depth-average viscosity given by

(11)

where e*is the strain-rate tensor and n is the Glen coefficient (Cuffey and Paterson, 2010). After scaling with the relationships given in Eqn (3) and the force balance (Eqn (1)) becomes

(12)

and the dimensionless membrane-stress operator, MMM (Eqn (10)), is then given

(13)

where is a measure of dimensionless viscosity given by

(14)

This system is often referred to as the shallow-shelf or membrane-stress approximation (MSA). In order to facilitate comparison with the expression for the SIA model, shown in Eqn (6), we can use Eqns (5) and (7), for ice flux and diffusivity, D, respectively, to write

(15)

which is the same as Eqn (6) when M = 0. Coupling the momentum balances with the continuity equation for steady state (Eqn (8)) gives the discretized matrix system in terms of velocities as

(16)

where BBB represents a matrix with diagonal elements β^-2 . Note that BBB and D are related quantities and the equations are now nonlinear, even for a linear rheology. The whole matrix represents the mass conservation equation and two momentum balance equations, and includes operators and values for both boundary and interior points (see Appendix).

5. Analytical Solution: Linear Viscosity And Sliding

In this section, we present an analytical solution that can be used as a simple test case for the numerical model. We assume the ice sheet is rectangular and the ice acts as a Newtonian viscous fluid (n = 1) under linear sliding:

(17)

where τ* is the stress tensor. Since n = 1, the depth-integrated viscosity, , is now constant and (Eqn (11)), thus = 1 (Eqn (13)) and is the dimensionless measure of viscosity. We assume in this case that the measure of basal sliding, β² , is a constant, β² = 1, and consider the case of constant ice thickness, H = 1. On using the operator M the scaled vertically integrated momentum equations with linear viscosity and sliding become

(18)

(19)

Consider a trial solution of the form

(20)

(21)

(22)

where û and are constants dependent on k_x and k_y and is the coefficient of slope. E and F are parameters in the surface function. Substituting these expressions back into the momentum equations (Eqns (18) and (19)) leads to expressions for û and :

(23)

(24)

The boundary value problem on a rectangle is given by

where

(25)

The boundary conditions are automatically satisfied on boundary A, and are satisfied on boundary B if we choose where . Similarly, the boundary conditions on C are also automatically satisfied. On boundary D, the strain rate, e_xy ¼ 0, is satisfied if and the expressions for û and can be used to deduce

(26)

Note that this expression for the normal stress arises as a consequence of the other boundary conditions and is not set independently. This is not physically realistic, but allows an exact analytical solution to be formulated.

5.1. Comparing the SIA and MSA using the exact analytical solution

The exact analytical solution can be used to verify the numerical method for the case of linear rheology and sliding. The test surface is shown in Figure 3 for an ice stream 100 km in length and 200 km wide. Note that k_x and k_y can be any integer multiples of =w ₁ and =w ₂, respectively. Using the principle of superposition, one could find a general solution to the problem by summing all solutions from m = 1 to m = α. However, in this case we are predominantly interested in a solution that satisfies both the equations and the boundary conditions, so we take a particular solution, k_x = 4π /w ₁ and k_y = π/w ₂. The accumulation is prescribed to give an ice sheet in balance,

(27)

Fig. 3. A trial dimensionless ice surface used for comparing analytical and numerical solutions, with parameter values shown in Table 1 and a grid size of 5 km.

Setting dimensionless β² 1 and H = 1 gives the diffusivity D = 1 everywhere (Eqn (7)). Comparing the velocities and diffusivities calculated using the MSA and the SIA with this value enables us to assess the accuracy of the inversion methods.

Figure 4 shows velocities and diffusivities calculated from the analytical solution, the SIA and the MSA for a resolution of 2.5 km. For the MSA, a high level of agreement between the analytical and numerical solutions was found for both longitudinal and transverse velocities and for the diffusivities. However, for the SIA there are noticeable differences between the analytical and numerical solutions. In particular, the y –velocities and the diffusivities are considerably overestimated in the ‘valley’ regions (in which there is zero slope in the x –direction).

Fig. 4. (a–c) x –velocities, (d–f) y –velocities and (g–i) diffusivities calculated using the analytical solution (a, d, g), the MSA (b, e, h) and the SIA (c, f, i) for the ice surface in Figure 3 and a grid resolution of 2.5 km. All variables are dimensionless and parameter values are shown in Table 1.

5.2. SIA and MSA error analysis

To assess the accuracy of the SIA and the MSA the percentage error between the analytical and numerical speed was calculated at each internal H gridpoint for both approximations (the boundary gridpoints could not be used because the SIA only returns velocities on the interior gridpoints). ‘Speed’ at each gridpoint is defined as the absolute velocity, so that both horizontal velocity components are accounted for . These errors in speed were then used to calculate the root-mean-square (rms) error of speed for a range of grid resolutions. The results for both the MSA and the SIA are shown in Figure 5.

Fig. 5. The root-mean-square of the percentage error for speed for the MSA (crosses) and the SIA (open circles). The MCL for linear rheology and sliding (Eqn (28)) is shown as a black dashed line.

It is clear that as the grid resolution increases (and grid size decreases), the rms error monotonically decreases for the MSA but decreases and subsequently increases for the SIA. For a grid resolution of 5 km the MSA error is >15 times smaller than the SIA error, and for even finer grids (resolution ≤2.5 km) it is up to two orders of magnitude smaller than the SIA and continues to diminish as resolution is increased.

Hindmarsh (2006a) also derived an expression for the membrane couping length, MCL^* , for a Newtonian fluid under linear sliding as

(28)

where u* is taken as the typical horizontal speed and η^* is the Newtonian viscosity (here if n = 1). For the exact solution for the case of linear rheology and sliding, we calculated the MCL as L = 7: 75 km (using ε = 0: 0067, u ^* = 100 m a^-1, η* = 8.99 × 10⁶ Pa a, ρ* = 917 kg m^-3). This is shown as a dashed line in Figure 5, where we note that the error of the SIA increases significantly close to the point at which the MCL is achieved by the grid spacing, while the error for the MSA continues to decrease. Thus, once the grid size is smaller than the coupling length, the SIA does not capture the dynamics whereas the MSA does. This dependency on grid size highlights the difficulty in using the SIA when high-resolution data are required.

Although the error for the SIA is much larger everywhere than the MSA error when the grid resolution is less than ~ 10 km, Figure 6 shows that the SIA error is not evenly distributed throughout the domain: the error is largest in the regions of zero slope in the x –direction. In these valley regions, super-concentration and very fast flow occur for the SIA but not for the MSA, indicating that it is in these areas that the membrane stresses appear to be important in resolving the ice flow and preventing such discrepancies. Note that the data are smooth and are not forcing the SIA solutions to create structure and wavelengths shorter than those found in the data; in other words, there is a fundamental problem with the SIA when used in balance-velocity calculations, arising from the prescription that velocity is aligned with surface slope.

Fig. 6. Percentage speed error for SIA (circles) and MSA (crosses) at y ¼ 0: 5 for a grid resolution of 5 km.

6. An Example of a Synthetic Ice Stream with Nonlinear Rheology

For the MSA method to be useful in calculating balance velocities and diffusivities for real ice sheets or ice-covered regions, the method will need to be efficient for more-realistic geometries and ice rheology. Thus, in this section, the method is tested for a synthetic ice surface with varying thickness, non-uniform diffusivity and nonlinear rheology (n = 3). An ice surface was generated by running an ice-sheet model (BASISM; Hindmarsh, 2009) forward to steady state with a prescribed accumulation and non-uniform diffusivity to create an ice stream. The forward model includes membrane stresses in the force balance and has a water-pressure boundary condition on the ice front and free slip on the other three sides. The resulting surface, diffusivities and steady-state speeds for a grid resolution of 2 km are shown in Figure 7. An accumulation rate of 0.3 m a^-1 was used and the basal elevation was prescribed as zero everywhere (thus s* = H ^*). These surface and accumulation profiles were then used in the balance-velocity inversion models for both the MSA and the SIA, with the same boundary conditions.

Fig. 7. (a) The prescribed synthetic ice surface generated by the forward model, (b) the prescribed diffusivity and (c) the steady-state ice speeds from the forward model. The grid resolution is 2 km.

The performance of the balance-velocity inversion methods can be calculated by comparing the balance speeds and retrieved diffusivities with the speeds from the steady-state forward run and the prescribed diffusivity, as shown in Figures 8 and 9, respectively. Figure 8b and d show that the SIA method is prone to over-channelization of fast flow in the ice-stream region when compared with the steady-state speeds from the forward model (Fig. 7c), the same effect as shown for the SIA velocities for the simplified analytical solution. Conversely, the same comparison for balance speeds calculated from the MSA (Fig. 8a and c) does not show these effects, and speeds are a better match to the steady-state speeds. In fact, the rms percentage speed error on comparing the speeds from both methods with the steady-state speed is over five times higher for the SIA than the MSA (7.86% and 1.38%, respectively). In terms of the retrieved diffusivity, Figure 9 shows that the diffusivity mismatch is also smaller for the MSA than the SIA, although the retrieved diffusivities are not as accurate as the retrieved speeds for either method (the rms errors are 6.17% and 8.92% for the MSA and SIA methods, respectively). Some edge effects are apparent at the ice divide and these contribute significantly to the high rms, but since the diffusivity is small here anyway, implying very little sliding, these errors do not significantly affect the ice velocities (Fig. 8). Overall, as for the analytical case, the MSA method provides a closer match than the SIA to both the prescribed diffusivities and speeds, highlighting the need to include these stresses in calculations for balance velocities.

Fig. 8. (a, b) Balance speeds for the surface shown in Figure 7a calculated using (a) the MSA inversion method and (b) the SIA inversion method. (c, d) The percentage error of the speeds for the MSA (c) and the SIA (d) methods, when compared with the steady-state speeds for the forward model (Fig. 7c).

Fig. 9. (a, b) Diffusivities for the surface shown in Figure 7a calculated using (a) the MSA inversion method and (b) the SIA inversion method. (c, d) Percentage diffusivity errors for the MSA (c) and the SIA (d) methods. The prescribed diffusivity is shown in Figure 7b.

7. MSA Balance Velocities in Real-World Cases

Our next step was to use real surface, bed and accumulation data from Greenland and Antarctica, with cases centred on ice divides (e.g. around Summit, Greenland) and lakes (e.g. around Ellsworth Subglacial Lake, Antarctica). Numerous attempts to use this direct method to retrieve diffusivities in real-world applications led to the system failing to converge. Since it is a nonlinear system, there is an issue of whether solutions exist, or, more likely, whether solutions exist with positive², which is a requirement for the second law of thermodynamics to be satisfied. We suppose that these problems are due to a combination of data errors and the fact that the solutions are prone to being noisy, which hampers or prevents convergence of the nonlinear equations. This prevents us from reproducing Figure 1b using the MSA rather than the SIA method for calculating balance velocities.

This is a subtly different issue from the usual non-uniqueness associated with data inversions, as it is related to the existence of non-negative solutions for the diffusivity with the nonlinear equation set. It might be profitable to investigate these further from a more mathematical point of view, but in terms of practical geophysics the use of balance velocities in an inverse framework (e.g. control methods (MacAyeal, 1992) or inverse Robin methods (Arthern and Gudmundsson, 2010)) may be a more practicable approach. In order to understand this further, we looked at idealized cases with non-smooth data as shown in the following section.

8. Linear Rheology with Steep Topography: The Effects of Non-Smooth Data

In order for the MSA balance-velocity method to be practical for the rugged topography of Antarctica, we need to investigate its robustness to steep topography and other rapid variations in the data. We create such a configuration using the analytical solution with linear rheology (n = 1) by summing exact solutions with different wavelengths, using the principle of superposition. This allows us to create a steep-sided trough, as shown in Figure 10a, for which surface and accumulation profiles with m = 1, 3,….21 21 are summed as follows:

(29)

Fig. 10. (a) The trial ice surface generated by summing solutions of different wavelength, k_xm (Eqn (29)), (b) the retrieved speeds and (c) diffusivities obtained using this surface and corresponding accumulation with the MSA inversion method. (d) The retrieved diffusivity across the flow at y ¼ 0: 5. The grid resolution is 5 km and the ice-stream width and length are 200 and 100 km, respectively.

The x – and y –velocities are obtained by summing the solutions for each k_xm , using this expression; uniform diffusivity, D = 1, uniform thickness, H = 1 and Glen index n = 1 were used. In this case, k_y = π/ 2w ₂, E = 25 and all other parameters are as shown in Table 1.

Table 1. Parameter descriptions and typical values for parameters used in the analytical solution (in non-dimensional units). The dimensional ice-stream width and length are shown in parentheses

The speeds and diffusivities retrieved using these surface and accumulation profiles with the MSA inversion method are shown in Figure 10. There are sharp gradients in the ice surface around the edges of the surface trough, and in these regions the retrieved speeds and diffusivities do not accurately match the analytical solution. In fact, these sharp changes cause oscillations in the retrieved diffusivity, D, in the transverse flow direction, as shown in Figure 10d. These problems may be indicative of issues that arise when attempting to calculate balance velocities for real datasets, as discussed in the preceding section. These problems are not unexpected, and might be due to numerical resolution of sharp gradients or more fundamental issues with the use of MSA in regions where there are sharp gradients.

9. Discussion and Conclusions

In this study we have described a method which uses a membrane-stress model and surface and accumulation data to invert for the basal diffusivity (a combination of sliding and rheological effects near the base). This is accomplished by incorporating membrane stresses into calculations for balance velocities, and in formulating this method we have explicitly shown how the omission of these membrane stresses can lead to unreliable results.

A numerical method for a vertically integrated model (based on the model described by Hindmarsh, 1997), which uses accumulation, elevation and thickness data to invert for horizontal velocities and diffusivities, has been described. To show the importance of membrane stresses for such calculations, we first calculated balance velocities using only the SIA for the whole of the grounded part of the Antarctic ice sheet, using the new Bedmap2 dataset (Fretwell and others, 2013). Comparison with measured velocity data (Rignot and others, 2011) demonstrated that at high resolution there is substantial over-channelization of ice flow in ice streams. Furthermore, we found that the length scale over which the membrane terms are expected to act (the MCL) varies widely over Antarctica, indicating that the surface smoothing that is frequently used with the SIA in order to account for the effects of membrane stresses (e.g. Bamber and others, 2000a; Le Brocq and others, 2006) should not be done uniformly over the domain. In addition, the MCL was found to be many tens of kilometres in some areas, especially in the fast-flowing ice streams, which is larger than the smoothing length often used. We note that in transient cases Williams and others, (2012) found that frontal effects can be passed many tens of kilometres upstream via the transmission of membrane stresses over distances considerably greater than the MCL. Thus we conclude that additional terms in the force balance, namely membrane stresses, are needed for inversions via the balance-velocity method.

The numerical method which includes the membrane-stress terms (MSA) has been verified using a simplified, exact analytical solution for the case of Newtonian viscosity and linear sliding. The SIA solution was computed for the same geometry and parameters, allowing a direct comparison of the two balance-velocity inversion methods. For the MSA presented in this paper, excellent agreement was found between numerical and analytical velocities and diffusivities. The average error in speed decreased monotonically as a function of increasing grid resolution, indicating no dependence of the solutions on grid spacing. However, for the SIA the error mismatch between the solutions was up to two orders of magnitude higher than the MSA error for a grid size of 2.5 km or less and, moreover, results were dependent on grid size, with the average error increasing once the resolution was greater than ˜10–12 km.

The length scale above which the SIA error begins to increase approximately matches the MCL. A key assumption in the SIA is that ice always flows downhill, causing ice to speed up too much and over-concentrate in narrow streams at the bottom of valleys. If the grid resolution is less than the MCL, the SIA is unable to capture the smoothing dynamics of membrane-stress coupling; instead the ice flow is sensitive to the smoothing provided by the numerical discretization. The incorporation of membrane stresses smooths ice flow to form broader ice streams that are not grid-size dependent and which accurately match the exact solution. Thus, whereas other studies use uniform computational smoothing techniques over the whole domain to account for the effects of membrane stresses (e.g. Bamber and others, 2000a; Testut and others, 2003), we have negated the need for such methods by resolving the system in a more physically realistic manner by explicitly including membrane stresses.

Balance velocities and diffusivities have also been calculated using both the SIA and MSA methods for an idealized ice surface generated from a forward model for the case of nonlinear rheology and non-uniform thickness and diffusivity. Diffusivity was prescribed in order to create an ice stream, and comparing the retrieved speeds and diffusivities with the prescribed or known equivalents showed, as for the simpler case, that the SIA method is less accurate than the MSA and results in over-channelized flow in ice streams. However, the retrieved diffusivities were not as accurate as the balance velocities (MSA error of 6.17% for D and 1.38% for speed). Since the MSA method does provide accurate speeds, it may be that these errors in diffusivity do not have a large effect on the ice flow. This appears to be the case for the area close to the ice divide, where the diffusivity is ˜50% wrong in places but the ice velocity is hardly affected. In any case, the MSA method is still considerably more accurate than the SIA method when inverting for basal properties in this instance.

We have shown that the MSA balance-velocity approach can be used for inverting surface data to give balance velocities and a measure of the basal sliding viscosity of an ice sheet in some cases. This method is relatively straightforward and retains the essential physics, possibly making it a suitable alternative to initializations by the SIA and spin-up techniques. The next stage is to apply this method to a real-world ice sheet, such as Antarctica. However, initial explorations indicate that while this method works well for idealized synthetic data, it appears to be sensitive to errors in the datasets and to large variations in surface topography, which we found can lead to oscillations in the retrieved diffusivity. These errors can lead to the method returning negative values of the diffusivity, which is clearly non-physical, and so, in practice, Antarctic or Greenland balance velocities cannot be calculated using this method. This implies that a smoothing or initialization with additional data, such as the horizontal velocities of Rignot and others, (2011) (Fig. 1a), is required, which implies relaxing the constraint on the flux divergence. This is a route we are investigating.

The issue of whether the full-Stokes equations should be used in initializations arises. Note that while the SIA produces non-smooth solutions from smooth data, such that the assumptions underlying the SIA are violated, this is not the case for the MSA. However, with non-smooth data, it seems problematic to achieve solutions with the MSA, and either the data should be smoothed, or a smoothing/inverse technique in the estimation of the basal viscosity needs to be used.

To conclude, incorporating membrane stresses into balance-velocity calculations in theory represents a significant and necessary improvement when using the SIA, and results for idealized cases are promising. To use this method as an inverse procedure, the numerical method needs to be made more robust, possibly by constraining with additional data, before being employed to initialize models of the major ice sheets of Antarctica and Greenland which can then be run forward in time to make predictions on the effects of future climate and oceanic changes on ice-sheet volume.