Interactive comment on “ Reformulating the full-Stokes ice sheet model for a more efficient computational solution ”

The manuscript presents a new formulation of the Stokes equations that includes only two horizontal velocities as independent variables. In addition to reduced number of independent variables (two instead of original four), this formulation is positive-definite that makes it more attractive from a numerical point of view. This formulation definitely has a merit, however, the manuscript needs substantial revisions in both presentation and the equation formulation before it could be published.


Introduction
The most general and accurate model currently used for the simulation of ice sheet dynamics is based on non-Newtonian Stokes flow (e.g., Greve and Blatter, 2009).At present, however, a full-Stokes model presents formidable challenges for large-scale modeling, although such models exist and are being used (e.g., Zwinger and Moore, 2009, implemented in the ELMER (http://www.csc.fi/english/pages/elmer)code package).As a consequence, there is considerable interest in various approximate models (e.g., the first order or Blatter-Pattyn approximation, and the shallow ice and shallow shelf approximations) that are more limited but computationally far cheaper (e.g., Pattyn et al., 2008).
Typically, a discretized Stokes model may be written in matrix form as where A = A T is a square, symmetric, positive-definite matrix representing the negative of the discrete nonlinear stress divergence operator in the momentum equations, u i is a vector of three-dimensional velocities, P is the pressure, G is the discrete gradient operator, and G T is the negative of the discrete divergence operator.The right hand side contains contributions from gravitational forces and boundary conditions.A matrix system in the form of Eq. ( 1) is known as a saddle point problem that typically arises, as in this case, from a constrained optimization problem.The system matrix on the left hand side of Eq. ( 1) is symmetric but indefinite, meaning that its eigenvalues are real but have both positive and negative values.There are three main difficulties in the solution of such problems.First, the matrix problem is quite large, with matrix rank of order 4N, where N is the number of cells in the mesh, i.e., three velocity components and the pressure.Second, large-scale saddle point problems are typically solved iteratively using Krylov subspace methods (conjugate gradient-type algorithms).Such methods tend to converge slowly and are prone to failure when applied to saddle point problems, so it is necessary to find and apply a good preconditioner to achieve reasonable convergence.In fact, there is a voluminous literature on appropriate methods for the numerical solution of saddle point problems (see Benzi et al. (2005), for example).Finally, in the finite element context, basis functions for the pressure and velocity have to be chosen carefully so that the discrete problem is well posed (this involves satisfying the so-called Brezzi-Babuska or inf-sup condition, see Brezzi and Fortin (1991) or Elman et al. (2005), for example).
In glaciology, these difficulties have typically been avoided by the use of an approximate Stokes model, the socalled first-order model, otherwise called the Blatter-Pattyn model, first introduced by Blatter (1995) and refined by Pattyn (2003).The Blatter-Pattyn model is obtained by invoking the small aspect ratio approximation, i.e., assuming that the ratio of the characteristic vertical and horizontal length scales in the ice sheet velocity field is small, thus neglecting the mixed horizontal-vertical stress tensor components.As a result, it becomes possible to vertically integrate the vertical momentum equation and the continuity equation to obtain pressure as a function of the vertical velocity, P = P (w), and the vertical velocity as a function of the horizontal velocity components, w = w(u (i) ) (see Pattyn, 2003).This allows the elimination of both the pressure and vertical velocity from the approximated Stokes model to obtain a reduced system in terms of the horizontal velocity components only, which may be expressed in matrix form as follows where the index (i) represents just the horizontal components, and Ã is a symmetric, positive-definite matrix of reduced rank (of order 2N) as compared to the system matrix in Eq. (1).In contrast to Eq. (1), the system corresponding to Eq. ( 2) is associated with the minimization of a positivedefinite functional and is therefore ideally suited to solution by Krylov subspace methods (Knoll and Keyes, 2004) or even by direct numerical optimization methods (Nocedal and Wright, 2006).As a result, the Blatter-Pattyn system is much easier to solve than the full Stokes system.However, the Blatter-Pattyn model is much more limited in application (e.g., see the discussion and results in Pattyn et al., 2008).This is because of the small aspect ratio approximation and, in addition, because of a further approximation implicitly built into the Blatter-Pattyn model, limiting it to small basal slopes (see Dukowicz et al., 2011, henceforth referred to as DPL11).To partially remedy this last problem, DPL11 introduced a second Lagrange multiplier, , to enforce tangential flow at the base.
In the present paper we make the observation that there is no need for Lagrange multipliers P and if one already has a velocity field that satisfies both continuity and the basal no-penetration boundary condition for use as a trial function in a variational formulation, in loose analogy with the Ritz method.We note that such a velocity field is available, at least in principle, from vertically integrating the continuity equation to obtain the vertical velocity in terms of horizontal velocities, w u (i) , as is done in the Blatter-Pattyn model.In Dukowicz et al. (2010) (henceforth referred to as DPL10) and in DPL11 it was shown that non-Newtonian Stokes flow, including boundary conditions, may be expressed as a constrained variational principle expressed in terms of an action functional, A S [u i ,P , ], whose arguments represent the functions with respect to which a stationary point is to be found.Eliminating the vertical velocity from the Stokes action, we obtain the "reformulated" Stokes action, as follows which, together with w = w u (i) forms a complete specification of the Stokes problem.Note the following properties: (a) the action A RS u (i) is exactly equivalent to the Stokes action, as indicated in Eq. ( 3), (b) since both Lagrange multipliers are zero, the reformulated action is positive-definite, just as in the Blatter-Pattyn model, and (c) this action leads to a matrix system of exactly the same form as Eq. ( 2).The resulting matrix system, therefore, has exactly the same beneficial properties as the Blatter-Pattyn system, except now without approximation.
It is interesting to note that Pattyn (2008) presents a reformulation of the full Stokes model that superficially resembles the Blatter-Pattyn model.However, this reformulation basically amounts to expressing the pressure P in terms of an alternative variable, the vertical stress component τ zz , and thus leads to an iteration scheme that is effectively equivalent to the solution of a system in the form of Eq. (1).One feature of the present reformulation is that it, in effect, leads to an integro-differential formulation of the Stokes problem.Integro-differential formulations have appeared previously in the glaciological literature (e.g., Van der Veen and Whillans, 1989;Hindmarsh, 1993).These early formulations appear to involve the vertical integral of the shear stress τ xz and are therefore different in substance and motivation from the present formulation.
In the remainder of the paper we review the variational formulation of the basic Stokes problem in Sect.2, in terms of an "action" functional, making the simplifying assumption that the ice sheet is in contact with and sliding along a rigid, fixed bed, as in DPL11 and elsewhere in the literature.As pointed out in DPL10 and DPL11, the action in a variational formulation completely determines the problem and is in fact the preferred starting point for a discretization of the problem.In Sect. 3 we generalize the basal boundary condition to allow for a moving basal surface and the possibility of mass flux across the surface, as at the base of a floating ice shelf.In Sect. 4 we obtain the reformulated Stokes action in two different versions.In Sect. 5 we illustrate the implementation of the present method, and thereby provide some justification for the claims of computational efficiency, by means of a relatively simple but nontrivial test problem involving the sliding of an ice sheet along an inclined plane.This test problem is particularly attractive because it also provides an analytical solution, which helps in the understanding of the present method and makes it easy to check the validity and accuracy of the numerical solution.For completeness, in Sect.6 we also obtain the corresponding reformulated Euler-Lagrange partial differential equations and boundary conditions.These equations may be of interest for comparison with the full Stokes system of equations, and possibly may also suggest other approximations, perhaps more accurate, to the Stokes model.Finally, in Sect.7 we summarize and draw some conclusions.

The basic Stokes model
We begin with the variational principle for the non-Newtonian ice sheet Stokes model whose action functional (see DPL11) is given by where u i ∈ {u,v,w} is the three dimensional velocity vector, g i is the gravitational acceleration vector (typically g i = (0,0,−g)), ρ is the ice density, assumed constant, and ε2 = εij εij is the second invariant of the full Stokes strain-rate tensor, where We define where is the Glen's law viscosity coefficient, typically used with exponent n = 3, and µ 0 (θ) is a temperature-dependent coefficient.As in DPL10 and DPL11, we illustrate the effect of basal stress forces by j (u) = −βu i u i n j 2, which represents a linear frictional sliding law with a constant coefficient, β; β ≥ 0. However, other frictional laws are easily accommodated as in Schoof (2010), for example.The three integrals above cover the entire ice sheet volume and its upper and basal surfaces, respectively.Here, x i ∈ {x,y,z} is the position vector, and n i is the outward-pointing unit vector at the ice sheet bounding surfaces.Note that Cartesian tensor notation is being used, and, where appropriate, the summation convention on repeated indices.In general, tensor indices are three-dimensional, i.e., i,j,••• ∈ {x,y,z}, except when an index appears in parentheses, in which case it denotes an index in the horizontal plane only, e.g., (i),(j ),••• ∈ {x,y}, so that, for example, u (i) ∈ {u,v}, u i u i = u 2 + v 2 + w 2 and u (i) u (i) = u 2 + v 2 .As mentioned previously, the functional, Eq. ( 4), represents a constrained minimization principle, with constraints enforced by three Lagrange multipliers, P , s , and b .The pressure P enforces incompressibility, and b enforces tangential flow at the base.In spite of its role as a Lagrange multiplier, pressure also has a physical role in the presence of gravity.For instance, in the case of static flow (maintained by confining walls, for example), the pressure satisfies a hydrostatic balance equation and therefore we need an upper surface boundary condition, P = P s , where P s is some known or somehow specified pressure at the upper surface.Very frequently, we have P s = 0 if there is no ice or water weighing down from above.In the general case, the pressure will also require a separate boundary condition at the upper surface, and P = 0 is appropriate if atmospheric pressure is negligible.This condition is enforced by s .Alternatively, in many cases it may be simpler to directly insert such a boundary condition into the matrix equation.
The variational principle states that the solution of this dynamical system in terms of the arguments, i.e., the velocity components u i , pressure P , and Lagrange multipliers s , b , is to be found at the stationary point of the action, obtained by setting the functional derivatives with respect to the arguments equal to zero, as follows This yields the following Euler-Lagrange equations: (a) A three-dimensional momentum equation, where σ ij = τ ij − P δ ij is the Cauchy stress tensor and τ ij = 2µ n (ε 2 )ε ij is the deviatoric stress tensor, and (b) the continuity equation for incompressible flow, In addition, the following boundary conditions are implied.At the upper surface S (s) , specified at any instant of time by z = z s (x,y,t), we have (c) stress-free boundary conditions and, using (10), we have deduced that A s = 0. Along a fixed basal surface S (b) , specified by z = z b (x,y), we have (d) frictional tangential sliding boundary conditions The variational principle actually yields an equation containing b − P , but this is easily eliminated, as in DPL11, to obtain Eq. ( 14).The unit normal vectors that appear here are defined as follows For clarity, we employ superscripts (s) and (b), and subscripts s and b, to indicate an upper surface or basal value, respectively, particularly in those cases where confusion is possible.For concreteness, we have assumed a simplified ice sheet configuration illustrated in Fig. 1 that is subject to boundary conditions Eqs. (11-14), namely, an upper surface entirely exposed to the atmosphere and a basal surface that is entirely in contact with and sliding along a rigid bed.Further, for the purpose of this paper we implicitly define the upper surface by the condition n (s) z > 0, and the basal surface by n (b) z < 0. We have chosen to use this commonly used configuration since there are great many possibilities and it is impossible to deal with them all.The Stokes model itself is of course entirely general.In the next Section we shall indicate how to generalize to a moving and possibly melting basal surface, as at the base of a floating ice shelf.

Generalizing the basal boundary condition
So far we have assumed a fixed and rigid basal surface specified by z = z b (x,y).In such a case the no-penetration condition, Eq. ( 13), is given by More generally, for a moving material surface (i.e., a Lagrangian surface with no inflowing or outflowing flux due to a gain or loss of mass crossing the surface) and specified by z = z b (x,y,t), we have In addition, assuming an outwardly directed flux of mass at the basal surface with a normal velocity of magnitude u n , which may be due to melting, ablation, etc., we obtain where the effective basal vertical velocity due to both the motion of the interface and an outflowing mass flux.In general, and in particular at the base of a floating ice shelf, we might expect that w (b) n = 0.For our present purpose we assume that it is a given quantity.In general, however, the velocity w (b) n is unknown and must be determined by the simultaneous solution of the ice sheet problem and the external environment.
Integrating the continuity equation, Eq. ( 10), in the vertical direction with Eq. ( 19) as the boundary condition, the vertical velocity is given by or, alternatively, using Leibniz's theorem, one obtains Either one of these expressions corresponds to the relation w = w u (i) referred to earlier.The choice between them will depend on which is preferable from the point of view of discretization.
We note that w (b) n will vanish along certain sections of the ice sheet basal surface (i.e., when the ice sheet is sliding in contact with a fixed and rigid bed) but may have nonzero values elsewhere.It is therefore to be considered as a general function of the horizontal position vector x (i) over the entire basal surface.Similarly, the friction coefficient β may be considered as a function of horizontal position over the entire basal surface, vanishing when the ice sheet is no longer in contact with the bed.This way, the basal surface integral in Eq. ( 4) may be extended over the entire basal surface without loss of generality.However, in general we might expect that β is zero when w (b) n is nonzero and vice versa, so that βw (b) n = 0.In the following, we shall assume this to be true, while leaving open the possibility of exceptions under unusual circumstances.

The reformulated action principle
As discussed in Sect. 1, the Lagrangian multipliers P and b are no longer needed if the vertical velocity given by Eqs. ( 20) or ( 21) is used in the action functional.This is because the three-dimensional velocity field, given by the horizontal velocity components and the vertical velocity from Eqs. ( 20) or ( 21), already satisfies the continuity equation, Eq. ( 10), and the correct basal boundary condition, Eq. ( 13).Furthermore, eliminating P also removes the need for s .Substituting this velocity field into Eqs.( 4) and ( 5), the variational principle now becomes a function of horizontal velocity only, where and w u (i) is given by either Eq. ( 20) or ( 21).The subscript RS stands for "Reformulated Stokes".Observe that ε2 RS is actually the same as ε2 since the velocity fields in the Stokes and reformulated Stokes cases are the same.In general, the term involving w (b) n vanishes in Eq. ( 24) because of our assumption that βw (b) n = 0.The action, Eq. ( 22), may be simplified somewhat.As shown in Appendix A, the gravitational work term in Eq. ( 22) is expressible as follows The last term on the right hand side is independent of u (i) ; as such, it does not participate in the variational principle and may be omitted.Substituting this into Eq.( 22), the action takes the following alternative form, It may be observed that both functionals (excluding gravitational terms since they are only responsible for the forcing) are positive-definite, in contrast to the standard Stokes functional.Therefore, the variational principle is now a true minimization problem subject to gravitational forcing, just as in the Blatter-Pattyn approximate model.Also, as noted before, this is a fully three-dimensional problem in only two variables, i.e., the two horizontal velocity components, again as in the Blatter-Pattyn model.Furthermore, all boundary conditions are automatically and correctly incorporated, including the basal no-penetration (or tangential flow) boundary condition.Note that these functionals are to be used jointly with Eq. ( 20) or ( 21), as emphasized in Eqs.(22)(23)(24), to obtain the complete three-dimensional velocity field.
This action, Eq. ( 26), (or alternatively, Eq. ( 22)) is the preferred starting point for a numerical solution of the problem.This is because the discretization of the variational principle applied to the action automatically yields a symmetric, positive-definite matrix problem, analogous to Eq. (2), which is optimal for an efficient numerical solution, as discussed earlier.However, one possible disadvantage of this reformulation is that the action contains higher order derivatives than in the standard case, which may impose additional continuity requirements on the approximating space.A discussion of the requirements for the approximation space, however, is beyond the scope of the present paper.

A simple test problem
In order to better understand the reformulated Stokes problem and to partially justify the claims of improved efficiency we consider a simple two-dimensional test problem that deals with the sliding of a slab of ice of uniform thickness, characterized by a constant viscosity µ, on a basal plane inclined at an angle θ, as illustrated in Fig. 2. The problem is very attractive because it provides a nontrivial analytic Stokes flow solution that is relevant to ice sheets.Moreover, the problem can be reduced from a two dimensional configuration to a one-dimensional problem by rotating the coordinate system counterclockwise by an angle θ, i.e., (x, z) → x ,z , to align it with the ice slab.In this case all variables will be functions solely of z since the problem is longitudinally isotropic and extends to infinity in lateral directions.In addition, we have velocity v = 0 since there is no forcing in the transverse direction.In spite of its simplicity and linearity, this problem is nevertheless still useful as a means of evaluating the computational properties of the reformulated model in comparison to the standard formulation.
The required coordinate transformation is given by and therefore ice velocities transform as follows Since ∂ ∂x = 0, we obtain Note that the basal surface is located at z = 0, the upper surface at z = H , and the unit normal vectors are given by n (s) (30)

The analytic Stokes flow solution
Let us nondimensionalize by introducing a velocity scale ρgH 2 µ, a length scale H , and a pressure scale ρgH .As a result, the problem is characterized by only two independent nondimensional parameters, the angle of inclination θ, and a basal friction parameter η = βH µ.In this section, therefore, we shall consider all variables to be nondimensional.In transformed coordinates, the nondimensional Stokes system of momentum equations, Eq. ( 9), becomes while the continuity equation, Eq. ( 10), is given by These equations may be combined to obtain a separate equation for each of the three variables, as follows Note that Eq. (34) represents hydrostatic balance; this again reinforces the conclusion that pressure requires a separate boundary condition at the surface.Boundary conditions at the stress-free surface z = 1 , Eqs. ( 11) and ( 12), are given by where for simplicity we have assumed that P s = 0. Actually, in this simple test problem the pressure boundary condition is superfluous since, making use of the continuity equation, Eq. ( 32), the last two equations of Eq. ( 35) already imply P = 0. Thus, simplifying Eq. ( 35), the upper surface boundary conditions become Similarly, the basal surface z = 0 boundary conditions, Eqs. ( 13), ( 14), become sinθ u − cosθ w = 0, These are not independent.Simplifying, we therefore obtain the remaining two basal boundary conditions, as follows Finally, the system consisting of Eqs. ( 33), ( 34), (36), and (38) may be solved to obtain This solution represents ice flowing parallel to the base with an upper surface velocity of magnitude (η + 2)sinθ 2η and a basal velocity of magnitude sinθ η.Since the velocity magnitude is proportional to sinθ, the ice ceases to flow when the slab is horizontal, and conversely, velocity reaches its maximum value when the slab is oriented vertically.In the absence of friction (η → 0) there is nothing to oppose gravity and velocity becomes infinite, while for an infinite friction parameter (η → ∞) surface velocity goes to sinθ 2 and basal velocity goes to zero.It is convenient to define a basal slip parameter γ as the ratio of the basal velocity to the upper surface velocity, i.e., γ = 2 (2 + η).Thus, in the absence of friction there is 100 % slip (γ = 1), while in the limit of very large friction there is no slip (γ = 0).

Variational formulations for the Stokes and the reformulated Stokes test problems
In the present case, the Stokes action, Eq. ( 4), per unit crosssectional area, may be written as follows This incorporates all boundary conditions.As noted above, since the surface pressure boundary condition is superfluous we might simply set s = 0, although this is not necessary.The corresponding action principle leads to a onedimensional Euler-Lagrange system of equations and boundary conditions, i.e., an ordinary differential equation boundary value problem for the three variables, u,w,P , and the basal boundary constant, b , that is entirely equivalent to the system of Eqs. ( 31), ( 32), ( 35), (37).
In a similar manner, the reformulated Stokes action, Eq. ( 22), becomes where corresponds to Eq. ( 21), although Eq. ( 20) could have been used.Note that in the present case we have w (b) = 0. Eq. ( 42) is retained in this form as a reminder that this is what is discretized in the general, multidimensional case.However, in the present case it is more convenient write it in its equivalent form, Substituting this in Eq. ( 41), we obtain The variation of action, Eq. ( 44), results in a single Euler-Lagrange equation, and boundary conditions, This agrees with the horizontal part of the Stokes system from Eqs. ( 33), (36), and ( 38), and therefore leads to exactly the same solution as given in Eq. ( 39).

Discretization of the standard Stokes action
Let us now introduce a uniform one-dimensional grid with cell width h = 1 N, where N is the number of cells.The cells are indexed by k ∈ {1,2,•••,N}, and cell nodes by k ∈ {0,1,2,•••,N}, such that z k = kh.Thus, cell k is bounded by node k on the right and node k − 1 on the left.We assume that discrete velocity values are located at nodes, resulting in a piecewise linear velocity distribution, as follows and similarly for the vertical velocity component w.Noting that pressure is specified at the upper surface, it is convenient to also assume a piecewise linear distribution for the pressure, analogous to Eq. ( 47), as follows and specifically set P N = 0, or else enforce this condition by using the Lagrange multiplier s .We have already noted that in general there are subtle issues in connection with the choice of basis functions for pressure and velocity in saddle point problems, and in the Stokes system in particular.In the present case, if we use the pressure distribution, Eq. ( 48), and do not set P N = 0, we obtain a singular problem.On the other hand, if instead we use a piecewise constant pressure distribution, then the system is well behaved whether we set the surface pressure to zero or not.Substituting Eqs. ( 47) and (48) into Eq.( 40), the discretized action becomes 49) The variational principle applied to Eq. ( 49) states that This corresponds to a matrix equation, where the subscript S refers to the standard Stokes system, originating from the action, Eq. ( 49), v is the vector of unknowns, such that 52) and b is the corresponding right hand side vector, given by The matrix is square and its column or row dimension is 3N + 5.In this case the matrix system, Eq. ( 51), is the form of a saddle point problem (Benzi et al., 2005), common to problems that arise from a constrained optimization.Since the matrix A S originates from a variational principle, we conclude that it is symmetric, i.e., A S = A T S , implying that its eigenvalues are real.However, as a saddle point problem, the matrix is indefinite and therefore it is characterized by both positive and negative eigenvalues.

