Royal Society Publishing

Numerical modelling and data assimilation of the Larsen B ice shelf, Antarctic Peninsula

Andreas Vieli, Antony J Payne, Zhijun Du, Andrew Shepherd


In this study, the flow and rheology of pre-collapse Larsen B ice shelf are investigated by using a combination of flow modelling and data assimilation. Observed shelf velocities from satellite interferometry are used to constrain an ice shelf model by using a data assimilation technique based on the control method. In particular, the ice rheology field and the velocities at the inland shelf boundary are simultaneously optimized to get a modelled flow and stress field that is consistent with the observed flow. The application to the Larsen B ice shelf shows that a strong weakening of the ice in the shear zones, mostly along the margins, is necessary to fit the observed shelf flow. This pattern of bands with weak ice is a very robust feature of the inversion, whereas the ice rheology within the main shelf body is found to be not well constrained. This suggests that these weak zones play a major role in the control of the flow of the Larsen B ice shelf and may be the key to understanding the observed pre-collapse thinning and acceleration of Larsen B. Regarding the sensitivity of the stress field to rheology, the consistency of the model with the observed flow seems crucial for any further analysis such as the application of fracture mechanics or perturbation model experiments.


1. Introduction

Ice shelves are floating and therefore any mass loss does not directly contribute to sea-level change. However, they have been found to play an important role for the dynamics of the Antarctic ice sheet, by acting as a coupling element between the ocean and the grounded inland ice (Rignot et al. 2002; Shepherd et al. 2002; Payne et al. 2004). The recent collapse of parts of the Larsen ice shelf in the Antarctic Peninsula led to accelerated ice discharge and thinning of the tributary glaciers (Rott et al. 2002; Scambos et al. 2004; Rignot et al. 2004). This catastrophic ice shelf collapse also illustrates how dynamic ice shelves are, and their high potential for rapid change under a changing climate. Understanding the dynamics of ice shelves and the controls on their flow is therefore crucial to our understanding of the dynamics of the Antarctic ice sheet and, consequently, our ability to predict near future sea-level change.

The ice shelves in the Antarctic Peninsula have undergone a gradual retreat over the last few decades, which is still ongoing (Vaughan & Doake 1996). The Larsen ice shelf (consisting of Larsen A, B and C) is the largest ice shelf in the Antarctic Peninsula and underwent drastic changes within the last decade. In 1995, Larsen A (Rott et al. 1996) and, in 2002, parts of Larsen B collapsed over a time-scale of days to weeks. Furthermore, a thinning of up to 2.4 ma−1 of the ice shelf has been detected by satellite altimetry between 1992 and 2001 (Shepherd et al. 2003). These changes are likely to be related to the strong warming that the Antarctic Peninsula has been experiencing since the 1950s (King 1994; Vaughan et al. 2001). Several studies (Rott et al. 1998; Scambos et al. 2000; MacAyeal et al. 2003; Hulbe et al. 2004) have focused on the mechanism that caused the ice shelf to collapse. They suggest that the increased availability of surface melt-water has a destabilizing effect on the shelf by enhancing crevasse fracture and, therefore, acting as a trigger for the collapse.

Although this mechanism provides a link between regional climate warming and the final ice shelf break-up, it does not explain the pre-collapse response dynamics of the shelf to a warming climate, such as the causes and effects of the gradual retreat and thinning. The observed thinning, for example, suggests that the Larsen ice shelf reacted to the climate warming long before the final collapse and may have made the ice shelf susceptible to crevasse fracture (Shepherd et al. 2003). It has also been speculated that the observed thinning is mainly related to enhanced melting at the base of the ice shelf due ocean warming.

This study aims to investigate the rheology and the controls on the flow of the pre-collapse Larsen B ice shelf, which is done by combining a numerical ice shelf model with observations from remote sensing using data assimilation techniques.

Previous modelling studies focused on the stress field of the Larsen B to investigate its instability due to crevasse fracture (Doake et al. 1998; Scambos et al. 2000). In these studies, the amount of observational data, used to tune the model, was limited and consisted mainly of a series of velocity measurements in one cross and one longitudinal profile (Rack et al. 2000). To match the observed flow the authors softened the ice rheology along the margins considerably. How much the somewhat arbitrary pattern of softening and the limited amount of observational data affected their results is difficult to judge.

An alternative approach to modelling the flow of the Larsen ice shelf has been undertaken by Sandhaeger (2003), who used more physical criteria by softening the ice in regions where strain rates exceeded a critical value. However, his approach was not focused on matching the modelled flow to observations.

Here, we use a more objective method to adjust our ice flow model to the observational data, which is based on data assimilation techniques. Our observational data consist of two sets of satellite interferometry velocities (InSAR-velocities; from Du, Shepherd, Payne & Vieli in preparation) and the data coverage is therefore much more extensive than that of previous studies.

