Rheology of the Ronne Ice Shelf, Antarctica, inferred from satellite radar interferometry data using an inverse control method

The Antarctic Ice Sheet is surrounded by large floating ice shelves that spread under their own weight into the ocean. Ice shelf rigidity depends on ice temperature and fabrics, and is influenced by ice flow and the delicate balance between bottom and surface accumulation. Here, we use an inverse control method to infer the rigidity of the Ronne Ice Shelf that best matches observations of ice velocity from satellite radar interferometry. Ice rigidity, or flow law parameter B, is shown to vary between 300 and 900 kPa a1/3. Ice is softer along the side margins due to frictional heating, and harder along the outflow of large glaciers, which advect cold continental ice. Melting at the bottom surface of the ice shelf increases its rigidity, while freezing decreases it. Accurate numerical modelling of ice shelf flow must account for this spatial variability in mechanical characteristics.


Introduction
[2] Antarctica is surrounded by ice shelves. Ice shelf flow is governed by many factors, which include: 1) ice spreading viscously under its own weight; 2) ice shelves drag along the margins of ice rises, and they drag over the seabed beneath ice rumples; 3) the absence of basal stress or friction at the bottom; 4) snow accumulation at the surface and ice melting and accretion at the bottom; 5) sea-water pressure along the ice front; and 6) the advection of cold continental ice from the interior. The shape and flow of an ice shelf results from the delicate balance between these factors. The rate of flow follows Glen's law, s 0 = B _ e À 2 3 _ 1 3 , where s 0 is the deviatoric stress, _ e the effective strain rate, _ the strain rate, B the flow law parameter and n the flow law constant.
[3] B depends on ice temperature and ice fabrics [Paterson, 1994] which are not well known. The level of uncertainty on B is high, as illustrated by Paterson [1994], where B varies by a factor 3 at T = 263 K on different ice shelves. For simplicity, ice shelf flow models typically use a uniform ice rigidity, which gives reasonable, first-order fits of ice velocity [MacAyeal et al., 1998].
[4] Here, we propose an inverse control method to infer the distribution of B, over a major ice shelf in Antarctica, the Ronne Ice Shelf. Constraints to the model are provided by the ice shelf geometry, precise and comprehensive ice velocity observations from satellite radar interferometry, and ice thickness and elevation data from BEDMAP [Lythe et al., 2001]. We discuss the impact of the results on ice shelf modelling studies.

Data
[5] Ice velocity data from Radarsat-1 interferometry [Joughin and Padman, 2003] are used to constrain the inverse control model. We limit our study to the Ronne Ice Shelf, setting an arbitrary boundary with the Filchner Ice Shelf (see Figure 1). The model boundaries are the Ronne ice front, the grounding lines determined from ERS-1/2 radar interferometry [Rignot and Thomas, 2002], grounding lines inferred from the Radarsat-1 Mosaic [Joughin and Padman, 2003] and the limited availability of interferometry data. No attempt was made to bridge areas of missing data, because our inverse control model is sensitive to discontinuities in ice velocity. Boundary conditions are: dynamic at the ice front and kinematic everywhere else.
[6] Ice shelf elevation and thickness data are provided by BEDMAP. From the bed elevation, surface elevation and ice thickness, and using 917 kg m À3 for the ice density, we calculate [Jenkins and Doake, 1991] an equivalent amount of air in the firn column of 14 m. Following MacAyeal et al.
[1998], we subtract 14 m from the ice thickness to account for the lower density of firn. In this manner, the ice shelf is treated as having a depth-averaged density of 917 kg m À3 . Ice thickness is limited to 200 m near the ice front. The uncertainty is 70 m in thickness.
[7] We compared this thickness map with that deduced from surface elevation by Bamber and Bindschadler [1997], using the same ice density, but a firn layer of 17.5 m. Model calculations using either maps produce results within 5% of each other. In the remainder of this study, we use BEDMAP thickness.

Model
[8] MacAyeal [1993] developed an inverse control method to determine ice viscosity when the velocity is constrained at the n boundaries. This approach was used with ice shelves where ice velocity is imposed at the ice front [Rommelaere and MacAyeal, 1997].
[9] Problems with this approach include: 1) The boundary conditions at an ice shelf front are dynamic in nature: the stress is constrained by seawater pressure; and 2) The viscosity depends on the effective strain rate _ e and the flow law parameter B. By mixing the two, the interpretation of the control method results becomes difficult.
[10] Here, a new inverse control method for ice shelves is proposed, with natural boundary conditions at the ice front and where the flow law parameter, B, is inferred instead of the viscosity. We present the details of the model in Appendix A. The inverse method finds the distribution B that minimizes the misfit J between InSAR and forward model velocities. Viscosity and flow law parameter are linked by MacAyeal [1989]: where (u,v) is the velocity vector in the x and y horizontal plane.
[11] Model validation follows Rommelaere and MacAyeal [1997]: we use a square ice shelf with a uniform B, except for a square area in the middle of the ice shelf where B is lowered by 50%. The corresponding velocity is calculated using the ice flow model. We then use the result as an observation, starting from a different initial distribution of B, and are able to invert B to within 5%.
[12] We use a mesh of 18,000 finite elements, homogeneously distributed across the Ronne Ice Shelf. The initial value B = B 0 = 2.1 Â 10 8 Pa s 1/3 = 665 kPa a 1/3 is chosen from MacAyeal et al. [1998]. The computation is stopped when J reaches a minimum and variations in J fall below 1%. Computations conducted using different initial values B 0 lead to similar inversion results.

