Introduction

Unsteady ice-shelf behavior evident within current mass-balance estimates, relict crevasse patterns, and ice-rumple configurations on the West Antarctic ice shelves may result primarily from ice-stream discharge fluctuations (Reference JezekJezek, 1984;Reference MacAyeal, Bindschadler, Shabtaie, Stephenson and Bentley MacAyeal and others, 1987;Reference Shabtaie and Bentley Shabtaie and Bentley, 1987;Shabtaie and others, unpublished). This hypothesis is motivated by the wide discharge variation evident within the West Antarctic ice-stream population (Reference McIntyreMcIntyre, 1986; Reference Shabtaie and BentleyShabtaie and Bentley, 1987), and implies that ice shelves may be modified by factors independent of atmospheric and oceanic climate.

Ice Streams B and C, which discharge into the Ross Ice Shelf (Fig. 1), exemplify the extreme stages of temporal ice-stream evolution that produce ice-shelf thickness and flow anomalies. Ice Stream B currently discharges more than is accumulated in its catchment (Reference Shabtaie and BentleyShabtaie and Bentley, 1987). Ice Stream C is a fossil of past drainage (Reference Bindschadler, MacAyeal, Stephenson, Veen van der and OerlemansBindschadler and others, 1987) and is distinguishable from stagnant inland ice only by its buried surface crevasses (Reference RoseRose, 1979). Unsteady ice-shelf behavior linked to Ice Streams B and C includes the current build-up of Crary Ice Rise (Reference MacAyeal, Bindschadler, Shabtaie, Stephenson and BentleyMacAyeal and others, 1987) and deviations between current flow lines and englacial debris streaks observed on the Ross Ice Shelf (Reference JezekJezek, 1984).

Initially, to characterize ice-stream-driven variability of West Antarctic ice shelves, a companion paper (Reference MacAyeal, Veen van der and OerlemansMacAyeal and Barcilon, 1988) describes ice-thickness and velocity anomalies on unconfined, one-dimensional ice tongues. Here we include the effects of coastal confinement and large-amplitude ice-stream forcing on two-dimensional ice-shelf fluctuation (the companion paper restricted analysis to small-amplitude changes). To accomplish this extension, we adopt numerical methods to solve the governing equations (the companion paper was based on van derReference Weertman Veen’s (1983) analytic treatment), and restrict our attention to an ideal, rectangular ice-shelf domain (Fig. 2) forced by a single ice stream. This restriction simplifies boundary conditions, which are poorly understood, and eliminates geometric effects associated with complex West Antarctic geography. Ice-stream discharge forcing is varied among a number of simulations to identify aspects of ice-stream dynamics critical to ice-shelf change.

Fig. 1. The eastern Ross Ice Shelf (polar stereographic projection). Black regions denote open water or ice-shelf rifts. Shaded regions denote floating ice shelf. Ice streams (both active and fossil) are designated by names A - F and are delineated by marginal crevasse bands (hatched). Features of the Ross Ice Shelf associated with possible ice-stream forcing include the present build-up of Crary Ice Rise, disagreement between relict extensions of the marginal crevasse bands and ice-shelf stream lines and the geographic co-ordination of the various ice rises with ice-stream outlets. The broken line represents the 400 m thickness contour.

Fig. 2. The ideal, rectangular ice-shelf domain (above) and time-average thickness (contoured in m) and velocity magnitude (in m/a) for test 2 (below). In all figures, the axis of symmetry is exploited to present contoured fields more efficiently.

Methods

Our numerical simulations are conducted using the simple, rectangular ice-shelf domain shown in Figure 2. Horizontal dimensions of the rectangular channel (300 km χ 150 km) and ice-stream entry (50 km width) represent typical Antarctic conditions. The ideal ice shelf studied here is smaller than the Ross or Filchner—Ronne Ice Shelves (approximately one-tenth of the surface area); however, the ideal geometry is selected to represent a typical unit cell consisting of a single ice stream and its ice-shelf effluent. Boundary conditions and rheological parameters defined below represent most of the essential features of ice-shelf/ice-stream interaction. For simplicity, and for lack of guidelines evident from field observations, we do not allow the ice-shelf boundaries to change with time nor allow grounding on the sea bed.