The control method is used here to assimilate the data in our model. Previously, applications of the control method in glaciology have mainly been to invert for the basal properties of ice streams (MacAyeal et al. 1995; Vieli & Payne 2003; Joughin et al. 2004). Rommelaere & MacAyeal (1997) applied the control method to the Ross ice shelf to investigate its large-scale rheology (i.e. the effective viscosity). In contrast to their application, in this study we simultaneously invert for the ice rheology and the velocities at the shelf's inland boundaries. We also consider the case of Glen's flow law. Furthermore, on Larsen B we deal with much smaller spatial scales and with larger rheological variations. Larour et al. (2004) used the control method in a similar way to this study; however, they focused on the propagation of rifts on a small section of the Ronne ice shelf.

2. Model and inversion method

The inversion procedure used here consists of a numerical forward ice shelf flow model and the control method used to adjust the forward model.

(a) Ice shelf model

The ice shelf model calculates the two-dimensional horizontal flow based on the leading-order, vertically averaged stress equilibrium equations (MacAyeal 1989):Embedded Image(2.1)Embedded Image(2.2)where u and v are the vertically averaged velocities in the two horizontal directions x and y, and h and s are the thickness and surface elevation of the shelf. The vertically averaged shelf density is given by Embedded Image and g is the acceleration due to gravity. The main assumption for the above equations is that the horizontal flow is independent of the vertical coordinate. Using Glen's flow law (Glen 1955) as constitutive equation, the effective viscosity ν is given byEmbedded Image(2.3)where τ is the effective stress, B is the vertically averaged rate factor (in Pa a1/n; Embedded Image, with A being the rate factor in a−1 Pan commonly used in the literature) and in general the flow law exponent is n=3.

At the ice shelf/inland boundaries, the two vertically averaged velocity components Embedded Image and Embedded Image are specified (kinematic boundary condition). At the shelf-front/ocean boundaries the depth averaged sea water pressure is specified (dynamic boundary condition).

The model equations (2.1) and (2.2) are solved numerically for the horizontal velocities u and v using the finite difference method and by iterating for the effective viscosity ν. Given the geometry of the ice shelf, the velocities along the shelf boundaries (inland boundaries) Embedded Image and Embedded Image and the rate factor-field B, the discretized model equations represent a linear system of equations for u and v using ν from the previous iteration.

Because, the finite-difference model requires simple rectangular geometry at the boundaries where the dynamic condition is specified, we artificially extend the ice shelf beyond the front into the domain of the ocean, but the thickness is reduced here to 1.16 m. However, at the edge of the ocean domain the dynamic boundary condition is used and is derived from the depth averaged force balance between ice and displaced ocean water. Artificially extending the shelf beyond its actual front is a typical way of dealing with the shelf-front/ocean boundary in finite-difference shelf models which require simple rectangular geometries (MacAyeal et al. 1996).

(b) Control method

Some of the parameters or quantities that are needed as input of the shelf flow model are unknown or not well constrained. However, information from observations of the modelled counterpart can be used to further constrain and adjust these model parameters. The aim here is to fit the modelled flow (u, v) to its observed counterpart (Embedded Image) by adjusting the model parameters. In this study, the adjustable model parameters (which are input parameters of the forward model) are the ice rheology field B and the two components of the velocities at the inland boundaries Embedded Image and Embedded Image. The misfit between observations and the model is measured by a cost function J which is taken as the least mean square errors between the observed velocity components and their modelled counterparts summarized over all grid points jEmbedded Image(2.4)The weighting functions Embedded Image and Embedded Image measure the uncertainty of the observed velocity components Embedded Image and Embedded Image and are here taken as the inverse of the observational errors Embedded Image and Embedded Image of the two velocity components and are assumed to be uniform over the whole domain. The additional penalty function term Embedded Image includes additional constraints to the problem and will be specified later.

Searching for the model parameter set that gives the best fit between model and observations is the same as minimizing the cost function J while ensuring that the model equations are satisfied. This minimization is done here using the control method which is based on the method of Lagrange multipliers and uses the conjugate gradient descending algorithm (Polak–Ribiere method; Press et al. 1989, p. 303). The details of the control method and inversion procedure used here are described in Vieli & Payne (in preparation) and for other glaciological applications of the control method we refer to MacAyeal et al. (1995), Vieli & Payne (2003) and Rommelaere & MacAyeal (1997).

In our finite difference model, we get two sets of model parameters consisting of the ice rheologies Embedded Image at each grid point within the shelf domain and the boundary velocities Embedded Image and Embedded Image along the inland margin of the shelf. An initial guess for these two model parameter sets has to be given to run the forward model.

(i) Inversion for the effective viscosity ν

On real ice shelves, in regions of high shearing rates Glen's flow law with n=3 may no longer be valid, due to very high stresses or partial decoupling by crevassing (rifts). At least the flow law exponent n is expected to be different and most likely higher. This contradicts the assumptions of Glen's flow law used when inverting for B. Alternatively, we can directly optimize for the effective viscosity ν as done by Rommelaere & MacAyeal (1997), which puts all the ice rheology into the effective viscosity and we no longer rely on Glen's flow law. In that case, solving the ice shelf forward model no longer requires any iteration because the discretized equation system is linear in u and v.

