Notation

1 Introduction

The stability of the West Antarctic ice sheet (WAIS) has been in question since the early 1970s when Reference WeertmanWeertman (1974) made the argument that ice sheets grounded below sea level are unstable and could collapse catastrophically in response to a modest increase in sea level. Characterizing the present and future stability of the WAIS is complicated by the presence of several ice streams. Their capacity to transport mass rapidly from the interior to the sea possibly increases the potential for rapid collapse of the WAIS. Understanding the past behavior of these ice streams should provide important clues to their role in WAIS stability.

Inter-ice-stream ridges are good places to look for evidence of past ice-stream activity, because a record of changes in ice flow at the edges of these ridges may be recorded in their geometry and internal stratigraphy. Siple Dome is the ridge between Ice Streams C and D in West Antarctica (Fig. 1). There is evidence for recent change at its boundaries. Reference Retzlaff and BentleyRetzlaff and Bentley (1993) used ice-penetrating radar to detect buried crevasses of the former shear margin and date the shut-down of Ice Stream C at about 130years ago. Jacobel and others (1996a) showed that internal layers are truncated beneath a linear topographic feature on the northeast flank of Siple Dome and suggested that this feature is a formerly active margin. It is possible that these or other past events have altered the ice flow at Siple Dome and have left a signal in the pattern of internal layers.

In this paper, we present evidence for migration of the Siple Dome divide over the past several thousand years at a rate of 0.05-0.50 m a⁻¹ northward toward Ice Stream D. Our results are based on analysis of the shape of internal layers obtained from radar measurements in the vicinity of the divide. Ice-flow models and inverse techniques are used to determine the rate of divide motion and its sensitivity to the pattern of accumulation near the divide. We do not attempt to identify the specific cause of divide migration in this analysis. Reference Nereson, Hindmarsh and RaymondNereson and others (1998) have found that this divide migration could be caused by either an increase (decrease) in elevation of the Ice Stream D (Ice Stream C) -side boundary of Siple Dome by less than 100 in or by a modest change in the accumulation pattern over the past several thousand years.

Fig. 1.

AVHRR image of Siple Dome and he Stream C and D. Radio echo-sounding(RES) profiles measured in 1994 and 1996 are shown by the white lines. Courtesy of T. Scambos,Nalional Snow and Ice Data Center, University of Colorado.

2. Measurements

The geometries of the surface, bed and internal layers in the vicinity of the Siple Dome summit were obtained from Global Positioning System (GPS) surveys and broad-hand, monopulse radio echo-sounding (RES) measurements made in 1994 and 1996 as part of a collaborative project among the University of Washington, St. Olaf College and the University of Colorado. Detailed RES measurements were made every 20-100 m along a 4 km profile across the summit using several frequency bands (dominant wavelength 25-80m in ice, Table 1). Low frequencies were chosen to reveal deep internal layers while higher frequencies were chosen to reveal shallower layers. Each RES measurement is a voltage time series that is filtered with a zero-phase, fourth-order Butterworth filter at 1-10 MHz to reduce high-frequency noise and low-frequency coupling within the radar system. The signal voltage is mapped to color. Sequential, equally spaced measurements of a profile are plotted together to reveal the shape of the bedrock and internal layer pattern. The internal layers arc caused by variations of electrical properties in the ice (Reference HarrisonHarrison, 1973; Moore and others, 1992). They are assumed to be isochrones, that is, former ice-sheet surfaces which have been buried and deformed over time by ice flow (Reference HammerHammer, 1980).

