1. Introduction
Radio-echo sounding (RES) has been used to investigate the subsurface properties (internal structure) of ice over large areas of the Antarctic and Greenland ice sheets. Internal layering from RES measurements (Fig. 1) provides a vertical profile in an ice sheet, which could previously be obtained only from drilled boreholes. The radar echoes arising from internal layers in ice are caused by the sudden changes in complex dielectric permittivity of layers in the ice sheets (Reference FujitaFujita and others, 1999). The mechanisms causing the sudden changes have been examined in detail, and changes in shallow layers have been attributed primarily to ice-density variations. Changes in deeper layers have mainly been assigned to changes in electrical conductivity due to acidic fallout from volcanic eruptions and/or changes in impurity concentration associated with climatic transitions (Reference HarrisonHarrison, 1973; Reference Paren and RobinParen and Robin, 1975; Reference CloughClough, 1977; Reference HammerHammer, 1980; Reference MillarMillar, 1981; Reference MooreMoore, 1988; Reference Fujita and MaeFujita and Mae, 1994; Reference FujitaFujita and others, 1999). Variation of ice-crystal orientation can also cause the changes in permittivity at the large depth of the ice sheet (Reference HarrisonHarrison, 1973; Reference Fujita and MaeFujita and Mae, 1994; Reference FujitaFujita and others, 1999).
Fig. 1. RES imagery showing internal layers along the airborne flight-line indicated by dashed line in Figure 2.
RES-detected internal layers are widely recognized as isochrones: former ice-sheet surfaces that have been buried and deformed by ice flow (Reference WhillansWhillans, 1976; Reference CloughClough, 1977; Reference HammerHammer, 1980; Reference Dahl-JensenDahl-Jensen and others, 1997; Reference Morse, Waddington and SteigMorse and others, 1998; Reference Nereson, Raymond, Jacobel and WaddingtonNereson and others, 2000). As such, they would contain information about the ice-sheet history, from which inferences can be made about some of the climate changes. For example, known ages of internal layers enabled past accumulations to be inferred from the horizontal variation in thickness between the layers on Taylor Dome, East Antarctica (Reference Morse, Waddington and SteigMorse and others, 1998), and around Summit, Greenland (DahlJensen and others, 1997; Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others, in press). Lacking time information, Reference Nereson, Raymond, Jacobel and WaddingtonNereson and others (2000) estimated the pattern of accumulation over Siple Dome, West Antarctica, from the observed pattern of internal layering.
The most common method of dating the internal layers is to combine them with annual-layer thickness measurements from the associated ice cores. A recently developed system (Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others, in press) traces the internal layers from RES imageries and then extends the ages of the layers at the Greenland Icecore Project (GRIP) ice-core site at the Summit of the Greenland ice sheet, based upon the available GRIP age–depth relationship (Reference JohnsenJohnsen and others, 1997), to the same layers along the airborne flight-lines. The dated internal layers mainly cover the northern areas of the Greenland ice sheet (W. Abdalati, unpublished information), providing a good reference to validate the numerical model calculation. In regions where the internal layers cannot be detected or cannot be extended from the dated ice cores (e.g. the ablation zone), numerical modeling may give a good approximation of the ice age at any ice-sheet depth.
In this paper, we apply the anisotropic-ice flowline model developed by Reference Wang and WarnerWang and Warner (1999) to a flowline through Swiss Camp (69.57° N, 49.28°W), West Greenland, to estimate the dates of internal layers obtained by the RES technique (Reference Fahnestock, Abdalati, Luo and GogineniFahnestock and others, in press). The available measurements of dated internal layers are used to validate the model calculation.
2. The model
The model used here is an anisotropic steady-state ice flow-line model (Reference Wang and WarnerWang and Warner, 1999). The stress–strain-rate relationship for anisotropic ice is characterized by an enhancement factor, defined as the ratio of tertiary strain rate for anisotropic ice to minimum strain rate for isotropic ice, based on the laboratory observations of ice deformation under combined compression and shear stresses (Reference Li, Jacka and BuddLi and others, 1996). The model has been described in detail by Reference Wang and WarnerWang and Warner (1999). Here we give a brief review of the flow relations used in the model to describe anisotropic-ice flow.
In this two-dimensional model, the horizontal transverse flow is neglected so that the flow is entirely constrained to the vertical compression and horizontal shear along the flowline in steady-state balance. If x denotes the direction of flow and z the vertical direction, the flow relations between strain rate and stress for components of shear and compression are
where 
is shear strain rate, is compressive strain rate, A0(T) is a temperature-dependent parameter, T is temperature, τo is octahedral shear stress, τxz is shear stress, 
is compressive factors of deviatoric shear and stress, compression G(λc) and component, F (λc) are enhancement respectively, and λc is a compression factor defined as
Laboratory experiments indicate that G(λc) and F(λc) can be simplified to be equal to E(λc) (Reference Li, Jacka and BuddLi and others, 1996) according to:
where Es and Ec are respective enhancement factors for shear or compression alone. While λc varies from 1 to 0 as the stress situation varies from purely confined compression (shear stress τxz = 0) to simple shear (compressive deviatoric stress enhancement factor E(λc) increases from 3 to 10.
In Equations (1) and (2), the temperature-dependent coefficient A0(T) is based on the laboratory experiment results (Reference Budd and JackaBudd and Jacka, 1989) tabulated by Reference Wang and WarnerWang and Warner (1998, table 2) which have similar values calculated from the Arrhenius relation (Reference PatersonPaterson, 1994). The shear stress is taken as the driving stress
in terms of the ice density ρ, the acceleration due to gravity g, the depth Z and the surface slope α. The octahedral shear stress is taken as
assuming the ice flow corresponds to a confined vertical compression stress combined with a horizontal shear stress.
After rearranging the above equations, a cubic equation for the shear strain rate, τ xz is obtained as
This equation involves vertical compressive strain rate, τ z which can be calculated from horizontal velocity based on the assumption that ice is incompressible.
Equation (7) is the stress–strain-rate relation for anisotropic-ice flow used in the flowline model and solved iteratively for the shear strain rate (see Reference Wang and WarnerWang and Warner, 1999).
3. Application Of The Model To The Flowline
Flowlines over the whole Greenland ice sheet were generated based on the assumption that the ice flows downslope in the direction perpendicular to the surface contours using 5 km gridded surface elevations derived from European Remote-sensing Satellite (ERS-1) radar altimeter data (Reference Zwally, Brenner and FuZwally and Brenner, 2001). The flowline studied in this paper (Fig. 2), from the ice-divide ridge through the Swiss Camp station northwestward to the coast, was chosen because an airborne radar survey ran very close to the flow-line. The internal layers were traced from RES imagery (Fig. 1) along the flight-line and were dated by extending the same layers from the Greenland Icecore Project (GRIP) ice-core site where dating was available, which were used to validate the model calculation. Furthermore, the previous filed observations, such as surface velocity (Reference HofmannHofmann, 1975), accumulation (Reference BensonBenson,1962) and temperature (Reference Mock and WeeksMock and Weeks, 1966), were available near this flowline and provided the input data to the model.
Fig. 2. Map of flowlines in West Greenland, with the contours of surface elevation in meters. Thick line shows the studied flowline through Swiss Camp station; dashed line shows the airborne flight-line; dotted line shows the Expéditions Glaciologiques Internationales au Groenland (EGIG) traverse-line surface velocity observations (Reference HofmannHofmann, 1975); solid circles show locations of the measured accumulation (Reference BensonBenson, 1962) and temperature (Reference Mock and WeeksMock and Weeks, 1966); and asterisks show the sites of the GRIP borehole and Swiss Camp station.
Table 1. Ages of the internal layers displayed in Figures 3a and 4
The Swiss Camp (69.57°N, 49.28°W), located in west-central Greenland, was established near the equilibrium line in 1990, estimated from surface balance measurements of the 1980s. The surface mass balance has been studied by means of global positioning system measurements. The modeling study along the flowline through Swiss Camp provides information about the internal dynamics of ice in the vicinity of the Camp.
Figure 3 shows the data profiles used as inputs to the flowline model. Swiss Camp station is 510 km from the ice divide along the direction of flow.
Fig. 3. Profiles of input data used in the model, with the location of Swiss Camp station indicated by asterisk. (a) Flowline topography (smoothed curves) used in the model, the measured topography from airborne flight (rough curves) and internal layers (thin lines) representing the ages in Table 1. (b) Surface accumulation rates: dots indicate previous measurements (Reference BensonBenson, 1962) at the locations shown in Figure 2, the dashed line is interpolated from 50 km grid database (Reference Zwally and GiovinettoZwally and Giovinetto, 2000) and the solid line is the trend of profile used in the model. (c) Surface temperatures: dots indicate previous measurements (Reference Mock and WeeksMock and Weeks, 1966) at the locations shown in Figure 2, and solid line is interpolated from 50 km grid database (personal communication from M. B. Giovinetto, 2001) used in the model. (d) Surface horizontal velocities: dots indicate the measurements (Reference HofmannHofmann, 1975) along EGIG traverse line shown in Figure 2, asterisk shows the observation at Swiss Camp site and the solid line is the profile used in the model.
The ice-sheet topography along the studied flowline is shown in Figure 3a. Surface elevations used in the model were taken from 5 km grid satellite radar altimeter data (Reference Zwally, Brenner and FuZwally and Brenner, 2001), and bedrock elevation by subtracting ice thickness (Reference Bamber, Layberry and GogineniBamber and others, in press) from the surface elevation. Surface and bedrock profiles obtained from RES measurements along the flight-line are plotted in Figure 3a in the section with the observed internal layers. Comparisons with the RES measurements show the irregular basal topography, but the profile from the 5 km grids appears reasonably to represent the smoothed topography incorporated into the model.
Twelve internal layers are traced from RES imagery (Fig. 1) and displayed in Figure 3a. The discontinuous internal layers 410–430km from the ice divide are due to unclear RES imagery in this section, but the pattern is then resumed beyond 430 km. The corresponding ages of the layers shown, obtained by tracing the continuous layers back to the dated GRIP core site, vary from 2339 to 12329 years (see Table 1).
Internal layers along the flowline seem to reflect the shape of the bedrock undulation with decreasing amplitude as the distance from the bedrock increases. This can be seen more clearly in Figure 1. Near the Greenland Summit along the ice-divide ridge, however, where the bedrock is smooth, the internal oscillations are caused by variations of the dynamic velocity fields rather than the bedrock undulations (Reference Dahl-JensenDahl-Jensen and others, 1997).
Surface accumulation rate along the flowline used in the model is a trend line shown in Figure 3b. This was determined based on the previous observations (Reference BensonBenson, 1962) at the locations near the flowline (see Fig. 2) and 50 km grid database (Reference Zwally and GiovinettoZwally and Giovinetto, 2000). The accumulation rate is zero at Swiss Camp station on the equilibrium line and is negative at the lower elevations where ablation exceeds precipitation.
The surface temperature profile, shown by the solid line in Figure 3c, was interpolated from a 50 km grid database (personal communication from M. B. Giovinetto, 2001), showing a good agreement with the previous observations (Reference Mock and WeeksMock and Weeks, 1966) near the flowline (see Fig. 2).
The surface velocity profile (Fig.3d) is assumed to be the observations along the EGIG traverse line (Reference HofmannHofmann, 1975), which runs close to the flowline studied in this paper (see Fig. 2), based on the assumption that surface horizontal velocities along the flowlines vary slowly along the surface elevation contour lines. Using laboratory-based flow relations results in an overestimate of surface velocity (~8 times too high) by integrating shear strain rate down to the bed-rock. This is because the model does not consider the reduction of enhancement and shear stress near the bedrock. This reduction may cause the reduced shear strain rates near the bedrock which have been found from several of the borehole inclination measurements in Antarctica (e.g. Reference Russell-Head and BuddRussell-Head and Budd, 1979; Reference EtheridgeEtheridge, 1989; Reference Morgan, Ommen, Elcheikh and JunMorgan and others, 1998). In Greenland it has been found at Dye 3 borehole (DahlJensen, 1985; Reference Dahl-Jensen and GundestrupDahl-Jensen and Gundestrup, 1987) that maximum shear strain rates occur at the bottom due to a high concentration of dust and other impurities, but at 200–257m depth above the bed the shear strain rates are almost constant along with the reduction of the enhancement factors. To compensate for the overestimated velocity without turning any parameter in the flow law, we adopted a reasonably simple scheme from Reference Wang and WarnerWang and Warner (1999) by terminating integration of shear strain rate near the bedrock to match the measured surface velocity.
Basal sliding velocity was estimated based on the study of the ice-sheet modeling along a flowline jointed with the central flowline of Jakobshavn Isbræ drainage basin (see Fig. 2) (Reference Funk, Echelmeyer and IkenFunk and others, 1994) using the Reference Huybrechts and WoldeHuybrechts and de Wolde (1999) equation
where a constant is basal shear stress and Z ¤ is the height above buoyancy.
All input data were interpolated so that the horizontal resolution was 1km. A rescaled vertical coordinate was used by subdividing the ice thickness into 100 evenly spaced bands. The basal temperature gradient of 0.02°Cm–1 (Reference Kostecka and WhillansKostecka and Whillans, 1988; Reference Funk, Echelmeyer and IkenFunk and others, 1994) was used in the calculation of the steady-state temperature.
4. Results And Discussion
4.1. Comparison of modeled isochrones with observed internal layers
The ages of the ice at any depth of the ice sheet along the flowline were determined by numerically integrating the velocity fields following the ice-particle trajectories from the surface, using the present-day surface accumulation rates and assuming a steady-state ice flow.
Several isochrones calculated from the model are displayed in Figure 4 to compare with the observed internal layers. Using the smoothed bedrock and a trend surface accumulation rate results in smoothed modeled isochrones which do not capture the high-frequency oscillations of the observed internal layers, but agree quite well on larger scales.
Fig. 4. Isochrones (solid lines) calculated from the model, and observed internal layers (dashed lines) traced from RES imagery (Fig. 1). Surface and bedrock elevations used in the model are plotted in smoothed curves, and the RES observations in rough curves.
The model simulating isotropic-ice flow has been run by replacing the flow relation for anisotropic ice with the relation for isotropic ice, i.e. using Equation (7) with E(λc) = 1 in the model calculations. Under all of the same conditions, the isochrones calculated from the model for isotropic-ice flow do not match the observed internal layers well (see Fig. 5 and 6a).
Fig. 5. Plot of the ages computed from the models against observed ages from the layer measurements. Open triangles and open circles indicate mean values of computed ages from the models for anisotropic and isotropic ice, respectively, along each measured layer. Bars show the range of variations between computed and observed ages along each layer. The dashed line indicates the ideal match situation. importance of including anisotropy in ice-flow modeling.
Fig. 6. (a) Relation of age and depth, at 200 km from the ice divide, generated from the models for anisotropic-ice flow (solid line) and using Equation (7) with E(λc) =1for isotropic-ice flow (dashed line), and from the observed internal layers (solid circles). (b) Surface accumulation rate used in the model (thick line) and reconstructed from the observed internal layer nearest the surface of the ice sheet with 2339 years (thin line).
Figure 5 gives a quantitative expression of the ages computed from the models for both anisotropic and isotropic ice with the ages obtained from GRIP cores for 12 observed internal layers. The modeled ages at the depths of the observed horizons are calculated as a function of x and shown in the bars. Themean values are obtained by averaging those modeled ages along each observed layer. The larger error bars in lower layers are caused by the bedrock undulation, which will be substantially reduced by using flight-line bedrock (rough curve in Fig. 4) in the models. Figure 5 shows the
4.2. Relation of age and depth, and reconstructed surface accumulation rates
The comparison of the relations of age and depth derived from the model results with the observations is shown in Figure 6a for a point 200 km from the ridge. All three profiles in Figure 6a closely agree for approximately the top one-third of the ice thickness. At greater depth, the modeled profile for anisotropic ice still follows the observations closely, in contrast with the obvious departure shown by the dashed line for isotropic ice, which indicates that isotropic ice is no longer appropriate when anisotropy of the ice-crystal fabric develops with increasing depth and that it is important to consider the effect of the ice-crystal anisotropy on the modeling of ice flow in the deeper parts of ice sheets.
Based on the assumption of steady-state ice flow, we reconstructed the accumulation rates for the internal layer nearest to the surface of the ice sheet with 2339 years. The age–depth relation given in Figure 6a and the modeled vertical strain were used in the calculation. The reconstructed accumulation rates are close to the present-day values (see Fig. 6b), confirming that this part of the Greenland ice sheet has not experiencedmajor changes in accumulation rate and ice-flow status.
4.3. Ice flow at Swiss Camp station at/near the equilibrium line
In this study, present-day surface accumulation rate, surface temperature, surface elevation, ice thickness and surface velocity were used as model inputs. The ice-flow regime through the depth was modeled along the flowline. The agreement between the modeled isochrones and the observed internal layers implies that the model incorporated with ice-crystal fabric anisotropy has given a proper description of the ice flow. Here, we use Swiss Camp site, at or near the equilibrium line, as an example of the modeling results. The mass balance of the Greenland ice sheet at Swiss Camp site will be studied based on these results and the surface observations.
Vertical profiles of several outputs from the Swiss Camp model are summarized in Figure 7. The parameters shown are: shear and vertical compressive strain rates, horizontal and vertical velocities, enhancement, stresses, temperature and age.
Fig. 7. Vertical profiles of several outputs from the model at Swiss Camp station, 510 km from the ice divide: (a) shear strain rate (ɛ̇xz ); (b) horizontal velocity (Vx ); (c) enhancement factor (E(λc)); (d) temperature (T); (e) vertical compressive strain rate (ɛ̇z ); (f) vertical velocity (Vz ); (g) shear stress (τxz ), compression stress deviator and octahedral shear stress (τo ); and (h) age.
Shear strain rates (ɛ̇xz) increase with depth, which is accompanied by an increased enhancement (E(λc)) from 3 to 10, as the range of stress regimes changes from predominantly vertical compressive stress (τz) in the upper layers to predominantly horizontal shear stress (τxz) near the bed. The transition layer from the compression stress dominance to the shear stress dominance is about 300 mdeep (Fig.7g).
It is important to note the following points:
(1) Since accumulation rate is zero at the site, the vertical velocities (Vz) represent the total horizontal advections.
(2) Temperature profile (T) shows that the bottom ice reaches the melting point, indicating the existence of basal sliding.
(3) The age–depth profile (Fig. 7h) shows that the ice near the bottom is >20 000 years old.
5. Conclusion
The application of an anisotropic steady-state ice flowline model to the flowline in West Greenland is presented in this paper. A close agreement is found between the isochrones generated from the model using present-day input data and observed internal layers with confirmed dates, and this strong agreement indicates that this part of the Greenland ice sheet is essentially in steady state. This part of the ice sheet has been close to its present form for a period of at least 12 000 years. This successful comparison shows that the model provides a good approximation of the flow of the ice sheet, and that the ice-crystal fabric anisotropy must be considered in order to accurately model deeper ice flow.
Acknowledgement
We wish to thank M. Beckley for the help with the graphics.