(ii) Penalty functions

In most applications, the model parameters are not completely unknown and rough estimations of maximum or minimum constraints are available. This information is accounted for in additional penalty function terms Embedded Image in the cost function (2.4). Here, we use three additional terms given byEmbedded Image(2.5)A short description of these terms is given below (for details see Vieli & Payne in preparation). The first term Embedded Image constrains the velocities at the inland boundaries (Embedded Image) and penalizes the discrepancy to the initial input values weighted by their uncertainty (discrepancy to the power of two).

The temperature within an ice shelf is not expected to be colder than the temperature of the inflowing ice Embedded Image which can be roughly estimated from mean annual air temperatures. Using the relation between temperature and rate factor (Paterson & Budd 1982), this gives an upper limit Embedded Image for the ice rheology which is accounted in the penalty function Embedded Image. It penalizes values of B that are above B=−18 °C (discrepancy to the power of four). There is no lower limit, to allow softer rheologies in shear or fracture zones within the shelf.

The ice rheology B is in general continuous and rather smooth in horizontal space. Thus, as a further constraint for the ice rheology field the penalty function Embedded Image is added which penalizes the ‘noisiness’ of the B-field by measuring the curvature of B in both horizontal directions of space at each grid point (Thacker 1988).

(c) Acceptance of inversions

The observational quantities are not perfectly known and are subject to some uncertainty. Therefore, the model counterparts do not have to perfectly fit the observations, they only have to fit them within the observational uncertainty. In general, we accept an inversion run when a minimum in J is found (no significant improvement with further iterations) and when the least mean square velocity error divided by the number of grid points with observations is clearly below the observational error. In general, after 22 iterations the inversions presented here fulfil these conditions. The optimization technique does not give us a unique solution for the inverted model parameters and does not guarantee that the global rather than a local minimum of J is found. This has to be kept in mind when interpreting the results.

(d) Observations and model input data

(i) Model setting and geometry

For the surface elevation, the RAMP digital-elevation model (Liu et al. 1999) has been used (corrected for the geoid) and the shelf base and thickness (figure 1a) have been derived from the flotation condition. The location of the grounding line has been estimated from a combination of the RAMP grounded ice-mask and direct mapping from the interferometric velocity analysis used in this study. A mask has been created and is read into the model to indicate whether a gridpoint is on the shelf (floating), is grounded, is a grounded inland boundary point or in the open ocean (see figure 1b). In general, a 2 km horizontal grid resolution is used in the model.

Figure 1

Map of Larsen B showing (a) surface elevation and (b) the mask used in the model indicating the ice shelf (dark grey), grounded inland ice (light grey), inland boundaries where velocities are prescribed (black) and open ocean (white). The dashed lines mark the two profiles used in later figures.

(ii) Velocity data

Two datasets of surface velocities derived from satellite interferometry (InSAR) are used in this study (details see Du, Shepherd, Payne & Vieli, in preparation). The first dataset is an InSAR mosaic, using images from three ascending and three descending orbits, and is from the period between October 1995 and February 1996 (referred as 1995/1996 data) and it provides the x- and y-components of the flow (figure 2).

Figure 2

Observed velocity field from the InSAR-mosaic of 1995/1996 derived from interferometry and projected on to the model grid (in m a−1, contour interval 50 m a−1). The black areas indicate shelf regions with no velocity data. The circles indicate the points on the inland boundary.

The second dataset only uses one descending orbit from the 9 November 1999 (referred as 1999 data) and therefore only gives the velocity component in the satellite look direction (figure 3), which is, however, almost parallel to the flow in the main part of the shelf. A rough comparison of the two datasets shows a maximum increase in speed of about 150 ma−1 from 1995/1996 to 1999.

Figure 3

Observed velocity component in satellite look direction from 1999, derived from interferometry (in m a−1, the contour interval is 50 m a−1). The satellite look direction is indicated by the two arrows. The black areas indicate shelf regions with no velocity data.

Both datasets have been corrected for the tides using a tide model (Rignot et al. 2000); however, it is the uncertainty in the tidal amplitude that dominates the error in the velocity calculations. For both datasets and all components, an uncertainty in the velocity of Embedded Image has been estimated. However, a comparison of the velocities between the 1995/1996 dataset and field measurements by GPS in profiles along and across Larsen B around the same time period (Rack et al. 2000) shows differences well below the estimated 40 ma−1 error.