The measurements show that Siple Dome is a two-di-mensional ridge overlying flat bedrock with ice thickness 1009 ± 7 m at the summit (Reference Raymond, Nereson, Gades, Jacobel and ScambosRaymond and others, 1995; Reference Raymond, Nereson, Gades, Jacobel and ScambosScambos and Nereson, 1995; Reference Jacobel, Scambos, Raymond and GadesJacobel and others, 1996b). figure 2a shows a 10 km cross-section of Siple Dome obtained from RES measurements with a wavelength of~80 m in ice. A detailed 4 km profile measured using a shorter wavelength (≈35 m) is shown in figure 2b. The internal layers are digitized for shape analysis by selecting the signal travel time associated with the maximum or minimum reflection amplitude in a hand-prescribed time window. Each digitized internal layer is then smoothed horizontally with either a 100 m box-car filter (for profiles with 20 m measurement spacing) or a zero-phase, low-pass filter routine with cut-off frequency at 5x10 ⁻⁴m⁻¹ (for long profiles with 100 m measurement spacing). Travel time is converted to depth using a ray-tracing program that includes a correction for the variation of density with depth (Reference WeertmanWeertman, 1993). Firn density is based on measurements of a 160 m core taken near the divide in 1994 (Reference Mayewski, Twickler and WhitlowMayewski and others, 1995).

Table. 1.

Summary of RES measurements

Figure 2 shows that the internal layers have a shape that varies with depth, with no relation to surface or bed topography. The layers are warped convex-up in a ≈2 km wide zone beneath the divide with a maximum upward displacement of about 50 m (Reference NeresonNereson and Raymond, 1996). The pattern is asymmetric. The depth to a given layer on the north flank of Siple Dome is greater than its depth to the south. Figure 3 shows that variation in shape is small along the divide ridge, so that the layer shapes are largely two-dimensional in a vertical plane perpendicular to the ridge (Reference Jacobel, Scambos, Raymond and GadesJacobel and others, 1996b). The positions of the apex of the curved layers are not vertically aligned; the axis is tilted toward the north side of Siple Dome by about 60????????? from vertical (fig. 2b; Reference NeresonNereson and Raymond, 1996).

Fig. 2.

(a) A 10 km radar profile of Siple Dome measured perpendicular to the divide at the summit. Measurements were made every 100m with a center frequency of 2 MHz. (b) Shallow part of a detailed radar profile where measurements were made every 20m with a center frequency of 5MHz. The data in both panels have been corrected for GPS-determined topography and hand pass filtered at 1-10 MHZ. The black band over the top 80-100 m of the profiles denotes the part of the returned RES signal obscured by the direct airwave. The apex of selected internal layers are marked with black dots in the lower panel. The solid line shows the trend of the apex axis.

Regardless of the cause of the local warping of the internal layers, the warping is likely a feature associated with the presence of the divide. Therefore, we interpret the tilt of the apex axis as direct evidence for divide migration. However, inferring the rate of migration is not direct. An ice-flow model predicts that a layer at 60? depth is about 10⁴ years old (Reference NeresonNereson and others, 1996). The horizontal displacement of its apex from the present divide position is about 700 m, which might suggest a migration rate of 0.07 m a⁻¹. However, this simple analysis fails to include the integrated effect of ice flow on the shapes of isochrones over time.

Fig. 3.

Three-dimensional representation of a 10x10x1km cube centered at the Siple Dome summit showing surface-elevation contours, bed topography and the shape of one prominent internal layer identified from RES measurements. Horizontal coordinates are in kilometers and denote distance “north” and “east”of the summit relative to the east west-trending divide ridge. True geographic north is 16.4° west of local “north”. Vertical coordinates are in meters above the WGS84 ellipsoid.

3 Flow Models

Our goal is to determine the rate of divide migration and its sensitivity to unknown parameters. This goal requires a forward model that predicts the shape of an isochrone layer resulting from its deformation since deposition arid includes a near-divide deformation process which produces up-warped isochrone layers.

Reference RaymondRaymond (1983) showed that the non-linear creep of ice leads to an anomalous flow pattern spanning a few ice thicknesses beneath an ice divide. Deep ice is relatively ”stiff” due to low deviatoric stress, so that downward-moving isochrones are predicted to drape over this stiff zone and become warped convex-up. A two-dimensional, steady-state, finite-element model (FEM) of this description of ice flow is used as a reference model against which the observed layer shapes are initially compared. We do not attempt to account for effects such as the linear behavior of ice under low shear stress (Reference MellorMellor and Testa, 1969;) Reference Doake and WolffDoake and WolfF, 1985; Reference AlleyAlley, 1992; Reference Nereson, Waddington and RaymondWaddington and others, 1996;) or the development of crystal fabrics on the deformation of isochrones. Rather, we consider the special divide flow regime predicted by an isotropic. non-linear flow law as one possible cause of warped isochrones.

