A numerical study of glacier advance over deforming till

. The advance of a glacier over a deforming sediment layer is analysed numerically. We treat this problem as a contact problem involving two slowly-deforming viscous bodies. The surface evolution of the two bodies, and of the contact interface between them, is followed through time. Using various different non-linear till rheologies, we show how the mode of advance depends on the relative effective viscosities of ice and till. Three modes of advances are observed: (1) overriding, where the glacier advances through ice deformation only and without deforming the sediment; (2) plug-ﬂow, where the sediment is strongly deformed, the ice moves forward as a block and a bulge is built in front of the glacier; and (3) mixed-ﬂow, where the glacier advances through both ice and sediment deformation. For the cases of both overriding and mixed-ﬂow, an inverse depth-age relationship within the ice is obtained. A series of model experiments show the contrast in effective viscosity between ice and till to be the single most important model parameter deﬁning the mode of advance and the resulting thickness distribution of the till. Our model experiments indicate that the thickness of the deforming till layer is greatest close to the glacier front. Measurements of till thickness taken in such locations may not be representative of deforming till thickness elsewhere. Given sufﬁciently large contrast in effective viscosity between ice and till, a sediment bulge is formed in front of the glacier. During glacier advance, the bulge quickly reaches a steady state form strongly resembling single-crested push moraines. Inspection of particle paths within the sediment bulge, shows that particles within the till travel at a different speed from the bulge itself, and the push moraine to advance as a form-conserving non-linear wave.


Introduction
Glacier length changes have been recorded by measurements of the snout position and their mass-balance (Oerlemans, 1989;Haeberli and Beniston, 1998;Oerlemans, 2001), by geomorphological observations of the glacier forefield deformation, as well as by using traces in the landscape and historical glacier records to interpret the glacier fluctuations in the past (e.g., Drewry, 1986;Benn and Evans, 2010;Haeberli and Beniston, 1998;Oerlemans, 2001).Observed length changes have been used to extract climate signals (e.g., Boulton, 1986;Oerlemans, 2005;Beedle et al., 2009).
Glacier flow can occur by several mechanisms: (1) ice deformation; (2) deformation of the glacier bed; (3) sliding at the ice-bed interface; or (4) a combination of these processes (e.g., Benn and Evans, 2010).A glacier advance by ice deformation only, i.e.where the ice is frozen to the bed, leads to an advance by "overriding", which means that the advancing glacier snout rolls its surface over the glacier forefield, potentially giving rise to inverse depth-age relationship within the ice behind the glacier front.For an advance by (2) or (3) only, the ice mass is moving as a block ("plug-flow").A glacier advance over a deformable till can lead to the formation of "push moraines", where the glacier appears to be bulldozing over and through the ice-marginal sediments (Boulton, 1986;Van der Wateren, 1995;Benn and Evans, 2010;Bennett, 2001;Motyka and Echelmeyer, 2003).
In contrast to the large number of observational studies of the advance of glaciers over a deforming substrate, there is a conspicuous lack of numerical studies.Although a variety of numerical models have been used for calculating temporal changes in front position (e.g., Greuell, 1992;Schlosser, 1997;Schmeits and Oerlemans, 1997;Wallinga and Van de Wal, 1998;Oerlemans, 2001;Smedt and Pattyn, 2003), none of these models attempt to describe the details of ice flow, or G. J.-M. C. Leysinger Vieli and G. H. Gudmundsson: Glacier advance over deforming till the mechanical interactions between ice and till in the vicinity of the glacier terminus.Here, assumptions commonly used in flow modelling of large ice masses, such as the shallow ice approximation (SIA) model (Hutter, 1983) do not apply.Furthermore, we are not aware of any numerical work where the time-dependent deformation of both the ice and the till is calculated in a fully coupled way using the full-set of the momentum equations.
In this paper we aim to answer the question of how the front of a glacier residing on a layer of sediment advances and to investigate the details of the flow at an advancing glacier front, by using a two dimensional finite-element model designed for this purpose.This study extends and complements a previous study of ours (Leysinger Vieli and Gudmundsson, 2004) focusing on the advance of a glacier over non-deforming substrate.
Our modelling approach is generic in its nature and we do not try to simulate the advance of any one particular glacier.Rather, we are interested in the overall characteristics of the flow of ice and till, and the mechanical interaction between the advancing glacier and the underlying till.Of special interest is the flow regime in the vicinity of the snout where the assumptions of most previous models break down.We aim at answering how till deformation is induced by the flow of a glacier, how the thickness distribution of till changes as a glacier advances over it, and where till deforms primarily through shearing as opposed to horizontal compression.We study numerically the formation of push moraines and suggest that push moraines may be considered to be a type of form-conserving non-linear till waves.