The velocities at the inland boundary around the shelf, that have to be prescribed in the ice shelf model, are here taken from the 1995/1996 dataset. In regions with no data, which occur often close to the grounding line, assumptions had to be made. Five regions with significant inflow (tributaries) have been identified from the InSAR-data and missing data in the grounding line region have been filled by interpolation between the shelf velocities and the velocities in the grounded tributaries. Elsewhere along the inland boundary, there is almost no inflow and the boundary velocity components Embedded Image and Embedded Image have been set to zero. As initial estimate for the optimization we always use the 1995/1996 data. To account for a higher uncertainty an error of Embedded Image has been assumed for the velocities at the inland boundaries. The optimization for the boundary velocity components Embedded Image and Embedded Image is restricted to the five regions of major inflow. Elsewhere, Embedded Image and Embedded Image are kept constant at zero.

3. Results

If not noted otherwise the results below refer to the geometry (front position) and velocity dataset from 1995/1996.

(a) Forward model run

First the forward ice shelf model is run with a uniform rate factor distribution of 0.41×106 Pa a1/3, which is equal to the averaged B value over the shelf from a typical inversion and corresponds to an ice temperature of −10 °C. This temperature is also about the pre 1981 mean annual surface temperature value for the shelf (Vaughan & Doake 1996). The modelled flow is shown in figure 4 and does not agree with the observed flow in either magnitude or pattern, despite the fact that the velocities at the inland boundaries are taken directly from the observations. Even for a lower uniform B, with a value corresponding to −5 °C, the observed flow is generally too slow and the pattern (shear bands) cannot be reproduced.

Figure 4

Modelled velocity magnitude (contours in m a−1, contour interval 50 m a−1) and direction (arrows) before optimization using a uniform B-field equal to the mean B-field from figure 8, corresponding to −10 °C. The circles indicate the points on the inland boundary where the velocities are prescribed.

(b) Optimization for B and boundary velocities

Using the observed velocity fields and a uniform initial rheology Embedded Image Pa a1/3 (corresponding to −15 °C), we optimize for the B-field and the inland boundary velocities (Embedded Image) simultaneously. Note that the optimization for the boundary velocity components Embedded Image and Embedded Image is restricted to the five regions of major inflow; elsewhere, Embedded Image and Embedded Image are kept constant at zero.

(i) Optimization for time slice 1995/1996

The velocity field after optimization, using the velocity dataset from 1995/1996, is shown in figure 5 and agrees well with the observed flow (figures 2 and 6), both in terms of magnitude and pattern. Ninety-three per cent of the velocity points are within the error uncertainty. The few points with larger discrepancies are mostly located near the grounding line where a higher uncertainty of the velocity data is expected.

Figure 5

Modelled velocity magnitude (contours in m a−1, contour interval 50 m a−1) and direction (arrows) after optimization for the rate factor B. The circles indicate the points on the inland boundary.

Figure 6

Optimized modelled velocities (u and v) against observed velocities (Embedded Image and Embedded Image) for the B-optimization case: (a) x-velocity components and (b) y-velocity components. The dashed line indicates the observational error.

The optimized velocities at the inland boundaries are mostly within the prescribed error estimate, but are in general slightly higher than the input values. For an additional model inversion without optimizing for the boundary velocities Embedded Image and Embedded Image the misfit between observed and modelled flow was significantly higher and not within our criteria of acceptance.

Figure 7 shows how the cost function J and its different terms evolves during the optimization process. The main contribution to the total cost function J is from the discrepancy between observed and modelled velocity fields. In this particular example, after about 22 iterations, J does not improve much further.

Figure 7

Optimization for B: evolution of the total cost function J (solid line with dots) and its individual terms (as labelled) with the number of iterations.

The inverted ice rheology field Embedded Image is shown in figure 8. Narrow bands of low B, located mostly along the margins, are obtained, while in the main body of the shelf, B is rather high and around the initial value of Embedded Image. The obtained soft bands in B correspond to the observed zones of high shearing, which can be seen clearly in the two cross profiles in figure 9. The absolute values in these soft zones go below the value for which a physical interpretation in terms of ice temperature can be made.

Figure 8

Optimized ice rheology B-field in 106 Pa a1/3 (contour interval Embedded Image Pa a1/3).

Figure 9

Optimization for B: (a) and (b) modelled (solid line), observed (solid line with dots) and modelled using a uniform B corresponding to −10 °C (dotted line) velocity magnitudes along profiles A and B (indicated in figure 1b). (c) and (d) Embedded Image for the corresponding cases; the thick black line indicates the location of the shelf (floating), with the thinner line at the ends indicating the location of the inland boundary points.

(ii) Optimized stress field

The first and second principal stresses (Embedded Image and Embedded Image) are shown in figure 10a,b. Positive values of both Embedded Image and Embedded Image correspond to extension and negative values of both to compression. Tensile stresses (positive principal stresses) are in general between 30 and 70 kPa and highest compressional and extensional stresses occur along the margins. The exceptions are the very weak zones, where the stresses are rather low (below 30 kPa). These low values in the very weak zones can also be seen in figure 11 showing the effective stress τ, which is a measure for the stress environment.

Figure 10

