Consistent approximations and boundary conditions for ice-sheet dynamics from a principle of least action

Abstract

The formulation of a physical problem in terms of a variational (or action) principle conveys significant advantages for the analytical formulation and numerical solution of that problem. One such problem is ice-sheet dynamics as described by non-Newtonian Stokes flow, for which the variational principle can be interpreted as stating that a measure of heat dissipation, due to internal deformation and boundary friction, plus the rate of loss of total potential energy is minimized under the constraint of incompressible flow. By carrying out low-aspect-ratio approximations to the Stokes flow problem within this variational principle, we obtain approximate dynamical equations and boundary conditions that are internally consistent and preserve the analytical structure of the full Stokes system. This also allows us to define an action principle for the popular first-order or ‘Blatter–Pattyn’ shallow-ice approximation that is distinct from the action principle for the Stokes problem yet preserves its most important properties and elucidates various details about this approximation. Further approximations within this new action functional yield the standard zero-order shallow-ice and shallow-shelf approximations, with their own action principles and boundary conditions. We emphasize the specification of boundary conditions, which are problematic to derive and implement consistently in approximate models but whose formulation is greatly simplified in a variational setting.

1. Introduction

Ice-sheet flow is typically modeled as an incompressible gravitationally forced Stokes flow, albeit for a non-Newtonian fluid with a power-law rheology. This is combined with an energy equation describing the evolution of internal temperature. Here we are concerned only with the dynamical part of the problem, the equations of incompressible, nonlinear Stokes flow, with the assumption that the internal temperature distribution is known. The Stokes problem for ice sheets has certain properties with important consequences for both the physical fidelity and numerical solvability of an ice-sheet model: namely, the system of incompressible Stokes equations is self-adjoint and is characterized by positive-definite dissipation (i.e. internal deformational heating). These properties can be traced to the existence of a variational principle (or an ‘action’ principle) that is physically meaningful for this problem: namely, a measure of ice-sheet dissipation plus the rate of loss of potential energy is minimized under the constraint of incompressible flow.

At present, the numerical solution of the full Stokes problem is not practical or routine for large-scale modeling of ice sheets. Various approximations are employed to produce simpler models that are easier and cheaper to solve. These include the ‘shallow-ice approximation’ (SIA; Hutter, 1983), the ‘shallow-shelf approximation’ (SSA; Morland, 1987; MacAyeal, 1989) and the ‘first-order’ approximation (Blatter, 1995; Pattyn, 2003; often referred to as a ‘higher-order’ approximation). Traditionally, these approximations are derived from the partial differential equations (PDEs) of the Stokes problem by means of a scaling analysis, and the resulting approximate PDEs are discretized using finite differences or finite elements. Unfortunately, there is no guarantee that the numerical model so obtained will preserve the above-mentioned mathematical and physical properties of the Stokes system. In particular, there is no guarantee that the associated system matrix will be symmetric or positive-definite, in which case numerical solution techniques will be less than optimal. Such problems may be avoided by deriving both approximations and discretizations within a variational framework, ensuring that the numerical model inherits the favorable properties of the underlying variational principle.

Variational or action principles are ubiquitous in physics. A famous example is Hamilton’s principle in classical mechanics and its many extensions to modern field theories. The existence of a variational principle is advantageous for both the analytical and numerical formulation of a problem for a number of reasons, including:

1. 1. A variational principle is a concise statement of a problem in terms of a single scalar quantity, the functional, that is unchanged in all coordinate systems. Thus, the dynamical equations obtained from the variation of the functional naturally incorporate metric terms associated with any given coordinate system (e.g. Hunke and Dukowicz, 2002). A coordinate system is not even required since an unstructured grid may be used.

2. 2. The functional involves lower-order partial derivatives that are easier to discretize than the higher-order derivatives in the corresponding PDEs.

3. 3. Given an action functional in discrete form, the variational process automatically yields symmetric, often positive-definite matrix systems that, as a result, are generally easier to solve efficiently (e.g. Newton–Krylov methods). In fact, the variational principle may often be cast as an optimization problem, the minimization of an action functional, in which case a number of efficient solution techniques are available (e.g. the Fletcher–Reeves algorithm).

4. 4. A variational principle provides a natural means of enforcing constraints by the method of Lagrange multipliers.

5. 5. Boundary conditions involving stresses are implicit in the variational principle, eliminating substantial practical difficulties in their discretization and implementation. Similarly, boundary conditions involving specified velocities may be directly incorporated into the variational functional so that all boundary conditions become included in the equations derived from the variational principle.

6. 6. Finally, it is advantageous to make approximations within the variational principle for problems that possess one, thereby preserving the internal consistency of the resulting approximations with respect to the original equations and boundary conditions. This method of obtaining approximations to physical problems has a long history. For example, Salmon (1983, 1985) has been a strong advocate of extracting approximations to ocean dynamics using Hamilton’s principle. In that case, the approximate equations preserve the conservation laws of the underlying exact equations provided the associated symmetries are not violated by the approximation.

In this paper, we focus primarily on the last two items above: deriving consistent ice-sheet dynamical approximations in a variational framework, and demonstrating how complex boundary conditions, commonly a source of practical difficulty in numerical ice-sheet models, can be generated and implemented simply and consistently. The PDEs obtained from the approximate variational principle (i.e. the Euler–Lagrange equations) and their discrete approximations preserve a positive-definite dissipation function, which is important for physical fidelity. The fact that the functional depends on just a single scalar function, an invariant of the strain-rate tensor, allows us to make approximations in only one place. As a result, the approximate equations are automatically internally consistent.

We stress the distinction between a variational principle and a variational formulation as used in the finite-element method, for example. A variational principle is obtained when a variational formulation is converted into a single scalar functional that has an extremum at the desired solution. Note that a variational formulation can be constructed for almost any problem, and this is what makes the finite-element method so powerful. The finite-element method is used to solve a wide variety of numerical problems. It is defined by means of a variational or ‘weak’ formulation: the problem is posed in integral form by integrating the residuals of PDEs and boundary conditions over a set of arbitrary test functions; the solution to the associated PDEs is determined if the integrals vanish for all such test functions. However, not all variational formulations can be converted to a variational principle, let alone a ‘natural’ variational principle. A natural variational principle is usually associated with significant physical or mathematical properties of the underlying problem, and its existence implies the advantages outlined above. The existence of a variational principle may be recognized in the application of the finite-element method, as noted below in connection with ice-sheet modeling, although it need not be. These distinctions, as they pertain to the finite-element method in particular, are lucidly explained by Zienkiewicz and Taylor (2000, ch. 3).

