Introduction

Glaciers and ice sheets on soft beds flow by some combination of ice creep, basal sliding and subglacial sediment deformation (Reference Alley,Alley, 1989). It remains unclear what processes control the partitioning between sliding and sediment deformation, and without new instruments and extensive subglacial measurements it is doubtful that this question can be settled. To this end, we have developed an instrument, the UBC ploughmeter, to sense the mechanical interaction between a glacier and its bed (Reference Fischer, and Clarke,Fischer and Clarke, 1994). The essential idea is that a steel rod, with its tip protruding into subglacial sediment, will bend elastically as the tip is dragged through the sediment by movement of the overlying ice Fig. 1a), much as a stick bends if it is used to inscribe lines on a sandy beach. The main body of the rod is embedded in ice, while the tip is instrumented with strain gauges that respond to bending in two mutually orthogonal directions Fig. 1b). Detailed information on the construction, calibration, installation and theory of this device is given in Reference Fischer, and Clarke,Fischer and Clarke (1994). For the purpose of the present paper it is only necessary to understand that the ploughmeter measures the bending moment acting near the rod tip and that this moment is influenced by the depth of penetration, the sliding rate and bed properties.

Fig. 1. (a) Schematic diagram of ploughmeter operation. (b) Arrangement of strain gauges near the tip of the steel rod.

For a homogeneous, fine-grained substrate, modelled as a linear viscous fluid having viscosity μ, the force per unit length F acting near the tip is approximated by

where U is the relative velocity between rod and sediment (here assumed to equal the basal sliding rate) and L is a characteristic length scale (Reference Fischer, and Clarke,Fischer and Clarke, 1994). For a heterogeneous or clast-rich sediment, a purely linear viscous assumption is inappropriate and the concepts of mean free path and clast-collision frequency suggest themselves. The mean free path is the typical travel distance between clast collisions; collision frequency is the inverse of the typical time interval between collisions. Potential complicating factors are temporal variations in sliding rate and spatial variations in bed granulometry.

Observations

During the 1991 summer field season, two ploughmeters were installed at the bed of Trapridge Glacier, Yukon Territory, Canada. The location and a detailed description of the Trapridge Glacier study area are given in Reference Clarke, and Blake,Clarke and Blake (1991). Figure 2 shows roughly 60 d of observation for ploughmeter 91PL02. During the course of these measurements, we collected data at 20 minute intervals. We estimate that the ploughmeter was inserted ∼0.15 m into the basal sediment. Figure 2a indicates the bending force applied to the tip of the ploughmeter, and Figure 2b indicates the azimuth of the force with respect to the internal coordinates of the ploughmeter. During the entire observation period, both the force record and the azimuth record display a jagged appearance with rapid and large fluctuations. Variations in azimuth of the applied force could result from rotation of the ploughmeter about its long axis, from temporal variations in the direction of glacier flow, or from translational motion in a direction perpendicular to the glacier flow. Because we believe that ploughmeters are firmly gripped by glacier ice (Reference Fischer, and Clarke,Fischer and Clarke, 1994), it is unlikely that they will undergo large and rapid back-and-forth rotations within a borehole Fig. 2b). To further interpret the data, we postulate that azimuth variations due to changing glacier flow directions are generally long-term (cycles of 1 d or more) and that variations due to cross-flow motion are generally short-term (cycles of on the order of hours or less). In our analysis, we used a Gaussian low-pass filter (cut-off frequency ≈ 0.074 d⁻¹) to extract the slowly varying component of the azimuth time series (dashed line in Fig. 2b) and attributed this slowly varying component to changes in the glacier flow direction. Any remaining high-frequency disturbance on the azimuth record indicates non-zero cross-flow motion. By aligning the slowly varying azimuth component with the glacier flow direction, we can decompose the net force values Fig. 2a) into down-flow and cross-flow components Fig. 3). To effect this decomposition we assume that the principal direction of ploughmeter motion is down-glacier.