Principal stress fields after optimization for B in kPa (contour interval 20 kPa). (a) First principal stress Embedded Image and (b) second principal stress Embedded Image. Positive values indicate extension and negative values compression.

Figure 11

Effective stress τ in kPa: (a) after optimization for B; (b) using a uniform Embedded Image Pa a1/3 (no optimization).

In contrast, in the non-optimized case with a uniform Embedded Image Pa a1/3 (figure 11b) these low stresses in the very weak zones cannot be seen and in general the stress field, particularly in the marginal zones, is much smoother compared to the optimized case. However, within the main body of the ice shelf the values of the τ-field are similar (between 50 and 60 kPa).

(iii) Sensitivity to initial ice rheology

To investigate the robustness of the above inversion results, several inversions with different uniform initial ice rheologies Embedded Image have been performed. The inversions for Embedded Image Pa a1/3 are shown in figures 12 and 13 and clearly show that the pattern of weak bands along the shear zones is very robust. In the main ice shelf body, however, the absolute value of the ice rheology B is not very well constrained. Although in some regions the B-values vary significantly, the stresses seem to be very similar (figure 13).

Figure 12

Sensitivity of optimization of B to different initial rheologies Embedded Image Pa a1/3.

Figure 13

Sensitivity of optimization of B to different initial rheologies Embedded Image: (a) and (b) optimized Embedded Image along profiles A and B (indicated in figure 1b) and (c) and (d) the corresponding optimized effective stress τ.

(iv) Optimization for time slice 1999

An additional optimization for B and the boundary velocity has been done using the InSAR dataset of 1999 as velocity observations (figure 14). Note that here the observed flow field is less well constrained because only the velocity magnitude along the look direction is available. This means the direction of flow (vector) is not fully known. This reduction of observational data makes it easier to adjust our model parameters to fit the observations. However, the resulting optimized B-field (figure 14) is very similar to the 1995/1996 case. Again, clear bands of low B are obtained, mostly along the shelf's margins, which correspond to the zones of high shearing. In the main shelf body, the values of B are high and fairly constant around the initial value. The fit between observed and optimized velocities along the satellite look direction (figure 15) is reasonably good. The exceptions are points from a part of the ice shelf at the left margin (marked with S in figure 3), where the model underestimates the observed flow significantly. A closer look at the observational velocities shows here velocities in satellite look direction of around 120 ma−1. However, the margin of grounded ice in this area is almost perpendicular to the look direction, which would suggest look-direction velocities that are close to zero, as obtained from the model results and observed in the 1995/1996 data. This discrepancy between modelled and observed flow indicates a possible offset in the 1999 InSAR velocity data and may require some reanalysis of the velocity data. However, the 1999 inversion confirms the pattern of weak zones found in the 1995/1996 inversion, although the shelf front has significantly retreated over that period.

Figure 14

Optimized B-field using the 1999 velocity dataset (contour interval Embedded Image Pa a1/3). Note the retreated shelf front position compared to 1995/1996.

Figure 15

Optimized against observed satellite look direction velocities from 1999 for the B-optimization case.

A closer look at the optimized boundary velocities shows that the mean boundary velocity magnitude (from the main tributaries) is about 40 ma−1 higher. This indicates that the flow in the tributaries also increased over that time period, a result which is confirmed by observations from satellites (Rignot et al. 2004).

(c) Inversion for effective viscosity ν and Embedded Image

The same model settings and the 1995/1996 dataset are now used to invert for ν, except that a uniform effective viscosity field Embedded Image Pa a is now prescribed at the beginning of the optimization and we simultaneously optimize for ν (rather than B) and the velocity at the inland boundary.

The general pattern of the optimized effective viscosity field, shown in figure 16, is very similar to the inverted B-fields from before. Bands of low ν (below 5×106 Pa a), corresponding to weak ice, are obtained in the shear zones and are mostly along the margins. In the main shelf body and in regions of major inflow high ν values (between 20 and 30×106 Pa a) are found. Figure 17 shows the velocities and ν values across the shelf for the two profiles A and B (see figure 1) and nicely illustrates these zones of very low ν across the flow located in the shear zones where the velocity gradients are highest.

Figure 16

Optimized effective viscosity ν in 106 Pa a (contour interval Embedded Image Pa a).

Figure 17

Optimization for ν: (a) and (b) modelled (solid line) and observed (solid lines with dots) velocity magnitudes along profiles A and B (indicated in figure 1b). (c) and (d) optimized ν in the corresponding profiles; the thick black line indicates the location of the shelf (floating) and the thin ends the location of the inland boundary points.

Additional inversions with various initial effective viscosities Embedded Image between 10 and 31×106 Pa a show that the pattern of bands with low ν in the shear zones is very robust and almost independent on the initial Embedded Image values (figure 18). The absolute values within these weak zones are also similar (see cross profiles in figure 19). In contrast, within the main shelf body (in between the weak zones) the effective viscosity is not very well constrained (figures 18 and 19).