Up-warped isochrones can also be caused by slow burial and low vertical velocities associated with a local minimum in the spatial accumulation pattern over the divide (as found by Reference Fisher, Koerner, Patersnn, Dansgaard and ReehFisher and others (1983) at Agassiz Ice Cap). This is a second possible process leading to up-warped isochrones.

The actual deformation of isochrones near ice divides may be caused by some combination of anomalous flow from non-linear ice deformation and from the spatial accumulation pattern. We estimate divide migration rate at Siple Dome separately for each deformation process. The range of the results is assumed to be indicative of the uncertainty which arises from our lack of knowledge about the exact nature of deformation near the Siple Dome ice divide.

The FEM provides the basis for constructing a kinematic flow model which is used to model the evolution of isochrone shapes for each case. The flow field in the kinematic model is shifted to simulate divide migration and adjusted to account for spatial variations in accumulation rate. The evolution of isochrone shapes is found by tracking particles in the resulting time-dependent flow field. The divide migration rate is then determined by the shift of the modeled flow field required to match the shape of internal layers observed at Siple Dome.

3.2. Kinematic representation of flow field

The FEM is computationally intensive and does not allow for time-dependent effects such as migration of the divide flow field. A kinematic model is constructed to match the spatial variation of the flow field from the finite-element model calculation. This kinematic flow model is computationally efficient and allows simulation of divide migration. It is time-dependent in the sense that it calculates the effect of a migrating divide on isochrone deformation. The model does not allow for varying rates of divide migration. The migration rate and accumulation pattern are assumed constant for all time. Therefore, the initial divide position is outside the model domain.

Fig. 4.

Comparison between isochrone shapes predicted by a steady-state ice-flow model with a nan-linear flow law and the shapes of internal layers observed using RES. (a) Modeled layers in gray and observed layers in black, (b) Model and data discrepancy plotted according to the gray scale to show the spatial pattern of the mismatch.

We assume that the ice thickness H is constant in the ± 4 km region around the divide, since radar measurements show that the ice thickness there varies by less than 2?. We also assume that H has not changed over time. The depth-averaged horizontal velocity u(x), the mass balance b(x) and the distance along flow from the divide x are related by

(1)

The depth variation of horizontal velocity can be written in terms of a shape function that varies with position x and height above the bed z in the vicinity of the summit

(2)

(3)

Combining Equations (2) and (3) and integrating over depth gives

(4)

Equations (1) and (2) then require

(5)