Time evolution of ice-shelf thickness, flow, and temperature are governed by stress equilibrium, mass conservation, and heat-flux continuity principles subject to appropriate boundary conditions and material properties. For the purpose of this study, we simplify the equations that express these principles by: (i) restricting solutions to leading-order balance conditions imposed by ice-shelf geometry (small aspect ratio implies depth-independent horizontal velocity; Reference Morland, Veen van der and OerlemansMorland, 1987;Reference Muszynski and Birchfield Muszynski and Birchfield, 1987; Reference MacAyeal, Veen van der and OerlemansMacAyeal and Barcilon, 1988), (ii) disregarding thermo-mechanical coupling through temperature-dependent flow-law parameters (Reference Williams and FlutterWilliams and Hutter, 1983), (iii) disregarding firn densification, (iv) adopting a simplified flow law that does not account for ice-crystal fabric, density, or chemical composition (such as salinity of basal ice), (v) applying constant surface accumulation and zero basal freezing, and (vi) restricting the lateral boundaries of the ice shelf to be fixed in time. We accordingly treat ρ (all variables are defined in Table I) as constant (917kg/m³) and express the flow-law constant B(ρ,θ) as a z-dependent function to approximate the effect of vertical temperature gradients in typical West Antarctic ice shelves. Horizontal and time derivatives of B(z), which arise generally as results of horizontal gradients in snow accumulation and basal melting not considered here (see Reference MacAyeal and ThomasMacAyeal and Thomas (1986) for discussion), are assumed to be zero.

Under these simplifications, the stress-equilibrium conditions in the three spatial directions may be written (derivations are presented inReference Morland, Veen van der and Oerlemans Morland (1987); see alsoReference MacAyeal and BarcilonMacAyeal and Barcilon (1988)):

The horizontal flow U = u(x,y,t)n_x + v(x,y,t)n_y is independent of z. The depth-average effective viscosity v ^z is strain-rate-dependent to account for power-law ice rheology,

where B ^z is the depth-average of B ^z , n is the dimensionless flow-law exponent (taken here to be 3; see discussion by Reference Doake and WolffDoake and Wolff (1985)), and the second invariant of the strain-rate e is (under the simplification of depth-independent horizontal velocity),

The mass-balance equation is written:

where the snow-accumulation rate M is assumed to be constant (0.25 m/a ice equivalent) and the basal melting rate is assumed to be zero.

Lateral boundary conditions on the stress-equilibrium equations consist of kinematic constraints on the ice-stream entry and other coastal boundaries:

and dynamic constraints at the ice front

where Γ′ is the contour of the ice-stream entry, and Γ is the contour of the ice-front (Fig. 2). The kinematic constraint on Γ′ treats ice-stream flux as an externally controlled forcing. In reality, this flux depends on ice-shelf configuration and could be parmeterized in more advanced studies. The no-slip condition applied at other coasts constitutes the most extreme condition of coastal restraint. The least extreme condition is treated in the companion paper (Reference MacAyeal and BarcilonMacAyeal and Barcilon, 1988). The dynamic constraint at the ice front implies that sea-water pressure alone acts on the ice front.

Lateral boundary conditions on the mass-continuity equation are expressed in terms of mass flux Q = UH and are,

and

The flux condition at the ice front implies an exact balance between forward advection and iceberg calving. This balance is required by our assumption of fixed ice-shelf geometry, and is supported by observation (Reference Jacobs, MacAyeal and ArdaiJacobs and others, 1986).

Equations (1)–(13) are solved using the finite-element method, an implicit time-step scheme (Reference MacAyeal and ThomasMacAyeal and Thomas, 1982;Reference MacAyeal, Veen van der and Oerlemans MacAyeal, 1987), and the method of weighted residuals (Reference Lapidus and PinderLapidus and Pinder, 1982). The method of weighted residuals requires that Equations (1), (2), and (6) be converted to integral constraints that incorporate all of the boundary conditions. We display this conversion for Equation (1) as follows: multiplication by an arbitrary weighting function R(x,y) and integration over the area of the ice shelf gives

