Introduction

Mount St Helens is an active volcano with a glacier in its crater. The eruption that began in September 2004 has been largely a quiescent one featuring the growth of a lava-dome complex (Reference Vallance, Schilling, Thompson, Messerich, Sherrod and ScottVallance and others, in press). Early in the eruption sequence, a ‘fin’ of solidified magma broke through the surface of the ∽1km² Crater Glacier and grew southward, intersecting the south crater wall in mid-November 2004 and cleaving Crater Glacier into eastern and western halves (Fig. 1; see also Reference Walder, Schilling, Vallance, Lahusen, Sherrod and ScottWalder and others, 2007a, Reference Walder, Schilling, Vallance and LaHusenb). The fin then expanded eastwards, squeezing the east crater glacier against the east crater wall until mid- to late-April 2005, when dome growth shifted to the west. The response of East Crater Glacier to the squeeze event and its aftermath was recorded by photography, photogrammetry and single-frequency, helicopter-deployed Global Positioning System (GPS) stations on the glacier. East Crater Glacier underwent a hitherto never-described style of deformation: strain rates were of a similar order of magnitude to those observed in glacier surges (Reference KambKamb and others, 1985), but here the causative stresses were oriented across flow, rather than along flow.

Fig. 1. Field setting of Mount St Helens crater glacier. (a) Mount St Helens in October of 2000, looking south. Crater Glacier is partially obscured on the east side of the 1980–86 lava dome by rock-fall debris, and on the west side merges with a rock-laden ice mass shed from the west crater wall (photograph courtesy of S.P. Schilling, USGS). (b) The crater of Mount St Helens on 10 April 2005. View to the south. The dashed line is the approximate centerline of East Crater Glacier. Locations of GPS stations are marked with an ‘x’ (photograph courtesy of J.J. Major, USGS).

As the squeeze event progressed, the thickness of ice in the upper part of East Crater Glacier became 2 to 3 times that in the lower part (Fig. 2), so if deformation were only by simple shear in the vertical, the difference in surface velocity between the upper and lower portions of East Crater Glacier should have been a factor of about 16–80 for a flow-law exponent n = 3 (Reference van der Veenvan der Veen, 1999). The observed velocity difference (Fig. 3), however, was only a factor of 3–4. Reference Walder, LaHusen, Vallance and SchillingWalder and others (2005) suggested an explanation: the lower third of the glacier, which was not squeezed, acted as a dam that pushed back against the ice upstream. This would lead to a large, longitudinal stress gradient. To assess this hypothesis, we use a two-dimensional, full-stress flowband model to simulate evolution of East Crater Glacier.

Fig. 2. Surface elevation profiles along the East Crater Glacier centerline in Figure 1. Line types denote the centerline elevation at various dates during the squeezing event. The fine lines represent the initial computational grid. The approximate location of the velocity measurement station ELE4 is labeled. Labels ‘a’ and ‘b’ mark the approximate upglacier and downglacier bounds of the domain squeezed by the growing lava dome.

Fig. 3. Horizontal speed of East Crater Glacier GPS stations. ICY4 and ICY5 were on the glacier while the lava-dome spine was expanding eastward. ELE4 was fortuitously placed on the glacier about the time that the spine stopped growing. HIE5 was on the glacier in mid-summer. Azimuth of motion for all stations was within 18^˚ of north. For comparison we show surface-speed data (after Reference Anderson, Walder, Anderson, Trabant and FountainAnderson and others (2005)) for a target on Kennicott Glacier, a temperate valley glacier in Alaska, during the year 2000. The record for Kennicott Glacier shows large-amplitude, commonly diurnal fluctuations not seen at East Crater Glacier.

The model is constrained by observations of surface elevation and surface motion, and is tuned by varying not only the stiffness of the ice but also the bulk density. Because independent data provide bounds on the rock-debris content and thus the effective density of the ice, we are able to constrain the rheological properties of temperate glacier ice containing a large fraction of rock debris. In effect, the combination of our flow model and the observations from Mount St Helens in early 2005 serve as a field-scale experiment on the deformation of debris-laden ice.

Field Setting