Discretization of the reformulated Stokes action
In a similar manner, the reformulated action, Eq. ( 41), may be discretized as follows where, from Eq. ( 42) we have and u z is given by the piecewise-linear distribution, Eq. ( 47).In the present one-dimensional problem F z = u z , and therefore in agreement with Eq. ( 43).Substituting this in Eq. ( 54), we obtain The discrete variational principle, as in Eq. ( 50), but this time with respect to u k only, yields the matrix equation where the subscript RS refers to the reformulated Stokes system arising from Eq. ( 57).The vector of unknowns v this time is given by and the right hand side vector b becomes The row or column dimension of this system is just N + 1.Since the matrix A RS arises from the action, Eq. ( 57), we conclude that it is symmetric, A RS = A T RS , and positive definite, implying that its eigenvalues are real and positive.

Numerical iterative solution
For illustration, let us consider the case θ = 18 • , η = 18, and h = 0.02, giving 50 cells in the vertical and a slip parameter γ of 10 %.The numerical solution from Eqs. ( 51) and (58) for the horizontal velocity component u z is shown in Fig. 3 in comparison with the exact solution from Eq. (39).
The numerical properties of a linear system are largely determined by the distribution of eigenvalues associated with the matrix, and in particular by the condition number.In the present example, with N = 50 and h = 0.02, we plot the distribution of the eigenvalues λ for the standard Stokes system, Eq. ( 51), in Fig. 4, and for the reformulated Stokes system, Eq. ( 58), in Fig. 5.Note the presence of both negative and positive eigenvalues in Fig. 4, indicating an indefinite matrix in the standard Stokes system, while all eigenvalues in Fig. 5 are positive, indicating a positive-definite matrix.
From the eigenvalue distribution, the condition number for the standard Stokes system is given by and for the reformulated Stokes system by Thus, the standard Stokes system is quite poorly conditioned, even in this simplified one-dimensional problem.One would  51).Note that negative eigenvalues have been greatly amplified for clarity.expect to see significant numerical errors in comparison with the reformulated Stokes system as one goes to higher resolutions, although this is not yet evident in Fig. 3. Large problems, particularly multi-dimensional ones, are solved by iterative methods.The iterative method of choice for symmetric systems, as in this case, is typically one of several possible Krylov subspace methods.The simpler Krylov methods typically require definite matrices.Indefinite systems, on the other hand, require special methods (Paige and Saunders, 1975;Fletcher, 1976).Furthermore, the convergence of Krylov subspace methods depends on the condition number (Saad, 2003).From this we could infer that the reformulated Stokes system might be much easier to solve than  the standard Stokes saddle-point system.We shall illustrate this by looking at the iterative convergence of the two formulations when using the conjugate gradient method (Saad, 2003), a prominent Krylov subspace method, even though this method is prone to break down for indefinite problems due to a possible division by zero or a near zero.The iteration is initiated with all unknowns set to unity.We plot the convergence history for the two different cases in Fig. 6 for a small problem, N = 50, and in Fig. 7 for a much larger problem, N = 1000.In all cases we plot the horizontal velocity error, i.e., the L2-norm of the difference between the horizontal velocity solution and the exact horizontal velocity from Eq. ( 39), as a function of the iteration number.The initial error is approximately 8 and 35 for the small and large problem, respectively.
The error at convergence (at about iteration 50) of the reformulated system is on the order of 10 −13 , while for the standard Stokes system the corresponding error (at about iteration 150) is on the order of 10 −4 .
In the absence of round-off error the conjugate gradient method gives the exact answer in N steps, where N is the order of the system.This can be seen in Figs. 6 and 7 for the reformulated Stokes case, where the method effectively terminates in 50 and 1000 steps, respectively.For very large problems, and particularly for multi-dimensional problems, it is not feasible to carry on the calculation for that many steps and the method becomes simply an iterative scheme whose convergence depends on the condition number.Indeed, viewed in this way we see that the reformulated Stokes system converges significantly faster than the standard Stokes system.Moreover, the convergence of the conjugate gradient method for the standard Stokes system can break down due to its being an indefinite system.We observe breakdowns in the vicinity of iteration 80 and 125 in Fig. 6 for the small problem, and in the vicinity of iteration 500, 700, and 1600 in Fig. 7 for the larger problem.Although the method recovers and convergence resumes, it does so with a larger error.In fact, for the larger problem, the reformulated system ends up with an error of order 10 −12 after about 1000 iterations, while the standard Stokes system continues to converge very slowly and eventually ends up with a much larger error of about 10 −3 after 4000 iterations.

