Introduction

Due to the limited accessibility of the subglacial and englacial environment, the determination of subsurface ice velocities and components of the strain rate and stress is a technically difficult and expensive task. Such information is of crucial importance for understanding basal processes, such as variable basal sliding and its consequences for the near-basal strain field.

A number of methods and techniques, mostly in boreholes, have been applied to probe the englacial velocity and strain fields. An often used method is the observation of ice deformation in an open borehole with inclinometry (Reference Gerrard, Perutz and RochGerrard and others, 1952; Reference SharpSharp, 1953; Reference RaymondRaymond, 1971; Reference Hooke and HansonHooke and Hanson, 1986; Reference Hooke, Pohjola, Jansson and KohlerHooke and others, 1992). This method was used to determine longitudinal strain in a chain of boreholes (Reference Shreve and SharpShreve and Sharp, 1970), to identify extrusion flow type velocity profiles (Reference Hooke, Holmlund and IversonHooke and others, 1987), to determine local sliding velocities (Reference Copland, Harbor, Minner and SharpCopland and others, 1997) and to investigate the effect of a cold margin on ice flow (Reference Moore, Iverson, Brugger, Cohen, Hooyer and JanssonMoore and others, 2011). Reference Harper, Humphrey and PfefferHarper and others (1998, Reference Harper, Humphrey, Pfeffer, Huzurbazar, Bahr and Welch2001) used arrays of boreholes to determine three-dimensional flow and strain fields. Various methods were applied to measure vertical strain rates in open boreholes (Reference PatersonPaterson, 1976; Reference Raymond, Rogers, Taylor and KociRaymond and others, 1994; Reference GudmundssonGudmundsson, 2002; Reference Hawley, Waddington, Morse, Dunbar and ZielinskiHawley and others, 2002; Reference Sugiyama and GudmundssonSugiyama and Gudmundsson, 2003). Reference Schwerzmann, Funk and BlatterSchwerzmann and others (2006) used a combined inclinometer/caliper probe to measure the velocity profile together with the closure, cross-sectional deformation and vertical normal strain of a borehole.

The installation of continuously measuring inclinometers in boreholes allows us to perform measurements of high temporal resolution (Reference Gudmundsson, Bauder, Lüthi, Fischer and FunkGudmundsson and others, 1999; Reference Amundson, Truffer and LüthiAmundson and others, 2006). This method continuously observes the tilt and azimuth of the axes of instruments that are inserted into the ice. Modeled tilt curves that are calculated from assumed or modeled flow fields were fitted to the measured tilt curves to interpret the measurements.

An in situ measurement of the stress-strain-rate relation of glacier ice can only be performed if all stress components and strain components are measured simultaneously at the same location within the ice (Reference Pfeffer, Humphrey, Amadei, Harper and WegmannPfeffer and others, 2000; Reference Marshall, Harper, Pfeffer and HumphreyMarshall and others, 2002) without additional assumptions on the flow field or numerical modeling. In this work we propose and discuss a technique that allows us to meet these requirements for measurements of strain rates. The practical application of this technique may not be trivial, and in fact, the analysis was carried out as a result of difficulties in an experiment on Rhonegletscher, Switzerland. The coupling of the instruments to the ice in boreholes in temperate ice may take several months and seems to be unstable even then.

For the evaluation of the inclinometer data we assume that the instrument behaves as a Lagrangian unit vector attached to the ice, which leads to a time-evolution equation for the unit vector. This equation is (to lowest order in the time-step) identical to the one Reference Gudmundsson, Bauder, Lüthi, Fischer and FunkGudmundsson and others (1999) solved by an explicit Euler method. The evolution equation for the unit vector can be solved analytically with two restrictive assumptions: (1) the flow field is well represented with the plane-strain first-order approximation (Reference BlatterBlatter, 1995) and (2) the strain components are constant over the measurement period.

