1. Introduction

Realistic predictions of the ice-sheet contribution to sea-level rise, both in the past and in the coming centuries, require numerical ice-sheet models that include both fast flow near the ice-sheet margins and more traditional sluggish ice-sheet-like behavior in the ice-sheet interior (Reference Lemke and SolomonLemke and others, 2007). The challenge is to do so in a straightforward and computationally feasible way.

Early attempts to include both fast and slow modes of ice flow in an ice-sheet model employed a mixture approach, a shotgun marriage between two very different thin film approximations: the shallow-ice approximation (SIA) and the shallow-stream approximation (e.g. Reference Marshall and ClarkeMarshall and Clarke, 1997; Reference Hulbe and MacAyealHulbe and MacAyeal, 1999). In the interior of the ice sheet, where the ice is primarily frozen to the bed or sliding very slowly, the SIA is assumed to hold (e.g. Reference HutterHutter, 1983; Reference FowlerFowler, 1997). In this limit, the pressure gradient is approximately balanced by gradients in vertical shear stresses, and adjacent columns of ice interact solely through the mass continuity equation. Hence, the velocity field is effectively two-dimensional (2-D) and can be determined independently for each column of ice. In contrast, in regions where rapid basal sliding occurs, lateral stretching becomes important, and the pressure gradient is primarily balanced by gradients in (depth-averaged) lateral stresses. This approximation, often called the shallow-shelf approximation or shelfystream approximation (SSA), was originally derived by Reference MacAyealMacAyeal (1989) to describe ice-stream flow and has proven to be successful in describing ice streams and , at least qualitatively, in explaining the acceleration of glaciers in response to a decrease in ice-shelf buttressing, a decrease in basal friction and increased basal melting (Reference MacAyealMacAyeal, 1989; Reference Schmeltz, Rignot, Dupont and MacAyealSchmeltz and others, 2002; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004; Reference Dupont and AlleyDupont and Alley, 2005). In the SSA the velocity is approximately constant with depth and , like the SIA, is effectively 2-D. However, the SSA neglects all vertical shear within the ice, an unrealistically severe assumption throughout the majority of an ice sheet.

In the mixture approach, the SSA and SIA are coupled together solely through the continuity equation, neglecting any transfer of momentum. This, combined with the need to a priori specify which portions of the ice sheet are sliding, has led researchers to abandon the SIA and SSA in favor of fully three-dimensional (3-D) models. These models either solve the full Stokes equations for all components of velocity and pressure simultaneously, or some long-wavelength approximation thereof (e.g. Reference BlatterBlatter, 1995; Reference Payne, Vieli, Shepherd, Wingham and RignotPayne and others, 2004; Reference Pattyn, Huyghe, De Brabander and De SmedtPattyn and others, 2006). Fully 3-D models are more general, but the greater generality is bought at a steep price: these models are computationally expensive and much more difficult to implement numerically. All computational challenges are magnified when considering transient ice-sheet behavior, especially over the millennial timescales of glacial and orbital cycles. To put things in perspective, since ice dynamics are highly sensitive to basal conditions and even under the best conditions these are poorly known, the increased accuracy of full Stokes models compared to more efficient reduced models may be illusory.

In this study, instead of wholly abandoning the SIA, we attempt to marry the SIA and SSA in an approximation whereby mass and momentum transfer is continuous across different flow regimes. More specifically, we seek an approximation that is more accurate than mixture models yet easier to implement and more computationally efficient than either full Stokes or so-called higher-order ‘Blatter’ models (Reference BlatterBlatter, 1995; Reference Pattyn, Huyghe, De Brabander and De SmedtPattyn and others, 2006). Our efforts in this direction are modest and build upon those of many others (Reference BuddBudd, 1970; Reference HutterHutter, 1981; Reference Kamb and EchelmeyerKamb and Echelmeyer, 1986; Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others, 1998; Reference HindmarshHindmarsh, 2004; Reference Pollard, DeConto, Hambrey, Christoffersen, Glasser and HubbardPollard and DeConto, 2005; Reference Bueler and BrownBueler and Brown, 2009). In particular, our model is similar in spirit to the model first proposed by Reference Kamb and EchelmeyerKamb and Echelmeyer (1986) and later classified as an ‘L1L2’ model by Reference HindmarshHindmarsh (2004). The key assumption we make is that the SIA and SSA are valid long-wavelength approximations within their respective domains of applicability and the primary issue is that of coupling the two limits together such that traction is continuous between different flow regimes.

The starting point of our investigation is a form of mechanics not often employed in glaciology. Instead of using the conventional Newtonian vectorial framework, we use an alternative, variational form of mechanics based on Hamilton’s principle (see, e.g., Reference Landau and LifshitzLandau and Lifshitz, 1976; Reference LanczosLanczos, 1986). This break in tradition has three reasons. The first reason is that the approximation we present is based on the Rayleigh–Ritz method which relies on the existence of a variational principle as its starting point. The second reason is more philosophical in nature. Although variational principles have played a significant role in several theoretical advances in ice-sheet dynamics (Reference FowlerFowler, 1979; Reference Picasso, Rappaz, Reist, Funk and BlatterPicasso and others, 2004; Reference SchoofSchoof, 2006a,Reference Schoofb, Reference Schoof2010) they are often presented as an ad hoc mathematical trick. This misses the deep physical connection between the variational principle, classical mechanics and ice-sheet energetics. The third reason that we introduce Hamilton’s principle is that there have been significant theoretical advances in the treatment of the dynamics and thermodynamics of complex, nonlinear fluids, especially the application of dual variational principles to derive bounds on behavior (Reference MorrisonMorrison, 1986; Reference Grmela and ÖttingerGrmela and Öttinger, 1997; Reference Öttinger and GrmelaÖttinger and Grmela, 1997). Applying these relatively new techniques to ice sheets requires an in-depth understanding of the connection between variational principles and the energetics of ice-sheet flow.