New Euler-Lagrange equations for the reformulated Stokes system
It is of interest to obtain the partial differential equations that characterize the reformulated Stokes system, if only to compare them with the standard Stokes system, Eqs.(9)(10)(11)(12)(13)(14).For this we need to derive the Euler-Lagrange equations associated with the reformulated Stokes action, which we do next.
Taking the variation of the action, Eq. ( 26), as detailed in DPL10, and making use of Eqs. ( 21) and ( 24), we obtain (63) Note that this is linear in the velocity perturbations δu (i) , δu (b) (i) , and implicitly in δu (s) (i) also.Recall that the variational principle, i.e., Eq. ( 8), implies that the variation of the action, Eq. ( 63), must vanish for arbitrary velocity perturbations.Therefore, Eq. ( 63) must now be manipulated into a form such that the integrands in the volume and surface integrals are linear functions of the velocity perturbations themselves.Since the velocity perturbations are arbitrary, the coefficients multiplying each of the velocity perturbations must vanish, and this gives the required set of Euler-Lagrange equations and also the associated natural boundary conditions.The manipulations required to put Eq. ( 63) into this form are rather complicated.We do this in Appendix A, and obtain where The associated free-stress upper surface boundary condition is and the generalized basal boundary condition becomes As noted earlier, we may set βw (b) n = 0 in Eq. ( 68) except possibly under unusual circumstances.These are the partial differential equations and boundary conditions that constitute the reformulated Stokes problem.The basal boundary conditions include sliding along a rigid bed as well as a generalized floating boundary condition that may, for example, include conditions at the base of an ice shelf.The above equations are very similar to the corresponding Blatter-Pattyn equations (see DPL10) except for extra terms, which we enclose in square brackets for emphasis.These extra terms, in effect, convert the Blatter-Pattyn model into the full-Stokes problem.