If U = (u,v) is a solution to Equation (14) for arbitrary weighting functions R(x,y), then it must be a solution of Equation (1) as well (this solution is called the weak solution because it can violate Equation (1) in ways that still allow the integral constraint expressed by Equation (14) to be satisfied). To incorporate the boundary conditions, Equation (14) is integrated by parts to yield

where dλ. is a vector-line element along the contour defining the entire ice-shelf boundary. The kinematic constraints given by Equations (7) and (8) imply that we are not free to vary the solution U along the ice-stream entry or other coasts in our attempts to minimize the residual dissatisfaction of Equation (15). We can thus choose R = 0 along such boundaries without loss of generality. On the ice front, where dynamic constraints apply, R remains an arbitrary non-zero function. The dynamic boundary conditions along Г given by Equations (9) and (10), and the restrictions R =0 along the rest of the boundary, allow the boundary integrals of Equation (15) to cancel:

This is the most compact form of the integral constraint to be solved by the finite-element method.

For a clearer physical interpretation, we express Equation (16) in an alternative, less compact form that corresponds to energy conservation. We undo one of the integrations by parts done previously:

and re-express Equation (16) as

If R is interpreted as a virtual velocity in the x-direction, Equation (18) represents the conservation of work associated with this virtual motion. The left-hand side represents viscous energy dissipation. The terms or. the right-hand side represent: work done against gravity by vertical mass redistribution, work done against sea-water pressure along the sloping ice-shelf base, and work done against sea-water pressure at the ice front.

Equations (2) and (6) may be converted to integral constraints in a manner similar to that used to derive Equation (18). The results are

where S(x,y) and P(x,y) are arbitrary weighting functions. As with R, S has been set to zero at boundaries where kinematic constraints are applied. P, in contrast, is an arbitrary non-zero function at all boundaries because ice thickness is free to vary in response to the prescribed mass flux.

To discretize Equations (18)–(20), the domain is partitioned into 5 km × 10 km rectangular elements, H, u, v, P, R, and S are approximated as piecewise linear functions, and v ^z is approximated as a piecewise constant function. The requirement that Equations (18)-(20) hold for arbitrary P,R, and S yields a system of non-linear algebraic equations. These equations are solved (seeReference MacAyeal and Thomas MacAyeal and Thomas, 1986) for values of H(t + Δt). u(t), and v(t) at locations corresponding to nodes of the finite-element mesh and, for values of v ^z , at locations corresponding to the centers of the finite elements.

The numerical procedure was tested through comparison with exact solutions and with solutions obtained using different numerical schemes, and by monitoring the mass and mechanical energy budgets. The spreading rate associated with free-slip boundary conditions at the ice-shelf sides, for example, was compared withReference Williams and FlutterWeertman’s (1957) exact solution. Computational short-cuts, such as the lag of horizontal velocity by Δt in the left-hand side of Equation (20), were tested through comparisons with more accurate (and computationally expensive) numerical schemes. In circumstances where direct comparison with analytic results could not be made, mass and mechanical energy budgets were monitored to ensure their balance. As a means of eliminating small-amplitude node-to-node oscillations in the simulated thickness, artificial diffusion was added to the mass-balance equation (with a diffusivity of 5.0 × 10⁻²m²/s). This diffusion did not significantly modify the mass balance, and is comparable to the truncation error of the other terms. Sensitivity to parameters not varied in the tests conducted here (such as the flow-law parameters) is examined elsewhere (paper in preparation by M.A. Lange and D.R. MacAyeal).

Experimental Design

Several numerical experiments were conducted to determine ice-shelf response to a variety of ice-stream forcing scenarios (the ice-stream forcing parameters are listed in Table II). In one class of experiments, ice-stream forcing was periodic to determine ice-shelf sensitivity to frequency and phase lag between grounding-line thickness and velocity fluctuations. In another class of experiments, ice-stream changes were impulsive to determine how rapidly ice-shelf conditions could equilibrate with sudden changes of inland-ice drainage patterns.

In simulations of periodic ice-stream discharge, ice-stream boundary conditions U = (u,v) and Q = (q_x,q_y) were specified to be sinusoidal functions of time with period T = 2π/w:

where

and are the time-average discharge velocity and flux, respectively, U ₀and Q ₀ are the amplitudes of the

velocity and flux variations, respectively, and Φ is the phase lag between flux and velocity. In all simulations of this type, and were 20 km³/year and 500 m/year (representing typical Antarctic conditions;Reference Shabtaie and Bentley Shabtaie and Bentley, 1987), respectively, and T and Φ were varied to determine sensitivity of the ice-shelf response. After starting from arbitrary initial conditions (usually in steady state with median ice-stream discharge), seven complete cycles of ice-stream variation were simulated to dissipate transients induced by the initial conditions. On the eighth and subsequent cycles, model data were sampled to measure ice-shelf response. The time-step size (Δt) was 1/40T.

In simulations of response to impulsive ice-stream behavior, ice-stream velocity and flux were suddenly accelerated from a state of rest or decelerated from a state of full activity (tests 8 and 9 of Table II). In each case, initial conditions were in steady state with initial ice-stream discharge and subsequent time evolution was simulated (using 10 year time steps) until conditions equilibrated with the new forcing.

Response Linearity

Non-linearity of the governing equations and flow law suggests that ice-shelf response may not exhibit the same frequency of oscillation as ice-stream forcing. To test this possibility, we compared power spectra of ice-thickness and velocity anomalies generated by tests 1–3 (Table II). These anomalies are defined by

where H ^t and U ^t are the long-term average thickness and velocity fields. Time series of h and ǀ u ǀ were sampled at each time step during the eighth cycle of a simulation eight cycles long, and at a variety of locations throughout the ice-shelf domain. Power spectra were generated from these time series for the frequency range 0 ≤ x ≤ 64π/T, using the method ofReference Press, Flannery, Teukolsky and Vetterling Press and others (1986, p. 428).

For all locations tested, approximately 95% of variance in h and ǀ u ǀ was confined within a narrow spectral peak at the forcing frequency. Side peaks at the first and second harmonics were found to be largest (having amplitude approximately 5% of the peak at the forcing frequency) in tests 6 and 7 and at locations near the ice-stream outlet. These harmonics produce minor temporal asymmetry in the thickness-anomaly cycles. Ice-shelf response is thus essentially linear in the frequency and amplitude ranges tested despite non-linearity of the governing equations. One can understand this result by comparing the time-average ice-shelf configuration shown in Figure 2 with the amplitudes of h and ǀ u ǀ shown in Figures 3, 4, and 5. Products of anomaly amplitudes are insignificant when compared with products of anomaly amplitude and time-average amplitude. In this circumstance, harmonics of the forcing frequency are not forced strongly by interactions between anomalies at the fundamental forcing frequency.

The apparent response linearity demonstrated empirically above affords two conveniences in conducting our experiments. First, anomalies may be described by their Fourier components at the forcing frequency w = 2π/T:

where A _h and A _u are the amplitudes of thickness and velocity anomalies, respectively, and Φ _h and Φ _u are the phases of thickness and velocity anomalies, respectively, measured with respect to the phase of ice-stream velocity (increasing Φ corresponds to increasing time lag). Secondly, the amplitudes of h and ǀ u ǀ (A _h and A _u ) are linearly related to the amplitudes of ice-stream forcing (U ₀ and Q ₀ ). If U ₀ and Q ₀ are doubled (as between tests 2 and 3), A _h and A _u will also double at all points throughout the domain. Ice-stream forcing amplitude thus does not need to be varied as an independent parameter. Results of stronger or weaker forcing may be extrapolated from tests 4−7.

Fig. 3. Amplitude (upper half of each diagram, m) and phase (lower half of each diagram. degrees) of ice-thickness anomalies for tests 2. 4. and 5 (middle, above, and below, respectively), Ice-stream forcing period is varied between 300 and 1200 years, but ice-stream phase lag (between thickness and velocity fluctuations) is zero in all three tests.

Fig. 4. Amplitude (upper half of each diagram, m/a) and phase (lower half of each diagram, degrees) of velocity anomalies for tests 2, 4, and 5 (middle, above, and below, respectively). Ice-stream forcing period is varied between 300 and 1200 years, but ice-stream phase lag (between thickness and velocity fluctuations) is zero in all three tests.

