Introduction

A recent work (MacAyeal and others, 1996) has shown that accuracy in modelling the velocity of the Ross Ice Shelf, Antarctica, is limited by lack of knowledge of the large-scale variation of rheological parameters in the ice shelf. Rheological parameters are generally incorporated in ice-shelf models through specification of an ice-column average of stress-dependent viscosity based on Glen’s flow law. This viscosity accounts for variations of temperature in the ice shell, but it does not account for density variation, chemical contaminants (e.g. sea salt), ice-crystal fabrics or other factors which influence ice creep (and which are not represented in Glen’s flow law). In addition to these shortcomings, it is not obvious a priori whether Glen’s law, which is laboratory-based, can represent the space and stress scales of deformation occurring in situ in the natural ice-shelf environment.

In the present study, we derive the spatial distribution of depth-averaged effective viscosity in the Ross Ice Shelf from velocity and thickness measurements under the assumption of isotropic ice. By the expression “effective viscosity”, we mean a spatially variable scalar that relates the local deviatoric stress to the strain rate. The spatial variation of this scalar is presumed to reflect the flow law which actually applies to the ice shelf as well as the variation of physical ice conditions such as temperature. Data used here have been described in MacAyeal and others (1996) and may be accessed at the World Data Center for Glaciology A, Boulder.

Our derivation of effective viscosity variations is motivated by three goals. The first is to see if the relationships between the strain rates, surface temperatures and effective viscosities are in agreement with Glen’s flow law. The second is to detect regions with strong basal melting or refreezing. The third is to determine whether temporal ice-stream discharge variability may have left a signature in the rheological properties of the ice shelf.

Notation

Control Method

The purpose of the control method is to determine the best fit between an observed field and its modelled counterpart through adjustment of model parameters. The resulting solution, i.e. values of the adjusted parameters, is not necessarily unique and could depend on the initial guess of the model parameters. In our study, there is one adjusted model parameter: the depth-averaged viscosity η (c,y). The misfit between observation and model used to constrain this adjusted parameter is evaluated by a least-squares definition of velocity error:

(1)

The equations that link the depth-averaged viscosity to the velocities, and which form the basis of the ice-shelf model, are the diagnostic equations of ice-shelf motion (Morland, 1987; MacAyeal, 1989):

(2)

These diagnostic equations, which form constraints on the minimization of J, come from a vertical integration of the Stokes’ equations in the ice shelf. The viscosity appearing in Equations (2) represents a vertically averaged quantity. The main underlying assumptions which justify Equations (2) are that acche smallnseleration terms and inertia are negligible, that horizontal velocities are independent of depth (because ofe ts of . we dthe aspect ratio in an ice shelf) and that ice is isotropic. No other assumption is made about the flow law (i.eo not assume that the viscosity follows Glen’s flow law). Equations (2), together with boundary conditions, describe what we refer to as the “forward model”.

A useful means to minimize J as defined in Equation (1) is to use a Lagrange multiplier vector (MacAyeal, 1993) to incorporate the constraints of Equations (2) into the definition of J as follows:

(3)

Note that the second and third integrals on the right-hand side are null when u and v are solutions of the forward model. Mathematical manipulation, primarily the use of the Green–Ostrogradski theorem on the null terms, allows us to differentiate J with respect to the free model parameter η , This derivative serves two purposes: its smallness determines when the minimum of J is approached, and it provides useful information for improving previous estimates of η (see next section). The expression for this derivative is obtained from the variation of J with respect to η .

(4)

Fig. 1. Flow chart of the control method.

The Lagrange multiplier vector components, referred to above, are obtained from the following equations:

(5)

with λ = μ = 0 specified on the boundaries of the study area (Γ). Observe that the model/observation misfit, u _d – u and v _d – v, serve as forcing to the equations which determine λ and μ. The similarity between the constraints on the Lagrange multiplier vector represented by Equations (5) (also called the adjoint equations), and the diagnostic-equations of the “forward model”, represented by Equations (2), is advantageous since the same kind of numerical algorithm can be used to compute either the velocity or the Lagrange multiplier vector components.