Since we anticipate that most glaciologists are unfamiliar with Hamilton’s principle, we start by providing an introduction to this formulation of mechanics. This allows us to introduce our notation and show how a variational principle for ice-sheet flow can be constructed solely from ice-sheet energetics. This treatment is analogous to the force-balance arguments by which the Stokes equations for a non-Newtonian flow are usually derived. The variational formulation of non-Newtonian Stokes flow is not new (see, e.g., Reference FinlaysonFinlayson, 1972) and those readers already familiar with or uninterested in this extensive background material may wish to skip ahead to section 2.3. It is here that we start to develop an approximation for large-scale ice-sheet flow that includes both lateral stresses and vertical shear stresses. In section 5 we use the approximation developed in the previous section to determine the steady-state profile of ice sheets for a range of basal conditions and Newtonian and plastic sliding laws.

To put this study in historical perspective, the approach we advocate has much in common with the spectral approach employed by Reference OerlemansOerlemans (1982) to approximate the solution of the temperature equation in early thermomechanical ice-sheet models. We differ from that much earlier study in that we apply the spectral approach directly to a variational principle. To our knowledge, the only instance where approximate equations of motion for ice-sheet flow were derived from an energy principle is the work of Reference MacAyeal and ThomasMacAyeal and Thomas (1982), where the equations governing ice-shelf flow were derived using d’Alembert’s principle, a close relative of Hamilton’s principle.

2. Hamilton’s Principle

2.1. Formulation of Hamilton’s principle

Let the position of marked fluid particles at time τ be given by

The coordinates a = (a, b, c) are continuous variables that ‘label’ individual fluid ‘particles’that are small compared to the macroscopic dimension of the ice sheet but large compared to the individual grain size of ice. The label coordinates allow fluid particles to be tracked over time and can be chosen in the most convenient way so long as labels are conserved. For reference, the coordinate system and the geometry of the ice sheet used are defined in Figure 1.

Fig. 1. The coordinate system and geometry considered in this paper. We choose the Cartesian coordinate system such that z = 0 corresponds to sea level and we measure the effective depth of the water, D, relative to the ice-sheet base.

We use the symbol τ as a shorthand for time t. However, ∂/∂τ implies differentiation with respect to time with label coordinates (a, b, c) held fixed, and ∂/∂t implies differentiation with spatial coordinates (x, y, z) held fixed. With this in mind, the velocity of fluid particles is by definition:

We also suppose that the inverse of x (a, b, c), y (a, b, c), z (a, b, c) exists such that we can map from a-space back to x-space and vice versa. Hence the Jacobian cannot be singular:

Following Salmon (Reference Salmon1983), we assign label coordinates such that equal areas in a-space contain equal mass. Denoting the density of ice associated with point (x, y, z) by ρ we have

The specific volume, α, is then by definition

(1)

If the density is constant, this implies that the specific volume is also constant. For constant α, application of ∂/∂τ to Equation (1) shows that

(2)

where we have used the Jacobian to switch from material to Eulerian coordinates. This, however, is just the incompressibility condition. Thus a constant Jacobian implies an incompressible fluid.

Hamilton’s principle states that of all possible trajectories of a mechanical system, the actual trajectory is that which results in a stationary value of the time integral of the kinetic energy plus the work done by/to the system

(3)

where T and W are the kinetic energy and work densities, δ stands for arbitrary independent variations in the position of fluid particles relative to fixed Cartesian coordinates and L is the Lagrangian functional (see, e.g., Reference LanczosLanczos, 1986). The integral is over all fluid particles. The equations of motion that result from rendering stationary the functional L are called the Euler–Lagrange equations and correspond to Newton’s laws of motion, i.e. conservation of mass and momentum. The problem is that for dissipative systems there may be no such functional L whose variation is given by Equation (3).

In the simplest case, all forces acting on the system are conservative and hence path-independent. In this case, L is called the ‘action’. Our interest is in situations which include non-conservative forces. We start by noting that regardless of whether the forces are conservative or not, the differential work δW (sometimes called the virtual work) can always be expressed as the product of a force times a differential displacement. Let the components of the force in the (x, y, z) directions be given by (X, Y, Z). The differential work is then

(4)

where the integral is over all particles within the fluid.

Now suppose that the work can be separated into a conservative portion V and a non-conservative portion W _d such that

(5)

where we have introduced the negative sign because work done by the system is negative. Since the conservative portion of the work depends only on the current configuration and not how it arrived at that configuration, we can write this portion of the work as a scalar work function V. However, the non-conservative work, δW _d, is path-dependent and can only be expressed as the product of a force (or stress) times a differential displacement. The total non-conservative portion of the work W _d in changing the configuration of the system of particles is obtained by adding up all of the contributions from incremental displacements of the particles as they are moved from the initial to the final configuration. For non-conservative, dissipative forces, the path dependence of the differential δW _d implies that there is no scalar work function and the total amount of work done depends on the particular path integral of deformation. Hamilton’s principle is then only applicable in the differential form of Equation (3).

2.2. A modified form of Hamilton’s principle

If we consider an interval of time that is sufficiently short that the geometry of the ice changes very little and hence the kinetic energy is nearly constant, we can approximate the time integral in Equation (3) using a Taylor series:

(6)

where we have used the fact that the variation operator δ commutes. We can express this as the time rate of change of the Lagrangian functional:

(7)

where the dot decoration represents differentiation with respect to time τ.

Following Reference SalmonSalmon (1998), we suppose that the conservative portion of the work V consists of two terms: (1) internal energy E, which is assumed to be only a function of the specific volume α; and (2) external potential energy φ, which only depends on the position of fluid particles. With these assumptions, the conservative portion of the work becomes

(8)

whence after differentiation with respect to time we have

(9)

If we define the pressure

(10)

the internal energy becomes

(11)

Further generalizations of the internal energy to account for entropy and/or for multiple fluid components (e.g. water and ice) are possible but are not the focus of this study. However, it is pleasing to note that the definition of pressure given by Equation (10) is consistent with the usual definition from thermodynamics (Reference SalmonSalmon, 1998).

Returning to the non-conservative portion of the work , we partition this into the rate of working done by viscous stresses and the rate of working done by frictional forces along the ice-sheet base . Although we have already determined that no work function exists for dissipative forces, if we denote incremental change in position of fluid particles relative to an instantaneous reference state with the vector δξ = (δx, δy, δz), then recalling Equation (4) the rate of working done by viscous stresses is given by