The Mount St Helens crater glacier formed following the catastrophic eruption of 1980 and the subsequent growth of a lava dome from 1980 to 1986. It is fed by copious avalanching of snow and rock debris from the crater walls. Given these hazardous conditions, our knowledge of the glacier’s development comes almost entirely from occasional visual inspection and photogrammetry. Comparison of sequential digital elevation models (DEMs) revealed that the crater-floor fill accumulated up to mid-1988 was, on average, about 60% rock debris by volume (Reference MillsMills, 1992). A similar exercise by Reference Schilling, Carrara, Thompson and IwatsuboSchilling and others (2004) showed that by September 2000, the average volumetric rock content of the entire glacier was down to 30%. Crevasses – indicators of flow – were first noticed on aerial photographs taken in September 1996. The general picture presented by these studies, along with field observations by an author (JSW) and other US Geological Survey scientists, is of crater-floor fill that grades upward from mostly rock – and clearly not glacier ice in a rheological sense – to ‘dirty’ firn and ice with debris mostly in distinct layers of decimeter-scale thickness.

For purposes of modeling, we wish to consider as ‘glacier’ only the portion of the crater-floor fill that behaves as a creeping solid, and to exclude the lowermost, rock-rich material, which probably behaves as a Coulomb-frictional material with ice in the interstices. Picking the effective glacier bed is admittedly an uncertain exercise. Our choice is to pick the glacier bed as the crater-floor surface defined by DEMs for 12 October and 12 November 1986. The rationale for this choice is presented by Reference Walder, Schilling, Vallance, Lahusen, Sherrod and ScottWalder and others (2007a; Reference Walder, Schilling, Vallance and LaHusen2007b). Average rock content of the glacier so-defined is about 15% by volume. We recognize that the creeping material may in fact extend below the 1986 surface, in which case the average rock content would be greater.

Flow Model

Model description

Our flowband model solves the full two-dimensional momentum-balance equations and includes terms to account for converging and diverging flow. The model is based on the Finite Volume Method (FVM) (Reference PatankarPatankar, 1980; Reference Versteeg and MalalasekeraVersteeg and Malalasekera, 1995) and is described fully by Reference PricePrice (2006). Here, we discuss relevant portions of the model including the governing equations, the general solution method, and the parameterizations that account for the effects of changes in flowband width.

For low-Reynolds-number flow of a viscous fluid, conservation of momentum in a Cartesian reference frame is expressed as

(1)

where x, y and z are the along-flow, across-flow, and vertical coordinate directions and repeat indices imply summation. The first term on the left-hand side of Equation (1) is the body force, the product of ice density ρ and the acceleration due to gravity g _i . The second term is the stress divergence, where σ_ij are components of stress and related to deviatoric stress Ƭ_ij and the mean compressive stress P through

(2)

The constitutive relation is written according to Nye’s generalization of Glen’s flow law,

(3)

(e.g. Reference van der Veenvan der Veen, 1999) where ∊_ij is strain rate,

(4)

ɳ is the effective viscosity,

(5)

and the u_i represent components of the velocity vector (henceforth called u, v and w) parallel to the x, y and z directions, respectively. In Equation (5), E is a scalar enhancement factor (equal to 1 for normal ice), B(T) is a temperature-dependent rate factor (here taken constant), n is the power-law exponent (here taken equal to 3), and ∊_e is the effective-strain rate, given by

(6)

The ice is assumed to be incompressible, in which case

(7)

We solve Equation (1) in a two-dimensional, boundary fitted, orthogonal, curvilinear-coordinate system. The transformation between this and a standard Cartesian coordinate system is discussed in detail by Reference PricePrice (2006). The model solutions discussed below (e.g. velocity fields) have been rotated into a Cartesian coordinate system.

General solution method

Integrating Equation (7) over a single finite volume (i.e. a single grid cell) gives

(8)

where ρ is density. The subscripts U, D, T, B, L and R refer, respectively, to the upstream, downstream, top, bottom, left and right faces of a single cell. A is the area of the relevant cell face: for example, A _U = W _UΔz _U, where W _U and Δz _U are the flowband width and volume height respectively, at the upstream cell face. Equation (8) states that the net mass flux into and out of a volume is equal to zero. A flowband is bounded by two flowlines and if, for the moment, we take the flowband width (and thus the finite-volume width) as constant, the second term in Equation (8) is zero. Using Equations (2)–(6) in Equation (1) with an estimated pressure field and solving computationally, we obtain estimated velocity fields u* and w*. Inserting these values into the remaining terms from Equation (8), we obtain

(9)