Fig. 5. Amplitude (upper half of each diagram, m/a) and phase (lower half of each diagram, degrees) of thickness anomalies for tests 6 and 7 (above and below, respectively). Phase lag between ice-stream velocity and volume flux (or equivalently, thickness) is varied between +45 and –45°, but forcing period is constant (600 years).

Response To Periodic Discharge

Seven experiments were conducted to determine response to periodic ice-stream discharge. Details of each experiment, described in Table II, were varied to determine sensitivity to ice-stream fluctuation frequency and phase lag between ice-stream velocity and flux (or, equivalently, between ice-stream velocity and thickness).

Ice-shelf sensitivity to ice-stream forcing period T is displayed by comparing tests 2, 4, and 5. Forcing period T was specified to be 300, 600, and 1200 years in the three tests, respectively. Phase lag between ice-stream volume flux and velocity was the same (Φ = 0) in each test. The three forcing periods sampled were chosen because they are multiples of the advective time-scale of the ice shelf (600 years, given by the ratio of ice-shelf length to ). Contour maps of the amplitude and phase of h and ǀ u ǀ are presented in Figures 3 and 4, respectively. Time-average thickness and velocity magnitudes for test 2 are displayed in Figure 2 for information.

Characteristics common to tests 2, 4, and 5 include: (i) a monotonie increase of Φ _h along the longitudinal axis of the ice shelf and towards the rear corners of the ice shelf, (ii) three local maxima of A _h (one situated 50–100 km down-stream of the outlet and two situated on either side of the outlet where rifts would occur in Nature), (iii) little variation of Ф _u throughout the domain, and (iv) two local maxima of A _u (one located at the outlet, the other at the ice front).

Differences among tests 1, 4, and 5 reflect ice-shelf sensitivity to the time-scale of ice-stream variation. As displayed by Figure 3, thickness anomalies are most severely trapped within the vicinity of the ice-stream outlet at the smallest forcing period. Phase lag increases with decreasing ice-stream forcing frequency. These tendencies suggest that ice-shelf dynamics act as a low-pass filter to ice-stream fluctuations. Thickness anomalies furthermore appear to be propagated advectively along the two-dimensional analog of the slow characteristic discussed byReference MacAyeal and Barcilon MacAyeal and Barcilon (1988). In contrast with h, ǀ u ǀ remains locked in phase with ice-stream fluctuations at all forcing frequencies. Velocity anomalies appear to propagate instantaneously along the two-dimensional analog of the fast characteristic (Reference MacAyeal, Veen van der and OerlemansMacAyeal and Barcilon, 1988).

Sensitivity to phase difference between ice-stream velocity and volume flux (or, equivalently, thickness) is displayed by tests 6 and 7 in which Φ was 45° and −45°, respectively. As shown in Figure 5 (see also Fig. 3), thickness-anomaly amplitude is dramatically increased by non-zero phase angle, Φ. Grounding-line thickness anomalies are amplified by the need to accommodate volume-flux extremes at times in the ice-stream forcing cycle when ice-stream velocity is not extreme.

Impulsive Ice-stream Forcing

Conditions observed on the Ross Ice Shelf, most notably the stagnant state of Ice Stream C, suggest that ice-stream discharge fluctuations may be episodic and sudden in Nature (as if governed by fast physical processes such as subglacial hydrology). To investigate ice-shelf adjustment to discontinuous changes in ice-stream discharge of large amplitude, we conducted two simulations (Table II) in which ice-stream discharge was suddenly introduced or eliminated. In both tests, the initial ice-shelf thickness and flow were in steady state with either active or stagnant ice-stream discharge. Subsequent evolution was monitored to describe how the thickness and flow equilibrate with new ice-stream activity.

Plots of instantaneous ice-thickness distribution and the fraction of adjustment F(x,y,t) are presented in Figures 6 and 7 at arbitrary time intervals (50, 250, and 500 years; we present data at 250 years because Reference Shabtaie and BentleyShabtaie and Bentley (1987) conjectured that Ice Stream C may have suddenly halted 250 years ago). The fraction of adjustment F(x,y,t), defined by

