Numerical simulations of the ice flow dynamics of George VI Ice Shelf, Antarctica

Abstract

A diagnostic, dynamic/thermodynamic ice-shelf model is applied to the George VI Ice Shelf, situated in the Bellinghausen Sea, Antarctica. The George VI Ice Shelf has a peculiar flow geometry which sets it apart from other ice shelves. Inflow occurs along the two longest, and almost parallel, sides, whereas outflow occurs on the two ice fronts that are relatively short and situated at opposite ends of the ice shelf. Two data sources were used to derive the ice thickness distribution: conventional radioecho sounding from the British Antarctic Survey was combined with thickness inferred from surface elevation obtained by the NASA GLAS satellite system assuming hydrostatic equilibrium. We simulate the present ice flow over the ice shelf that results from the ice thickness distribution, the inflow at the grounding line and the flow rate factor. The high spatial resolution of the ice thickness distribution leads to very detailed simulations. The flow field has some extraordinary elements (e.g. the stagnation point characteristics resulting from the unusual ice-shelf geometry).

Introduction

Ice shelves are elements sensitive to the climate and govern, to a large extent, the mass loss of the inland ice which they bound. It is therefore important to understand their dynamics, especially how they react to changes influencing their mass and temperature such as may accompany regional climate change. Here, we present an analysis directed towards understanding the dynamics of the George VI Ice Shelf, located on the Antarctic Peninsula. Study of this particular ice shelf is motivated by regional warming in the Antarctic Peninsula observed by Comiso (2000) and Vaughan and others (2001). In previous years, the Larsen Ice Shelf (Rott and others, 2002) and Wordie Ice Shelf (Doake and Vaughan, 1991a, b) have reacted to this warming by disintegrating. The present George VI Ice Shelf may thus represent an ice shelf that is in the process of disintegration (Vaughan and Doake, 1996; Vaughan, 1998), and the characterization of its present flow regime is a primary objective of this study.

The George VI Ice Shelf, situated in the Bellinghausen Sea between Alexander Island and Palmer Land, is an ice shelf of 22 019 km² located at 71–74º S, 66–76º W. This ice shelf is unusual because it has two ice fronts, a northern one from which icebergs calve into Marguerite Bay and a southern ice front in the Ronne Entrance as shown in Figure 1. There are glaciers from both Alexander Island and Palmer Land draining into the ice shelf. The glaciers on Alexander Island are generally small, 9–16 km long and 6–9 km wide, whereas the glaciers flowing from Palmer Land (Fig. 1) have wide entrance gates led by Goodenough Glacier, the main contributor of inflow. The dominant inflow from Palmer Land leads to a flow distribution in which the ice masses show stagnation points: ice moving towards the opposite grounding line, then turning just a few kilometres ahead of the grounding line along Alexander Island (Pearson and Rose, 1983; Reynolds and Hambrey, 1988). Ice-free surfaces of Alexander Island contribute small wind-blown particles of loose sedimentary rock to the surface of the ice shelf, which triggers extensive surface melting in the summer season (Wager, 1972; Pedley and others, 1988; Ridley, 1993). The wintertime snow deposition is therefore not transformed into firn and the ice shelf consists of completely consolidated ice. We will take some advantage of this below.

Another peculiarity of the George VI Ice Shelf is several lakes which connect with the ocean below the ice shelf: socalled tidal lakes (Moutonnée Lake, Ablation Lake, Hobbs Pool) discussed by Reynolds (1981) and Pedley and others (1988). The ice flow forms minor calving margins that can occasionally calve small icebergs into the lakes. These are, however, not included in our numerical simulations.

As shown by Lucchitta and Rosanova (1998), the northern ice front of the George VI Ice Shelf retreated remarkably from 1974 to 1995. MODIS (moderate-resolution imaging spectroradiometer) imagery shows that in early 2002 a rift was created almost parallel to the northern ice front, forming a new iceberg. In another MODIS image, acquired on 19 April 2006, the iceberg has still not calved and thus the northern ice front has not retreated during the past 4 years. The southern ice front, on the other hand, seems to be quite stable in time. The numerous, although small, ice rises and the Eklund Islands act as pinning points, making the position of the southern ice front more stable. A possible disintegration of the ice shelf by an ice-front retreat mechanism is thus expected to occur from the northern ice front towards the southern ice front.