Fig. 2. Data from ploughmeter 91PL02. (a) Force record indicating force applied to the tip of the ploughmeter. (b) J Azimuth of the force with respect to the internal coordinates of the ploughmeter. Dashed line represents the slowly varying azimuth component (see text for details).

Fig. 3. Decomposition into down-flow (upper trace) and cross-flow (lower trace) force components of the data from ploughmeter 91PL02 shown in Figure 2.

Interpretation

The jagged appearance of the force record Fig. 2a) and the lack of an apparent correlation between the records of down-flow and cross-flow force components Fig. 3) suggest that ploughmeter 91PL02 is interacting with a clast-rich area of the glacier bed. The basis for this suggestion lies in the idea that the spikes in the data record are the result of the ploughmeter colliding with individual clasts as it is dragged through basal material. An alternative interpretation is that the spikes reflect interactions with groups or “bridges” of grains. During shearing of a granular material, grain alignment can become sufficiently coaxial that the resulting bridges support anomalously high shear and normal stresses (Reference Hooke, and Iversion,Hooke and Iverson, 1995). However, as these bridges fail by crushing of the particles, there is a tendency for the material to develop a self-similar size distribution in which particles of the same size are separated from each other (Reference Sammis,, King and Biegel,Sammis and others, 1987). In such a spatial arrangement of particles, intergranular stresses are more uniformly distributed and the tendency for mechanical heterogeneity is minimized (Reference Iverson,, Hooyer, and Hooke.Iverson and others, 1995). Hence, because the particle-size distribution for Trapridge sediment was shown to approximate a fractal distribution (Reference Fischer,Fischer, 1995) we prefer the interpretation of a ploughmeter colliding with individual clasts.

To support our idea we simulated collisions between a ploughmeter and individual clasts as the ploughmeter moves through a heterogeneous subglacial sediment. We developed a simple numerical model to describe the response of a ploughmeter to such forcing by interaction with a synthetic till. For simplicity, we only Considered nearest neighbour collisions. The model till was constructed by filling a volume with spheres of different sizes. The spheres represent Clasts in a matrix of fine-grained solids and water-filled voids. We assumed a porosity of 0.3 and a clast-size distribution that is based on the granulometry of a basal till from Trapridge Glacier (Reference Clarke,Clarke, 1987). Only clasts of the seven largest size classes were included in the model, corresponding to spheres of diameters in the range 4–32 mm (Table 1). Grains of smaller size were assumed in form the matrix of fines. To achieve a random yet spatially fractal distribution of clasts within the till, the spheres were arranged inside the Volume using a method similar to the way a dense concrete is created. As a first step, the largest spheres were randomly assembled. Subsequently, spheres of the next smaller size class were placed into the voids left between the larger ones. This second stage was then followed by adding yet smaller spheres which fill the gaps between the first and second sets of spheres. By recursively placing smaller spheres in between the larger spheres, this filling process was continued untill the smallest size class had been reached. The resulting cascades of smaller and smaller spheres filling the gaps between larger ones is similar to the packing of circles in the fractal Apollonian gasket Fig. 4).

Fig. 4. Cross-sectional cut through the model till. The area shown is 0.1 m by 0.1 m in size. The light shaded regions represent the matrix of fine-grained solids and water-filled voids. Inset shows a mathematically perfect Apollonian gasket (adapted from Reference Mandelbrot,Mandelbrot, 1983, p. 170).

The transverse motion of a vertical cylinder through the volume represents the ploughmeter being dragged through subglacial material Fig. 5). The cylinder was taken to have a radius of 16 mm (Table 2) which is consistent with the cross-sectional dimension of the ploughmeter (Reference Fischer, and Clarke,Fischer and Clarke, 1994). To account for the insertion depth of the ploughmeter into basal sediment, the height of the model till volume was assumed to be 0.15 m (Table 2). We view the spheres as being suspended in a viscous medium. Thus, we treat the synthetic till as a solid-liquid dispersion and calculate the forces on the cylinder as it moves through this dispersion. Two parts contribute to these forces. The drag on the cylinder is calculated using Equation (1). This force contribution has only a component in the down-flow direction and, assuming that the glacier slides at a steady rate, is constant with time. The calculation of the forces that arise from collisions with the spheres is based on Stokes’ law