Steepest-Descent Algorithm

Numerically, the spatial distribution of viscosity is described by discrete values η _i,j at specified grid points. This means that the space-dependent function to be adjusted η (x,y) reduces to a multi-parameter variable. An initial guess of η _i,j is needed to initiate the algorithm of minimization (sec Fig. 1). Using this initial guess, Equations (2) (with appropriate kinematic and dynamic boundary conditions) are solved to obtain the velocity components u and v. The finite-difference methodology used to solve Equations (2) is described by Rommelaere and Ritz (1996). Once u and v are determined, the misfit model/observation velocity is evaluated and used to compute λ and μ from Equations (5). Finally,

Fig. 2. (a) True viscosity field in MPa a of the idealized ice shelf used to evaluate accuracy, (b) Reconstructed viscosity field obtained by control method.

once λ and μ are obtained, Equation (4) is evaluated to determine whether the minimization condition δJ = 0 has been achieved or a new viscosity field that will reduce the model/observation misfit is to be found by searching along the direction of the gradient of J, e.g.

(6)

The problem is now reduced to a simple univariate minimization with respect to the scalar Δα (Δα is positive and has the dimension of a viscosity). This minimization is performed by a method combining golden section search and successive parabolic interpolation. The algorithm is run until a satisfactory convergence on the modelled velocity components is achieved.

Accuracy

It will never be possible to exactly fit the observed velocity field with our model, for two reasons. First, the assumptions used to justify Equations (2) may be wrong. An example of such a problem, in our case, is the assumption of isotropic ice which may be wrong in places where ice has experienced strong shear stress and has reoriented its crystalline axes. A second reason preventing an exact fit is error introduced in the measurements. In our study, these errors are estimated to be about 30 m a⁻¹ (Thomas and others, 1984).

To assess the influence of measurement error on the significance of the derived viscosity field, we conducted a series of “identical twin” experiments involving an idealized, rectangular ice shelf of 610 km × 710 km (Fig. 2), for which we prescribe both the ice thickness and viscosity. In the idealized ice-shelf geometry the ice front is located on one boundary of the domain. On the three other sides, the velocity is set equal to zero. The ice thickness decreases linearly from 850 m on the upstream boundary to 250 m at the ice front (lower boundary).

We contaminate the solution of the forward problem u and v with random observational noise of roughly the same amplitude as noise in the RIGGS data. This “contaminated” solution is then used as the “observation” from which η _e is estimated using the steepest-descent algorithm described previously. Differences between the prescribed η and the estimated η _e provide an experimental quantification of the

Fig. 3. Fig. 3 (a) “Exact” velocity magnitude in m a⁻¹ computed with the true viscosity (Fig. 2a) and a model which solves Equations (2), (b) “Measured” velocity magnitude in m a⁻¹ in which the exact velocity field (a) is contaminated with a white Gaussian noise.

uncertainly of η due to random observational error in u and v (see figs 2 and 3).

The steepest-descent algorithm is run from an initial spatially uniform viscosity of 25 MPa a (789 TPa s). Our results (Fig. 2b) show how noise introduced to the velocity input is passed to the estimated viscosity field. Variations of about 10% in the area where the viscosity should be constant are observed. Nevertheless, it is encouraging to see that the main features of the viscosity are detected (Fig. 2b). The identical-twin experiment suggests that the relative error is less than 20% everywhere. This level of precision and the detectability of main features are sufficient to study viscosity variations on real ice shelves in a qualitative manner.

Rheology of the Central Part of the Ross ICE Shelf