The Model

A three-dimensional, thermodynamical consistent dynamic model for ice-shelf flow was developed by Weis (2001), using a general, continuum mechanical approach. Here, we solve diagnostically the elliptic boundary-value problem for the horizontal velocity in the shallow-shelf approximation (SSA). The deformation of ice is described by Glen’s flow law (e.g. Paterson, 1994):

(1)

where D = sym grad v is the strain-rate tensor (symmetric part of the gradient of the velocity v = (v _x, v_y , v _z )), t^D the Cauchy stress deviator, σ = (tr(t^D)²/2)^1/2 the effective stress, n the stress exponent, T the absolute temperature, A(T) the flow-rate factor and E the flow enhancement factor.

Fig. 1. Map of the George VI Ice Shelf with an arrow indicating the inflow gate. Background image: Mosaic of Antarctica (Scambos and others, 2007).

Insertion of the flow law (1) into the SSA limit of the horizontal force balance yields the elliptic differential equations for the horizontal velocity (v _x, V_y ) (see Weis and others, 1999):

(2)

where the coordinates x and y span the horizontal plane, is the depth-integrated viscosity, ϱ = (ρ_sw — ρ_ice )/ρs _w is the relative density (ρ_ice , the density of meteoric ice, is 910 kg m⁻³ and ρ_sw , the density of sea water, is 1028 kg m⁻³), H is the ice thickness and g the acceleration due to gravity (9.81 ms⁻²). The depth-integrated viscosity is given by

(3)

where z is positive upwards, E_s = E ^1/n is the stress enhancement factor, B(T) = (A(T))^−1/n the associated rate factor, h_s and h_b are the positions of the free surface and the base (ice–ocean interface) respectively, ice thickness H = h_s — h_b , and

(4)

is the effective strain rate (second invariant of the strain-rate tensor).

The elliptic system of differential equations (2) is subject to two different types of boundary conditions: (i) inflow of ice along the grounding line from the adjacent inland ice, and (ii) a vertically integrated stress boundary condition at the front edge. Further, the ice-thickness distribution H(x, y) must be described.

The equations are solved using the commercial software COMSOL (www.comsol.com) using the finite-element (FE) technique. COMSOL isa high-performance FE solver for stationary and non-stationary non-linear systems with an optimized user interface.

Simulation Set-up: Input Quantities

Since there is no information about the vertical ice temperature profile (either in the form of direct information as a temperature profile or indirect information as a function of bottom melting or properties of the ice core), the flow rate factor B (see Equation (3)) was assumed to be constant. The product of stress enhancement factor and flow rate factor was set to BE_s = 0.95 x 10⁸Pa s, a value close to that shown to lead to good results in other ice shelves, such as the Brunt Ice Shelf and Fimbulisen (e.g. Bohlander and others, 2004; Humbert, 2006, 2007; Scambos and others, 2007). Since the George VI Ice Shelf has no firn layer, the flow rate factor was increased slightly to account for an intrinsic firn correction that is included in the value from the Brunt Ice Shelf and Fimbulisen. Note that the ice temperature profile is influenced strongly by the processes at the ice–ocean interface. The magnitude of the flow rate factor is dependent upon whether an ice shelf exhibits basal melting or basal freezing. We assumed here that the George VI Ice Shelf is similar to the Brunt Ice Shelf and Fimbulisen and that it is dominated by basal melting. The FE mesh was triangulated with a maximum triangle size of 4000 m². The number of elements was 107 856.

Fig. 2. Ice thickness distribution.

Boundary: grounding line and ice front