where S is non-zero because, in general, our initial estimate for the velocity field will not satisfy continuity. To satisfy continuity the mass source (or sink) S at each volume is eliminated through an iterative pressure-correction method (Reference PatankarPatankar, 1980). A non-zero mass source defines a pressure perturbation that improves upon the estimated pressure and velocity fields. Through Equation (9), the updated velocity field leads to a further improvement in the estimate for the mass source (i.e. one with a smaller magnitude) and a further improvement on the estimated pressure perturbation. Simultaneously, the updated velocity field is used to update the effective-viscosity through Equations (4) and (5). Iterations continue until the solution has converged.

The converged velocity field is used to predict the shape of the free surface (and thus the domain geometry) at a future time step. Changes in domain geometry and the redistribution of mass within the (x, z) plane are accounted for when regridding the finite-volume mesh at the start of each time step.

Changes in flowband width

Our model accounts for the kinematic and dynamic effects of a changing flowband width, W = W(x, t ). During any single time step, spatial changes in flowband width affect the velocity field directly through continuity, and indirectly through the dependence of the effective-viscosity on strain rate. The direct effect is accounted for automatically when specifying W(x, t ) through Equation (8). The indirect effect is

accounted for by including a non-zero, transverse-normal strain rate in Equation (6), parameterized as

(10)

which can be derived from considerations of continuity along a flowband (Reference WaddingtonWaddington, 1982).

To maintain global continuity, special considerations are required when changing the flowband width from one time step to another. While regridding of the domain accounts for mass redistribution in the (x, z) plane (along the glacier centerline), flowband narrowing (or widening) over time leads to mass convergence towards (or divergence from) this plane. When the flowband width changes from one time step to another, this additional mass must be introduced (or removed) from the model domain. The solution method discussed above offers a natural means for doing this: in any volume affected by a time change in flowband width, the mass source term in Equation (9) is augmented by the additional amount

(11)

where ∂W/∂t is the rate of flowband narrowing (∂W/∂t < 0) or widening (∂W/∂t > 0) and ΔxΔz is the area of the volume face through which that mass passes. The negative sign in Equation (11) ensures that, at the glacier centerline, convergence due to flowband narrowing is treated as a mass source and that divergence due to flowband widening is treated as a mass sink. When ∂W/∂t 6≠ 0, mass converging on or leaving the plane of the glacier centerline is automatically redistributed by the flow field so that global continuity is maintained. As with along-flow changes in width, temporal changes in width add an additional strain rate component to the right-hand side of Equation (6). This term is given by

(12)

Results

We explored the sensitivity of modeled surface elevation and velocity fields to choices for bulk density ρ and scalar enhancement factor E. Bulk density obviously depends upon debris content, and it seems likely that E does too, with E decreasing as rock content increases. Our goal is not to fit observations exactly: capturing all of the observed surface detail would require an impractical number of grid cells and/or excessive tuning of the model. Instead, we aim to constrain the most reasonable values for E and ρ based on the overall trends in the misfit between the model and the observations. We tested the model for 0:01 ≤ E ≤1 and 918 ≤ ρ ≤ 1700 (kgm^–3) (representing debris content ranging from 0 to ∽50% by volume). While we present only a subset of all model results, additional results support the main points discussed below.

Misfit between the modeled and observed surface elevations for a density of 1100 kgm^–3 and enhancement factors of 1, 0.2 and 0.05 is shown in Figure 4. Note that ρ = 1100 kgm^–3 corresponds to about 15% rock debris by volume, the value estimated for the actual glacier as defined by the 1986 DEMs. While the misfit for all models is similar for the first ∽50 days, clear trends emerge by day 107. Models with E = 1 lead to a mass deficit in upstream regions and a mass surplus in downstream regions, while models with E = 0.05 show the opposite trend. For E = 0.2, the maximum misfit generally falls within 10m throughout the model run. When comparing modeled and measured surface velocities at ELE4 for the same simulations (Fig. 5), the reason for trends in surface elevation misfit in Figure 4 becomes obvious. For E = 1, horizontal velocities in the region of squeezing are too large and ice moves from the bulge towards the terminus too quickly. Conversely, for E = 0.05, velocities in the region of squeezing are too small and too much ice remains in the upper portion of the glacier after the squeezing stops. The model with E = 0.2 provides a good fit to the velocity at ELE4, and thus we adopt the case ρ = 1100 kgm^–3, E = 0.2 as a benchmark for further discussion. The decline in horizontal velocity over time, which is shown by all models, is a reflection of the decline in squeezing rate.