Figure 18

Sensitivity of optimization of ν to different initial effective viscosities Embedded Image Pa a.

Figure 19

Sensitivity of optimization to different initial effective viscosities Embedded Image: (a) and (b) optimized ν along profiles A and B (indicated in figure 1b), and (c) and (d) the corresponding optimized effective stress τ.

The stress field corresponding to the optimized flow in the main shelf body is quite robust to different initial conditions and values for τ between 50 and 100 kPa are obtained (figure 19d). In contrast, in the marginal regions where the major inflow occurs, the resulting stresses seem sensitive to the initial Embedded Image chosen (figure 19c) and the effective stress sometimes rises to 400 kPa; values which are clearly above the critical stress for ice failure.

In the ν-inversions shown, a penalty function is added to the cost function to avoid ν values that are significantly higher than 30×106 Pa a and therewith to avoid unrealistically high stresses. Inversions without such a penalty show very similar results in terms of the general pattern of the ν-field and stress field; however, the maximum values of ν are higher in the marginal regions where high inflow occurs, and the maximum stress values there are even higher than the already high values obtained for the inversion with a constraint for ν.

4. Discussion

(a) Adjusting the model

Using a uniform ice rheology, we are not able to get the observed flow pattern. Using the relation between B and temperature a vertically integrated ice rheology field could be estimated from the annual mean air temperatures and the melting point temperature of ice at the shelf base. However, this B-field would be rather smooth and close to the uniform case. Softening of the shelf-inland margins by some rule is a better way to adjust the model and has previously been done in applications to the Larsen B (Doake et al. 1998; Scambos et al. 2000). However, the amount of softening needed is difficult to guess a priori and our example of Larsen B shows that in several locations the softening clearly does not follow the shelf-inland margin. In the case of almost continuous data coverage, the control method used here is therefore a more objective way to fit our model to the data and results in a continuous flow and stress field that is consistent with its observed flow. The differences in terms of flow and in particular stress fields between the optimized and the non-optimized uniform B-field case, underlines the importance of a proper adjustment of an ice shelf model to observations, before any further analysis such as the application of fracture mechanics.

(i) Inverted ice rheology

Bands of soft rheology ice have been identified by optimizing the flow model against observed velocities from InSAR. They agree with the shear zones which delineate the fast flowing shelf ice from the stagnant floating ice or grounded ice. In general, the weak zones follow the flow trajectories of the shelf, starting at the margins of the main tributaries (figures 5 and 8). This indicates that the pattern of a few narrow fast flowing shelf tributaries strongly affects the main shelf flow by inducing these weak zones and possibly advect them down stream as suggested by modelling work of Sandhaeger (2003). However, advection of weak ice from the margins in the confluence area of the two main ice shelf arms into the middle of the shelf would then be expected, but cannot be seen in our inverted ice rheology field (figure 8). This could be a consequence of the fact that both shelf arms already flow with comparable high speeds and therefore almost no lateral shearing between the two occurs in the confluence area. Thus, the flow is very insensitive to the ice rheology there, and the weak zone, although in reality present, cannot be detected. This interpretation is consistent with the finding from the sensitivity study on the initial ice rheology above, that the rheology within the main body of the shelf is not well constrained by the observed flow. Due to convergence and hence compressional flow, there may also be some ‘healing’ process going on, meaning a transformation from weak, most likely fractured, to harder and more compact ice.

The location of the bands of weak ice also agree with visible features on the InSAR images and zones in which coherence is lost in the interferometric analysis, which occurs typically in zones of very high deformation or fracturing of ice (see zones of missing data within the ice shelf in figure 2).

In the weak zones, the inversion for B shows values far above the upper temperature limit of temperate ice (0 °C) which corresponds to 0.2×106 Pa a1/3. This indicates, that to allow such high shearing, additional weakening of the ice other than strain heating is needed. This could be through crevassing or changes in fabric (Jacka & Budd 1989), which are both a result of high strain rates. Heavily crevassed areas are weakened because they are no longer mechanically coupled over the full ice thickness but only over a fraction of it. Crevassing may also lead to stiffer ice by cooling through pooling of cold air in crevasses in winter (Harrison et al. 1998). These bands and zones of weak rheology can be compared to the shear margins of Antarctic ice streams for which a strong rheological softening (factors up to 10) is needed to explain the observed flow fields (Echelmeyer et al. 1994; Whillans & Van der Veen 1997). Extensive crevassing can locally lead to decoupling and the ice can strictly speaking no longer be treated as a continuum. Regarding the grid resolution of 2 km, we are not able to resolve these discontinuities, but include them by softening the rheology averaged over one grid cell.

(ii) Inverted rheology pattern

