Introduction

Understanding the response of glaciers to global climate change represents a key problem for society. Solving this problem requires a full appreciation of the processes operating beneath glaciers, particularly of spatial and temporal variations in subglacial water pressure. However, our ability to record such information is seriously restricted by the inaccessibility of subglacial environments. At bedrock-based glaciers, subglacial water is commonly routed via cavities linked by Nye channels(Reference Nye, Whalley, Jones and GoldNye,1973)incised in to the substrate: so-called linked-cavity drainage systems(Reference WalderWalder, 1986). Previous calculations of water flow through linked cavities have been limited in two ways. First, they have failed to account for realistic channel geometry and roughness, despite these features controlling the turbulence and velocity of waterflow within those channels. Second, calculationsofwater pressure have, up until now, been based on theoretical considerations of spatially averaged ice pressure but neglecting the hydraulics of the actual basal flow pathways involved. This yields general pressure fields and makes a detailed analysis of subglacial hydrodynamicsacross entire drainage networks impractical.Recent advances influid dynamic models (see, e.g., Reference Lane, Bradbrook, Richards, Biron and RoyLane and others, 1999; Reference Ma, Ashworth, Best, Elliot, Ingham and Whit-combeMa and others, 2002) make it now possible to calculate the hydrodynamics of these systems, using accurate channel boundary conditions. Such modelling can provide information on water velocity, pressure and erosion, and may facilitate further work on the acquisition of solutes by subglacial water. Such information is fundamental to investigations of glacier sliding, the interpretation of glacial hydrological field measurements, and to calculations of fluvial erosion in such environments.

The immediate foreground of Glacier de Tsanfleuron, Switzerland, is an ideal field site for measuring water flow in Nye channels (see Fig. 1a nd discussionbelow). By combining knowledge of the field site (Reference Hubbard, Siegert and MccarrollHubbard and others, 2000) with new measurements of the recently exposed subglacial geomorphology (including sub-mm-scale roughness) and numerical modelling of subglacial fluid dynamics, a full representation of water flow through a Nye channel has been developed for the first time.

Fig. 1. (a) Map showing the margin of Glacier de Tsanfleuron, and the area of the recently exposed limestone glacier forefield. (b) Photograph of the Nye channel measured by an EDM showing the propeller current meter suspended by a tripod (taken from the ice margin looking across the limestone platform).

Objectives

In this paper, we aim to establish a fully quantitative understanding of the hydrodynamics of a short reach of subglacial Nye channel. The objectives of the project were as follows:

• (1) To obtain an accurate digital elevation model (DEM) of a real channelized subglacial environment in the fore-field of Glacier de Tsanfleuron.

• (2) To determine the velocity and turbulence of water through this system by routing specific fluxes of melt-water to channels at the head of the DEM area.

• (3) To calibrate a computational fluid dynamics (CFD) model by reconstructing these measurements of water flow through channels.

• (4) To make simple adjustments to this calibrated CFD model to account for the effect of overlying ice, and determine from these water fluxes velocities and pressures within a linked-cavity environment.

Field Data Collection

Fieldwork was carried out between 24 August and 4 September 2001 on the immediate forefield of Glacier de Tsanfleuron, located northwest of Sion in Switzerland. The glacier has a surface area of ~4 km₂ and, based on recent global positioning system (GPS) survey and DEM data, extends from -2420 to ~2850ma.s.l. The glacier’s retreat since its Little Ice Age position, which has been rapid (-5 ma–₁) in recent years, has resulted in the exposure of the former subglacial bed with an extraordinarily fresh subglacial morphology (Fig. 1). Nye channels, linked cavities and striations are evident throughout the platform (Reference Sharp, Gemmell and TisonSharp and others, 1989; Reference Hubbard, Siegert and MccarrollHubbard and others, 2000). A 3 m long section of a Nye channel located at the glacier margin was selected for investigation (Fig. 1). This location benefited from being recently exposed and therefore relatively unaltered by subaerial erosion (Reference Hubbard and HubbardHubbard and Hubbard, 1998). It was also proximal to a melt-water source, which could easily be routed into the channel.