The existence of a variational principle for the Stokes problem (both linear and nonlinear) has been known for several decades (e.g. Bird, 1960; Johnson, 1960). Early applications in glaciology of the principle formulated by Johnson (1960) may be found in Fowler (1981) and Oakberg (1981). Variational principles in ice-sheet modeling involving the various approximate models, particularly in the finite-element context, are to be found in Colinge and Rappaz (1999), Glowinski and Rappaz (2003) and Rappaz and Reist (2005) for example, and most recently in the work of Schoof (2006, 2010). The main use of the variational principle in these works has been to study the mathematical issues of existence and uniqueness, and to obtain error estimates for the weak formulation of the numerical problem. The present work differs in that it focuses solely on the use of variational principles for ice-sheet modeling, and particularly on the derivation of ice-sheet approximations and boundary conditions by use of an action functional. Our premise is that the variational principle is fundamental and that it should form the basis for both the analytical and numerical study of ice sheets. We aim to strengthen the justification for the commonly used approximations and, ultimately, to improve their numerical formulation and implementation. Our aim is to derive approximate variational principles for ice-sheet modeling successively from known, more accurate variational principles starting from the Stokes variational principle. This again differs from earlier work where the approximate variational principle is deduced from existing Blatter–Pattyn equations rather than from a more fundamental variational principle.

We begin the discussion with a review of the basic Stokes flow model for ice-sheet dynamics (section 2). The variational principle associated with the Stokes flow model is then introduced, and we discuss its properties, particularly with respect to boundary conditions (section 3). This is followed by deriving several common ice-sheet approximations, their boundary conditions and associated variational principles (section 4). Finally, we summarize the results and present perspectives for future work (section 5), particularly as they apply to numerical modeling. However, the discussion of numerical implementation is deferred to a later publication.

2. Stokes Flow Dynamics

2.1. Governing equations

We use Cartesian tensor notation and, where appropriate, the summation convention on repeated indices. Let x_i = (x, y, z) represent the components of a Cartesian position vector, u_i = (u, v, w) the corresponding velocity components, g_i the components of the gravitational acceleration, and ρ the density, assumed constant. Also, Cartesian indices may sometimes be labeled with coordinate labels, i.e. i, j, k, … ∈ {x, y, z}. For convenience we assume horizontal orientation so that g_i = (0, 0, − g).

The Stokes flow equation expressing conservation of momentum is

(1)

where τ_ij is a symmetric deviatoric stress tensor and P = − σ_kk /3 is the mean compressive stress (or in the case of u_i = 0, the hydrostatic pressure). The deviatoric stress tensor and the pressure make up the total stress tensor, σ_ij given by

(2)

where δ_ij is the Kronecker delta. The flow is assumed to be incompressible,

(3)

The deviatoric stress tensor is expressed in terms of the strain rate by a non-Newtonian constitutive relation:

(4)

where is the strain-rate tensor, given by

(5)

The effective viscosity coefficient is a function of , the second invariant of the strain-rate tensor, given by

(6)

Definitions of the second invariant found in the literature can differ by a constant (e.g. the left-hand side of Equation (6) often includes a factor of 2). In glaciology the constitutive relation is generally given by Glen’s flow law (Paterson, 1994), for which the effective viscosity is expressed as

(7)

a highly nonlinear, positive-definite scalar function of the strain rate and temperature, θ. Here μ ₀(θ) is a temperature-dependent prefactor where the temperature is provided by the solution of a separate energy equation. For present purposes the temperature is assumed to be a known function of position, θ = θ (x, y, z), independent of time. In ice-sheet studies, the exponent n is usually taken equal to 3 (for reference, a Newtonian fluid corresponds to n = 1). It is also assumed that the location of the upper and lower surfaces, z(s) = z _s(x, y) and z(b) = z _b(x, y), is specified. Thus, the Stokes system given by Equations (1–7) together with appropriate boundary conditions forms a closed, time-independent system for determining the velocities u_i and the pressure P.

2.2. Boundary conditions

A set of boundary conditions must be specified to uniquely determine velocities and pressure. The ice-sheet configuration may be entirely general and may include a calving front, for example. However, there is a great variety of possible boundary conditions and ice-sheet configurations. We do not attempt to be exhaustive. Instead, we discuss a representative subset of boundary conditions and surface boundary configurations. The treatment of other cases not specifically described will hopefully be obvious from those that are considered. In order to simplify the presentation, we assume that the ice sheet is isolated and bounded by two surfaces: an upper surface that is in contact with the atmosphere or partially submerged and whose outward unit normal vector has an upward-pointing vertical component , and a basal (bottom) surface that is in contact with the bed (or possibly submerged in part and not in contact with the bed) and whose outward unit normal vector has a downward-pointing vertical component . Following standard practice, we also assume that the slope of the two bounding surfaces is small,

although this assumption is not necessary for much of the development. The relevant outward unit vectors are defined in Appendix A.

The part of the upper surface that is exposed only to the atmosphere is taken to be stress-free, assuming that any applied surface stresses (e.g. air pressure and wind stress) are negligible. Thus, the appropriate upper surface boundary condition is given by

(8)

where n_j is the outward unit vector normal to the surface (denoted in Appendix A by ). Basal boundary conditions are both more complicated and not as well established. The simplest case occurs over that part of the basal surface where the ice sheet is grounded and frozen to a solid bed, in which case a no-slip condition applies:

(9)

These boundary conditions apply simultaneously but on different parts of the boundary. The simple but relatively common situation when Equation (8) applies to the entire upper surface and Equation (9) applies to the entire bottom surface completely determines the Stokes flow problem, for example. In general, the specification of three independent conditions at each point of the boundary, such as the three vector components of the velocity in Equation (9), is required to determine the problem. This becomes clear when considering the variational formulation of the problem in section 3.

More general boundary conditions are required locally if the basal surface is partially submerged in water, or if sliding occurs, as when the ice–bed interface is thawed (or deforming plastically). Here we consider only two out of many possible situations to illustrate the additional complications that may arise. In preparation, let us note that any vector may be decomposed into normal and tangential components, , such that , where the normal component is given by , and the tangential vector is given by .

The first case occurs when part of the ice-sheet upper surface is immersed in water. A related boundary condition concerning a calving front in the two-dimensional (2-D) depth-integrated shallow-shelf approximation is discussed separately in section 4.3.2. The water exerts a pressure of magnitude p _w normal to the surface, directed inward, while stress-free conditions apply in the tangential direction. Thus, the normal stress force is given by (σ_ijn_j )^⊥ = (n_kσ_kjn_j )n_i = − p _w n_i , while the tangential force vanishes, (σ_iin_j )^∥ = σ_ijn_j − (n_kσ_kjn_j )n_i = 0. Combining these two conditions, we obtain

(10)

for an upper surface immersed in water, as opposed to Equation (8) for a free surface or a surface exposed to the atmosphere. Note that Equation (10) may also be applied to the floating part of the basal surface by using the appropriate basal unit normal vectors. To distinguish between these two possibilities we use the notation when referring to an immersed boundary at the surface, while for the basal boundary case we write , where and are the externally applied pressures.

