Applicability of the Shallow Ice Approximation inferred from model inter-comparison using various glacier geometries

Applicability of the Shallow Ice Approximation inferred from model inter-comparison using various glacier geometries M. Schäfer, O. Gagliardini, F. Pattyn, and E. Le Meur LGGE, CNRS, UJF-Grenoble, BP 96, 38402 Saint-Martin d’Hères Cedex, France Laboratoire de Glaciologie, ULB, Bruxelles, Belgium Received: 30 May 2008 – Accepted: 10 June 2008 – Published: 14 July 2008 Correspondence to: M. Schäfer (smartina@gmx.de) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
The Shallow Ice Approximation (SIA) has been widely used to simplify the Stokes equations, allowing the construction of models that can simulate the flow of large ice masses like the entirerty of the Antarctic or Greenland ice-sheets.The use of SIA is well founded as long as the horizontal gradient of the bedrock and surface topographies are much smaller than the vertical ones (Hutter, 1983), which is the case in a large part of Antarctica and Greenland.The SIA presents, besides the convenience of easiness in coding, the huge advantage of being very efficient in terms of computing resource, both from the memory and time CPU consumptions.
These advantages certainly explain why the SIA has also been applied to smaller ice masses like mountain glaciers (Le Meur and Vincent, 2003;Hubbard, 2000).With the help of a 2-D model, Leysinger Vieli and Gudmundsson (2004)  the snout position of a valley glacier is relatively well reproduced by a SIA model.This is confirmed in the present work with a 3-D model.Greuell (1992) states also that longitudinal stresses have no big influences on the response of an alpine glacier (the Hintereisferner in Austria).
Recent papers have started to contemplate the applicability of SIA to smaller ice masses (Leysinger Vieli and Gudmundsson, 2004;Le Meur et al., 2004;Hindmarsh, 2004).The Ice Sheet Model Inter-comparison Project for Higher-Order ice-sheet Models (Pattyn and Payne, 2006) has clearly shown that SIA is not longer valid to compute a diagnostic velocity field for a given geometry, especially for high aspect ratios (Pattyn et al., 2008).With the SIA, the velocity is only function of the local surface slope and elevation, and no stress transfer from one grid point to the other is allowed.In cases where the bedrock topography or the basal sliding conditions change within grid points over typical wavelengths that are much smaller than the ice-thickness, as is the case for mountain glaciers, one can then easily understand why the SIA is not longer valid.
In the continuation of these studies, this paper presents several original timedependent test cases for which the SIA is compared to two other models which incorporate longitudinal stresses.The synthetic glacier geometries have been chosen to mimic real glaciers, but in a simplified way.Therefore, small scale variation of the basal topography or sliding condition are neglected, and the difference between the models can only be attributed to large scale phenomena like the transmission of longitudinal stresses within grid points or the transmission of lateral shearing in a valley.The paper is organised as follows: we first recall the equations governing the flow of an isothermal non-linear viscous glacier (Sect.2).Then the different models are presented on the basis of the simplification of the mentioned equations, and some technical points are discussed (Sect.3).The next section presents the three glacier geometries that will be used to compare the three models (Sect.4).The main section discusses the results obtained for both diagnostic and time-dependent simulations for the three models applied to the three glacier geometries (Sect.5).The CPU time of each model is also discussed in this section.