The margins of the George VI Ice Shelf were manually detected from a MODIS Mosaic of Antarctica (MOA) image (formed using images acquired between 20 November 2003 and 29 February 2004 (Scambos and others, 2007) in combination with Ice, Cloud and land Elevation Satellite Geoscience Laser Altimeter System (ICESatGLAS) data. The sharp change in ice surface slope at the transition from ice sheet to ice shelf is visible in the ICESat data and thus gives additional information about the position of the grounding line. The domain begins in Marguerite Bay and extends eastwards as far as Spaatz Island, Ronne Entrance, where the ice rise defines the natural end of the George VI Ice Shelf. The southern ice front consists of two parts, intersected from the De Atley Island ice rise. Several ice rises were introduced to the model: the Eklund Islands and ice rises nearby, one east of the mouth of Goodenough Glacier and one close to the mouth of Millet Glacier. In particular, the ice rises along the southern ice front have a strong influence on the flow field.

Ice thickness

In general, the flow of ice shelves is found to be highly sensitive to the ice thickness (Humbert and others, 2005; Humbert and Pritchard, 2006). Numerical simulations therefore require very precise ice thicknesses. Several ice thickness measurements were performed by the British Antarctic Survey (BAS) in the years 1966–86. In the 1984/85 summer season, traverses perpendicular to the regional geology were run in which seismic measurements of the ice thickness were made. These datasets are available through the BEDMAP database (Lythe and others, 2000). Additional thickness data can be obtained by using ice surface elevation data. Since the whole body of ice consists of completely consolidated ice, the density–depth function is therefore a constant and the transformation from ice surface elevation to ice thickness is simple and precise (see Introduction). This has been carried out using ICESat Laser 2A data. Together with the BAS ice-thickness dataset, an ice thickness distribution was developed (Fig. 2). The thinning towards the ice front occurs with a relatively smooth gradient, whereas the thickness gradient from the thickest area to the southern ice front is much greater. We show below how the thickness gradient, which acts as a driving force for the ice flow, is reflected in the ice velocity distribution.

Fig. 3. Simulated horizontal speed distribution assuming that there is no inflow. The dashed line indicates the position of the flow divide.

Inflow velocity

There are very few field measurements of the ice velocity. Pearson and Rose (1983) show some ice velocity vectors on the ice shelf between the grounding line at Goodenough Glacier and the Coal Nunatak (72–72.3° S, 67.2–68.2° W) and another few points on the ice shelf close to Carse Point (70–70.3° S, 69–68.3° W). Rignot (2006) determined interferometrically derived ice velocities from RADARSAT-1 data. An approximate value of this velocity distribution was used as the inflow velocity for the model. The inflow direction was determined from the flow patterns in the MOA image, as far as they could be detected. Along the grounding line of Palmer Island, the inflow at the glacier entrance gates was set to 100 m a⁻¹ for Meiklejohn, Millet, Bertram, Ryder and Armstrong glaciers. The highest inflow velocity was 500 m a⁻¹ at the mouth of Goodenough Glacier and at an inflow gate at 73° S, 68° W (arrow in Fig. 1 ). The inflow between these two was set to ∼100 m a⁻¹ with slightly varying bearings. The speed of the glaciers of Alexander Island was set to 220 m a⁻¹ at Mercury Glacier and 100 m a⁻¹ at Grotto Glacier. These two prove to be the greatest contributors.

Results of the Numerical Simulations

Using the above set-up, the horizontal flow field was calculated. The results are shown in Figures 3 and 4. In order to distinguish the effect of the ice thickness gradient and the inflow, we show in Figure 3 the computed horizontal flow field in the hypothetical case that there is no inflow at the grounding line.

Assuming no inflow along both grounding lines, the flow field shows a strong difference in horizontal velocities along the two ice fronts (Fig. 3). The southern ice front reaches speeds up to 327 m a−¹ and is considerably faster than the northern ice front where the ice moves extremely slowly at a speed of at most 25 m a−¹. The position of the flow divide is shown in Figure 3 as a dashed curve. To the west of it the ice moves towards the northern ice front; to the east the ice moves towards the southern ice front.

Figure 4 shows the horizontal velocity distribution using the inflow discussed above. The speed at the southern ice front reaches 380ma−¹, whereas the northern ice front shows ice velocities only as high as 141 ma−¹. Unfortunately, the dataset by Rignot (2006) does not show speeds directly along the northern ice front. However, speeds below 100 ma∼¹ between Meiklejohn and Bertram Glaciers provide a good match with the simulated speeds. The narrowest part of the ice shelf is a waist between the inflow of Ryder and Armstrong Glaciers. The measured speeds in this area vary from a few tenths of metres per annum to ∼10 m a−¹. The simulated speeds of at most 15 m a−¹ agree very well with the measured speeds. A further consequence of the inflow is that the position of the flow divide has moved towards the west and therefore the northern ice front.

Fig. 4. Simulated horizontal speed distribution. The inset displays the direction of the flow in an area which is dominated by a stagnation point. The dashed line indicates the position of the flow divide. The colour coding is the same as in Figure 3.

The inset of Figure 4 displays an area which shows a strong stagnation point. The ice moves from Goodenough Glacier towards Alexander Island. A strong change in the direction of the ice movement occurs just a few kilometres ahead of the grounding line along Alexander Island. Melt events at the ice surface in summer produce melt lakes that form in the valleys between flowlines. MODIS images acquired in summer (e.g. 11 January 2005) show the pattern of the ice movement as dark, water-filled lines. The pattern found in the MODIS imagery agrees with the general flow field of the simulations. The computed flow field along the grounding line of Alexander Island agrees well with the flowline map presented by Pearson and Rose (1983), derived from oblique aerial photography and satellite imagery.

Pearson and Rose (1983) show a velocity cross-section measured in the region of the ice rise east of Goodenough Glacier. The speeds in this cross-section are about 15 m a−¹. The ice rises usually reduce the ice speeds of the surrounding ice shelf, acting as obstacles. Since the maps of Pearson and Rose (1983) do not show the ice rises at all, a possible explanation for the difference in the speeds between the simulation and the observation might be a change in the dynamic state of the ice shelf in this area (e.g. by formation of a grounded area).

Comparison between the simulations with and without inflow shows that the influx increases the speeds at the northern ice front dramatically more (560%) than the speeds at the southern ice front (16%). The dynamics of the ice along the southern ice front are thus driven by the ice front, whereas the effect of the inflow velocity is stronger along the northern ice front.

The flow between the Eklund Islands and the grounding line is almost decoupled from the flow of the interior of the ice shelf. A decoupled flow leads to a strong velocity and therefore strain gradient, which typically forms crevasses. These crevasses are highly visible in MODIS images and thus verify the pattern of the calculated flow field. A similar behaviour was found at the Ross Ice Shelf (Humbert and Pritchard, 2006).

Between De Atley Island and the Eklund Islands the ice is moving with speeds of up to 270 m a−¹. However, the ice mass between De Atley Island and Spaatz Island is moving slowly (60 m a−¹), demonstrating a generally decoupled behaviour. This pattern is found in simulations with and without inflow along the grounding line; the low value of flow is therefore assumed to be due to a small thickness gradient.

Conclusions

A new ice-thickness distribution (Fig. 2) was derived using seismic data and ice thicknesses calculated assuming hydrostatic equilibrium from ICESat GLAS surface altimetry data. The hypothetical case of no inflow along the grounding line was studied (Fig. 3), demonstrating that the ice dynamics along the northern ice front are dominated by the inflow whereas the movement of the ice along the southern ice front is driven by the thickness gradient. Simulations taking interferometric synthetic aperture radar (InSAR) inflow velocities into account (Fig. 4) result in speeds along the southern ice front of 380 m a−¹. The northern ice front, however, is considerably slower-moving, reaching velocities of only 141 m a−¹. The computed flow field shows stagnation points (inset, Fig. 4) that are also represented in visible satellite imaginary (e.g. MODIS) and field measurements (Pearson and Rose, 1983; Reynolds and Hambrey, 1988).