In the second case, we consider a more general basal boundary condition: an ice sheet sliding in direct contact with a solid bed such that a linear frictional sliding law applies. There are many other possibilities for a frictional sliding law, including those that involve plasticity (e.g. Schoof, 2006, 2010). The basal boundary condition is then composed of two separate parts: a condition normal to the surface, and a condition in the tangential plane. The condition normal to the surface requires that the ice sheet not penetrate the bed:

(11)

and therefore . That is, velocity normal to the bed is zero and the basal velocity is tangential. This specifies one component of the boundary condition. The tangential condition specifies the surface frictional force, in this case given by a linear function of velocity,

(12)

where β is the coefficient of the basal drag law (note that β ≥ 0 as required by the Clausius–Duhem inequality). Since the vector lies in the tangential plane, it is sufficient to specify any two of its components, which together with Equation (11) will provide three independent conditions. In this particular case it is advantageous to use the two horizontal components since the lower surface of a grounded ice sheet is assumed to be near horizontal.

For clarity, let us write some of these conditions in full in rectangular Cartesian coordinates. The submerged upper surface boundary condition (Equation (10)) is then

(13)

(14)

(15)

where the components of the upper surface normal vector are defined in Appendix A. The stress-free upper surface boundary condition (Equation (8)) is obtained by simply setting the upper surface hydrostatic water pressure, equal to zero.

At the bed, the no-penetration condition (Equation (11)) may be explicitly written out as

(16)

The second part, involving Equation (12), is more complicated. Using Equation (2), the tangential friction force is expressible in terms of just the deviatoric stress:

(17)

where the magnitude of the normal component of the deviatoric stress force is given by τ_n = n_kτ_kjn_j . Note that the tangential components of the frictional sliding stress force do not depend on pressure, which is to be expected. As indicated previously, it is sufficient to specify only the horizontal components of the frictional boundary condition; the vertical component is determined by the fact that the force is tangential to the basal surface. Using Equation (17), the horizontal components of the friction force are therefore given by

(18)

(19)

where the normal stress per unit area at the base is and we have indicated that we are explicitly dealing with the tangential basal velocity vector by the use of superscript (b). Thus, the complete basal boundary condition is given by Equations (16), (18) and (19). For completeness, noting that and making use of Equation (17), the vertical component of the friction force is given by

(20)

For later use, in order to discuss more general specified or applied surface force boundary conditions, it is helpful to introduce a tensor ς_ij , which is given by

(21)

for the specific case of the linear basal frictional sliding force (Equation (12)), and by

(22)

for a submerged surface boundary condition given by Equation (10).

2.3. Energy conversion budget

It is useful at this point to review the energy conversion budget associated with the Stokes system because it has a direct connection to the variational principle introduced in section 3. We stipulate that there are no applied surface traction forces and no body forces other than gravity, consistent with our assumption of an isolated ice sheet. After contracting Equation (1) with u_i , integrating by parts over the entire volume and applying Gauss’ theorem, we obtain

(23)

where V is the integration volume containing the entire ice sheet, S is its surface area, dS_j = n_j dS is the directed element of surface area pointing outward, dS is the area of a surface element, and the frictional surface stress ς_ij has been defined previously. Here we use g_iu_i = −gw, where w is the vertical component of velocity. Note that the rate of change of total internal energy (or pressure work) has been eliminated from Equation (23) because the flow is incompressible. The first term on the right-hand side of Equation (23) is the rate of change of total gravitational potential energy, defined as

(24)

The second term is the dissipation rate, or rate of heat generation by internal deformation, defined as

(25)

Since the coefficient of viscosity is positive-definite as required by the Clausius–Duhem inequality , the frictional heating always dissipates energy, . This is an important physical property that should be preserved in all approximations and discretizations of the full ice-sheet Stokes system. The last term on the right-hand side is the frictional dissipation on the boundaries by surface stresses, defined as

(26)

where ς_ij is defined by Equation (21) in the present case, and it is positive-definite because β ≥ 0 as noted earlier. Thus, both forms of frictional dissipation are positive-definite. Indeed, we expect dissipation to be positive-definite in general.

The physical content of Equation (23) is simply that changes in the gravitational potential energy, , of an isolated ice sheet are converted to heat due to internal deformation and surface friction, , such that

(27)

When boundary conditions are given by Equations (8) and/or (9) and there is no frictional dissipation on the boundary, the potential energy change is balanced by internal dissipation alone . When other boundary conditions are involved, such as Equations (11) and (12), then and potential energy changes are balanced by the total dissipation as in Equation (27).

3. A Variational Principle for Non-Newtonian Stokes Flow Dynamics

The Stokes problem for ice-sheet dynamics is usually formulated and solved as a system of PDEs together with the specification of the rheology and the associated boundary conditions. As we show next, it is possible to formulate the underlying physics of this problem in terms of an action principle that also fully determines the problem. In fact, such a formulation can be considered as the more fundamental of the two.

3.1. The action principle and the Euler–Lagrange equations

Elements of a variational principle for the Stokes flow of an incompressible non-Newtonian fluid already exist in the literature (Bird, 1960; Johnson, 1960). These have been our starting point for constructing a variational principle suitable for ice-sheet modeling. Consider the following functional, which we call the ‘action’ in loose analogy to other problems in physics:

(28)

where, in general,

(29)

and which, in the case of Glen’s law rheology, is given by

(30)

(We thank C. Schoof for making us aware of this compact method of incorporating an arbitrary strain-rate-dependent viscosity into the functional. This method may also be found elsewhere (e.g. Bird, 1960; Glowinski and Rappaz, 2003)). Note that the variable s in Equation (29) corresponds to in Equation (7). Here Σ _j (u) is a specified velocity-dependent surface vector that incorporates surface stress boundary conditions in the functional. It is defined by the requirement that ∂Σ _j/∂u_i = ς_ij . The integrals comprising the functional are taken over the volume of the entire ice sheet and its boundary, as appropriate. The arguments on the left-hand side of Equation (28) indicate that the variation is to be taken with respect to the components of u_i and P. Note that the functional maps the velocity and pressure (P) fields into a single scalar quantity (i.e. the sum of the integrals on the right-hand side).

The action (Equation (28)) makes it clear that P is not a physical pressure in the interior of the domain but is actually a Lagrange multiplier that enforces the incompressibility constraint, ∂u_i/∂x_i = 0. However, to be consistent with its treatment in section 2, we set physical pressure boundary conditions so that at least on the boundaries P is equivalent to pressure. These boundary conditions are contained in Σ _j .