(12)

while the rate of working due to frictional forces is given by

(13)

In Equations (12) and (13), τ_ij represents the viscous stresses, τ _bi represents the basal tractions due to friction in the ith direction, ∂Ω_b represents the basal surface area and summation over repeated indices is implied.

To compute the rate of working, it is convenient to define ‘generalized potentials’ sometimes called dissipation functions (originally proposed by Lord Rayleigh in 1873). These dissipation functions yield the appropriate variational rate of working and , but due to the path dependence of work are not simply related to their corresponding work functions except under very restrictive circumstances. We denote the generalized potential for viscous dissipation by Г and that for friction by Π.

Assembling the previous components together, we now have the functional

(14)

needed for an integral variational principle. Unfortunately, it is only the variation of the dissipation functions Γ and Π that has a clear physical meaning. To apply the variational principle to ice sheets, we must consider explicit forms for each of the terms in Equation (14). We do this next.

2.2.1. Internal energy

Recalling Equation (10), the rate of change of internal energy is a linear function of the rate of change of the specific volume,

where the constant of proportionality P is defined to be

The integral over all particles of the rate of change of internal energy is then simply

(15)

where we have used the Jacobian to transform the integral from a-space to x-space.

2.2.2. Potential energy

The gravitational potential energy φ of the ice is related to the height z of fluid particles and can then be written as

The integral over mass in label-space coordinates is then simply

(16)

Differentiating Equation (16) with respect to τ we find

(17)

where we have again used the Jacobian to transform the integral from a-space to x-space.

The above expression does not account for any changes to the gravitational potential energy due to the influence of the ice-sheet bed (or sides). To account for the fact that these may exert an upward force on the ice sheet, we introduce a potential energy at the bed of the form:

(18)

where r = b − z is the bed penetration distance and ∂Ω_b represents the ice-sheet bed surface area. (One could postulate more general relationships where the potential barrier is a function of all particle coordinates.) The potential could account for a variety of ice-sheet–bed interactions including the change in gravitational potential of water beneath the ice if the ice is floating, or elastic potential energy if the ice is in contact with bedrock or till.

Differentiating Equation (18) with respect to τ and making use of the fact that

we find

(19)

where the Greek index v runs over the horizontal coordinates x and y. If the potential V _b is a very rapidly increasing function of the penetration distance, then the bed can be approximated as ‘rigid’ and we may impose the kinematic no-penetration constraint

(20)

or, more succinctly, u_in_i = 0, where n_i is the outward-directed normal vector. Introducing a Lagrange multiplier K to enforce the constraint, the potential energy becomes

(21)

The last term in Equation (21) accounts for the kinematic constraint that along a rigid ice-sheet bed the vertical velocity is related to the horizontal velocity through the kinematic condition u_in_i = 0. We shall use this throughout, but alternative specifications of the ice-sheet–bed interaction potential are possible.

2.2.3. Viscous dissipation

Defining the instantaneous strain-rate tensor in the usual way

(22)

and recalling the previous discussion of non-conservative forces, we now introduce a generalized viscous potential Γ, defined such that

(23)

It is often assumed that the rheology of ice is approximately independent of the first and third strain-rate invariants, so we assume Γ depends exclusively on the second strain-rate invariant (Reference PatersonPaterson, 1994),

We can either view Equation (23) as the definition of the dissipation function or adopt a more thermodynamic viewpoint where viscous stresses τ_ij are defined by Equation (23).

The variation δΓ is straightforward to compute using the calculus of variations:

(24)

which shows that the variation of Γ yields the appropriate differential rate of work needed to apply Hamilton’s principle. Although the dissipation function itself does not have a convenient physical interpretation, its variation can be interpreted as the differential rate of working. Thus Γ can be used as a functional that yields the appropriate virtual rate of working even though the path dependence of dissipation precludes a convenient interpretation in terms of the rate of dissipation.

Since the rheology of ice is often assumed to have a power-law dependency on σ (Paterson, Reference Paterson1994) we adopt a dissipation function of the form

(25)

with n usually taken to be ∼3 (Hutter, Reference Hutter1983; Paterson, Reference Paterson1994). We can write the dissipation function more conveniently in terms of viscous stresses and strain rates:

(26)

This shows that Γ is proportional to the rate of dissipation with the proportionality determined by the flow-law exponent n. For later use, we also define the effective viscosity

(27)

such that the viscous stress tensor can be written .

2.2.4. Frictional dissipation

The last term we need is the rate of working by frictional forces. This can be treated in a similar way to the rate of work done by viscous deformation of the ice, defining the frictional dissipation function II, this time a function of the sliding speed tangent to the ice-sheet bed u _b. Assuming a power-law parameterization, we can write the frictional dissipation function as a functional of the sliding speed u _b:

(28)

where γ is a (small) constant parameter which maintains the positive definiteness of the dissipation function, p is a friction exponent and β is a spatially varying coefficient that may depend on, for example, pore pressure, sediment grain size, etc.

Assuming velocity components at the bed are independent, the traction at the ice-sheet base due to friction in the ith direction, τ _bi , is related to Π by

(29)

where u _bi is the ith velocity component at the bed. However, the velocity components may not be independent at the bed.For instance, for a rigid bed where w = u_xb_x + u_yb_y we would have

(30)

where we have muddied our notation by using Greek subscripts on b to indicate differentiation with respect to x and y while Greek subscripts appended to u are used to indicate (x, y) components of u (not derivatives). For a rigid bed, the basal traction would be given by:

(31)

which it should be noticed is tangent to the bed. Regardless of this complication, for p = 1, Equation (29) reduces to a Newtonian friction:

For p = 1/2 we have:

For small γ ² this approximates a Coulomb-plastic material (Schoof, Reference Schoof2006b). Intermediate values of p correspond to a sliding law similar to that proposed by Budd and others (Reference Budd, Keage and Blundy1979).

2.3. Lagrangian formulation of ice-sheet flow

Assembling the individual components from section 2.2, the Lagrangian functional can be written explicitly:

(32)

This is identical to the variational principle first proposed by Johnson (Reference Johnson1960) and used by Fowler (Reference Fowler1979) for the first time in a glaciological context.