where Ξ(χ) = ]^{χ. z)(lz β(χ) = jj /;(.r)d;r.Thc v lirai velocity field w(:r, s) is found from incompressibilky

(6)

which gives

(7)

Equations (5) and (7) describe a continuous two-dimensional flow field which includes a special near-divide flow regime. Tracking the motion of ice particles in this flow field through time predicts the shape of isochrone horizons. In the simplest case, £(x, z) = 1, b(x) = b₀. and Equations (5) and (7) reduce to familiar expressions used to derive the “Nye” time-scale (e.g. Reference PatersonPaterson, 1994, p. 277): u(x,z) = (b₀,/H) x and w(x,z) = (-b₀H)z.

To simulate divide migration, the flow field is moved in the same reference frame as the divide, so that the horizontal velocity in the divide reference frame u* (x*, z*) is

(8)

where M = m/b(0) is a constant divide-migration rate m scaled to a reference accumulation rate b(0) which is taken to be 0.10-0.15ma⁻¹ ice equivalent representing a range about the modern rate of 0.13 m a⁻¹ (personal communication from K. Kreutz), x* = x — x,x_div and z* = z. With this construction, x* = 0 is always the divide position. The expression for vertical velocity w{x,z) is unchanged so that w*(x*, z*) = w(x*, z*). We drop the asterisks from now on, assuming we are always in the moving divide reference frame and x = O is the divide position.

3.3. Case 1: non-linear divide deformation

The first model parameterization assumes that internal layer warping beneath the divide is caused only by the anomalous flow regime predicted near ice divides from a non-linear ice-flow law. The shape functions Q and ζ,??ν in Equation (3) are described by an eighth-order Chebyshev Polynomial fit to the FEM How field at x=10km and. X=0 km, respectively. The partitioning function φ(x) is chosen so that the horizontal velocity u(x, z) given by Equation (5) matches u(x, z) predicted by the FEM. We find that φ(x) is well approximated by the continuously differentiable, even function:

(9)

where l₁ =775m and l₂ =3740m represent scale lengths which define the width of the divide zone. This partitioning function illustrates that the characteristic width of this zone associated with divide flow is about two ice thicknesses (~2 km) (Reference RaymondRaymond, 1983; Reference HvidbergHvidberg, 1996). Figure 5 shows the horizontal velocity shape functions from the FEM and front Equations (3) and (9). Despite the simplifications involved in defining {x,z), the kinematic model reproduces the layer shapes predicted by the finite-element reference model with an rms error of 5m.

Fig. 5.

Horizontal velocity shape functions φ(x,z)from the finite-element model calculation (dashed line) and from Equations (3) and (9) (solid line) for x – 0,0.5,1.0,3.0, 5.0 and 10.0 km.

A simple spatial variation in accumulation b(x) is allowed, because the asymmetry of the observed layers implies some accumulation gradient across the divide independent of the layer warping. We describe a hinge-like accumulation pattern with two gradient parameters G _s and G_n which describe the accumulation gradient on the south (x < 0) and north (x > 0) sides of the divide, respectively:

(10)

where b(0) is a reference accumulation rate and x — x/H. Since we assume that the ice thickness H is constant, we scale the velocity pattern to the accumulation rate so that Equation (1) holds. velocity Equations (5) and (7), with shape functions given by Equations (3) and (9) and parameters given by Equations (8) and (10), define the kinematic model for isochrone deformation due to a non-linear flow-law. There are three non-dimensional free parameters, M, G_s, and G_n which we vary and find the combination which best fits the observed layer shapes.

3.4. Case 2: accumulation minimum at divide

The second parameterization assumes that internal layer warping is caused solely by a local accumulation minimum over the divide. We construct a uniform flow field with no inherent special divide flow zone by setting φ(x)= 0 in Equation (3) thereby removing the spatial variation in the velocity shape function [x, z) so that

(11)

We further define a new accumulation pattern that allows a local minimum over the divide. This pattern could be caused by wind scouring at the divide and redeposition elsewhere (Reference Fisher, Koerner, Patersnn, Dansgaard and ReehFisher and others, 1983) or by a deposition pattern which causes less snow to be deposited over the divide relative to die Hanks. The snowfall at Siple Dome may be associated with meso-scale cyclones which track from the north or northeast side of Siple Dome (Reference BromwichBromwich, 1988; Reference Carrasco, Bromwich and LiuCarrasco and others, 1997). Numerical models suggest that winter surface winds are likely from the north to northeast (Reference Bromwich and ParishBromwich and others, 1994). Wind perpendicular to the cast-west-trending divide could cause scouring and redeposition on the lee side. Orographie effects may cause more deposition on the north flank of Siple Dome and a precipitation "shadow" on the south flank. We simulate these effects by superimposing one period of a cosine curve on the best-fitting hinge-like accumulation pattern from case 1:

(12)

where G*_s, and G*_n emerge from case 1 since the far-field layer shapes suggest an accumulation gradient independent of the divide zone.

Using this second parameterization, the model flow field u(x, z) and w(x, z) is defined by Equations (5), (7) and (8) with the shape function given by Equation (11) and accumulation pattern given by Equation (12). We vary three non-dimensional parameters to find the best fit to the observations: the scaled amplitude A and wavelength A of the superimposed cosine accumulation curve and the scaled divide migration rate M.

4 Choosing The Best Model

Because we do not know a priori the accumulation history of Siple Dome or the age of any layer, our goal is to match layer shapes rather than any particular isochrone. We must therefore decide which modeled layer we compare to a given observed layer in a way that makes a unique one-to-one assignment of observed to modeled layers. Options include comparing layers with the same height (1) at the present divide x = 0, (2) at their shallowest point, (3) at either edge of the domain or (4) when averaged over the domain. We want to emphasize the shape of layers near the divide to determine the divide migration rate. Options (1) and (2) are therefore poor choices, because they emphasize the mismatch between the model and the data at edges of the domain where we expect the effects of a migrating divide to be smallest. Option (3) presents problems when the layers are asymmetric. We choose option (4) and compare modeled and observed layers with the same average height z. The shape S of each layer is defined as its height variation relative to its average height:

(13)

For a discrete set of digitized layers, each layer is denoted by the index j and identified by ils average depth Z_j. Points along the horizontal coordinate are denoted by i The shape of each layer is quantified by measuring the elevation S _ij at each position i along layer j so that

(14)

To compare modeled (S^m _ij) and observed (S ^d _ij) layer shapes (data), we define a mismatch index J;

(15)

The term e (S^m _ij) is the expected combined error from the model and the data for a particular layer j. A weighting function fciy for point î on layer j is chosen to give increased weight to the residual in the divide area. The number Nj is the total number of points i along a layer j and L — 30 is the total number of layers, T = J] ? Nj is the total number of points sampled in the domain and p is one less than the number of free parameters in the model. In this application, each layer is sampled every 200 m, and Nj depends on the length of the RES profile (3 7 km).

