Field equations

The problem to be solved is the two-dimensional 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 parallel-sided slab geometry and use a Cartesian coordinate system with the x axis along the bed pointing down-glacier, and the z axis normal to the bed pointing upward. The constitutive relationship used is

where ἑ_ij and τ_ij are the components of the strain-rate and deviatoric-stress tensors, respectively, and τ_e is the effective stress. The rate factor A and the flow-law exponent n are material parameters. We use the commonly accepted value of n = 3 for the flow-law exponent, and use a rate factor of A = 10 MPa–₃ 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–₃ a–₁ which was used by Reference GudmundssonGudmundsson (1999) in a three-dimensional numerical modeling of Unteraargletscher, Bernese Alps, Switzerland.

Numerical scheme and model description

Equations (1–3) were solved with a finite-element method that was coded for this study. Figure 1a shows the finite-element mesh used in the computation. The geometry of the 6 km central part of the finite-element grid is a simplified longitudinal cross-section 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 strain-rate 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.

Fig. 1. Finite-element mesh used in this study. (a) All elements and coordinates used for the longitudinal cross-section model of a conceptional glacier with its head at x=0m. (b) Enlargement of the mesh near the glacier bed. Shaded region indicates the 1m thick subglacial thin layer used to introduce basal sliding.

The discretization error was estimated by comparing numerical results with the analytical solution for plane-slab 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 sub-basal layer within the central part of the finite-element 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 sub-grid bed roughness. Spatially varying basal slipperiness is generated by giving different viscosity values to some of the elements of the sub-basal layer. A no-slip boundary condition is applied for the bottom nodes of the sub-basal layer (Fig. 1b).

Ignoring the effects of longitudinal stresses within the sub-basal layer, the deformational velocity of the layer is related to the viscosity of the sub-basal 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 sub-basal 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 finite-element mesh to evolve with time.

Stress-free conditions are imposed at the surface, and atmospheric pressure is neglected (τ_ijn _j = 0, where nj is a surface-normal 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).

Transmission of basal sliding to the surface

Figure 2 shows model results for which the rate factors of the elements in the sub-basal layer are set to 25 MPa–₁ a–₁ for a slippery zone within 7000 < x < 9000 m, and 5 x 10–₄ MPa–₁ a–₁ 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 down-glacier, 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 small-amplitude 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).

Fig. 2. Computed fields of (a) horizontal velocity, (b) vertical velocity, (c) longitudinal deviatoric stress component τxx, and (d) vertical strain rate. Only the central domain of the modeled field (5000 <x< 11000 m) is shown. Contour lines were established using the method of triangulation with linear interpolation.

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₁ a₁ at 7000 < x < 9000 m). We define the transmission rate of horizontal speed as

Fig. 3. (a) Longitudinal distribution of the horizontal velocity obtained for the sliding condition described in the text. Surface (solid circle) and basal velocities (open circle) were computed at nodes and midpoints of the finite-element mesh. Solid line indicates the surface velocity for the non-slip condition. (b) Horizontal surface velocities at x = 8000 plotted against the basal velocity u_b at x = 8000 and the length of the slippery zone L. The surface velocity is indicated for each experiment with a unit of m a–₁ and contoured with solid lines. The shaded contours show the transmission rate of horizontal velocity f_u.

where Δu_s is the difference in horizontal surface velocity between non-slip 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 sub-basal 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–₁ 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 non-linearity 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 sub-basal 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

Fig. 4. (a) Longitudinal distribution of the surface (solid circle) and basal velocities (open circle) obtained for the condition given by Equation (5). The solid line is the surface velocity for the non-slip condition. (b) Longitudinal gradients of the horizontal surface velocity at x =8000 plotted against basal velocity gradient du^b/dx at x = 8000 and the length the slippery zone L. The surface velocity gradient is denoted each experiment with a unit of 10^–3 a^–1 and contoured with lines. The shaded contours are the transmission rates of horizontal velocity gradient f^grad.

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 water-pressure 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 con­solid 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 down-glacier 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 sub-basal 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⁰ 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 surface-melting rate decreases, the basal velocity starts to decrease and the zone moves down-glacier during the period 1800–0000h. Eventually, this down-glacier propagation of the slippery zone replaces the compressive flow at the study site with a weaker extending flow (0300–0600h).

Fig. 5. Diurnal evolutions of (a) surface and (b) basal velocity distributions computed with the basal condition scenario described in the text.

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.

Fig. 6. (a) Diurnal variations of (a) horizontal surface velocity and (b) vertical strain rate. Open circles with the dashed line are model results, and the solid lines are measurement data taken on Unteraargletscher on 25 August 2001. Vertical strain rate is computed as the mean for the upper 200 m of the glacier, and measurements were taken from a 174 m deep borehole and plotted with the scale magnified 10 times.

The fact that the computed strain rate accounts for only about 10% of the measured borehole-depth 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 strain-rate 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 τ₀ 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.1 MPa so that ice becomes less viscous under a low-stress condition, and applied the basal conditions identical to those used for Figure 2. Computed flow velocities, stress and strain-rate 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 flow-law 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.

Fig. 7. Computed fields of (a) horizontal velocity, (b) vertical velocity, (c) longitudinal deviatoric stress component τ_xx, and (d) vertical strain rate. Modified flow law (Equation (7)) was used with the basal conditions identical to those in Figure 2.

The studied velocity variation accompanied by a strain-rate anomaly is similar in nature to the “mini-surges” 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 fast-flow 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–₁ that is higher than our modeling result in Figure 6b.