An electro-optical distance meter (EDM) was used to generate a 3 x 1 m D EM of a section of proglacial Nye channel at a grid spacing of < 5 cm and with a vertical resolution of 1mm. Within this domain, the Nye channel was approximately 0.5 m wide and 0.5 m deep and near-rectangular in shape. This resulted in a raw topographic dataset of approximately 4000 points, of which approximately half were located below the surveyed water level. Water was routed into the former subglacial bedrock from the nearby ice margin. The water velocity through the channel within the study area was then measured with a propeller current meter using a 5 min integration period at each measurement point. A long integration period is needed to produce a true mean velocity value where measurement noise is minimized, and the result is a value for mean downstream velocity in the x and y Cartesian directions (vertical velocities are not readily determined with this technique). Finally, bed roughness was determined according to Strickler’s grain-size law (see, e.g., Reference Chanson, distriChanson, 1999, p. 84), as the channel substrate consisted of relatively large and angular gravel. Velocity measurements were taken at the head of the reach to establish both the volume flux into the section and the inlet velocity profile, and in the middle of the domain in order to validate the model independently. The water surface profile through the reach was also surveyed using the EDM. Over the 3 day study period, water surface elevations within the channel varied by only 2–3 cm.

Numerical Modelling

The field data provide the necessary topography, boundary condition, parameterization and validation information to fully configure a three-dimensional finite-volume CFD model previously developed at the University of Bristol (Reference Stoesser, Bates, Falconer, Lin, Harris and WilsonStoesser and Bates, 2002; Reference Stoesser, Wilson, Bates and DittrichStoesser and others, in press). The model calculates hydrodynamics for a general three-dimensional geometry, discretized by the finite-volume method on curvilinear coordinates. The Reynolds averaged Navier–Stokes (RANS) equations are solved with the continuity equation written as

(1)

and the three momentum equations written generally as

(2)

where U_i (i = u, v, w) is the velocity averaged over the time t, x_i (i = x, y, z) is the spatial geometrical scale, p is the water density, P is the pressure, δ is the Kronecker ij delta, u' are the instantaneous velocity fluctuations in time during the time-step dt, when U is subtracted, and the subscripts i and j refer to the two velocity components relevant to each momentum equation. The second term on the right side of the equation is the Reynolds stress term, which needs to be solved by a turbulence model. Initially, a two-equation k-ε (model (Reference RodiRodi, 1980) was tried. As expected, however, this proved unstable because the topographic complexity of the domain induced very steep gradients in turbulent kinetic energy (k). These gradients were sufficient to prevent convergence of the iterative scheme used in the numerical model. Instead, a simpler one-equation turbulence model, the Smagorinsky scheme (Reference SmagorinskySmagorinsky, 1963), was selected. This approach is based on Prandtl’s mixing-length model (Reference PrandtlPrandtl, 1925) and is most commonly used as a sub-grid-scale closure in large eddy simulations (Reference Rodi, Ferziger, Breuer and PourquieRodi and others, 1977), although it can also be used in RANS modelling as employed here. Lastly, the SIMPLEC methods (semi-implicit method for pressure-linked equations) were used for pressure velocity corrections, and the water level was specified as a rigid no-stress lid. At the centre of each finite-volume cell, the model therefore predicts pressure and the three components of the mean velocity vector U = (u, v, w) averaged over the time-step in the x, y and z Cartesian directions, respectively.