Ice flow models
The modelling of glaciers and ice sheets has a large range of applications, such us ice cores dating or simulating the role of large ice masses on the climatic system.Mountain glaciers are important not only because they act as good climatic indicators but also for their possible economic impacts and the role they can play in natural catastrophes.Whatever the objectives, numerical flow models seem to be the best way to capture the complexity of glacier evolution.As depicted in Fig. 1, an ice flow model solves the momentum and mass conservation equations for given initial and boundary conditions.The evolution of a glacier, i.e. the evolution of its surface geometry, is then fully determined by the ice flow and the climatic forcing.In what follows, we will only consider isothermal glaciers, so that the climatic forcing reduces to a mass balance distribution at the surface of the glacier.
A right-handed (O, x, y, z) Cartesian coordinate system as depicted in Fig. 2 is used.The surface and bedrock are defined as z=S(x, y, t) and z=B(x, y, t) respectively, but the bedrock will be assumed to be a fixed boundary so that z=B(x, y).The icethickness is given by H=S−B.In what follows, the ice is considered as an isothermal non-linear viscous incompressible material.
Due to ice incompressibility, the mass conservation equation can be expressed as: Here, v x , v y , v z denote the components along the three directions of space of the velocity vector v .The stress tensor τ splits into a deviatoric part τ and the isotropic pressure p, where δ i j is the Kronecker symbol.The constitutive relation, that links deviatoric stresses to strain rates ˙ , is a power-law, referred to as the Glen's flow law in glaciology: where A(T ) is the deformation rate factor (hereafter reduced to A for the isothermal ice body, see Table 2 for its adopted value) and τ is the second invariant of the deviatoric stress tensor which is defined as: The strain-rate components ˙ i j are defined from the velocity components as: The force balance (quasi static equilibrium) in the three directions of space lead to the Stokes equations: where σ i stands for τ i i , g is the norm of the gravitational acceleration vector and ρ is the glacier ice density (see Table 2).
In temperate glaciers, in addition to the ice deformation, basal sliding contributes to the ice motion, but will be omitted in this study.At the stress-free upper surface (air-ice interface) the following kinematic equation applies for all points such that z=S(x, y, t): where a is the mass balance function, considered as a vertical flux.These equations finally lead to a complex system of five partial differential equations with five independent unknowns (three terms of the velocity vector, one for isotropic pressure and another for the free surface elevation).Generally no analytical solution is possible.Different models thus use different simplifications and numerical methods to solve this set of equations.

Description of the three different models
In this inter-comparison study, three different models with different degrees of complexity will be used, i.e. a simple zero-order SIA model (Sch äfer and Le Meur, 2007), a higher-order (HO) model (Pattyn, 2003) and a full-Stokes (FS) model using the finite element code Elmer (Elmer Manuals, CSC, 2006).In the interest of brevity and because it did not systematically succeed in converging, results from the HO model are omited for some simulations.The aim of this section is to present in detail these three different models.

SIA model
In the SIA formulation, important simplifications are introduced in the equations presented above from a scale analysis by which the orders of magnitude of the various variables are assessed.All variables are expressed as the product of a characteristic value and a dimensionless quantity.Following the scale analysis from Hutter (1983) and Morland (1984), the dimensionless variables are expanded in (perturbation) power Introduction

Conclusions References
Tables Figures

Back Close
Full series of the aspect ratio as small parameter, assuming the aspect ratio is defined as and expressing the shallowness of the ice-sheet or glacier.
[H] and [L] are respectively a characteristic horizontal and vertical dimension of the studied ice-body.The characteristical horizontal and vertical velocities ([V L ] and respectively [V H ]) scale in the same way as the corresponding geometrical quantities, However, as mentioned by Baral et al. (2001), the scaling of the velocities in Eq. ( 9) is valid only for land based ice masses, but not for ice shelves.
All variables in the equations presented above are replaced by their corresponding power series.In the zeroth-order SIA case only the remaining O-th-order terms are kept.This implies that in the stress tensor all components vanish besides τ xz and τ yz .σ x , σ y , σ z , τ xy τ xz , τ yz (10) A rigorous derivation has been developped by Hutter (1983) and can also be found in Baral et al. (2001).It is clear that <1 is a prerequisite for the SIA to apply and that the smaller , the more accurate is the approximation.Moreover, Morland (1984) shows that the steeper the bedrock slope, the less valid the SIA.Specifically, the SIA (0-th order) leads to where v b is the contribution of basal sliding if any, v ⊥ =(v x , v y ) stands for the two horizontal velocity components and ∇ ⊥ defines the horizontal components of the gradient operator ∇ (i.e.∇ ⊥ •=(∂•/∂x, ∂•/∂y)).Integration of the horizontal velocities from the ice bottom B to the upper free surface S gives the horizontal fluxes as Vertical integration of the mass conservation Eq. (1) from B to S assuming the kinematic boundary condition (7) leads to a transport equation (Greve, 1997): Assuming ∂B/∂t=0, this transport equation is transformed into a diffusion equation: where D is analogous to a diffusivity which reads when sliding is neglected.
Most of the time this equation is treated numerically with a semi-implicit (Sch äfer and Le Meur, 2007) and in some cases with an over-implicit scheme (Hindmarsh, 2001) after being discretized according to a finite-difference method on a staggered regular 50 m×50 m grid.