where H _i(x,y) and H _f (x,y) are the initial and final steady-state thickness profiles, is displayed to show the progress of ice-shelf equilibration. In both tests, ice thickness adjusts most quickly along the longitudinal axis of the ice shelf. Thickness in the rear corners adjacent to the ice-stream outlet initially overshoots the final distribution and subsequently displays damped oscillations.

Comparison of the two experiments (Figs 6 and 7) reveals asymmetry between the two adjustment processes. Ice thickness in the switch-on experiment approaches its new steady-state distribution more rapidly than in the switch-off experiment. Fifty years after ice-stream switch-on, for example, F exceeds 0.9 at the ice-stream outlet; whereas, 50 years after switch-off, it remains below 0.5. This asymmetry can be explained by the differences in the advective time-scale (time required for an ice parcel to traverse the longitudinal axis of the ice shelf) between the two tests. This time-scale is reduced in the switch-on experiment because of the high velocity imposed at the ice-stream outlet.

Fig. 6. Adjustment parameter F (upper half of each diagram, non-dimensional units) and instantaneous thickness (lower half of each diagram, m) at various times after sudden introduction of ice-stream discharge (test 8, switch-on). Contours of F(x,y,t) in the up-stream corner of the ice shelf are not drawn because of crowding and small-scale detail.

Fig. 7. Adjustment parameter F (upper half of each diagram. non-dimensional units) and instantaneous thickness (lower half of each diagram, m) at various times after sudden elimination of ice-stream discharge (test 9. switch-off). Contours of F(x,y,t) in the up-stream corner of the ice shelf are not drawn because of crowding and small-scale detail.

Conclusion

Ice-shelf thickness and flow anomalies induced by ice-stream transience have been explored through time-dependent numerical simulations in which ice-stream discharge is prescribed as an external forcing. Thickness anomalies generated by periodic ice-stream discharge fluctuations exhibit: (i) amplitude trapping near the ice-stream outlet, with a decay scale determined by the ice-stream forcing time-scale, and (ii) appreciable phase lag across the ice-shelf domain. The first effect suggests that ice-shelf dynamics impose an effective low-pass filter to periodic ice-stream behavior. Ice-stream fluctuations with short time-scales (measured in relation to the advection time-scale of the ice shelf, which is approximately 600 years) produce little measurable fluctuation down-stream. Conversely, the influence of long time-scale fluctuations is widespread. Thickness anomalies are sensitive to phase lag in the timing of ice-stream velocity and thickness fluctuations. If this lag is zero, for example, ice-shelf thickness anomalies exhibit localized maxima down-stream of the grounding line. These down-stream maxima may correspond to points where ice rumples or rises would be generated in response to stochastic ice-stream variability.

Velocity anomalies generated by periodic ice-stream discharge are not trapped within the vicinity of the ice-stream outlet, and propagate instantaneously to most parts of the ice shelf. Typically, velocity anomalies exhibit two maxima: one at the outlet, and the other at the ice front. The first is due to the direct influence of the ice stream. The second is due to increased longitudinal spreading rates associated with thickness anomalies generated on the ice shelf.

The companion study (Reference MacAyeal and BarcilonMacAyeal and Barcilon, 1988) demonstrated that two distinct characteristics serve to transmit the effects of ice-stream fluctuation: a result which is borne out by the present study. Coastal confinement and two-dimensionality do not significantly alter the distinction between slow- and fast-propagation pathways for information dispersal. Thickness anomalies propagate on the slow path because they are constrained to move with the horizontal velocity of individual ice columns. Velocity anomalies, conversely, propagate instantaneously (or with seismic speeds) throughout the ice-shelf environment because they are governed by an elliptic boundary-value operator.

Acknowledgements

This research (Contribution No. 73 of the Alfred-Wegener-Institut für Polar- und Meeresforschung) was supported by the U.S. National Science Foundation (DPP 8509451), the Alfred-Wegener-Institut für Polar- und Meeresforschung, and the NATO Scientific Affairs Division (27–854/85). We thank G. Platzman, L.W. Morland, and an anonymous referee for comments on the manuscript.