Equation (28) is analogous to the functional given by Bird (1960) except for the presence of the last term on the right-hand side, which incorporates more general boundary conditions. The great utility of Equation (28) derives from the fact that the entire dynamical problem discussed in section 2, including the various boundary conditions, is converted into a different problem, namely, that of finding an extremum of the action in terms of the velocity and pressure fields. In mathematical terms, the Stokes problem of section 2 is equivalent to the following statement for the extremum of the action:

(31)

for arbitrary variations of the velocity and pressure, δu_i and δP respectively, except that δu_i may be subject to certain restrictions at the boundaries (e.g. Dirichlet conditions). Strictly speaking, because of the incompressibility constraint, this involves finding a stationary point, i.e. a so-called ‘saddle point’, rather than an extremum. However, this would become an extremum problem if the space of admissible functions were limited to divergence-free velocity fields.

To show that Equations (28–31) are equivalent to the Stokes system of section 2, we take the first variation of the action as indicated in Equation (31), as follows:

(32)

The stress, σ_ij , is given in terms of the deviatoric stress and the pressure by Equation (2), and the deviatoric stress is given by Equation (4). As usual in variational manipulations, we have integrated by parts and applied Gauss’s theorem. For clarity, the derivation of Equation (32) is carried out in more detail in Appendix B. We conclude that a condition for the action to be stationary (i.e. in order to have for arbitrary variations δu_i , δP) is that the following Euler–Lagrange equations have to be satisfied:

(33)

(34)

These equations are identical to the Stokes system of equations. Thus, the Stokes PDEs are implicitly contained in the functional, as are many of the boundary conditions, as we show next.

3.2. Boundary conditions

The second requirement for the stationarity of Equation (31) is the condition,

(35)

which must apply along all bounding surfaces. This boundary condition is very general and includes all the boundary conditions discussed in section 2.2, as well as others. There are potentially five possibilities for satisfying Equation (35):

1. 1. δu_i is orthogonal to the vector, (σ_ij − ς_ij )n_j . This case may be eliminated from consideration because it is impossible to satisfy for an arbitrary vector, δu_i .

2. 2. δu_i = 0. This case corresponds to a specified velocity on the boundary (a Dirichlet-type boundary condition; δu_i = 0 because no variation in velocity is then possible).

3. 3. (σ_ij − ς_ij )n_j = 0. This case corresponds to a specified boundary stress force (a Neumann-type boundary condition). Neumann boundary conditions, implicit in the functional, are an example of a natural boundary condition associated with the action principle. The stress-free case σ_ijn_j = 0 when ς_ij = 0, is a particular example. (Courant (1943) defines natural boundary conditions as those that imply the vanishing of boundary terms in the first variation of the action functional for a ‘free problem’. A free problem is one in which the functions admissible in the functional are not subject to boundary constraints.)

4. 4. and (σ_ijn_j )^∥ − (ς_ijn_j )^∥ = 0. This is a mixed boundary condition in which a Dirichlet boundary condition in the normal direction is combined with a Neumann boundary condition in the tangential direction. This case is quite common and occurs, as noted, in the case of basal sliding (e.g. Equations (16), (18) and (19)).

5. 5. and (σ_ijn_j )^⊥ − (ς_ijn_j )^⊥ = 0. This is another mixed boundary condition in which a Dirichlet boundary condition in the tangential direction is combined with a Neumann boundary condition in the normal direction. Although an actual possibility, this condition does not appear to occur in practice.

Neumann- and Dirichlet-type boundary conditions behave quite differently in a variational formulation. Neumann boundary conditions are already implicit in the functional, either as natural conditions or as those specified by the surface integral in the action, and need not be explicitly specified outside the action principle. In this case, the boundary velocities are unknowns to be determined by the Euler–Lagrange equations. Thus, the system of equations resulting from applying the variational principle already incorporates all Neumann boundary conditions. This is a valuable property of the variational principle, particularly in the discrete case, because the alternative, in which boundary conditions must be specified and implemented independently in the partial differential equations, can be problematic.

Dirichlet boundary conditions, on the other hand, must be explicitly incorporated into the functional before the variation is performed (thus implying δu_i = 0). That is, only internal velocities are determined by the Euler–Lagrange equations (i.e. only those velocities not specified on the boundary). This is difficult to do analytically but much easier in the discrete case. That is, it is much easier to choose basis functions that already satisfy the Dirichlet boundary conditions, particularly local basis functions as in the finite-element method.

The same considerations apply to mixed boundary conditions. Thus, in case (4) the normal component of velocity must be specified on the boundary within the functional, whereas the tangential stress boundary conditions are already present. However, it may be more difficult to incorporate just the normal component of the velocity into the functional rather than the entire velocity vector. One way of doing this may be to use Equation (16) to eliminate the vertical velocity wherever it occurs at the basal boundary, and then to perform the variation only with respect to the horizontal velocities. However, such a procedure may be practical only in the discrete case, but, on the other hand, in the discrete case the basal slopes in Equation (16) may be ill-defined. Another possibility is to add the no-penetration boundary condition to the boundary integral in Equation (28) as a constraint with its own Lagrange multiplier. Using these techniques, the matrix equations obtained from a discrete variational principle will incorporate all boundary conditions and will automatically be symmetric.

We have yet to specify the surface stress vector Σ _j (u) that appears in Equation (28). For simplicity, we assume that the surface stress vector is composed of just the two contributions already considered, , where is the contribution from the immersed boundary and is the contribution from frictional sliding. For the immersed surface boundary condition (Equation (22)), we have

(36)

which implies that

(37)

and for the case of the linear sliding law (Equation (21)), we have

(38)

which gives

(39)

with similar expressions for other types of sliding laws. Note that Equations (37) and (39) may apply simultaneously but over different parts of the basal boundary, as mentioned previously, i.e. Equation (37) would apply over an immersed floating part of the basal boundary, while Equation (39) would apply over the part of the basal boundary that is sliding in contact with the bed.

3.3. Relation to dissipation

The Stokes action (Equation (28)) and the associated action principle (Equation (31)) are a statement regarding the relationship between the rate of heat dissipation by the ice sheet and the loss of total potential energy. Consider the case of Glen’s law and the linear sliding law (Equation (12)), with no external forcing. Taking account of Equations (24–26), and restricting the admissible velocities to those that satisfy the continuity equation , which are denoted by , we may write the action for the case of Glen’s law as

(40)

where . Therefore, the Stokes action principle may be interpreted as requiring the minimization of a linear combination of total internal dissipation, total frictional boundary dissipation, and the rate of change of total potential energy. Noting that are concave functionals of velocity, and is directly proportional to the mean vertical velocity, this means that the ice sheet slumps or flows under gravity such that potential energy is converted directly into heat (section 2.3). Also, the velocity field adjusts itself to minimize total dissipation, given by the two terms in brackets in Equation (40) plus the rate of decrease of total potential energy. These properties are also expected to hold for an arbitrary rheology and for arbitrary frictional sliding laws. In the case of arbitrary rheology, applying the mean value theorem to Equation (29) twice, we obtain