Fig. 4. Misfit (elevation_model – elevation_observed) for elapsed model times corresponding to dates of DEMs. Three enhancement factors were considered for ρ = 1100 kgm^–3, our best estimate of bulk density: E = 0.2 (thick solid line), E = 0.05 (dashed line), and E = 1 (dotted line). Results are also shown for pure ice (ρ = 918 kgm^–3) with E = 0.2 (upper thin solid line) and ice with ρ = 1400 kgm^–3 and E = 0.2. The gray-shaded region encompasses an elevation misfit of 5m.

Fig. 5. Model fit to observed horizontal velocities at ELE4. Gray dots are observations from Reference Walder, Schilling, Vallance and LaHusenWalder and others (2007b). Solid, dashed and dotted lines are as in Figure 4. The lower dash-dash-dot line is for the end member of stiff, low-density ice (E = 0.05 and ρ = 918 kgm^–3) and the upper dash-dash-dot line is for the end member of soft, high-density ice (E = 1 and ρ = 1400 kgm^–3).

Consider next the effect of varying ρ within the model. With E = 0.2, surface elevation misfit (Fig. 4) for ρ = 918 kgm^–3 (pure ice) and ρ = 1400 kgm^–3 (30% rock debris by volume) is worse than for the estimated actual density of 1100 kgm^–3, but still generally falls within the range of 10 m. The same variation in density yields an acceptable fit between measured and modeled surface velocity (Fig. 5).

While the end-member models tested here – stiff, ‘heavy’ ice (E = 0.05 and ρ = 1400 kgm^–3; dashed line ‘c’ in Fig. 6) and soft, ‘light’ ice (E = 1 and ρ = 918 kgm^–3; dotted line ‘a’ in Fig. 6) – result in a fit to the observations that is worse than our favored model, they highlight an important tradeoff between density and enhancement factor. The observations can be fit by assuming either stiff, dense ice or soft, less dense ice because a similar velocity field (not shown) results in either case. If there were no constraint on density other than a lower bound (the density of pure ice), we would find that observations could be fit equally well by assuming combinations of enhancement factor and density ranging from 0.1–0.5 and 918–1700 kgm^–3, respectively. In fact, however, we do have an independent constraint on bulk density and can therefore use the sensitivity tests to narrow the range of uncertainty in the enhancement factor. Our best estimate for density, as noted above, is 1100 kgm^–3, corresponding to about 15% rock debris by volume; the upper bound is ρ = 1400 kgm^–3, corresponding to about 30% rock debris by volume – the average debris content of all the crater-floor fill accumulated since 1980. If the density were 1400 kgm^–3, the enhancement factor would have to be close to 0.1 to match surface elevation and velocity as closely as our benchmark case (ρ = 1100 kgm^–3, E = 0.2).

Fig. 6. Misfit between modeled and observed elevations on day 193. Line types represent model results for relatively stiff (E = 0.05, dashed) and relatively soft (E = 1, dotted) ice and for densities of (a) 918, (b) 1100 and (c) 1400 kgm^–3. The black solid line shows the misfit for the benchmark values of density and enhancement.

Observed and computed surface profiles on several dates are shown in Figure 7 for our benchmark case and for the case of pure glacier ice with no enhancement (E = 1). The pure ice case gives a surface elevation profile with a systematic error – too low in the upper reach and too high in the lower reach – that worsens with time.

Fig. 7. Model fit to observed surface elevations (open circles). Model results for our benchmark case (ρ = 1100 kgm^–3, our best estimate of actual density, and E = 0.2) are given by the solid lines. Model results for ‘normal’ ice values of ρ = 918 kgm^–3 and E = 1 are shown by dashed lines. The seemingly good fit at early times for the model based on ‘normal’ ice is best interpreted with caution as the actual bulk density used in this calculation is certainly incorrect. The dotted line represents the initial surface elevation profile on 3 January 2005.