4.1. The expected error

The combined error term Cj represents errors in layer shape S{Xi, Zj, z) and arises from errors in determining the elevation z of a point at, r, on a given modeled or observed layer and in determining the average height of the laver ζ,. For we define

(16)

where all values on the right hand side correspond to layer j and subscripts m and d denote modeled and observed values, respectively.

Both the model and the data contribute to the total error. We do not include errors due to using possibly wrong assumptions in the flow model (flow law, temperature distribution, etc.), because instead we consider two distinct descriptions of internal layer deformation. However, there is an error in z_m associated with simplifications used to define the kinematic model. Based on a comparison to the FRM, these model simplifications correspond to △z_m =1 -5 m, depending on the height of the layer. We assume no error in average depth so that △z_m = 0.

Main sources of errors in the data include (1) picking errors: accurately defining an internal layer from the radar data, (2) processing errors: filtering and smoothing the internal layers, and (3) physics errors: converting the signal travel lime to depth. Part of the error △z_d is due to random noise associated with picking a layer from the raw radar data which contributes to (1). This error is significantly reduced during the smoothing process and amounts to less than 1 m. Another contribution to Δz_d is uncertainty associated with horizontal space between radar measurements and how well this spacing is known: △x_d We somewhat arbitrarily set △x_d to be half of the measurement spacing: 10-50 m. The contribution to Az,\ is small (<2m) because the layer shapes vary only slightly with position. The error in average depth of a layer Az,\ is systematic for a given layer and arises from identifying a specific return time associated with an internal layer from a rather broad reflection signal and from converting the signal travel lime to depth. Both of these depend on signal frequency and we take the error to amount to about one-tenth of the wavelength of the re-flected pulse: about 2- 10 m. Based on these estimates, the total combined error e_j, From Equation (16) ranges from 3 to 6 m, depending on the average layer height and signal frequency.

4.2. The weighting function

The simplest weighting scheme is to prescribe equal weight to all of the residuals by setting w_ij = 1. However, we are primarily interested in determining the rate of divide migration and information about divide position is contained in the shape of the "divide bump". We therefore want to give greater weight to the residuals associated wit h this feature so that residuals occurring where the divide signature is most pronounced are given more weight than those where the signature is weak.

We define a weighting function ω_ij which is proportional to the component of the modeled layer shape that is caused by deformation processes associated with the presence of the divide. This divide shape component β_ij is obtained by subtracting the layer shape predicted under no special divide flow or divide accumulation pattern from the modeled layer shape S^m _ij.This removes the shape associated with bed topography and large-scale accumulation patterns. The resulting shape contains only the signature from the divide flow field. For Siple Dome layers, this adjustment amounts to removing any linear trend present in the modeled layer shape. The divide signature shape β_ij is generally positive because the divide deformation processes considered here produce up-warped isochrones. We do not consider models where the divide deformation produces down-warped isochrones and negative β_ij values. The divide signature shape is normalized so that