(41)

where

(42)

and where is a nondimensional constant for any realization of a Stokes problem. We conclude that with a positive proportionality constant that may be rheology- and problem-dependent. In special cases, need not be problem-dependent (i.e. for a Newtonian fluid, and for a Glen’s law rheology). Similar results may be demonstrated for an arbitrary frictional sliding law, which will be dissipative in general. This relationship between dissipation and total potential energy provides important insight into the physical content of the Stokes system. One way of preserving this content in numerical models is to ensure that these models are consistent with (or derived from) the variational principle (Equation (28)).

4. Approximations to Stokes Flow for Ice-Sheet Dynamics

The numerical solution of the full three-dimensional (3-D) nonlinear Stokes system can be quite difficult, and although such models exist (e.g. Zwinger and others, 2007; Pattyn, 2008) they are computationally expensive. However, several approximations are available that are sufficiently accurate for many applications and are much cheaper to solve. We next show that these approximations can be derived and thus formulated within a variational framework, similar to that introduced in section 3 for the Stokes system, suggesting the possibility of new and improved numerical solution methods as a result.

For ice-sheet modeling, the fundamental approximation to the Stokes system of equations depends on the lowaspect-ratio assumption, δ = d/L ≪ 1, where d and L are characteristic vertical and horizontal length scales respectively, and d is of the order of the ice-sheet thickness. We refer to an approximation as ‘first-order’ if it is accurate to (i.e. the error is ). The model introduced by Blatter (1995) is indeed first-order accurate as shown by Schoof and Hindmarsh (2010). Two other models of practical importance, the SIA and the SSA, start from the first-order approximation and are derived by introducing further approximations. Because the first-order approximation is primary, we discuss it first and in most detail.

Approximate ice-sheet models have traditionally been derived by applying scaling arguments to the full Stokes system of section 2.1 and the associated boundary conditions of section 2.2. Such scaling arguments are nontrivial. Schoof and Hindmarsh (2010) show that Blatter’s (1995) original scaling argument used in deriving the first-order approximation is incorrect, although the correct set of first-order equations was actually obtained. They provide the correct scaling argument, which turns out to be quite complex. In Appendix C we show that Blatter’s (1995) derivation of the first-order approximation requires two separate and independent assumptions: a low-aspect-ratio approximation to the strain-rate tensor, and the seemingly unjustified neglect of two components of the corresponding deviatoric stress tensor. However, the latter assumption turns out to be justifiable by the physical requirement of a positive-definite dissipation rate. Remarkably, only the first assumption is required if the low-aspect-ratio approximation is made within an action principle. The second of Blatter’s assumptions is a natural by-product of using the action principle and is obtained because positive-definite dissipation is implicit in the action principle. This illustrates the consistency of approximations when derived from an action principle.

Historically, the first-order approximation was derived in two steps. In the first, Blatter (1995) used a scaling argument based on δ = d/L ≪ 1 to approximate the strain-rate and deviatoric stress tensors, but otherwise retained the form of the Stokes system as given by Equations (1–7) in four independent unknowns: u, v, w, P. However, he recognized that the vertical velocity, w and pressure, P in this approximate system are expressible in terms of the horizontal velocities, u, v by a simple quadrature, thus simplifying the system. Nevertheless, the model was expressed in non-standard form and the proposed solution method was unduly complicated, involving eight unknown variables. Pattyn (2003) took the next logical step and gave the problem its current formulation, thus simplifying the system to a set of two equations for the horizontal velocities alone (e.g. equations (15) and (16) in Pattyn, 2003; equations (2.1a,b) in Schoof, 2010). The solution is completed when w and P are recovered once the system for the horizontal velocity is solved. We refer to the reduced system for the horizontal velocities, together with the related equations for the vertical velocity and pressure, as the ‘Blatter–Pattyn’ (BP) model. As noted earlier, the other approximations, SIA and SSA, follow from further approximations to the BP model. The hierarchical evolution of these approximations is illustrated in Figure 1.

Fig. 1. Evolution of ice-sheet model approximations from the Stokes action (with the respective action functionals indicated). Shaded ovals designate an intermediate model.

4.1. Low-aspect-ratio scaling

Let us assume the existence of typical horizontal and vertical velocity scales, U, W respectively, and corresponding velocity gradient length scales L, d, D that characterize ice-sheet dynamics:

(43)

The vertical length scale d is assumed to be of the order of the ice-sheet thickness, while the magnitude of the vertical length scale D depends greatly on whether sliding is slow or fast. Thus, D is very large when sliding is fast and vertical shear is small, and relatively small when sliding is slow or when bed friction is important, as when the ice sheet is frozen to the bed. These three length scales define two non-dimensional parameters,

(44)

which determine the domain of validity of the various approximations. Schoof and Hindmarsh (2010) use an alternative but closely related definition for λ. Their parameter is λ _SH = λ ^1/n , where λ is the parameter defined in Equation (44) and n is the power-law index in Glen’s law. Assuming that tumbling or slumping motion dominates ice-sheet flow, the continuity equation (Equation (3)) further implies that U/L ≈ W/d. Thus, the scaling parameter δ may alternatively be expressed in terms of velocity scales,

(45)

As noted previously, the various approximations to Stokes flow are essentially governed by approximations to the strain-rate invariant in the Stokes action. Based on the above scaling, Schoof and Hindmarsh (2010) show that the Stokes strain-rate invariant (Equation (6)) may be approximated by

(46)

This may be considered a precursor to the strain-rate invariant that will later characterize the Blatter–Pattyn approximation, so the subscript ‘pBP’ stands for ‘proto-Blatter–Pattyn’. The omission of the horizontal gradients of vertical velocity that are present in the Stokes strain-rate invariant is characteristic of the Blatter–Pattyn approximation. According to Equation (46), the Blatter–Pattyn strain-rate invariant, and therefore the Blatter–Pattyn approximation itself, is at most first-order accurate in the low-aspect-ratio approximation. It should be noted, however, that this approximation selects a particular preferred direction, in this case the vertical direction. This means that the associated models are no longer rotationally invariant.

The two additional approximations derived from the Blatter–Pattyn approximation, the SIA and the SSA, depend on the relative magnitude of the various terms in Equation (46). The SIA, valid for the slow sliding regime, is obtained when the shear terms dominate, λ ≫ 1, and the SSA, valid in the fast sliding regime, is obtained when the horizontal gradient terms dominate, λ ≪ 1. Schoof and Hindmarsh (2010) carried out an extensive scale analysis using the method of matched asymptotic expansions to characterize the various approximations. They identify two regimes in which the SIA is valid: (1) the parameter range 1 ≪ λ ≤ δ ⁻¹ in which vertical shearing is large and dominates the mass flux, and (2) the parameter range λ ≫ δ ⁻¹, in which vertical shearing is sufficiently large to dominate over sliding. The SSA is also divided into two regimes: (1) the parameter range δⁿ ≤ λ ≪ 1 in which horizontal stresses are significant compared to friction but friction is not necessarily weak, and (2) the parameter range λ ≪ δⁿ ≪ 1, in which horizontal stresses are largely dominant or friction is weak. The domain of validity of the various approximations given by this scaling is illustrated in Figure 2.