The pattern of weak zones is a very robust result of our inversions. It seems independent of the initial choice of B, and is found for both time slices 1995/1996 and 1999. Furthermore, inverting for ν instead of B produces a very similar pattern of weak zones. The robustness of the pattern suggests that the ice shelf flow is dominantly controlled by these weak margins and the rheology within the main body is secondary. This has two main implications on the dynamics of the Larsen B ice shelf. First, the flow of an ice shelf with such weak lateral zones is more sensitive to changes in its tributary flow, because the main ice shelf body can almost be pushed forward as a block. A corollary to this is that a perturbation in flow at the ice shelf front (due to a retreat) can propagate easily upshelf to the grounded tributaries and therefore affect the tributary flow. This may have consequences for the interpretation of the tributary-glacier acceleration observed by DeAngelis & Skvarca (2003) and Rignot et al. (2004). Second, changes in rheology in these weak zones are expected to have a much higher impact on the flow than the same rheological change within the main shelf body. The weak zones are also the zones of high shearing and therefore strain heating may take place, but also extensive crevassing is expected to occur there. Surface melt-water that gathers in crevasses is likely to deepen them further, but also warms up the ice around the crevasses, which both lead to a further rheological softening of the ice. In these zones of high shear, some crevasses may open up through the whole depth and form rifts. A study by Larour et al. (2004) indicates that melting of ice melange trapped in these rifts or crevasses would significantly increase the propagation rates of the rifts. This implies that these weak zones are also more sensitive to enhanced surface melt. Overall this suggests that the weak shear zones may play an important role in the observed thinning and acceleration taking place on Larsen B. Also, these weak zones may be important for the dynamics of many smaller ice shelves with such existing shear zones, but not for larger ones such as the Ross, the Ronne–Filchner and the Amery ice shelf.

(iii) Stress field

The tensile stresses (positive principal stresses) obtained, for both the B- and the ν-inversions, are over most regions of the ice shelf between 30 and 70 kPa, which is within the tensile values needed for single crevasses to exist (for multiple crevasses larger stresses are needed). This is in contrast to the modelling results of Scambos et al. (2000), which found tensile stresses below 20 kPa over most regions of the shelf and very high values just along the highly softened margins. Our generally higher stresses would therefore imply that active crevassing could have been more widespread than suggested from this previous model study. In this paper, we do not want to explore the details of the discrepancy in the modelled stress fields, but add two comments here. First Scambos et al. (2000) adjusted their model on the basis of only one cross and one longitudinal profile of velocity measurements. Indeed, their modelled flow field is, in some areas away from the profiles, distinctly different from the observed InSAR-field used here (compare figure 2 in this paper to fig. 8 in Scambos et al. (2000)). Second, their model adjustment method, of strongly softening a very narrow zone along the margin, may have affected their model results, as indicated by the concentration of the high stresses in these soft margins (fig. 10 in Scambos et al. (2000)).

In the main part of the shelf, the stress fields obtained in both the B- and the ν-inversions are similar and do not depend on the initial rheological conditions employed. However in marginal regions, close to the high inflow from tributaries, the stresses from the ν-inversion case can be distinctly higher than those of the B-inversion and the final inversion shows some dependency on the initial Embedded Image value. This lack of constraint on the stresses in these regions in ν-inversion is discussed in more detail in the next paragraph and may partially explain the differences between the stress fields generated by the B- and the ν-inversions. However, these features result in a degree of uncertainty that has to be kept in mind in any further analysis of the stress-fields.

The result of very low stress values within the zones of weak rheology occurs for both the B-inversion case and the ν-inversion case, and indicates that this feature is not solely a result of the chosen rheology of Glen's flow law. We interpret these low stresses as follows. Very weak and soft ice cannot resist high stresses and also stress transfer through them is difficult. These very weak zones therefore seem to act to decouple the shelf ice, and consequently the stresses cannot be fully supported and have to be transferred elsewhere.

(iv) Inversion for ν

As mentioned before, Glen's flow law may no longer be valid in zones of high shearing so the interpretations of the B field above have to be taken with care. The inversion for ν does not rely on Glen's flow law and puts the constitutive equation into the effective viscosity ν. However, interpreting the absolute values of ν is more difficult, because unlike B, they cannot directly be related to the temperature. Also, the ν values close to the inflowing tributaries are not well constrained. Depending on the initial Embedded Image, in such regions very high ν-values occurred, which lead locally to unrealistically high stress values (figures 18 and 19). Figure 20a shows the inverted effective viscosity ν against the modelled effective stress τ for each grid point (note the logarithmic scale) to give some information on the constitutive relation hidden in ν. For very low effective viscosity values, which correspond to the very weak zones, ν seems to decrease with decreasing τ. In the range of expected general stresses however, there seems to be no clear relationship. The fact that no values of ν are found above 30×106 Pa a is a result of applying the penalty in the cost function for ν-values above.

Figure 20

Effective viscosity ν against effective stress τ both in logarithmic scale for (a) ν-optimization case and (b) the B-optimization case. The black line indicates the relationship between τ and ν using Glen's flow law (n=3) for a uniform B corresponding to −18 °C.