(17)

and

(18)

With this definition of the weighting function, residuals far from the predicted divide signature receive little weighting, and layers with a more pronounced divide signature are weighted more heavily than layers with no signature such as those found near the surface and bed. Because the weighting function is normalized, models which produce a mismatch index J < 1 match the observations to within the weighted errors (a 68% confidence level) and models corresponding to, J < 2 fit l be data to within twice the weighted errors (a 95% confidence level). Over- or underestimation of the errors would change the value of the minimum, J (e.g. from 1 to 2) but not its dependence on the model parameters. Since we have only three free parameters for each case, we can explore the entire parameter space and find the shape of the three-dimensional parameter volume that defines the best fit (minimum J).)

Our results arc not sensitive to a reasonable choice of weighting functions. Prescribing equal weight to all points along a given layer by choosing w_ij == 1 does not change the value of the three best-fitting parameters for either case; though the overall value of the mismatch J increases slightly for case 1.

5 Results

Figure 6 shows a comparison between the shape of the observed layers and the hivers predicted by the kinematic model with the combination of parameters that produces the smallest J for case 1 where the divide bump is caused by a non-linear ice-flow law. The small area of larger residuals likely indicates that the simple hinge-like spatial pal-tern of accumulation allowed in the case J model is not sufficiently flexible to allow an exact match between model and observations. The real accumulation pattern is probably more complex. Figure 7 shows how the mismatch parameter J depends on the three model parameters G_s, G_n. The models with parameters within the J < 1 contour fit the data to within the expected errors. Given this threshold, the inferred scaled migration rate M = 2.0-3.5 corresponds to a divide-migration rate M = 0.2-0.5 ma⁻¹ northward toward Ice Stream D when we take b(0)= 0.10-0,15 m a⁻¹. This is three to seven times faster than the migration rate inferred from the angle of the apex axis and a hypothetical depth-age scale. The inferred migration rate is relatively insensitive to the accumulation gradient south of the divide G _s (Fig. 7). Northward divide migration (toward Ice Stream D) smears out layer shapes to the south, leaving a “wake” of relatively flat layers and obscuring information about accumulation patterns there (paper in preparation by E. Waddington and others). Therefore, G_s is not well con-strained and ranges from about -0.03 to 0.02. The minimization indicates a relatively strong positive accumulation gradient north of the divide with G _n =0.05-0.07.

Fig. 6.

Best-fit case assuming the internal layer warping is caused primarily by non-linear ice flow. M=3, G_n = 0.06, G_s = -0.01. (a) Modeled layers in gray and observed layers in black, lb) Model and data discrepancy plotted according to the gray scale la show the spatial pattern of the mismatch.

Fig. 7.

Mismatch index J as a function of each model parameter for case 1: non-linear ice-flow law. Each panel is a cross-section of J taken through the minimum at model parameter values G_s = -0.01, Gs = 0.06 and M = 3, respectively. Contour lines for J are shown in intervals of 0.5.

For case 2, where the divide bump is caused by a local accumulation low, we choose G* _n = 0.06 and G*_s = 0 in Equation (12) as derived in case 1. We also choose p = 4 in Equation (15) to reflect two additional degrees of freedom, because we allow a more complicated accumulation pattern and Lise results from case 1 in Equation (12). Figure 8 shows a comparison between the observed layers and the predicted layers for the combination of A, λ and M which give the minimum J. Figure 9 shows the shape of J. Overall, the accumulation-low model produces a slightly better fit (J < 1) to the observations than the non-linear flow-law model, probably because of the curvature allowed in the prescribed accumulation pattern. However, this slightly improved fit is beneath the level of expected errors and does not imply that case 2 is more plausible than case 1. The predicted migration rate is slower than ease 1, with M ranging from 0.5 to 2.0 at the J 1 level· The predicted migration rate Mis largely independent of the amplitude A or the wavelength λ of the accumulation pattern. The predicted amplitude A of the accumulation low ranges from 2 to 6% of b(0).The predicted wavelength λ of the accumulation feature varies between 2.5 and 6.5 ice thicknesses (2.5-6.5km), increasing with amplitude A. This range of λ is expected, since the warped feature observed in the layer shapes spans about two to four ice thicknesses, suggesting λ ≈4-8. The accumulation patterns predicted by this model are reasonable possibilities near the divide.