It would appear that because ice is very nearly incompressible, we should be able to neglect the first term in Equation (32). This is not the case. A perfectly incompressible fluid would imply a constant thermodynamic pressure. Instead, the pressure in Equation (32) can be thought of as a Lagrange multiplier, the presence of which enforces the approximation of a constant Jacobian, or more intuitively, incompressibility (Lanczos, Reference Lanczos1986).

2.4. Flowline Lagrangian formulation of ice-sheet flow

For simplicity, in the calculations that follow, we restrict our attention to a flowline and define the two velocity components (u_x , u_z ) ≡ (u, w) such that the horizontal velocity u is in the x direction and the vertical velocity w is in the z direction (Fig. 1). Using Equation (26) to rewrite the viscous dissipation function in terms of the viscous stresses τ_ik and rearranging, we have

(33)

The viscous stresses themselves should be understood to be nonlinear functions of the velocities (u, w) and we have appended a subscript to remind ourselves that this Lagrangian is now an approximation of the full 3-D Lagrangian.

The Lagrangian can be simplified slightly by insisting that Equation (33) vanish for arbitrary δP, yielding the incompressibility condition

Next, making use of the incompressibility condition in the form

to eliminate τ_zz from Equation (33):

(34)

We have labelled the new Lagrangian to emphasize that the new Lagrangian differs from the original Lagrangian. We have also introduced the symbol P _h to denote the fact that in changing the Lagrangian, the Lagrange multiplier, P _h, need not have the same value as the previous Lagrange multiplier and in fact now corresponds to the hydrostatic component of the pressure. The work done by the viscous pressure drop is accounted for by the factor of 2 multiplying τ_xx .

Although probably less familiar to most readers, upon invoking Hamilton’s principle, Equations (33) and (34) are equivalent to the vectorial form of the full Stokes equations for a non-Newtonian fluid (see Appendix A). An advantage of the formulation just presented is that the variation of each of the terms in Equations (32) and (33) is related to ice-sheet energetics providing an intuitive physical interpretation. Thus this approach has many of the advantages of the force-budget technique in terms of interpreting the underlying equations (Van der Veen, Reference Van der Veen1999).

For future use, we also note that the normal and tangent vectors to the bed and surface are given by

(35)

while those for the surface are

(36)

3. Approximations

3.1. Non-dimensionalization and the long-wavelength approximation

