Manufactured analytical solutions for isothermal full-Stokes ice sheet models

Introduction Conclusions References


Introduction
Model verification is crucial in developing a numerical model.The ice-sheet modeling community has been using two tools to verify models, comparison of numerically computed solutions to analytical solutions when possible, and intercomparison, that is, measuring differences between various models' results on the sets of simplified geometry benchmark tests.
For shallow-ice approximation (SIA) models, the simplified geometry tests as well as the results of intercomparison of different SIA models can be found in (Huybrechts et al., 1995).As for the exact solutions for SIA equations, two techniques have been used to generate analytical solutions, the similarity reduction technique (an approach that identifies equations for which the solution depends on certain groupings of the independent variables rather than depending on each of the independent variables separately (Nye, 2000;Halfar, 1981Halfar, , 1983;;Bueler et al., 2005) and the manufactured solutions technique (an approach that chooses a reasonable "solution" function, for example, a velocity-field and pressure, substitutes them into the Stokes equations, and determines the body force necessary to make the chosen functions into actual solutions (Bueler et al., 2005(Bueler et al., , 2007;;Bueler and Brown, 2006).
For higher-order models and full-Stokes models, the simplified geometry tests and the results of intercomparison of different models can be found in (Pattyn et al., 2008) and an analysis on the CPU performance of the tests can be found in (Gagliardini and Zwinger, 2008).As for the exact solutions, mathematical work has mainly focused on the flow of linear media, and quasi-analytical solutions have been found for the first-order approximation equations for Published by Copernicus Publications on behalf of the European Geosciences Union.
All the above solutions give physical insight into the flow processes; however, they cannot be easily used to benchmark the numerical solutions.For example, Gudmundsson in (Gudmundsson, 2003) obtained the three-dimensional solution of the linearized zeroth-order problem for a linear viscous medium.To use this solution for benchmarking numerical ice sheet models, the exact error estimate must be known (Raymond and Gudmundsson, 2005).
In this paper, we present the detailed construction of a manufactured exact solution to time-dependent and steadystate isothermal full-Stokes ice sheet problems.The solutions are constructed for three-dimensional (3-D) full-Stokes and two-dimensional (2-D) flowline ice sheet models with variable viscosity.The construction is done by choosing for the specified ice surface and bed the velocity distributions that satisfy the mass conservation equation and the kinematic boundary conditions, and by then calculating the required force distribution that makes the chosen velocities and pressure into exact solutions of the conservation of momentum equation and its boundary conditions.In the appendices we give the formulas that can be used to calculate the compensatory stress terms for the momentum equation in the 2-D and 3-D full-Stokes models and supplement to the manuscript contains a fortran 77 code to calculate stress terms for the 2-D model.
The steady-state solutions constructed in this paper are variations of the benchmark experiments A and B in (Pattyn et al., 2008).However, by substituting different ice surface and bed geometry into the derived formulas, analytical solutions for different geometries can also be constructed.
The boundary conditions can be specified as essential Dirichlet conditions or as periodic boundary conditions.By changing a parameter value, the analytical solutions allow modelers to investigate their solutions for a range of aspect ratios as well as for different, frozen or sliding, basal conditions.Finally, the analytical solutions may help the modelers to estimate the numerical error in the case when the effect of the boundary conditions are eliminated, that is, when the exact solutions values are specified as inflow and outflow boundary conditions.

Model equations
We consider an ice sheet model in the Cartesian coordinates x = ( x, ỹ, z) with the domain 0 ≤ x ≤ L, 0 ≤ ỹ ≤ L, b( x, ỹ) ≤ z ≤ s( x, ỹ, t), where t is time, s( x, ỹ, t) defines the surface and b( x, ỹ) defines the base of the glacier.
Bed elevation b( x, ỹ) and accumulation rate ã are time independent, while surface elevation s( x, ỹ, t) can change with time.The solution is the velocity vector ṽ = ( ũ, ṽ, w) and ice pressure p. Dimensional variables in this work are denoted with a tilde and non-dimensional variables without.
The field equations for the isothermal ice sheet model consist of the conservation of mass and the conservation of momentum: where ρ is the ice density, g is the gravitational acceleration, μ is the effective viscosity , B is a temperature-independent rate factor, and n is the stress exponent.

Boundary conditions
The model is time-dependent in the usual sense that the ice sheet geometry evolves according to a mass continuity equation.We assume that the ice has a hard bed, ∂b ∂t = 0.The kinematic boundary conditions applied at the upper and lower surfaces of the ice mass are as follows: The stress-free boundary conditions at the upper surface s( x, ỹ, t) are defined as: For the frozen-based grounded ice, the boundary conditions at the bed b( x, ỹ) can be specified as Dirichlet conditions: For the ice with sliding bed, the shear stresses may be specified at the bed b( x, ỹ) as Robin conditions: and β2 is the friction coefficient.
Along the glacier's upstream and downstream boundaries, either periodic where f = ũ, ṽ, w, p, or Dirichlet boundary conditions may be specified.
Here we assume that functions fexact are known.

Dimensionless equations
To non-dimensionalize variables, we choose the following typical values: Z -the mean thickness of the ice-sheet, Lthe length of ice-sheet, U -a typical velocity in the horizontal direction, W -a typical velocity in the vertical direction, P -the mean pressure, A -the mean accumulation/ablation rate, and introduce the following non-dimensional variables (variables without tilde): To further simplify the equations, we introduce the aspect ratio parameter δ: and require that scale factors L, U , W , and P satisfy the following relationships: The nondimensional steady-state conservation of mass and momentum equations are then as follows: The stress-free boundary conditions at the upper surface s(x,y,t) become as follows: To satisfy the 2-D version of the kinematic boundary conditions (13-14), we assume that in the interior of the domain, where s(x,t) > b(x), the vertical velocity w is From ( 21), it follows that If we substitute (22) into the incompressibility Eq. ( 9), we obtain the following equation containing only variable u and its derivatives: Equation ( 23) is a first-order quasi-linear partial differential equation with two independent variables (x and z) and one dependent variable (u).The system of ordinary differential equations is called the characteristic system of Eq. ( 23).If we can find two particular independent solutions of this system, which are called the integrals of system (24), in the form where c 1 and c 2 are arbitrary constants, then the general solution of Eq. ( 23) can be written as where θ is an arbitrary function of two variables.With Eq. ( 26) solved for φ, the general solution can be written in the form where ϑ is an arbitrary function of one variable.Thus, to solve Eq. ( 23), we have to find integrals φ and ψ of the system (24).The first integral of the system (24) can be found by solving equation Equation ( 28) can be re-written as follows: We multiply both sides of Eq. ( 29) by s − b and recognize that the left side of the equation is now the following product rule, (s . After replacing the left side of the equation with this product rule, we obtain: Equation ( 30) has a solution where c 1 is a constant.
The second integral of the system (24) can be found by solving equation After multiplying both sides of Eq. (32) by 1 s−b , the equation can be transformed into: Equation ( 33) has a solution Thus, the general solution of Eq. ( 23) can be written as where θ is an arbitrary function of two variables.With Eq. ( 35) solved for u, the general solution can be written in the form  where ϑ is an arbitrary function of one variable.
The formula (36) shows that the functions satisfying the kinematic boundary conditions (13-14) and the conservation of mass Eq. ( 9), derived under assumption (21), depend on the form of the function ϑ and ice surface and bed curves.
To generate a solution similar to the benchmark experiment B in (Pattyn et al., 2008) and to keep the mathematics simple, choose function ϑ as follows: where λ, c x , and c b are constants.The first term on the righthand side of (37) may be considered as component of velocity associated with internal deformation, and c b as the basal sliding velocity coefficient.
Then the velocity field satisfying the 2-D versions of the kinematic boundary conditions (13-14) and the conservation of mass Eq. ( 9) is: For a zero-accumulation ( ȧ = 0) steady-state ( ∂s ∂t = 0) flow with frozen bed (c b = 0), the horizontal velocity scaled to the surface velocity can be written as a function of ice scaled This expression shows that the horizontal velocity from internal deformation increases with power λ of ice depth.For λ = 4 this is consistent with lamellar flow (der Veen, 1999) as shown in Fig. 1.
As can be seen from ( 38) for a zero-accumulation ( ȧ = 0) steady-state ( ∂s ∂t = 0) flow, if λ > 0, then These expressions show that c b can be interpreted as the iceflux due to sliding flow and c x can be interpreted as the iceflux due to deformation flow.
In addition to velocities, the ice pressure function should also be constructed.
The manufactured solution for the ice pressure can be chosen, for example, as in Pattyn's higher-order model (Pattyn, 2003): or in nondimensional form: The constructed velocity and pressure functions do not necessarily satisfy the conservation of momentum Eqs.(10)(11)(12) or the surface and basal boundary conditions (15-17) and (18-20).To make the constructed velocity and pressure The Cryosphere, 4, 285-311, 2010 www.the-cryosphere.net/4/285/2010/functions into exact solutions of these equations, we substitute them into the equations and calculate the right-hand side functions that match these solutions.This can be done when a specific surface s(x,t) and bed b(x) are chosen.Equations (38-39) also satisfy prognostic equation describing the change of local ice thickness h(x,t) = s(x,t) − b(x) in space: Equations (38-39) and ( 41) are solutions of flow with a general surface s(x,t) and bed b(x).Below are specific solutions for a particular case of an ice surface and a sinusoidal bed, similar to the benchmark experiment B in (Pattyn et al., 2008).

A manufactured solution for a time-dependent flow with a sinusoidal bed
To generate a particular solution, assume a flow with zero accumulation/ablation rate, ȧ = 0, a sinusoidal bed defined as in (Pattyn et al., 2008), and an ice surface that changes from a linear sloping surface to the one that is draped over the topography of the bed: where Constant c t shows how fast the ice surface changes (relative to the value of the ice bed) at the beginning of the test: For a flow down an infinite plane with a mean inclination tan(α), periodic boundary conditions for a function f are defined as follows: f (0,z+tan(α)) = f (1,z) and the analytical solutions ( 38), ( 39), ( 41) satisfy these conditions for geometry (43-44).
Appendix A contains the formulas and supplement to this manuscript contains a simple fortran 77 code that can be used to calculate the exact solutions and compensatory stress terms for the momentum equation in the 2-D flowline model.The code dumps the generated solutions to specified files.All input data are specified in file parameter2d.h in the supplement.
Values of flow parameters and constants are chosen from (Pattyn et al., 2008) and are given in Table (1), the starting linear slope of the ice surface α = 0.5 • .The length scale of the domain is chosen 80 km, which results in aspect ratio δ = 1/80.
Constants of the test are chosen as follows: coefficient in (45) c t = 10 −6 and coefficients in (38) c x = 10 −6 , c b = 10 −6 , and λ = 4.This experiment can be considered as an ice-stream flow over a bumpy bed.The values of constants c x , c b , and c t chosen to generate a reasonable dimensional values of the flow functions which calculated from nondimensional values using formulas ( 6) and ( 8).For the flow parameters from (Pattyn et al., 2008), the value of the nondimensional ice flux will be around 10 −6 : A. Sargent, J. L. Fastook: Analytical solutions for isothermal full-Stokes ice sheet models The choice of parameter c t was dictated by the typical scale value of the time parameter: Velocity is shown in km/yr and pressure in MPa.

Figure (
2) shows the bed (44) and the time-evolution of the ice surface (43) (left graph) and the time-evolution of the norm of the surface velocity (right graph) over 14-year period.The ice surface changes from a linear sloping surface to the surface draped over the topography of the bed.Ice thickness is spatially uniform when the steady-state solution is reached.The surface velocity at the beginning is anti-correlated with the ice thickness -it is larger over the bump than over the trough.At the steady-state, the surface velocity is spatially uniform and does not depend on the bed topography.
Figure (3) shows the horizontal velocity, vertical velocity, and pressure at the beginning (left graphs) and at the time when the steady-state solution is reached (right graphs).
Figure ( 4) shows the compensatory horizontal and vertical stress terms in the conservation of momentum equation at the beginning (left) and at the time when the steady-state solution is reached (right).At the beginning both stress terms have largest values over the bump.At the steady-state solution, the stress terms are zeroes almost everywhere except a small surface layer over the bump.

A steady-state manufactured solution for a flow with a linear sloping surface and a sinusoidal bed
To generate a steady-state solution, assume that in (43) the function γ (t) = 0, that is, a linear sloping surface and a sinusoidal bed are defined similar to the ones of the benchmark experiment B in (Pattyn et al., 2008): If we substitute the above functions for bed and surface into (38-39), then the corresponding steady-state flow's velocities are as follows: Choice of coefficient c b = 0 generates frozen bed flow with zero basal velocities, while c b = 0 generates flow with a sliding bed.
As can be seen from (48-49), if λ > 0 then The last expression shows the conservation of mass flux, q = hu = c x = constant.This anti-correlated relationship between surface horizontal velocity and ice thickness is consistent with the simulation of the smallest length scale L = 5km Experiment B in (Pattyn et al., 2008), by all full-Stokes models.Figure ( 5) shows the horizontal and vertical velocity, ice pressure, and the norm of the surface velocity corresponding to the flow with a linear sloping surface with a slope α = 0.5 • and a frozen sinusoidal bed (c b = 0).The constants in (48) are chosen as c x = 10 −6 and λ = 4.0 and the aspect ratio δ = 1/80.

Analytical manufactured solutions of the 3-D isothermal full-Stokes ice-flow model
Assume as in the 2-D case that in the interior of the domain, s(x,y,t) > b(x,y), the vertical velocity w is: If we substitute (51) into the incompressibility Eq. ( 9), we obtain the following equation containing only variables u, v and their derivatives: Equation ( 52) is a first-order quasi-linear partial differential equation with three independent variables (x, y, and z) and two dependent variables (u and v) of type: Similar to the 2-D flowline manufactured solutions, we choose velocity u(x,y,z,t) as the following function: or where γ 1 , λ 1 , c x , c bx are constants, d(x,y,z,t) = s−z s−b is scaled ice depth, and h(x,y,t) = s − b is ice thickness.
Then the derivatives of function u(x,y,z,t) are www.the-cryosphere.net/4/285/2010/The Cryosphere, 4, 285-311, 2010 A. Sargent, J. L. Fastook: Analytical solutions for isothermal full-Stokes ice sheet models The graphs show the horizontal velocity u, the vertical velocity w (both in km/yr), and the ice pressure p (in kPa) at the beginning (left) and at the time when the steady-state solution is reached (right).At the beginning, the horizontal velocity is anti-correlated with ice thickness: it is larger over the bump than over the trough.At the steady-state, the horizontal velocity is spatially uniform and increases from the bed to the surface with power of λ = 4 of the ice thickness.At the beginning, the vertical velocity is largest at the bed where ice shearing is the largest.At the steady-state, the vertical velocity is almost uniform in every vertical slide.This is consistent with ice-stream flow.The ice pressure is proportional to the ice thickness.At the steady-state, it is equal zero at the ice surface.
The   The characteristic system of Eq. ( 56) is as follows: Two independent particular solutions of this system can be found by solving the equations: along the y = 1/4 slide at the beginning (left) and at the time when the steady-state is reached (right).At the beginning, velocity has two local maximums, over the bump and over the bed where the bed changes the most.At the steady-state position, velocity is spatially uniform and proportional to the ice thickness.+ v 2 + w 2 1/2 (right) change along y = 1/4 slide.Ice surface and the norm of the surface velocity are shown every 1.5 years over the 14-year period, green curves are the initial values and red curves are the final values.At the beginning, velocity has two local maximums, over the bump and over the bed where the bed changes the most.At the steady-state position, the norm of the surface velocity is spatially uniform.

Equation (58) has a solution
where c 1 is a constant.
Equation ( 59) can be re-written as follows: where a(y) is an unknown function.
Substituting Eq. ( 62) into Eq.( 61), we obtain an equation for a: which has a solution: Substituting ( 63) into Eq.( 62), we obtain Then, the general solution of Eq. ( 56) can be written as where θ is an arbitrary function of two variables.With Eq. ( 65) solved for v, the general solution can be written in the form where ϑ is an arbitrary function of one variable.
If we assume again that function ϑ in ( 66) is of the form where λ 2 , c y , and c by are constants, then functions ( 54), ( 50), and (66) satisfying the mass balance equation and the kinematic boundary conditions are as follows: The manufactured solution for the ice pressure can be chosen again as in Pattyn's higher-order model (Pattyn, 2003): or in nondimensional form: where non-dimensional ice viscosity . (72) The constructed velocities satisfy the surface and bed kinematic boundary conditions (13-14) and the mass conservation Eq. ( 9).However, the constructed velocities and pressure do not necessarily satisfy the conservation of momentum equations and the basal and surface boundary conditions.To make the constructed functions into exact solutions of these equations, we substitute them into those equations and calculate the right-hand side functions which accommodate the solutions.This can be done when specific surface s(x,y,t) and bed b(x,y) are chosen.
The constructed solutions do not satisfy ice-sheet evolution equation describing the change of local ice thickness h(x,y,t) = s(x,y,t) − b(x,y) in space: To make the constructed functions into exact solutions of Eq. ( 73), the equation can be modified by adding to the righthand side of the equation a compensatory term.

A time-dependent analytical solution for a flow with a sinusoidal bed
To generate a particular solution, assume a flow with a sinusoidal bed defined similar to the bed in the benchmark experiment A in (Pattyn et al., 2008) and an ice surface that changes from a linear sloping surface to the one that is draped over the bed: where η(x,y) = 1 2 sin(2π x)sin(2πy),γ (t) = 1 − e −c t t , c t is a constant.
To calculate integral in (69), substitute functions (74-75) for bed and surface into the integral in (69).Since it is difficult to calculate the integral analytically for general constants γ 1 and λ 1 , particular values, for example, γ 1 = 1 and λ 1 = 1, can be chosen.
If we substitute the calculated integral and functions (74-75) for bed and surface into (68-70), we obtain the following formulas for velocities: For a flow down an infinite plane with a mean inclination tan(α), periodic boundary conditions for a function f are defined as follows: The constructed solutions (77-79), (71) satisfy periodic boundary conditions only in the horizontal x-direction and do not satisfy periodic boundary conditions in the horizontal y-direction for all values of the input parameters.To satisfy periodic boundary conditions in all lateral directions, the accumulation-ablation rate may be chosen as follows: Appendix 4.1.1contains the formulas that can be used to calculate the compensatory stress terms for the momentum equation.For the 3-D ice-stream flow over a bumpy bed experiment, the parameters of the flow are chosen as follows: the horizontal domain is chosen 80 km×80km which results in aspect ratio δ = 1/80, the starting linear slope of the ice surface α = 0.5 • , sliding bed parameters c bx = c by = 10 −8 , and the remaining constants in ( 77) and ( 78) c x = c y = 10 −6 , λ 2 = 4, and c t = 10 −6 .As in 2-D case, all graphs are given for the dimensional values of variables which are calculated from non-dimensional values using formulas (6).
Figure (6) shows the bed (75) and the ice surface (74) at the time zero and at the time when the steady-state solution is reached.Ice flow is from left to right.The ice surface changes from a linear sloping surface to the surface draped over the topography of the bed.Ice thickness is spatially uniform when the steady-state solution is reached.
Figure (7) shows the horizontal and vertical velocity at the beginning.At the steady-state, the horizontal velocity field is smoothed out, both x-and y-horizontal velocities are almost spatially uniform (≈ 46 km/yr).
Figure ( 8) shows the vertical velocity and pressure at the beginning and at the time when the steady-state solution is reached.
Figure ( 9) shows the norm of the velocity along the y = 1/4 slide at the beginning and at the time when the steadystate is reached.Figure (10) shows change over time of the ice surface elevation and the norm of the surface velocity along the y = 1/4 slide.At the beginning, velocity has two local maximums, over the bump and over the bed where the bed changes the most.At the steady-state position, the norm of the velocity is spatially uniform and at each vertical slide is increasing with ice thickness.

A steady-state analytical solution for a flow with a linear sloping surface and a sinusoidal bed
To generate a steady-state solution, assume that in (74) the function γ (t) = 0, that is, a linear sloping surface and a sloping sinusoidal bed are defined as in the benchmark experiment A in (Pattyn et al., 2008).
The coefficients are α = 0.5 • , λ 2 = 2.25, c x = c y = 1, c bx = c by = 0, δ = 1/80, and accumulation rate ȧ = π sin(4πx) 4 . All functions, the surface horizontal x-and y-velocities, the vertical z-velocities as well as the surface ice pressure, for this steady-state experiment are very similar to the corresponding graphs in Figs. ( 7) and ( 8) of the time-dependent experiment at the beginning time.

Conclusions
The detailed constructions of manufactured exact solutions to 3-D and 2-D flowline time-dependent and steady-state isothermal full-Stokes ice sheet problems are presented.The solutions are valid for non-linear Glen-type flow.The construction of exact solutions done by using manufactured solution technique (Bueler et al., 2007) while the suggested experiments follow directly from ice sheet intercomparison (Pattyn et al., 2008).
The steady-state solutions, constructed in this paper, are variations of the benchmark experiments A and B in (Pattyn et al., 2008).However, by substituting different ice surface and bed geometry formulas into the derived formulas, analytical solutions for different geometries can also be constructed.
Although artificially constructed, the solutions may be useful for testing numerical methods.They offer several ben-efits to potential ice sheet modelers.By changing a parameter value, the analytical solutions will allow the modelers to investigate their algorithms for a different range of aspect ratios as well as for different, frozen or sliding, basal boundaries.The lateral boundary conditions can be specified as periodic boundary conditions or as essential Dirichlet conditions.Specifying Dirichlet conditions, when the exact solutions are specified as inflow and outflow boundary conditions, allows the modelers to check the model accuracy in the inside of the problem domain.

Appendix A Calculation of compensatory stress functions in the 2-D flowline full-Stokes diagnostic equations A1 Compensatory terms in diagnostic equations and in the boundary conditions
The constructed velocities (38-39) satisfy the 2-D versions of the surface and bed kinematic boundary conditions (13-14) and the mass conservation Eq. ( 9) but do not necessarily satisfy the conservation of momentum Eqs.(10-12) and its basal and surface boundary conditions (15-17) and (18-20).Following (Bueler et al., 2007), we introduce compensatory stresses x and z in the conservation of momentum equations to make the chosen velocity and pressure functions into exact solutions of the equations.
To make the chosen velocities satisfy the boundary conditions, we introduce compensatory terms υ x ,υ z ,τ b , and τ z in the boundary conditions.
A. Sargent, J. L. Fastook: Analytical solutions for isothermal full-Stokes ice sheet models At the upper surface s(x,t), the boundary conditions are as follows: 1

A2 Calculation of derivatives
Calculation of the compensatory stress terms requires calculation of derivatives of the exact solutions ( 38), (39), and (41).To simplify calculation of the derivatives, we re-write these functions as follows: and the second derivatives are where, for a surface (43) and a sinusoidal bed (44), If we name the expression For further calculations we need the following derivatives: ) Substituting (A9-A10) and (A13-A14) into (A1-A2), (A3-A4), and (A5-A6) generate formulas for compensatory terms x , z ,υ x ,υ z ,τ x , and τ z .If constant λ in (37) is chosen so that λ > 2, then the calculation of the second derivatives is well defined.

Appendix B Calculation of compensatory stress functions in 3-D full-Stokes diagnostic equations B1 Compensatory terms in diagnostic equations and in the boundary conditions
The constructed velocities (68-70) satisfy the surface and bed kinematic boundary conditions (13-14) and the mass conservation Eq. ( 9).They do not necessarily satisfy the conservation of momentum equations and its basal and surface boundary conditions.Following (Bueler et al., 2007), we introduce compensatory stresses x , y , and z in the conservation of momentum equations to make the chosen velocity functions into exact solutions of the equations.

B2 Calculation of derivatives
Calculation of the compensatory stress terms requires calculation of derivatives of the exact solutions (68-70, 71).To simplify calculation of the derivatives, we re-write these functions as follows: A. Sargent, J. L. Fastook: Analytical solutions for isothermal full-Stokes ice sheet models   For further calculations we need the following derivatives: Substituting (B10-B14) and (B17) into (B1-B3), (B4-B6), and (B7-B9) generate formulas for compensatory terms x , y , z ,υ x ,υ y ,υ z ,τ x ,τ y , and τ z .If constant λ 2 in (B11) is chosen so that λ 2 > 2, then calculation of the velocities' first and second derivatives is well defined.

Fig. 1 .
Fig. 1. 2-D flowline steady-state manufactured solution (coefficient λ = 4): horizontal component of velocity scaled to the surface velocity as a function of dimensionless thickness.Horizontal velocity increases with the fourth power of ice thickness.Most shearing ice concentrated near the glacier base which is similar to lamellar flow.
Fig. 2. 2-D flowline time dependent experiment -ice stream flow over bumpy bed; the steady bed and transformation over time of the ice surface (left) and transformation over time of the norm of the surface velocity (right).Ice surface and the norm of the surface velocity are shown every 1.5 years over the 14-year period, green curves are the initial values and red curves are the final values.

Fig. 3 .
Fig. 3. 2-D flowline time dependent experiment -ice stream flow over bumpy bed.The graphs show the horizontal velocity u, the vertical velocity w (both in km/yr), and the ice pressure p (in kPa) at the beginning (left) and at the time when the steady-state solution is reached (right).At the beginning, the horizontal velocity is anti-correlated with ice thickness: it is larger over the bump than over the trough.At the steady-state, the horizontal velocity is spatially uniform and increases from the bed to the surface with power of λ = 4 of the ice thickness.At the beginning, the vertical velocity is largest at the bed where ice shearing is the largest.At the steady-state, the vertical velocity is almost uniform in every vertical slide.This is consistent with ice-stream flow.The ice pressure is proportional to the ice thickness.At the steady-state, it is equal zero at the ice surface.

-Fig. 4 .Fig. 5 .
Fig. 4. 2-D flowline time dependent experiment -ice stream flow over bumpy bed.The graphs show the compensatory horizontal x and vertical z stress terms (in kJ) at the beginning (left) and at the time when the steady-state solution is reached (right).At the beginning the stress terms have the largest values over the bump.At the steady-state solution, they are zeroes almost everywhere except a small surface layer over the bump.

Fig. 6 .Fig. 7 .
Fig. 6. 3-D time dependent experiment -ice flow over a bumpy bed.The graphs show the bed and ice surface at the beginning (left) and at the steady state (right).All distances are scaled.Ice flow is from left to right.The ice surface changes from a linear sloping surface to the surface draped over the topography of the bed.Ice thickness is spatially uniform when the steady-state solution is reached.

Fig. 8 .
Fig. 8. 3-D time dependent experiment; the graphs show the ice surface vertical velocity w (in km/yr) and the ice surface pressure p (in kPa) at the beginning (left) and at the time when the steady-state solution is reached (right).
Fig. 9. 3-D time-dependent experiment.The graphs show the the norm of the velocity u 2+ v 2 + w 2 1/2 along the y = 1/4 slide at the beginning (left) and at the time when the steady-state is reached (right).At the beginning, velocity has two local maximums, over the bump and over the bed where the bed changes the most.At the steady-state position, velocity is spatially uniform and proportional to the ice thickness.

Fig. 10 .
Fig. 10.3-D time dependent experiment; the ice surface elevation (left) and the norm of the surface velocity u 2+ v 2 + w 2 1/2 (right) change along y = 1/4 slide.Ice surface and the norm of the surface velocity are shown every 1.5 years over the 14-year period, green curves are the initial values and red curves are the final values.At the beginning, velocity has two local maximums, over the bump and over the bed where the bed changes the most.At the steady-state position, the norm of the surface velocity is spatially uniform.

Table 1 .
Constants for the benchmark experiments.
To make the chosen velocities satisfy the boundary conditions, we introduce compensatory terms υ x ,υ y ,υ z ,τ x ,τ y , and τ z in the boundary conditions.At the upper surface s(x,y,t), the boundary conditions are as follows: