Time-dependent basal stress conditions beneath Black Rapids Glacier, Alaska, USA, inferred from measurements of ice deformation and surface motion

Abstract

Observations of surface motion and ice deformation from 2002–03 were used to infer mean stress fields in a cross-section of Black Rapids Glacier, Alaska, USA, over seasonal timescales. Basal shear stresses in a well-defined zone north of the center line (orographic left) were approximately 7% and 16% lower in spring and summer, respectively, than in winter. Correspondingly higher stresses were found near the margins. These changes in the basal shear stress distribution were sufficiently large to cause mean surface velocities to be 1.2 and 1.5 times larger in spring and summer than in winter. These results were inferred with a simple inverse finite-element flow model that can successfully reproduce bulk surface velocities and tiltmeter data. Stress redistribution between the well-defined zone and the margins may also occur over much shorter time periods as a result of rapidly changing basal conditions (ice–bed decoupling or enhanced till deformation), thereby causing large variations in surface velocity and strongly influencing the glacier’s net motion during summer.

1. Introduction

The interaction between a glacier and its bed is an important aspect of glacier dynamics, but remains poorly understood. For glaciers underlain by till, basal motion occurs primarily by sliding along the ice–till interface and/or deformation of the till. Both basal sliding and sediment deformation are facilitated by high basal water pressures.

Experiments conducted from a tunnel beneath Engabreen, Norway, demonstrated that, although elevated water pressure helps to initiate till deformation, the ice can decouple from the till when high water pressure is sustained for a period of a few hours (Iverson and others, 2003). Both till failure and ice–till decoupling locally limit basal shear stresses; these stresses must be redistributed if the glacier is to remain at stress equilibrium (Truffer and others, 2001).

Non-uniform basal motion, which can be considered evidence for stress redistribution, was first observed by Raymond (1971) during a study of borehole inclinometry in a cross-section of Athabasca Glacier, Canada. Kavanaugh and Clarke (2001) and Mair and others (2003) later demonstrated that basal stresses can be rapidly redistributed during short-lived motion events in summer. However, neither of these studies quantitatively described stresses along the glacier bed, nor did they indicate how much these events influence the net motion of the glacier. In this study we expand on the conclusions of Kavanaugh and Clarke (2001), Truffer and others (2001) and Mair and others (2003) by inferring time-dependent basal stresses in a cross-section of Black Rapids Glacier, Alaska, USA, from observations of surface motion and ice deformation during 2002–03. To this end we use a simple inverse approach with a finite-element (FE) model to calculate ice velocities. Basal stresses are computed from the modeled velocity gradients.

This paper begins with descriptions of the field methods and data collected, followed by a summary of the modeling effort. The model results, which are based on the data, are presented and discussed in the context of till mechanics and ice–till coupling.

2. Field Setting and Methods

Black Rapids Glacier is a 40 km long surge-type glacier located in the central Alaska Range (Fig. 1). It last surged in 1936–37 and, based on the current position of loop moraines and on surge cycles of similar glaciers in the region, is thought to have a surge cycle of approximately 75–100 years. Currently there are no signs of an impending surge. It lies on the Denali Fault, a major tectonic feature (Post, 1969), along which a magnitude 7.9 earthquake occurred on 3 November 2002.

Fig. 1. Map of Black Rapids Glacier showing distance from the headwall (circles), the locations of the drilling transect (diamonds), the GPS (global positioning system) base station (small star northwest of the drilling transect) and the Denali Fault (dashed line).

The glacier has been studied since 1970 by the US Geological Survey, the University of Washington and the University of Alaska Fairbanks (e.g. Heinrichs and others, 1996). Annual average velocities and mass balance have been measured throughout this period. A seismic study (Nolan and Echelmeyer, 1999a, b) and two drilling projects (Truffer and others, 1999; Truffer and Harrison, 2006) have revealed a till layer beneath the glacier that is locally up to 7 m thick. These studies focused on a region of the glacier located 16 km from the glacier headwall, where the magnitude and variation of seasonal and annual average velocities are the highest on the glacier (Heinrichs and others, 1996).

In May 2002, we deployed instruments in six boreholes along a roughly north-to-south transverse cross-section (Figs 1 and 2) located near the previous study sites. The boreholes, which were drilled to a diameter of approximately 15 cm, were equipped with three to five instrument packages (48 cm long, 9 cm diameter stainless-steel or aluminum casings) suspended from a shared cable. Each package contained a microprocessor (BasicStamp 2, Parallax Inc.) and one dual-axis electrolytic tiltmeter with a signal conditioning board (Frederiks Co). Those near the bottom were also equipped with a pressure transducer (MSI) and a few had three-axis magnetometers (Precision Navigation) to help resolve tiltmeter orientation. The microprocessors directly measured pulse-width signals output by the tiltmeters; pressure transducers and magnetometers were read by the microprocessors through analog-to-digital converters and serial ports, respectively. Measurements were made every 6 hours. A data logger on the surface switched a relay to provide power to all instruments in the borehole via a common power line. The processors sequentially turned on, computed the average of five measurements from each instrument, digitized and added an identification number to the signal, and transmitted the digitized information to the surface via a common data line. A Campbell 21X data logger recorded the signal and stored it in a storage module. At one borehole (S1) the storage module was corrupted, and we lost several months of data.

Fig. 2. A cross-section profile of Black Rapids Glacier at the drilling transect. Dotted lines indicate boreholes, and crosses indicate tiltmeter locations. At least one pressure transducer was installed in each borehole.

The data collection and transmission methods were very reliable. Signal digitization within the instrument packages ensured that measurements were not affected by issues like data logger temperature variations or cable stretching. We can be confident that, when data were transmitted, they were transmitted accurately.

The nominal range of the tiltmeters is ±15°; they were calibrated in the laboratory to ±20° and appear to be well behaved over this extended range. The resolution of the tiltmeter readings is better than 0.0001° over the range of ±20^°, though the accuracy, which depends on the accuracy of the instruments used for calibration, is closer to ±0.1°. The pressure transducers were calibrated in the field as the instruments were lowered into the boreholes; the resolution of the transducer readings, which was limited by the analog-to-digital converter, was 0.22 m.

All instrument records that contain less than 3 weeks of data have been removed from this analysis since we are primarily investigating seasonal changes in the flow field. Unfortunately, many instruments, especially those located near the bed, were lost due to water leakage within the first few weeks of deployment. This includes most of the pressure transducers and magnetometers. Furthermore, the few magnetometer data that were collected show no temporal variation in the strength of the magnetic field components (along the magnetometers' axes);deformation of the ice should change the orientation of the magnetometers, thus affecting the magnetic field readings along the magnetometers' axes. One explanation for this lack of variation is that the magnetometers were adversely affected by the electronics in the instrument packages. All of the magnetometer records were neglected in the analysis.

To facilitate borehole closure and coupling of the instrument packages to the ice, we pumped some water out of the boreholes before inserting the instruments. The borehole tops were packed with snow to reduce the risk of water draining through the boreholes and affecting the instrument records. Although it is difficult to confidently comment on borehole closure rates, synchronous diurnal fluctuations in tilt observed by all of the tiltmeters and the correlation of some short-lived tilt events suggest that the instrument packages were well coupled to the ice shortly after deployment.

The borehole positions were surveyed with GPS (global positioning system) equipment during the drilling campaign in April and May 2002, and again in September 2002 and May 2003. Additionally, a GPS station was placed on the ice surface near the N1 borehole to record position twice daily throughout the summer (25 April–7 August 2002; days 115–219). The GPS data were differentially corrected to a base station located on the valley wall approximately 2 km from the drilling transect (Fig. 1). The error in the GPS measurements is about ±0.005 m, as previously estimated by a study of the baseline between two benchmarks located near the study site. Tilting of the GPS station at the N1 borehole also contributed some error to the twice-daily velocity measurements. In calculating the error in those velocities we therefore conservatively estimate that the uncertainty in the antenna position is ±0.01 m.

We supplemented the measurements of ice motion and water pressure with daily average air temperature at Gulkana Glacier, available from the US Geological Survey (Fig. 3; R. March, unpublished data). The Gulkana Glacier meteorological station is about 200 m lower and 60 km east of our study site. Although the station is located on the south flank of the Alaska Range and may observe different weather patterns than Black Rapids Glacier, it is the closest station and the only one nearby located at a similar elevation.

Fig. 3. Daily average air temperature at the Gulkana Glacier meteorological station (R. March, unpublished data). The dotted curve represents the 40 year daily average temperatures.

3. Observations

3.1. Ice deformation

Tiltmeters measure tilt from horizontal for two mutually perpendicular axes. The tiltmeter data are transformed from the tiltmeter frame to a map frame with three rotations (Fig. 4), enabling calculation of tilt from vertical, θ, and rotation angle, ψ (Blake, 1992). It is impossible to resolve the azimuth, ϕ, from the tiltmeter data alone. Magnetometers were employed to determine this angle but they did not function properly (see section 2). Tiltmeter records are labeled by borehole name and height above the bed (in meters).

Fig. 4. Tiltmeter orientation in a map frame: are longitudinal, transverse and vertical up directions, respectively. The angles θ and ψ can be resolved by rotating the tiltmeter from the coordinate system to the coordinate system.

The tiltmeter records generally indicate steady deformation over periods of weeks to a year (Fig. 5). Superimposed on these general trends are large fluctuations and/or steps in tilt angle. Small gaps in the records are due to failure during data transmission; large gaps are the result of the tiltmeter turning off and restarting at a later date. This could be due to a loss of battery power, though we are unable to provide satisfactory explanations for many of the data gaps.

Fig. 5. Tilt angle as a function of time. The solid curves represent tiltmeter data; the dashed curves are model-derived synthetic tilt curves (see section 4.1). Dotted curves indicate time periods during which data were omitted during generation of the synthetic tilt curves. Tiltmeter records are labeled by borehole name and height above the bed (in meters).

N2-11 (Fig. 5a) was the only tiltmeter in the N2 borehole to operate for more than a few days. Its tilt rates were considerably higher than those observed by other tiltmeters. The negative initial slope of the tilt curve suggests that the tiltmeter was initially tilted up-glacier; vertical shearing probably caused the top of the tiltmeter to move towards vertical and the tiltmeter to achieve a minimum tilt angle as it passed through the vertical transverse plane.

Two tiltmeter records from the N1 borehole exceeded 100 days. N1-51 (Fig. 5b) experienced steady tilt rates throughout the year, with slightly higher rates in winter than in summer. N1-11 (Fig. 5c) exhibited very low tilt rates during summer, punctuated by steps in tilt of up to 1.5°.

Only limited tilt data are available from the CEN borehole. The short and sparse record from CEN-6.5 (Fig. 5d) indicates a decrease in tilt rate as summer progressed.

Tiltmeters in the S1 borehole exhibited highly variable behavior. S1-51 (Fig. 5e) experienced large fluctuations in tilt angle during early summer, including two large events on days 154 and 180. The tilt rate was constant between days 192 and 241. The S1-6 tilt record (Fig. 5f) spanned days 148–252. In general, its tilt angle steadily decreased during summer; superimposed on this trend is a positive departure between days 164 and 194. S1-6 was presumably tilted up-glacier. Two tiltmeter records from the S1 borehole were neglected from this analysis because their tilt angles often greatly exceeded design specifications.

Tilt data from the S2 borehole are also very limited: only S2-6 (Fig. 5g) yielded reliable data. It demonstrated steady tilt rates throughout summer, with a large drop in tilt on day 180. The tilt rates before and after this event were nearly identical.

Tiltmeters in the S3 borehole (Fig. 5h–k) provided the longest and most complete records. These tiltmeters were located farther from the bed than the tiltmeters in the other boreholes because they became stuck during deployment, as indicated by a drop in cable tension and a stabilization of the water-pressure readings at an elevation of 70.6 m above the bed. The uppermost tiltmeter, S3-170.6 (Fig. 5h), exhibited large tilt variations in summer. During winter (days 271–426), its tilt rate remained steady and maintained a slope similar to the mean slope of the summer data. In late winter its tilt rate rapidly increased. S3-120.6 (Fig. 5i) exhibited the steadiest deformation of all the tiltmeters. The tilt angle of S3-75.6 (Fig. 5j) oscillated during the first few weeks of summer and then slowly and steadily decreased. The lowest tiltmeter, S3-70.6 (Fig. 5k), demonstrated slow, steady changes in tilt throughout the study period, except for a large departure from the general trend between days 172 and 200. Tiltmeters S3-75.6 and S3-70.6 were, presumably, initially tilted up-glacier.

Finally, all of the tiltmeters exhibited small diurnal fluctuations in tilt during summer (Fig. 6) that cannot be seen at the scale of an entire tilt curve. These fluctuations are very similar to diurnal tilt fluctuations that have been observed by tiltmeters installed in till beneath other glaciers (Iverson and others, 1999; Kavanaugh and Clarke, 2006). In our measurements the maximum and minimum daily tilt typically occur at 0600 h and 1800 h, respectively, throughout the summer, though it should be noted that measurements were only obtained every 6 hours. The lack of drift in the timing of the maximum and minimum daily tilt indicates that the fluctuations cannot be attributed to the principal components of Earth tides. Furthermore, the fluctuations cease abruptly on or before day 270 (27 September 2002).

Fig. 6. Examples of diurnal fluctuations in tilt angle from three boreholes. All tilt curves showed in-phase diurnal variations during summer that stopped abruptly on or before day 270.

3.2. Surface velocity

The mean spring (days 120–136, 30 April to 16 May 2002) and summer (days 136–257, 16 May to 14 September 2002) surface velocities were approximately 1.2 and 1.5 times larger, respectively, than the mean winter (days 257–490, 14 September 2002 to 5 May 2003) surface velocities, across the drilling transect (Fig. 7). This is typical of velocities on Black Rapids Glacier (Heinrichs and others, 1996).

Fig. 7. Mean spring (30 April to 16 May 2002; squares), summer (16 May to 14 September 2002; circles) and winter (14 September 2002 to 5 May 2003; triangles) surface velocities across the drilling transect. The curves are model results (see section 4.1) for various relative contributions of the velocity and tiltmeter root-mean-square errors; the velocity data are not adequately reproduced when the inverse model depends only on tiltmeter data (not shown). There is only one curve for the spring velocity profile because no tiltmeter data were available for that time period, and so the inverse model depended solely on velocity data. The thick curves represent the solutions we used in our analysis: the velocity and tiltmeter errors were weighted by 0.25 and 0.75, respectively. The error in the velocity measurements is smaller than the line thickness.

Daily surface velocities near the N1 borehole during summer are presented in Figure 8. Velocities were between 0.15 and 0.2 m d⁻¹ (55 and 73 m a⁻¹) prior to the annual spring speed-up, which began on day 140. The maximum spring speed-up velocity of 0.55 m d⁻¹ (201 m a⁻¹) was reached on day 146. Velocity peaks were either small but long-lived (days 146, 153, 160 and 170) or large but shortlived (days 180, 196, 201 and 214). Similar velocity variations during summer have been observed previously on Black Rapids Glacier (Nolan and Echelmeyer, 1999a; Truffer and others, 2001) and elsewhere (e.g. Willis, 1995). We do not know vertical velocities during this period because the GPS station was resting on the ice surface.

Fig. 8. Twelve-hour mean surface velocity (solid curve) compared to water level measured in the till (dotted curve) at the N1 borehole during summer. Velocities are given in m a⁻¹ for comparison with model results.

3.3. Water pressure

Although pressure transducers were installed in each borehole, records are only available from two boreholes. Data from this study were therefore supplemented with water-pressure measurements from a pressure transducer installed in the till near the N1 borehole (Truffer and Harrison, 2005).

The N1-till record (Fig. 9a) is lengthy, but it is also sparse due to wireless data transfer methods (Harrison and others, 2004). Major peaks in water pressure were observed on days 146, 198 and 249, and a large drop in pressure occurred around day 175. Owing to the sparsity of the dataset, it is difficult to comment on magnitudes and rates of water-pressure fluctuations.

Fig. 9. Piezometric surface for (a) N1, (b) S1 and (c) S3 boreholes. The N1 pressure data were obtained from two pressure transducers installed in the till near the N1 borehole. The dashed lines represent the ice-overburden pressure at each borehole. The data in (a) were collected at irregular intervals due to wireless data transfer methods (Harrison and others, 2004).

The S1 record (Fig. 9b) lasted from day 148 to day 170. Water pressure was generally high during this period, with major drops in pressure occurring on days 149, 156 and 163. Diurnal fluctuations in water pressure began on day 165.

The most complete record is from the S3 borehole (Fig. 9c). It indicates steady water pressure until day 160, followed by large diurnal fluctuations throughout the summer. Large drops in water pressure occurred on days 163 and 173; major peaks were observed on days 165, 170 and 180. Water pressure was generally low between days 190 and 200; in late summer it gradually increased until it reached a nearly constant value that was sustained throughout the winter.

4. Modeling Approach

We obtained velocity and tiltmeter data in an attempt to infer temporal and spatial variations in the basal shear stresses beneath Black Rapids Glacier. The water-pressure data were collected to help interpret the changes in basal conditions. However, the velocity and tiltmeter records, when treated individually, can be explained by a variety of stress fields. Gudmundsson and others (1999) addressed this issue by acknowledging that velocity and tiltmeter data obtained concurrently are the consequence of the same flow field. By parameterizing the velocity field and adjusting its parameters, they were able to determine the velocity profile (including basal motion) that was best able to reproduce the surface velocity and multiple tiltmeter datasets from a borehole in Unteraargletscher, Switzerland. We adapt their approach to a two-dimensional (2-D) finite-element flow model; the model solution includes the velocity field and its gradients, which are used to determine the stress field and basal shear stresses.

4.1. Model description

Ice flow is modeled through a transverse cross-section (Fig. 2) 16 km from the glacier headwall (Fig. 1). The model is simplified by assuming that there is (1) no extension or compression along the flowline and (2) no vertical or transverse flow. The first assumption appears valid over seasonal timescales, during which longitudinal strain rates at the glacier surface are on the order of 10⁻³ a⁻¹ (Nolan, 2003). This turns out to be more than one order of magnitude smaller than the shear strain rates estimated by our FE model at the location of the uppermost tiltmeter (S3–170.6) and should therefore have little effect on the model solutions, especially due to the non-linearity of the flow law (see Equation (2) below). The second assumption requires a flat transverse surface profile, a straight channel and little influence from tributary glaciers, all of which are reasonable approximations for this site. These assumptions reduce the ice-flow equations to a non-linear Poisson equation:

(1)

where η is the stress-dependent viscosity, ρ is ice density, α is the mean out-of-plane surface slope, g is acceleration due to gravity, u is the out-of-plane velocity, y is the transverse direction and z points upward and is perpendicular to the mean glacier slope. A power-law rheology (Glen’s flow law) is used to describe the non-linear viscosity of ice:

(2)

where A is the flow-law parameter, n is an empirical constant, is the second invariant of the strain-rate tensor and = 10⁻¹⁵ a⁻² is a finite viscosity parameter used to prevent infinite viscosity at low stresses. We set α = 1.8°, A = 0.1 a⁻¹ bar⁻³ (3.17 × 10⁻²⁴ s⁻¹ Pa⁻³) and n = 3. The surface slope, α, was estimated from GPS measurements. The flow-law parameter, A, was determined by adjusting its value until a flow model with no sliding along the bed produced the lowest measured velocities at this site during the past 30 years (Heinrichs and others, 1996).

The model domain is derived from a radio-echo sounding (RES) profile (Gades, 1998) and borehole depths. The errors in the RES data are on the order of ±5 m, but the close agreement between our borehole depths, known to better than ±0.1 m, and the RES data suggests that the bed profile is known to better than ±5 m. The surface is assumed to be stress-free, which results in a Neumann boundary condition ( = 0, where is the outward normal vector), and we assume a Dirichlet condition (u = f(y)) for the bed.

We seek a basal velocity function, f(y), that minimizes the error between model results (velocities and model-derived synthetic tilt curves, described below) and measured surface velocities and tiltmeter data. In reality, basal velocity is a function of the stress field, bed roughness and other physical characteristics of the bed, which would require a mixed (Robin) boundary condition or, in the case of till deformation, a stress boundary condition. We adopt the Dirichlet boundary condition in order to derive the basal velocity distribution, and ultimately the basal stress distribution, without attempting to make a conclusive statement about the nature of the actual boundary condition. Note that, for the purposes of an inversion, it is irrelevant which basal boundary condition is assumed. Existence and uniqueness of the forward problem (i.e. solving Equations (1) and (2) with the stated boundary conditions) were proven by Colinge and Rappaz (1999); thus, the velocity field that best reproduces all the data is uniquely determined regardless of the boundary condition implemented in the model. We therefore choose to use a Dirichlet boundary condition along the bed because it provides for fast convergence and is easier to implement than mixed or stress boundary conditions.

Synthetic tilt curves are generated according to the procedures outlined in Gudmundsson and others (1999). The modeled velocity gradients (∂u/∂y and ∂u/∂z), initial tilt angles, θ _o, and initial azimuths, ϕ ₀, are required to generate synthetic tilt curves. We use the initial tilt angle from a given tiltmeter record and the modeled velocity gradients at the tiltmeter’s location, and iteratively search for the value of ϕ _o that minimizes the root-mean-square (rms) error between the synthetic and measured tilt curves.

Determining the basal velocity distribution from measurements of ice deformation and surface velocity is an inverse problem. A basal velocity distribution can be found that reproduces the data exactly; potentially there are many such solutions. Due to errors in the data and uncertainties in the physical model, a solution that reproduces the data exactly would probably be unrealistic. In general inverse theory, a preferred solution is one that fits the data to within a given error and also minimizes some undesired quality of the unknown boundary condition (Parker, 1994). In the problem at hand, a reasonable solution might be one that forces the basal velocity distribution to be a smooth function with small first derivatives (Truffer, 2004);the data are then fitted to within a given error by minimizing a norm that reflects the smoothness of the function.

Applying general inverse theory to a highly non-linear, 2-D FE model with an irregular boundary is non-trivial and computationally expensive. To simplify the problem, we seek a smooth basal velocity function that can be approximated by a fourth-order polynomial, with the additional assumption that the velocity equals zero at the margins. Although this choice of basal velocity function is somewhat arbitrary, it appears to be the simplest that can satisfactorily reproduce our measurements without restricting the location or magnitude of basal velocity maxima and minima. (We also tried higher-order polynomials and found the results to be very similar.) We are forced to specify the velocity at the margins because the actual velocity there is unknown, and leaving it unspecified results in very unrealistic solutions. However, the velocity specified at the margins appears to have very little effect on the model solution, especially near the glacier center line (see appendix A in Amundson, 2006). The basal shear stress distribution, which we ultimately desire, can be computed from the model solution.

Ideally the flow field would be determined as a continuous function of time. This is impossible with our velocity data because they are of insufficient temporal resolution to constrain the inverse model over short time periods. Furthermore, all our measurements were restricted to a transverse cross-section, thus forcing us to neglect longitudinal strain. This approximation is adequate for investigating strain rates over periods of months or longer in our study area (discussed above). However, over shorter time periods such as during the spring speed-up (Nolan, 2003) and in summer (when velocity variations are large), longitudinal strain may be important. Our model also does not take into account viscoelastic effects that may be associated with short-lived motion events.

The data were therefore divided into spring (days 120–136, 30 April to 16 May 2002), summer (days 136–257, 16 May to 14 September 2002) and winter (days 257–490, 14 September 2002 to 5 May 2003), corresponding to trips to the field that occurred in spring and autumn 2002, and spring 2003. This division roughly agrees with changes from above-freezing to below-freezing air temperatures at nearby Gulkana Glacier (Fig. 3). In this manner we limit our investigations to seasonal changes in the mean stress field. Acknowledging these limitations, we will, however, qualitatively discuss short-term stress variations by considering tilt deviation from the synthetic tilt curves.

The coefficients of the seasonal basal velocity functions were determined simultaneously using a multidimensional unconstrained non-linear optimization (Nelder–Mead method, implemented with Matlab’s fminsearch) to minimize the error between data and model results. We coupled the summer and winter models by maintaining continuity in the tilt angle and azimuth of the synthetic tilt curves at the summer/winter transition; the spring model was run separately because it was based entirely on surface velocity data.

The modeling approach is summarized as follows:

1. 1. Guess the coefficients of the summer and winter basal velocity functions. Only three coefficients are required for each, since we are requiring the basal velocity functions to be fourth-order polynomials with zero velocity at the margins.

2. 2. Solve the spring, summer and winter FE models. Export strain rates and surface velocities.

3. 3. Calculate the percentage rms error between modeled surface velocities and observations (later referred to as the velocity error).

4. 4. For the summer and winter models, generate synthetic tilt curves for each tiltmeter record using the modeled velocity gradients and compute the percentage rms error between synthetic and measured tilt curves. The percentage rms error of each tiltmeter is weighted by the number of data points in the respective set to give the longer records more importance. Sum the weighted rms errors and divide by the total number of data points in all the records (later referred to as the tiltmeter error).

5. 5. Compute the total percentage error using some combination of the velocity and tiltmeter errors. Since it is not clear how to do this, we considered several possibilities (see section 4.2).

6. 6. Iteratively search for the coefficients of the basal velocity functions (summer and winter) that minimize the total percentage error using the Nelder–Mead method.

4.2. Model results

The model captured the seasonal tilt rates and surface velocities well (see Figs 5 and 7), although the agreement between synthetic and measured tilt curves could partly be attributed to the unknown initial azimuth of the synthetic tilt curves. Modeled velocity fields are shown in Figure 10a–c; the corresponding octahedral stress fields were determined from the modeled velocity gradients and are shown in Figure 10d–f. The results suggest that basal shear stresses in a zone about 500 m north of the deepest point in the channel were approximately 0.07 bar and 0.16 bar (7% and 16%) lower in spring and summer, respectively, than in winter (Fig. 11), and that basal shear stresses near the margins were correspondingly higher in summer than in winter (thus allowing the glacier to remain in stress equilibrium). This stress redistribution is qualitatively similar to that predicted from forward models assuming a Coulomb-friction rheology for the till (Truffer and others, 2001). However, it should be noted that the ‘spring’ period was a short period that immediately preceded the spring speed-up, an event that is typically associated with longitudinal extension (e.g. Nolan, 2003). Thus, the assumption of no longitudinal strain may be invalid during this period.

Fig. 10. Model results (shown in order of increasing velocity): (a–c) mean winter, spring and summer velocity fields (m a⁻¹); (d–f) mean winter, spring and summer octahedral stress fields (×100 kPa). Flow is into the plane.

Fig. 11. Stress difference plots. The winter octahedral stress field is subtracted from (a) the spring octahedral stress field and (b) the summer octahedral stress field. The shaded regions indicate areas where stresses were lower in spring and summer than in winter. The values given on the contours are in units of 100 kPa; note that the plots have different contour intervals.

The model results were largely independent of the relative weights assigned to the velocity and tiltmeter errors in the inverse model (see Fig. 7), although surface velocities had to be included in the analysis to ensure a good fit to all the data. Removing tiltmeter data from the analysis did not significantly affect the solutions. This provides tentative support for inverse models based solely on surface data, although in our case the agreement between the synthetic and measured tilt curves greatly increases our confidence in the model results. Hereafter we have weighted the velocity error by 0.25 and the tiltmeter error by 0.75, as this provided the best fit for both velocity and tiltmeter data.

Solution uniqueness cannot be proven for most nonlinear inverse problems (Parker, 1994); it is therefore difficult to determine whether our model solutions represent global or local error minima. However, the coefficients of the basal velocity functions are at least locally well constrained (Amundson, 2006).

5. Discussion

5.1. Rapid motion events

Seasonal velocity variations at Black Rapids Glacier appear to be the result of basal stress redistribution between a well-defined region north of the center line and the margins. This well-defined region coincides with the area where Nolan and Echelmeyer (1999a, b) observed the till becoming seismically transparent during lake drainages and associated rapid motion events. The bed there seems to be highly susceptible to changes in water flux, maybe due to the several lakes along the northern margin of the glacier that might drain subglacially and significantly influence the subglacial drainage system.

The influence of the drainage events on the basal stresses can be investigated with our FE model by making a few assumptions:

1. 1. The velocity spike on day 180 (Fig. 8) can be attributed to a lake drainage event since it is similar in magnitude and duration to the velocity events that Nolan and Echelmeyer (1999a) and Truffer and others (2001) related to lake drainages. These events occur every summer.

2. 2. The velocity prior to day 180 is very close to the mean velocity during the following winter, and therefore the basal velocity immediately preceding the motion event can be approximated by the previously determined winter basal velocity distribution.

3. 3. Water from lake drainages is routed through the region near the N1 borehole.

4. 4. The water influx from the drainage events is sufficiently large to locally decrease the normal stresses imparted to the bed by the glacier (Nolan and Echelmeyer, 1999b). This results in ice–bed decoupling, which can occur at pressures beneath the ice overburden pressure (e.g. Iverson and others, 1999, 2003).

Thus, we numerically decouple the ice from the bed near the N1 borehole (by setting the local basal shear stresses equal to zero) and specify the basal velocity along the rest of the bed according to the previously determined winter basal velocity distribution. By adjusting the length and position of the decoupled zone until the modeled velocity at the N1 borehole agrees with the velocity observed on day 180, we can estimate the basal shear stress distribution during a lake drainage event.

Under these assumptions, the section of the boundary that must be decoupled in order to achieve the velocity on day 180 is very close to the length of the bed over which mean basal shear stresses are lower in spring and summer than in winter (Fig. 12). The same surface velocity can be obtained by shifting the region of the bed that is decoupled from the ice by a few hundred meters in either direction; increasing or decreasing its length by more than 100 m significantly reduces the agreement between the modeled velocity at the N1 borehole and our observations. Although this model was not well constrained (we used a surface velocity measurement at just one point) and did not take into account additional till deformation or basal sliding that could result from increasing the shear stresses adjacent to the decoupled zone (Mair and others, 2003), it does reinforce the contention that large increases in surface velocity (over a variety of timescales) are due to stress transfer from a well-defined zone near the N1 borehole toward the margins.

Fig. 12. Model-derived mean summer and winter basal shear stresses (dashed and dotted curves, respectively) and a possible basal shear stress distribution of a decoupling event (solid curve). The spring basal shear stress distribution is not shown because it is difficult to see the difference between the winter and spring stress distributions at this scale. The negative shear stresses and the large jump in stress observed at the zero basal shear stress/specified velocity transition are due to a mathematical singularity (Hutter and Olunloyo, 1980).

Observations of till deformation (Truffer and Harrison, 2006) agree with this model in that little till deformation appears to have occurred near the N1 borehole during the short-lived motion events. However, enhanced till deformation, which occurred almost exclusively during the spring speed-up, would also have produced stress redistribution toward the margins (Truffer and others, 2001). The spring speed-up and the short-lived motion events that we attribute to lake drainages (Fig. 8, days 142–163, 180, 196, 201 and 214) significantly influenced the mean velocities during summer; the mean velocity at the N1 borehole during these periods was 192 m a⁻¹, while the mean velocity during the rest of the summer was 46 m a⁻¹, about 7 m a⁻¹ slower than the mean velocity during winter 2002/03. These rapid motion events, which we suggest are caused by rapid stress redistribution toward the margins, are thus responsible for the overall larger ice displacement in summer than in winter.

5.2. Diurnal tilt fluctuations

Tiltmeters react to a strain field, which cannot be determined uniquely from the tiltmeter measurements without additional data or assumptions. It is thus difficult to quantitatively interpret the diurnal fluctuations in tilt angle observed by all of the tiltmeters (Fig. 6), as we are unable to neglect any of the strain-rate components over short timescales (see section 4.1). By considering the effect of each strain-rate component on the tiltmeter readings we can, however, attempt to qualitatively explain the diurnal fluctuations in tilt. For example, diurnal fluctuations in shearing across the glacier would cause some tiltmeters to achieve maximum daily tilt at the same time that others are achieving minimum daily tilt, except in the very unlikely event that all tiltmeters were oriented to the same side of the vertical longitudinal plane. Thus, shearing across the glacier is an unlikely mechanism for the diurnal tilt fluctuations. Fluctuations in shearing along the flowline, for which we have no information, could partly account for the diurnal tilt fluctuations. However, the negative tilt rates during the diurnal tilt fluctuations (for tiltmeters oriented down-glacier) cannot be explained by diurnal fluctuations in shearing along the flowline, as that would imply up-glacier deformation. This suggests that longitudinal strain is also important.

The diurnal tilt fluctuations are consistent with observations and a mechanism proposed by Sugiyama and Gudmundsson (2003). They obtained detailed observations of vertical strain in Unteraargletscher, from which they concluded that increased water influx to the bed during the day causes the glacier to accelerate; the acceleration is amplified up-glacier because the drainage efficiency is lower there, thus causing longitudinal compression. At night, when water influx decreases, the longitudinal flow speeds become more uniform and motion is dominated by shearing. This mechanism could cause tiltmeters to achieve maximum daily tilt in the morning and minimum daily tilt in the evening, similar to our observations (Fig. 6). However, the diurnal increases in tilt angle exhibited by tiltmeters oriented up-glacier (e.g. Fig. 6b) suggest that some extension occurred at night. This does not contradict the mechanism proposed by Sugiyama and Gudmundsson (2003); compression could result in vertical thickening that is larger down-glacier, thus causing subsequent extension. Further supporting this mechanism on Black Rapids Glacier is the observation that the diurnal fluctuations in tilt angle abruptly cease on or before day 270, which probably coincides with the onset of winter conditions (Fig. 3) and decreased water influx.

5.3. Non-steady deformation

In our model we have assumed that ice deforms as a continuous medium. However, most tiltmeters experienced short-lived tilt events that are probably indicative of nonsteady deformation. Some of these events were localized and had little or no correlation with other tiltmeter, velocity or water-pressure records. For example, the large fluctuations in tilt angle during summer observed by tiltmeters N1-11, S1-51 and S3-170.6 (Fig. 5c, e and h) are poorly correlated with any of the instrument records. There were, however, at least two instances in which tilt events were well correlated between two instruments located in different boreholes but were not observed by more proximal tiltmeters: the large changes in tilt exhibited by S1–51 and S2-6 on day 180 (Fig. 5e and g) and the slope reversals between days 170 and 180 in the S1-6 and S3-70.6 tilt curves (Fig. 5f and k). Although the lack of correlation between most tilt events could be due to poor coupling between the tiltmeters and the ice, the strong synchronicity of the diurnal fluctuations in tilt (Fig. 6) and the strong correlation between the S1-6 and S3-70.6 tilt curves suggest that the tiltmeters were well coupled to the ice.

Other tiltmeter records exhibit similar behavior that cannot be explained by steady-state flow models. For example, tiltmeters installed near Jakobshavn Isbræ, Greenland, indicated overthrusting at depth (Lüthi and others, 2003), and some (but not all) tiltmeters installed in a single borehole in Unteraargletscher, demonstrated quasi-repetitive tilt oscillations (Gudmundsson and others, 1999). More knowledge of small-scale glacier mechanics and the coupling of tiltmeters to ice will be required if we are to be able to interpret all the fluctuations in tilt rate observed by a tiltmeter. Also, it remains unclear how much non-steady deformation, such as faulting, influences the net motion of a glacier.

6. Conclusions

A well-defined zone approximately 500 m north of the deepest point in the channel of Black Rapids Glacier experienced mean basal shear stresses in spring and summer 2002 that were 7% and 16% lower, respectively, than those observed during winter 2002/03; similarly higher shear stresses were found near the margins. The seasonal changes in the basal shear stress distribution were sufficient to cause mean surface velocities in spring and summer to be approximately 1.2 and 1.5 times greater than the mean surface velocities during the following winter. These results were inferred with an inverse FE model that incorporates measurements of ice deformation and surface velocity from a transverse cross-section of the glacier.

The well-defined zone also corresponds to a region where seismic anomalies were observed during the drainage of marginal lakes (Nolan and Echelmeyer, 1999a, b);the bed in that region may be particularly susceptible to changes in water flux. The high surface velocities associated with the drainage events can be reproduced using the FE model by decoupling the ice from the bed (setting the basal shear stress equal to zero) along this zone. The velocity field during the drainage events, predicted by our model, suggests that the short-lived motion events can have a very large effect on the seasonal stress distributions. This is supported by the observation that surface velocities during most of the summer (excluding the days during which a motion event occurred) were very similar to the mean surface velocities in winter.

Our modeling effort assumed steady-state deformation and did not allow for viscoelastic effects or short-lived fluctuations in basal conditions. Many of our tiltmeter data, as well as observations from other recent studies (e.g. Gudmundsson and others, 1999; Bindschadler and others, 2003; Lüthi and others, 2003; Elsberg and others, 2004; Kavanaugh and Clarke, 2006), demonstrate that glacier stresses can vary over small regions and change on timescales that cannot be addressed under these assumptions. The amount of motion that is due to fracturing and other local phenomena needs to be further investigated in order to increase our understanding of both short- and longterm velocity variations.