Introduction
Temperate glaciers are known to change their flow velocity over a timescale of a week, a day and in some cases hours (Reference WillisWillis, 1995). Such relatively shortterm fluctuations are understood to be controlled by changes in subglacial conditions, particularly by subglacial waterpressure fluctuations (Reference Iken and BindschadlerIken and Bindschadler, 1986; Reference JanssonJansson, 1995). The transmission of the basal variability in flow speed to the surface is complex and at first sight counterintuitive (Reference BaliseBalise and Raymond, 1985), and the influence of temporal changes in basal lubrication on surface velocities and glacial deformation difficult to predict in detail.
In this study, we use a numerical finiteelement model to determine how a nonlinear viscous medium (ice) moving over a bed with spatially nonuniform basal lubrication causes flow perturbations. The aim of this work is not to comprehensively study the general aspect of the transmission of basal slipperiness perturbations to the surface of glaciers, but rather to find basal conditions that can reproduce the particular diurnal fluctuations in surface flow velocity and vertical strain observed during the ablation season in 2001 on Lauteraargletscher, Bernese Alps, Switzerland.
Methods
Field equations
The problem to be solved is the twodimensional flow of incompressible viscous material governed by Stokes equations:
where u = (u, w) is the velocity vector, p is the hydrostatic pressure, μ is the viscosity, p is the density of ice, and g = (g_{x}, g_{z}) is the gravity vector. We employ parallelsided slab geometry and use a Cartesian coordinate system with the x axis along the bed pointing downglacier, and the z axis normal to the bed pointing upward. The constitutive relationship used is
where ἑ_{ij} and τ_{ij} are the components of the strainrate and deviatoricstress tensors, respectively, and τ_{e} is the effective stress. The rate factor A and the flowlaw exponent n are material parameters. We use the commonly accepted value of n = 3 for the flowlaw exponent, and use a rate factor of A = 10 MPa–_{3} a_{–1}, which gives modeled surface velocities similar to those observed. Taking into account the reduction in shear stress due to the valley shape by using a shape factor of 0.5 (Reference NyeNye, 1965), this value of the rate factor compares favorably with 75 MPa–_{3} a–_{1} which was used by Reference GudmundssonGudmundsson (1999) in a threedimensional numerical modeling of Unteraargletscher, Bernese Alps, Switzerland.
Numerical scheme and model description
Equations (1–3) were solved with a finiteelement method that was coded for this study. Figure 1a shows the finiteelement mesh used in the computation. The geometry of the 6 km central part of the finiteelement grid is a simplified longitudinal crosssection of Lauteraargletscher. We assumed an ice thickness of 400 m and a surface slope of 4° (Reference Funk, Gudmundsson and HermannFunk and others, 1995; Reference BauderBauder, 2001). To minimize boundary effects, the numerical grid was extended an additional 6 km upstream and 6 km downstream from this central section. We used Galerkin’s method with quadratic shape functions (Reference ZienkiewiczZienkiewicz, 1977) to compute the velocity and pressure fields. Starting from the solution of a linearly viscous flow, a new effective viscosity distribution was calculated from the previously determined strainrate field, and the calculation repeated until the horizontal velocity field converged within 10–5 ma–1. The computation gives the velocity components at each node and midpoint, and the maximum resolution of the velocity field is 250 m longitudinally and 25 m vertically.
The discretization error was estimated by comparing numerical results with the analytical solution for planeslab flow for no slip. The largest error of 1.2% occurred in the horizontal surface velocity. By comparing the numerical and the analytical vertical velocity profiles, we found that the main source of error was near the bed, where the largest spatial gradients exist.
Boundary conditions
Basal motion is simulated by introducing a thin subbasal layer within the central part of the finiteelement model. The viscosity of this layer can be regarded as the effective viscosity of subglacial till, or, which in the case of Lauteraargletscher is more appropriate, as being related to the form drag generated by subgrid bed roughness. Spatially varying basal slipperiness is generated by giving different viscosity values to some of the elements of the subbasal layer. A noslip boundary condition is applied for the bottom nodes of the subbasal layer (Fig. 1b).
Ignoring the effects of longitudinal stresses within the subbasal layer, the deformational velocity of the layer is related to the viscosity of the subbasal layer, its thickness and the shear stress at the top of the layer. This gives rise to a sliding law of the form
where ub is the basal sliding velocity, d is the layer thickness, τb is the basal shear traction, and A' and n' are the parameters of the layer. We assume n' = 1. This method of introducing basal motion through a subbasal layer has been used before (e.g. Reference Vieli, Funk and BlatterVieli and others, 2000). Because longitudinal stress gradients will cause the thickness of the layer to change with time, this method is not suitable for transient calculations. Here we are interested in obtaining snapshots of the velocity field for a given glacier geometry, so we do not allow the surface or other boundaries within the finiteelement mesh to evolve with time.
Stressfree conditions are imposed at the surface, and atmospheric pressure is neglected (τ_{ij}n _{j} = 0, where nj is a surfacenormal vector component). The lower end of the modeled domain (x = 16 000 m) was connected to the upper end (x = 0 m) so that the glacier is assumed to be infinitely long (periodical boundary condition).
Results and Discussion
Transmission of basal sliding to the surface
Figure 2 shows model results for which the rate factors of the elements in the subbasal layer are set to 25 MPa–_{1} a–_{1} for a slippery zone within 7000 < x < 9000 m, and 5 x 10–_{4} MPa–_{1} a–_{1} elsewhere. Introducing this slippery zone in the model increased the surface horizontal velocities directly above the center of the zone by about 3 0% (Fig. 2a). With increasing distance up and downglacier, the amplitude of the velocity perturbation decreases. For example, five ice thicknesses away from the boundaries of the slippery zone (x = 5000 and 11000 m), the increase in surface velocities is down to 8%. The calculated vertical velocity field shows extending and compressive flows above the boundaries (Fig. 2b). Furthermore, tensile and compressive deviatoric stresses occur throughout the entire glacier but are strongest near the surface (Fig. 2c). Because the deviatoric stress is related to the vertical strain rate by Equations (1) and (3), this stress field generates vertical strain rate as shown in Figure 2d. The computed vertical strain rate near the surface helps to explain the observed vertical strain fluctuations discussed in the next subsection. Analytical solutions for smallamplitude slipperiness perturbations show a similar increase in deviatoric stress and strain with distance away from the source of the disturbance (i.e. the slippery zone) (Reference VonmoosVonmoos, 1999).
Figure 3a shows the longitudinal distribution of basal and surface horizontal velocities calculated from the same basal condition as used in Figure 2 (A' = 25MPa_{1} a_{1} at 7000 < x < 9000 m). We define the transmission rate of horizontal speed as
where Δu_{s} is the difference in horizontal surface velocity between nonslip and slip conditions and u_{b} is the basal velocity at the same position (Fig. 3a). Calculations were performed for various sets of viscosity of the subbasal layer and different lengths of the slippery zone (L). Results are summarized in Figure 3b, which shows the horizontal surface flow speed u_{s} and transmission rate of horizontal velocity f_{u} at x = 8000 m as functions of L and the basal velocity u_{b} As expected, the surface velocity increases with L. The transmission rate of horizontal speed f_{u} also increases with L, and becomes larger than 0.75 for u_{b} > 5 ma–_{1} when L = 6000 m. In general, f_{u} varies with L similar to an analysis for a linear medium (Reference BaliseBalise and Raymond, 1985), although our calculated values of f_{u} are slightly smaller. This discrepancy might be caused by an effective material softening close to the glacier bed that introduces nonlinearity into the flow law. Reference Blatter, Clarke and ColingeBlatter and others (1998) used the boundary condition of zero basal shear traction (free slip) in their numerical model, which resulted in much higher transmission rates than ours.
In order to examine the transfer of a velocity gradient from the bed to the surface, we assumed the rate factor of the subbasal layer to vary linearly with x within a region of length L in further numerical experiments. Figure 4a shows the resulting surface and basal horizontal velocity distribution for A' given by
Both the longitudinal gradients of basal and surface velocities were calculated at x = 8000 m, as well as the transmission rate of horizontal velocity gradient defined as
Figure 4b shows a contour plot of the surface velocity gradient du^{s}/dx and f^{grad} for a range of lengths L and basal velocity gradients du^{b}/dx. The surface velocity gradient increases as the basal gradient and L increase, but the transmission rate is rather small. In this case, the maximum value of f^{grad} is about 0.5, and almost no transmission occurs for L 5 3000 m.
Diurnal fluctuation in the surface flow speed and vertical strain
We now compare our modeling predictions to measurements of diurnal velocity and strain fluctuations made on Lauteraargletscher from 22 to 27 August 2001. Lauteraargletscher, a tributary of Unteraargletscher in the Bernese Alps of Switzerland, is a temperate valley glacier that is about 5km longand1 kmwide. Details of the measurements are reported by Reference Sugiyama and GudmundssonSugiyama and Gudmundsson (2003), but the main findings can be summarized as follows.

(A) Surface flow velocity at the center of the glacier varied diurnally, with an approximate minimum speed of 30 ma^{–1} in the morning and a maximum of60 ma^{–1} in the afternoon.

(B) From 24 to 26 August, diurnal fluctuations in vertical strain occurred in a 174 m deep borehole where the total ice thickness was about 400 m. The depth of the borehole increased by about 30 mm from morning to evening, and then decreased again by a similar amount towards the next morning. This vertical stretching/compression corresponds to a strain of 1.7610^{–4}.
Furthermore:

(a) The surface velocity and vertical strain measurements were made within the ablation area with the snowline about1 km upstream.

(b) Subglacial water pressure at the study site fluctuated diurnally throughout the summer, with an average pressure decreasing from July to August.

(c) Surface velocity varied diurnally at least as far as ±1.5km up and downstream from the study site, and the amplitude of the diurnal variation increased in the upstream direction.
If we assume no transverse component in the flow velocity, observation (B) suggests a strongly compressive flow at the study site in the daytime. This compressive flow could be due to the comparatively larger increase in flow velocities in the upstream direction in the early morning (c). The fact that the surface flow velocity was well correlated with the subglacial water pressure during the study period suggests that the surface velocity variation can be attributed to the diurnal waterpressure variations. The decrease in the average water pressure in the ablation season (b) seems to be related to the development of the subglacial drainage system beneath the study site. Because the drainage system starts to develop in the lower reach of the glacier and then progresses upward during the ablation season, water pressure fluctuates with larger amplitude and causes greater basal sliding in the upper reach. Therefore, we ascribe the observed vertical strain to the longitudinal gradient of basal velocity.
From Figure 4b, one can determine the basal velocity sliding gradient necessary to produce the observed vertical strain. However, the observed strain rate was about 0.1a^{–1}, which is too large to be generated by the velocity gradients in Figure 4b. This suggests that the perturbation of slipperiness was of even larger than the one used to produce Figure 4b.
for The observed diurnal variation in flow speed and its consolid nection to basal pressure variations suggests that the glacier became partly decoupled from the bed during daytime. This may be due to increased basal cavity formation or due to the failure of a thin till layer. If such a decoupled zone begins to develop upstream from the measurement site and then increases in spatial extent and sliding magnitude with moving downglacier in the course of the day, large deviatoric stresses and vertical strain rates near the surface can result that are consistent with the results in Figure 2c and d. Such a basal condition can be simulated by prescribing a slippery zone as described in the previous section. The sliding magnitude, spatial size and position of the slippery zone were varied to fit the observed velocity and strain rate at the study site as well as possible. Because the observed diurnal changes in thickness are only a few centimeters at most, changes in thickness and surface slopes can be ignored in the modeling calculations, and the diurnal variations can be analyzed without changing the model geometry.
We introduced a slippery zone upstream of the observation site by changing the rate factor of the subbasal layer as explained above. By subsequently changing the magnitude, size and position of the slippery zone, the observed spatial and temporal pattern in vertical strain rate could be reproduced qualitatively. Figure 5a and b show surface and basal velocity distributions calculated for the determined slippery zones. The conditions of the slippery zone at each time were determined by adjusting the distribution of A^{0} to fit the observed surface velocity at the study site. The temporal evolution of the slippery zone is explained as follows: The slippery zone starts to develop upstream of the observation site at about 0800 h when the surface melting starts. Then it increases in extent and magnitude, and forms a sliding boundary beneath the study site at 1500 h. In the evening, when the surfacemelting rate decreases, the basal velocity starts to decrease and the zone moves downglacier during the period 1800–0000h. Eventually, this downglacier propagation of the slippery zone replaces the compressive flow at the study site with a weaker extending flow (0300–0600h).
Figure 6 compares the field measurements with the computed temporal variations of horizontal surface velocity and vertical strain rate in the upper 200 m of the ice at the study site (x = 8000 m). The model reproduces the observed temporal variation in velocity and vertical strain on 25 August 2001. Even though the computed strain rates are still considerably smaller than the measured ones, they are an order of magnitude larger than those obtained in Figure 4b.
The fact that the computed strain rate accounts for only about 10% of the measured boreholedepth changes shows that some important aspects of the flow are not described adequately by the model. To obtain larger strain rates, the model could assume higher flow speeds in the upper region, but this assumption would contradict observed surface velocities. Another possibility is that the flow law used is incorrect. In the present model, stress concentration near the surface in Figure 2c does not generate clear strainrate concentration in Figure 2d. This is because Glen’s flow law (Equation (3)) gives much higher effective viscosity as the effective stress decreases near the glacier surface. To assess the influence of effective viscosity near the surface, another experiment was carried out with modified Glen’s flow law,
Here, as elsewhere, a constant value τ_{0} has been introduced in the flow law to avoid mathematical singularity where the effective viscosity is infinite when the effective stress vanishes (Reference HutterHutter, 1983). We took τ_{0} = 0.1 MPa so that ice becomes less viscous under a lowstress condition, and applied the basal conditions identical to those used for Figure 2. Computed flow velocities, stress and strainrate fields are shown in Figure 7 and compared with Figure 2. It is interesting to see the vertical strain increases near the surface while the deviatoric stress decreases (Fig. 7c and d). The results suggest that the modeled vertical strain will be improved by Equation (7) but will still be insufficient due to the reduction in the stress. Further study on the flowlaw exponent or flow law itself is required to interpret the field data adequately. The model can also be improved by including the average longitudinal strain and transverse stresses over the area in question.
The studied velocity variation accompanied by a strainrate anomaly is similar in nature to the “minisurges” or propagation of enhanced motion waves, observed in Variegated Glacier, Alaska, U.S.A. (Reference Raymond and MaloneRaymond and Malone, 1986; Reference Kamb and EngelhardtKamb and Engelhardt, 1987). Reference BaliseBalise and Raymond (1985) applied their analysis of a linear viscous flow to interpret the measurements on the surface by a propagation of basal velocity anomaly that was much sharper than the one we used for Figure 5b. They also suggested that the vertical straining of ice might explain uplift events, which are increases in surface elevation during a fastflow period. In order to apply our model to uplift events, quantitative discussion described above is crucial because the uplift rates measured at Variegated Glacier (Reference Kamb and EngelhardtKamb and Engelhardt, 1987), or at Unteraargletscher by Reference Iken, Röthlisberger, Flotron and HaeberliIken and others (1983), are in the same order of ∼10_{–4} day–_{1} that is higher than our modeling result in Figure 6b.
Summary and Conclusion
We used a finiteelement model with Glen’s flow law to study the effects of a local perturbation in basal slipperiness on the surface velocity and the englacial strainrate distribution. The model results were compared to recent measurements of diurnal strainrate variations on Lauteraargletscher. Our model of downglacier propagation of a slippery zone was able to qualitatively explain the observed temporal variation in flow velocities and vertical strain rates. However, quantitative comparison between model predictions and field measurements showed that calculated vertical strain rates were about an order of magnitude too small. The reason for this discrepancy remains unclear. A possible explanation is that the modeled ice near the surface is too stiff, which may arise from ignoring the average longitudinal strain or because the flowlaw exponent of ice varies with stress.
Acknowledgements
We are grateful to R. Naruse of the Institute of Low Temperature Science, Hokkaido University, for his support for this study. All the members of Section of Glaciology, VAWETHZ, Zürich, assisted us in the fieldwork at the glacier. We also thank L. Tarasov, H. Blatter, R. C. A. Hindmarch, D. Vaughan and an anonymous reviewer for their beneficial comments on the manuscript. The Swiss National Science Foundation through grant No. 2100063770.00 and the Inoue Scientific Field Study Foundation funded this research.