1 Introduction
Numerical modelling of glaciers and ice sheets generally involves a number of simplifications with respect to the physics of the ice mass. A common simplification is the so-called shallow-ice approximation (Reference HutterHutter, 1983), which applies to ice masses of a small aspect ratio (characteristic ice thickness divided by the characteristic horizontal extension) and characterized by relatively small surface and bedrock slopes. Ice sheets belong to this category. For glaciers as well as key areas in ice sheets, such as ice divides, grounding zones and areas with a pronounced basal topographic relief, the shallow-ice approximation does not hold anymore. A rigorous treatment of the force balance in an ice mass (the so-called higher-order solution) is less straightforward, but a number of numerical schemes have already been developed to calculate the velocity and stress field according to higher-order approximations. Among these schemes are those of Reference Herterich, Van derVeen and OerlemansHerterich (1987), Reference Dahl-JensenDahl-Jensen (1989), Reference Van der Veen and WhillansVan der Veen and Whillans (1989), Reference BlatterBlatter (1995) and Reference Colinge and RappazColinge and Rappaz (1999).
Few higher-order models include the time dependence, i.e. the evolution of the geometry of the glacier as a reaction to surface boundary conditions such as surface mass balance and its perturbations. Amongst these are the finite-element model of Reference HvidbergHvidberg (1996) and the finite-difference model of Reference MayerMayer (1996). Since those models already demand a significant computational effort, a time dependence, which iterates the aforementioned scheme, makes a proper solution less straightforward.
The higher-order model presented here is a dynamic flowline model that predicts the ice-thickness distribution along a fixed flowline in space in response to environmental conditions. This response is obtained by calculating at a given moment the two-dimensional flow regime and the temperature distribution, determined by the glacier geometry and its boundary conditions. A previously developed flowline model described in Reference PattynPattyn (1996), which is also a higher-order model, was based on the model of Reference Van der VeenVan der Veen (1989), i.e. by integrating the field equations from the base of the ice mass to the surface. However, that model still suffered from numerical instabilities when the grid size was reduced, and remained time-consuming as well. Reference BlatterBlatter (1995) applied a similar technique to solve the field equations in glaciers, i.e. by integrating them from the base of the glacier to the surface, but using the method of lines. The boundary-value problem is solved iteratively, starting from zeroth-order estimates of the basal shear stress. His method was further improved using a single-shooting Newton iteration (Reference Colinge and BlatterColinge and Blatter, 1998).
Another approach for computing the velocity field in glaciers is based on the method presented by Reference Herterich, Van derVeen and OerlemansHerterich (1987), i.e. by formulating the field equations as a function of the horizontal velocity. A finite-element scheme according to this approach is presented by Reference Colinge and RappazColinge and Rappaz (1999) . The model presented here follows a similar approach: the momentum equations are reformulated as a function of velocities by introducing the constitutive equation so that a differential equation as a function of the horizontal velocity is obtained, which is solved by means of the conjugate gradient method. A further series of iterations is then necessary to determine the velocity field accurately. Compared to the model of Reference PattynPattyn (1996), transverse strain rates and stresses are now included properly to account for convergence and divergence of the ice flow (e.g. Reference ReehReeh, 1988). The model is compared to results from Reference Blatter, Clarke and ColingeBlatter and others (1998) from Haut Glacier d’Arolla, Switzerland. It is further compared with the Van der Veen method (Reference PattynPattyn, 1996) and the shallow-ice approximation according to the specifications in the European Ice-Sheet Modelling Initiative (EISMINT) benchmark for ice-sheet models (Reference Huybrechts and PayneHuybrechts and others, 1996) . Finally, it is investigated how a glacier, simulated with a higher-order model, reacts to sudden and/or continuous perturbations in surface mass balance.
2 Field Equations
Consider a Cartesian coordinate system (x, z) with the x axis along the flowline, parallel to the geoid, and the z axis pointing vertically upward (z = 0 at sea level). The only non-zero velocity components are thus u, w (horizontal and vertical velocity, respectively), while the horizontal transverse velocity v = 0. Since the flowline geometry considers transverse velocity gradients, it follows that ∂v/∂y = (u/ω)(∂ω/∂x) where ω is the width of the drainage basin taken perpendicular to the flowline. It should be noted, however, that this assumption is valid for the axisymmetric case or when the flowline exhibits a low curvature. The equations for conservation of mass and linear momentum then become
where ρ is the ice density, g the gravitational constant, and τij the stress components. The values for these and all other constants used in this paper are given in Table 1. The parameter ∊ 2 adds complexity to the model: for ∊ 2 = 1, vertical resistive stresses and their gradients are included, so that the full system of equations is solved. This model will be further referred to as the full-system model. Ignoring this term (∊ 2 = 0) corresponds to the model description of Reference BlatterBlatter (1995), further referred to as the basic model. A recent review of asymptotic theories of large-scale motion in ice sheets by Reference Baral, Hutter and GreveBaral and others (2001) showed that the basic model description corresponds to an incomplete second-order approximation.
Neglecting atmospheric pressure, an expression for τzz is obtained by integrating Equation (3); from the surface s to a height z in the ice body
The second term of the righthand side of Equation (4) represents the vertical resistive stress Rzz, which is responsible for the so-called bridging effect (Reference Van der Veen and WhillansVan der Veen and Whillans 1989). Deviatoric stresses are related to the stress components by 
for i = x, y, z, and σxz = τxz. Making use of Equation (4) and inserting the deviatoric stresses, Equation (2) becomes
The constitutive equation governing the creep of polycrystalline ice and relating the deviatoric stresses to the strain rates is taken as a Glen-type flow law with exponent n = 3 (Reference PatersonPaterson, 1994),
where 
is the second invariant of the strain-rate tensor, defined by 
, and μ is the effective viscosity. 
is a small number (10−30 a−1) to validate Glen’s flow law in cases where 
equals zero, and singularity might occur. The use of such a small number does not influence the numerical outcome of the model. The flow-law rate factor A(T
*) is a function of temperature, where T
* is the ice temperature (K) corrected for pressure melting, i.e. T
* = T + βP, where P is the ice pressure (P = σxx
+ σyy + ∊
2
Rzz
− βg(s − z)) Following Reference HookeHooke (1981), A(T
*) is set to
Here Q is the activation energy for ice creep, R is the universal gas constant and m is an enhancement factor (or tuning parameter). Experiments carried out with the widely used “two-step” Arrhenius relationship to determine the flow-law rate factor A(T*) (e.g. Reference Payne and DongelmansPayne and Dougelmans, 1997) failed, probably due to this discontinuity. Iteration methods usually rely explicitly or implicitly on gradients being smooth as well as the function, which is not the case with a “two-step” Arrhenius relationship. However, using the Arrhenius equation (Equation (7)), which is continuous over the whole range of ice temperatures, allows for a proper solution of the iteration scheme.
By definition, strain rates are related to velocity gradients by 
Applied to the flowline geometry and making use of the conservation of mass (Equation (1)), one obtains
Combining the flow-law equation (Equation (6)) and the horizontal stress-field equation (Equation (5)), and replacing the strain-rate components by velocity gradients using Equations (8–11), one obtains
Where
Expanding Equation (12) and rearranging terms leads to
An expression for the vertical velocity w is obtained through vertical integration of the incompressibility condition (Equation (1)) from the base of the glacier to a height z,
The thermodynamic equation for the ice flow along a flow-line is written as
where k i and c p are the thermal conductivity and the specific heat capacity, respectively, and τ is the effective stress or the second invariant of the stress tensor. The heat transfer is thus a result of horizontal and vertical diffusion, horizontal and vertical advection, and internal friction due to deformational heating.
Using a kinematic boundary condition at the upper and lower surface of the ice mass (see below), the mass-balance equation (Equation (1)) is integrated along the vertical in order to obtain an expression for the change of local ice thickness in space
where 
is the depth-averaged horizontal velocity (m a−1), H is the ice thickness (m), Ms
is the surface mass balance (m a−1 ice equivalent) and M
b is the melting rate at the base of the glacier (m a−1).
In the experiments described below, a comparison is made with a zeroth-order model (according to the shallow-ice approximation). Ignoring terms involving ∊ 2 as well as normal stress components (τxx) in Equation (2), it is straightforward to derive the velocity u(z), which is then solely dependent on the local geometry, i.e. ice thickness and surface gradient. The horizontal velocity in the grounded ice sheet according to the shallow-ice approximation is determined as (e.g. Reference Huybrechts and OerlemansHuybrechts and Oerlemans, 1988)
where u b is the basal velocity.
3 Boundary Conditions
Boundary conditions to the ice mass are a symmetric ice divide (∂s/∂x = ∂H/∂x = 0) in the case of an ice sheet, or H = 0 in the case of a valley glacier, zero ice thickness at the downstream edge of the glacier, a stress-free surface and a basal sliding function at the bottom of the glacier, which is a function of the basal drag τ b,
where A s = 2.0 × 10−7 N−2 m5 a−1 is a sliding parameter and N the effective pressure at the base of the ice mass (Reference Van der Veen, Van derVeen and OerlemansVan der Veen, 1987). The kinematic boundary condition at the lower ice surface reads
The basal drag is defined as the sum of all basal resistive forces and written as (Reference Van der Veen and WhillansVan der Veen and Whillans, 1989)
It should be noted that the basal drag is linked to the driving stress τ d = −ρgH∇s, corrected for longitudinal pushes and pulls and the vertical resistive stress. This expression is obtained after vertically integrating Equation (5) from the base of the ice mass to the surface:
The stress-free surface implies that
where Rzz (s) = τxz (s)∂s/∂x. Written in terms of velocity gradients, this results in
Boundary conditions to the thermodynamic equation (Equation (16)) follow from the mean annual air temperature at the surface. At the base, the temperature gradient is defined as
where G is the geothermal heat flux. However, a constant geothermal heat flux is not strictly correct; it might be more realistic to consider heat conduction in the bedrock, so that G = −kr(∂T/∂z), where k r is the thermal conductivity in the underlying bedrock. In view of the experiments carried out hereafter, the choice of basal thermal boundary condition will not alter the results. Furthermore, a constant value of G is a common boundary condition in ice-sheet and glacier models, and facilitates future model intercomparison. The basal temperature in the ice mass is kept at the pressure-melting point whenever it is reached, and the basal melt rate M b is calculated as
where L is the specific latent heat of fusion.
4 Numerical Solution
For numerical convenience, a dimensionless vertical coordinate is introduced to account for ice-thickness variations along the flowline, and which is defined as ζ = (s − z)/H, so that ζ = 0 at the upper surface and ζ = 1 at the base of the ice mass. The coordinate transformation maps (x, z) → (ξ, ζ), where ξ ≡ x is the transformed horizontal coordinate. Details on the transformation of the global velocity term (Equation (14)), the viscosity term (Equation (13)) and the surface boundary condition (Equation (24)) are given in Appendix A. After this transformation, all equations are solved on an irregularly spaced grid in ξ and ζ. Central differences are used to compute first- and second-order gradients. At the boundaries, upstream and downstream differences are employed (see Appendix B for more details on the procedure).
41 Velocity field
The finite-difference forms of the transformed differential equations for u (Equations (52) and (53)) are written as a set of linear equations with u(ξ, ζ) as unknowns. In matrix notation this becomes
where ℓ is the iteration number. Starting from a zeroth-order estimate of the horizontal velocity field u 0, a new estimation u 1 of the velocity is obtained by solving the set of linear equations. For a numerical grid of Nξ = 17, Nζ = 21, the matrix A consists of (17 × 21)2 = 127 449 nodes. However, only 2642 nodes contain non-zero elements, clustered around the diagonal, which is 2% of A.
A solution to the linear system of Equation (27) is accomplished using the sparse-matrix algorithms of Reference Press, Teukolsky, Vetterling and FlanneryPress and others (1992), which are based on the conjugate gradient method. Although the coding of sparse matrices is rather complicated, they are far more efficient in terms of computation time than point-relaxation algorithms on the full or even parts of the matrix. Because of the non-linear nature of Equation (52), A and b contain three parameters that are still a function of u, i.e. the viscosity term μ, the vertical velocity w, and the vertical resistive stress Rzz, which have to be determined in an iterative fashion (the iterative determination of both w and Rzz applies only to the full-system model for ∊2 = 1). The successive substitution method or Picard iteration was used for this purpose. In order to optimize the rate of convergence, a relaxation formula was added based on the unstable manifold correction (Reference Hindmarsh and PayneHindmarsh and Payne, 1996). Therefore, Equation (27) is written as
where u * is the velocity estimate obtained with the conjugate gradient method. Consider an iterative solution of a non-linear equation which generates a series of approximate solutions u ℓ+1, u ℓ,…, being updated by a series of correction vectors c ℓ+1, c ℓ,…, such that u ℓ+1 = u ℓ + c ℓ Since the correction vector is defined as c ℓ = u * − u ℓ, the Picard iteration would simply update the velocity u ℓ+1 with this correction vector, so that u ℓ+1 = u *. If eℓ+1, eℓ,… is taken as the error in the solution vector u ℓ+1, u ℓ,…, then we can state that (eℓ+1, eℓ,…) = α(c ℓ+1, c ℓ,…). Assuming that the decay is on a straight line in the correction space, we obtain (Reference Hindmarsh and PayneHindmarsh and Payne, 1996)
A new update of the velocity vector is obtained by
where the modified correction vector becomes 
. The direction θ between successive correction vectors is computed as
where the norms refer to the L 2 norm. Whenever this angle is close to 0 or π the subspace iteration might become applicable. Using this subspace relaxation algorithm, the solution vector properly converges, which is not the case with the Picard iteration method.
42 Glacier evolution and thermodynamics
The continuity equation (Equation (17)) is reformulated as a diffusion equation for ice thickness H, i.e.
where 
. Equation (30) results in a tridiagonal system of Nx
equations, and is solved using the tridiagonal algorithm of Reference Press, Teukolsky, Vetterling and FlanneryPress and others (1992).
The thermodynamic equation (Equation (54)) is solved implicitly in the vertical, giving rise to a tridiagonal system of Nζ equations which is solved using the tridiagonal algorithm of Reference Press, Teukolsky, Vetterling and FlanneryPress and others (1992). A two-point upstream difference notation was employed for the horizontal, while central differences were used in the vertical. The horizontal implicit terms are found by iteration of this scheme. Only a few iterations are necessary to obtain a good convergence.
5 Model Comparison and Validation
51 Haut Glacier d’Arolla
One way to validate the sequence of numerical schemes is to compare the model — or part of the model — with other studies. For this purpose, the longitudinal surface and bedrock profiles of Haut Glacier d’Arolla were digitized from Reference Blatter, Clarke and ColingeBlatter and others (1998), in order to compare our basic model (∊ 2 = 0) with the Blatter model (Fig. 1a). The longitudinal profile of this glacier has a very simple geometry, hence the resulting stress field is not influenced by geometrical perturbations such as the presence of a steep icefall. For this model run we keep the glacier geometry fixed. A zero basal velocity was considered, and ω was kept equal to 1 along the whole flowline domain. The flow-law rate factor A(T*) was taken constant over the whole model domain, and equals A(T *) = 10−16 Pa −n a−1, a value corresponding to an isotherm ice mass of around −2°C (Reference PatersonPaterson, 1994). The model was numerically solved on a regular grid in both the horizontal and vertical dimensions, with 50 equidistant layers in the vertical. A series of experiments were carried out, each with different grid resolutions, i.e. 380, 200, 100, 50, 20 and 10 m, respectively. To obtain a glacier geometry at these different resolutions, the longitudinal bedrock and surface profiles were resampled using a cubic spline to assure smooth surface slopes even at the highest model resolution. The resulting basal drag and driving stress for the 20 m resolution experiment are displayed in Figure 1c, and can be compared to those of Reference Blatter, Clarke and ColingeBlatter and others (1998, p.463, fig. 13). Small differences in driving stress compared to the results of Reference Blatter, Clarke and ColingeBlatter and others (1998) are due to the accuracy of the digitizing process, since the driving stress depends solely on the local glacier geometry. For the model results, however, the comparison between the driving stress and the basal drag is more relevant.
Fig. 1 (a) Longitudinal profile of Haut Glacier d’Arolla, taken from Reference Blatter, Clarke and ColingeBlatter and others (1998) (solid line ).The dashed line is the calculated steady-state surface profile. The glacier is in both cases 3.8 km long (b) The difference between the basal drag and the driving stress τ
b
− τ
d
(solid line) vs the vertically integrated longitudinal stress gradient 
(dotted line). (c) Basal drag τ
b
(solid line), basal shear stress τ
xz
(dotted line) and driving stress τd (dashed line). The horizontal model resolution is 20 m.
The basal drag is a smoothed version of the driving stress (Fig. 1c), and displays a pattern that is very similar to the profile obtained by Reference Blatter, Clarke and ColingeBlatter and others (1998). The driving stress oscillates between 0.7 and 2.2 bar, and becomes zero at both ends of the longitudinal profile. Differences between the τ d and τ b profiles are as high as 0.5 bar (25% of τ d). The basal-drag signal is thus characterized by a lower amplitude and phase shifts compared to the profile of the driving stress. The direction and magnitude of these phase shifts are similar to those observed in Blatter’s Arolla profile. Such phase shifts are due to local pushes and pulls as a result of bedrock and surface perturbations. Reference Mayer and HuybrechtsMayer and Huybrechts (1999) identified phase shifts between π/4 and π/2 on Ekströmisen, Antarctica, suggestive of a topographical control on surface undulations. However, it will be shown below that the differences between driving stress and basal drag — in both magnitude and phase shift — gradually disappear when the free surface is allowed to react on surface mass-balance conditions and a steady-state glacier profile is reached.
The model consistency is verified by Equation (22), which states that the basal drag, as defined in Equation (21), must balance the difference between the driving stress and the gradient of the vertical mean longitudinal stress deviator (for ∊2 = 0). Both curves, i.e. τ b calculated from Equations (21) and (22), respectively, coincide well for the Arolla experiment (Fig. 1b). A second verification involves the integration of Equation (22) along the whole flowline. For the basic model this becomes
Since the ice thickness H as well as the stress deviators equal zero at both ends of the flowline (at x = 0 and x = L), it follows from Equation (31) that 
. Hence, the average of the basal drag over the entire bed equals the average of the basal driving stress over the entire bed and thus is an invariant that only depends on the geometry of the glacier. This relation therefore yields an accuracy test for the numerical solution (Reference Blatter, Clarke and ColingeBlatter and others, 1998). It can be seen that this invariant depends almost linearly on the model resolution, which confirms Blatter’s findings (Fig. 2). Moreover, the invariant 
holds for all resolutions, contrary to results of Reference Blatter, Clarke and ColingeBlatter and others (1998). The difference between the mean basal drag and the mean basal driving stress, especially at larger Δx, is due to a round-off error specified by the iteration scheme. This error is < 1% at grid sizes smaller than 100 m. According to Blatter’s findings, however, the basal drag remains constant over the whole range of grid resolutions (Reference Blatter, Clarke and ColingeBlatter and others, 1998). This discrepancy might be due to the difference in resampling of the glacier geometry at different grid resolutions.
Fig. 2 Effect of grid resolution on the calculated longitudinal average of the basal drag 
(solid line) and of the driving stress 
(dashed line) for the profile of Haut Glacier d’Arolla.
52 EISMINT benchmarks
The computer code of the zeroth-order model (according to the shallow-ice approximation), the basic model (∊ 2 = 0) and the full-system model (∊ 2 = 1) is tested against the EISMINT benchmarks for ice sheets with simplified geometry (Reference Huybrechts and PayneHuybrechts and others, 1996). Furthermore, the full-system model is compared to the force-budget method of Reference Van der VeenVan der Veen (1989) and Reference Van der Veen and WhillansVan der Veen and Whillans (1989), which was previously used in Reference PattynPattyn (1996). The EISMINT moving-margin experiment was chosen for the intercomparison of the different model codes, especially in relation to the time-dependent solution. Since for these experiments aspect ratios are small, all models should more or less give the same result. Moreover, the full-system model should give exactly the same result as the Van der Veen model of Reference PattynPattyn (1996).
The EISMINT benchmark consists of a flowline of 16 gridpoints, regularly spaced at intervals of 50 km, so that the model domain stretches from x = 0 (a symmetric ice divide) to x = 750 km (zero ice thickness). In the vertical, 50 equally spaced layers are considered. In order to produce an axisymmetric ice-sheet profile, the width parameter ω was chosen to vary linearly along the flowline, with ω = 0 at the ice divide and ω = 750 km at the last gridpoint. The absolute value of ω is not important, as long as its gradient is correct, which is in this case a constant value over the whole model domain. The climatic boundary conditions at the surface are determined as (Reference Huybrechts and PayneHuybrechts and others, 1996)
where d is the distance (in km) from the ice divide, a = 10−2 m a−1 km−1 is the mass-balance gradient, and R el = 450 km is the distance at which the surface mass balance changes from positive to negative values.
A calculation over a time period of 50 000 years with a time-step of 20 years takes 20.98 central processing unit (CPU) seconds on a Compaq-Alpha server GS140 for the basic model solution starting from zero ice thickness over the whole model domain, and 8.32 CPU seconds when starting from an already calculated ice-sheet profile. For the zeroth-order model it takes about 2.87 CPU seconds in both cases. The evolution of the ice surface as well as the temperature field is exactly the same for the full-system model solution and the Van der Veen model. Both deviate slightly from the shallow-ice approximation, but this deviation, except for the ice divide, is only 1–2% of the ice thickness and global velocity field, as is to be expected for such a large-scale ice-sheet experiment. This confirms the validity of the shallow-ice approximation in large ice sheets resting on a flat bedrock. At the ice divide, the difference in ice thickness with the shallow-ice approximation is somewhat larger (4–5 m), as longitudinal stress gradients are not neglected in the basic model.
6 Glacier Evolution
61 Steady-state solutions
Starting from the present longitudinal profile of Haut Glacier d’Arolla (Fig. 1a, solid line), the glacier is now allowed to react to a given time-independent surface mass-balance distribution defined as
where a = 10−2 a−1 is the mass-balance gradient, estimated from Reference Willis, Arnold, Sharp, Bonvin, Hubbard, Lane, Chandler and RichardsWillis and others (1996), and ELA = 2800 m the equilibrium-line altitude, chosen in such a way that the steady-state front position coincides with the present observed position of the glacier terminus (Fig. 1a, dashed line). The width of the drainage basin ω was also determined from the maps by Reference Willis, Arnold, Sharp, Bonvin, Hubbard, Lane, Chandler and RichardsWillis and others (1996). Figure 3 displays the basal stress field and surface conditions of Haut Glacier d’Arolla in steady state. Here, the relation between the driving stress and the basal drag is quite striking. Compared to the fixed-geometry calculations of Figure 1, the large oscillations in the driving stress are apparently smoothed out when the free surface is allowed to react to mass-balance conditions at the surface. Furthermore, the basal-drag profile is almost similar to the driving stress, which means that longitudinal pushes and pulls tend to cancel themselves out with surface and bedrock slopes. Not only does the basal drag tend to equal the driving stress in amplitude, but the phase shift between both disappears as well. The reason for this behaviour lies in the surface slopes (Fig. 3a), since they readjust to the stress field in two ways. First, the steady-state surface slopes are smoother than the observed ones, except at the beginning and end of the longitudinal glacier profile. Second, a phase shift appears in the surface slope profiles when the basic model is compared with the zeroth-order model. This difference between the steady-state surface profile and the observed profile might be attributed to the fact that the present glacier surface is not in steady state, but also to the fact that the applied surface mass-balance function is a smoothed function and might deviate from reality, and to assumptions in the model description that are different from the real world (isothermal ice, no basal sliding, etc.). It should be added that the longitudinal profile of Haut Glacier d’Arolla has a relatively simple geometry and lacks, for instance, a pronounced icefall, where compression and extension of ice flow would eventually alter the basal-drag profile compared to the driving stress.
Fig. 3 (a) Surface slopes along the flowline of Haut Glacier d’Arolla of the steady-state profile according to the basic model solution (solid line), according to the shallow-ice approximation (dotted line), and according to the present observations in Reference Blatter, Clarke and ColingeBlatter and others (1998) (dashed line). (b) Basal drag τ b (solid line), basal shear stress τ xz (dotted line), driving stress τ d (long-dashed line) and basal longitudinal stress σ xx (short-dashed line). The horizontal model resolution is 50 m.
62 Response to perturbations at the surface
Glaciers advance and retreat in response to changes in mass balance induced by climate. Such front migrations do not occur instantaneously, but with a given time lag in response to the climatic perturbation. The time lag is due to the downstream mass transport of the glacier and governed by the deformational and flow properties of the glacier. A theoretical concept that translates this mass transfer is found in the theory of kinematic waves, and was applied to glaciers by Reference NyeNye (1960). Subsequent validation and comparison of the kinematic wave theory with numerical glacier models has been done by Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others (1989) and Reference Van de Wal and OerlemansVan de Wal and Oerlemans (1995). In this paper, we will limit our investigation to whether the model complexity is responsible for different response patterns to sudden perturbations in the surface mass balance.
Starting from the steady-state glacier conditions, a sudden change in surface mass balance was applied, i.e. a perturbation of +0.5 m in accumulation over the entire glacier surface, lasting for 1 year. The glacier reacts to such a perturbation by transporting the excess mass from the head of the glacier towards the front, where the ice accumulates and the front considerably thickens. After the perturbation, the ice thickness starts to decrease in the accumulation area and increases in the ablation zone. Such a behaviour is expected from the linear-wave theory and has been demonstrated by numerical model experiments (Reference Van de Wal and OerlemansVan de Wal and Oerlemans, 1995). Our experiments confirm this behaviour, though the basic (higherorder) model transports the ice faster towards the front, so that the ice thickness at the terminus is higher compared to the zeroth-order model (Fig. 4). The pronounced frontal thickening is also confined to a smaller area. Due to the accelerated downstream transport of ice, a steady-state condition after the perturbation is reached sooner than with the zeroth-order model (Fig. 4). For both models the terminus position did not change throughout the experiment.
Fig. 4 Change in ice thickness after a sudden surface mass-balance increase of 0.5 m a−1 along the whole flowline, lasting lyear. (a) and (c) show the results for the basic higher-order model, and (b) and (d) for the shallow-ice approximation. (a, b) Response after lyear (solid line), 5 years (long-dashed line) 10 years (short-dashed line) and 20 years (dotted line). (c, d) Response after 20 years (solid line), 50 years (long-dashed line), 100 years (short-dashed line) and 150 years (dotted line).
A better insight into these mechanisms is obtained by forcing the model with a sine-wave ice-thickness perturbation in the accumulation area:
where H 1 is the new ice thickness and H 0 refers to the ice thickness before the perturbation was applied. The perturbation in ice thickness starts at the head of the glacier. At t = 0, a local maximum in ice-thickness perturbation is thus situated at x = λ/2. The downstream propagation of the wave is monitored by determining for each position along the flowline the time needed for the maximum ice-thickness perturbation to arrive, which is a measure for the phase velocity of the wave. The results for A = 1, 2, 4 and 6 m, respectively, obtained by both model configurations are shown in Figure 5. λ was taken 1000 m, thus covering the accumulation area at t = 0, and which corresponds to 25% of the total glacier length. According to the shallow-ice approximation, the wave travels almost at constant speed towards the front (dashed lines in Fig. 5). The situation is slightly different for the basic model (solid lines in Fig. 5), where the wave travels more slowly in the accumulation area, accelerates in the ablation area and reaches the terminus before the zeroth-order model. Due to this acceleration, and confirming the findings of the previous experiment, more mass is stored near the front, so that at higher amplitudes (A = 6 m) the model glacier is more sensitive to a front advance, contrary to the shallow-ice approximation. This front migration therefore disturbs the velocity signature for A = 6 m in Figure 5, where the wave is retarded as mass is spread out in the horizontal and not only in the vertical.
Fig. 5 Time needed for an ice-thickness perturbation (as given by Equation (34)) to travel down the glacier The curve shows the time of the maximum ice thickness reached at a specific site on the glacier A denotes the amplitude of the perturbation. Results given for the basic model (solid line) and the shallow-ice approximation (dashed line). The irregular nature of some of the curves is due to the interpolation procedure which locates the maximum ice-thickness perturbation in the (discrete) output files. They are by no means numerical artefacts produced by the model itself.
A more detailed view of the glacier behaviour is given in Figure 6 for the experiments with A = 4 and 6. Here, the local ice-thickness variation at four positions along the flowline is given (x = 0.53, 0.80, 0.96 and 1.00 scaled distance). Although the perturbation travels more slowly in the upstream part according to the basic higher-order model, a steady state is reached sooner compared to the shallow-ice approximation. Close to the ice front, we notice that the response signature is disturbed due to the non-linearity of the model equations. A similar behaviour is noticed in other experiments as well. With higher values for the amplitude A, for instance, the pattern becomes even more complex. Altering the grid-size resolution or the error cut-off value in the subspace relaxation algorithm leads to the same result, which means that the observed behaviour is not a numerical artifact.
Fig. 6 Response of the glacier surface at particular locations (x = 0.53, 0.80, 0.96 and 1.00) to a sudden perturbation (as given by Equation (34)) for the basic model (A = 4.0 m (thick solid lines) and A = 6.0 m (long-dashed lines)) and according to the shallow-ice approximation (A = 4.0 m (thin solid lines) and A = 6.0 m (short-dashed lines)). The irregular nature of some of the curves is explained in the text.
The above experiment was repeated for waves with different wavelengths as well as with a glacier of simplified geometry. Those experiments confirm the above-described glacier response.
63 Response to climatic variations
The previous series of experiments demonstrated how a sudden increase in ice accumulation at the head of the glacier affects the response of this ice mass. The response is generally fast and, with one exception, does not influence the length of the glacier. When studying glaciers over a longer time-span of several hundreds of years, significant glacier-length variations do occur. This series of experiments is designed to investigate the reaction of a glacier described by a higherorder model to continuous perturbations in glacier mass balance on longer time-scales, which involves a distinct pattern in glacier advance and retreat. We applied, therefore, a sinusoidal mass-balance forcing in time:
where T is the time period, t time (years) and A the mass-balance perturbation amplitude. The reaction of the glacier to the mass-balance perturbation (Equation (35)) can be expressed in a variety of measures. One measure is the response time lag to the climatic perturbation, obtained from a correlation between the forcing series (M s (x, t)) and the response time series (H(x, t)) at different time lags. The lag corresponding to the maximum correlation then gives the time lag between both time series (Fig. 7). A second measure is the response range, defined as the maximum minus the minimum ice thickness during the investigation period. Starting from steady-state conditions, the model of Haut Glacier d’Arolla was forced according to Equation (35) with A = 0.1 m a−1 and a time period T = 500 years, for a total period of 2000 years. This experiment was repeated with the shallow-ice approximation as well. The lagged response of the local ice thickness follows approximately 40 years (in the accumulation zone) to 90 years (near the terminus) after the mass-balance perturbation. The difference between both model types is largest in the central part of the glacier and is of the order of 5 years (a 10% difference between the shallow-ice approximation and the basic model, the latter showing a larger time lag). However, both advance and retreat of the glacier occur faster with the higher-order model, but a local maximum ice thickness is reached at a later moment in time. Variations in range also increase towards the terminus. For this parameter we did not observe a significant difference between both model types.
Fig. 7 Glacier response time lag (a) and range (b) to a perturbation in surface mass balance according to Equation (35) (A = 0.1 m a−1;. T = 500 years). Response curves are given for the basic higher-order model without (solid line) and with basal sliding (long-dashed line) and for the shallow-ice approximation without (short-dashed line) and with (dotted line) basal sliding.
A second experiment consisted of introducing basal sliding over the whole model domain, according to Equation (19). For the zeroth-order model, the basal drag was replaced by the driving stress τ d. Despite significant differences in model description, the overall behaviour is in accordance with the previous non-sliding experiments (Fig. 7). Moreover, the difference in response pattern due to the introduction of basal sliding is of minor importance compared to the use of higherorder physics (2% vs the aforementioned 10%). One reason might be that the basal-drag pattern is rather similar to the driving-stress pattern in steady-state conditions (Fig. 3). Since these stresses directly influence the basal sliding conditions, both models should therefore react in a similar way.
7 Conclusions
In this paper, a new higher-order numerical glacier model is developed, which is numerically stable and fast and copes with small horizontal grid sizes. Moreover, the possibility of using an irregularly spaced grid allows for high-resolution modelling adapted to the available input-data distribution. The numerical model is compared with (i) the Blatter-model solution on the longitudinal profile of Haut Glacier d’Arolla (Reference Blatter, Clarke and ColingeBlatter and others, 1998), (ii) the Van der Veen model (Reference Van der VeenVan der Veen, 1989; Reference PattynPattyn, 1996) and (iii) the shallow-ice approximation of the EISMINT benchmarks (Reference Huybrechts and PayneHuybrechts and others, 1996).
The calculation of the stress field in the longitudinal profile of Haut Glacier d’Arolla (fixed geometry experiment) demonstrates that the basal-drag profile has a smaller amplitude than the driving-stress profile and that a pronounced phase shift occurs between the two curves. Furthermore, the average of the basal drag over the entire bed equals the average of the driving stress over the entire bed 
and thus is an invariant that solely depends on the glacier geometry. This invariant holds for any grid resolution.
When the glacier surface is allowed to evolve in response to a surface mass-balance distribution, the steady-state basal-drag and driving-stress profiles tend to fall together. The phase shift, necessary to balance the longitudinal stresses, disappears in the stress field but becomes apparent in the surface gradients. This might indicate that the steady-state glacier profile (as a reaction to the surface mass-balance conditions) is less influenced by the higher-order physics than it is by the surface mass-balance distribution. This transient glacier behaviour was investigated by applying surface mass-balance and ice-thickness perturbations in the accumulation area and to analyze the mass redistribution. Such a perturbation is transported towards the glacier front, where a significant thickening appears. The zeroth-order model transports this excess mass at a more or less constant speed towards the front, while the higher-order model transports this mass at a slower rate in the accumulation area, but faster in the ablation zone, so that the front position is more sensitive to migration than the zeroth-order model. The non-linear behaviour in mass transport is responsible for a 10% difference in time-lagged response to the climatic signal compared to the response of the zeroth-order model.
Acknowledgements
This paper forms a contribution to the Belgian Research Programme on the Antarctic (Federal Office for Scientific, Technical and Cultural Affairs), contracts A4/DD/E03 and EV/03/08A. The author is indebted to H. Blatter and D. Dahl-Jensen for their critical remarks as well as to the scientific editor, R. Greve.
Appendix A Coordinate Transformation
Due to the coordinate transformation (x, z) → (ξ, ζ), whereby ζ = (s − z)/H and ξ = x, the function derivatives transform to (Reference LliboutryLliboutry, 1987)
where ∂ξ/∂x = 1, ∂ξ/∂z = 0, ∂ζ/∂x = a, ∂ζ/∂z = −1/H. The function derivatives can thus be rewritten as
Where
so that the incompressibility condition (Equation (1)) becomes
Equation (14) then becomes after coordinate transformation, by making use of Equations (A3–A9),
Where
The strain are rewritten as
Equation (A11) can be decomposed further by replacing the terms involving ∂ω/∂ζ with the incompressibility condition of Equation (A10). Rearranging terms finally leads to
Similarly, the upper boundary condition (Equation (24)) is transformed to
where c = ∊ 2(∂s/∂x)2 − 1. The thermodynamic equation (Equation (16)) is transformed as follows:
Appendix B Finite Difference Approximation
Equations (A17-A19) are solved on an irregular grid in ξ and ζ. Central difference approximations on an irregular grid of first and second derivatives are governed by (e.g. Reference Payne and DongelmansPayne and Dongelmans, 1997)
The horizontal gradients ∂f/∂ξ and ∂2f/∂2ξ are obtained from Equation (B1) by replacing ζ k with ζi. An expression for ∂2f/∂ξ∂ζ is obtained by combining the finite-difference approximation of ∂f/∂ζ and ∂f/∂ζ At the surface boundaries i = Nξ and k = N ζ, first derivatives are computed using upstream differences:
At the surface boundaries i = 1 and k = 1, downstream differences are formulated by replacing (k − 1) → (k + 1) and (k − 2) → (k + 2) in Equation (B2).With the exception of the surface boundaries, central differences are used everywhere.
