Introduction
A glacier-dammed marginal lake fills until the water reaches a certain threshold, at which the lake water connects to the subglacial drainage system, allowing it to drain. The drainage is of unstable nature (Reference NyeNye, 1976) and may lead to potentially hazardous flooding in the glacier outlet stream. The phenomenon, termed ‘jökulhlaup’ (e.g. Reference BjörnssonBjörnsson, 1992), drives large amounts of water into the basal environment of the adjacent glacier within a short time, causing a strong perturbation of subglacial water flow. Increased water pressure and temporary water storage at the glacier bed locally enhance basal motion, which is reflected in increased surface flow speeds (e.g. Reference IkenIken, 1981; Reference Sugiyama and GudmundssonSugiyama and Gudmundsson, 2004). Reference Anderson, Walder, Anderson, Trabant and FountainAnderson and others (2005) emphasized the coherence between the speed-up of Kennicott Glacier, Alaska, and the drainage of Hidden Creek Lake (HCL), whereas Reference WalderWalder and others (2005, Reference Walder2006) discussed the response of the ice dam in the vicinity of HCL to its filling and drainage.
Gornergletscher, Switzerland, is a valley glacier located in the southern Swiss Alps. At the place where it meets Grenzgletscher, an even larger tributary, it dams the ice-marginal lake, Gornersee (Fig. 1). Gornersee originates from meltwater retention during spring and usually drains in early summer. The bottom of Gornersee is essentially ice, except along the eastern lake border, where the lake floor is formed by a steep slope of morainal debris plunging westward into the lake. Five drainage events of Gornersee in past years (2004–08) have been observed. Results of the impact of individual Gornersee jökulhlaups on the ice flow and hydrology of Gornergletscher have been published by Reference Huss, Bauder, Werder, Funk and HockHuss and others (2007) and Reference Sugiyama, Bauder, Zahno and FunkSugiyama and others (2007a, Reference Sugiyama, Bauder, Huss, Riesen and Funk2008). The results from passive seismic surveys are described by Reference Walter, Deichmann and FunkWalter and others (2008). Reference Werder and FunkWerder and Funk (2009) and Reference Werder, Loye and FunkWerder and others (2009) give details of extensive dye-tracing experiments performed during the observed Gornersee drainage events.
(a) Map of Gornersee (grey) and Gornergletscher with contours of ice surface elevation (solid, 50 m intervals), glacier margin (thick solid), location of the tachymeter and nearby GPS reference station (triangle TS), boreholes bh430 and bh6 (squares) and ice-flow measurement locations (crosses). The dashed box indicates the geographic coverage of Figures 4–7. (b) Zoom panel of the box in (a) (ice dam), indicating the measurement locations in the ice-dam area. For identification, we label each marker with a number composed of the survey year (first digit: 4 or 6) and an individual integer (second digit). A ‘g’ in front of the integer indicates a GPS-equipped marker.
In this paper, we further examine the interaction between draining lake and damming glacier using field measurements and numerical modelling. We measured ice-flow velocities of Gornergletscher in the vicinity and down-glacier of the Gornersee marginal lake during the 2004 and 2006 drainage events and studied the observed changes in flow speed and direction during these distinct events. To reproduce and interpret the observed flow changes we used a full-Stokes three-dimensional (3-D) numerical ice-flow model, which is able to treat local variations in basal motion, as well as the effect of lake water pressure acting on the ice. To detect and corroborate the mechanisms responsible for the observed flow changes, we carried out different flow simulations with adjustments in the model boundary conditions.
Field Methods
We measured the ice flow of Gornergletscher by surveying the displacements of eight poles (markers) installed in the ice. Five markers (41–45 in 2004, 61–65 in 2006) were located in the vicinity of Gornersee (henceforth ‘ice dam’; Fig. 1b). These markers were equipped with optical survey prisms. Every hour an automated tachymeter estimated distance and both horizontal and vertical angles from the prisms’ reflectance patterns. This tachymeter was situated on the rock ledge at the northern margin of the glacier (Fig. 1a; marked ‘TS’). We fixed multiple reference prisms onto rock outcrops along the glacier margin to carry out an atmospheric correction for the angle and distance readings. The relative azimuth and distance data were then converted to absolute coordinate positions. For the optical survey, the position errors are dependent on the tachymeter orientation (Table 1). The error is larger in the down-glacier (east–west) direction than in the cross-glacier (north–south) direction, due to increased influence of vertical tilt errors of the tachymeter device in this axis. The absolute position error of 35 cm in the down-glacier direction is large, but diminishes for repeated measurements. Within 3–4 days, when N = 50 position estimates are acquired, it decreases to ∼5 cm. For our studies, N is in the range 50–350, hence the error becomes even smaller. Also, the ice displacement at the key markers 41 and 42 (61 and 62) is approximately aligned with the axis of minor errors. Three markers (46–48 in 2004, 66–68 in 2006) were also installed ∼1–2.5 km down-glacier of Gornersee (Fig. 1a). These stakes moved down-glacier, i.e. principally west, and differential GPS receivers (Leica 500) at these marker locations were used to overcome the low accuracy of the tachymeter in the down-glacier direction. The GPS reference station was installed close to the tachymeter. We recorded 1 hour sessions of L1 + L2 phase signals every 3 hours at a sampling rate of 1/30 Hz. We processed the GPS data in static differential mode relative to the reference station, using commercial software provided by Leica. The error of the position estimates using GPS is ±1 cm in the horizontal, and independent of orientation.
Absolute position errors (cm) for the optical ice-flow measurements (tachymeter) with standard deviation, SD, and standard error from the mean,
for N = {50, 150, 3001} repeated measurements
Note that the data of an individual marker were acquired by either GPS or tachymeter and not the result of a combined survey. However, we compared trajectories of poles installed side by side and surveyed by either GPS or tachymeter. These datasets were generally coherent, inspiring confidence in the trajectories obtained by means of the optical survey. The horizontal flow velocity, v, of an individual stake was computed as the least-squares fit of a first-degree polynomial x(t) = v · t, with x(t) = xi (ti ) the vector of N measured positions in the time interval Δt = tN − t 1. Thus, the velocity corresponds to the slope of a linear fit to the stake’s displacement trajectory in the time interval Δt.
In addition to the geodetic measurement, we measured basal water pressure in boreholes drilled to the glacier bed using a hot-water drilling system. In each borehole, a vibrating-wire pressure sensor (Geokon Model 4500) connected to a logger at the glacier surface was operated at a sampling interval of 10 min. We show data from borehole bh430 which was used in 2004, and from borehole bh6, used in 2006; their respective locations are depicted in Figure 1. Also, the lake water level was recorded using a Keller PAA-36W pressure transducer with a 10 min sampling interval.
Numerical Model
We solve the balance equations for incompressible, isothermal steady Stokes flow as
where ρ = 900 kg m−3 is the density of ice, g is the gravitational acceleration and t is the deviatoric stress tensor with tr(t) = 0. The fields to solve for are velocity, v, and pressure, p. The constitutive relation adopted for the stress is a power law of the form
where
is the strain-rate tensor, n is the power-law exponent and A is the constant rate factor. The problem given by Equations (1–3) was solved by adapting the variable viscosity, η(v), through a Newton–Raphson iteration scheme until convergence of the velocities was achieved, starting from an initially uniform viscosity distribution or a previously determined solution. A penalty method was used for the incompressibility constraint. The model is implemented in the commercial finite-elements software MSC.Marc (MSC Software Corporation, 2005) and has been used and validated in previous studies (Reference GudmundssonGudmundsson, 1994, Reference Gudmundsson1999; Reference HelbingHelbing, 2005; Reference Raymond and GudmundssonRaymond and Gudmundsson, 2005). In all computations we used n = 3 and A = 75 Mpa−3 a−1 (Reference GudmundssonGudmundsson, 1999).
From a digital elevation model (DEM), acquired in September 2004 using airborne photogrammetry, and radio-echo sounding measurements of the bed topography (Reference Sugiyama, Bauder, Huss, Riesen and FunkSugiyama and others, 2008), we generated a mesh representing the Gornergletscher geometry. The model domain includes the confluence area and the lower tongue of Gornergletscher. The domain does not cover the entire glacier system with the Grenzgletscher and Gornergletscher tributaries, due to a lack of information about bed geometry upstream of the confluence. On the two-dimensional (2-D) transverse sections where the model boundary traverses the Grenzgletscher and Gornergletscher tributaries, we prescribed fixed flow velocities as inflow boundary conditions at the element nodes. These flow velocities were derived from 2-D cross-section flow fields, separately computed on identical geometries as the inflow faces, using the implementation of 2-D glacier flow in a cross section of Reference Sugiyama, Bauder, Zahno and FunkSugiyama and others (2007a) with equivalent parameters, n and A. The 2-D flow fields at the inflow boundaries were tuned so the computed surface flow-speed profiles matched nearby measurements of annual surface flow velocities. Note that the model of Reference Sugiyama, Bauder, Zahno and FunkSugiyama and others (2007a) was solved using a finite-differences grid. The nodal velocities input to the finite-element model were interpolated from the separately computed flow fields, then translated to the 3-D model faces, aligned to the mean slope of the glacier surface at the transverse section, and oriented in the main flow direction. The slopes were determined from the DEM.
In the vertical, the mesh consisted of six layers of elements (3330 total) of decreasing thickness towards the bottom. A Dirichlet condition, v = 0, was applied at the bottom nodes of the sixth element layer and at all marginal nodes. The glacier bed is represented by the plane consisting of the upper faces of the sixth-layer elements. Note that this sixth layer of elements has no physical meaning otherwise, but we refer to it as the ‘bed layer’ from now on. We applied a parameterization of basal motion, utilizing deformation of the bed-layer elements. In this way, the basal velocity, u b, at the glacier bed, i.e. at an upper face of an element of the bed layer, corresponds to
where d is the element thickness, η is the element viscosity and τ b is the basal shear stress. We used a Newtonian rheology with η = constant for the bed-layer elements and substituted C = d/η into Equation (4), so the linear relation for the basal boundary condition was incorporated as
Since we only solved for velocities and computed stresses separately, we controlled the simulation of the influence of the lake drainage directly through enhancement or reduction of basal velocities. Thus, we increased or decreased C in subsets of the bed-layer elements to locally soften or harden the bed-layer rheology, as the bed-layer elements viscosity, η, is inversely proportional to C. The thickness of the bed layer (elements) is small compared to the overburden ice thickness, but a moderate thickness was needed when large deformations of the elements occurred. We used a uniform thickness of d = 10 m, which is ∼3% of the maximum ice thickness in the central confluence. Another description of this parameterization of basal motion is given by Reference Vieli, Funk and BlatterVieli and others (2000). The glacier surface was treated as a free surface with t · n = 0, where n is the surface normal. On the upper element faces within the lake domain, we described distributed pressure loads in order to simulate the effect of the lake water pressure acting on the ice surface covered by the lake water. The load normal to an element face was calculated as the pressure p f = ρ w g (h l − h f), where ρ w = 103 kg m−3 is the density of water, h f is the mean elevation of the upper element face f, and h l is the lake water-level elevation.
Results 1: Observations
We observed perturbations of the ice flow within the ice dam of Gornergletscher during two lake drainage events (in 2006 and 2004; Figs 2 and 3, respectively). Each drainage event was analysed in a discrete sequence of three successive time intervals: (1) before, (2) during and (3) towards the end of the drainage. Figures 2a and 3a show the evolution of the lake water level and borehole water level, as measured during the drainages of Gornersee. Vector plots of the flow velocities measured during the three intervals are shown in Figures 2b and 3b (ice dam) and in Figures 2c and 3c (down-glacier). Note that in the following we consider the unexpected drainage event first, and discuss the drainage of 2006 before that of 2004.
(a) Measured lake water level (dashed) and water level in borehole bh6 (solid) during the drainage of Gornersee in 2006. The drainage event is subdivided into the three intervals 1–3, as indicated by the shaded bars and vertical lines (thin solid). (b) Vector plot of horizontal flow velocities observed in the ice dam during the three intervals 1–3. (c) Vector plot of horizontal flow speeds observed down-glacier of the ice dam during the three intervals 1–3.
Same as Figure2, but for 2004.
2006 drainage event
On 5 July 2006, the lake reached a volume of ∼4 × 106 m3 and had a water level that corresponded to overflow of the ice dam at the westernmost corner of Gornersee. The overflowing water discharged into a small moulin located close to borehole bh6 (Fig. 1). The superficial lake outflow filled the moulin and balanced the filling of the lake. The lake water level did not increase further until 7 July, when it started to decrease. Between the lake and the drainage moulin, the lake water carved a gorge, 5–10 m wide, 100 m long and 30 m deep into the ice by the time the lake was fully drained. During the drainage period (>20 days) the rate of decrease of the lake water level was approximately constant. We consistently observed water draining into the moulin shaft, so we assume that the lake drained completely via overflow into the moulin. It appears that no other significant subglacial connection existed. The borehole water level in bh6 shows diurnal fluctuations with an amplitude >80 m and with peak magnitude increasing from 25 June to 5 July. The period of high water level without fluctuations between 5 and 7 July coincided with a constant lake water-level height, then the borehole water level dropped suddenly on 7 July. This is consistentwith the beginning of the lake water-level decrease and initiation of the drainage of Gornersee. From the interval 1 to 2 (denoted as ‘1 → 2’) we observed a clockwise rotation of the flow vector and a slight increase in speed at markers 61 and 62. In 2 → 3 the rotation of the flow direction continued and flow speed increased again noticeably. Note that at marker 61 the rotation of flow direction proceeds gradually during 1 → 2 → 3, whereas at marker 62 the incremental rotation in 2 → 3 is more pronounced. For 2 → 3 the measured flow velocities at 63–65 are insignificantly larger than in interval 1. Down-glacier of the ice dam, markers g67 and g68 did not show any significant changes in flow direction. Marker g66 showed a counterclockwise rotation in 1 → 2 and a slight increase in speed in 2 → 3 (Fig. 2c).
2004 drainage event
In 2004, the drainage began in a similar way to the 2006 event, i.e. following the overflow of the ice dam on 1 July (Reference Sugiyama, Bauder, Weiss and FunkSugiyama and others, 2007b). However, the superficial outflow ceased shortly afterwards and by the early morning of 3 July the lake water level decreased noticeably. Most of the lake water drained within 5 days (2–7 July). We suspect that the lake water invaded the subglacial environment of Gornergletscher at the southwestern lake border (east of markers 41 and 42; Fig. 1b). Located close to bh6, two excavated englacial channels were observed during inspection of the lake floor after drainage (Reference Sugiyama, Bauder, Huss, Riesen and FunkSugiyama and others, 2008). Unfortunately, during the evening of 5 July, ice blocks broke the cable to the pressure transducer that was recording the lake water level, which resulted in a data stoppage (Fig. 3a). Borehole bh430 in the central confluence indicates a clear perturbation of the subglacial conditions starting on 2 July, when the amplitude of diurnal water-level variations of ∼150 m abruptly decreased to ∼25 m. The water level remained consistently high during most of intervals 2 and 3, while the lake water level went down. Then the borehole water level abruptly dropped below 200 m on 7 July, apparently marking the termination of the drainage. The diurnal oscillations in basal water level recovered again from 9 July onwards. Figure 3b shows the measured flow velocities on the ice dam during the 2004 lake drainage. Note that markers 41 and 42 were located further away from Gornersee than markers 61 and 62 (2006). In 1 → 2, the ice-flow speed increased considerably and the flow direction rotated counterclockwise (down-glacier) at 41 and 42. In the same interval, the flow speeds strongly increased at markers g44 and 45 (Fig. 3b). In 2 → 3, the sense of rotation of the ice-flow direction reversed and the ice changed to flow towards Gornersee at markers 41–45. At markers 42, 43 and g44, flow speeds were highest in interval 3, whereas at 41 and 45 flow speeds decreased again after a maximum in interval 2. Down-glacier of the ice dam the measured flow velocities increased by 20–70% during interval 2 at all three locations, g46, g47 and 48 (Fig. 3c). At markers g46 and g47, the increase in flow velocities was accompanied by a deflection of the ice-flow direction northwards and a subsequent return during interval 3. In interval 3 the flow speed decreased to the same magnitude as that in interval 1 at all three markers.
Results 2: Modelling
We conducted preliminary tests in order to clarify the model sensitivity to the applied boundary conditions and to estimate errors due to numerical discretization. We varied the prescribed inflow boundary velocity distributions by both magnitude and orientation and affirmed that the computed flow regimes of Gornergletscher would not be compromised by these boundary conditions. We then constructed an additional mesh with higher vertical resolution (6105 elements, 11 layers) and verified the results of surface flow velocities. The model used is low-resolution but produces admissible horizontal surface flow regimes with minimal computation times. We performed many model runs (>2000) to systematically investigate how the free manipulation of C at different local patches of the bed layer affects the ice-flow regime. Since the basal boundary condition (Equation (5)) is linear, the variation in C operates on the magnitude of velocities, and changes in flow direction are caused mainly by influences of the geometry of Gornergletscher. Thus, the model behaviour remained manageable.
Initial flow conditions
For the subsequent simulations, an adequate flow field of Gornergletscher was required as the initial configuration. To obtain the initial flow conditions, we optimized the computed flow solution, based on all available surface flow speed observations for the period April–September 2006. The optimization was achieved by manually adjusting C in local subsets of basal elements. Varying C is an alteration of the model boundary conditions and produces systematic changes of the orientation and magnitude of the residual error vectors of the mismatch between observed and modelled flow velocities. The fitness of the distribution of C and the computed flow field is best estimated in terms of systematic errors and we seek a representative flow field where the length of the sum of normalized error vectors is minimal. The sum of normalized error vectors (SNE) is computed as
where
is the measured horizontal flow velocity,
is the computed equivalent and
is the magnitude of the measured flow speed at the ith stake location of a total N stakes. For a random distribution of error vectors, R
i
, the length R = |R| approaches zero when N → ∞ (Reference GudmundssonGudmundsson, 1999). The initial flow-field solution and the settings of C of the bed-layer elements are given in Figure 4. For this flow field, R is 0.003. For comparison purposes, R is 0.112 for a distribution with uniform C = 11.5 MPa−1 a−1 m. The computed flow field is in good agreement with the observed one. The major discrepancies occur in the ice-dam area on the side of Grenzgletscher tributary. Reference GudmundssonGudmundsson (1999) calculated flow fields for the confluence of Unteraargletscher and obtained similar residual discrepancies between modelled and observed flow regimes in the confluence area, though the error vectors point towards the junction and not down-glacier. His model systematically overestimated the down-glacier component of the flow, whereas in our model it is underestimated. Although Gudmundsson’s model did not incorporate basal motion, our results did not change when we excluded basal motion (C = 0 ≡ 10−8). However, the residual differences are not obstructive to the intended simulations. We provide additional remarks on these discrepancies in the discussion section.
(a) Comparison of measured (black) and optimum calculated (grey) summer horizontal flow velocities at a set of N = 30 stakes. The model domain outline is indicated by a dashed curve. (b) View from above of the bed-layer elements and the spatial settings of the constant C used for the flow-field computation shown in (a). The black dots correspond to the velocity measurement markers as in (a), and the outline of Gornersee is indicated.
Modelling the ice-flow changes
Based on the initial flow field, we reproduced the flow velocities observed during the three intervals, 1–3, of the two drainage events by considering the following two scenarios:
Scenario 1: The distribution of C presented in Figure 4b remained unchanged. We applied a load of lake water pressure corresponding to the highest lake water level reached in interval 1, then reduced the lake water load in drainage intervals 2 and 3, according to the decrease in lake water level measured during each period.
Scenario 2: The initial situation of interval 1 remained as in scenario 1, except that the lake water pressure was not considered at all. Flow-field solutions for the two subsequent intervals, 2 and 3, were obtained by readjusting C locally.
Note that in both scenarios the computed flow changes result from the sequence of snapshots of flow regimes which conform to different prescribed boundary conditions of lake water pressure or various C values. In the following we discuss the plausibility of the applied changes in the boundary conditions.
Scenarios 1 and 2 for the drainage in 2006
In scenario 1, when the lake water load is released, at markers 61 and 62 and, to a lesser extent, at marker 63, the ice-flow direction rotates towards Gornersee during 1 → 2 → 3 (Fig. 5a). Meanwhile, flow velocity magnitudes increase by ∼2–5% at markers 62 and 63. When comparing the computed flow velocities according to scenario 1 (Fig. 5a) with the observations (Fig. 2b) we find agreement on the clockwise flow azimuth rotation for markers 61 and 62. However, in the model, the increase of the flow velocities at markers 61 and 62 during 1 → 2 → 3 is much less pronounced.
(a) Computed flow velocities from scenario 1 for the 2006 lake drainage. The velocity of interval 1 includes the full lake water pressure, which is removed in two steps, of −6 and −18 m in intervals 2 and 3, respectively. (b) Computed flow velocities from scenario 2 for the lake drainage in 2006. (c, d) Relative changes applied to the initial C-configuration (Fig. 4b) to reproduce the velocities of intervals 2 and 3 from scenario 2, given in (b). The changes are expressed in the factor of increase or decrease of C. The greyscale referenced ‘+1’ and ‘+5’ indicates that C was set from zero (∼10−8) to 1 or 5.
In the results of scenario 2, the enhanced flow speeds at markers 61 and 62 in interval 3 can be better captured when changing the C-distribution of the bed-layer elements in the vicinity of Gornersee. The changes applied to C are depicted in Figure 5c and d, where the greyscale refers to the factor by which the value of C was increased or decreased relative to the initial configuration of Figure 4b. The increase of C (from 0 to 1 MPa−1 a−1 m) in the lake area for interval 2 (Fig. 5c) and the increase of C by a factor of 1.6 along the southwestern border of Gornersee and by 2.0 in the central confluence (Fig. 5d) gradually enhance the flow velocities at markers 61 and 62 and match the observed speed increase during the sequence 1 → 2 → 3. However, a change of flow direction, i.e. the rotation (clockwise) towards Gornersee, could hardly be introduced. We reduced C to half its value in the broad area down-glacier of the confluence area (Fig. 5c and d). The lowering of C relaxes the longitudinal coupling between the ice flow below the confluence and the ice-dam area and decreases the down-glacier component of the ice flow in the ice-dam area, hence the slight up-glacier rotation noticeable in 1 → 2 → 3. The decrease of C below the confluence results in lower velocities at markers g68 and g67 (intervals 2 and 3), which is inconsistent with the observations (Fig. 5b* cf. Fig. 2c). It was not possible to obtain a distribution of C where the modelled flow velocities were fully consistent in both the ice-dam area and down-glacier of the confluence, because the rotation of flow directions in the ice-dam area is less influenced by local enhancement/lowering of C. It was possible to modify the flow directions using the aforementioned (non-local) longitudinal-coupling effect, but this is barely supported by the observations.
Scenarios 1 and 2 for the drainage in 2004
The results from scenarios 1 and 2 corresponding to the drainage event of 2004 are given in Figure 6. In scenario 1, where we considered only the release of lake water pressure acting on the ice surface, the impact is minor at markers 43–45 and also at markers 41 and 42 (Fig. 6a). This is because the locations of markers 41 and 42 were ∼100 m to the west of their equivalents (61 and 62) in 2006. At this distance from the lake, the influence of the lake water pressure becomes small. The changes in flow speed and direction observed in 1 → 2 → 3 (Fig. 3b) cannot be reproduced in scenario 1. The markers are situated beyond the noticeable influence of the lake water pressure. This indicates that counterclockwise excursion of the flow direction from interval 1 to 2, observed at markers 41–45, does not relate to the influence of lake water pressure.
(a) Computed flow velocities from scenario 1 for the 2004 lake drainage. The velocity of interval 1 includes the full lake water pressure, which is removed in two steps, of −14 and −10 m in intervals 2 and 3, respectively. (b) Computed flow velocities from scenario 2 for the lake drainage in 2004. (c, d) Relative changes applied to the initial C-configuration (Fig. 4b) to reproduce the velocities of intervals 2 and 3 from scenario 2 given in (b). The changes are expressed in the factor of increase or decrease of the multiplier C. The greyscale referenced to ‘+4’ indicates that the multiplier C was set from zero (∼10−8) to 4.
In scenario 2 (Fig. 6b*) we captured the dominant speed increase of 1 → 2 in the modelled flow velocities, as featured at markers 41, 42, g44, 45, g46, g47 and 48. The changes applied to C for 1 → 2 are given in Figure 6c. On the southwestern border of Gornersee, C was increased by a factor of 1.2 and down-glacier of the ice dam by 2.5–3.5. The area of enhanced basal motion thus incorporates the southwestern border of Gornersee and extends northwest into the central confluence and northern part of Gornergletscher, then down-glacier to the lower tongue. Though we have only shown three marker locations situated along the centre flowline, other measurements confirm that the observed speed-up of 1 → 2 (Fig. 3c) takes place across the whole width of the lower tongue.
In interval 3 we restored C to its initial setting in the area down-glacier of the ice dam (Fig. 6d). Note that this would result in flow speeds at markers 41 and 42 similar to interval 1 (Fig. 6b*). To be consistent with the observations of interval 3, we must retain enhanced flow velocities while the iceflowdirection rotates up-glacier (clockwise) at these markers. This was achieved by further enhancement of C in the lake vicinity and along the inflow tributaries of Grenzgletscher and Gornergletscher (Fig. 6d).
Discussion
Mismatch of initial flow field in the ice-dam area
In the computation of the initial flow-field configuration, the largest discrepancies between modelled and measured flow velocities occurred in the ice-dam area. Everywhere else the agreement between measured and observed velocities was good. We observed that the orientation of flow velocities in the ice-dam area was significantly better matched when (1) C was considerably increased in the lower reaches of Gornergletscher and/or (2) the vertical resolution of the model was increased. The effect of (1) was an intensified longitudinal coupling between the ice-dam area and the lower tongue, so the ice was pulled down-glacier. We know that on Gornergletscher the ice viscosity, and thus the flow regime, is influenced by a massive core of cold ice which is advected down-glacier along the centre flowline of Grenzgletscher (Reference Eisen, Bauder, Lüthi, Riesen and FunkEisen and others, 2009). This was not considered in our model, but we speculate that cold rigid ice along the glacier central flow produces a drag on the ice-dam area. An increased vertical resolution (11-layer mesh) further reduced the mismatch between measured and computed flow velocities in the dam area but was not considered, due to reduced computational efficiency.
Scenario modelling
We considered two mechanisms by which the drainage of an ice-marginal lake affects the ice flow of the adjacent glacier: the load and release of lake water pressure arising from the filling and drainage; and enhanced basal motion caused by the subglacial drainage. We observed the simultaneous change in flow directions and increase of flow velocities at markers 61 and 62 during the 2006 lake drainage. In the modelling of these flow changes, when considering only the effect of decreasing lake water pressure, we were not able to capture the observed speed increase, whereas by modifying the C-distribution and thus influencing basal conditions, the observed rotation of the ice flow could hardly be reproduced. In Figure 7, we display a third scenario for 2006, in which the release of lake water pressure is combined with adjustments to C. In this case, we did not need to modify C below the confluence to reproduce the observed up-glacier (clockwise) rotation, as this event is now attributed to the impact of the lake water pressure. Aside from the two initial days of drainage, neither the borehole water-level measurements in bh6 and other boreholes (e.g. Reference Werder and FunkWerder and Funk, 2009), nor the monotonic flow velocities observed below the ice dam (Fig. 2c) indicate the spreading of a pronounced subglacial perturbation. The subglacial conditions were probably affected close to the region where the moulin shaft reached the glacier bed, which we estimate to be somewhere between bh6 and marker g66 (Fig. 1). In the combined scenario, the essential changes of C do correspond to this location, so we assume that enhanced basal motion only played a role local to the ice-dam area in the case of the 2006 lake drainage. We showed that the rotation of the ice flow due to decreasing lake water pressure yields the correct changes in flow direction and is a controlling mechanism affecting the ice flow in the ice-dam area. For the Gornersee drainage in 2006, the observed flow changes found in 1 → 2 → 3 are thus caused by the release of the lake water pressure acting on the ice and locally enhanced basal motion. Due to the slow, superficial drainage of 2006 and the absence of a pronounced subglacial perturbation, the effect of the release of lake water pressure could be identified.
(a) Computed flow velocities for a combined scenario in 2006. The velocity of interval 1 includes the full lake water pressure, which is removed in two steps, of −6 and −18 m in intervals 2 and 3, respectively. (b, c) Relative changes applied to the initial C-configuration (Fig. 4b) to reproduce the velocities of intervals 2 and 3. The greyscale referenced to ‘+1’ and ‘+5’ indicates that C was set from zero (∼10 − 8) to 1 or 5.
In interval 2 of 2004 the borehole water level in bh430 was high and then decreased monotonically between 5 and 7 July, dropping sharply on 7 July. Also, in the ice dam basal water pressure was high before 6 July (Reference Sugiyama, Bauder, Huss, Riesen and FunkSugiyama and others, 2008, fig. 7). In interval 2 the lake water presumably overcharged the subglacial drainage system of Gornergletscher and is certainly responsible for disturbed subglacial conditions, inducing broadly enhanced basal motion. For 1 → 2 the observations of strongly increased flow velocities in both the ice-dam area and the lower tongue can be corroborated with the plausible modelling results (Fig. 6b and c). We assume the lake water influenced the ice flow during the drainage of 2004 in a similar manner to that in 2006, as the maximum lake water level was similar at the two drainage events. Nevertheless, the results indicate that markers 41 and 42 were less influenced by the lake water due to the increased distance from Gornersee and also because the impact of enhanced basal motion was overwhelming. The reversal in flow direction with an almost complete switch in flow direction from west to (south)east was observed in 2 → 3 (Fig. 3b), i.e. at g44 and 45. This feature is absent in the model results. Also, we do not think the further enhancements applied to C in Figure 6d provide a plausible interpretation for 2 → 3.
In the current model, we could not find an appropriate configuration of the constant, C, to reproduce these changes. The linear relation in the basal boundary condition may be insufficient. Errors in the bed topography may also play a role. Referring to the same part of the ice dam and the drainage in 2004, Reference Sugiyama, Bauder, Weiss and FunkSugiyama and others (2007b) mentioned that changes in longitudinal stress can produce reversal of the ice motion, but changes in basal conditions (basal motion and separation, sediment deformation) did not provide a sufficient explanation. Sugiyama and others also mention that the reverse motion correlates with the increase or decrease in lake discharge, which is closely connected to the lake water level. We conjecture that the flow change 2 → 3 with the unexplained observations at markers 43–45 may be influenced by further factors, due to the effect of lake water pressure, not captured in the current model. When the lake fills to its maximum possible size as in both 2006 and 2004, the lake water pressure can perturb the ice flow at a distance of up to 300 m from Gornersee. The lake emptied significantly faster in 2004 than in 2006. In 2004, the lake drained in <5 days, whereas in 2006 the drainage took place over 20 days. The rate of decrease of the lake water level, the change of stress conditions in the ice dam and fracturing of the ice further influence the ice in the lake vicinity. It may be worthwhile to analyse the ice flow as part of a transient study considering a time-dependent, possibly elastic, response of the ice to explain the presence (absence) of the reversal in flow direction at markers 43–45 (63–65) in 2004 (2006) with variable loading/unloading test cases.
Conclusion
We studied ice-flow changes on Gornergletscher during two drainage events of Gornersee in 2004 and 2006. We used a 3-D numerical flow model and produced scenarios of ice-flow changes through adjusting the boundary conditions of spatially non-uniform basal motion and loading/unloading of lake water on the ice covered by the lake.
For each drainage event, we compared the results of the numerical experiments to our observations and identified the major mechanism generating the flow observations. The slow lake drainage of 2006 had a minor, local impact on the basal conditions, and the ice-flow changes were dominated by the release of the lake water load. In 2004, the draining lake water rapidly entered the subglacial environment and perturbed the basal conditions, locally enhancing basal motion of Gornergletscher.
The lake water pressure acting on the ice surface of the adjacent ice dam has not, until now, been considered to have a relevant influence on the ice flow during a glacial lake drainage event. We demonstrated that changes in lake water-pressure load do have a notable and varying effect on the ice flow. The lake water pressure also affects the local stress conditions in the ice dam. Further experiments investigating the response of the ice dam to variable time-dependent loads and enhanced constitutive properties of the ice (i.e. elasticity, damage) may provide further information on the role of the ice dam.
Acknowledgements
This research was funded by Swiss National Science Foundation grants 200021-103882/1 and 200020-111892/1. We are grateful to the International Foundation of High Altitude Research at Gornergrat (HFSJG) for providing facilities and support. We thank M. Huss, M. Lüthi and an anonymous reviewer for comments which considerably improved the manuscript. Thanks are extended to S. Braun-Clarke for proofreading and H.A. Fricker (Scientific Editor).