Results
[13] Our results show large spatial variations in B (Figure 1). The minimum value is 300 kPa a 1/3 (55% of the initial value) along the shear margins of the ice shelf, especially along d'Orville Coast, Berkner Island, Evans and Rutford Ice Streams. This value corresponds to an ice shelf with a depth-averaged temperature of À6.7°C (Table 1). The maximum value is 900 kPa a 1/3 (135% of the initial value) at the northern extremities of Korff and Henry ice rises (near the areas marked B and C in Figure 1). This value corresponds to an equivalent temperature of À31°C. B remains close to B 0 south of Korff and Henry ice rises. The ice shelf is more rigid in the wake of Evans Ice Stream towards the ice front, between 78°S and 76°S latitude (areas marked G H and D) and in the wake of stagnant ice between Korff and Henry ice rises (areas A and F). Areas of soft ice appear in the northwestern side of Korff and northeastern side of Henry ice rises and in the middle of the ice front region (near areas B and C and at area E).
[14] The misfit between InSAR velocity, V obs , (Figure 2a) and model velocity, V, (Figure 2b) is shown in Figure 2c. We also show the misfit between InSAR and the initial model velocity V o (Figure 2d) calculated with B 0 . The improvement in velocity is significant: the average misfit is reduced by a factor 4 from ±200 m a À1 to ±50 m a À1 . The flow on the western part of the ice shelf, downstream of Evans Ice Stream and along d'Orville Coast, is reproduced significantly better. The initial misfit in this area is 300 m a À1 . Strong softening is inferred along d'Orville Coast with the inverse method, which increases ice velocity to better fit the InSAR observations.
[15] To assess the sensitivity of B to ice thickness, we carry out another computation with a test ice shelf similar in shape to the Ronne Ice Shelf. The computation indicates that an increase of 1 m in ice thickness leads to a propor-   Figure 2. a: InSAR velocity V obs (m a À1 ) from Joughin and Padman [2003]. b: Forward model velocity V (m a À1 ) after control method. c: Misfit V À V obs between InSAR velocity and model velocity using a control method. d: Misfit V 0 À V obs without control method (B = B 0 ). tional increase of B by 1 kPa a 1/3 . This result suggests that the 70 m uncertainty in ice thickness translates into a 70 kPa a 1/3 uncertainty in B. This is one order of magnitude less than our overall variability in B across the ice shelf (about 700 kPa a 1/3 ), which increases our level of confidence in the results.