Modelling the advance of a glacier
For isotropic ice Glen's flow law is commonly used as a constitutive law (Glen, 1955).For subglacial till various different rheological models have been proposed (e.g., Fowler, 2002;Clarke, 2005).The rheology of the till has been suggested to behave like a linear or slightly non-linear viscous fluid (Alley et al., 1986(Alley et al., , 1987;;Boulton and Hindmarsh, 1987;Alley, 1989;Alley et al., 1989;MacAyeal, 1989MacAyeal, , 1992)), where the strain rates ˙ are related to the effective stress τ through with a stress exponent m in the range of 1≤m≤5 (e.g., Kamb, 2001), such as m = 1.33 as fitted by Boulton and Hindmarsh (1987) for the till beneath Breidamerkurjökull or m = 3.04 as optimised by Gudmundsson (2007) for the subglacial till of Rutford Ice Stream.It has also been pointed out that since till is a granular medium its rheological behaviour should fall within the range commonly exhibited by such materials (Iverson et al., 1998;Tulaczyk et al., 2000a,b;Truffer et al., 2000;Kamb, 2001;Kavanaugh and Clarke, 2006).Therefore, it should show Coulomb-plastic rheology (termed "treiboplastic" by Kamb, 2001), with a yield stress which is controlled by intergranular friction and which depends only slightly, if at all, on the strain rate (Kamb, 1991;Iverson et al., 1998).
A perfectly plastic till rheology can be obtained by letting the flow exponent m in Eq. ( 1) go to infinity.For stresses below the yield stress (τ o ), deformation is negligible and for stresses exceeding the yield stress the deformation would be infinitely large.As a result, the stress in the sediment never exceeds the yield stress (Paterson;1994;Van der Veen;1999).Kamb (2001) suggests that till rheology, as observed for ice streams B, C and D (now Whillans, Kamb and Bindschadler), is best described as "imperfect treiboplastic", because of an appreciable strain rate (˙ ) dependence on the stress (σ ), which he expresses in a highly non-linear flow law with an exponent m ≈ 40 ± 20.Here we use the term near perfect plastic to refer to such a highly non-linear rheology.
As pointed out by Fowler (2002Fowler ( , 2003) ) a rigid/perfectly plastic material can deform over long time-scales as if viscous.More recently the controversy is dealt with by investigating constitutive models of subglacial sediments which allow for different deformation under different conditions (e.g., Sane et al., 2008;Altuhafi et al., 2009).Here we consider cases corresponding to both non-linear viscous and near perfect plastic rheology.

Model description
A time-dependent two-dimensional numerical plane-strain model is used to study the advance of a glacier over a layer of deforming till.The model treats the glacier and the underlying till as two independent slowly-deforming viscous bodies.The surface evolutions of both bodies are followed with time and contacts between the bodies are automatically determined.In each time step, nodes that enter into or out of contact are found and the boundary conditions are changed accordingly.Unstructured gridding and automated remeshing are employed to limit element distortion.The full set of the equilibrium equations is solved, allowing for accurate estimates of stress and strain.

Ice rheology
Glen's flow law (Glen, 1955;Steineman, 1958) is used to describe the relationship between deviatoric stresses and strain rates in the ice as well as in the sediment layer.Here ˙ ij are the components of the strain rate tensor, and σ (d) ij are the components of the deviatoric stress tensor The Cryosphere, 4, 359-372, 2010 www.the-cryosphere.net/4/359/2010/where σ ij are the components of the stress tensor and δ ij is the Kronecker delta.The effective stress τ is defined by The rate factor A and the stress exponent n are material parameters.For the ice body the flow law exponent is chosen as n = 3.Values for other parameters are given below.All parameters are made non-dimensional by the use of appropriate scalings (see Sect. 2.1.4and Appendix A).