In the following sections, we give some technical information on the inclinometer/magnetometer probes and their calibration. The mathematical theory of strain measurements with inclinometer probes is outlined and discussed, and we present an analysis of time series of zenith and azimuth angles of probe axes to determine the strain components if only one sensor is available and show some examples of time series observed in Alpine glaciers. In the final section (Conclusions and prospects), we discuss the possibilities for strain measurements and propose experiments to further develop the method.

Inclinometer Probe

The inclinometer probe consists of a cylindrical aluminum casing of 46mm diameter and 63 mm length. A two-axis inclinometer of type SCA100T-D01 (VTI Technologies, http://www.vti.fi/en/products/inclinometers) with a nominal resolution of 0.044 mrad, and a three-axis magnetometer of type MicroMag3 (PNI Corporation, http://www.pnicorp.com/products/products) with a nominal resolution of 015 μT allow us to measure the absolute direction of the probe axis in space and, furthermore, to determine the absolute orientation of an orthogonal set of vectors that are fixed with respect to the casing. The limiting factor for the accuracy is given by the magnetometers and the large uncertainty in the orientation of the Earth's magnetic field at depth in the ice.

Calibration and processing of data

The combined inclinometer/magnetometer probes provide the strength of the magnetic field along three axes and the strength of the gravity field along two axes,

(1)

where R_i are the magnetic and gravity field strengths, S_i are the readings of the sensor output, a_i are the slopes of linear functions and b_i are the offsets at zero field strength. The ten unknown parameters, a_i and b_i , must be determined by a calibration procedure.

Ideally, the calibration should be performed under laboratory conditions. If, however, the absolute value and the orientation of the magnetic field at the site of the measurement is not known, the problem can still be treated with the same number of unknown parameters. In this case, we simply assume the R_i to be normalized fields, and the parameters depend on the absolute field strengths. For the magnetic field this means that the calibration has to be repeated on site, as the geomagnetic field there may be different from that in the laboratory.

Given a parameterization of the pure rotations in three dimensions, three parameters completely describe the orientation of the instrument with respect to a fixed coordinate system at a given measurement. For N measurements, there are 3N unknown orientation parameters, i.e. a total of 10 + 3N unknowns (including instrument parameters). Every measurement results in five measured values of magnetic and gravity components along the corresponding axes. Thus, for N measurements, we obtain 5N readings. For N = 5 the problem is closed, for N > 5 measurements, the problem is overdetermined, and the solution can be optimized by a least-squares optimization procedure, yielding both ten instrument parameters and 3N orientation parameters of best fit. In detail, we need to minimize the sum

(2)

where ϕ_j is the set of orientation parameters at measurement j, and R_i(ϕ_j) are the normalized fields projected to the instrument’s coordinate system. The minimization with respect to ϕ_j, a_i and b_i is best done using the Gauss–Newton algorithm (e.g. Reference BjörkBjörk, 1996).

The calibration was done in several steps. First, a laboratory calibration measurement was used to determine the four gravimeter parameters. Once the instrument was at the field site, measurements close to the borehole were carried out to calibrate the magnetic sensors. The data processing of the sensors in the boreholes can be done using the same optimization procedure as for the calibration with the instrument parameters from the calibration kept constant and optimizing only the orientation parameters. This allows us to calculate for every measurement an orthogonal matrix, Q, that describes the orientation of the instrument with respect to a reference frame fixed in space, given by the gravity and geomagnetic fields. Q contains as rows three unit vectors each describing the orientation of one of three axes of the instrument with respect to the external reference frame. The angles θ and ϕ, that are used in the next sections, are the zenith (i.e. tilt from vertical) and azimuth angles of the row of Q which corresponds to the longitudinal axis of the instrument. It would be desirable to do a similar time-series analysis for the whole orthogonal system; however, it is not yet clear how similar theoretical evolution curves for the whole system can be obtained. Details of the underlying mathematics are given in Appendix A.

Mathematical Theory

A thin cylindrical probe inserted in the ice is assumed to move in the same way as a unit vector given by the direction of the image of a Lagrangian coordinate increment with respect to deformation. Using the notation of Reference HutterHutter (1983), we denote by dX_A a vector increment in the reference (Lagrangian) configuration, and by dx_i-(X _A) the corresponding vector in the present (Eulerian) configuration. Using the Einstein summation convention, we define the unit vector through its components, e_i , as

(3)

so it points in the direction of the vectorial quantity dx. The temporal evolution of the unit vector is given by

(4)

where

(5)

denotes the components of the velocity gradient. That is, the evolution of e_i is given by the projection of L_ije_j on a plane perpendicular to the vector e_i . A detailed derivation of Eqn (4) is given in Appendix B.

For practical purposes, we assume Cartesian coordinates with the base unit vectors n₁, n₂ and n₃, such that the unit vector e = e₁ n₁ + e2 n₂ + e₃ n₃. Furthermore, we use spherical coordinates with 0 ≤ θ < π, ≤ ϕ < 2π such that

(6)

Thus, the time derivative of the unit vector, e, has the form

(7)

where the vectors

(8)

are orthogonal to vector e and orthogonal to each other,

(9)

Scalar multiplication of Eqn (4) with a and b, respectively, yields

(10)

and

(11)

The second part on the right-hand side of Eqn (4), e_i(e_jL_jk e_k), does not have any influence on and as it is colinear with e, and thus orthogonal to both vectors a and b. Equations (10) and (11) are two linear equations for the nine components L_ij of the velocity gradient.

The analytical solution can be used for a regression with the measured time series of zenith angles. Interestingly, the azimuth angle does not occur in the solution for θ, except for its initial value, ϕ₀, which can be used as a fitting parameter in the regression.

In principle it is possible to give an analytical solution including L₁₂ ≠ 0 as well but this is much more complicated than Eqn (16), which makes it much less useful for purposes of field data analysis.

Analysis of time Series

One absolute mono-axial (linear) inclinometer/magnetometer yields one value for and at one time, and thus only two equations for the nine unknowns, L_ij . Three linear independent unit vectors, e_k k = 1, 2, 3 with measured angles θ_k and ϕ_k yield six equations. One additional spherical inclinometer/magnetometer, tracking the rigid rotation of the velocity field, would yield three more equations, and would thus allow one to close the system of equations.

However, in many situations interesting for glaciologists, fewer than nine components L_ij are of interest. A suite of assumptions and approximations may reduce the number of unknowns. Some assumptions may not be trivial, though, and require careful testing to estimate the possibly imposed errors.

For a typical channel-like valley glacier it is reasonable to assume that in a coordinate system with x₁-direction along the flow and x₃ vertically (parallel to gravity), the transversal velocity component and its derivatives are comparably small, thus

(12)

Along a central flowline, where the transverse velocity vanishes and the along-flow velocity reaches its transversal maximum, the derivatives with respect to x ₂ become small. Thus, the first-order plane-strain approximation for incompressible ice, L₃₃ = −L₁₁, with the only non-vanishing components L₁₃ and L₃₃ reduces Eqns (10) and (11) to

(13)

(14)

In principle, Eqns (13) and (14) are linear equations for L₁₃ and L₃₃, and thus could be used to calculate the temporal evolution of the velocity-gradient components. The analysis depends on the calculation of time derivatives of the measured angles, θ and ϕ. Depending on how noisy the data are, this procedure is more or less unstable and requires a smoothing of the data. Even then, the evaluation of Eqns (13) and (14) resulted in temporally drifting values for L₁₃ and L₃₃, that should result in a temporal drift of the surface velocity, which could not be observed. Apart from that, L₁₃ and L₃₃ depend explicitly on the evolution of ϕ. The measurement of this angle requires the use of magnetometers, which are much more susceptible to environmental perturbations and thus are less accurate than gravity-based inclinometry. Furthermore, Eqns (13) and (14) are only valid in a coordinate systemwith the x1-axis oriented parallel to the flowdirection. The azimuth angle measured by an instrument fixed to the ice is difficult to orient absolutely in space, and thus to align with respect to the ice flow direction.

With the assumption that L₁₃ and L₃₃ are constant over the period considered, Eqns (13) and (14) can be solved analytically for the time evolution of the zenith angle, θ = θ(t),

(15)

with

(16)

A detailed derivation of Eqns (15) and (16) is given in Appendix C. A numerical evaluation of this equation reproduces the same curves as the Euler scheme used by Reference Gudmundsson, Bauder, Lüthi, Fischer and FunkGudmundsson and others (1999) to a high accuracy (not shown here in detail).

In Figure 1, two examples for possible tilt curves according to solution Eqn (15) are given. There are always two asymptotic values, π/2 for L₃₃ t < 0, and

Fig. 1.

Examples of possible tilt evolutions according to Eqn (15). Solid curve: L₁₃ = 3L₃₃, θ₀ = 0. Dashed curve: L₁₃ = 2L₃₃, θ₀ = 0.2, ϕ₀ = π/2.

(17)

on the half-axis with L₃₃t > 0. Thus, if the motion is dominated by shear, L₁₃ » L₁₃, both asymptotes are close to π/2, whereas for L₁₃ ≈ L₃₃ significantly lower asymptotes are possible. If ϕ₀ = 0 or ϕ₀ = π (depending on the sign of L₃₃) or θ ₀ = 0, the tilt reaches a minimum of θ = 0 and has a kink at the minimum (solid line in Fig. 1). In all other cases the minimum is greater than zero and is passed smoothly; in these cases e never reaches a vertical position.

Results

Field observations

Here we demonstrate the application of Eqns (15) and (16) to fit tilt data from field experiments. We consider field data from Rhonegletscher and Gornergletscher, both situated in the Swiss Alps. The Rhonegletscher data stem from a chain of five inclinometer/magnetometer probes that was installed in a borehole in summer 2009. The data were read every 5 min, except during winter2010/11, when the measurement interval was extended to 8 hours in order to cope with power supply problems. The instruments are labeled A, B, C, D and E, with corresponding approximate positions of 4, 23, 40, 58 and 75 m above the glacier bed (ice thickness 119 m). The borehole is situated in temperate ice, and was permanently water-filled during the experiment. The data measured by sensor E are extremely noisy, probably due to bad coupling of the sensor to the ice. Due to large shear rates, sensor A ran out of the measurement range of the instrument within 1 year. For comparison, we also show measurements obtained with the same type of instrument installed in summer 2006 in a borehole in Gornergletscher. InstrumentG1 is located 127 m above the bed (ice thickness 271 m) in cold ice of —1°C above a basal temperate layer of ~70m thickness, and can be assumed to be coupled tightly to the ice. Instrument G2 was inserted at 38 m above the bed in the temperate ice. In Figure 2 the measured tilt evolutions are displayed together with least-square fits obtained with Eqns (15) and (16).

Fig. 2.

Evolution of tilt angle, θ, on Rhonegletscher and Gornergletscher. The thin lines are the fitted curves according to Eqns (15) and (16). (a) Rhonegletscher, sensor B, L₁₃ = 5.55 × 10^-9 s^-1, L₃₃ = 3.40 × 10^-9 s^-1. (b) Rhonegletscher, sensor C, L₁₃ = 2.73 × 10^-9 s^-1, L₃₃ = 8.18 × 10^-10 s^-1. (c) Rhonegletscher, sensor D, L₁₃ = 1.39 × 10^-8 s^-1, L₃₃ = 1.10 × 10^-7 s^-1. (d) Gornergletscher, sensor G1, L₁₃ = 2.69 × 10^-9 s^-1, L₃₃ = 1.86 × 10^-9 s^-1.

Daily variations

At some times, the long-term evolution of the zenith angles is modulated by daily variations. These short-term variations cannot be captured by analytical tilt curves that are optimized for long periods of measurements with constant strain rates. The lowest (sensor A) and the second lowest (sensor B) instruments (Fig. 3a and b) show variations in shear strain in phase, whereas for the next highest instrument (sensor C) they are in antiphase (Fig. 3c). Interestingly, the shear rates change sign during one cycle. These variations may be caused by variations in basal sliding; however, due to possibly insufficient coupling of the instruments to the ice in the open and partially water-filled borehole, an influence of the daily variation in the water flux through the borehole cannot be excluded. Further interpretation of these patterns lies beyond the scope of this paper. Figure 3d demonstrates an example for diurnal variations in the tilt evolution inside a temperate basal layer in Gornergletscher.

Fig. 3.

Diurnal variations of tilt angles, θ, on Rhonegletscher and Gornergletscher. The scatter of the single points indicates the angular resolution of the sensor around 0.1mrad. (a) Rhonegletscher, sensor A. (b) Rhonegletscher, sensor B. (c) Rhonegletscher, sensor C. (d) Gornergletscher, sensor G2.

Conclusions and Prospects

The analysis of the tilt evolution observed with the inclinometer/magnetometer probes is based on the assumption that a thin cylindrical probe inserted into glacier ice follows the same movement as a unit vector, which is given as the direction of the image of a Lagrangian coordinate increment with respect to deformation. The time-evolution equation for this unit vector is defined by the flow field, or more precisely, by its spatial derivatives. For simple flow fields (simple shear with some vertical and horizontal stretching), explicit analytical solutions to the evolution equations in polar coordinates were used to fit field data.

The above findings allow us to calculate theoretical tilt curves, which should be recorded by an instrument sensitive to inclination that is fixed within the material, and can therefore be helpful for calibration/validation purposes or for the reconstruction of a local flow field from tilt measurements. Of course, one could try to avoid the mathematical difficulties of solving the evolution equations analytically and just rely on a numerical integration of the coupled system. The analytical solution has, however, the great advantage that we can write down the evolution of the tilt angle, θ, without simultaneously solving the equation for the azimuth angle; instead, the evolution of θ is parameterized by the initial azimuth angle, ϕ₀.

The explicit time evolution, θ(t), has been used to fit data from inclinometer tilt experiments on Rhonegletscher and Gornergletscher, showing its ability to reproduce the measured curves. With the assumption of incompressible ice flow and the first-order plane-strain approximation, strain- rate measurements can be performed with much simpler instruments, only measuring the zenith angle of the probe axis.

Tilt measurements with high temporal resolution show examples of diurnal variations of strain rates. This phenomenon has been observed in different glaciers, such as the entirely temperate Rhonegletscher, the polythermal Gornergletscher and in the polythermal region of Jakobhavns Isbræ, Greenland ice sheet (Reference Marshall, Harper, Pfeffer and HumphreyLüthi and others, 2002). Reference Amundson, Truffer and LüthiAmundson and others (2006) measured similar features on Black Rapids Glacier, a large surge-type glacier in the central Alaska Range, USA. In these examples, the patterns in variations are similar in shape and magnitude. The mechanism causing these variations and the temporal coincidence with diurnal variations in englacial water pressure suggest that variations in sliding may explain the variations in the strain field. The partial reversibility of deformations perhaps indicates diurnal variations of normal strain (Reference GudmundssonGudmundsson, 2002) or elastic reactions of the ice base to variable basal water pressure (Reference Sugiyama, Bauder, Weiss and FunkSugiyama and others, 2007).

An unsolved question concerns the time evolution of an orthogonal system of unit vectors fixed to a rigid body inside the material. If all unit vectors evolve as described above, the orthonormality condition may be violated in the presence of a shear motion. One therefore has to define a way to couple the rigid body to the deforming material such that orthonormality of the coordinate system fixed to it is kept. This problem can probably be addressed by modeling the motion of a rigid or hard elastic body in a deforming viscous fluid, which, however, lies beyond the scope of this paper.

One suggestion that can be drawn from the above considerations is a possible in situ measurement of the local viscosity of glacier ice. A system of three thin cylindrical inclinometer/magnetometer probes with independent orientations and the same type probe with spherical shape to measure the rigid rotation components, allows us to close the system of equations for the nine velocity-gradient components. Combined with a probe with adequate pressure sensors to measure the stress components as suggested by Reference Pfeffer, Humphrey, Amadei, Harper and WegmannPfeffer and others (2000), this yields, in principle, the complete information to compute the ice viscosity