Fig. 8.

Best-fit case assuming the internal layer warping is caused by a local low in accumulation over the divide. M = 1, A = 0.04, λ = 3. (a) Modeled layers in gray and observed layers in black. (b) Model and data discrepancy plotted according to the gray scale to show the spatial pattern of the mismatch.

The lower predicted migration rates for case 2 (accumulation low) arise from the fact that a combination of A and λ alone can be found to produce a sufficient match to the general shape of the observed layers with no divide migration. Only a slight shift of the flow field (small migration rate) is required to bring the modeled layers into agreement with observations. For case 1 (non-linear flow), the deformation field which produces the divide bumps is largely fixed by our assumptions about ice dynamics, and a more significant migration of the flow field is required to produce agreement between the model and the data.

6 Discussion

It is possible that the real situation at the divide is best described by a combination of a special flow regime associated with non-linear flow and a complex accumulation pattern. Assuming that case 1 and case 2 are equally likely and, using the J 1 mismatch level as a reasonable threshold, we estimate that the Siple Dome divide is migrating northward toward Ire Stream D at about 0.05-0.50 m a⁻¹, corresponding to M = 05-3.5 and b(0) =0.10 -0.15 m a⁻¹. If we have over-(under)-estimated the errors, then the predicted range would be smaller (larger). The simple estimation of the migration rate from the angle of the layer apex axis and a presumed depth-age scale produces a value in the low end of this range at about 0.07 m a⁻¹. The predicted rate is only slightly sensitive to the presumed mechanism which produces the internal layer warping beneath the divide. The non-linear divide-flow mechanism predicts higher migration rates than the accumulation minimum mechanism. However, this analysis does not distinguish which mechanism is operating at Siple Dome or reveal any information about the ice-flow law. Resolution of these questions will require information from ice cores and measurements of the depth variation of strain.

Fig. 9.

Mismatch index J as a function of each model parameter for case 2: local accumulation low. Each panel shows a cross-section of J taken through the minimum at parameter values λ= 3, A = 0.04 and M=1 respectively. Contain lines for J are shown at intervals of 0.5.

A large south north accumulation gradient of 5 10x10⁻⁶, a⁻¹ is a robust feature of our analysis in the 2-4km zone north of the divide. This result is independent of the migration rate and is consistent with both mechanisms of bump formation considered here. The predicted accumulation gradient would be artificially high if the divide zone calculated using the FEM is too wide (giving low vertical velocities over a large zone). However, the accumulation pattern predicted for case 2 (no divide zone from non-linear flow law) does not require a smaller gradient. Therefore, the large predicted accumulation gradient near the divide is a requirement of the data, not an artifact of the calculation. Measurements from shallow cores and snow pits in the vicinity of the Siple Dome summit indicate a regional-scale south-north accumulation gradient over the divide which is 10 30% of our prediction (personal communications from K. Kreutz and K. Taylor). We expect that the strong-accumulation gradient predicted here is a divide-local effect and must decrease with distance from the divide. The available field sampling is too sparse to detect the localized pattern suggested by this analysis.

We do not allow for time variation of the migration rate. If past changes in migration rate or accumulation pattern have occurred, then the discrepancy between the model and the data would be depth-dependent. A model which fits the shallower (younger) layers would not fit deeper (older) layers. Such a pattern is not clearly evident in figures 6 and 8, suggesting that either the migration rate and accumulation pattern have been relatively constant in time or that isochrone shapes are not sensitive lo past changes that have occurred.