To compare, figure 20b shows the same plot (ν against τ) for the B-optimization run, which uses Glen's flow law as constitutive equation. For low stresses, a similar relationship as for the ν-case is found, but for higher stresses and at the upper limit of ν we can make out some inverse relationship between τ and ν, which is simply a result of Glen's flow and having an upper limit for B. For a constant B-field, all the points would be on one line with the same slope (see equation (2.3)) as the line indicated in figure 20b. However, because of the large variations in B, except for the points at the upper limit of B (corresponding to about −18 °C), it is difficult to make out a relationship between τ and ν, despite the fact that they are related through Glen's flow law (equation (2.3)).

In marginal regions with high inflow from tributaries, the inversions with lower initial Embedded Image show generally lower ν-values, and more reasonable maximum stresses are obtained (figure 19). Thus, the somehow arbitrarily chosen upper limit of ν at 30×106 Pa a should maybe be chosen lower. The rather large variations in such regions with high inflow in terms of stresses between different ν-inversions indicate that with no further constraints for ν, the modelled stress field for the Larsen B ice shelf in these regions is not well constrained. This uncertainty has to be taken into account in any further analysis of the stress fields.

Despite the difficulties in the inverted stress fields, the general pattern of weak zones (low ν) is a robust result of the ν-inversion case.

(v) Uncertainties and limitations

As illustrated by the sensitivity studies the presented inversion does not give one unique solution for the ice rheology field, and therefore the inverted values are not necessarily the true values, but rather our best estimates with some uncertainty. It is in the nature of the flow of ice shelves that variations in rheology in some parts of the shelf do not affect the overall flow field and therefore the rheology in these regions is not well constrained from the observed flow. Examples of such regions, that are identified by the sensitivity inversion runs, are the stagnant shelf areas between the margins and weak zones, or the main shelf body between the weak zones (figure 12). However, the pattern of weak bands is a very robust result of the inversion.

It is the velocity data that constrain our model, thus the uncertainty of our inversion result is strongly affected by the quality, but also the coverage of the observational data. There are a few areas with no data, mostly close to the grounding line and consequently the inverted rheology there is less certain. However, these regions are small (a few grid points wide) and we do not expect them to dramatically affect our inversion result.

As an example of how errors or poor quality in the data can affect the inversion, we refer to the optimization using the 1999 data (figure 14). The appearance of a second soft band on the left top side of the shelf in the 1999 inversion is expected to be a result of the suspected offset in the velocity data.

Another limitation of our model is the spatial grid resolution. Here, we use a 2 km resolution, thus small features such as ice rises are at the edge of being resolved in the model. For the application to Larsen B, probably more crucial are the tributary glaciers feeding the ice shelf, which are sometimes only a few kilometres wide. Most important for the flow of the ice shelf are the fast tributaries, which are in general also the widest and are captured in our model by at least two or more grid points. Inversion runs with a larger grid size of 4 km showed very similar results and therefore suggest that our results are independent of the chosen grid resolution.

5. Concluding summary

In this study, we present a working algorithm to assimilate InSAR velocities to a shelf flow model by simultaneously optimizing for the ice rheology and the velocities at the inland boundaries. The application to the Larsen B ice shelf shows that a strong softening of the ice in the shear zones, mostly along the margins, is necessary to fit the observed shelf flow. This pattern of bands with weak ice is a very robust feature of the inversion, whereas the ice rheology within the main shelf body is found to be not well constrained. This leads to the conclusion that these weak zones play a major role in the control of the flow of the Larsen B ice shelf and may be the key to understanding the observed pre-collapse thinning and acceleration of Larsen B, and the future behaviour of the whole Larsen ice shelf. This idea is further supported by the fact that the weak shear zones are expected to be more vulnerable to increased surface melting under a warming climate. The existence of these shear margins suggests that stress changes generated by changes in the position of the ice front are transmitted very effectively to the grounding line and to the ice shelf's tributary glaciers.

The inverted rheology field is not unique and also depends on the coverage and quality of the data. However, in our application the rheology pattern obtained is very robust. The inversion procedure presented here is a more objective method than ad-hoc softening of the shelf's shear zones. And importantly, the assimilation procedure results in flow and stress fields that are firstly, continuous, and secondly consistent with the observed flow from InSAR. The stress field of Larsen B is found to be sensitive to the imposed rheology and its pattern. Consistency with the observed flow is therefore crucial for any further analysis, such as the application of fracture mechanics or perturbation model experiments.


We thank the two anonymous reviewers for their useful comments. This research was funded by the Centre for Polar Observation and Modelling (CPOM) a Centre of Excellence of the UK Natural Environmental Research Council (CPOM core project 2.2: Observation and modelling of the Larsen ice shelf instability).


  • One contribution of 14 to a Discussion Meeting Issue ‘Evolution of the Antarctic Ice Sheet: new understanding and challenges’.


View Abstract