Higher-order model
The model used here is the model developed by Pattyn (2003).Two approximations are applied on the initial set of equations presented above (Sect.2, Eqs. 1 to 7).The major assumption consists of applying the hydrostatic approximation on the momentum equation.With this assumption, the third equation of the momentum set of Eqs. ( 6) reduces to and σ z can be directly obtained by integration.Assuming a zero normal stress at the surface, the first two equations can be rewritten as follows: (20) In this model, the Glen's law is expressed with the help of the second invariant of the strain rate tensor ˙ defined as The effective viscosity, defined as where ˙ 0 is a small artificial number in order to prevent viscosity from becoming infinite in case of zero strain-rate.
The second assumption of the model is that the horizontal gradients of the vertical velocity are assumed to be small compared to the vertical gradient of the horizontal velocities: With this second assumption, the strain-rates ˙ xz and ˙ yz reduce to Combining the Glen's flow law ( 22), the horizontal stress field Eqs. ( 20), the strain rate components and the equation of mass conservation (1) gives By replacing the strain-rates by their expression function of the velocity assuming Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
Eq. ( 24), these equations read: where the effective viscosity is given by Solving these two linear equations leads to the horizontal velocity components.These equations are treated as a pair of coupled linear equations with v x and v y as unknowns and η is approximated from a previous iteration.The problem is solved in two steps: (i) solving a linear set of equations using an estimate of the horizontal velocity field and (ii) iterating the nonlinear part by updating η with new estimates of v x and v y .
The vertical velocity is obtained by integration of the mass conservation equation.
To obtain the new surface elevation the kinematic Eq. ( 7) is rewritten as a diffusion equation, the resulting sparse system is solved using a conjugate gradient method.For further details see Pattyn (2003).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

Full-Stokes model
The numerical solution of the FS equations is obtained using the Finite Element Method based on the code Elmer (Elmer Manuals, CSC, 2006).In the present case, the free surface Eq. ( 7) and the Stokes Eq. ( 6) are coupled and solved iteratively using an implicit scheme during the increment time step.

Free surface equation
The non-integrated Eq. ( 7) is written in its discrete variational form: where S is developed as S(x, y, t)=φ i (x, y)S i (t) and where Ψ is a test function.S i stands for the discrete value of S at the i -th node of the domain.Stabilisation is obtained by applying the Stabilised Method (Franca et al., 1992).The vertical redistribution of the mesh nodes according to the moving boundary free surface is done by solving a linear elasticity equation.

Stokes equation
The mass conservation Eq. ( 1) and the Stokes Eqs. ( 6) are written in their variational form A vector-like test function Φ and the weight functions ψ i are introduced, the left-hand side term in the momentum equation has been integrated by parts and reformulated by applying the Green's theorem.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The deviatoric stress tensor τ is expressed in terms of the strain-rate tensor ˙ with the help of inversion of the power law (3).The non-Newtonian stress-strain relation introduces non-linearities into the system implying the application of an iteration scheme.The numerical solution of this system of equations is obtained either by the Stabilised (Franca and Frey, 1992) or the Residual Free Bubble method (Baiocchi et al., 1993).
For more details about the numeric in Elmer and the way the Stokes and free surface equations are solved, the reader can refer to Le Meur et al. ( 2004); Zwinger et al. (2007); Gagliardini and Zwinger (2008).

The different glacier geometries
Three simple-shaped synthetic glaciers are tested.

Flattened half-sphere
The first one, a flattened half-sphere on an inclined ramp (along the y-direction), is the one that was already used by Le Meur et al. (2004) and Sch äfer and Le Meur (2007).This half-sphere has a radius R of 500 m and is flattened with a factor 0.3 (see the corresponding profile on Fig. 3).Its center is at the origin.To asses the validity of the SIA for increasing bedrock slopes, different slopes p from 0 to 0.3 are used.
where all dimensions are expressed in meter.

Conic glacier
The second glacier shape is chosen in analogy to a real case glacier, the photo and profile on Fig. 4).The flat zone (radius R=400 m) in the centre represents in a very simple way the crater which is initially set free of ice and where the prescribed mass balance equals zero.
The ice-thickness is chosen to be uniform (40 m) over the whole glacier, except over the upper and bottom parts where it is assumed to reduce linearly to zero from respectively 5500 to 5800 m and from 4880 m to 4800 m.
The bedrock slope p is varied from 0.3 to 0.8 as well.The real Cotopaxi glacier has a mean slope of 0.55.The Cotopaxi geometry used is where all dimensions are expressed in meter.

