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 (Reference PatersonPaterson, 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. Reference SchoofSchoof, 2006, Reference Schoof2010). 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.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 (Reference BirdBird, 1960; Reference JohnsonJohnson, 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. Reference BirdBird, 1960; Reference Glowinski and RappazGlowinski 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 (Reference Bird1960) 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. (Reference CourantCourant (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.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. Reference Schoof and HindmarshSchoof 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, Reference Schoof and HindmarshSchoof 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. Reference Schoof and HindmarshSchoof 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 Reference BlatterBlatter (1995) and Reference PattynPattyn (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 (Reference PattynPattyn, 2003; Reference SchoofSchoof, 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, Reference Colinge and RappazColinge and Rappaz (1999) and Reference Glowinski and RappazGlowinski 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. Reference Rappaz and ReistRappaz 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. Reference SchoofSchoof (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 (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003; Reference SchoofSchoof, 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 Reference SchoofSchoof (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 Reference PattynPattyn (2003) and Reference SchoofSchoof (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 Reference SchoofSchoof (2010) in connection with more general frictional sliding stress laws, and it corresponds to those found in Reference Van der Veen and WhillansVan der Veen and Whillans (1989) and Reference PattynPattyn (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. Reference HuybrechtsHuybrechts, 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 Reference MacAyealMacAyeal (1989) and appendix A of Reference SchoofSchoof (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.