which describes the total drag of a sphere of radius R moving through a viscous fluid at velocity U (Reference Stokes,Stokes, 1851). The idea is that after impacting, the cylinder deflects the spheres to the side, thereby pushing them through the viscous medium. Because impacts are likely to be oblique with respect to the direction of motion of the cylinder, these forces have down-flow and cross-flow components.

Fig. 5. Schematic diagram of “clast collision” model. The vertical cylinder represents the ploughmeter moving through Subglacial sediment.

Instead of considering the motion of a cylinder through the fluid we change the frame of reference and analyze the uniform flow around the cylinder located at the origin (Fig.6). Fluid flow is assumed to be directed along the x axis. The spheres in the model till move at the fluid velocity and follow a trajectory that corresponds to a streamline. Collision between the cylinder and a given sphere occurs when the streamline is such that the sphere is not carried far enough around to clear the cylinder. If the velocity of the fluid at infinity is U then the velocity field at any point P can be expressed as (e.g. Reference Eskinazi,Eskinazi, 1962, p. 288)

where R_C is the radius of the cylinder, r is the length of the position vector for point P, and θ denotes the angle between r and the direction of fluid flow at infinity Fig. 6). The total velocity in the flow field is

Fig. 6. Uniform fluid flow around a cylinder (adapted from Reference Eskinazi,Eskinazi, l962, p.288).

When the cylinder and the sphere arc in contact, r is replaced by the sum of the radii of the cylinder and the sphere. It is the radial velocity component υ_r , that causes the sphere to be deflected from its original streamline. In order to calculate the force components in the down-flow and cross-flow directions we therefore have to decompose υ_r , into its x and y components Fig. 7),

Hence, substituting Equation [3a] into Equation (4) yields

where R_S denotes the radius of the sphere. With reference to Figure 7 we recognize that cosθ = x/(R_C + R_S ) and sinθ = y/(R_C + R_S ) and thus obtain two equations that describe the motion of the sphere around the cylinder along the x and y directions, respectively,

The temporal evolution of the force components in the down-flow and cross-flow directions that arise as the sphere is being deflected from its original streamline can now be easily calculated by substituting Equations (6) into Equation (2). However, these force calculations only account for interactions that take place while the cylinder and the sphere are in contact, and neglect any forces that are exchanged hydrodynamically before the collision has actually occurred. This “pre-collision” hydrodynamic interaction was simulated by ramping up the force components proportionally to the distance by which the streamline is deflected from horizontal Fig. 6).

Fig. 7. Velocity component on a sphere colliding with a cylinder.

In an attempt to describe the viscous behaviour of the synthetic till, it is useful to treat the solid liquid dispersion as an effective medium to which an effective viscosity can be ascribed. If a small number of particles is added to a liquid, its apparent viscosity increases as described by the relation (Reference Einstein,Einstein, 1906, Reference Einstein,1911)

where c is the volume fraction of the added particles (c ≪ 1), μ is the apparent viscosity of the dispersion, and μ ₀ that of the suspending medium. Following the principles of effective medium theory (Reference Van de Ven,Van de Ven, 1989, p. 539), we computed an effective viscosity for the model till by successive substitution. Algorithmically, this amounts to the following procedure. We start with the matrix of fine-grained solids and water-filled voids and treat this initial dispersion as a homogeneous material to which we ascribe a background viscosityμ _B. Subsequently, we add the spheres of the smallest size class to this material and calculate the increase in apparent viscosity according to Equation (7). We then regard the resulting “effective” dispersion with this new effective viscosity as the homogeneous material for the Succeeding substitution Step. This process continues until the spheres of the largest size class have been added to the dispersion. The associated final effective viscosity is what we might refer to as bulk viscosity of the synthetic till and is analogous to the apparent viscosity of subglacial sediments inferred from in situ or laboratory measurements (Reference Humphrey,, Kamb,, Fahnestock, and Engelhardt,Humphrey and others, 1993; Reference Fischer, and Clarke,Fischer and Clarke, 1994; Reference Porter,, Murray, and Dowdeswell,Porter and others, 1997). However, when analyzing the drag on a sphere of a given size class we substitute the viscosity that was calculated for this size class into Equation (2). The reason for not substituting the bulk viscosity is that the sphere moves in a “cage” formed by neighbouring spheres of the next larger size class. Because the typical distance of these neighbouring spheres is 3–4 times larger than the size of the sphere under consideration (Table 1), its motion is only affected by spheres of the same and smaller size classes. The foregoing thinking represents a standard application of effective medium theory (Reference Van de Ven,Van de Ven, 1989, p. 539) and is not original to this work.

As a first simulation Fig. 8), we substituted a background viscosity μ_B = 1.0 × 10¹⁰ Pa s and a fluid velocity U = 0.04 m d⁻¹ into our model equations. This background viscosity corresponds to a bulk viscosity of the synthetic till μ = 2.00 × 10¹⁰ Pa s which approximates values inferred the Trapridge sediment Reference Fischer, and Clarke,Fischer and Clarke, 1994). In our analysis, the fluid velocity is equivalent to the translational velocity of the cylinder (U = U_C ) and represents the differential velocity between ploughmeter and sediment which can he seen as a proxy for the glacier sliding rate. A value of 0.04 m d⁻¹ is a typical basal sliding velocity measured at Trapridge Glacier (Reference Blake,, Fischer, and Clarke,Blake and others, 1994).

Comparison of Figures 3 and 8 shows that the records generated using our model display many of the features seen in the real data. However, we note the lack of “broad” spikes in the synthetically generated records of down-flow and cross-flow force components Fig. 8). The absence of broad spikes could reflect the fact that clasts larger than 32 mm in diameter were not included in our model till (Table 1). Even though the largest clasts in a basal till sample from Trapridge Glacier (Reference Clarke,Clarke, 1987) were of this size, evidence from the ablation till in the forefield of the glacier indicates that the basal sediment layer contains clasts as large as cobbles and boulders. The upper limit of clast diameter of 32 mm, as determined by sedimentological size analysis (Reference Clarke,Clarke, 1987). is likely to be artificial and results from a bias towards smaller clasts during the sampling procedure. The sampling involved filling a 0.008 m³ pail with “typical till”, excluding material too large to shovel (Reference Clarke,Clarke, 1987). Upon closer inspection of Figure 3, we see that collisions with larger clasts seem to have occurred mainly during the interval 30 August-22 September, and that broad spikes are mostly absent during the remaining time period. For this reason we base the analysis in the following section mainly on the second half of this force record.

Fig. 8. Synthetically generated records of down-flow (upper trace) and cross flow (lower trace) force components. Note the similarity to Figure 3.

Estimation of Basal Sliding Rate and Sediment Viscosity

In the above section we have shown that the jagged appearance of the force and azimuth record of ploughmeter 91PL02 Fig. 2) can be interpreted as the ploughmeter interacting with a clast-rich sediment bed whereby the spikes in the data record indicate collisions with clasts. If we assume that the bed granulometry is spatially homogeneous and that temporal variations in sliding rate are negligible, the velocity with which the ploughmeter moves through the sediment is proportional to the rate of collision with clasts. In addition, assuming a linear viscous rheology for the subglacial sediment layer, a change in till viscosity and ploughmeter velocity results in a proportional change in the amplitude of the force records (Equations (1) and (2)). Therefore, by comparing the frequencies and the total power contained in the measured and synthetically generated force records (Figs 3 and 8), we are able to estimate the glacier sliding rate and an effective sediment viscosity.

Suppose that ploughmeter-sediment interactions yield a response r(t). The frequency content of this response can be represented by the power spectrum S(f) given by

where R(f) is the Fourier transform of r(t) and angular brackets denote time averaging. If r(t) is sampled at N points over a time period T, the total power of this time series is given by

A reasonable starting-point is to generate synthetic records of the force components in the down-flow and cross-flow directions for some plausible values of till viscosity μ ^s and glacier sliding rate U _s (e.g. Fig. 8). We denote these simulated times series by

and . We can now relate these simulated time series to those that were obtained from our field measurements Fig. 3). Presumably the values of till viscosity and sliding rate are different from those used in the simulation, hence

where μ ^m and U ^m are the values that we seek to determine. Note that a different sliding rate also causes the time axis to be rescaled because of a change in the rate of collision. Introducing the dimensionless ratios α = U ^m/U ^s and β = μ ^m/μ ^s gives

Consequently, the power spectra of

and are respectively related to those for the simulated time series and by

Here, we have made use of the Fourier transform similarity and addition theorems (e.g. Reference Bracewell,Bracewell, 1986, p. 122). Further, the power in the time series

and can be related to the power in the simulated time series as follows:

Note that the power is independent of the sampling interval (see Equation (9)). The implication of Equation (13) is that a simple multiplicative factor α²ß² can be used to adjust the power levels of the simulated time series

and to match the power levels of the observed time series and . Combined, Equations (12) and (13) suggest a procedure for determining α and β which leads to an evaluation of both till viscosity and sliding rate.

For both the data of ploughmeter 91PL02 Fig. 3) and the synthetically generated data Fig. 8), the records of the cross-flow force component were analyzed by taking the Fourier transform (Reference Press,, Teukolsky,, Vetterling, and Flannery,Press and others, 1992, p.490), and the power spectral densities were estimated using the periodogram method (Reference Press,, Teukolsky,, Vetterling, and Flannery,Press and others, 1992, p. 531). If plotted on log-log scales, we see that both spectral energy functions,

for the 91PL01 data Fig. 9) and for the synthetic data Fig. 10), show good power-law dependence on frequency over large parts of the spectra. Thus, the power spectra are of the form

which suggests that the rate of change of power is independent of scale. Considering that the ploughmeter interacted with a sediment having a self-similar particle-size distribution (Reference Fischer,Fischer, 1995), we are not surprised to find that the frequency-amplitude statistics of the force records are fractal. The least-squares log-log slopes from Equation (14) on the quasi-linear sections of the spectra give b = 2.7 for both cases. Because of this power-law scaling behaviour, we can represent the two power spectra

and in log-log space as two straight lines and, after defining γ = αβ rewrite Equation (12) as

where b ^m and b ^s are the slopes and a ^m and a ^s the intercepts of the power spectra for the measured and synthetic data series, respectively. The product γ = αβ can be evaluated by comparing the power levels of

with (Equation (13)). Now, Equation (15) can be solved to yield a value for α. This operation corresponds to shifting one of the power spectra with respect to the other, and we use this approach to determine the till viscosity and sliding rate that w ill achieve optimal overlap of the two power spectra.

Fig. 9. Power spectral density function for the record of the cross-flow force component shown in Figure 3. The dashed line represents Equation (14) with b ^m = 2.7.

Fig. 10. Power spectral density function for the record of the cross-flow force component shown in Figure 8. The dashed line represents Equation (14) with b^s = 2.7.

Table 3 summarizes the results of our simulations. Calculations were performed with a variety of different input values of the background viscosity and the velocity of the cylinder to demonstrate the robustness of the model. The velocity U _C was assumed to take on values of 20 to 80 mm d⁻¹, the latter of which corresponds to the surface velocity of Trapridge Glacier. Viscosity values were taken to span an order of magnitude and were in the range μ _B = 0.5−5.0 × 10¹⁰ Pa s. In all cases, we obtained consistent results for the predicted glacier sliding rate and till viscosity. The glacier sliding velocity of ∼45 mm d⁻¹ is in excellent agreement with direct measurements of sliding at the base of Trapridge Glacier using a “drag spool” (Reference Blake,, Fischer, and Clarke,Blake and others, 1994). Similarly, the value of the sediment viscosity of ∼2.0 × 10¹⁰ Pa s agrees well with estimates that had previously been obtained for Trapridge sediment from a theoretical analysis of the forces on a ploughmeter (Reference Fischer, and Clarke,Fischer and Clarke, 1994).