Summary and conclusions
We have presented a reformulation of the full Stokes problem for ice sheets that converts it from the standard constrained minimization formulation in six variables (u,v,w,P , s , b ) into an unconstrained minimization in only two variables (u,v).This not only reduces the size of the problem but makes the problem much more tractable numerically.Instead of the original indefinite saddle point problem we obtain a positive-definite minimization or optimization problem that is directly amenable to a number of efficient solution techniques.In this respect, the reformulated problem is similar to the first-order or Blatter-Pattyn approximation, but without the associated approximation errors.An important byproduct of the present formulation is the fact all boundary conditions are already incorporated in the action functional, thereby avoiding many problematic issues with the implementation of boundary conditions in practical models.As an aside, note that this work provides a further example of the usefulness of the fundamental action principle for ice sheets presented in DPL10 and DPL11.
On the negative side, the new system matrix is less sparse and may impose additional continuity requirements on the approximating space, as can be seen from the presence of integrals and (effectively) fourth-order horizontal velocity derivatives in Eq. ( 66).Note, however, that due to the nonlinearity of the problem in general, it might be expected that the JFNK (Jacobian-Free-Newton-Krylov) method of Knoll and Keyes (2004) will be the preferred solution method, in which case only the functional, Eq. ( 22) or (26), is required (i.e., the system matrix is never actually formed) and so only second-order horizontal velocity derivatives are needed.
We have noted many advantageous properties of the reformulated Stokes system compared to the standard Stokes system.We illustrate some of these properties by means of a simple linear two-dimensional ice sheet problem that is reducible to a one-dimensional representation.This simplified problem demonstrates better conditioning and convergence for the reformulated system compared to the standard Stokes system.This is encouraging for the application of the present method to more general problems.At this point it is not possible to conclude how computational costs will compare.This question is beyond the scope of the present paper and can only be answered when the method is implemented and evaluated in realistic, three-dimensional problems.In short, the proposed reformulation isn't likely to solve every problem with full-Stokes modeling, but it is hoped it will ameliorate many of them and will lead to new directions in ice sheet modeling.