The usual recipe for making approximations in fluid mechanics is to first non-dimensionalize the governing vector equations and then try to simplify the resulting system of equations by neglecting small terms. The same recipe applies to the Lagrangian formulation of ice-sheet mechanics. We define characteristic ice-sheet thickness H ₀, horizontal length L ₀ and scale velocities (u, w) with the characteristic speeds (U ₀, U ₀ H ₀ /L ₀₎ such that

The tilde perched on top of symbols indicates dimensionless variables. With these definitions we can define a characteristic viscosity scale

We scale the pressure by ρgH ₀ and stresses with unit μ ₀ U ₀/L ₀. We also scale the sliding coefficients β and γ such that

The aspect ratio of the ice sheet is

There are two choices for the velocity scale. We can either choose a leading-order balance between pressure gradients and gradients in lateral stresses or a leading-order balance between the pressure gradient and gradients in vertical shear stress. We arbitrarily choose the former balance, which selects the velocity scale U ₀ = ρgH ₀ L ₀/μ ₀. This choice provides an upper bound on the velocity scale; the actual velocity within the ice sheet may actually be much smaller than order unity and this will have consequences for the terms maintained in the long-wavelength limit.

Lightening the notation by dropping the tilde decoration, the non-dimensional versions of Equations (22) and (27) are

(37)

(38)

in which μ′(θ) represents any (non-dimensional) temperature dependence of the viscosity. Using the above expressions, Equation (34) in non-dimensional form becomes

(39)

where

(40)

is a non-dimensional number that measures how easy it is to shear a column of ice relative to the underlying substrate. As we shall see, this dimensionless number, which we call the basal traction number, plays an important role in controlling the transition from laminar flow to plug flow and hence from the SIA to the SSA.

The long-wavelength approximation corresponds to the condition λ ≪ 1. Dropping terms of order λ ², the non-dimensional stresses and viscosity become, respectively,

(41)

(42)

Examining Equation (39), it is clear that

(43)

and makes a much smaller contribution to the Lagrangian than the remaining terms. After dropping the above term, we find

(44)

or, in terms of velocities,

(45)

where the basal slip velocity u _b is equal to u to order λ ². We have affixed the subscript λ to to remind ourselves that this expression is now only an approximation of the true Lagrangian, valid for features that have a length scale that is large compared to the ice thickness. We have selectively maintained some terms in the expansion that would be present if we had instead chosen a dominant balance between vertical shear stresses and the pressure gradient. See Schoof and Hindmarsh (Reference Schoof and Hindmarsh2010) for further discussion of the two-parameter expansion that this implies.

3.2. Elimination of vertical velocity

It would appear from Equation (45) that we need to solve for both u and w. However, variation with respect to the Lagrange multiplier K yields the no-penetration condition

(46)

providing a relationship between u and w along the ice-sheet bed. Now the only place that w appears is within the first integral in Equation (45):

(47)

Using the product rule followed by the divergence theorem, this becomes

(48)

where the surface integral is over the entire ice-sheet surface and bed and n_z is the z component of the outward normal vector. Using the calculus of variation, we insist that the volume integral vanish for arbitrary variation of w, from which we conclude that the (non-dimensional) hydrostatic pressure P must be

(49)

Making use of the fact that along the ice-sheet bed to consistent order this implies that P = h, w = ub_x and n_z = −1, the integral over the ice-sheet bed becomes:

(50)

Substituting back, Equation (44) becomes

(51)

Using the incompressibility condition again, the non-dimensional viscosity can also be written solely in terms of u:

(52)

whence we see that now depends exclusively on u. This is the variational formulation for the flowline equivalent of Blatter’s model (Blatter, Reference Blatter1995; Colinge and Rappaz, Reference Colinge and Rappaz1999; Glowinski and Rappaz, Reference Glowinski and Rappaz2003; Schoof, Reference Schoof2010). Our task is to find an approximation for u.

4. The Rayleigh–Ritz Approximation

Equation (52) corresponds to a variational formulation of the higher-order ‘Blatter’ model (Blatter, Reference Blatter1995). The general problem of finding an expression for the continuous function u(x, z) that renders Equation (51) stationary for a given ice-sheet geometry is formidable. Instead, we seek a physically based approximation for u using the Rayleigh–Ritz method. It is here that we depart from previous studies. The object of the Rayleigh–Ritz method is to replace the problem of finding stationary points of functionals with the (much easier) problem of finding stationary points of functions of several variables. We propose expanding the horizontal velocity u as a linear combination of known functions

The functions ψ_i satisfy the appropriate boundary conditions, but are otherwise arbitrary, and the parameters a_i are, for the moment, undetermined. By the usual methods of differential calculus, to find the stationary value of L we need to find the coefficients a_i such that

This results in a set of simultaneous (nonlinear) equations for the coefficients a_i that can be solved using standard methods. As the number of coefficients, a_i , increases, provided a suitable set of complete basis functions has been used we expect the approximation to become increasingly accurate. The key in the Rayleigh–Ritz approximation is to find a set of basis functions that approximate the ‘true’ solution with a small number of terms.

Most previous studies have focused on finding approximations for the stress within the ice and then inferring what this means for the velocity structure. The Rayleigh–Ritz method reverses the procedure and we instead try to approximate the velocity and then infer what this means for the stress within the ice. (In the end, the two approaches must be equivalent. However, sometimes one approach is easier or more intuitive than the other.) With this in mind, our goal is to find a semi-analytic approximation for the velocity u. We shall try expanding u using the partial basis set:

(53)

This expansion has the advantage that only the first term in the expansion is non-zero at the ice-sheet base. Thus, the first term accounts for the contribution of sliding to the horizontal velocity; the remaining terms in the sum account for internal deformation within the ice. It is also apparent that the first term in the above expansion corresponds to the plug flow limit of the SSA while the n + 1 th term corresponds to the SIA limit. More generally, we suspect that u(z) will be smooth and thereby well approximated by a small number of terms in a polynomial expansion. This hints that a simple approach to marrying the SIA and SSA is to assume an approximation that is a combination of SIA and SSA such that

(54)

In this approximation we assume that the ‘true’ velocity profile anywhere in the ice sheet can be approximated by a linear combination of plug flow plus laminar deformation.

For a ₀ ≠ 0 we can write the expansion as

(55)

with

(56)

We have introduced the factor of n+1 in the denominator of the second term to simplify later calculations. The variable ∊(x) is for now arbitrary, but we shall see that it corresponds to the spatially variable component of the basal traction number ∊₀. If ∊(x) ≫ 1, the velocity is, to a good approximation, laminar and the velocity is proportional to that of the SIA. Conversely, if ∊(x) ≪ 1, the velocity is nearly plug-like and we recover the SSA. As hinted at earlier, improved accuracy of this approximation can be obtained by including a larger number of terms in the expansion. In this case it may be more appropriate to use a better-conditioned set of basis functions such as Chebyshev or Jacobi polynomials (Boyd, Reference Boyd2001).

It would seem that for the above two-term Rayleigh–Ritz expansion we still need to solve for both the functions U(x) and ∊(x). However, we have the additional constraint that u must also satisfy the surface and basal boundary conditions. Furthermore, since we have already made the long-wavelength lubrication theory approximation, we need only satisfy the boundary conditions approximately to the same order as the original differential equation. As we shall show, this simplifies the dynamics such that the momentum balance is essentially a one-dimensional (1-D) problem and the dependence of u(x, z) on z can be determined analytically.

4.1. Matching the surface and basal boundary conditions

The boundary conditions that u must satisfy are buried in Equation (51), and extracting them requires some manipulation. Noting that from Equations (35) and (36) the outward normal vectors to the ice-sheet surface, n _s, and base n _b, to order λ ² are

we apply the divergence theorem to the first term in Equation (51):

The integral over the ice-sheet base that arises from the divergence theorem will cancel the equivalent and opposite signed term in the integral over the ice-sheet base in Equation (51). Since the hydrostatic pressure at the ice-sheet surface is zero, the integral over the ice-sheet surface that would also normally arise is also zero.

Applying the divergence theorem to the next group of terms

(57)

the Lagrangian becomes

(58)

where the integral over ∂Ω_s denotes the area of the ice-sheet surface.

Since the integrals over the surface and base of the ice sheet must vanish for arbitrary variations δu, application of the calculus of variations with the aid of Equation (41) shows that the horizontal velocity must satisfy the following surface and basal boundary conditions:

(59a)

(59b)

Rewriting these in terms of U and neglecting terms of λ ² and higher, the approximate boundary conditions that U must satisfy are

(60a)

(60b)

The above derivation ignores the boundary terms that arise along the ice-sheet margin. Since these terms do not play a meaningful role in the following analysis and their inclusion requires the cumbersome chore of writing an additional integral along the ice-sheet calving front with every application of the divergence theorem, we defer this calculation to Appendix B.

Substituting Equation (56) into the above two boundary conditions and using the fact that

(61)

we find that the surface boundary equation, Equation (60a), is automatically satisfied by the expansion. To satisfy the basal boundary condition, Equation (60b), ∊(x) must be given by

(62)

where μ _b ≡ μ(b) is the viscosity evaluated at the ice-sheet bed (itself a function of U) and we have defined the effective sliding coefficient.κ,

(63)

As anticipated,∊(x) represents the spatial variability of the basal traction number. With these results, the z dependence of u in Equation (56) is defined analytically, although implicitly determined, and the 2-D problem is transformed into the 1-D problem of finding an expression for U(x). The approximation for the horizontal velocity can be written as u(x; e), where the parameter ∊ now determines the vertical structure of the velocity field.

4.2. Matched approximation with analytic z dependence

We are now in a position to write down an approximation that marries the SIA and SSA. To do this we start by combining Equations (55) and (56) to form an explicit expression,

(64)

where ∊(x) is defined by Equation (62). This expression for the velocity can then be substituted into either Equation (51) or Equation (58). Substituting the expression for velocity into Equation (51) (i.e. prior to application of the divergence theorem) provides an approximation based on the weak integral form that would be appropriate for a finite-element implementation. Substituting the velocity into Equation (58) (i.e. after the divergence theorem is invoked) provides an approximation that is valid at each point within the ice sheet. In this strong form we have the additional requirement that velocity is at least twice differentiable. In the absence of discontinuities in the velocity field, these two approaches should yield equivalent (but numerically different) approximations. In practice, the numerical approximations also differ, in that for the former we often discretize U before Hamilton’s principle is invoked whereas for the latter we usually discretize U after the variation is taken.

4.2.1. Approximation 1

In the first approximation, we substitute Equation (55) into Equation (51) to find

(65)

The integral over the ice-sheet bed is now independent of z, allowing us to drop the subscript, reminding us that the integral is supposed to be performed along the ice-sheet base. We can use the calculus of variations to take the variation of Equation (65) with respect to U to find

(66)

This must hold for arbitrary variation δU. Each of the terms in Equation (66) has a straightforward interpretation in terms of the energy balance of the ice sheet: (1) the first term represents the rate of work done against gravity by redistributing mass vertically; (2) the second and third terms represent the viscous energy dissipation; (3) the fourth term represents the energy dissipation by basal friction; and (4) the fifth term represents the rate of work done to push a column of ice up (or down) a basal slope.

If we discretize U using a (finite) set of basis functions (splines, piecewise polynomials, etc.) then we arrive at the form of conservation of momentum that would be used in a Galerkin finite or spectral element model (FEM). In finite-element terminology the variation δU is called a weighting function or test function and is arbitrary (Reference HughesHughes, 1987). If ψ = 1, then we can perform the required integration over the ice-sheet thickness to find

(67)

with depth-averaged viscosity denoted by . This is the flowline equivalent to the original finite-element formulation for the SSA (Reference MacAyealMacAyeal, 1989; see Equation (B2) in Appendix B).

For more general ψ ≠ 1, Equation (66) represents a generalization of the SSA to include vertical shear stresses. Unfortunately, because the viscosity is also a nonlinear function of U, analytic vertical integration is not usually possible and we need to resort to a numerical quadrature scheme. This difficulty notwithstanding, any finite-element formulation designed to solve the SSA should be easily modified with a subroutine that performs the vertical quadrature. This is the only change required to augment a SSA model to include an improved treatment of vertical shear stresses.

4.2.2. Approximation 2

As an alternative starting point, we can begin with Equation (58). Again, using the calculus of variations to take the variation with respect to δU we have

(68)

where we have used Equations (55), and (60a) and (60b). For this to vanish for arbitrary δU implies

(69)

After integrating over the ice thickness, again with the aid of Equations (60a) and (60b), we find the Euler–Lagrange equation corresponding to conservation of momentum for the two-term matched approximation:

(70)

With the exception of the z-dependent weighting factor ψ, the approximation is equivalent to the SSA (Reference MacAyealMacAyeal, 1989). The depth dependence of ψ has the effect of down-weighting the magnitude of lateral stresses. This is because a fraction of the gravitational potential energy is being dissipated by vertical shearing within a column of ice as well as lateral stretching; the larger the value of ∊ ₀, the larger the fraction of gravitational potential energy that is dissipated by vertical shearing.

Equation (70) also reveals an important limit of our approximation: we need to differentiate the basal sliding parameter β twice. This suggests some care must be taken in applying Equation (70) to regions where there is a rapid or stepwise change in β, say in the transition from grounded to floating ice, because this will lead to large horizontal gradients in ψ. In this case, neither of the approximations that we spliced together is valid. Thus a reasonable expectation is that the approximation is appropriate for situations where basal properties vary over a length scale that is large compared to the ice thickness. The alternative form, Equation (67), is not immune to this weakness. However, Equation (67) only requires that β ² be differentiated once. For this reason, the weak form may have several advantages over the strong form when there are significant variations in basal properties.

To summarize, having started with an expansion for u that combines the shallow-ice and shallow-stream approximations, we have found an approximation that stitches these two limits together smoothly using the basal boundary conditions. We show in Appendix C that Equation (70) analytically reduces to the SIA or SSA depending on the size of the quantity ∊ ₀ λ ⁻².

5. Numerical Examples

To provide a more concrete illustration of the approximation developed in the previous section, we seek steady-state solutions of Equation (70) for a range of basal friction parameters. This requires solving the nonlinear differential equation along with the continuity equation

(71)

along with the continuity equation

(72)

In Equation (72), is the accumulation/ablation rate and is the depth-averaged horizontal velocity. At x = 0 we assume an ice-divide symmetry boundary condition such that

At the (fixed) ice-sheet margin, x = L, we specify a large, negative, accumulation rate, as would be true if there were a very large calving rate at the margin. This physical condition forces the ice thickness at the ice-sheet margin to be zero without the need to (over-)specify an additional boundary condition for Equation (72). The velocity boundary condition that we use is the calving-front boundary condition discussed in Appendix B:

where D is the depth of the water. For those examples where D is zero we have (see Appendix B)

To facilitate comparison with analytical steady-state ice-sheet profiles, we take the bed to be flat with constant accumulation rate .

Equation (71) was discretized using a staggered grid for U and h. Centered finite differences were used in the interior with second-order, one-sided differences at the end points. We used the trapezoidal rule with ten evenly spaced points for the vertical integration. We experimented with using a larger number of nodes, but found no difference in any of the results using either 50 or 100 nodes for the vertical integration, suggesting error is dominated by the horizontal discretization. The continuity equation was discretized using a third-order upwind stencil with an implicit backward Euler step for the time derivative. We solved the nonlinear system of equations using Newton–Kantorovich iteration (e.g. Boyd, Reference Boyd2001, p.526–530) with a relative tolerance limit of 0.01%.

To find steady states we marched the continuity equation forward in time until the percent change in ice-sheet volume was <0.01% over 100 years. The numerical algorithm was probed for errors by checking to make sure we could reproduce analytic steady-state ice-shelf profiles for positive, negative and constant accumulation rates and analytic steady-state SIA profiles for constant accumulation and no basal sliding. All the following results were computed using a 100 km long ice sheet with grid spacing, Δx, of 5 km. We could, alternatively (or additionally), use Equation (65) as the basis of a finite-element implementation. None of the results depend greatly on the numerical scheme. Parameters used for all calculations are shown in Table 1.

Table 1. A list of the parameter values used for the numerical experiments

Figures 2 and 3 show the steady-state profiles and velocity contours for four different values of the β ² coefficient for linear (p = 1) and pseudo-plastic (p = 1/2, γ = 10⁻⁸ m s⁻¹) sliding laws, respectively. The analytic, SIA steady-state profile is also shown as a dashed curve (Paterson, Reference Paterson1994, p.243, equation 10). For a linear sliding law with β = 10¹⁴ Pa s m⁻¹, the numerical steady-state solution is indistinguishable from the analytic solution. As β ² decreases, the velocities increase and the steady-state profile begins to diverge from the SIA. For β = 10⁹ Pa s m⁻¹, the velocity profile has transitioned completely to plug flow with a velocity profile that is uniform with depth. The steady-state profiles with a pseudo-plastic sliding law show the same qualitative pattern as the Newtonian sliding law. High values of β ² are indistinguishable from the SIA steady-state solution. As β ² decreases, increased sliding results in a transition to a plug-like velocity profile. Both Newtonian and pseudo-plastic sliding laws show plug-like flow near the ice divide, where the strain rates are very small. (Note that the units of β ² for a linear versus plastic sliding law are different. Thus the profiles for linear and plastic profiles are not directly comparable.)

Fig. 2. Steady-state ice-sheet profiles and velocities for a linear, Newtonian sliding law (p = 1) and four different values of β ². The dashed curve shows the analytic, SIA steady-state ice-sheet profile. In the top panel, the numerical steady-state solution is indistinguishable from the analytic SIA solution. As β ² decreases, the velocities increase and the steady-state profile begins to diverge from the SIA. In the bottom panel, the velocity profile has transitioned completely to a profile where the velocity is uniform with depth.

Fig. 3. Steady-state ice-sheet profiles and velocities for a plastic sliding law (p = 1/2, γ= 10⁻⁸ m s⁻¹) and four different values of β ². The dashed curve shows the analytic, SIA steady-state ice-sheet profile. We seethe same qualitative pattern as in Figure 2 where high values of β ² lead to steady-state profiles that are indistinguishable from the SIA steady-state solution. As β ² decreases, increased sliding results in a transition to a more plug-flow-like velocity profile. Note that the units of β ² have changed relative to the Newtonian sliding law shown in Figure 2.

To better illustrate the transition from a strong bed (large β ²) to a weak bed (small β ²), we imposed a spatially varying sliding coefficient for a linear (p = 1) sliding law of the form,

(73)

β ₁ represents the strong bed, β ₂ the weak bed, and the position and width of the transition is determined by x ₀ and ℓ ₁, respectively. This is similar to the study of marine ice-sheet transition zones studied by Pattyn and others (Reference Pattyn, Huyghe, De Brabander and De Smedt2006) using a higher-order model. Figure 4a shows an example where β ₁ = 10¹⁶ Pa s m⁻¹, β ₂ = 10¹⁰ Pa s m⁻¹, x ₀ = 50 km and ℓ ₁ = 5 km. This places the transition at the center of the ice-sheet domain. (Note that, in contrast to Figure 2, the velocity is contoured with linear rather than logarithmic increments.) In contrast to previous profiles, the ice sheet now has a low surface slope and convex ice-sheet-like profile in the interior, but switches to a convex ice-streamlike profile near the margins. The shape is reminiscent of an ice sheet with a lobe that forms where the ice flows out from a region where the bed is hard crystalline rock into a region of soft till, as perhaps occurred with the Laurentide ice sheet (e.g. the Michigan lobe).

Fig. 4. The transition from a strong to a weak bed for two different cases, each with a linear (p = 1) sliding law. (a) Example of a steady-state, grounded ice sheet with an imposed change in sliding coefficient from a strong to a weak bed as defined by Equation (73). Note the change in convexity as the ice-sheet transition from a strong to a weak bed generates an ice-sheet profile with a lobe-like feature. (b) Example of a steady-state marine ice sheet, grounded 250 m below sea level, with an attached ice shelf. The sliding law used for this experiment assumes that the basal traction is tapered near the grounding line according to Equation (74).

Figure 4b shows an example of a steady-state marine ice sheet with an ice shelf (again with linear, p = 1, sliding law). The ice-sheet bed is grounded 250 m below sea level. Because the basal sliding coefficient enters into Equation (70) in differentiated form, we insist that β ² be continuous across the floating–grounded transition. To avoid a discontinuous basal traction, we taper the sliding coefficient near the grounding line such that

(74)

For the calculations shown here, β ₀ = 10¹⁶ Pa s m⁻¹ is the friction coefficient away from the grounding line, ℓ ₂ = 9 km is the width of the transition zone and x _g is the grounding line location, which moves as determined by buoyancy. Imposing a transition zone smooths the break in slope that would otherwise be present at the transition from grounded to floating ice. Tapering the sliding coefficient leads to more stable numerical results and would be physical if, for example, pore pressure increased substantially towards the grounding line. Despite the cavalier treatment of grounding line migration, the model is at least able to reproduce an ice sheet with an ice shelf attached. For the steady states computed, we used the steady-state profile shown in Figure 4a as the initial condition. As pointed out by others, the position of the grounding line may be sensitive to our numerical scheme (Vieli and Payne, Reference Vieli and Payne2005).

6. Discussion

Formally, the SIA and SSA are long-wavelength, asymptotic limits that correspond to small aspect ratios and large or small values of the basal traction number respectively. All models, ours included, must reduce to these two special cases in the limit of large wavelengths and appropriate values of basal shear stresses. Where different approximations and models (analytic and numeric) differ is in the region where (1) basal shear stress is intermediate to these two limits and (2) wavelengths become less than or equal to the ice-sheet thickness. In the long-wavelength limit, we need only worry about the transition from SIA to SSA.

A convenient measure of the degree to which the flow is laminar or plug-like is the slip ratio, defined as the ratio of the basal velocity to the surface velocity. A slip ratio close to 1 indicates that the velocity profile is nearly constant with depth (plug-like, SSA limit). In contrast, a value of the slip ratio close to 0 indicates that the ice is stuck to its bed and the velocity–depth profile is nearly laminar (SIA limit). Figure 5 shows the non-dimensional slip ratio plotted against the non-dimensional basal traction number for each point in the steady-state profiles shown in Figures 2–4. Since both the slip ratio and basal traction number are dimensionless, this curve is independent of ice-sheet geometry, temperature, basal sliding law, etc. Notice that a smooth curve connects regions with slip ratios that asymptotically approach a value of 1 at low values of ∊(x) to slip ratios that asymptotically approach 0 at high values of ∊(x). Different models will approximate the ‘true’ transition curve to greater (or lesser) accuracy. Mixture models represent the transition as a step function with an abrupt, instantaneous step from a slip ratio of 1 to a slip ratio of 0. This leads us to suggest that curves like Figure 2 and their generalization to include the ratio of wavelength to ice thickness can be beneficial for model intercomparisons. The suite of experiments of, say, the Ice-Sheet Model Intercomparison Project (ISMIP; Pattyn and others, Reference Pattyn2008) could be simplified by comparing the scaling laws implied by transition curves such as Figure 2. This may provide a better metric to quantify model differences than those currently used.

Fig. 5. Scaling of the slip ratio (the ratio of the basal to surface velocities) as a function of the non-dimensional basal traction number. A slip ratio close to 1 indicates that the ice is freely slipping over the bed with little vertical shear, whereas a slip ratio close to 0 indicates that the majority of the deformation is due to vertical shearing within a column of ice. We see that small values of basal traction number asymptotically approach a slip value of 1, corresponding to the shallow-stream approximation. In contrast, large values of the basal traction number asymptotically approach a slip value of 0, corresponding to the SIA. The two limiting behaviors are joined by a smooth curve.

At this point, we return to an earlier conjecture, in which we postulated that the primary limitation of ice-sheet models was not the ability to resolve 3-D stresses, but our uncertainty in basal parameters. From Figure 2 it is clear that an uncertainty in the basal traction number (itself a function of ice thickness, sliding coefficient and viscosity) will lead to an uncertainty in position on the transition curve. Calculating the transition curve to a much greater accuracy than the uncertainty in the basal traction number is then unlikely to increase model accuracy. We view the approximation presented here and the analogous approximations used by others as complementary to, rather than competing with, higher-order and full Stokes models, such as those proposed by Blatter (Reference Blatter1995) and Pattyn (Reference Pattyn2003). For many modeling applications, the full Stokes and/or higher-order approach are too computationally intensive to use for pal eo-climate/ice-sheet modeling simulations. For these long-timescale simulations, an intermediate complexity approximation such as the one presented may be an acceptable compromise, especially since the accuracy of the model is severely limited by poor knowledge of basal hydrology and /or basal sliding.

The approximations presented here (Equations (70) and (66)) are similar to the heuristic approximation developed by Kamb and Echelmeyer (Reference Kamb and Echelmeyer1986) and used by, amongst others, Hindmarsh (Reference Hindmarsh2004) and Reference Pollard, DeConto, Hambrey, Christoffersen, Glasser and HubbardPollard and DeConto (2005). Both approximations reduce to the SIA and SSA under appropriate limits. The two approaches, however, differ in two significant ways. Our approximation uses a single nonlinear partial differential equation for the velocity. In contrast, Reference Pollard, DeConto, Hambrey, Christoffersen, Glasser and HubbardPollard and DeConto (2005) iterate on both the viscosity and the solution of two different stress-balance equations. Whether one method is more accurate or advantageous is a question that requires further study and we do not pursue it further here. The advantage of our approach is that any ice-sheet model designed to solve the SSA can easily be generalized to include vertical shear stresses by adding a depth integration subroutine to an existing SSA model. The accuracy of our approximation can be improved by adding additional terms to the series expansion for u. In essence, we propose that higher-order and/or full Stokes models should consider using a spectral expansion in the vertical direction. Since the velocity profile is usually very smooth, a small number of terms are necessary. The two-term approximation that is the focus of this study is a special case of this more general expansion.

7. Conclusions

We have shown that a modified version of Hamilton’s principle can be used to systematically develop variational principles for the mechanics of slow-moving ice sheets, provided we can neglect the kinetic energy. This approach has the advantage that each of the terms in the equation of motion has a well-defined meaning in terms of ice-sheet energetics. Using this approach, we have developed an approximation, based on a two-term Rayleigh–Ritz expansion of the velocity field, that gently and continuously marries the shallow-ice approximation to the shallow-stream approximation. This approximation transitions smoothly between the two end-member cases of the SIA and SSA and analytically reduces to the SIA and SSA solutions, respectively, depending on whether a non-dimensional number called the basal traction number is large or small. Numerical examples show that the steady-state ice-sheet profiles generated by our approximation look realistic.

A separate question is whether Hamilton’s principle is a useful alternative to the vectorial formulation of mechanics. Although the approximation we developed could be derived (possibly in fewer steps) using more conventional methods, it illustrates the potential that Hamilton’s principle provides to find new approximations and the ease with which one can find both the variational principle and the corresponding Euler–Lagrange equation from the variational principle. But it is not necessary to use Hamilton’s principle and we could have obtained the same results starting from the vectorial formulation of mechanics. Although Hamilton’s principle appears to provide nothing new and merely regurgitates the well-known Stokes equations, the advantage that it provides is that we can systematically construct variational principles from the underlying energetics. Thus the elegant theoretical work done by Schoof (Reference Schoof2006a,Reference Schoofb) in deriving a variational principle for the shallow-stream approximation can be viewed as a special case of Hamilton’s principle. A final advantage, which we have not pursued here, is that the non-Newtonian power-law rheology of ice permits complementary variational principles (Finlayson, Reference Finlayson1972) which may be useful in finding strong upper and lower bounds on the response of ice sheets to perturbations. The ability to systematically bound ice-sheet behavior may be the single greatest advantage of Hamilton’s principle. Nonetheless, as we alluded to earlier, the appreciation of Hamilton’s principle may, in the end, primarily be a matter of taste.