The input data for the Ross Ice Shelf viscosity analysis include a digitized map of ice thickness (Bentley and others 1979; MacAyeal and others, 1987) and the measured horizontal surface velocities (Thomas and MacAyeal, 1982; Thomas and others, 1984; MacAyeal and others, 1987). The ice thickness has heen reduced by 14 m (estimated using the ice-density profile of station J9 on the Ross Ice Shelf) to account for the amount of air in the ice and firn. The ice front is set to its measured position of 1972. One can notice on the maps of the Ross Ice Shelf (figs 4 and 5) that data east of Roosevelt Island have not been considered. We find that the dynamics of this excluded region strongly depends on boundary conditions along the sides of Roosevelt Island. Preliminary study revealed that we lacked sufficient understanding of these boundary conditions to reduce the misfit between observed and modelled velocities to a level consistent with misfit achieved elsewhere. The domain (Γ) on which we shall focus is thus restricted to the central part of the ice shelf.

The study domain (Γ) is a 267 000 km² area of the Ross Ice Shelf. The spatial resolution of the finite-difference grid is 6.822 km, so the numerical domain consists of 5738 active grid points. The performance index J is minimized by selecting the 5738 values η _i,j at the i, jth grid point, using the steepest-descent algorithm described above. The initial

Fig. 4. (a) Comparison of the observed and modelled velocity magnitude at each RIGGS station of the domain. (b) Measured velocity magnitude (ma⁻¹) in the central part of the Ross Ice Shelf, (c) Best fit of the velocity magnitude (m a⁻¹) obtained by the control method.

guess of η is constant at 20 MPa a (631 TPa s). Several different initial guesses were tried, and the final viscosity pattern (Fig. 5) does not seem to vary significantly as a result.

Figure 4 shows a comparison between the magnitude of the observed velocity field and the best fit obtained by the control method. The general flow pattern seems to be well reproduced, and what appears to be noise in the measurements is spatially filtered by the forward model to yield a smoother model velocity field. There are several RIGGS stations (Thomas and others, 1984) in the domain where the point-by-point comparison of modelled and observed flow shows substantial misfit. Apparently, this disagreement appears more pronounced for the velocity direction than for the velocity magnitude. This illustrates the fact that viscosity variations can easily speed up or slow the ice flow, but its influence on the direction of flow is far less important. The velocity direction seems to be determined almost entirely by the geometry and the ice-thickness field of the ice shelf.

The viscosity pattern does not show much variation, compared with what has been measured in grounded ice sheets or in laboratories. The minimum value here is only about one-third of the maximum value, and the area average η of the ice shelf is about 30 MPa a (946 TPa s). According to the most commonly used parameters of Glen’s flow law with an exponent equal to three, the average viscosity we derive corresponds to an isothermal ice at −25°C with a second strain-rate invariant of about 10⁻³ a⁻¹ (which is the order of magnitude of that observed on the Ross Ice Shelf). The viscosity determined by the control method is related to the vertically averaged temperature–depth profile. The “isotherm-equivalent” temperature estimated above should logically lie between the basal and the surface temperatures of the ice shelf. The value of −25°C is very close to the surface temperature, and this suggests that ice-shelf models including a rigorous treatment of temperature variation with depth and Glen’s flow law would predict a softer and thus faster Ross Ice Shelf than observed.

Discussion

An ice shelf that is freezing at its base should appear less viscous than where its base is melting, because basal freezing tends to warm the temperature—depth profile and to entrap salt which softens ice. For similar reasons, an ice shelf that has strong surface accumulation should appear more viscous. Most of the variations in our estimate of η are probably due to such effects, However, a quantitative evaluation of basal melting or refreezing is not possible from ice-flow considerations alone. To further interpret our results, knowledge of temperature profiles and information about ice-crystal fabric is also necessary.

Several regions identified from our estimate of η are of special interest (Fig. 5). Downstream of Crary Ice Rise, there is an area of soft ice that experiences relatively low stresses. This suggests a very strong temperature effect due to basal freezing. Such a phenomenon has already been mentioned for instance by MacAyeal and Thomas (1986) and by Bindschadler and others (1990). Indeed, many glaciological measurements have been carried out in this area, and it is generally agreed that Crary Ice Rise is in the transition zone between basal melting (upstream) and basal accretion (downstream). Furthermore, many crevasses are observed in this zone, and this may also soften the ice.