Valley glacier
The last glacier is also chosen close to a real case glacier, a typical mountain valley glacier (ex.Glacier d'Argenti ère in the Mont Blanc range).This glacier is composed of an ablation valley tongue (denoted V in the formulas below) and a cirque (C in the formulas).y is the direction of the ice-flow (the origin is somewhere below the snout), x is orthogonal to this direction with x=0 at the centre of the glacier.At the interface (y=5000 m) of these two parts of the glacier, the icethickness is given by the parameter E 0 and is varied from 100 m to 700 m in order to modify the aspect ratio of the glacier (aspect ratios ranging from 0.1 to 0.7).
In the valley part, the bedrock is the superposition of a linear down-sloping function in the direction of ice flow (slope 0.1) and a parabolic function in x (orthogonal to the direction of flow) with a flat zone in the middle (again, all the dimensions are meter), Figures

Back Close Full Screen / Esc
Printer-friendly Version Interactive Discussion see formulas below and Fig. 5. with In the upper part of the glacier the bedrock is given by a 3rd degree function where the different parameters are defined by the following boundary conditions: For the upper part (cirque) this initial surface reads: whereas for the lower part (valley), the following quadratic reduction of the ice-thickness is chosen; in which the parameters are chosen for the surface function to be a one time differentiable function at the interface between the two parts of the glacier and for the snout to be positioned at y=1500 m.These different surfaces (bedrock, ice) are depicted in Fig. 5. Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

Simulations and results
Simulations with different degrees of complexity are carried out, that is diagnostic and prognostic simulations, as well as simulations with and without mass balance.These different types of simulations are performed to assess separately the influence of the viscous ice deformation and the effect of mass balance.An overview of the different simulations is given in Table 1.For all the simulations, the same numerical values as given in Table 2 are adopted.They are typical values for an alpine glacier (see for instance Lliboutry, 1965).

Diagnostic simulations
In the diagnostic simulations, velocity distributions of the different models are compared for given and fixed geometries.In all simulations an overestimation of the velocities over the whole ice-body by the SIA model was observed.In most cases, there was a good agreement between HO and FS velocities.These observations are illustrated in the following paragraphs for the different glaciers.

Flattened half-sphere glacier
In Fig. 6, the norm of the surface velocity vector is shown over the symmetry axis for the different models in the case of a bedrock slope of 0.3.An overestimation with the SIA appears clearly whereas the HO model leads to a smaller amplitude underestimation.In the case of the flattened half-sphere, there was no real difference regarding this over(under)estimation seen for different bedrock slopes.Probably this is due to the fact that the aspect ratio of the ice-body itself is already so large that the bedrock slope does not play an important role anymore.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

Conic glacier
In the case of the conic glacier, the dependence of the validity of the SIA on the bedrock slope is clearly observable.In Fig. 7, the ratio between SIA and FS velocities is shown for the radial component with different bedrock slopes (the almost uniform velocities in the region where the ice-thickness is 40 m are here considered).The overestimation of the velocities produced by the SIA model, increases with increasing bedrock slope.

Valley glacier
For this geometry, the over-implicit scheme has been used for the SIA model, allowing a better convergence with bigger time-steps (for diagnostic as well as prognostic simulations).
In the same way, a degradation of the SIA with the aspect ratio of the ice body is detectable in the case of the valley glacier.In Fig. 8, the evolution of the norm of the surface velocity vector along the symmetry axis (x=0) is shown for two different ice-thicknesses.
As for the other geometries, the HO model leads to an underestimation of velocities (hardly visible in the figure due to the scale).
In addition, this type of geometry indicates that models not resolving the FS equations not only disagree in terms of amplitude, but there can be also velocity components that are not distinguished at all with some models.This is illustrated in Fig. 9, where the vertical profiles of the different velocity components in one specific point, not located on the symmetry axis, are shown.The velocity in the x-direction is not distinguished with the SIA model.In the SIA the ice-flow occurs always along the highest slope, the velocity in the x-direction equals zero, as d S/dx=0.This reiterates the weakness of local theories, where stress is not transmitted between grid points.Introduction

Conclusions References
Tables Figures

Back Close
Full No experiences with sliding are presented here.Nevertheless several simulations with different Weertman-type sliding laws (Weertman, 1964) have been conducted.The main result of these simulations is that qualitatively the agreement between the SIA and the FS model is in the most cases better than without sliding.This result is surprising as a degradation of the SIA is expected as the spatial transmission of the stress gradients which is neglected in the SIA becomes more important (Gudmundsson, 2003).
In the chosen simulations the better agreement can be explained by the fact that the total velocity is dominated by the sliding velocity compared to the deformation velocity and very good agreement for the basal velocities is achieved between both models.

Prognostic simulations
In the prognostic simulations we compare the final velocity distributions of the different models as well as the deformation of the surface under the ice-flow and mass balance (mass balance is applied only in Sect.5.3).As in the diagnostic case, in all simulations the SIA model overestimates velocities, but this overestimation decreases with time.
This decrease of the overestimation can be explained by a negative feedback.Initial (diagnostic) velocities are overestimated within the SIA which leads to extra deformation of the ice.This in return gives a thinner and flatter glacier compared with the FS approach which is assumed to be close to reality.As SIA-velocities are proportional to ice-thickness and surface slope (raised to the 4th and 3rd powers, respectively), they rapidly reduce in the following time steps.The glacier slows down until its velocity field is coherent with its geometry.For this geometry, no steady state simulation is possible under zero mass balance.The simulations were thus conducted over a period of 50 years, a period long enough to observe different behaviours of the different models and short enough to deal with reasonable calculation time.In terms of velocity the same observations as in the diagnostic case are made, i.e. overestimation by the SIA model and small underestimation by the HO model.However, the longer the simulation, the better is the agreement between the three models.Figure 10 showing the final vertical profile of the norm of the velocity vector at point (x, y)=(300, −250) for the prognostic simulation (left of figure) as well as for a diagnostic simulation (right of figure) illustrates the better agreement in the case of a prognostic simulation.In both cases velocities are normalized to the corresponding FS surface velocities.The initial ratio between the SIA and the FS velocities decreases from about 4 (diagnostic) to nearly 1 (at the end of the prognostic simulations).
The observations made on the surface deformation are in agreement with those observed with velocities.The SIA glacier is too thin and too long (because has undergone too much deformation), whereas the HO glacier is not deformed enough (Fig. 11, thick lines).

Conic glacier
In the same way as for the flattened half-sphere glacier, the simulation time was fixed to 50 years in the case of simulations where no mass balance is applied.As for the velocities, a degradation of the agreement for the terminus position between the SIA and the FS models is observed with increasing bedrock slopes (Fig. 12 model is rather good.The previously mentioned negative feedback is investigated in more details with this glacier.Figure 13 shows the absolute differences for the surface altitude (dotted line) and the vertical velocity (full line) between the SIA and the FS model as a function of time during the 50 years of simulation.The difference in velocity is decreasing asymptotically to a small value close to zero, whereas the difference for the surface altitude increases asymptotically to a given value.Thus, for a sufficient long simulation, the SIA velocities get close to FS velocities and the discrepancy in terms of ice-thickness stays constant, as the error on the ice-thickness represents the cumulative errors on velocities.
The simulation duration needed to achieve these asymptotic values depends on the bedrock slope, the larger the slope the faster the discrepancy becomes constant (not shown).

Valley glacier
Some simulations with the HO model are missing as the required simulation time appeared too long.Simulations with the valley glacier confirm the observations made with the other glacier types.
Given the negative feedback exposed by the previous simulations, it appeared worthwhile to pursue by changing the simulation time as a function of the initial ice-thickness, i.e. the parameter E 0 .The simulation duration was determined according to a criterion based on surface deformation during the simulation.The different simulations with the parameter E 0 are run until the same amount of deformation is achieved, i.e. when the ice-thickness at the middle of the junction between valley and cirque part (0 m, 5000 m) has reduced to 65% of its initial value.The only exception is for E 0 =100 m for which the simulation was stopped after 100 years because of too small an amount of deformation.
As can be seen from Fig. 14 model eventually appears with large initial thicknesses.At the end of the prognostic simulation, the figure shows a convex shape for the FS model, whereas the surface becomes concave with the SIA.

Prognostic simulations where mass balance is included
Prognostic simulation have also been carried out with prescribed mass balance fields.
In the case of the flattened half-sphere glacier, a spherical mass balance distribution a S was chosen with 5 m of ice in the centre (shifted downhill by 250 m compared to the centre of the ice-body).This distribution remains constant in time.
In the case of the conic glacier, a linear dependency with altitude was chosen with a change in the regression coefficient at the equilibrium line (ELA=5100 m).However for simplicity knowing the relatively thin ice-thickness the mass balance distribution a C was kept constant in time and calculated according to the initial surface altitude.This distribution is motivated by observations made at another nearby volcano glacier (Antizana).
In the case of the valley glacier a linear dependency with altitude was chosen as well for the mass balance distribution a V (with a coefficient of 0.004 m −1 as observed on the Glacier d'Argenti ère in France).
However, the equilibrium line ELA has been artificially shifted well above the glacier surface in order to prevent, in a simple way, unrealistic ice settlements over the bare 577 Introduction

Conclusions References
Tables Figures

Back Close
Full rock slopes around the upper cirque.For convenience the mass balance was calculated once again accordingly to the initial surface altitude.Apart from the case of the valley glacier, all simulations achieve a steady state surface.For the valley glacier the same simulation time as with zero mass balance was chosen because this duration implies the same amount of deformation.
In all simulations a good agreement for the snout position was observed (see Figs. 11 and 12), but not necessarily for the entire glacier geometry.The good agreement for the snout position is explained probably just by the fact that the geometries reduce easily to a 2-D geometry (i.e. a problem with two dimensions in space).In 2-D the snout position in a steady state is fully determined by the mass balance distribution (which is kept constant in all conducted experiences) and the upper border of the glacier, as the total surface mass balance integrated on the glacier surface is equal to zero.
In the case of the valley glacier an interesting point regarding time evolution has been brought up.With the FS model the glacier systematically retreats from the very start.With the SIA model for large thicknesses, the glacier first re-advances due to high deformation rates, and then retreats because of the negative mass balance field.In some situations the overestimation of deformation of the SIA can overcome the effect of a mass balance field and completely change the behaviour.

Comparison of CPU time
As one example for comparison of CPU time, the CPU time necessary for a simulation with the different models was investigated with the flattened half-sphere (slope 0.2).This work does not intend to be an exhaustive study on the CPU requirements with the different models.Instead, the different models as they currently exist are used, i.e. without putting any emphasis on optimisation of the codes, the convergence parameters, initial conditions, time steps, spatial grid steps, etc.Nevertheless, all these parameters have been varied in order to get an idea on the effect on the total CPU time.The aim of this section is thus just to get in a fast and simple way an overall idea of the differences between the models in terms of CPU time.Memory size requirements are not covered 578 Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion in this study.Two different scenarios have been chose; a diagnostic simulation and a prognostic simulation of 10 years without sliding and with zero mass balance.

Diagnostic velocity field
Some interesting observations with the FS model have been made for the CPU time, which (i) does not increase linearly with calculation precision, (ii) increases nearly linearly with the number of grid nodes and (iii) is very sensitive to the initial conditions.In the case of the HO model, CPU time increases more than linearly with the number of grid-points, but initial conditions are less important.
Taking into account the different parameter sets that are not detailed here, a ratio of O(10 000) was found between the FS model and the SIA model, the ratio between the FS model and the HO model was only between O(10) and O(100) and the one between the HO model and the SIA solution was O(100).Specifically, 0.06 s have been needed for the SIA model, between 17 and 28 s for the HO model (depending on different convergence parameters and initial conditions chosen) and between 4 and 18 min for the FS model (also depending on different parameters and initial conditions).

Prognostic simulation
No linear dependency on the time-step was found both for the FS and HO models.Increasing of the time-step by a factor of 2 led to an increase of a factor of 1.3 in calculation time for the FS model and to a factor of 2.3 for the HO model.
In this case the computing time was between 0.22 and 0.39 s for the SIA model, between 30 and 108 min for the HO model and between 2.5 and 4 h for the FS model (in all cases for different parameter sets chosen).
Taking into account the different simulations, the ratio for the CPU time between the FS model and the HO model was now less than 10, the ratio between FS and SIA 579 Introduction

Conclusions References
Tables Figures

Back Close
Full remains O(10 000) and the one between the HO model and the SIA increases to nearly O(10 000) as well.
The large difference between the results on diagnostic and prognostic simulations is explained by the high sensitivity of the FS model to initial conditions.The calculation of the diagnostic velocity field takes disproportionally more time compared to the next time steps; whereas the HO model takes nearly the same time for all time steps -the initial one included.This explains the different ratios for the two types of simulations.
Finally we can conclude that the saving of CPU time with an SIA model is important, but taking a HO model instead of a FS model will not necessarily bring an advantage in CPU time.The latter case has to be investigated further for each given application.

Conclusions
The inter-comparison here uses three different models, a Shallow Ice Approximation (SIA) model, a higher-order (HO) model and a full-Stokes (FS) one.Different synthetic shaped glaciers have been used with different aspect ratios and bedrock slopes.
SIA models are frequently used for modeling large ice masses because of their huge advantages in terms of CPU time and implementation/development facilities.However, this approximation is only valid for ice bodies with a small aspect ratio and a weak bedrock slope.In the SIA, longitudinal stresses (σ x,y,z ) are neglected which prevents stress transmission within grid points and therefore makes the theory local.Glacier evolution at a given point of the glacier thus only depends upon the local topography without including any interactions from the neighboring points.This approximation is certainly not very realistic, especially in the case of a non regular surface or bed topography.Horizontal shear stress (τ xy ) is neglected as well in this theory, which is especially not justified in the case of a glacier deeply channelled in a valley, where the importance of such stresses is obvious and for example clearly visible with the bending of the Forbes bands.
The main conclusion is that deformation is clearly over estimated by the SIA.This Introduction

Conclusions References
Tables Figures

Back Close
Full can be observed on velocity fields as well as on evolved glacier geometries.This overestimation is partly a consequence of the local character of the model due to the lack of longitudinal coupling.In contrast, the HO model leads to a small underestimation.It is also obvious that the smaller the aspect ratio, the better the SIA results as for instance also shown in Le Meur et al. (2004).
More specifically, the following observations are made with the different simulation types: in the diagnostic simulations a huge disagreement between the SIA and the other models has been observed.This disagreement increases with the aspect ratio (valley glacier) and the bedrock slope (conic glacier).In the prognostic simulations the agreement for velocities is much better because of a negative feedback obtained with the SIA.Nevertheless, the overestimation of deformation is still noticeable from the time-dependent geometry evolution.
Mass-balance essentially affects the snout position and a very good agreement between the different models was obtained for this observable, however, a disagreement in terms of surface geometry remains.In some cases, even the glacier behaviour can change, e.g. an initial advance rapidly followed by a retreat with the SIA model contrary to the FS model which continuously retreats from the start.
In terms of comparison of CPU time, the advantage of the SIA model has clearly been shown.But the switch between a HO and a FS model does not necessarily bring a big advantage in term of calculation time.Thus, the SIA seems to be a good approximation in the case of prognostic simulations with small aspect ratios and small bedrock slopes.The mass balance does not have a real effect on this issue but can nevertheless affect the global behaviour.
As a next step, it would be interesting to make some more comparisons on prognostic simulations, not only by comparing final velocities and surfaces, but by comparing the dynamics over the entire simulation time.
remarks on the influence of sliding , bottom part of the figure).The HO model simulations are sometimes missing because of problems for the model to converge.Nevertheless, when available, the agreement with the FS , one additional observation can be made with this kind of glacier shape.Not only velocities and deformation are regularly overestimated over the whole ice-body, but a change in the transverse surface between the FS and SIA 576

Fig. 11 .Fig. 12 .Fig. 13 .
Fig. 11.Final glacier surface for prognostic simulations without (thick lines) and with (thin lines) mass balance.In the case of the simulation with no mass balance, the simulation ended after 50 years and in the case of the simulation with mass balance, steady state was achieved (after about 100 years).The mass balance function is detailed in the text.

Fig. 14 .
Fig. 14.Valley glacier, final ice surfaces after 100 years (left) and 1.4 years (right) of simulation for initial ice-thicknesses E 0 of respectively 100 m and 700 m.Inlets represents zooms on the surface with appropriate modification of the vertical scales for a better display.
have shown that Figures

Table 1 .
Overview of simulations.D (P) stands for Diagnostic (Prognostic).Simulations with are missing for the HO model.

Table 2 .
Numerical values of the parameters adopted for the all the simulations.