Fig. 2. The domain of validity of the Blatter–Pattyn approximation (BP: shaded and stippled, δ ≪ 1) and its subordinate shallow-ice (SIA: stippled, right-hand side, λ ≫ 1) and shallow-shelf (SSA: stippled, left-hand side, λ ≪ 1) approximations. The dashed lines delineate regimes 1 and 2 discussed in the text.

In addition, the specified frictional sliding stress vector in Equation (39) may be approximated,

(47)

where u _H indicates that only horizontal velocity components are involved. Also, as in this equation, and unless specifically indicated otherwise, an index in parentheses henceforth indicates a horizontal index, i.e. (i), (j), (k), … ∈ { x, y}.

4.2. The first-order or Blatter–Pattyn approximation

Using Equations (46) and (47) in Equation (28), we obtain a first-order accurate approximate Stokes action,

(48)

that supports a set of first-order accurate Euler equations and boundary conditions. However, it is possible and advantageous to reformulate this problem further following the precedent of Blatter (1995) and Pattyn (2003).

4.2.1. An action principle

Recall that parameter P is a Lagrange multiplier used to enforce compliance with the continuity equation (Equation (3)), thus ensuring that continuity is satisfied at the extremum. We are therefore justified in using the continuity equation to replace certain terms in Equation (48) to obtain the modified action,

(49)

where

(50)

is an effective strain-rate invariant obtained from Equation (46) by eliminating the vertical velocity by means of the continuity equation. We have also used Equation (37) to eliminate . The reason for the BP subscript notation will become clear shortly. Note that we have introduced a new Lagrange multiplier, P ^*, that again enforces compliance with the continuity equation at the extremum. Although the actions and are quite distinct, they yield the same Euler equations and boundary conditions and therefore they coincide at the extremum. However, it will be more convenient to start from rather than in deriving the Blatter–Pattyn action. In Equation (49) the surface integral is intended to apply over the entire surface of the ice sheet, i.e. both the upper and basal surfaces. Recall that dS_j = n_j dS, and note that we refrain from using subscripts or superscripts to designate the various parts of the upper or basal surface. In the following, however, for clarity we need to subdivide the ice-sheet surface into different parts.

It is convenient to rewrite Equation (49) as

(51)

where we decompose the surface integral into contributions from three different parts:

1. (a) The upper surface whose integrand involves , which may be either zero for the part of the surface that is above water and assumed stressfree, or nonzero for the part that is immersed. There are two basal surface contributions:

2. (b) A portion of the basal surface, where the no-penetration boundary condition applies and the integrand involves . Note that the contribution due to P ^* on this portion of the surface vanishes due to the no-penetration boundary condition (Equation (11)). This at least partly incorporates the no-penetration boundary condition into the functional according to the discussion in section 3.2.

3. (c) A portion of the basal surface, , whose integrand contains . This is the part of the basal surface that is submerged but not in contact with the bed, i.e. the ‘floating’ part of the basal surface.

Equation (51) may now be decomposed as:

(52)

where

(53)

(54)

We use the notation dS_z = n_z dS = dxdy for the projection of the surface element area dSj onto the horizontal plane. Similarly, dS _(i) = n _(i)dS represents the projection of the surface element area onto a vertically oriented plane. Thus, dS_x = n_x dS is the projection onto the y-z plane, and dS_y = n_y dS is the projection onto the x-z plane. Note that is the only part of that depends on vertical velocity. This decomposition is crucial for the derivation of the Blatter–Pattyn approximation and is sufficient to determine P ^*. Taking the variation of with respect to the vertical velocity, we obtain

(55)

The action principle applied to this component of the variation yields a vertical Euler equation,

(56)

which represents simple hydrostatic balance, and also two boundary conditions: an upper surface (z = z _s) boundary condition,

(57)

that specifies either an immersed or a stress-free surface boundary condition on different parts of the upper surface, and a pressure boundary condition on the floating part of the basal surface (z = z _b),

(58)

where we assume that has been specified independently. However, this presents an inconsistency. Note that Equation (56) is a first-order ODE that cannot support two boundary conditions. This is a price one pays for the lowaspect-ratio approximation, i.e. we have lost higher-order vertical pressure gradients, as is also the case in other forms of the hydrostatic approximation. Following standard Blatter–Pattyn practice, we choose to use Equation (57) as the boundary condition for the integration of Equation (56). Integrating Equation (56) from the upper surface down, we obtain P ^* as a function of the coordinates:

(59)

This implies that basal pressure is hydrostatic, i.e.

(60)

where H = (z _s − z _b) is the ice-sheet thickness. However, note that P ^* (x, y, z _b) need not be equal to on . For consistency, therefore, the applied basal pressure must be hydrostatic in the Blatter–Pattyn approximation, which is not unrealistic, or else it is necessary to eliminate the specification of pressure as a boundary condition on the floating part of the basal surface, which is probably more practical. Incidentally, substituting these results in Equation (54), we observe that at the extremum.

Using these results in Equation (53), we finally obtain

(61)

a component of the action that depends on horizontal velocities only, and which, therefore, we call the Blatter–Pattyn action (i.e. becomes ). This form of the action is remarkable because the only surface contribution arises from the part of the basal surface that involves sliding, . Note that it incorporates all boundary conditions discussed earlier except for any Dirichlet boundary conditions on horizontal velocities. However, a no-slip basal Dirichlet boundary condition may be obtained in the limit of a very large basal drag coefficient, β. As shown below in section 4.2.2, this functional yields the ‘standard’ Blatter–Pattyn model and the associated boundary conditions (Pattyn, 2003; Schoof, 2010).

It is worth noting that the nonlinear part of the functional (Equation (61)) is positive-definite. The resulting Euler–Lagrange equations are therefore positive-definite and self-adjoint, and the corresponding matrix in the discrete case is positive-definite and symmetric. The linear part of the functional specifies the source or right-hand-side terms. Unlike in the Stokes case, the action principle in this case does correspond to a minimization problem.

Essentially the same functional as Equation (61) has appeared in the ice-sheet modeling literature in the context of finite elements, albeit for much simplified Blatter–Pattyn models in most cases. Thus, Colinge and Rappaz (1999) and Glowinski and Rappaz (2003) consider a 2-D functional (x-z plane) in the absence of the surface term in Equation (61), i.e. corresponding to free slip at the bed, and bounded by horizontal upper and basal surfaces. Rappaz and Reist (2005) extended this functional to three dimensions and allowed for variable upper and basal surface elevations, but were still limited to a free slip boundary condition at the bed. Schoof (2010) incorporated much more general boundary conditions such that his functional in a generalized sense might be viewed as containing Equation (61) as a special case. In all these cases, however, the functional was derived from a weak formulation, given the partial differential equations and boundary conditions, rather than directly from the Stokes action as is done here.