Till rheology
The till is treated in a formally similar manner as the ice, that is as an incompressible non-linear viscous medium (Glen's flow law, Eq. 2).The rate factor and the stress exponent of the till are denoted by A and m, respectively.The till density was set to ρ s = 2200 kg m −3 .It is one of the interesting facts of glaciological research that despite decades of work on this subject no general consensus has emerged on the most appropriate numerical values for the modelling parameter m.Values for m ranging from 1 to infinity are commonly used (e.g., Iverson et al., 1998;Kamb, 2001;Clarke, 2005).Here we bracket the range of possible rheologies by using a wide range of values for the rate factor (A ) and the stress exponent (m).
The till rheology is determined by a number of factors and processes that do not concern us here.We do not consider in detail how the particular numerical values of the rheological parameters of the till are related to grain size, pore water pressure, etc. which is an important question in itself but entirely unrelated to our modelling work.

Model geometry and mass balance distribution
The general model geometry is depicted schematically in Fig. 1.The glacier is everywhere in contact with a deformable till layer.The till layer extends sufficiently far to the front of the initial glacier geometry to ensure that while advancing the glacier remains in contact with the till.Initially, the thickness of the till is uniform and set to about 1/20th of the maximum ice thickness.Where the ice is in contact with the till no differential motion between ice and till is allowed.No mass is added to, or extracted from, the sediments in the course of the calculations.In all our experiments the geometry of the base of the till layer is kept fixed, with a steep upper part (45 • ) and a lower part sloping at 5 • as shown in Fig. 1.Our choice for the initial geometry and the slope of the base of the till layer does not affect our final conclusions.
We chose a height-dependent mass balance function such as where ELA is the equilibrium line altitude, located at a height of 0.5 dimensionless units (defined below) and z s is the surface height of the glacier.Above the ELA a positive mass balance gradient a acc of 0.0117 (non-dimensional units) is used.Note that this mass-balance distribution includes altitude-mass balance feedback in the accumulation area only.The use of this particular form of the surface massbalance distribution is motivated by the requirement that the overall mass balance must remain positive in order to ensure continued advance.Setting the mass-balance to zero for z s ≤ELA simplifies our analysis of the movements of material particles in the vicinity of the snout.Due to this chosen mass-balance distribution no equilibrium is reached and the glacier is in a perpetual state of advance.An initial glacier geometry was generated for a hard bed by running a SIA model.Once the glacier had advanced to a total length of three times the size of the accumulation area, we switched from using the SIA model and did the rest of the run using the full-system model.Our conclusions are not affected by the details of the starting geometry, as the geometry obtained from the SIA model is, except for the frontal region, similar to the one calculated for the FS model (see Leysinger Vieli and Gudmundsson, 2004) and the differences at the front are quickly adjusted by the different flowfield obtained by the FS model (see e.g.adjustment for the initial surface marked with "b" in Fig. 2a).The use of SIA in the initial modelling stage is based on practical considerations and a desire to limit computational times involved.

Normalisation
We introduce non-dimensional variables for the space coordinates {X,Z}, the components of the velocity vector {U,W }, the time {T } and the effective stress

Numerical solution technique
The model calculations were performed with the commercial finite-element (FE) program MARC (2000) tailored for use in glaciology by Gudmundsson (1999).The code solves the full set of momentum equations and has been used extensively for flow modelling of alpine and grounded tidewater glaciers (Gudmundsson, 1999;Vieli et al., 2001;Leysinger Vieli and Gudmundsson, 2004).We refer to this model as the full system model (FS).The transient evolution of the surface and the ice-sediment interface is followed using a mixed Lagrange-Euler approach.As the surface evolves, automated remeshing is used when needed to limit mesh distortion (Leysinger and Gudmundsson, 2000;Leysinger Vieli, 2003).The remeshing strategy (Leysinger Vieli, 2003) uses a cubic spline interpolation to reposition nodal points along all singular surfaces.It has been verified that the remeshing algorithm preserves the total volume accurately and that the model results are not affected by smaller time steps or more refined mesh.The FE code has been verified by comparing its numerical results with transient analytical solutions (Gudmundsson, 1997;Adalgeirsdóttir et al., 2000;Raymond et al., 2003;Gudmundsson, 2007) and by comparing the time-dependent response for a linear flow law (Newtonian flow) with the analytical solution (ISMIP-HOM Benchmark experiment F, Pattyn et al., 2008).

Model parameters
We consider two end-member cases of different modes of advance where (a) the till is much stiffer than the ice and all the deformation takes place in the ice, and (b) the opposite situation where ice is much stiffer than the till and all deformation takes place within the till layer.These different limiting cases of modes of advance can be generated in our numerical model simply by varying the contrast in effective viscosity between ice and till.After normalisation, the results of the numerical runs depend only on the ratio of the rheological parameters.We keep the rate factor A and the stress exponent n fixed and vary the corresponding rheological parameters A and m of the till.The non-dimensional rate factor for the ice is A = 2 (Eq.A4), this follows directly from the scalings used, when we set n = 3.For a soft, highly deforming sediment the ratio between the rate factor A in the ice and A in the sediment is small ( 1) and for a highly rigid sediment this ratio is large ( 1).For the till, values of both m = 3 and m = 40 were used.The value m = 40 is used to simulate a nearly perfect plastic rheology as suggested by Kamb (2001).
To be able to get similar effective viscosities (η = 1/2A τ m−1 ) in the sediment at given stresses with different m we define B for the till as and rewrite Glen's flow law (Eq.2) as Note that the material parameter B has the dimensions of stress, and that increasing B has the effect of decreasing the strain rates for a given state of stress.We refer to B as the hardness factor while A is the rate factor.For m→∞, Eq. ( 7) shows that the rheology of till corresponds to a perfectly plastic rheology with B being the yield stress τ o .Hence, by adopting a high stress exponent such as m = 40, the sediment rheology goes from a non-linear viscous to a near perfect plastic one.For such a high stress exponent the till will only deform significantly in regions where the effective stress τ is approaching B. The values used in the model experiments are shown in Table 1.

Results
As already mentioned, the model is formulated in dimensionless variables and therefore the results shown are all dimensionless.One can easily go from non-dimensional variables to dimensional ones.Our non-dimensional vertical space coordinate Z extends from 0 to 1, as the extent of the grid in z has been chosen to be the same as the thickness [h] of the infinitely extended inclined ice slab (see Appendix A).For example, choosing a rate factor of A = 2.06×10 −15 s −1 (kPa) −3 for temperate ice (Hubbard et al., 1998;Gudmundsson, 1999;Albrecht et al., 2000), a mean thickness [h] = 1000 m for the ice slab and a slope of α = 5 • results in a mean horizontal surface velocity of For this example one non-dimensionless unit in length, velocity and in time, corresponds to 1000 m, 1.5665 × 10 4 m a −1 and 0.0638 years, respectively.Note that since our nondimensional vertical space coordinate Z extends from 0 to 1, the mean dimensionless thickness of the modelled glacier [H ] is approximately 0.2 Z (e.g.see Figs. 3,10,13 and 14), which corresponds to a mean glacier thickness [h o ] = 200 m for our example.

Experiment A: relatively stiff and moderately non-linear till
We start our discussion of the numerical results by considering a case where the till is relatively stiff in comparison to the ice.We use B = 8, which is a high value as compared to the mean effective basal stress of unity, and m = 3.This experiment illustrates how glaciers advance over stiff till, and serves as a reference case for the set of experiments described later where the advance gives rise to significant till deformation.
The front advance is shown in Fig. 2a as a series of snapshots depicting the longitudinal surface profiles at different times.The surface of the till is also plotted in Fig. 2a for the same times, but, as expected, no change in till thickness is visible due to its relative stiffness in this experiment.
The surface profile at t = 0, marked with "b", in Fig. 2a corresponds to the last profile obtained with the SIA model, and the starting profile of the full-system (FS) model.As Fig. 2a shows, the surface profile of the SIA model is markedly different from the one obtained with the FS model.The SIA model calculates velocity from local slope and thickness, and predicts infinitely large velocities for vertical slopes.The FS model includes the effects of horizontal deviatoric stresses on the flow and gives a much more accurate description of the flow field in the vicinity of glacier terminus.Figure 2a shows how the surface slope initially increases, finally resulting in an overhanging snout.Figure 2b-d show how various material particles (dots in the figures) within the ice move as the front advances under no-melt conditions.In Fig. 2b positions of selected material particles at the surface, at the base, and within the glacier are shown at different points in time.Due to the relative stiffness of the till, all material particles at the ice-sediment interface stay fixed at their respective locations, and the ice-sediment interface remains parallel to the bed (see Fig. 2b-d).Material particles within the ice, however, gradually move down-slope eventually passing the position of the initial glacier front.As the front advances, and owing to our assumption of nomelt, ice particles that previously were located at the glacier surface come in contact with the till surface.This leads to overfolding within the ice of the glacier front near the icesediment interface, and forms a new basal layer of ice giving rise to depth-age inversion.
Figure 3 shows vertical profiles of the non-dimensional horizontal velocity U and effective stress T eff at three different positions along the glacier for the overriding case.At each location the effective stress generally increases with depth.Note that at the snout the effective stress reaches its maximum value at the interface between the ice and the till, and not at the base of the till as is the case further upstream.

Experiment B: relatively soft and moderately non-linear till
We now consider the situation where the till is considerably softer than the ice by setting B = 0.35.In this experiment the glacier advances primarily through deformation of the underlying sediment and without any significant deformation of the ice.In Fig. 4a, glacier and sediment surface profiles are shown at different times.Clearly visible is the formation of a sediment bulge in front of the glacier.The amplitude of the bulge as a function of distance is shown in Fig. 4b.The geometry and the positions of a number of material points within the ice and at the till interface are shown in Fig. 4c-e and the corresponding profiles in Fig. 4a are labelled c, d, and e.Further details of the sediment deformation are given in Figs. 5 and  6.
One of the most interesting results of the experiment is the formation of the sediment bulge in front of the glacier snout.After an initial period of rapid growth the amplitude of the bulge approaches a steady state value (see Fig. 4a).Following the movement of material points within the till (see Figs. 5 and 6) reveals that till particles within the bulge travel at a different speed to that of the bulge itself.Thus, the bulge is an example of a propagating till wave.Till particles get incorporated into the till wave for a limited period of time, and are then left behind as the wave travels further.The distribution of horizontal and vertical velocity, and of the effective stress is shown in Fig. 7a, b, and c, respectively.This figure can be compared directly to Fig. 8a, b, and c depicting results from Experiment A. Because of the relative stiffness of the ice, the velocity field of the glacier is close to being constant everywhere, and almost the entire forward motion is due to till deformation (Fig. 7a).The vertical velocity distribution within the till bulge (Fig. 7b) shows how the part of the bulge in contact with the ice is being pressed downward, resulting in negative vertical velocities, while the down-glacier side of the bulge moves upwards.The topmost layer shows the strongest deformation due to the uplifting of sediment needed to form the bulge (b-d).After a while the bulge is stable in amplitude and shape, and only a small uplift of the sediment takes place at the down-facing side of the bulge (c).Three former vertical straight lines (solid) illustrate how such lines deform with time.the ice and the till.This is the required stress condition for the stress fields of two bodies in contact.The continuity of the stress field seen in Figs.7c and 8c is a demonstration of the accuracy of the numerical solution.
The velocity distribution within the sediment bulge is shown in detail with velocity vectors in Fig. 9.As the figure shows, the ice in the frontal area presses the till downwards and forwards at the same time.Note that the ice velocity is larger than that of any of the material particles in the till bulge.With time the glacier therefore overrides the till particles previously within the bulge, while at the same time, due to the upward motion in the frontal region, new till particles are incorporated into the bulge.
From the ratio of basal to surface velocity in Fig. 10b we see that at the glacier front the flow is dominated by sediment deformation (U b /U s ≈1), and that the contribution of the ice to the total deformation increases in the upstream direction.The increase in till deformation towards the glacier front is also seen in the vertical profiles of the horizontal velocity in Fig. 10c-e and in Table 2.The observed distribution of the effective stress in the three vertical profiles is similar to the one obtained for Experiment A with the exception of the frontal profile, where the highest shear stress is found at the sediment base (Fig. 10e) and not at the ice-sediment interface.

Experiment C: moderately stiff and moderately non-linear till
We set B = 2 −1/3 ≈0.79, giving an equal rate factor within the ice and the till.The advance now follows through a mixture of both ice and till deformation (see Figs. 11a-e and 12).In this experiment, flow features that were only seen either in Experiment A or B act in combination.As the ice advances over the till, material particles previously at the surface come in contact with the till and form a new basal ice layer.At the same time the till deforms giving rise to a propagating till wave similar to that seen in Experiment B (Fig. 11a-e).
Detailed view of the till deformation is given Fig. 12.The advance of the glacier causes shearing within the till, combined with a horizontal compression in the snout region.As in Experiment B a sediment bulge is formed, but due to the larger stiffness of the till the bulge is about 8 times smaller.
The ratio of basal to surface velocity as a function of distance is shown in Fig. 13b.The ratio increases sharply towards the snout implying that measurements of till deformation and surface velocity in that area are not reliable estimates of the glacier wide contribution of basal motion to the forward motion of a glacier.The effective stress distribution shown in Fig. 13c-e also reveals that the magnitude and distribution of stresses within the till in the terminus area (Fig. 13e) differ strongly from those found further up-glacier.In the terminus area, the stresses in the till decrease with depth rather than increasing as is the case elsewhere.

Non-linear viscous vs. near perfect plastic rheology
In the above experiments glacier advance over deforming till was modelled using a stress exponent m = 3 (Eq.7).Here we repeat those experiments for a near perfect plastic till rheology by choosing a flow law exponent m = 40 (Kamb, 2001).We performed a number of experiments using parameter sets that produced both relatively stiff and moderately soft till rheologies.
The main effect of increasing the value of m is to sharpen the contrast between those regions where the effective stress, τ , in the till is either similar to or smaller than the hardness factor B (see Eq. 7).Where τ >B the till deforms readily, but hardly at all where τ <B.We found in all our experiments that the highest effective stresses in the till were reached in the terminus area.One of the consequences of increasing the value of m from 3 to 40, while keeping B fixed, was therefore to raise the ratio between basal motion and surface velocity in those areas where τ >B (see Fig. 14).For example, for B = 0.42 (Fig. 14c-e) at approx.1/20[H ] behind the glacier front 62% of the horizontal deformation is due to till deformation, 33% at 2/3[H ] and no till deformation is taking place at only 5 times the mean ice thickness, [H ], behind the front (Table 2).
By increasing the hardness factor B from 0.42 to 0.44, the value used in Experiment B described above, the point of transition from non-deforming to deforming till migrates further towards the snout.The contribution of sediment deformation for B = 0.44 to the total horizontal deformation decreased to 56% at 1/20[H ], 23% at 2/3[H ] and 0% at 5[H ] behind the glacier front (Table 2).
The till thickness generally increases towards the terminus.This can be seen clearly in Fig. 14c-e where the till thickness is indicated by a horizontal dotted line (note  different vertical scales).Note that since the initial prescribed till thickness was everywhere uniform, this thickening of the till in the terminus area is solely caused by the flow of the glacier and the mechanical interaction between the ice and the till.Note also that since we do not describe in the model any sediment movements due to action of water and/or erosion, all sediment deformation is solely due to shearing through the till column and horizontal compression towards the terminus.The formation of a sediment bulge in front of the snout was also observed for m = 40 just as it was for m = 3.Although the details of the form and the shape of the bulge change as modelling parameters such as m and B are varied, the genesis of this bulge, hence, does not depend on the exact parameter values chosen and is, in this sense, a robust modelling feature.

Discussion
Our model was not designed at the outset to investigate the formation of push moraines.As illustrated by the findings of a study of the large-scale push moraine of Taku Glacier, Alaska, by Kuriger et al. (2006), the shape and form of the  sediment bulge forming in front of the glacier terminus, however, bears strong similarities to those of push moraines.Furthermore, we find the value of the till hardness factor B to be the single most important modelling parameter affecting the mode of till deformation.This agrees favourably with Kuriger et al.'s (2006) finding that decreased effective frictional resistance of the sediment is the key parameter affecting rates of internal sediment deformation.Numerous descriptions of the morphology of push/ squeeze moraines and conceptual models of their formation are found in the literature (e.g., Price, 1970;Kälin, 1971;Meyer, 1983;Van der Wateren, 1995;Winkler and Nesje, 1999;Bennett, 2001;Evans and Hiemstra, 2005;Kuriger et al., 2006;Evans, 2009).To our knowledge the model presented here is, however, the first numerical model of glacier advance over till that produces a feature closely resembling push moraines as they are observed.
The process of advance producing a push moraine is often observed to be a seasonal process with advance in late winter while the till is frozen and squeezing out of water-soaked till from beneath the glacier during the early summer melt period (e.g., Price, 1970;Boulton, 1986;Evans and Hiem-   stra, 2005).By varying the hardness of the sediment seasonally our model should be able to reproduce such seasonal push/squeeze moraines.Truffer et al. (2009) inferred from field measurements and numerical modelling that the mode of glacier advance changes on seasonal timescales on Taku Glacier, Alaska, from internal deformation (winter) to plugflow (summer), depending only on small changes in the basal boundary condition.
The numerical modelling experiments show how a glacier advancing over a soft till causes shearing within the till that, depending on the value of the stress exponent m, gives rise to spatially highly variable rates of basal motion (see Figs. 13b  and 14b).In all our modelling experiments, the highest ratio of basal motion to surface velocity was found in the terminus area.Measurements of basal motion in that area are therefore potentially not representative of a glacier as a whole.An example for such a measurement site is the tunnel system dug in the ice about 1 to 2 m above the glacier sole and approximately 24 to 44 m behind the terminus of Breidamerkurjökull, Iceland, described by Boulton (1979).By assuming no slip at the glacier sole, the displacement due to sediment deformation at this site was between 95% and 80% of the forward movement of the glacier (Boulton and Hindmarsh, 1987).If we assume an average thickness of 400 m for Breidamerkurjökull (Björnsson et al., 2001; using slope as in our model), then the position behind the terminus is about 1/16 to 1/9 of the typical mean ice thickness [h o ] behind the glacier front.Their measured sediment displacement, especially at the higher end, is in the same range as obtained in this study for the plug-flow case (Experiment B).Mechanical interaction between ice and till alone is sufficient to give rise to till thickening with time in the terminus area.The thickening is caused by combined action of shearing through the till column and horizontal compression towards the terminus.No other potential modes of sediment transport, i.e. suspended sediment transport or entrainment of sediment into the ice, were described in the model.Including erosion and/or hydraulic transport would modify longterm evolution of the sediment thickness (e.g., Motyka et al., 2006).Similarly, the moraine structure in front of the terminus would be modified by pro-glacial fluvial and eolian processes.
If sufficiently soft, the till is extruded from underneath the glacier towards the terminus.This happens through direct vertical mechanical action of the ice on the underlying sediments.Such moraine forming by extrusion flow has been observed in nature and discussed in the literature (e.g., Price, 1970;Sharp, 1984;Evans and Hiemstra, 2005).The shape and the position of the bulge as depicted for example in Figs. 5 and 6 resembles push/squeeze moraines (e.g., Price, 1970;Evans and Hiemstra, 2005;Evans, 2009).What has been touched on by Van der Wateren (1995) but possibly not been fully realised since, and is also suggested by our modelling, is the possibility that the push moraines are a form of non-linear shape-conserving waves where horizontal movement of individual till particles is insignificant in comparison to the total distance travelled by the moraine.We also emphasise that the formation of the push moraine does not require the terminus to be in a steady state for any significant period of time.In fact in our modelling approach the terminus is never in steady state but advances continuously.
Where the till is too stiff to be significantly deformed, the glacier advances through "overriding" or "mixed-flow".In nature examples of this type of advance can be found in advancing glaciers in sub-/polar regions where melting is negligible (e.g.Glaciers on Axel Heiberg Island, Canadian Arctic Archipelago and in the Dry Valleys, Antarctica).A possible example of a glacier advance by overriding or by mixedflow is Crusoe Glacier from Axel Heiberg Island where the folding structure revealed in a photo taken of the west front by Alean (2008)  cier advance calculated with a altitude-dependent mass balance distribution over the whole ed in Leysinger Vieli and Gudmundsson (2004) and shown using the same scale as in Fig. 2a.41 Fig. 15.Glacier advance calculated with a altitude-dependent mass balance distribution over the whole glacier as used in Leysinger Vieli and Gudmundsson (2004) and shown using the same scale as in Fig. 2a.example of advance by overriding are active rock glaciers, which have insignificant mass loss through surface ablation.On the other hand, for glaciers with surface melting in alpine environments a realistic shape of an advancing glacier front would be less steep.Previous model calculations for the evolution of advancing glaciers residing on a bedrock with a realistic mass balance distribution, done by Leysinger Vieli and Gudmundsson (2004), indeed show that the front geometry is less curved than the one with zero mass balance (Fig. 15).Although incomplete, one would still expect to find the overfolded ice layers within the glacier front for the overriding and the mixed-flow case, as widely seen in terminal ice cliffs (e.g.Hook and Hudleston, 1978;Benn and Evans, 2010).Due to ablation, the surface layers are continously truncated, and therefore the completeness of the folding structure depends very much on the prevailing melt and/or dry calving rate.
Raising the flow exponent m from 3 to 40 changes the vertical distribution of the horizontal velocities within the sediment.For this near perfect plastic rheology the velocity increase is largest close to the base of the sediment layer (Fig. 14d), whereas for m = 3 the increase is approximately linear over depth.The main effect of increasing the value of the stress factor, however, is to increase the relative contrast in the contribution of basal motion to the overall forward motion between different areas.The result is a more pronounced spatial variation in distribution of basal velocities.

Conclusions
By modelling numerically the advance of a glacier over a deforming sediment layer as a contact problem involving two viscous bodies, we have shown how push moraines can be formed.During glacier advance, till particles enter and then subsequently leave the push moraine.After an initial phase the size and the shape of the push moraine does not change with time.The formation of the push moraine is a "robust" modelling result in the sense that its genesis does not depend on a fortuitous choice of parameters but occurs in the model for a wide range of parameter values.
Although we only modelled a glacier in its advancing stage, the formation of the moraine and analysis of the trajectories of the till particles suggest that such a push moraine would also form during a stationary stage or even during a retreat phase.
Despite the model parameters describing the till rheology (m and B), as well as the initial till thickness distribution, to be spatially uniform, resulting distribution of basal motion is highly non-uniform.Highest rates of basal motion, as compared to surface velocity, are found in the vicinity of the terminus area.Rates of basal motion measured close to a glacier terminus are not reliable estimates of basal motion a few ice thicknesses further up-glacier.

Fig. 1 . 27 Fig. 1 .
Fig. 1.Schematic view of the glacier considered in the two-layer model.The ice body on a sediment layer (grey) which extends into the glacier forefield.The inclination of th 45 • and of 5 • for the lower part.figure

Fig. 2 .Fig. 2 .
Fig. 2. Experiment A: Relatively stiff and moderately non-linear till (m=3, B=8): (a) are shown at time intervals of approximately two dimensionless time units.Selected fro surface shapes marked (b), (c) and (d) are shown in more detail in subfigures (b-d).(b the front showing the surface, four selected internal layers and the ice-sediment interface times T equal to 0 (b), 10 (c) and 20 (d).The dots indicate material particles that throughout the three pictures.The connecting lines illustrate how material layers (fine lines) with initial shape as in (b) become deformed with time.The arrow marks the initi glacier front.

Fig. 3 .Fig. 3 .
Fig. 3. Experiment A: Relatively stiff and moderately non-linear till (m=3, B=8): (b-d) of the non-dimensional horizontal velocity U (lines with crosses) and effective stress T eff are shown at a dimensionless time T =15.4 for the three locations A, B and C shown in ( vertical scale changes for each profile.

eFig. 4 .
Fig. 4. Experiment B: Relatively soft and moderately non-linear till (m = 3, B = 0.35): (a) several glacier surfaces are shown at approximately two dimensionless time unit intervals.Selected front positions and surface shapes marked (c), (d) and (e) are shown in more detail in subfigures (c-e).(b) shows the maximum height H b of the sediment bulge as a function of its horizontal position.Note how, after the initial build up phase, the bulge reaches a constant thickness.(c-e)snapshots of the front showing the surface, four selected internal layers and the ice-sediment interface at dimensionless times T equal to 0 (c), 16 (d) and 32 (e).The dots indicate material particles and can be followed throughout the three pictures.The connecting lines illustrate how material layers (fine solid and dashed lines) with initial shape as in (c) become deformed with time.

Fig. 5 .Fig. 5 .
Fig. 5. Experiment B: Relatively soft and moderately non-linear till (m=3, B=0.35): Deformation of originally parallel internal sediment layers are shown for dimensionless times T equal to 0 (a), 8 (b) and 16 (c).The dots represent material particles throughout the pictures.The topmost layer shows the strongest deformation due to the uplifting of sediment needed to form the bulge (b-d).After a while the bulge is stable in amplitude and shape, and only a small uplift of the sediment takes place at the down-facing side of the bulge (c).Three former vertical straight lines (solid) illustrate how such lines deform with time.

Fig. 6 .
Fig. 6.Experiment B: Relatively soft and moderately non-linear till (m = 3, B = 0.35): Three snapshots (a-c) at intervals of two dimensionless time units T show the evolution of five internal sediment layers which are parallel at the time T = 15 (a), when the sediment bulge had an almost stable slope.Sediment uplift takes place at the transition from the down-facing bulge slope to the undisturbed sediment layer (b-c).A vertical straight line in (a) is followed throughout (b-c).

Fig. 7 .
Fig. 7. Experiment B: Relatively soft and moderately non-linear till (m = 3, B = 0.35): Snapshots at dimensionless time T = 15 showing the distribution of horizontal velocity (a), vertical velocity (b) and effective stress (c) in the frontal part of the glacier.In (b) the dashed line represents the zero velocity contour line and the symbol + means positive vertical velocities, that is upward movement.

Fig. 8 .Fig. 8 .
Fig. 8. Experiment A: Relatively stiff and moderately non-linear till (m=3, B=8): Snap mensionless time T =15.4 showing the distribution of horizontal velocity (a), vertical veloc effective stress (c) in the frontal part of the glacier.Note that the effective stress in the decreasing distal of the glacier front.34 Fig. 8. Experiment A: Relatively stiff and moderately non-linear till (m = 3, B = 8): Snapshots at dimensionless time T = 15.4 showing the distribution of horizontal velocity (a), vertical velocity (b) and effective stress (c) in the frontal part of the glacier.Note that the effective stress in the sediment is decreasing distal of the glacier front.

Fig. 9 .
Fig. 9. Experiment B: Relatively soft and moderately non-linear till (m = 3, B = 0.35): Velocity vectors of the sediment bulge are shown at the dimensionless time T = 15.Velocity vectors are mostly parallel to the bed.Sediment uplift takes place at the downfacing side of the bulge.

Fig. 10 .
Fig. 10.Experiment B: Relatively soft and moderately non-linear till (m=3, B=0.35 profiles of the non-dimensional horizontal velocity U (lines with crosses) and effective with dots) are shown at a dimensionless time T =15 for the three locations A, B and Note that the vertical scale changes for each profile.(b) shows the ratio between the ho at the ice-sediment interface U b to the one at the glacier surface U s .

Fig. 10 .
Fig. 10.Experiment B: Relatively soft and moderately non-linear till (m = 3, B = 0.35): (c-e) vertical profiles of the non-dimensional horizontal velocity U (lines with crosses) and effective stress T eff (lines with dots) are shown at a dimensionless time T = 15 for the three locations A, B and C shown in (a).Note that the vertical scale changes for each profile.(b) shows the ratio between the horizontal velocity at the ice-sediment interface U b to the one at the glacier surface U s .

GFig. 11 .
Fig. 11.Experiment C: Moderately stiff and moderately non-linear till (m = 3, B = 0.79): (a) glacier surfaces are shown at time intervals of approximately two dimensionless time units.Selected front positions marked c, d and e are shown in more detail in subfigures (c-e).(b) shows the maximum height H b of the sediment bulge with its horizontal position.Note that the bulge is significantly smaller than in Fig. 4. (c-e) snapshots of the front showing the surface, four selected internal layers and the sediment surface at dimensionless times T equal to 0 (c), 8 (d) and 14 (e).The dots indicate material particles that can be followed throughout the three pictures.The connecting lines illustrate how material layers (fine solid and dashed lines) with initial shape as in (c) become deformed with time.The arrow marks the initial position of the glacier front.

Fig. 12 .
Fig. 12. Experiment C: Moderately stiff and moderately non-linear till (m=3, B=0.79): Sed which were originally parallel to the bed (a) are deformed and form a wave to keep the sed in position (b-c).The snapshots were taken at dimensionless times T equal to 0 (a), 8 (b) Three vertical straight lines in (a) are followed throughout (b-c). 38

Fig. 12 .
Fig. 12. Experiment C: Moderately stiff and moderately non-linear till (m = 3, B = 0.79): Sediment layers which were originally parallel to the bed (a) are deformed and form a wave to keep the sediment bulge in position (b-c).The snapshots were taken at dimensionless times T equal to 0 (a), 8 (b) and 14 (c).Three vertical straight lines in (a) are followed throughout (b-c).

Fig. 13 .
Fig. 13.Experiment C: Moderately stiff and moderately non-linear till (m = 3, B = 0.79): (c-e) vertical profiles of the non-dimensional horizontal velocity U (lines with crosses) and effective stress T eff (lines with dots) are shown at a dimensionless time T = 15 for the three locations A, B and C shown in (a).Note that the vertical scale changes for each profile.(b) shows the ratio between the horizontal velocity at the ice-sediment interface U b to the one at the glacier surface U s .

Fig. 14 .Fig. 14 .
Fig. 14.Moderately stiff and highly non-linear till (m=40, B=0.42): (c-e) vertica horizontal velocity (lines with crosses) and effective stress (lines with dots) are shown at time T =15 for the three locations A, B and C shown in (a).(b) shows the ratio betwee velocity at the ice-sediment interface U b to the one at the glacier surface U s .
could be interpreted as the inverse layering obtained from overfolding within the glacier front.Another www.the-cryosphere.net/4 where [h] is a typical mean thickness of a glacier, [u s ] a typical mean horizontal surface velocity and [τ d ] = ρg[h]sin(α s ) the driving stress assuming the same surface slope α s = 5 • as the [h] corresponds to the thickness of the ice slab and [u s ] to the analytical solution for the surface velocity (Eq.A1).The non-dimensional ice thickness H and mean ice thickness [H ] are defined as h o /[h] and [h o ]/[h], respectively, where h o is the local and [h o ] the mean ice thickness of the observed or modelled glacier.

Table 1 .
Flow parameters used in the sediment layer for the different experiments.

Table 2 .
Contribution of sediment deformation at positions A, B, C in percentage to the total deformation.