1. Introduction

A significant issue in ice-sheet modeling concerns the basal sliding boundary condition. This boundary condition is complicated by the fact that it involves two different types of conditions: a no-penetration condition in the normal direction, which is a Dirichlet boundary condition, and specification of the sliding stress in the tangential direction, which is essentially a Neumann condition. Further, it may be complicated in practice to accurately specify the required normal and tangential directions. In the past, this issue has been avoided, particularly in connection with approximate ice-sheet models, by assuming small basal slopes (e.g. Reference HutterHutter, 1983; Reference Van der Veen and WhillansVan der Veen and Whillans, 1989). However, basal slopes can be significant in practice. For reference, we note that a preliminary ∼12 000 km², high-resolution (∼1 km) basal topography dataset of the Jakobshavn region in Greenland (https://www.cresis.ku.edu/data/greenland) shows that ∼14% of the bed area has slopes in the range 5–10⁸, ∼6% has slopes of >10° (this includes the majority of the subglacial trench containing Jakobshavn Isbræ), and maximum slopes are as large as ∼43°. Thus, errors associated with a small slope approximation can be expected to impact ice-sheet simulations, and in particular those that aim to exploit high-resolution datasets of bed topography.

Various attempts have been made to account for the effects of finite bed slope in approximate models (e.g. the shallow-ice approximation (SIA)) that implicitly make the low-aspect-ratio approximation δ = d/L ≪ 1 (where d is a vertical scale on the order of the ice-sheet thickness and L is a characteristic horizontal length scale) and that also assume small bed slopes, |∂z _b/∂x| ≪ 1 (where z _b is the vertical coordinate of the bed surface). For example, Reference MorlandMorland (2000) and Reference SchoofSchoof (2003) consider shallow-ice dynamics driven by upper-surface perturbations, which permit local flow on much smaller length scales over ‘bumps’ in bottom topography. These models require explicit integration of the equations over the short topography length scales. Reference SchoofSchoof (2003), on the other hand, considers the problem of parameterizing the effects of significant basal topography for use in the original approximate models that depend on the small-basal-slope assumption.

In this paper, we consider the specification of arbitrary-slope basal boundary conditions in large-scale, three-dimensional (3-D) models incorporating non-Newtonian Stokes flow dynamics. These models, which do not involve the low-aspect-ratio approximation, are becoming feasible with the increasing availability of supercomputing capabilities (e.g. Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007). Although the analytic specification of various types of boundary conditions is relatively straightforward, this is not so in the case of discrete models, particularly when the boundary is not smooth as, for example, when it has to conform to a discrete grid. In this regard it is noteworthy that non-Newtonian ice-sheet or glacier dynamics may be described in terms of a variational principle. A variational formulation of this kind has many advantages over a standard partial differential equation (PDE) formulation (Reference Dukowicz, Price and LipscombDukowicz and others, 2010, henceforth referred to as DPL. The present paper is a companion paper to DPL. It is intended to be read independently; however, access to the earlier paper may be helpful). In particular, a variational formulation has considerable advantages for the specification of boundary conditions, especially when applied to numerical modeling. Many boundary conditions, Neumann boundary conditions in particular, are inherently contained within the variational principle, and these so-called natural boundary conditions need not be specified separately, i.e. outside the variational principle. Such boundary conditions are particularly advantageous from the numerical point of view because it is only necessary to discretize the variational functional and not the boundary conditions themselves. The discretization of the functional is relatively straightforward because it is a scalar quantity that, in general, involves lower-order derivatives than the differential equations or the boundary conditions. Dirichlet boundary conditions, on the other hand, typically are not natural boundary conditions and must be specified separately from the variational principle. In particular, the no-penetration basal boundary condition suffers from this problem. There are various possible ways of implementing Dirichlet boundary conditions, some briefly discussed by DPL. For example, in the finite-element method this is commonly accomplished by the use of penalty functions. Many of these methods either are approximate or specify such boundary conditions outside the variational principle. Our preference is to incorporate all boundary conditions into the functional in order to preserve the advantages of a variational formulation. It is indeed possible to do this for Dirichlet boundary conditions, and specifically for the exact no-penetration boundary condition, by means of the method of Lagrange multipliers (i.e. by treating the Dirichlet boundary condition as an additional constraint). This transforms the no-penetration condition into a natural boundary condition with all the above-mentioned advantages as a byproduct. Further, no small-slope approximation is required.

Frequently, however, a simpler and less expensive approximate model of Stokes flow may be appropriate. In particular, the commonly used first-order or Blatter–Pattyn approximate model (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003) is accurate within the low-aspect-ratio approximation of the Stokes system (Reference Schoof and HindmarshSchoof and Hindmarsh, 2010; hereafter referred to as SH). However, the specification of discrete boundary conditions is again problematic, just as in the Stokes model. Furthermore, the associated basal boundary conditions implicitly assume a small bed slope, |∂z _b/∂x| ≪ 1, as in the SIA. Fortunately, the Blatter–Pattyn model also possesses a variational principle (see Reference SchoofSchoof, 2010; DPL). This allows us to apply the Lagrange multiplier method to incorporate the no-penetration boundary condition into the functional, resulting in boundary correction terms that are second-order accurate with respect to the basal topographic slope. This modification, which we call the ‘modified Blatter–Pattyn model’, improves the accuracy of the standard Blatter–Pattyn model for finite basal slopes. It does not, however, remove the error entirely.

2. The Basic Stokes Model

We begin with the variational principle for the basic Stokes model whose action functional is given by

where u_i ∊ {u, v, w} is the velocity vector, x_i ∊ {x, y, z} is the position vector, g_i is the gravitational acceleration vector, typically g_i = (0, 0, −g), ρ is the ice density, assumed constant, and is the second invariant of the full Stokes strain-rate tensor,

giving

We define

where

is the effective Glen’s law viscosity, typically used with exponent n = 3, and μ ₀(θ) is a temperature-dependent coefficient. The variable, P, in Equation (1), although commonly identified as the isotropic pressure, is actually a Lagrange multiplier for enforcing incompressibility, ∂u_i /∂x_i = 0. The volume integral includes the entire ice sheet. Note that we are using Cartesian tensor notation and, where appropriate, the summation convention on repeated indices. In general, tensor indices are 3-D, i.e. i,j, … ∊ {x, y, z}, except that an index appearing in parentheses will denote the horizontal plane, i.e. (i), (j), … ∊ {x, y}. For example, u_iu_i = u ² + v ² + w ² and u _(i) u _(i) = u ² + v ².

The variational principle states that the dynamical system is to be found at the extremum of the action, Equation (1), with respect to functional derivatives relative to velocity, u_i , and pressure, P, as follows:

This yields the following Euler–Lagrange equations:

(a) a 3-D momentum equation,

where σ_ij = τ_ij – Pδ_ij is the Cauchy stress tensor and is the deviatoric stress tensor,

(b) the continuity equation for incompressible flow,

and

(c) stress-free boundary conditions over the entire surface, S, that surrounds the ice-sheet volume, V,

DPL provide details of the derivation.

As shown by Equation (8), the action, Equation (1), incorporates only stress-free boundary conditions. For present purposes we assume that the surface, S, is composed of an upper surface, S ^(s), given by z = z _s(x, y), and a basal surface, S ^(b), given by z = z _b(x, y), that is in contact with the bed: S = S ^(s) + S ^(b). The corresponding outward-pointing unit normal vectors are given by

1. (i) Upper surface

   (9)

   

2. (ii) Basal surface

The action must be modified to include basal sliding, with or without frictional slip, and we do this in two steps. First, we incorporate a frictional stress force in the functional, as follows:

Note that the surface integral covers only the basal surface, S ^(b); it is assumed that the frictional force acts only on S ^(b) so that there is no contribution from the stress-free upper surface. For simplicity, as in DPL, we assume that Σ _j (u) = −βu_iu_in_j/2, which represents a linear frictional sliding law with constant coefficient, β (β ≥ 0). Other frictional laws may be easily accommodated, as in Reference SchoofSchoof (2010), for example. The variational principle applied to Equation (11) gives the following boundary conditions:

(d) a stress-free condition on S ^(s):

and

(e) a frictional condition on S ^(b):

In actuality, the basal condition involves sliding along the bed. This means that we must require the no-penetration condition,

in addition to Equation (13), which then becomes a tangential boundary condition. Note that Equation (14) is a Dirichlet boundary condition since it specifies the normal velocity, while (Equations 12) and (13) are Neumann-type conditions since they specify velocity gradients.

As discussed by DPL, a variational principle, such as one based on Equation (1) or (11), incorporates Neumann boundary conditions as natural boundary conditions but not necessarily Dirichlet boundary conditions. Leaving aside penalty-function methods, Dirichlet boundary conditions may be incorporated in the functional in two ways. First, velocities may be restricted to those that satisfy the Dirichlet boundary conditions. In practice, this option is feasible only in the discrete case where velocities are expanded in local basis functions that satisfy the specific Dirichlet boundary conditions. In the ice-sheet problems that we are considering, all boundary conditions are Neumann conditions except for the basal no-penetration boundary condition, which specifies only one component of the basal velocity. This situation may be problematic to specify in terms of local basis functions. The second alternative is to enforce Dirichlet boundary conditions as constraints by means of Lagrange multipliers. This alternative is attractive because variational principles are well suited to the use of this method, including the situation when only one component of velocity is specified, and also because all boundary conditions then become natural boundary conditions that are inherently part of the variational principle. As a result, boundary conditions become easier to implement in the discrete case. Accordingly, we pursue the second alternative here.

To accommodate the second option, we modify the action (Equation (11)) as follows:

where the first term in the surface integral has been added to enforce the no-penetration constraint at the basal surface, u_in_i = 0, by means of the Lagrange multiplier, Λ. As before, taking the variation of Equation (15) with respect to velocity components u_i and the Lagrange multipliers, P and Λ, we obtain the same Euler–Lagrange equations and stress-free boundary condition on the upper surface, S ^(s), as before, except that at the basal surface, S ^(b), we have

(f) the no-penetration condition,

and

(g) a constrained basal stress condition,

So far, we have mostly omitted the superscript notation denoting upper and basal surfaces when discussing boundary conditions, taking for granted that this should be clear from the context. Note that we have added a new variable, Λ, and an additional boundary equation, Equation (16). However, this new variable appears only locally at the basal boundary and, as a result, the two basal boundary equations (16) and (17) are easily solved for the combination Λ – P, as follows:

Substituting this in Equation (17), we have

This is exactly the tangential stress sliding boundary condition for the Stokes system written out in detail by DPL, except that now it is incorporated in the action as a natural boundary condition. The basal slope is not restricted. Note that Equation (19) specifies only the tangential components; the normal component of this equation vanishes because of Equation (16). Thus, Equations (16) and (17) are an alternative (and equivalent) way of specifying the basal frictional sliding boundary condition.

3. A Modified Blatter–Pattyn Model

The commonly used first-order or Blatter–Pattyn approximate model (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003) is specified by the following equations:

where in the case of Glen’s law the effective Blatter–Pattyn viscosity, μ _BP, is given by

and is the approximate Blatter–Pattyn strain-rate invariant:

As before, the boundary conditions used in the Blatter–Pattyn model are a stress-free condition at the upper surface and a frictional sliding condition at the bed. These are given by

(a) Stress-free upper surface:

(b) Basal frictional sliding:

for the simple case of a linear drag law. This boundary condition approximates the true horizontal components of tangential frictional sliding stress at the bed only in the small slope limit |n_x |, |n_y | ≪ 1, as we show later. Thus, the Blatter–Pattyn model is valid only in the small slope limit. We shall refer to the model that incorporates the above boundary conditions as the standard Blatter–Pattyn model. The more general form of the boundary conditions containing exact unit normal vectors, i.e. Equations (24–27), is found in Reference SchoofSchoof (2010), for example. However, in practice the model is most commonly used with the Reference Van der Veen and WhillansVan der Veen and Whillans (1989) boundary conditions, obtained by further approximating the unit normal vectors by means of the small-slope approximation, |∂z _s/∂x _(i)|, |∂z _b/∂x _(i)| ≪ 1 which results in

Given the surface and basal boundary conditions, Equations (20) and (21) determine the horizontal velocities, and this determines the vertical velocity by means of the continuity equation:

DPL provide more details.

Just as in the Stokes case, the standard Blatter–Pattyn model is obtainable from a variational principle, this time associated with the action

The action principle associated with Equation (31) is exactly equivalent to the equations of the standard Blatter–Pattyn model, Equations (20–27). Essentially the same functional as Equation (31) has appeared in finite-element ice-sheet formulations such as Reference SchoofSchoof (2010), where it is inferred from the known Blatter–Pattyn equations. Compared with the Stokes model, the Blatter–Pattyn model is much simpler because it only involves horizontal velocities. Furthermore, Equation (31) is positive-definite and therefore leads to a symmetric, positive-definite (although nonlinear) matrix problem. Unlike in the Stokes case, the action principle now corresponds to an actual minimization problem (rather than a saddle-point problem). As mentioned earlier, presupposing that basal slopes are , SH have shown that the standard Blatter–Pattyn model is first-order accurate in the low-aspect-ratio approximation to the Stokes model; that is, it is accurate to O(b). This implies that the Blatter– Pattyn action, Equation (31), and the basal boundary condition, Equations (26) and (27), are also first-orderaccurate low-aspect-ratio approximations.

In this regard, DPL have shown that the Blatter–Pattyn basal boundary condition, Equations (26) and (27), approximates the correct tangential boundary condition, Equation (19), only in the limit of small basal slopes, i.e. when Equation (29) is valid. Note that the low-aspect-ratio and small-basal-slope approximations are related but do not necessarily imply one another. It might be expected that a sufficiently large basal slope will cause the velocity field to locally violate the low-aspect-ratio approximation, which requires

However, this does not mean that a small basal slope, one that satisfies Equation (29), is required, and in fact a moderately large basal slope may be compatible with the low-aspect-ratio approximation.

The action, Equation (31), was derived directly from the Stokes action, Equation (11), by DPL. In this derivation the no-penetration basal boundary condition was only partially taken into account, and the small-basal-slope assumption was not specifically invoked but was presumably implied. This suggests that explicitly incorporating the no-penetration boundary condition into the Blatter–Pattyn action might be beneficial. Therefore, noting the modifications that were necessary in the Stokes model, the standard Blatter–Pattyn action, Equation (31), may be modified as follows:

where we use the subscript MBP to indicate that this differs from the standard Blatter–Pattyn action in two ways. First, just as in the Stokes case, there is a term in the surface integral imposing the no-penetration condition, and second, the term specifying frictional stress involves the full 3-D velocity rather than just the horizontal velocity. As a result, the variation is to be taken with respect to the 3-D velocity. Note that the strain invariant has the same form, Equation (23), as in the standard Blatter–Pattyn model, and we indicate this by using the notation

Taking the variation of the modified action, Equation (33), the Euler equations (20) and (21) and the stress-free surface boundary condition, Equations (24) and (25), are unchanged. (Note that a more general upper surface boundary condition can be found in DPL.) The basal boundary condition, on the other hand, now becomes

and

where

is an effective Blatter–Pattyn stress quasi-tensor. Just as in the Stokes case, we may solve for A by contracting Equation (35) with to obtain

Substituting this in Equation (35) and expanding, the horizontal components become

where we have made use of the fact that to eliminate . This is exactly the same as the correct tangential component of the sliding boundary condition derived by DPL for a generic low-aspect-ratio approximate model, except that now it appears as a natural boundary condition in the modified Blatter–Pattyn variational principle. It is clear that these equations reduce to the standard Blatter–Pattyn boundary condition, Equations (26) and (27) in the limit of very small basal slopes, as mentioned earlier. In addition, the modified Blatter–Pattyn action provides the correct vertical velocity at the basal surface, in contrast to the standard Blatter–Pattyn action, which specifies only the horizontal velocities. This simplifies the determination of internal vertical velocities directly from the continuity equation, rather than by the use of Equation (30).

4. Analysis of the Modified Blatter–Pattyn Model in the Slow Sliding Regime

We now address the extent to which the modified boundary condition, Equations (38) and (39), improves the accuracy of the Blatter–Pattyn model compared with the standard boundary condition, Equations (26) and (27), the models being otherwise identical. We first simplify by limiting ourselves to the case of two-dimensional flow in the x-z plane with constant viscosity, so that, for example, μ _BP = μ, a constant. The analysis can be done most conveniently for an essentially horizontal ice sheet in the slow sliding regime, i.e. where the term ∂u/∂z dominates all other components in the strain-rate invariant, Equation (2). There are two advantages of doing this. First, the analysis is straightforward and can be carried out analytically, and second, this simplification can be applied to all three models, including the reference Stokes model, so that a proper evaluation of the two approximations can be made. Note that this is the regime in which the SIA applies, although we do not make the assumption of small basal slope, [∂z _b/∂x] ≪ 1, that is commonly made in this case.

From the action principles given earlier, in this regime the relevant equations for all three cases are given by

together with the boundary condition at the stress-free surface,

However, the basal boundary condition is different for each model. We have

1. (a) the standard Blatter–Pattyn model (SBP),

   (42)

   

   which is given by Equation (26),

2. (b) the modified Blatter–Pattyn model (MBP),

   (43)

   

   which is given by Equation (38), and

3. (c) the Stokes model (S),

   (44)

   

   which is implied by Equation (19).

In each vertical column this system amounts to a simple, one-dimensional ordinary differential equation with constant coefficients. Its solution is given by

where

and use has been made of the fact that . Note that we must have |n_x | ≤ 1, i.e. a maximum basal slope of 90°, far in excess of realistic basal slopes. We therefore limit consideration to slope inclinations of <45°. Also note that, consistent with this limit, |f(n_x )| ≤ 1 in all cases. Given the horizontal velocity, Equation (45), the vertical velocity, w(x,z), may be obtained from Equation (30).

It is important to note that the horizontal variation of the velocity in Equation (45) is implied by the horizontal variation of z _s and z _b (therefore of n_x also). Thus, it can always be arranged for the horizontal gradients of these quantities to be sufficiently small that the assumption regarding the dominance of ∂u/∂z is valid. That is, |∂u/ ∂z| ≫ |∂w/∂x|, |∂w/∂z|, where we make use of the fact that |∂u/∂x| = |∂w/∂z| from the continuity equation. Note that this implies the standard small aspect ratio assumption, which only requires |∂u/∂z| ≫|∂w/∂x| as in Equation (32), but not vice versa.

We are now in a position to evaluate the error associated with the two Blatter–Pattyn models relative to the Stokes model, at least in the regime where the SIA is valid. From Equation (45), the basal (or sliding) velocity is given by

where H = z _s – z _b is the ice-sheet thickness. It is convenient to use this velocity as a metric to evaluate the accuracy of the two approximate models as a function of the basal slope. We define a normalized basal horizontal velocity, using Equation (47), as follows:

where

is the basal velocity for a flat horizontal bed (n_x = 0). In Figure 1 we plot the normalized sliding velocity from Equation (48) as a function of the inclination of the basal surface (°), given by . In this regime the basal velocity error of the standard Blatter–Pattyn model relative to the reference Stokes model is quite large, except at very small basal slopes. The modified Blatter– Pattyn model cuts this error by approximately a factor of two for slope inclinations smaller than 45°.

Fig. 1. Normalized basal velocity in the slow sliding regime

From Equation (45), the upper surface velocity is given by

A requirement for being in the slow sliding regime is that ∂u/∂z dominates all other terms in the strain-rate invariant. Alternatively, this may be interpreted as requiring that |u _s| ≫ |u _b|, or equivalently, that

This is always satisfied provided

because, as we remarked earlier, |f(n_x )| ≤ 1. Since ice-sheet viscosity is very high, this is a requirement that β must be large for the ice sheet to find itself in the slow sliding regime.

5. Conclusions

We have shown how to extend the variational principle for the non-Newtonian Stokes model of ice-sheet dynamics so that all boundary conditions, including various forms of the sliding boundary condition, are incorporated as natural boundary conditions in the action functional. This substantially simplifies the formulation of accurate and consistent discrete models that are amenable to efficient numerical solution.

We have applied the same methodology to the approximate Blatter–Pattyn model, normally valid only in the limit of small basal slopes. Comparing the accuracy of the standard and modified Blatter–Pattyn models relative to the Stokes model in the slow sliding regime, we find that the error of the standard Blatter–Pattyn model can be large, and this error increases with increasing basal slope. The modified Blatter–Pattyn model cuts this error approximately in half, extending the range of validity of this model to somewhat larger values of basal slope at negligible additional cost. It does not, however, eliminate that part of the Blatter–Pattyn error associated with the small-aspect-ratio approximation. Although the analysis presented here applies only to the slow sliding regime, it suggests that similar results might be found in other regimes (e.g. fast sliding). This suggestion, however, remains to be verified numerically.