Deriving the Euler-Lagrange equations for the reformulated Stokes problem A1 Preliminaries
We shall be making frequent use of the following two results: We have introduced dummy variables z ,z as a reminder that variables other than x,y may be present.A useful special case is given when g z ,x = 1, as follows

A3 Derivations leading to the Euler-Lagrange equations
There now remain two terms in Eq. ( 63) that need to be manipulated into the required form, namely, The first term, I 1 , is by far the most complicated and we shall deal with it first.To do this we shall temporarily assume that the vertical velocity is an independent variable, as in the standard Stokes model, and therefore retain a three-dimensional velocity in the form u i ∈ u (i) ,w .However, from Eqs. ( 19) and ( 21), and noting that δw where The basal surface integral I 13 is the simplest; where, making use of Eq. (A8), it may be rewritten as follows (A12) The upper surface integral I 12 is more complicated.It may be expanded and rewritten as follows

Fig. 1 .
Fig. 1.A schematic diagram of the simplified ice sheet configuration considered in this paper.

Fig. 2 .
Fig. 2. The simple sliding ice sheet test problem configuration, showing the transformed (rotated) coordinate system, x , z .

Fig. 3 .Fig. 4 .
Fig. 3. Horizontal velocity component u z for the case of basal inclination θ = 18 • and basal friction parameter η = 18, corresponding to a basal slip parameter γ = 10%.The exact solution, from Eq. (39), is shown by the solid line.Discrete points from numerical solutions with h = 0.02 for the standard Stokes and the reformulated Stokes cases are shown dotted.The two cases cannot be distinguished visually.

Fig. 6 .
Fig.6.Conjugate gradient convergence history for the reformulated Stokes system (solid line) and the standard Stokes system (dashed line) for a small problem (N = 50, h = 0.02).

Fig. 7 .
Fig. 7. Conjugate gradient convergence history for the reformulated Stokes system (solid line) and the standard Stokes system (dashed line) for a relatively large problem (N = 1000, h = 0.001).
. (23) and Eqs.(20) or (21) define ε2 RS and w u (i) , respectively.The two-dimensional vectors composed of the horizontal components of the unit vectors at the boundaries, i.e., n , are defined in the Appendix.Thus, the Euler-Lagrange equations are given by +h z ,x,b(x) ∂b(x) ∂x − h z ,x,a (x) ∂a (x) ∂x .