Skip to main content Accessibility help


  • Access
  • Cited by 3


      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Numerical simulation of three-dimensional velocity fields in pressurized and non-pressurized Nye channels
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Numerical simulation of three-dimensional velocity fields in pressurized and non-pressurized Nye channels
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Numerical simulation of three-dimensional velocity fields in pressurized and non-pressurized Nye channels
        Available formats
Export citation


Channels incised into bedrock, or Nye channels, often form an important component of subglacial drainage at temperate glaciers, and their structure exerts control over patterns and rates of (a) channel erosion, (b) water flow-velocity and (c) water pressure. The latter, in turn, exerts a strong control over basal traction and, thus, ice dynamics. In order to investigate these controls, it is necessary to quantify detailed flow processes in subglacial Nye channels. However, it is effectively impossible to acquire such measurements from fully pressurized, subglacial channels. To solve this problem, we here apply a three-dimensional, finite-volume solution of the Reynolds averaged Navier– Stokes (RANS) equations with a one-equation mixing-length turbulence closure to simulate flow in a 3 m long section of an active Nye channel located in the immediate foreground of Glacier de Tsanfleuron, Switzerland. Numerical model output permits high-resolution visualization of water flow through the channel reach, and enables evaluation of the experimental manipulation of the pressure field adopted across the overlying ice lid. This yields an increased theoretical understanding of the hydraulic behaviour of Nye channels, and, in the future, of their effect on glacier drainage, geomorphology and ice dynamics.


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(Nye,1973)incised in to the substrate: so-called linked-cavity drainage systems(Walder, 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., Lane and others, 1999; Ma 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 (Hubbard 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).


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 km2 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–1) 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 (Sharp and others, 1989; Hubbard 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 (Hubbard 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., Chanson, 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 (Stoesser and Bates, 2002; Stoesser 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


and the three momentum equations written generally as


where Ui (i = u, v, w) is the velocity averaged over the time t, xi (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 (Rodi, 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 (Smagorinsky, 1963), was selected. This approach is based on Prandtl’s mixing-length model (Prandtl, 1925) and is most commonly used as a sub-grid-scale closure in large eddy simulations (Rodi 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., Hervouet 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., Chanson, 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.


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.


Funding for this project was provided by U.K. Natural Environment Research Council grant NER/A/S/2000/01144 to P.D.B. and M.J.S. B.P.H. acknowledges receipt of a grant from the University of Wales at Aberystwyth research fund which was used to support this work.


Chanson, distri, H. 1999. The hydraulics of open channel flow: an introduction. London, Arnold.
Hervouet, J.-M. and van Haren, L. 1996. Recent advances in numerical methods for fluid flows. In Anderson, M. G., Walling, D.E. and Bates, P. D., eds. Floodplain processes. Chichester, John Wiley and Sons, 183–214.
Hubbard, B. and Hubbard, A. 1998. Bedrock surface roughness and the distribution of subglacially precipitated carbonate deposits: implications for formation at Glacier de Tsanfleuron, Switzerland. Earth Surf. Processes Landforms,23(3), 261–270.
Hubbard, B., Siegert, M. J. and Mccarroll, D. 2000. Spectral roughness of glaciated bedrock geomorphic surfaces: implications for glacier sliding. J. Geophys. Res.,105(B9), 21, 295–21, 303.
Lane, S. N., Bradbrook, K. F., Richards, K. S., Biron, P.A. and Roy, A. G. 1999. The application of computational fluid dynamics to natural river channels: three-dimensional versus two-dimensional approaches. Geomorpholog y,29(1–2), 1–20.
Ma, L., Ashworth, P. J., Best, J. L., Elliot, L., Ingham, D. B. and Whit-combe, L. J. 2002. Computational fluid dynamics and the physical modelling of an upland urban river. Geomorpholog y,44(3–4), 375–391.
Nye, J. F. 1973. The motion of ice past obstacles. In Whalley, E., Jones, S.J. and Gold, L.W., eds. Physics and chemistry of ice. Ottawa, Ont., Royal Society of Canada, 387–394.
Prandtl, L. 1925. Bericht uber Untersuchungenzur ausgebildeten Turbulenz. Z. Angew. Math. Mech., 5, 136–139.
Rodi, W. 1980. Turbulence models and their application in hydraulics. Delft, International Association for Hydraulic Research. (IAHR Special Publication.)
Rodi, W., Ferziger, J. H., Breuer, M. and Pourquie, M. 1997. Status of large eddy simulation: results of a workshop. J. Fluids Eng.,119(2), 248–262.
Sharp, M. J., Gemmell, J. C. and Tison, J.-L. 1989. Structure and stability of the former subglacial drainage system of Glacier de Transfleuron, Switzerland. Earth Surf. Processes Landforms, 14, 119–134.
Smagorinsky, J. 1963. General circulation experiments with the primitive equations. Mon. Weather Rev., 91, 99–164.
Stoesser, T. and Bates, P. D. 2002. Compound channel now using an implicit time scheme and a hybrid turbulence model. In Falconer, R. A., Lin, B., Harris, E. L. and Wilson, C. A. M. E., eds. Hydroinformatics 2002. Proceedings of the Fifth International Conference on Hydroinformatics, 1–5 July 2002, Cardiff, U. K. London, IWA Publishing, 45–51.
Stoesser, T., Wilson, C. A.M. E., Bates, P. D. and Dittrich, A. In press. Development and application of a 3-D numerical tool for the design of an integrated hydraulic and ecological sensitive river scheme.J. Hydroinformatics.
Walder, J.S. 1986. Hydraulics of subglacial cavities.J. Glaciol.,32(112), 439–445.