To explain why the downstream portion of East Crater Glacier remains largely stagnant during the squeezing event while the upstream portion undergoes dramatic changes, Reference Walder, LaHusen, Vallance and SchillingWalder and others (2005) suggested that ice in the downstream portion of the glacier, which was not squeezed, acted as a dam against the ice upstream. We augment their qualitative discussion by using the model to understand quantitatively why a dam would form. In Figure 8 we show the horizontal velocity, horizontal strain rate, and horizontal deviatoric stress from the model on day 50 (21 February 2005), for our benchmark model values. The horizontal axis in Figure 8 is approximately centered at the transition from very rapid to no lateral squeezing. Moreover, there is an abrupt decrease in flowband width just downglacier of this location, where East Crater Glacier passes through a narrow gap between the 1980–86 lava dome and the east crater wall (Fig. 1). According to the model, this combination of factors results in highly compressive longitudinal-strain rates and corresponding highly compressive longitudinal stresses. This, in turn, leads to a large negative stress gradient across this width transition, and for several hundred meters downstream, which results in a force pushing back upstream, resisting the ice flow. Across this width transition, drag at the glacier bed, drag against the valley sidewalls, and longitudinal-stress gradients are of approximately equal importance in resisting the flow. In reality, the transition from a narrowing to a steady flowband width may be more gradual than shown here, in which case the overall effect on the strain and stress fields, while less dramatic than shown in Figure 8, would be the same.

Fig. 8. Horizontal velocity, strain rate and stress on day 50 (21 February 2005). (a) Horizontal velocity (contour interval = 0.1 md^–1); (b) longitudinal strain rate (contour interval = 10^–3 d^–1) and (c) longitudinal deviatoric stress (contour interval = 60 kPa). The downstream limit of squeezing is at x ∽ 650 m.

Discussion

If flow of the debris-laden East Crater Glacier at Mount St Helens is properly characterized by Glen’s flow law with n = 3, then the ice must be stiffer than debris-free ice by a factor of 5 to 10, with the stiffness increasing as the rock-debris content rises. This inference from modeling contrasts with the discussion by Reference Jacka, Donoghue, Li, Budd and AndersenJacka and others (2003), who reviewed laboratory data and concluded that there is no relation between deformation rate and debris (sand) content for debris fractions up to 15% by volume. (For sand/silt fractions in excess of 50%, there does seem to be a clear correlation between stiffness and debris content (Reference Mangold, Allemand, Duval, Geraud and ThomasMangold and others, 2002)). We speculate that the difference between our conclusion and that of Reference Jacka, Donoghue, Li, Budd and AndersenJacka and others (2003) reflects the radical difference between a laboratory experiment and the natural ‘ice-squeezing’ experiment performed by the Mount St Helens lava dome, and in particular the associated scale effects. In the laboratory, because the grain size of the ice is likely to be comparable to the grain size of the debris, the presence of thin water films at ice–debris interfaces probably facilitates regelation and other deformation mechanisms. Rock debris within East Crater Glacier, in contrast, spans the grain-size range from silt to boulders, with the grain size of the particulate material commonly much greater than that of the ice; deformation mechanisms that may be effective at laboratory scale may play little role at field scale.

A cautionary note regarding the above results must be included here. Matching a modeled evolution of the surface shape with data required that the modeled length of glacier being squeezed was about 100m longer than suggested by observations. One possible reason for this discrepancy is that the region of squeezing needs to be slightly longer in the model to correct for three-dimensional effects not captured in our two-dimensional model. We assume a simple flowband geometry in both the region being squeezed and in the region for which the glacier width does not change over time. In reality, both the glacier width and the squeezing rate vary along flow. Our need to adjust the length of the squeeze zone may reflect our over-simplified parameterization of spatial variability of squeezing rate. Another more speculative possibility is that matching the observed surface elevation profiles does in fact require mass convergence additional to that implied by surface observations of the motion of the dome–glacier contact. For example, magma may have intruded beneath the west margin of the glacier or beneath the glacier itself, slightly downstream from the new lava dome.

By simplifying a three-dimensional flow problem to two dimensions, we have likely missed some details governing the deformation of East Crater Glacier during growth of the new lava dome. Nevertheless, it is satisfying that the simple approach taken here provides a good fit to the observations with only a minor amount of model tuning. Without the use of a full-stress flow model, this would likely not be the case. Thickening of the upper glacier during growth of the new dome is a straightforward consequence of continuity but the redistribution of this thickened ice is strongly dependent on longitudinal stresses.

Conclusions

Flow modeling constrained by surface velocity data and independent estimates for glacier bulk density demonstrates that the debris-laden ice of the East Crater Glacier is significantly stiffer than debris-free glacier ice. Our favored model, in which the ice contains ∽15% rock debris by volume, requires a flow-enhancement factor of 0.2 (that is, ice stiffer than normal by a factor of 5) to fit the data. If the ice contains 30% rock debris by volume (the upper bound), an enhancement factor of 0.1 (that is, ice stiffer than normal by a factor of 10) is indicated. The model further makes evident the existence of strong longitudinal-stress gradients in the glacier while it was being squeezed by the growing lava dome.