1. Introduction
In the late 1950's and early 1960's, ideas on the response of valley glaciers to climate change were strongly influenced by the theoretical concept of kinematic waves (Reference Lighthill and WhithamLighthill and Whitham, 1955). They defined kinematic waves as a type of wave motion which exists in any continuum as a consequence of a conservation law and a coupling between discharge, concentration, and position. Application of this theory to glacier response was developed in a series of papers by (Reference WeertmanWeertman 1957, Reference Weertman1958) and (Reference NyeNye 1958, Reference Nye1960, Reference Nye1963). It was realised by these workers that the vertically integrated continuity equation for ice mass, together with a relation between ice discharge and (at least) ice thickness, supports the existence of kinematic waxes. With strong simplifications in the governing equations, some interesting results were obtained concerning the effect of mass-balance perturbations on the shape of a valley glacier. Notably, it was concluded that kinematic waves arc unstable in regions of decreasing ice velocity in the direction of flow (lower part of glacier) and stable in the region of increasing ice velocity (upper part of glacier).
We believe that some of the results have led to misinterpretation of features observed in the field. For instance, it has been claimed (Reference NyeNye, 1963) that the high level of trim-lines in the lower reaches of a glacier basin is in agreement with the sudden increase of ice thickness due to an unstable kinematic wave. This is only true because diffusion is neglected in the kinematic-wave solution. Including diffusion (which is realistic) yields a much smaller increase in the ablation area. A more logical explanation for the high level of trim-lines in the lower reaches of a glacier basin is therefore the slow retreat due to changes in the mass balance since the “Little Ice Age”. It is also argued (Reference NyeNye, 1960) that kinematic waxes are the fundamental reason why glaciers are such sensitive indicators of climatic change. This point has never been demonstrated in a convincing way from existing data. In our opinion glaciers are sensitive indicators of climatic change because a small change in climate may yield a large change in ablation, yielding an adjustment of the glacier volume (length). This can be demonstrated with simple continuity models.
On the other hand, the use of simple continuity-models as a tool to study the response of glaciers to climate change has been criticised (e.g. Reference PatersonPaterson, 1981, p. 268). It is a fact, however, that the simple continuity model has richer dynamics than the linearised-wave equation used by both Nye and Weertman. In view of this, we think that a comparison of linear-wave theory with the performance of a numerical “continuity model” is useful.
Before we start this comparison some comments on historical observations of kinematic waves are necessary. The classic example is the kinematic wave on the Mer de Glace from 1891 to 1899, as described by Reference LliboutryLliboutry (1965) and Reference Lliboutry and ReynaudLliboutry and Reynaud (1981). Figure 1 shows a kinematic wave with an amplitude of 2–6 m, a wavelength of about 3 km, and a wave velocity of 800 m a−1. The surface velocity itself increased from 125 to 155 m a−1, simultaneously over the entire ablation area. In spite of the limited accuracy of the survey technique at that time, this seems to be significant. Yet this is not what one would expect. Velocity should decrease after the passage of a kinematic wave at a specific point, since the ice thickness decreases locally and the surface slope decreases as well (ice thickness increases in the flow direction after the passage of the kinematic wave). Furthermore, an increase of 25% (from 125 to 155 m a−1) seems a bit too much for an increase in ice thickness of only a few metres. The simultaneous change in velocity presented in Reference Lliboutry and ReynaudLliboutry and Reynaud (1981) rather suggests die occurrence of a surge-type event with increased basal sliding, enabling velocity to increase over large areas of a glacier by about 25%. Although a few more examples of kinematic waves could be discussed, available data are limited. The major reasons for the scarcity of data on kinematic waves on glaciers are, in our opinion, the complex adjustment of a glacier to seasonal and long-term fluctuations in its mass balance, the rapid diffusion of kinematic waves and the lack of accurate field observations. All together these points make differentiation of the field observations between kinematic waves and variations in basal sliding very difficult.
Fig. 1. Change in the mean surface elevation of Mer de Glace, France, along four cross-profiles over a period of 9a. The broken line corresponds to a wave velocity of 800 m a−1 (from Reference LliboutryLliboutry (1958)).
In addition to a comparison with the linear theory, this paper attempts to present some quantitative insight into kinematic-wave velocities, and associated changes in ice thickness and ice velocity. In particular, we consider experiments using a “model glacier” with simple geometry and using the Hintereisferner, a valley glacier in Austria. Model experiments are presented which show the sensitivity of the kinematic waxes to variations in amplitude and wavelength.
We use a numerical ice-flow model in which the mass flux is directly coupled to the basal shear stress (e.g. Reference BindschadlerBindschadler, 1981; Reference KrussKruss, 1984; Reference OerlemansOerlemans. 1986; Reference Huybreehts, de Nooze, Decleir and OerlemansHuybreehts and others, 1989; Reference Stroeven, van de Wal, Oerlemans and OerlemansStroeven and others, 1989; Reference GreuellGreuell, 1992). We will not consider models that deal in one way or another with longitudinal stress gradients (e.g. Reference Budd and JenssenBudd and Jensen, 1975; Reference Shoemaker and MorlandShoemaker and Morland, 1984; Reference Van der Veen, Van der Veen and OerlemansVan der Veen, 1987), although these would be required if one wanted to simulate certain strongly localised dynamic features.
2. A Brief Review of Linear-Wave Theory
So-called continuity models are based on the vertically integrated mass-conservation equation together with a simple flow law for the vertical mean ice velocity. In the one-dimensional case:
Here H is ice thickness, t time, M specific balance, x distance along the flowline, U vertically averaged velocity, parallel to the bed, F “driving stress”. [C d, n, C s, m] a set of flow parameters, ρ ice density (910 kg m −3), g gravitational acceleration (9.8 ms−2) and h x surface slope. Subscripts d and s refer to contributions from internal deformation and sliding.
The kinematic-wave equation for glaciers, essentially based on linearisation of Equation (1) (Reference NyeNye, 1958: Reference WeertmanWeertman. 1958). reads:
C 0 is the kinematic-wave velocity, and D 0 the diffusivity of the kinematic waves. The linear theory assumes a reference state upon which small, independent perturbations of mass flux, ice thickness, and surface slope occur. The subscripts 0 and 1 refer to the reference and perturbed states, respectively. The reference state is normally, interpreted as an equilibrium state. Equation (4) can be solved if C 0 and D 0 are known functions of x. This is of course a valuable approach for gaining physical insight into the full non-linear problem. But Equation (4) also shows the limitation of this theory, since in reality Co and D 0 are not simple functions of x and in fact are only-known if the full ice dynamics arc included. The wave velocity as derived by Reference NyeNye (1958) i is given by:
Useful qualitative insight can be gained by prescribing C 0 and D 0 as a function of x. In Figure 2 the well-known results of Reference NyeNye (1960) are presented for the following formulation of C 0, with D 0 = 0:
Here r is a positive constant, and x is now the scaled glacier length (ranging from 0 at the glacier head to 1 at the terminus).
Fig. 2. (a) Time evolution of ice thickness in an idealised glacier following a sudden uniform increase in accumulation rate. The upper part of the glacier responds stably; the lower part responds unstably, until the kinematic wave from x = 0.5 arrives (from Reference NyeNye (1960)). The time of observation is scaled with the mean velocity gradient over the ablation area. The glacier length is also scaled (0 at the glacier head, 1 at the glacier front in the equilibrium state). (b) Time evolution of ice thickness following an addition of uniform layer of ice. The temporary instability of the lower half is relieved by the arrival of the kinematic wave generated at x = 0.5 (from Reference NyeNye (1960)).
The results presented in Figure 2, in which diffusion is neglected, will be compared with results from a simple How model that includes the effects of diffusion. The kinematic-wave velocity should be 3–6 times the horizontal surface velocity, because n and m are generally estimated to be 3–4 and 2, respectively. We use a model with n = 3 and no basal sliding. Therefore, a kinematic-wave velocity of 5 times the surface velocity is expected from the linear-wave theory.
3. Experiments with a Numerical Model
The numerical model solves Equations (1)–(3) on a grid. Flow is prescribed along a flowline with a constant width. Basal sliding is neglected (C s = 0). The flow parameter for deformation (C d) is set to 5 × 10−17 m6 Ν−1 a−1 for all experiments (including the Hintereisferner experiment). The specific balance is prescribed as:
where a is the balance gradient, h E the equilibrium-line altitude, and M max an upper limit (values of these parameters in all experiments reported here: 0.01 m m−1, 675 m and 1.25 m of water equivalent, respectively). A grid-point spacing of 0.1 km is used. Because of explicit time integration, this requires a time step of 0.05 a to maintain stability. It may be noted that all results presented are independent of the applied grid-point distance. This means that numerical diffusion is negligible. The steady-state properties of the resulting model glacier arc shown in Figure 3.
Fig. 3. Characteristics of the model glacier in equilibrium, (a) Bedrock and ice-surface elevation as a function of the scaled length, (b) Mass balance as a function of height E = 675 m, M* = 1.2 m. (c) Velocity as a function of the scaled length, (d) Gradient of the velocity as a function of the scaled length.
A: experiments with simple geometry (slope and width of bed constant)
Starting with the modelled glacier in a state of equilibrium, we have imposed a sudden and uniform instantaneous perturbation of the mass balance with a scaled amplitude H 1 equal to 1. The model calculates the time evolution of the ice thickness after the perturbation. The result shown in Figure 4a can be compared with Figure 2a because both are scaled by the velocity profile (r). The results obtained with the flow model are quantitatively different from the linear theory. There are three factors that can account for the observed discrepancies. First, the glacier length is kept constant in the linear theory, but this is not the case in the numerical model. Secondly, the linear theory neglects diffusion, which is not the case in the numerical model. Thirdly, in the linear theory ∂u/∂x equals r over the ablation area (and −r over the accumulation area) whereas in the model, ∂u/∂x decreases linearly over most of the ablation area (Fig. 3d). Overall the model yields a smaller increase in ice thickness, but the relatively greater increases in ablation area, particularly in the terminus area, are qualitatively comparable to those predicted from theory. The lower increase in ice thickness can be understood when we realise that the mass (flux in the model will increase as a result of the perturbation in the mass balance, yielding an adjustment of the modelled glacier length. A less pronounced increase in ice thickness can be observed in Figure 4a compared to Figure 2a as a result of the changing glacier length. However the most pronounced increase is found near the tongue of the glacier Since ∂u/∂x reaches a minimum value at the tongue, overruling the effect of diffusion in the upper ablation area. The role of diffusion is considered in more detail for a block of ice later on.)
Fig. 4. (a) Modelled increase in ice thickness due to an instantaneous increase in the mass balance. The increase in ice thickness is scaled with the magnitude of the disturbance. Length and time scaling is similar to Figure 2. (b) Modelled increase in ice thickness due to an instantaneous increase in the ice thickness. The scaling is similar to Figure 2.
A second comparison between the linear theory and the model experiments is presented in Figures 4b and 2b for a uniform increase in ice thickness (H 1) over the entire glacier (H 1 << H 0). The result depends only slightly on the magnitude of the perturbation. However, the quantitative result near the glacier front in Figure 4b should be regarded with some scepticism since it is affected by a relatively large truncation error. Nevertheless, a strong increase in ice thickness can be observed in the model results in the ablation zone. But the increase-in ice thickness is continuous in time (Fig. 4a and b), in contrast to the linear theory in which the increase in ice thickness shows a temporary instability. This instability is relieved by the arrival of the kinematic wave. The differences between theory and model results can be explained by the same arguments as stated in the previous paragraph. In a paper by Reference BindschadlerBindschadler (1982) kinematic waves were simulated with a similar flow model, using a different numerical approach and keeping the glacier length fixed. The resulting kinematic waves are similar to those in Figure 4a and b for central parts of the glacier.
Having drawn the comparison between the model and the linear theory, we will now give attention to quantitative aspects of the simulation of kinematic waves with a simple flow model in order to show the sensitivity of the results to the formulation of the imposed disturbances The dependence of the wave velocity on the wavelength of the disturbance (λ), the amplitude of the disturbance (A), and the mass-balance gradient will be presented. The same datum state as used in the previous experiments and presented in Figure 3a–d is considered. This equilibrium state is disturbed by an instantaneous perturbation of the ice thickness, a bump centred at the equilibrium line and described by
in which x
e is the grid point at the equilibrium line and x the grid point for which the perturbation is calculated. The amplitude (A) is in metres, and wavelength (λ) is dimensionless alter scaling by the glacier length (L). The subscripts 0 and 1 again indicate the equilibrium state and the perturbed state, respectively. Note that this formulation prescribes a positive perturbation of the ice thickness, since x ranges from 
, with a maximum perturbation at the equilibrium line. To study the influence of the wavelength and the amplitude we change only one variable (A, λ) at a time. The kinematic wave can be characterised by the velocity of the maximum perturbation of the ice thickness.
In Figure 5a the velocity of the kinematic wave is presented for different amplitudes (note that A « H 0) and a constant wavelength of half the glacier length. The velocity is scaled with the vertically averaged horizontal velocity at the equilibrium line in the equilibrium state. Larger wave velocities are observed for larger amplitudes which can be understood by the larger ice thickness and sleeper slopes (wavelength is constant), and hence greater ice velocities as well as greater wave velocities. Obviously the wave velocity increases towards the front of the glacier as |∂u/∂x| increases. This means that in the model glacier the expected reduction of the wave velocity due to diffusion is overruled by the increased wave velocity (wave velocity is proportional to the magnitude of |∂u/∂x|, Or (C 0) as a result of the increasing |∂u/∂x| towards the front of the glacier. Furthermore, we can observe than a doubling of the amplitude (A = 10 m) or a halving of the amplitude (A = 2.5 m) yields a roughly similar change in wave velocity. D 0 (= UH/hx ) can be evaluated from the steady-state conditions (Fig. 3). At the equilibrium line, U 0d ≈ 32 m a−1, H ≈ 164 m, and h, x ≈ 0.109, so D 0 is approximately 0.48 × 105 m2 a−1. The ablation area is about 2.8 km long, so ∂D 0/∂x averages about −17 m a−1 over the ablation area. For comparison C0 ≈ 160 m a−1 (Equation (5)), or an order of magnitude larger than the gradient in D 0. Using these estimates in the linear theory yields a kinematic wave speed, C 0 − ∂D 0/∂x (Equation (4)), of approximately 177 m a−1, or ∼5.5 times the surface speed. However, the scaled velocity calculated with the model (Fig. 5a) is typically 6–10 times the surface velocity over a large part of the ablation area, or somewhat larger than expected from the linear theory.
Fig. 5. (a) The scaled velocity of the kinematic wave for various amplitudes (A in metres) of the disturbance for λ = 0.5 L. The velocity is scaled by the velocity at the equilibrium line in the reference slate, (b) The scaled velocity of the kinematic wave for various wavelengths (λ scaled by the glacier length) of the disturbance (A = 5 m). (c) The scaled velocity of the kinematic wave for various mass-balance gradients (a) for a disturbance defined by A = 5 m and λ = 0.5. (d) As in (c), but the velocity is scaled by dividing by U∂u/∂x instead of U only.
A similar experiment is presented in Figure 5b. Here the kinematic-wave velocity is calculated for three different values of the wavelength, while the amplitude-is kept constant at 5 m. This experiment is somewhat more complicated lo understand because the larger increase in ice thickness for a longer wavelength, compared with a perturbation with a small wavelength, yields a higher diffusivity, bin on the other hand steeper slopes for smaller wavelengths yield a higher diffusivity. The net result is. however, an increasing diffusivity for longer wavelengths. If we consider the experiment in more detail, we may note that for a wavelength of 1λ, the tongue of the glacier is disturbed at the outset, as the front end of the perturbation is at the terminus (Equation (8)). The thickness at the tongue will, in this case, immediately begin to increase due to diffusion. Together with a transport of the wave in the flow direction, a greater velocity of movement of the locus of maximum disturbed ice thickness towards the tongue is expected. Following this line of reasoning, it is easy to understand that the wave velocities for λ = 1L are higher than for 
, although increasing wave velocities are still observed towards the margin. Doubling the wavelength from 0.5 to 1L yields a larger change in the wave velocity than halving the wavelength to 
(Fig. 5b).
The reason for the relative insensitivity of the wave velocity to perturbations in amplitude or wavelength is that the velocity gradient is comparable in the different model experiments presented so far.
By changing the mass-balance gradient, new equilibrium states can be calculated with a different ∂u/∂x profile in the flow direction. The new equilibrium states with different mass-balance gradients were perturbed with a wave with an amplitude of 5 m, and a wavelength of half the glacier. A smaller value of ∂u/∂x, due to a smaller mass-balance gradient (a = 0.005), reduces the kinematic-wave velocity as expected from linear-wave theory as can be observed in Figure 5c. A larger value of ∂u/∂x (a = 0.02) increases the kinematic-wave velocity. Scaling the velocity of the kinematic wave by dividing the kinematic-wave velocity by ∂u/∂x in the equilibrium state) yields equal wave velocities for the three different model glaciers, as shown in Figure 5d. This means that the kinematic-wave velocity scales with the ∂u/∂x This may seem trivial, but the difference with the linear theory is that the velocity is calculated in the numerical model and not prescribed.
To eliminate the influence of the gradient in the horizontal velocity, a few experiments were conducted for a block of ice with constant thickness and surface slope. This resulted in a steady stale with ∂u/∂x = 0 over the entire block length. The boundary condition at the outflow border, ∂H/∂x is constant, is time-independent, The kinematic-wave velocity is presented in Figure 6. The wave velocity is constant in the middle of the block because ∂u/∂x is negligible. Close to the margin, ∂u/∂x increases and the wave velocity increases accordingly. This result is in agreement with Bindsehadler (1982) who presented a constant wave velocity for kinematic waves in central parts of a uniform slab of ice with fixed length. One should, however, realise that diffusion cannot be neglected in general, since diffusion reduces the longitudinal velocity gradient and therefore the wave velocity if there is a velocity gradient. The wave velocity is independent of the wavelength and amplitude as long as λ ≤ 1 and A ≤ 10. the values used in the previous model glacier experiments.
Fig. 6. The velocity of a kinematic wave for λ = 0.5L and A = 5m in a block of ice (Η = constant) with ∂u/∂x = 0 initially. The velocity is scaled by the velocity in the reference state.
B: Hintereisferner experiment
In order to get some insight into how kinematic waves can he observed in a real glacier, a second set of experiments is presented. Kinematic waves were simulated on a glacier with the geometry and mass balance of Hintereisferner, Austria. Here variations in the width of the glacier and undulations in the bedrock may affect the results. A detailed description of this Hintereisferner model, and results from a Simulation of historical glacier variations were presented by Reference GreuellGreuell (1992). His model is nearly similar to the one used in the present paper because it also solves Equations (1)–(3). As a start, an equilibrium state is calculated, resembling the 1987 extent of the glacier (Fig. 7a–d).
Fig. 7. Steady-state characteristics of Hintereisferner as obtained from the numerical model, (a) Bedrock and ice surface elevation as a function of the distance from the glacier head. The glacier length is 7.4 km. (b) Mass-balance forcing as a function of height (Reference GreuellGreuell, 1992). (c) Vertically averaged velocity in the flow direction as a function of the distance from the bead, (d) As (c), for the gradient of the velocity.
Figure 8 shows the total ice volume after a perturbation of the ice thickness with a wave of 5 m, having a wavelength of half the glacier length and being centred around the equilibrium line. As soon as the wave front reaches the glacier terminus the volume decreases. This occurs after only a few years. Equilibrium is reached after about 70 a. Here the response time is defined as the time required to reach (1–1/e) of the volume change due to the disturbance in ice thickness. This response time is comparable to the one found by Reference GreuellGreuell (1992) for a perturbation in mass balance.
Fig. 8. Volume as a function of time after an instantaneous perturbation imposed on the glacier as it appeared in 1990 (A = 5 m, λ = 0.5L.).
The change in ice thickness at specific times after the onset of the perturbation is presented in Figure 9a. The maximum change moves down-glacier with a velocity 6 times the velocity at the equilibrium line. This is similar to the velocity of the kinematic waves in the model-glacier experiments discussed earlier. The geometrical boundary conditions create a rather stable zone about 5 km from the head of the glacier. Equilibrium is restored starting at the glacier head. A temporary increase of 200 m can be observed in the glacier length.
Fig. 9. (a) Increase in ice thickness at various times (t in years) after the onset of the perturbation (A = 5 m, λ = 0.5 L) as a function of the distance from the glacier head. The vertical line 11.) denotes the glacier length in the equilibrium state, (b) As in (a), now for a perturbation of 0.5 m in accumulation, lasting 1 a, over the entire glacier, (c) Time evolution of the ice thickness (as in Figure 9b) 20–100 a after onset of the perturbation.
In another experiment, an instantaneous uniform increase of 0.5 m in the specific balance during I a was imposed on the glacier. The insensitive zone is now centred around 4 km (Fig. 9b). Equilibrium is restored again from the glacier head down to the glacier front. The transition to a new equilibrium is shown in Figure 9b and c. No increase in glacier length is observed.
The effect of this perturbation on the velocity is shown in Figure 10. In the accumulation zone the velocity increases only slightly due to the counteracting effect of increased ice thickness and reduced surface slope (Fig. 10a). In the ablation zone the small changes in ice thickness and surface slope increase the local velocity typically by 10% (Fig. 10b). At the glacier terminus a very large increase is observed, but the large truncation error reduces the accuracy of the result here. The 10% increase in the velocity can be compared with observations at the Hintereisferner over the period 1916–21, as presented by Reference LliboutryLliboutry (1965). Changes in the surface velocity by a factor of 10, as Lliboutry presented, cannot be explained by a kinematic wave. The small changes in ice thickness and surface slope during the passage of the wave are insufficient to explain an order-of-magnitude increase in surface velocity Kinematic waves cause only a small increase in the local ice velocity, depending on position along the flowline.
Fig. 10. (a) The velocity scaled with the steady-state velocity at various times at two positions in the accumulation zone. The distance from the head of the glacier is denoted as x (in/km). (b) As in (a). for the ablation zone. Note the different scale along the vertical axis.
4. Conclusions
According to our numerical flow model, kinematic waves typically move with a velocity of 6–8 times the surface velocity, whereas linear theory predicts that they will move with a velocity of 5–5.5 times the surface velocity. This difference results from the increasing gradient in the horizontal velocity in the ablation zone in the flow model. This also explains the acceleration of the waves in the ablation zone. For a typical valley glacier the effect of diffusion is less important than the gradient in the horizontal velocity, as demonstrated by Figures 7 and 8.
Observing kinematic waves in the field is difficult due to the rather small changes (10%) in the local ice thickness and velocity (Figs 9 and 10). Furthermore, time series of both parameters must be observed, since ice thickness and velocity are coupled. Observations of ice thickness or ice velocity alone cannot discriminate between processes related to deformation and changes in sliding. Synchronous changes in ice velocity over a glacier are not an indication of kinematic waves.
No length variations can be observed (Figs 9b and c) for an instantaneous increase (during 1 a) of 0.5 m in the accumulation over the entire glacier. This suggests that observing only front variations is probably insufficient to detect kinematic waves in the field.
We believe that observations of so-called kinematic waves are often associated with variations in basal sliding (and not related directly to perturbations in ice thickness).
The numerical experiments discussed here have show that flowline models with local coupling of velocity and thickness/slope simulate kinematic waves well. These models are therefore suitable for the simulation of glacier fluctuations on a lime-scale of more than a few years. We believe that the success of such simulations is determined largely by the accuracy with which the mass-balance history can be reconstructed or formulated.
Acknowledgements
We thank W. Greuell lor providing us with the source code of his Hintereisferner model. We appreciated the useful comments on this work made by A. Stroeven and the discussions in the ice-and-climate group at our Institute. We also thank the reviewers for their thorough and helpful comments. Financial support was provided by the Netherlands Organization for Scientific Research (NWO).