The insensitivity of isochrone shapes to past flow changes arises because isochrone deformation is an integrated process. Layer shapes lose evidence of old flow regimes over time. It is therefore difficult to determine how long the divide has been migrating. We assume that it has been migrating for t>>H/b r 10⁴ a). However, the actual onset of migration may be relatively recent and the layer shapes no longer contain onset information. It is not correct to use the age of the oldest (deepest) detectable isochrone with a displaced apex as the age of divide-migration onset, because the moving flow field affects the deformation history of all isochrones, regardless of their age. However, a deep layer can be used to place an upper bound on the date of migration onset. If the divide had been in one location for a long time and then started in migrate, we would expect to see remnant warping at the old divide position in the deeper layers where the warping is expected to be most pronounced. It would take ~10⁴ years for the changing flow field associated with divide migration to smear out well-developed warped layers (paper in preparation by E. Wad-dington and others). Suppose the divide was once at the present apex of the deepest detected layer, about 700 m south of the divide. Then the divide has been migrating toward Ice Stream D for at least 700/m years, or about 1.5-15 ka.

Migration of the ice divide suggests that non-steady processes are affecting the geometry of Siple Dome. Candidates include changes in the elevation or position of the bounding ice streams or a change in the spatial accumulation pattern. Reference Nereson, Hindmarsh and RaymondNereson and others (1998) have explored these possibilities and concluded that relatively small changes in either candidate could produce the rate of observed divide motion. A further consequence of the Siple Dome divide migration and the associated effects on isochrone shapes is that pre-Holocene ice is thickest about 0.5 km south of the present ice divide. This information has been used to select the exact site for a deep ice core to be obtained for paleoclimate analysis.

The migration rate is sufficiently slow that the near-divide warping is maintained in the internal layer pattern. This indicates that the divide has been within a few ice thicknesses of its present position for the past 10³ 10⁴ years and suggests that there have been no major asynchronous changes in the configuration of the bounding ice streams or in the geometry of Siple Dome in that time. Significant asymmetric forcing at the boundaries of Siple Dome would cause rapid motion of the divide ι(Reference HindmarshHindmarsh, 1996) and prevent the development of a near-divide isochrone bump.

Since the development of a divide bump takes time, its presence at Siple Dome raises the question about the age of Siple Dome itself. There are two distinct hypotheses. The first is that Siple Dome has existed in about its present state for the past 10⁴ or more years. This suggests that Siple Dome has been a flow center for most of the Holocene and supports the theory that ice streams were present during the early stages of retreat of the Ross Ice Shelf (Reference Denton, Bockheim, Wilson and StuiverDenton and others, 1989). The second hypothesis is that Siple Dome formed more recently as a result of ice-stream initiation and thinning of the WAIS.

The analysis presented here assumes the first hypothesis: constant ice thickness. Given this assumption and a time-averaged accumulation rate of 0.10 m a⁻¹ ice equivalent, the time-scale to create a well-developed warp in the internal layers by either an accumulation low or by ice-dynamic processes is related to the thickness/accumulation rate time-scale (inverse of the characteristic vertical strain rate), about 10⁴ years. Without considering the smearing effects of divide migration, it would take at least 5 ka to deform isochrones to the extent we observe at Siple Dome. The first hypothesis is thus consistent with the data.

However, we cannot eliminate the second hypothesis. Possible rapid thinning of WAIS and the formation of Siple Dome would be associated with large vertical strain rates which increase the speed of ice deformation and reduce the time required to create internal layer warping beneath the divide. If Siple Dome formed by rapid synchronous thinning at its boundaries as a result of ice-stream initiation, then it is possible to form isochrone bumps beneath the divide at all depths in as little as 10 years. For example, suppose Siple Dome was once twice as thick as at present and thinned to its present thickness in about 10� years. Assuming no change in accumulation rate, the vertical strain rate during thinning would be roughly 10⁻³ a ⁻¹, which is about one order of magnitude larger than the vertical strain rate for the steady-state case. The time to form a divide bump, which is related to the vertical strain rate, would correspondingly be reduced by about a factor of ten. Measurements of age vs depth and total gas content from the deep ice core will help resolve this question. If Siple Dome has thinned rapidly in the recent past, then the thickness of layers deposited during and prior to thinning would be smaller (due to larger vertical strain rates). Thus, these layers would appear in the present ice sheet at shallower depths than if no thinning occurred. Accumulation-rate estimates from layer-thickness profiles which assume no dome thinning would be artificially low.