The Blatter–Pattyn action (Equation (61)) is used to determine horizontal velocities. As mentioned earlier, the variation of the action (Equation (49)) with respect to P ^* implies that the continuity equation (Equation (3)) must be satisfied. Integrating the continuity equation vertically upwards from the bed, we obtain

(62)

where we use the non-penetration boundary condition (Equation (16)) together with Leibniz’s rule to eliminate the basal vertical velocity. This establishes the vertical velocity once the horizontal velocities are known, at least over that part of the domain that lies directly above the portion of the basal surface denoted by where the basal vertical velocity is given by the non-penetration boundary condition. However, one also needs a basal vertical velocity over the remaining, ‘floating’ part of the basal surface, denoted by ; this can only be provided by coupling to an underlying ocean model and is not easily available. These issues are avoided if the basal surface is assumed to be everywhere grounded, i.e. if , as it is in the standard Blatter–Pattyn model (Blatter, 1995; Pattyn, 2003; Schoof, 2010) and in the remainder of this paper.

The derivation and use of to determine horizontal velocities implies a complete decoupling of and , or equivalently, of the horizontal velocity problem from the hydrostatic pressure problem. However, these two components of the total action are potentially coupled by the no-penetration condition (Equation (16)) that links the vertical and horizontal velocity components at the basal boundary. Although we have already taken this boundary condition at least partly into account in Equation (51), it is not clear that this is sufficient. Such an omission could lead to an error in the basal frictional sliding boundary condition in the standard Blatter–Pattyn model. However, we show in section 4.2.3 below that any such error is actually second-order and therefore negligible.

4.2.2. Model equations and boundary conditions

Following the procedure in section 3, the Blatter–Pattyn action (Equation (61)) yields the Euler–Lagrange equations for the horizontal velocities,

(63)

where we define an effective Blatter–Pattyn stress ‘tensor’ as

(64)

and where is the viscosity (Equation (7)) defined in terms of the Blatter–Pattyn strain-rate invariant, Equation (50). Note that 〈τ _(i)j 〉_BP is the same as the quantity Σ _ij in Schoof (2010). Aside from Dirichlet conditions, the boundary conditions implied by the action are given by

1. (a) on the upper surface, S _s,

   (65)

   

   and

2. (b) on a frictionally sliding basal surface, ,

   (66)

   

   and it is implicitly assumed that are the horizontal components of a tangential velocity vector. Recall that we have assumed .

Writing out these equations in Cartesian coordinates, the Euler–Lagrange equations (Equation (63)) become

(67)

(68)

and the boundary condition on the upper surface, case (a) as given by Equation (65), becomes

(69)

(70)

The above equations are essentially identical to the corresponding equations found in Pattyn (2003) and Schoof (2010) for example, except for the inclusion of a hydrostatic pressure gradient on the right-hand side of the Euler equations due to a possible immersed surface boundary.

The basal frictional sliding condition, case (b) as given by Equation (66), becomes

(71)

(72)

This boundary condition appears in Schoof (2010) in connection with more general frictional sliding stress laws, and it corresponds to those found in Van der Veen and Whillans (1989) and Pattyn (2003) by replacing the components of the basal unit vector with their small-slope approximation (Appendix A, Equation (A5)).

4.2.3. Analysis of the basal boundary condition

It is instructive to compare the above Blatter–Pattyn boundary condition with the corresponding condition as given by Equations (18) and (19) in the Stokes system. However, it is more convenient to start from the first-order ‘fully coupled’ Blatter–Pattyn model whose action is and whose effective deviatoric stress tensor 〈τ_ij 〉 is given by Equation (C3) in Appendix C. We wish to compare the horizontal components of the tangential stress force associated with this tensor with those implied by the standard Blatter–Pattyn boundary condition (Equations (71) and (72)). These components are given by

(73)

where

(74)

Making use of the fact that and substituting into Equation (73) yields the correct tangential boundary condition:

(75)

(76)

Recalling our assumption from Appendix A that

(77)

Equations (75) and (76) coincide with the Blatter–Pattyn basal boundary condition (Equations (71) and (72)) to leading order. That is, the Blatter–Pattyn basal boundary condition represents a first-order accurate approximation of the Stokes tangential sliding boundary condition (Equations (18) and (19)) in the small basal slope limit.

In principle, it is possible for basal slopes to be relatively large even while the low-aspect-ratio assumption holds in a global sense. However, noting that the spatial scale definitions (Equation (43)) also hold locally implies that the low-aspect-ratio approximation, δ = d/L ≪ 1, will fail locally for short-wavelength, large-slope basal topography. Thus, the Blatter–Pattyn model requires both assumptions to be valid simultaneously.

4.3. Shallow-ice and shallow-shelf approximations

4.3.1. The shallow-ice approximation

In the slow sliding regime we expect the ice sheet to be frozen to the base or the basal drag to be very large so that vertical velocity gradients are dominant. In other words,

(78)

and hence we may approximate the Blatter–Pattyn strain-rate invariant (Equation (50)) by

(79)

Using this in the Blatter–Pattyn action, the SIA action becomes

(80)

Ignoring complications due to submerged boundaries , the first variation of Equation (80) with respect to the horizontal velocities yields the following Euler–Lagrange equations:

(81)

where . The boundary condition at the stress-free surface is

(82)

while the basal sliding boundary condition is given by

(83)

Observe that vertical velocity gradients depend directly on the frictional sliding coefficient β. This implies that the coefficient β used in the SIA approximation must be sufficiently large to satisfy Equation (78). As in the case of the Blatter–Pattyn model, solution of the shallow-ice system (Equations (81–83)) must be completed by obtaining the pressure and vertical velocity. Thus, the SIA is a system of uncoupled one-dimensional PDEs and associated boundary conditions at each point in the horizontal plane (i.e. in each vertical column), which is particularly easy to solve. This at least partly explains why this approximation has been extensively used in ice-sheet modeling (e.g. Huybrechts, 1994).

4.3.2. The shallow-shelf approximation

In the fast sliding regime, we expect the ice sheet to be only weakly affected by the basal boundary conditions, meaning that the vertical velocity gradients are relatively weak. This implies that

(84)

so we may approximate the Blatter–Pattyn strain-rate invariant (Equation (50)) by

(85)

Using this in place of in the Blatter–Pattyn action (Equation (61)), a fully 3-D action functional for the SSA becomes

(86)