The model discretization was constructed using dedicated mesh generation software that constructs a smoothly varying planform discretization of the domain as a curvilinear grid. The vertical resolution of the mesh was then discretized within the code as a sigma transformation (see, e.g., Reference Hervouet, van Haren, Anderson, Walling and BatesHervouet and Van Haren, 1996) to give a normalized thickness of vertical layers irrespective of bottom topographic variations. In order to capture correctly the expected complex flow patterns, a relatively high mesh resolution was constructed for the 360.5m channel with 74625620 computationalpoints in the down-channel(x), cross-channel (y) and vertical (z) directions, respectively. For each point on the lowermost plane of the grid, a topographic height was assigned based on weighted nearest-neighbour interpolation of the four closest points in the high-resolution topographic survey. This gave a total of 33288 computational cells with approximate size 0.0460.0260.03 m (see Fig. 2). Figure 2 clearly shows the topographic complexity of the domain, even over such a short reach.

Fig. 2. Finite-volume discretization construct ed for a proglacial Nye channel shown in (a)two-dimensional planform and (b) full three-dimensional views. Note in the three-dimensional view the high topographic complexity of the study reach.

For the free-surface simulations, boundary conditions for the model consisted of the upstream inflow discharge and velocity distribution determined from the measured data and the downstream water surface elevation. Given the small scale ofthe reach and the near-absence of dynamic changes in depth, the surface boundary was discretized as a rigid no-stress lid according to the surveyed water surface profile. Finally, channel side and bed boundaries were treated as impermeable walls with the no-slip condition applied. For the bottom boundary, a Manning friction coefficient (see, e.g., Reference Chanson, distriChanson, 1999, p. 82–83) was used to represent energy losses due to bed roughness, and, given the relatively large clast size observed in the field, a value of 0.03 was specified for the whole domain.

The above boundary conditions were then used to drive steady-state solutions of Equations (1) and (2). Initial conditions for each simulation were zero velocity within the domain and hydrostatic pressure. Convergence was deemed to have occurred when all variables showed variations of 50.01% between iterations, which was achieved in approximately 500 iterations.

Once the free-surface model had been constructed and run successfully, its boundary conditions were changed to represent pressurized flow within the channelas if it were still subglacial. The ice lid was assumed to be at the same position as the free-water surface and was treated analogous to the channel wall boundaries. Simulations were then rerun, keeping all other aspects of the model structure constant.

Model Results

Figure 3 shows a comparison between measured and simulated downstream velocity profiles over the flow depth as recorded in the centre of the model domain (point A in Fig. 2). Specifically we compare the product of the time-averaged velocity vector components in the x and y Cartesian directions with the mean velocity as determined by a propeller current meter. Both modelled and measured values of velocity represent non-fluctuating mean fields in time and are integrated over approximately the same spatial volume (cell size of 0.0460.0260.03m in the model compared to a propeller diameter of ∼0.03m) and hence are broadly comparable. The section-average velocity measured using dye tracing also matched well with the corresponding average of the predicted velocity field. Thus,it is clear from Figure 3 that the model provides an acceptable match to the observed data and can therefore be used to make initial process inferences concerning the fluid dynamics of water in pressurized and non-pressurized Nye channels.

Fig. 3. Comparison of measured and predicted vertical velocity profiles (the product of mean velocities u and v in the x and y Cartesian directions) taken in the centre of the model domain (point A in Fig. 2a).

Visualization of the full three-dimensional velocity field produced by each model is inherently difficult. Therefore, Figure 4 shows the three-dimensional velocity vectors predictedby the model a long a two-dimensional slice in the horizontal plane through the domain at z = 0.80 m, or approximately 12 cm below the water surface. This represents the mean velocity vector U averagedover the time-step dt. Figure 4 shows the model to be capable of predicting complex recirculations present within the domain. These are manifest by strong vortices with vertically oriented axes. In the full three-dimensional velocity vector plot, it is also clear that there are regions of strong upwelling and down-welling within the flow. The flow is thus highly turbulent and complex, even over this relatively short section.

Fig. 4. Full three-dimensional time-averaged vector field U obtained from the free-surface model. Note the complex recirculations present.

