Spatial variability in the water content and rheology of temperate glaciers: Glacier de Tsanfleuron, Switzerland

Abstract

The physical properties of eight ice cores recovered from along a flowline at Glacier de Tsanfleuron, Switzerland, have led to the identification of three distinctive internal zones. We use variations in the bulk ionic chemistry of these zones to approximate their relative liquid-water concentrations and ice viscosities. Results suggest that relative bulk water concentration and ice softness vary by over an order of magnitude between the zones. Implications of this variability for predictions of the glacier’s response to climate change are evaluated by incorporating these relative softnesses into a multi-layered (two-dimensional) model of ice flow. Model output is compared with that from an identical model constrained with a spatially uniform ice viscosity under advance and retreat modelling scenarios. The former scenario is used to tune viscosity by growing a glacier to its present long-section geometry, resulting in best-fit ice hardness values of 1.2 a¹ bar^–3 for the englacial ice in the multi-layered model and 7.0 a^–1 bar^–3 for all of the ice in the single-layered model. Both result in close approximations to the current long profile, yielding rms deviationsbetween measured and modelled ice thicknesses that are 5 5 m. In contrast, a single-layered model constrained with a hardness of 1.2 a^–1 bar^–3 overestimates the current measured long-section area by 31%, having a rms ice-thickness error of 15.0 m. Under the retreat modelling scenario, which gauges the response of the glacier to an imposed 75 m rise in equilibrium-line altitude (ELA), the multi-layered model predicts a long-section area reduction of 78%. This contrasts with a reduction of 64% for the single-layered model (hardness = 7.0 a^–1 bar^–3) and 85% for the single-layered model (hardness = 1.2 a^–1 bar^–3).

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., Payne 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 Duval (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 km², 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 (Tison and Lorrain, 1987; Hubbard and Sharp, 1995) and subglacialcarbonate crusts exposed on its proglacial bedrock surface by retreat over the past ∼150 years (Hallet and others, 1978; Lemmens and others, 1983; Souchez and Lemmens, 1985; Sharp and others, 1990; Hubbard and Hubbard, 1998; Hubbard and others, 2000b).

Fig. 1. 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

Hubbard 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 10₂ g m–₃. 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 10₄ g m_–3 (Table 1).

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 (Hubbard and others, 2000a; Tison 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.

Fig. 2. 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⁺ , Ca^{2 +} , Mg ^{2 +} , NO₃ ²–, C l ^– and SO₄ ²– concentrations by ion chromatography (Dionex DX-100). Laboratory procedures are described elsewhere (Hubbard 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–¹ and on the order of ±100% at concentrations of ~1 μeq L–¹. HCO^{3 –} concentrations were determined by charge balance to an accuracy of ± 1 0% (Fairchild and others, 1999; Hubbard 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 (Nye and Frank, 1973; Glen and others, 1977; Wolff and Paren, 1984). Following Nye (1991) and Mader (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 M_b and the temperature depression 9 of the ice below its pressure-melting point:

(1)

where the constant K (in K m ³ g–¹) is given by the slope of a

(2)

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, M_a. Thus, the original bulk ionic concentration of the ice prior to sampling, M_b , is given by

(3)

In order to calculate the concentration of solute acquired during sampling, M_a, we apply acquisition rates derived from laboratory dissolution experiments to our laboratory sampling conditions. Following Lerman (1979), the concentration C of a solute in solution is expressed as function of time t as

(4)

where Cs (g cm–₃) is the ultimate (steady-state) concentration of the material in solution, C0 (g cm–₃) 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. Brown 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 (m²) per unit volume of water (m³); m–¹) of the sediment concerned, yielding, for Ca^{2 +} ,

^*The glacier’s actual long-section area, measured in 2001, is 116780 m².

(5)

and

(6)

Here, we use these authors’ acquisition rates for Ca^{2 +} because measured bulk solute at Glacier de Tsanfleuron is dominated by the dissolution of CaCO₃ (yielding the principal ions Ca^{2 +} and HCO₃ ²–). We double the predicted Ca^{2 +} concentration to account for the balancing HCO₃ ^– .

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 Hooke and Iverson (1995). Specific surface area by mass, SSA_m (m² g–¹), 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 SSA_m is multiplied by the measured concentration of that sediment in the ice facies (gm^{– 3}) to yield the SSA (m–¹). Results of this process yielded SSA_m values of 1.58 m² g–¹ for the englacial-facies debris, and 0.555 m² 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 10² gm–³), 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 10⁴ gm–³ . These calculations yield an overall sediment concentration for the 15 m thick clear-facies ice layer of 7.67 6 10² gm–³ . The corresponding SSA_m 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 m² g–¹. Resulting values of SSA are 3.00 6 10² m–¹ for the englacial facies, 4.70 6 10² m–¹ for the clear facies, and 2.62 6 10⁴ m–¹ for the dispersed facies (Table 1). Substituting these values into Equations (5) and (6) yields values of C0 and C_s of 0.14 and 0.44 gm–³ for the englacial facies, 0.22 and 0.69 gm–³ for the clear facies, and 12.2 and 38.4 g m–³ 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 C₀ and C_s derived above, into Equation0 (4) yields values of the solute concentration acquired during sampling C of 0.50 g m–³ for the englacial facies, 0.79 g m– ³ for the clear facies and 43.94 g m–³ 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 (M_b) of 0.68 g m–³ for the englacial facies, 1.23 g m–³ for the basal clear facies and 7.34 g m–³ for the basal dispersed facies (Table 1). Substituting these values of M_b 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–¹ bar–³ ) in Glen’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 (A_LZ) and the BZ ice (A_BZ) relative to that of the UZ ice (A_UZ) in addition to empirical tuning. Duval (1977) reported a strong linear correlation between measured water content W and effective strain rate, yielding:

(7)

Thus, the relative hardnesses of the clear facies (A_LZ) and dispersed facies (A_BZ) scale linearly with reconstructed values of W' to become 1.80 A_UZ and 10.74A_UZ, 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 (Blatter, 1995), and therefore includes the effects of longitudinal or normal deviatoric stresses, is described fully in Hubbard 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 × A_UZ) 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; Paterson, 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 A_UZ = 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 = A^UZ =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).

Table 2. Summary of two-dimensional modelling results, expressed in terms of long-section ice thicknesses or total long-section areas

Fig. 3. 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 A_UZ =1.2; (b) single-layer rheology model with A = 7.0; (c) single-layer rheology model with A = 1.2.

Fig. 4. Advance scenario modelling ice-thickness deviations from measured values: (a) multi-layer rheology model with A_UZ = 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.

Fig. 5. 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 A_UZ =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. A_UZ = 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.