Modeling the dynamic response of a crater glacier to lava-dome emplacement: Mount St Helens, Washington, USA

Abstract

The debris-rich glacier that grew in the crater of Mount St Helens after the volcano’s cataclysmic 1980 eruption was split in two by a new lava dome in 2004. For nearly six months, the eastern part of the glacier was squeezed against the crater wall as the lava dome expanded. Glacier thickness nearly doubled locally and surface speed increased substantially. As squeezing slowed and then stopped, surface speed fell and ice was redistributed downglacier. This sequence of events, which amounts to a field-scale experiment on the deformation of debris-rich ice at high strain rates, was interpreted using a two-dimensional flowband model. The best match between modeled and observed glacier surface motion, both vertical and horizontal, requires ice that is about 5 times stiffer and 1.2 times denser than normal, temperate ice. Results also indicate that lateral squeezing, and by inference lava-dome growth adjacent to the glacier, likely slowed over a period of about 30 days rather than stopping abruptly. This finding is supported by geodetic data documenting dome growth.

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 (Vallance 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 Walder and others, 2007a, b). 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 (Kamb 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 (van der Veen, 1999). The observed velocity difference (Fig. 3), however, was only a factor of 3–4. Walder 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 Anderson 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 (Mills, 1992). A similar exercise by Schilling 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 Walder and others (2007a; 2007b). 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) (Patankar, 1980; Versteeg and Malalasekera, 1995) and is described fully by Price (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. van 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 Price (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 (Patankar, 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 (Waddington, 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)

Resistance from valley sidewalls

To approximate the effects of drag against valley sidewalls, the body force (the first term on the left-hand side of Equation (1)) is multiplied by a shape factor, F _s ≤ 1. We use shape factor definitions from Nye (1965) assuming an elliptic-shaped channel, a width equal to our specified flowband width and a specified ice thickness.

Observations (Walder and others, 2007a; b) indicate that as the squeezing event progressed, the surface of East Crater Glacier bulged markedly and that a significant fraction of the glacier thickness was not in contact with either the lava dome or the east crater wall. This effect has to be accounted for to avoid overestimating the effects of side drag. On the basis of cross sections in Walder and others (2007a), we estimate that side drag was effective over a vertical distance of about 50 m greater than the pre-eruption, centerline ice thickness.

Initial conditions

Profiles of the initial glacier surface and bed elevations are necessary to define the initial domain geometry. We use surface elevation data from photogrammetrically-derived Digital Elevation Models (DEMs), as described in Walder and others (2007b). Our analysis is limited to a period of time for which the glacier can be adequately described by a two-dimensional flow model. We choose an initial glacier surface elevation profile, along the approximate centerline shown in Figure 1, from 3 January 2005 (model day 0), by which time East Crater Glacier was clearly separated from its western counterpart, and compression along the upper portion of the glacier was oriented almost entirely across-flow. Because the DEM for 3 January 2005 covers only the upstream ∽700m of the glacier centerline, we must estimate a surface elevation profile for the downstream ∽600m of the glacier centerline. For the region along flow from ∽900–1300m we use the next available surface elevation profile, that from 19 April 2005 (model day 107). A comparison of this and other profiles shows that there is very little elevation change along the stagnant, lower portion of the glacier. For the region of the centerline from 700–900m along flow, we have manually interpolated between the 3 January and 19 April profiles guided by the surface shape in subsequent elevation profiles. Bed elevations along the same profile are held constant and are based on surface elevation of the crater floor in October/ November 1986, as discussed above. Figure 2 shows these longitudinal elevation profiles as well as the finite volume grid at the start of the model run.

One other necessary initial condition is an estimate for the effective viscosity of the ice, which we take from Equation (5) assuming temperate ice and a value of 0.001 d^–1 for ɛ_{e .}

Boundary conditions

Deformation of the upper part of East Crater Glacier was clearly dominated first by the compression applied by the expanding lava dome, and afterwards by relaxation of the greatly thickened ice. In contrast, the speed of GPS station ICY4 (Fig. 3), downglacier of the squeeze zone (Fig. 1b), was in fact similar to what one would estimate from a balance-velocity argument (Walder and others, 2007b), and the surface elevation of the lower part of the glacier changed very little during the squeeze event. These observations motivate some simplifying assumptions with respect to boundary conditions. First, we assume that over the 193 day period considered here, the glacier’s response to traditional mass balance forcing is negligible relative to its response to the squeezing event; in other words, the thickening rate is much larger than anticipated surface elevation changes due to accumulation or ablation. This is justifiable because the average glacier thickening rate, from a total ice accumulation of about 80×10⁶m³ averaged over an area of about 1 km² in 20 years, is about 4ma^–1; in comparison, East Crater Glacier thickened at an average rate of about 0.6md^–1 during the squeeze event. Second, we assume that the entire glacier would be relatively stagnant over this same period of time in the absence of lava-dome intrusion. We therefore simply treat the terminus of the glacier as a zero-flux (u = 0) boundary. The upstream end of the glacier, which was pinned against the crater wall, is also treated as a zero flux boundary.

Walder and others (2007b) interpreted East Crater Glacier’s striking lack of both a spring speed-up and diurnal velocity fluctuations as evidence for absence of basal sliding. This unusual situation probably reflects absence of a drainage system conveying water along the glacier bed and instead a situation in which meltwater percolates into the thick rock-avalanche deposits underlying the glacier. We therefore assign a zero-slip boundary condition at the glacier bed.

The glacier surface, which is specified as stress free, evolves over time in response to the squeezing event. Its evolution serves as our primary constraint for matching the model output to the observations.

Finally, we must specify the rate at which the contact between the lava dome and the glacier moved, as this effect constitutes a mass source in our numerical scheme. Relevant data from Walder and others (2007b) are as follows: (i) Measurements using sequential DEMs show that, on average, the contact over the upper ∽550m of the glacier migrated at ∽1–2md^–1 between 3 January and 19 April 2005 (model days 0 to 107). (ii) A GPS station on the expanding dome moved eastwards at an average rate of ∽0.8md^–1 on model days 109 to 111. (iii) Comparison of DEMs for 19 April 2005, 13 May 2005 (model day 131), and 15 June 2005 (model day 164) shows that the contact moved no more than ∽5m after model day 107. A consistent interpretation requires that eastward dome growth must have greatly slowed after day 111, and motivates the following simple parameterization of squeezing rate: (1) From its upstream limit to 200m downstream (marked ‘a’ in Fig. 3), the rate of flowband narrowing increases linearly downglacier. (2) From 200m to some distance D downstream (marked ‘b’ in Fig. 3), the rate of flowband narrowing is held constant at 2md^–1 from 3 January 2005 (model day 0) through 3 March 2005 (model day 60), and 1md^–1 from then until 22 April 2005 (model day 110). (3) Downstream of D, there is no squeezing and the flowband width is held constant. (4) The squeezing rate after 22 April (model day 110) decreases exponentially in time over a period of ∽60 days with a characteristic timescale of ∽7 days. The net displacement of the glacier–dome contact over this period of time is ∽7m, which is in good agreement with observations that suggest ≤5m of motion. Observations suggest that D = 550 m, but we find computationally that mass conservation requires D to be closer to 650 m. That is, with D = 550 m, not enough mass converges on the glacier centerline near x = 600m to build the observed surface bulge (this finding is independent of the sensitivity tests discussed below). We discuss the possible reasons for this discrepancy below.

Other observational constraints

We use DEM-derived surface elevation profiles (Fig. 2) as primary targets for flow modeling that steps forward in time. Surface elevation has an uncertainty of a few decimeters (Walder and others, 2007a) except where the glacier is highly fractured – meaning the uppermost 200–300m long reach of the glacier – and we do not attempt to match the surface elevations there. As a secondary check, we compare modeled and measured surface velocities on the glacier during the squeezing event. Surface velocities were derived from GPS-measured positions of stations deployed on the glacier surface. The helicopter-deployed GPS stations were designed for volcano monitoring during the eruption and were available for glacier monitoring only sporadically. Fortuitously, the record at one station, ELE4, encompassed the period of time during which squeezing of the glacier slowed and ultimately stopped. Here, we focus on fitting the observed velocities at this station because (1) it is centrally located within the region of squeezing; (2) it contains the longest record of glacier velocity and (3) it requires that we match not only the glacier velocity but also the change in glacier velocity over time.

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 Walder 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, Walder 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 Jacka 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 (Mangold and others, 2002)). We speculate that the difference between our conclusion and that of Jacka 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.