Figures 5 and 6 show comparisons of u (the component of the mean velocity vector U in the x Cartesian direction averaged over time-step dt) and mean pressure P for the free-surface and ice-covered models. These diagrams show that imposition of an ice lid leads, as one would expect, to a change in both the pattern and magnitude of flow quantities. Given the volumes of data produced by the model, it is only possible to provide a snapshot of the total predicted field. However, the differences evident from comparison of Figures 5 and 6 are typical. In addition, the model is able to predict values for such quantities as fluid total viscosity (the sum of molecular and turbulent viscosity) and turbulent kinetic energy. Clearly, this gives an extremely powerful tool with which to commence the investigation of fluid dynamic processes within Nye channels.

Fig. 5. Comparison of predicted downstream time-averaged velocity distribution (vector component u in the x direction with units m s^–1) for a vertical slice taken at approximately x = 2.4 m through the (a) free-surface and (b) ice-covered flow models. Note the different scales.

Fig. 6. Comparison of predicted time-averaged pressure butions (in N) for a vertical slice taken at approximately x = 2.4 m through the (a) free-surface and (b) ice-covered flow models. Note the different scales.

In both pressurized and free-surface cases, the water flow appears complex. By comparison, however, water flow is more ordered in the “ice-covered”experiment (Fig. 5 and 6). In other words, gradients in water pressure and velocity orthogona lto the general line of water flow are lower in the ice-covered situation (Fig. 5 and 6).

The results are testament to the utility of our method in future glaciological research. Simulations were run on a 450 MHz Pentium 3 with 256 Mb of RAM and took 5 5 min to complete. This suggests that quite substantial sections of Nye channel, say several tens of metres, can be simulated using such an approach. Given the ability of the model to cope with complex topography, such as cavities and tributary junctions, it is likely that in the future we may be able to gain considerable insight into subglacial water, sediment and chemical fluxes in proglacial and subglacial channel networks using the approach detailed in this paper.

Conclusions

A 3 m long section of a recently exposed Nye channel, in the immediate forefield of Glacier de Tsanfleuron, was surveyed using an electro-optical distance meter, from which a 3×1m DEM was created. Water was routed into the channel, and the flow measured using a propeller current meter and dye tracing. The DEM was then used to create a 33288 cell mesh used as a topographic boundary condition for a CFD model, which was run to steady state.

The CFD model’s prediction of free-surface water flow through the Nye channel matched well with measurements from flow-meter and dye-tracing experiments. This match between model results and measurements allows us to have confidence that the model is appropriate, and that it produces useful process information with respect to water flow through Nye channels in both free-surface and pressurized-lid (subglacial) environments.

The results display complex circulation patterns within the water flow, including strong vortices with vertical axes. Upwelling and downwelling of water is predicted across the 3 m section of the model domain, signalling that the flow is highly complex and turbulent over this short distance. The pattern of water flow is changed markedly when a pressurized (fully subglacial) channel is modelled. In general, pressurized-lid simulation showed less complicated water flow than in the free-surface case. In both experiments, however, the water flow was found to be intricate. Hence, the application of CFD modelling to water flow through real Nye channels has allowed us, for the first time, to quantify hydraulic properties precisely and at a high spatial resolution. To date, such hydraulic data have largely been unavailable or assumed. CFD model output includes water-flow velocity, turbulence and the response of both to variations in input—critically, under both atmospheric and pressurized conditions. While the former has the capacity to provide insights into subglacial flow near the frontal margins of glaciers and in proglacial areas, the latter can allow the reconstruction of subglacial water flow beneath thicker ice. The quantitative reconstruction of such flow properties will have a direct significance for water transfer beneath hard-bedded glaciers and the relationship between that water transfer and motion of the overlying ice.

These results show how water flow through Nye channels can be predicted. Each model simulation took approximately 5 min to perform. Clearly, the technique demonstrated here can, in future, be used to model longer and multiple interconnected Nye channels. In this way, the water flow through large sections of a hard-bedded subglacial environment can be modelled. In addition, the CFD model approach allows for the prediction of turbulent kinetic energies and water temperatures. These predictions, in conjunction with the results presented here, may in future be used to understand Nye-channel erosion rates and solute acquisition and transport within sub- and proglacial channelized systems.