Introduction
Models of ice-mass motion are constrained by, amongst other things, imperfect knowledge of spatial variations in ice viscosity. Although the rheological effects oftemperature variations and (relatively soft) Pleistocene-age basal ice are adequately catered for in large, polythermal ice masses (see, e.g., Reference PaynePayne and others, 2000), spatial variations in rheologically important properties such as water content and crystal orientation remain unknown, particularly in temperate glaciers. Current models of temperate glacier flow are therefore forced to assume ice viscosity is spatially uniform, potentially disregarding important spatial variations in the response of the ice to applied stresses.
Here, we reconstruct the liquid-water content of Glacier de Tsanfleuron, Switzerland, from the bulk ionic composition of eight ice cores recovered from a flowline transect along the glacier. These water-content data are used to determine the relative viscosities of the ice, following Reference DuvalDuval (1977), and the resulting layered rheology is used to constrain a two-dimensional model of ice flow at the glacier. This model and one incorporating rheologically homogeneous ice are tuned by the current glacier’s geometry and run under advance and retreat scenarios. Finally, outputs from the two models are compared and the potential importance of rheological variations for predicting the geometrical response of temperate glaciers, such as Glacier de Tsanfleuron, to climate change is evaluated.
Field Site and Methods
Glacier de Tsanfleuron has a surface area of ∽3.5 km2, and flows over macro-porous Cretaceous and Tertiary limestone from ∽2950 to ∽2450ma.s.l. (Fig. 1). The glacier has been widely studied, particularly in terms of its basal ice layers (Reference Souchez and LemmensTison and Lorrain, 1987; Hubbard and Sharp, 1995) and subglacialcarbonate crusts exposed on its proglacial bedrock surface by retreat over the past ∼150 years (Reference Hallet, Lorrain and SouchezHallet and others, 1978; Reference Lemmens, Lorrain and HarenLemmens and others, 1983; Reference PayneSouchez and Lemmens, 1985; Reference PatersonSharp and others, 1990; Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and Reference Hubbard and HubbardHubbard, 1998; Reference Hubbard, Siegert and McCarrollHubbard and others, 2000b).
Glacier de Tsanfleuron, Switzerland, with core locations marked as closed circles and labelled by the last two digits of the year drilled, followed by a hyphen and the order drilled.
Ice-facies character and sedimentology
Reference Hubbard and SharpHubbard and Sharp (1995) sampled ice at marginal exposures and in subglacial cavities at Glacier de Tsanfleuron, identifying three principal ice facies, termed englacial, clear and dispersed (Table 1). Englacial ice comprises the bulk of the glacier’s volume and is formed by firnification near the glacier surface, and is commonly referred to as “glacier ice”. Three samples of this ice yielded a mean debris content of 1.90 6 102 g m–3. Clear- and dispersed-facies ice both acquire their distinctive character as a result of processes operating near the glacier bed. Clear-facies ice occurs as a massive layer, some metres thick, that contains dispersed debris clots towards its base, and from which virtually all gas bubbles have been expelled. The debris present towards the base of the facies is considered to have been introduced by deformation-related mixing with the underlying debris-rich (dispersed-facies) ice. Dispersed-facies ice is universally present as a layer -1.5 m thick immediately above the ice–bed contact. Analysis of 26 debris samples from this facies indicated that it contains relatively coarse debris that is dispersed throughout the ice at a mean concentration of 4.65 6 104 g m–3 (Table 1).
Summary of Glacier de Tsanfleuron ice properties (ice facies and column abbreviations are explained in the text)
Ice cores
Ice classification
We have recovered eight vertical ice cores from along a flow-line at Glacier de Tsanfleuron (Fig. 1). The ice forming these cores was analyzed for its ionic composition, its isotopic composition, its gas content and composition, its debris content and texture, and its crystal size and orientation (Reference Hubbard, Tison, Janssens and SpiroHubbard and others, 2000a; Reference Sharp, Tison and FierensTison and Hubbard, 2000). Visual logging at a vertical resolution of 10 mm allowed all of the ice to be categorized into one of three zones: an upper zone (UZ), a lower zone (LZ) or a basal zone (BZ). The UZ is formed predominantly of ice containing high concentrations of fine (5 1 mm diameter) bubbles, and closely corresponds to the englacial-facies ice described above. The L Z is defined by a complete absence of any fine bubbles, and is formed predominantly of ice containing either no bubbles or widely dispersed large bubbles. This zone equates directly with the clear facies described above. The BZ is formed exclusively of bubble-free and debris-rich ice, corresponding to the dispersed facies described above.
View of the long section and its constituent ice zones (UZ, upper zone; LZ, lower zone; BZ, basal zone) used for two-dimensional modelling. The properties of the ice zones are explained in the text and in Table 1.
Major-ion concentration
Two hundred and fifty-three ice samples recovered from the three zones as represented in all the ice cores were analyzed for Na+ , K+ , Ca2 + , Mg 2 + , NO3 2–, C l – and SO4 2– concentrations by ion chromatography (Dionex DX-100). Laboratory procedures are described elsewhere (Reference Hubbard and HubbardHubbard and others, 2000a), yielding an analytical precision that was typically ± 5 % at concentrations 4 5 0 μeq L–1 , ± 3 0% at concentrations <10 μeq L–1 and on the order of ±100% at concentrations of ~1 μeq L–1. HCO3 – concentrations were determined by charge balance to an accuracy of ± 1 0% (Reference FairchildFairchild and others, 1999; Reference Hubbard and HubbardHubbard and others, 2000a).
Calculation of Ice Liquid-Water Content
Most of the liquid water present at the scale of individual crystals within temperate glaciers exists within veins present at triple-grain junctions (Reference Nye and C.|FrankNye and Frank, 1973; Reference Glen, Homer and ParenGlen and others, 1977; Wolff and Paren, 1984). Following Reference NyeNye (1991) and Reference MaderMader (1992), the bulk vein-water concentration in ice (the volume of water per unit volume of ice) W can be approximated from its bulk ionic concentration Mb and the temperature depression 9 of the ice below its pressure-melting point:
where the constant K (in K m 3 g–1) is given by the slope of a
where the values are taken relative to the englacial ice of the upper zone, here denoted by the subscript “UZ”.
Correction for solute acquisition during sampling
Debris-bearing ice was sampled for chemical analysis by melting the solid sample and allowing the resulting melt-water to drain through a filter. This process necessarily involves the meltwater contacting the released debris, typically for a period of ~1 hour in the present study (where core samples were melted in the laboratory rather than in the field). Measured bulk solute concentration Mm therefore includes solute that has been acquired from the dissolution of incorporated debris during sampling, Ma. Thus, the original bulk ionic concentration of the ice prior to sampling, Mb , is given by
In order to calculate the concentration of solute acquired during sampling, Ma, we apply acquisition rates derived from laboratory dissolution experiments to our laboratory sampling conditions. Following Reference LermanLerman (1979), the concentration C of a solute in solution is expressed as function of time t as
where Cs (g cm–3) is the ultimate (steady-state) concentration of the material in solution, C0 (g cm–3) is the concentration of solute released by rapid, surface-exchange reactions during the early stages (5 180 s) of water–rock contact, and the constant k is the reaction rate parameter, dictated by the precise solvent–mineral combination involved. Reference Brown and HubbardBrown and others (2001) carried out a suite of low-temperature laboratory dissolution experiments to calculate the values of C0, Cs and k for contact between sediments and glacial melt-waters. Results indicated that C0 and Cs are strongly related to the specific surface area SSA (the total surface area of sediment (m2) per unit volume of water (m3); m–1) of the sediment concerned, yielding, for Ca2 + ,
*The glacier’s actual long-section area, measured in 2001, is 116780 m2.
and
Here, we use these authors’ acquisition rates for Ca2 + because measured bulk solute at Glacier de Tsanfleuron is dominated by the dissolution of CaCO3 (yielding the principal ions Ca2 + and HCO3 2–). We double the predicted Ca2 + concentration to account for the balancing HCO3 – .
We determine SSA from the grain-size distribution of the sediment entrained within the ice by calculating the number of particles within each measured size range N and the mean (spherical equivalent) diameter D of those grains following the methods of Reference Hooke, B. and IversonHooke and Iverson (1995). Specific surface area by mass, SSAm (m2 g–1), which is a material property of the sediment, is then calculated as the number of particles of each diameter per unit mass of sediment multiplied by the surface area of each particle. The resulting SSAm is multiplied by the measured concentration of that sediment in the ice facies (gm– 3) to yield the SSA (m–1). Results of this process yielded SSAm values of 1.58 m2 g–1 for the englacial-facies debris, and 0.555 m2 g – 1 for the dispersed-facies debris. Observations within basal cavities indicate that the concentration of sediment incorporated in clear-facies ice is generally similar to that of the overlying englacial facies, except within the basal ~1m which contains sediment at a concentration that is intermediate between that above and that of the underlying dispersed-facies ice. We therefore calculate the net sediment concentration of a 15 m thick clear-facies/LZ ice layer by assuming that: (a) basal sediment has been introduced into the uppermost 14 m of the layer at a concentration equal to that of the existing englacial sediment (yielding a concentration of 3.8 6 102 gm–3), and (b) the basal 1m of the layer contains a sediment concentration that decreases linearly from 50% of that of the underlying dispersed-facies ice at its base to that of the overlying clear-facies ice at its top, yielding a concentration, in this basal 1 m of clear-facies ice, of 1.20 6 104 gm–3 . These calculations yield an overall sediment concentration for the 15 m thick clear-facies ice layer of 7.67 6 102 gm–3 . The corresponding SSAm of this clear-facies debris (scaled such that the debris located within the basal 1 m is similar to that of the underlying dispersed facies and the rest of the debris is similar to that of the overlying englacial facies) is 0.613 m2 g–1. Resulting values of SSA are 3.00 6 102 m–1 for the englacial facies, 4.70 6 102 m–1 for the clear facies, and 2.62 6 104 m–1 for the dispersed facies (Table 1). Substituting these values into Equations (5) and (6) yields values of C0 and Cs of 0.14 and 0.44 gm–3 for the englacial facies, 0.22 and 0.69 gm–3 for the clear facies, and 12.2 and 38.4 g m–3 for the dispersed facies, respectively.
While sampling typically took ∽1 hour, for most of that time the sample was only partially melted, allowing correspondingly partial rock–water contact. We therefore adopt a value of t of 20 min (1200 s). Substituting this, along with the values for C0 and Cs derived above, into Equation0 (4) yields values of the solute concentration acquired during sampling C of 0.50 g m–3 for the englacial facies, 0.79 g m– 3 for the clear facies and 43.94 g m–3 for the dispersed facies. Substituting these values of M a into Equation (3) yields values of the bulk ionic concentration in the initial ice sampled (Mb) of 0.68 g m–3 for the englacial facies, 1.23 g m–3 for the basal clear facies and 7.34 g m–3 for the basal dispersed facies (Table 1). Substituting these values of Mb into Equation (2) gives W' = 1.80 in the clear-facies ice, and W' = 10.74 in the dispersed-facies ice (Table 1).
Modelling
Inmostglacierflowmodels, the ice softness value A (a–1 bar–3 ) in Reference GlenGlen’s (1955) flow law is held at a spatially uniform value that is tuned to the ice-thickness profile or velocity field of the ice mass concerned. Since our data indicate that the water content of Glacier de Tsanfleuron varies systematically between different ice zones, we constrain the hardness of the LZ ice (ALZ) and the BZ ice (ABZ) relative to that of the UZ ice (AUZ) in addition to empirical tuning. Reference DuvalDuval (1977) reported a strong linear correlation between measured water content W and effective strain rate, yielding:
Thus, the relative hardnesses of the clear facies (ALZ) and dispersed facies (ABZ) scale linearly with reconstructed values of W' to become 1.80 AUZ and 10.74AUZ, respectively.
We associate these relative viscosities with an approximation of the spatial extent of their host facies as revealed in our ice cores, and substitute the resulting layers into a two-dimensional (flowline) model of ice flow at Glacier de Tsanfleuron. To reconstruct these layers from our ice-core data, we consider the UZ to be composed entirely of englacial ice, the LZ entirely of clear-facies ice, and the BZ entirely of dispersed-facies ice. The model, which is based on a finite-difference first-order solution of the ice-flow equations (Reference BlatterBlatter, 1995), and therefore includes the effects of longitudinal or normal deviatoric stresses, is described fully in Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others (1998). The model is based on 50 m long finite differences and includes 40 stacked layers, each of which can be ascribed a unique hardness. We therefore define the following three layers: (i) a BZ with a relative hardness value of 10.74 (i.e. 10.74 × AUZ) that occupies the lowermost model layer (i.e. 2.5% of ice thickness); (ii) a LZ with a relative hardness value of 1.80 that occupies the next six layers above the BZ (i.e.15.0% of ice thickness); and (iii) a UZ with a relative viscosity of 1 that occupies the remaining 33 layers, extending 82.5% of the glacier’s thickness from the top of the LZ to the surface. These layers, which are modelled as constant proportions of ice thickness irrespective of the glacier’s geometry, are illustrated for the glacier’s current long section in Figure 2. Results of this multi-layered rheology model are compared with those of a single-layer rheology model run through two scenarios, a simulated glacier advance and a simulated glacier retreat.
Glacier advance
Under the advance scenario, the glacier is grown from zero thickness to its current surface profile by imposing an equilibrium-line altitude (ELA) of 2775 ma.s.l. on the ice-free valley profile. Each model configuration was run under a range of ice hardnesses from 0.10 a–1bar–3 (a value of 0.21 a–1bar–3 is considered to be a “standard” hardness of temperate glacier ice; Reference PatersonPaterson, 1994,p.96) through to 10 a–1bar–3 in order to provide the best-fit matches between the resulting modelled surface profile and that measured in 2001. This was achieved with AUZ = 1.2 a–1bar–3 for the multi-layer rheology model and A =7.0 a–1bar–3 for the single-layer rheology model. For intercomparison, the single-layer rheology model was also run with a hardness similar to the best-fit hardness calculated for the englacial ice in the layered model (i.e. A = AUZ =1.2 a–1bar–3).
Results of the model runs, presented in Figures 3 and 4 and summarized in Table 2, indicate that the multi-layer rheology model matches the actual surface profile closely, while the single-layer rheology model with A = 7 a–1bar–3 matches for most of the glacier’s length, but slightly underestimates thickness in the glacier’s accumulation area.These experiments yield rms ice-thickness deviations from those measured of 3.9 and 4.5 m, respectively. In contrast, the single-layer rheology model constrained with A = 1.2 a–1 bar–3 overestimates ice thickness substantially along the entire glacier’s length, yielding an areal excess of 31% of the current glacier’s long-section area (rms error =15.0 m).
Summary of two-dimensional modelling results, expressed in terms of long-section ice thicknesses or total long-section areas
Advance scenario modelling results. The current measured glacier surface profile is given as a solid line, and the modelled steady-state profile is given as a dashed line: (a) multi-layer rheology model with AUZ =1.2; (b) single-layer rheology model with A = 7.0; (c) single-layer rheology model with A = 1.2.
Advance scenario modelling ice-thickness deviations from measured values: (a) multi-layer rheology model with AUZ = 1.2; (b) single-layer rheology model with A = 7.0; (c) single-layer rheology model with A = 1.2.
Glacier retreat
Under the retreat scenario, an ELA rise of 75 m (from 2775 m a.s.l. to 2850ma.s.l.) is imposed on the current glacier profile, and each of the three models described above is run to a steady-state response. Results, presented in Figure 5 and Table 2, indicate major glacier retreat and thinning. In the case of the multi-layered rheology model, the glacier retreats ∼1750m to a residual length of ∼700m, shrinking to 22% of its 2001 area. In the case of the single-layer rheology with A = 7.0 a–1bar–3, the glacier retreats rather less, by ∼1200m to a residual length of ∼1250m, or to 36% of its 2001 area. Intercomparison of these responses in terms of long-section area reveals corresponding major differences in the model predictions. The long-section area of the glacier, predicted using a single-layer rheology with A = 7.0 a–1bar–3, is 65% greater than that predicted by the multi-layer rheology model. In contrast, the single-layer rheology with A = 1.2 a–1bar–3 results in a dramatically smaller residual long-section area, having retreated ∼1900m to a residual length of only 550 m.This is 31% smaller than the multi-layer rheology long-section area, and 42% smaller than that resulting from the single-layer rheology with A = 7.0 a–1 bar–3 model.
Retreat scenario modelling results. The current measured glacier surface profile is given as a solid line, and the modelled steady-state profile is given as a dashed line: (a) multi-layer rheology model with AUZ =1.2; (b) single-layer rheology model with A =7.0; (c) single-layer rheology model with A =1.2.
Summary and Discussion
Reconstructing the ionic composition of the three principal ice facies present at Glacier de Tsanfleuron allows their relative liquid-water content and softness to be approximated. If both these variables are scaled to a value of 1.0 in englacialfacies ice, their values in clear-facies ice and dispersed-facies ice become 1.80 and 10.74, respectively. Thus, other controls being equal, the basal ice layer present at Glacier de Tsanfleuron is approximately 10 times softer (or less viscous) than standard englacial ice that forms the bulk of the glacier. Introducing these relative viscosities into a layered numerical model of ice flow results in markedly different predicted responses to ELA perturbations relative to those of an otherwise similar model incorporating a spatially uniform ice viscosity. For example, increasing the ELA by 75 m results in modelled steady-state glacier long sections that reduce to 15–36% of its current long-section area. Importantly, these data indicate that if the same ice hardness term (A =1.2 a–1bar–3) were used as the basis for a multi-layered and a single-layered model then the latter would overestimate retreat by 31% relative to the former. The physical explanation for this overestimation probably lies in the relative stiffness of the single-layered glacier and its inability to extend downslope under a low surface gradient. In contrast, if optimized ice hardnesses are adopted, i.e. AUZ = 1.2 a–1 bar–3 in the multi-layered model and A = 7.0 a–1bar–3 in the single-layer model, the latter underestimates retreat by 65% relative to the former. In this case, the single-layered glacier is relatively soft and is correspondingly able to extend further downslope under a lower surface gradient.
This study includes several approximations and assumptions. For example, the expression of ice layering in our model is relatively crude, being modelled as uniform proportions of ice thickness irrespective of location along the glacier’s long section or the glacier’s overall geometry. We have also focused solely on the rheological role of liquid-water content, to the neglect of variations in other physical properties. This research consequently represents a first approximation of the magnitude and significance of systematic spatial variations in the physical structure and rheology of temperate glaciers. However, our results clearly indicate that such variations may be important for modelling predictions of the response of temperate ice masses to future climate change.
Acknowledgements
We thank L. Janssens, R.D. Lorrain, A. Khazendar and R. Hubbard for fieldwork assistance. This work was funded by U.K. Natural Environment Research Council grant GR9/2026. J.-L.Tison is Research Associate at the Fonds National de la Recherche Scientifique, Belgium.