Discussion
[16] The presence of strong bipolar singularities of B at the extremities of Korff and Henry ice rise (points B and C) suggests that the stagnant ice lying between them (area A) has a strong influence on ice flow over a large area (areas A, B, C, F and E). This stagnant ice is inferred to be very rigid by the inverse control method. Surface temperature [Giovinetto et al., 1990;Comiso, 1999] reaches a minimum, which could account for the increased rigidity. Sandhäger et al. [2004] suggest that ice in this location is 90% meteoric ice. The layers of slush associated with accretion of marine ice might also insulate the ice shelf from warm ocean water, further decreasing the bulk temperature of the ice shelf.
[17] The singularities at the extremity of these ice rises are located in areas of intense shearing between stagnant ice (area A) and fast ice from Rutford Ice Stream and Carlson Inlet (singularity B) and Foundation Ice Stream (singularity C) (see Figure 2a). This shearing creates an extensive crevasse field in the wake of the ice rise that is visible in Radarsat imagery (Figure 1). Loss of cohesiveness in the field of rifts can explain the low rigidity of ice immediately north of B and east of C. The same type of inconsistency appears along Hemmen Ice Rise, where a similar field of rifts is present.
[18] Residual misfit remains along the shear margins of d'Orville Coast and Berkner Island. This problem arises from the coarseness of the finite element mesh used in our calculation. The grounding zone of the ice shelf is about 10 km wide (Rignot et al. [2000]) whereas the typical size of our elements is 5 km. This limitation is imposed by the computational load of the inversion, but could be reduced by using a finer mesh.
[19] The lack of strong variations in B south of Korff and Henry ice rises can be explained by the nature of the inverse control method. The inverse control method is constrained at the grounding lines. Misfit J is nil at the grounded boundaries of the domain where the forward model velocity is constrained by InSAR. This implies that variations of B at each iteration are small near grounding lines.
[20] Softening along the ice margins averages 55%. On Evans Ice Stream and Rutford Ice Stream, softening along the southern shear margins is larger: this is probably due to the sideward rotation of the ice streams, which increases the stress on the southern margins, and releases it on the northern margins. Our average softening coefficient is consistent with Echelmeyer et al. [1994] along Siple Coast ice stream shear margins, and with MacAyeal et al. [1998] in his forward modelling of Hemmen Ice Rise area. Along the shear margins, viscous heating increases ice temperature, which causes softening. The enhancement factor for ice undergoing shear deformation parallel to the basal plane can reach 19 (62% of softening [Paterson, 1994]).
[21] Near grounding lines, bottom melting is significant (Figure 3), as shown by Rignot and Jacobs [2002], which implies a thinning of the ice shelf and a decrease in the depth-averaged temperature: B should therefore increase. Furthermore, advection of cold ice from the ice streams increases ice rigidity. We observe this effect in the wake of Carlson Inlet, Rutford and Foundation Ice Streams (areas G, H and I). Rigidity downstream of Rutford Ice Stream reaches 750 kPa a 1/3 , equivalent to a depth-averaged temperature of À28°C (Table 1). This analysis is consistent with the conclusions of Huybrechts and De Wolde [1999] that the mean ice shelf temperature is primarily controlled by advection of cold continental ice.
[22] In contrast, an area with intense bottom basal freezing should translate into lower values of B: this is observed around the northeast extremity of Henry Ice Rise, where intense bottom freezing (Figure 3) translates into low rigidity values (area C).

Conclusions
[23] The flow law parameter B for the entire Ronne Ice Shelf inferred from a control method applied to satellite observations of ice velocity is far from being homogeneous. Variations in B range from 300 kPa a 1/3 to 900 kPa a 1/3 . It is difficult to replicate velocity observations using a constant ice rigidity. Rigidity is indeed not only affected by air temperature, but by the full thermal regime of the ice shelf and change in ice fabrics. In effect, ice rigidity is much larger than could be estimated from surface temperature alone: the difference between surface temperature of the Ronne Ice Shelf and its effective temperature inferred from the values of B is of the order of 10 K in the wake of ice streams. This result illustrates the importance of a careful selection of B in numerical simulations of ice shelf flow. The control method used here promises more realistic models of ice shelf flow, and in turn more realistic models of ice sheet flow evolution.

Appendix A
[24] We present the modifications to the inverse model of MacAyeal [1993]. The misfit J, augmented with a Lagrange multiplier vector, becomes after integration, and taking into account dynamic boundary conditions: rgH l @z s @x þ m @z s @y dxdy À Z Z G @l @x 2nH 2 @u @x þ @v @y þ @l @y nH @u @y þ @v @x dxdy À Z Z G @m @y 2nH 2 @v @y þ @u @x þ @m @x nH @u @y þ @v @x dxdy þ Z g 2lnH 2 @u @x þ @v @y n x þ lnH 2 @u @y þ @v @x n y dl þ Z g 2mnH 2 @v @y þ @u @x n y þ mnH 2 @u @y þ @v @x n x dl where G is the integration domain, (u, v) the forward model velocity, (u d , v d ) the InSAR velocity, (l, m) the adjoint vector, n the viscosity, H the ice thickness, g the ice front boundary, r the ice density, g the gravitational acceleration and (n x , n y ) the normal vector to the ice front boundary.
[25] The differentiation of J 0 brings additional term dJ 0 g due to the dynamic boundary conditions: 2h 2 @l @x þ @m @y n x þ h @l @y þ @m @x n y @udl þ Z g 2h @l @x þ 2 @m @y n y þ h @l @y þ @m @x n x @vdl [26] This term brings different boundary conditions for the adjoint equations along the ice front: Z g 2h 2 @l @x þ @m @y n x þ h @l @y þ @m @x n y dl ¼ 0 Z g 2h @l @x þ 2 @m @y n y þ h @l @y þ @m @x n x dl ¼ 0 For the grounding line, the boundary conditions for the adjoint equations remain l = 0 and m = 0.
[27] For MacAyeal [1993], optimization relies on the computation of dJ/dn, gradient of the misfit with respect to the viscosity. Using (2), dJ/dB is deduced from dJ/dn. We then carry out the remainder of our optimization as was done by MacAyeal [1993].