which implies the existence of an effective SSA viscosity, . However, a 3-D model based on this action is not actually used in practice. Instead, it is possible to approximate and simplify it further.

Equation (84) implies that horizontal velocities are nearly independent of depth. To make use of this fact, we expand the velocity in a Taylor’s series about the vertical centroid (the mid-depth position in a column) and obtain

(87)

where is the vertical average of the horizontal velocity, and H(x, y) = z _s − z _b is the thickness of the ice sheet. Note that δ, λ≪ 1 in this approximation, making the truncation error essentially second-order. Following the procedure of Appendix B and expanding the action (Equation (86)) about vertical velocity averages, we obtain

(88)

where is the vertically constant strain-rate invariant obtained by substituting the vertical velocity average into Equation (85), A is the projection of the ice-sheet volume on the horizontal plane, and dA = dxdy, dV = dzdA. We have eliminated integration over depth since the integrand in the volume integral is depth-independent. For simplicity, in this section we drop the overbar notation and the special horizontal index notation, so, for example, u_i = (u, v) refers to the depth-averaged horizontal velocity. The Euler–Lagrange equations associated with Equation (88) become

(89)

where

(90)

is the effective 2-D stress tensor. Analogous equations appear in MacAyeal (1989) and appendix A of Schoof (2006).

Observe, however, that we are missing lateral boundary conditions for this 2-D problem. These are not needed provided H(x, y) = 0 at a lateral boundary, as we have implicitly assumed until now. However, when H(x, y) > 0, the ice sheet terminates in a vertical ice cliff and a lateral boundary condition is required. Assume we know the location of the lateral boundary given by the function B(x, y) = 0, a closed curve in the horizontal plane. We are free to choose the sign of this function such that the interior of the domain is specified by B(x, y) ≤ 0; the outward unit normal vector to the ice-sheet boundary is then given by The boundary condition at the cliff face is either stress-free (i.e. the cliff is still part of the ‘upper surface’) and given by

(91)

or else the cliff face is immersed in water and feels a resisting pressure (we ignore the complications of a grounding line, if it exists). In the latter case, the magnitude of the normal stress, directed inward, is given by

(92)

where is the opposing pressure at the cliff face, and there is no resistance in the tangential direction, which is specified by

(93)

If the cliff face is only partially immersed then Equation (91) will apply on the part of the face above water and Equations (92) and (93) apply on the submerged part. In this case, Equation (88) is modified as

(94)

where the area integration is over the entire horizontal domain of the ice sheet, dS is an area element on the cliff face, while the boundary segment corresponds to the immersed part of the cliff area where Equations (92) and (93) apply. The Euler–Lagrange equations are the same as in Equation (89), while the implied boundary conditions are given by

(95)

and these exactly correspond to Equations (91–93), where is the part of the cliff face that is not submerged, and S denotes the total cliff face area.

5. Discussion

The existence of an action principle for a non-Newtonian Stokes flow model for ice sheets has enabled us to derive a number of approximate models and the associated approximate action functionals directly and consistently from the fundamental Stokes action functional. Although these approximations already exist in the literature, there are important advantages to deriving them from a variational principle. The very existence of an action principle is important, particularly from the numerical point of view. As previously noted, derivation of the correct approximate equations and boundary conditions by means of a scaling analysis of the PDEs is non-trivial. However, the corresponding derivation from a variational principle is relatively simple and automatically produces consistent equations and boundary conditions. Given the existence of a variational principle, the approximate system preserves desirable properties of the exact system such as a positive-definite dissipation function and self-adjointness of the operator associated with the system of equations.

The principal approximation derived in this way is the first-order or Blatter–Pattyn approximation. We have shown that the derivation of the standard Blatter–Pattyn model involves a decoupling or splitting of the action functional into two parts. In the case of basal sliding, this decoupling can be justified only when basal slopes are very small. In other words, the standard Blatter–Pattyn model is valid only when both the low-aspect-ratio and small-basal-slope assumptions are valid. It follows that this also applies to the two subordinate approximations, the SIA and SSA.

The existence of a variational principle for ice-sheet dynamics provides a basis for new and potentially more efficient discretizations and improved numerical solution techniques. Here we briefly mention some of these implications, deferring a detailed discussion to a later publication. A potential solution technique, for example, is an analog of the classical Rayleigh–Ritz method in which the velocity is expressed as an expansion in trial functions that satisfy the incompressibility constraint and any specified velocity boundary conditions (Oakberg, 1981). The action functional would then be minimized with respect to the expansion coefficients. Note that the trial functions need not satisfy the stress boundary conditions since these are implicitly built into the functional. Alternatively, and more practically since such trial functions are not generally available, one can treat some of these problems (e.g. the Blatter–Pattyn model) as a problem in nonlinear optimization, i.e. the minimization of the action functional, and apply solution techniques developed specifically for such problems (e.g. Nocedal and Wright, 2006). The Stokes system, however, is not a minimization problem but involves the solution of a linear symmetric but indefinite matrix problem at each iteration (a so-called ‘saddle point’ problem; e.g. Benzi and others, 2006). Even here the existence of an action functional helps to ensure a symmetric and consistent matrix discretization. Nonlinear matrix problems (in this case arising from a nonlinear discrete action functional) are most often solved by a Newton–Raphson iteration method. For the Blatter–Pattyn model (and its approximations) the Newton–Raphson method involves the solution of a symmetric positive-definite matrix problem, which may be done efficiently using a preconditioned conjugate gradient algorithm. The hierarchy of approximations in Figure 1 suggests that the various approximations themselves might be used as physics-based preconditioners for the solution of problems at the next higher level of approximation. Since preconditioners for the conjugate gradient method must also be symmetric and positive-definite, variational principles for the Blatter–Pattyn model and the approximate models derived from it become very useful since they automatically produce preconditioners that possess these properties.

In addition, the existence of a variational principle has important implications for the actual discretization of the PDEs, as well as of the boundary conditions. Of particular significance is the fact that a properly formulated discrete action principle automatically includes the correct specification of boundary conditions. This is particularly important when the basal topography is not smooth (e.g. if it is described by discontinuous step functions), in which case even the form of the boundary conditions becomes indeterminate. Furthermore, there is no need to make a low-aspect-ratio approximation for the normal vectors at the boundary, as in the Blatter–Pattyn model, since these are automatically implied by the variational principle, thereby improving accuracy.

Finally, the self-adjoint property of a discrete forward model obtained from a variational principle, as described here, should facilitate the formulation of discrete inverse models. Using optimization techniques, these models can be applied to observations (e.g. ice-sheet geometry and surface velocity) to obtain best estimates for initial and boundary conditions and for other model parameters such as the value of the basal sliding coefficient, β (e.g. MacAyeal, 1993), or the distribution of the flow-law rate factor, μ ₀, or ice viscosity, μ (e.g. Vieli and others, 2006).