The variation of η downstream of Steershead Ice Rise is a bit more difficult to explain. We suspect that the signal of high viscosity is accentuated by local and very weak ice-shelf grounding (e.g. ice rumples) that is not accounted for in the forward model, In addition, measurement coverage near Steershead Ice Rise is far less dense than near Crary Ice Rise. However, our results suggest a difference in the influence of the two ice rises on the large-scale flow, and this difference is unexplained. Mean surface temperature is colder near Steershead Ice Rise (−28°) than near Crary Ice Rise (−25.5°C), but this is not sufficient to interpret the viscosity variations. Part of the difference could be explained

Fig. 5. Fig 5. Depth-averaged viscosity of the central Ross Ice Shelf reconstructed by the control method (MPa a). Three stream lines (dashed lines) are plotted to show the possible origin of the west-central region’s maximum of η (178° E, 10° S).

by temperature effects of basal melting. In places where the water column is relatively thin, typically close to the grounding line, heat can be provided to the ice shelf by vigorous tidal mixing (MacAyeal, 1985; Jenkins, 1991), and this initiates the melting of these areas. Tidal resonance phenomena at specific sites can amplify tidal effects: this may be the case for Steershead Ice Rise.

The ice-front region is melting according to oceanographic observations and should thus be composed of stiffer ice. Our results show the opposite, however. Ice may be softened near the ice front because of crevasses and cracks induced by flexural-gravity wave activity.

Finally, the most prominent and presumably the most interesting feature of the overall pattern is the maximum of viscosity in the west-central part of the ice shelf. Since this location is not tied to coastal features such as ice rises, we suggest that the high η may be a property of local ice that is advected with the flow. There is no mention of such an “anomaly” of viscosity in the literature, but it may be associated with the significant systematic bias MacAyeal and Thomas (1986) found between their model and the observations in the southwestern part of the Ross Ice Shelf. Indeed Equations (2) which determine ice-shelf flow as a function of ice-thickness and boundary conditions are elliptic in two directions, and are thus non-local: a difference of viscosity may induce velocity differences 300–400 km upstream or downstream and no differences locally. Estimated flowlines suggest that this high η signal may be associated with ice which came from Ice Stream A or B. Moreover, according to the mid-shelf location of the viscosity maximum, this ice would have been discharged from the ice stream at least 1000 years ago. The signal of high η in the west-central ice shelf is not necessarily due to theological properties of ice. Buoyancy anomalies, which would induce changes in the driving stress, could lead to the same flow conditions and fool our method into suggesting a high η . In any case, either a decrease of the driving stress or an increase of viscosity appears in the west-central ice shelf. Casassa (1993), with the help of AVHRR satellite imagery, has suggested that discharge of Ice Stream A increased relative to that of Ice Stream B about 1000 years ago. The η maximum in our results may thus be an additional sign of the same phenomenon.

Concluding Remarks

We have developed a method to detect large-scale viscosity variations in ice shelves. Its application on the Ross lee Shelf shows that the viscosity is of the same order of magnitude as predicted by Glen’s flow law with an exponent equal to three and a temperature-dependent coefficient comparable to that associated with the observed surface temperature. For very simplified modelling approaches of the Ross Ice Shelf, a constant viscosity of about 30 MPa a (946 TPa s) is able to catch the main characteristics of large-scale flow (such an approach is of course not recommended for simulations of ice-shelf evolution under changing climates).

Our method has the potential to reveal the local influence of unusual basal melting or freezing, and where ice viscosity is influenced by ice-shelf history. For instance, the west-central part of the Ross Ice Shelf possesses stiffer ice, and this is possibly linked with changes in discharge of Ice Stream A or B at least 1000 years ago. Finally, we recommend that this method be applied more extensively on the other large ice shelves of Antarctica to explain differences in their dynamics.