Concluding Discussion

In this paper, we have demonstrated that if the bed granulometry is known, the collision frequency and collision force as indicated by a ploughmeter can be used to obtain estimates of the basal sliding velocity and effective sediment viscosity. Our analysis assumed temporal variations in sliding rate to be negligible, a linear viscous rheology of the basal material and a spatially uniform granulometry. In the case of Trapridge Glacier we were able to confirm our velocity estimate by comparing it with direct measurements. Support for our viscosity estimate is given by previously derived values for Trapridge sediment.

In our analysis we have assumed that the differential velocity between ploughmeter and sediment is equivalent to the basal sliding rate of the glacier. Because glacier sliding is commonly defined as the motion between the base of the ice and the top of the bed, we have essentially neglected the contribution of the basal motion that arises from bed deformation. Therefore, the question of how a deforming sediment layer affects the ploughmeter record has to be addressed. If we assume that part of the basal velocity is due to bed deformation, sediments near the ice-bed interface have a velocity component in the down-flow direction. Thus, the differential velocity between ploughmeter and surrounding sediment is generally less compared to the case where the ploughmeter is moving through a non-deforming material. Accordingly, for this case, where the ploughmeter is moving through an already shearing material, both the frequency of collision and the force of collision are reduced. We outline the implications for our model by comparing the mode of motion in which 100% of the basal velocity is due to glacier sliding to that in which the basal velocity is partitioned between 50% glacier sliding (Reference Blake,, Fischer, and Clarke,Blake and others, 1994) and 50% bed deformation in a 0.15 m thick deformable layer with a shear strain rate that is constant with depth. To obtain similar frequencies and power levels in the force records generated with the two modes of basal motion, the basal velocity of the 50% glacier-sliding/50% bed-deformation mode has to be larger by approximately one-quarter than that of the 100% glacier-sliding mode. However, because half of this increased basal velocity is due to bed deformation, the sliding velocity is effectively reduced by about one-third. As a result, by neglecting sediment deformation in our analysis we are likely to overestimate the basal sliding rate.

We compared our velocity estimate with direct measurements of sliding using a drag spool. This device is suspended within the borehole close to the base of the glacier, and continuously measures the length of string paid out to an anchor in the bed (Reference Bracewell,Blake and others, 1994). Thus, we have defined “sliding” in an operational manner, as the motion between the spool and the anchor. We recognize that this definition is not equivalent to that of “true” basal sliding and that measurements taken with a drag spool will place only an upper limit on glacier sliding, because the anchor is placed within deformable sediment; any deformation of the sediment lying between the anchor and the ice-sediment contact will introduce additional string extension. With an anchor insertion depth of ∼0.18 m, a typical thickness of the deformable layer underneath Trapridge Glacier of ∼0.5 m (Reference Blake,, Fischer, and Clarke,Blake and Clarke, 1989) and the assumption that the sediment velocity varies linearly with depth, about 5–10% of the measured displacement is due to sediment deformation (Reference Fischer,Fischer, 1995). Hence, both methods of estimating the glacier-sliding rate yield values that tend to be somewhat on the high side.

Complications to our analysis arise if the bed granulometry is unknown and the clast-collision frequency is observed to vary with time. In this case, two contrasting interpretations which can be regarded as end-member cases present themselves: (1) the sliding rate is constant and bed granulometry spatially inhomogeneous; (2) the bed granulometry is spatially homogeneous and the sliding rate varies temporally. In the latter case, sliding rate is proportional to the clast-collision frequency. Suppose that sliding initially occurs at some constant rate U. If at some subsequent time the sliding velocity changed to some new but constant value U′ = αU, the time axis of the ploughmeter response would become rescaled as t′ = t/α (see Equations (10) and (11)). As a result, the new response has a power spectrum α² S(αf) (Equation (12)). Thus, the relative change in sliding rate α = U′/U produces a proportional rescaling of the frequency axis. It may therefore prove possible to use ploughmeters to observe seasonal variations in the basal sliding rate.