## Abstract

The problem of forecasting the future behaviour of the Antarctic ice sheet is considered. We describe a method for optimizing this forecast by combining a model of ice sheet flow with observations. Under certain assumptions, a linearized model of glacial flow can be combined with observations of the thickness change, snow accumulation, and ice-flow, to forecast the Antarctic contribution to sea-level rise. Numerical simulations show that this approach can potentially be used to test whether changes observed in Antarctica are consistent with the natural forcing of a stable ice sheet by snowfall fluctuations. To make predictions under less restrictive assumptions, improvements in models of ice flow are needed. Some of the challenges that this prediction problem poses are highlighted, and potentially useful approaches drawn from numerical weather prediction are discussed.

## 1. Introduction

Better predictions of sea-level rise are needed, in part to plan defences that protect towns and cities from coastal flooding. Engineering projects of this kind take many decades to plan and build, and may be designed to operate for centuries, so it is important to predict how sea-level might change over the coming century. Many processes can drive sea-level change: thermal expansion of ocean waters, storm surges, and ground water extraction are important. Many bodies of ice, such as Greenland, mountain glaciers, and ice caps, may contribute. However, the Antarctic ice sheet is perhaps the most difficult to predict. This paper focuses exclusively on the contribution from Antarctica.

There are many analogies between predicting the behaviour of the Antarctic ice sheet over the coming decades and the daily weather forecasts provided by centres of numerical weather prediction. The two essential components to a successful weather forecast are the observation network and the atmospheric model (Bennett 2002; Kalnay 2002). The atmospheric observation network monitors pressure, temperature, and other key variables across the globe. The model allows the equations that describe the atmosphere to be integrated numerically, forwards in time. It is only by including a model that there is any capacity for prediction: even the simplest extrapolation of observed change is a model, of a kind. In numerical weather forecasting, the procedures that combine observations and model to provide forecasts have come to be known as data assimilation: the idea being that the skill of the predictive model is strengthened by ingesting data from the observation network. Many of the methods of data assimilation have been adapted to oceanography (Wunsch 1996; Bennett 2002). Here, we describe work towards developing and testing similar methods for glaciology. More specifically, our aim is towards forecasting variations in the size of the Antarctic ice sheet.

More and more observations are being made of Antarctica. Satellites have so improved our knowledge of changes occurring on this continent that, for the first time, an operational observation network of the shape and flow of Antarctic ice is effectively in place. Numerical models of this flow have been used to make predictions of future behaviour (Huybrechts & De Wolde 1999; Church *et al*. 2001), but not in a way that optimally incorporates the many observations that are being recorded. It seems a good time to ask whether the methods of data assimilation can predict accurately how much ice will be delivered to the ocean from Antarctica during the coming century.

In this paper, §2 provides some background to the problem of forecasting the behaviour of the Antarctic ice sheet, and some of the reasons why this has proved challenging. Section 3 describes the close relationship between hypothesis testing, and the combination of data and models through data assimilation. Our intention here is to clarify which problems we can address using the current models of ice sheet flow, and which problems will require further development of those models. Section 4 briefly describes some of the observations that are being made of the Antarctic ice sheet. Sections 5 and 6 describe a linearized ice sheet model, and how to compute the spatio-temporal covariance of its state vector, when a spatially and temporally random pattern of snow accumulation is specified. Section 7 describes a simulation of a forecast of the ice sheet, and examines the information provided to this forecast by the various glaciological observations available. To complete the paper, §8 considers some possibilities for advancing, drawing further analogies with the discipline of numerical weather prediction.

## 2. Background

Climate models predict that warming of the Earth's surface will occur during the twenty-first century, because of increases in atmospheric carbon dioxide concentration (IPCC 2001). This has focussed attention on how such a warming might affect the ice sheets, especially their contribution to global sea-level (Church *et al*. 2001). Ice sheets are sensitive to external changes in climate, but they are also subject to internal instabilities that might allow them to change size of their own accord (MacAyeal 1992). Both types of variability need to be considered.

There is no strong consensus among glaciologists about how quickly the ice sheets can change size, either in response to internally generated variability, or to external forcing by climate (Vaughan & Spouge 2002). Most Antarctic ice is found in East Antarctica, which is widely considered more stable than West Antarctica. Indeed, aside from several glaciers near the coast that are thinning (Wingham *et al*. 2006), the interior ice of East Antarctica does not seem to be changing thickness in a significant way. Radar altimeters have revealed only slight overall thickening during the decade 1992–2003, at a rate that may well be explained by recent snowfall fluctuations (Davis *et al*. 2005).

It has been proposed that the West Antarctic ice sheet might react quickly to climatic warming because it rests below sea-level at its base (Mercer 1978; Oppenheimer 1998). This allows sections to go afloat as *ice-shelves* if they become thin enough. Many ice-shelves are confined by embayments or shoals, which prevent the ice from spreading as rapidly as it would otherwise. Some of the restraining force is transmitted to ice upstream that has yet to go afloat. Removal of the ice-shelf by warming would therefore remove a braking force from these glaciers. Mercer (1978) suggested this could lead to acceleration, thinning and runaway floatation upstream. Alternatively, the grounded glaciers might stabilize after a modest adjustment (Hindmarsh & Le Meur 2001). The conditions needed for stabilization are: (i) additional drag from the glacier bed or lateral margin must restore the balance of forces to the upstream ice; while (ii) the resulting changes in ice-transport and strain-thinning effectively cancel each other out, allowing the thickness to remain steady at the new point of floatation. The relative likelihood of runaway retreat verses stabilization is still being debated (Oppenheimer 1998; Hindmarsh & Le Meur 2001; Van der Veen 2002).

Observations that coastal thinning is actually occurring in parts of West Antarctica, in particular the Pine Island and Thwaites glaciers (Wingham *et al*. 1998; Rignot *et al*. 2002), has added a further impetus to predict the behaviour of this ice sheet. Radar altimetry has shown that a region approximately twice the size of Great Britain thinned by an average of about 10 cm per year over the period 1992–1996 (Wingham *et al*. 1998). At the coast, the thinning rate exceeds a metre per year, with regions of fast flow thinning fastest (Shepherd *et al*. 2001). This is apparently because a recent acceleration of the glaciers has placed them out of balance (Rignot *et al*. 2002; Joughin *et al*. 2003). Thinning at the oceanic margin may be implicated in these changes (Payne *et al*. 2004). In addition, several ice shelves of the Antarctica Peninsula have disintegrated. Glaciers feeding the Larsen B ice shelf have accelerated markedly since its collapse (Scambos *et al*. 2004; Rignot *et al*. 2004*a*). Predicting how these changes will evolve is an important and difficult task.

The climate has natural fluctuations. One approach to determining flood risk would be to assume these are stationary in the statistical sense (so that any particular sequence of events is equally likely to occur at any time). Inferences about what might happen in future could then be drawn from how often similar events have happened in the past. The ice sheets contain a record of past snowfall, from ice cores, and from internal layering revealed by ice-penetrating radar. These might be used to drive a statistical model of ice sheet behaviour. The problem with this approach is that it ignores the unusual warming that is anticipated for the twenty-first century, estimates of which vary widely (IPCC 2001). The predicted regional warming for Antarctica is even more uncertain. Nevertheless, the assumption that changes in the ice sheet will be statistically similar to the past is immediately questionable. Because of this, our forecasts of the ice sheet must ultimately be founded upon a physical understanding of the processes involved, rather than on a purely statistical analysis of the historical record.

The limits of our physical understanding are highlighted by a number of questions that have yet to be answered conclusively, and which have a bearing on the stability of West Antarctica. How does an ice sheet that rests on the ocean floor respond if some part of it thins and goes afloat? Do upstream glaciers always discharge ice more rapidly to the ocean after the disintegration of confined ice-shelves? Does a glacier terminating in deep water inevitably lose mass to icebergs more quickly than a shallow terminus, so that ice can be removed rapidly from deep subglacial trenches? Can warm ocean water penetrate to the ice front via those trenches, removing ice rapidly by melting the underside of ice-shelves? Can ice be lost suddenly from Antarctica through subglacial thermal or hydraulic instabilities of the type that make valley glaciers surge? Can the sediment fail beneath large sections of Antarctica, so that ice that was previously fixed in place at the bed begins to slide rapidly? Any serious forecast of the global impact of West Antarctic ice will have to predict in some detail the consequence for sea-level (or lack of consequence) implicit in each of these questions. The answers are under vigorous pursuit by glaciologists, but we are still not in a position to make such a comprehensive forecast.

Some of our difficulty in forecasting the future of West Antarctica stems from the behaviour that depends disproportionately upon whether some important threshold is crossed. Examples of such thresholds include ice sheet floatation, ice-shelf disintegration, the onset of basal melting, or the yielding of sediments. Mathematical models of the ice sheet that respond linearly to climatic forcing are inappropriate for such systems. Accurate forecasting depends upon predicting when sensitive thresholds are about to be crossed, and how the subsequent response, which may never have been observed, will unfold.

Aside from nonlinear effects, there is uncertainty concerning how the floatation of ice, or the production of icebergs at the margins, should be specified mathematically. This can determine whether the system of equations used to predict the evolution of the ice sheet is stable, unstable, or neutrally stable (i.e. indifferent to small changes; neither returning to previous conditions, nor running away from them). For systems that are unstable, wildly different outcomes might result from slight perturbations to the initial configuration or forcing: the ice sheet might vanish in a runaway retreat, or advance fully to continental margins, chasing nothing but subtle variations of the weather. Even if a stable configuration can eventually be regained, such behaviour could place limits upon the predictability of the ice sheet.

Further complications may arise when poorly measured parameters determine the stability of the ice sheet, or its sensitivity to external climatic forcing. Ideally, we would know, for each of the glaciers that drain Antarctica, the shape and material properties of the underlying rock or sediment, and how the temperature varies through the ice; we would know how much drag is experienced by the glacier from below and from the sides, and how that drag varies with the speed of sliding. At present we can only speculate upon some of these.

## 3. Data assimilation as hypothesis testing

In the shadow of the many difficulties outlined above, can we really claim to have a physical model of the ice that we trust sufficiently to perform a forecast of its future behaviour? Perhaps not yet, but it is worth pointing out that our models are really hypotheses about how the world works, and that there is a very close relationship between data assimilation and hypothesis testing (Bennett 2002). The relationship is so close that even to ask the question of whether the models are good enough forces us to undertake a combined analysis of data and models that is essentially identical to the procedure needed to perform a forecast using those models. Thus, provided that diligence and scepticism are applied to evaluating the success of the forecasts, the business of evaluating the current models, and the business of forecasting on the assumption that those models are accurate, amount to the same thing.

To make a start, we shall examine a hypothesis H: that (i) fluctuations in snowfall statistically similar to the past will be more important in driving changes in the size of the Antarctic ice sheet than twenty-first century warming; (ii) no important thresholds will be crossed during the next century; (iii) no runaway collapse through marginal instability is underway or imminent; and (iv) in these circumstances, a stable, linearized model of the ice flow, along with our best estimates of the important parameters, can provide adequate predictions. This is an optimistic hypothesis. We can test it by seeing whether the thickness changes observed in Antarctica lie within the range of behaviour exhibited by a stable, linearized model of ice sheet flow, forced by natural variability similar to that observed in the ice core record.

Using the linearized ice sheet model, and the available observations, we can forecast the evolution of the ice sheet under the hypothesis H. Since the forcing by weather and climate contains a component that is essentially random, and since the observations contain errors, there will be an envelope of possible outcomes, rather than a single deterministic evolution. While new observations remain within this forecast envelope there is little reason to doubt the hypothesis. If new observations depart significantly from the forecast envelope this will force us to reject the hypothesis, along with some of our optimism.

We should be clear about what can be accomplished using the models that are presently available, and what cannot. We can test whether observed changes in the Antarctic ice sheet are consistent with a slow response to natural variations in snowfall, of the kind that our current models predict. Failure to reject this hypothesis might save us from becoming unduly concerned about changes that lie inside the bounds of such natural variation. Rejection might alert us if we are being complacent about the time-scales on which the ice sheet can respond, or in the ability of our models to make predictions. Rejection in some regions, but not in others, might clarify which of the changes being monitored in Antarctica are most significant, or hint at how the models can be improved.

Using presently available models, we cannot accurately quantify the likelihood of a nonlinear, or unstable collapse of a substantial portion of the ice sheet. We cannot accurately evaluate the risk of coastal flooding posed by Antarctica at levels of significance needed to plan flood defences (e.g. protection against flood risks exceeding one part in a thousand probability). Nor can we accurately predict the financial costs of flooding in the coming century and beyond. To answer these questions, we will need to make forecasts under less restrictive hypotheses than we consider in this paper. Nevertheless, we are beginning to develop the tools that will be needed to do this.

The real value of the simple forecasting scheme described in this paper is perhaps that it can be built upon: even data assimilation schemes for highly nonlinear systems may benefit from having a linearized algorithm available (Kalnay 2002). With this in mind, the method outlined here should be considered as a step towards an operational system for forecasting the contribution of Antarctica to sea level.

The rest of this paper describes a method for predicting the behaviour of the ice sheet under the hypothesis H. The approach is based upon the combination of modelling and observations through data assimilation. The scheme has been developed independently for forecasting ice sheets, and is not a direct application of techniques developed for numerical weather prediction. Nevertheless, it does make use of methods from control theory very similar to those that lie behind numerical weather prediction schemes.

## 4. The glaciological observation network

Nowadays, a wealth of glaciological observations is available from Antarctica. These can be classified into two categories: (i) supporting observations; and (ii) assimilated observations. The distinction is a practical one. In making any individual forecast, supporting observations are treated as though they are free of errors (this assumption might later be relaxed through Monte-Carlo simulations). The assimilated observations are at the focus of our forecasting scheme. They directly inform us about the future state of the ice sheet. To optimize the forecast, we are forced to specify their uncertainties in some detail.

Over recent decades, Earth observation satellites have made valuable contributions to both classes of observation. Extra information has come from airborne surveys, and from the numerous and varied observations recorded by field parties. Laboratory experiments have augmented these by investigating processes related to the deformation and sliding of ice.

The supporting observations include the approximate ice thickness , and bed elevation . These vary as functions of position ** x** within the ice-covered area. These observations are analogous to the orography for numerical atmospheric models, or bathymetry for ocean models. Satellite radar altimetry has provided accurate maps of the elevation of the ice sheet (Bamber & Bindschadler 1997; Liu

*et al*. 1999). Visible and radar images (Ferigno

*et al*. 2000; Jezek

*et al*. 2002) show coastlines, and the location of major features. Maps of ice thickness can be derived from airborne radar, along with field based radar and seismic surveying (Lythe

*et al*. 2001). Floating ice can be identified by its vertical tidal motion, revealed by phase-sensitive radar (Rignot 1998). Supporting observations also include results from laboratory experiments, particularly those that inform on the rheology of the ice, or how it slides. Numerous other measurements are essential in the development of ice sheet models generally.

For numerical weather prediction, the assimilated observations have traditionally been the atmospheric temperature, pressure, and other key meteorological variables, recorded at 6 h intervals by weather stations distributed across the globe. Nowadays, the meteorological observation network includes a wider variety of measurements, many from satellites and aircraft.

To forecast the ice sheet, we need to assimilate the key glaciological observations into our model of ice flow, just as the meteorological observations are assimilated into a model of the atmospheric flow. We shall assimilate observations of , the precise rate of change of ice thickness, , the rate of snow accumulation, and , the rate at which surface material is drawn downwards by the glacial flow. These key glaciological variables are functions of position ** x** and time

*t*.

To begin with, we will ignore changes in the density of upper layers of snow. These can alter the thickness of the ice sheet (Arthern & Wingham 1998), so this effect will need to be considered eventually; but for now, we regard the ice sheet as incompressible. Conservation of mass and continuity are then expressed by the equation(4.1)

The rate of thickness change (*r*) across most of Antarctica has been mapped using radar altimeters that were designed primarily for operation over the ocean (Wingham *et al*. 1998; Davis *et al*. 2005). This time-series of thickness change will be extended further by ICESat, a laser altimeter, and Cryosat, an advanced radar altimeter. These satellites are more suited to operation over undulating ice surfaces. Their near-polar orbits and improved tracking of the ice surface will allow thickness changes to be measured in virtually all regions of Antarctica. Thickness changes have also been recorded by airborne surveys (Krabill *et al*. 2000; Rignot *et al*. 2004*a*,*b*). On broad spatial scales, time varying gravity anomalies, as measured by the GRACE satellite, can also provide information about changes in the mass of ice within a region (Wahr *et al*. 2000).

Our knowledge of *a*, the rate at which snow accumulates on the ice sheet, relies heavily upon field observations, including ice-cores, snowpits, snow-stakes (Vaughan *et al*. 1999; Giovinetto & Zwally 2001). The major difficulty is interpolating between sparse field observations to provide a map of the accumulation rate. This interpolation can be assisted by ground penetrating radar, which reveals the depth to which snow, deposited in previous years, has become buried (e.g. Arcone *et al*. 2004). Satellite observations of thermally emitted microwaves can also contribute useful information (Arthern *et al*. 2006). Numerical modelling of the atmosphere provides useful information about spatio-temporal fluctuations in accumulation rate (Genthon & Krinner 2001; van Lipzig *et al*. 2002; Bromwich *et al*. 2006).

Subsequent to its deposition upon the ice sheet, snow and ice are transported by glacial flow. Because this spreading flow is divergent, material at the surface is drawn downwards. Mathematically, the vertical component of the material velocity at the surface (i.e. the rate of drawdown *D*) is represented by(4.2)i.e. the horizontal divergence of the volumetric flux carried by the glacial flow. The flux is defined as the depth integrated velocity(4.3)

The rate of drawdown, or divergence (*D*) cannot be measured directly over broad areas, so it must be derived from other measurements. There are several methods for estimating the velocity (*v*) of the ice flow. Global positioning systems allow the motions of the stakes to be determined accurately, and stakes may be revisited over many years to provide an average velocity. Alternatively, the motion of natural features, such as crevasses, can be tracked in satellite images. In regions of rapid flow, crevasses can be identified in images taken several years apart, so the average velocity over many years can be obtained. Slower moving ice has fewer features to identify, but even here, the ice motion can be tracked from space using the technique of Interferometric Synthetic Aperture Radar. By analysing the phase and amplitude of echoes collected by an imaging radar system, the velocity of slow moving, featureless sections of the ice sheet has been determined (Joughin & Tulaczyk 2002; Rignot & Thomas 2002; Joughin & Bamber 2005). These methods measure velocity at the surface of the glacier. However, unless the glacier is afloat, or flows over extremely soft sediment, the ice velocity will decrease with depth, especially near the bed where shearing is concentrated. If the ice at the glacier bed is frozen in place, the velocity will decrease to zero there. To infer the ice flux from the surface velocity (), and the ice thickness (*H*), assumptions must be made about the velocity profile through the ice column.

More complete descriptions of the measurements that make up the glaciological observation network may be found elsewhere in this volume. Next, we turn our attention to the glaciological model that forms the second essential component of our data assimilation scheme.

## 5. A linearized ice sheet model

A common approximation of glaciological modelling is the *shallow ice approximation* (Hutter 1983). This treats the ice sheet as a shallow film that flows and spreads under its own weight. Assuming further that all resistance is provided at the glacier bed, the speed of the ice is determined by the local gradient of the surface elevation *s* and the local ice thickness *H*. Under this approximation the ice flux is(5.1)(5.2)(5.3)where is the density of ice and *g* is the acceleration due to gravity. The stress exponent *n*, and the depth-varying, temperature-dependent, coefficient derive from the Glen relationship relating strain rate *n* to deviatoric stress *σ*,It is commonly assumed that some effective value of *A* applies throughout the ice column, so that *c* is constant. Appendix A summarizes a linearized model (Hindmarsh 1997) of how the ice thickness responds to changes in accumulation rate. The model has the form,(5.4)where represents small fluctuations in ice thickness, driven by fluctuations in accumulation rate. In this linearized model, *θ* is a linear, second-order differential operator, providing a damping feedback to the thickness fluctuations that represents the change in flow. For notational convenience, we introduce and to represent fluctuations in the rate of thickness change, and the rate of drawdown, respectively.

## 6. Covariance of the ice sheet state

We can use equation (5.4) to investigate how statistical fluctuations in the rate of thickness change and flow divergence relate to each other, in space and in time, and how each relates to the accumulation rate fluctuations . If we have some independent information about the magnitude, length-scale, and time-scale, of snowfall fluctuations, this can provide prior information about the most likely state of the ice sheet. This prior information is available to us even before we have introduced any observations of *r*, *a* and *D*.

The art of data assimilation lies in weighing the information provided by the new measurements against the prior information already available, updating our estimate of the state of the ice sheet accordingly. Assuming that measurement errors and fluctuations , and have Gaussian distributions, the procedures for doing this are straightforward, as long as we can estimate the covariance of the errors, and the covariance of the fluctuations. Here, we consider the calculation of how , and covary in space and time.

Our hypothesis excludes internally generated instabilities that are associated with floatation, the onset or shutdown of sliding, or thermo-viscous instabilities. In our model, all changes in the form and flow of the ice sheet are driven randomly by the external climate via the fluctuating accumulation rate. (More generally, we can extend our hypothesis to allow for some underlying ‘background’ change that is described as a linear combination of known functions; this is described further in §8.)

The accumulation rate fluctuations cause the thickness of the ice sheet to vary, but these changes are damped by a negative feedback. As the ice thickens, the flow of ice accelerates, which tends to restore the thickness to lower values. The time-scale of this equilibration is set by the viscosity of the ice, and is expected to be several thousand years. This means that changes in flow on short time-scales are relatively small: rapid accumulation rate fluctuations will be highly correlated with rapid fluctuations in the rate of thickness change. On much longer time-scales, the thickness fluctuations are effectively damped out, so that slow variations in accumulation rate should correlate well with slow changes in ice flow.

The question of how the thickness change and accumulation rate changes are correlated under the shallow-ice approximation is dealt with in more detail in an earlier paper (Arthern & Hindmarsh 2003). There, a dynamical model for the thickness fluctuations is introduced, which describes the evolution of the state vector(6.1)when accumulation rate fluctuations are described by a continuous Markov model of the form(6.2)The parameter *T* represents the time-scale of correlations in accumulation rate. The random forcing is a Gaussian white noise process, uncorrelated in time, but correlated spatially over some distance characteristic of the meteorology of the area.

For the purposes of numerical computations, we define grid locations, and represent spatially varying fields as vectors. Discretized forms of equations (5.4) and (6.2) can then be combined as(6.3)with(6.4)where is the discretized equivalent of the operator *θ* and is a vector of spatially covarying white noise processes such that , and the matrix ** Q** is chosen with time-scale

*T*, length-scale

*L*, and variance , of fluctuations in accumulation rate chosen to approximate the meteorology of the local area.

Fluctuating quantities , , and can be written as linear operations on as follows:(6.5)(6.6)so the general covariance with can be written , where .

Following Bryson & Ho (1975, pp. 328–334), and assuming that all fluctuating quantities are statistically stationary in time, can be written(6.7)where is obtained by solving the Lyapunov equation(6.8)where is a matrix (the ‘transition matrix’ in the terminology of Bryson & Ho 1975) defined by(6.9)

## 7. Simulation of an ice sheet forecast

A full description of our estimation algorithm has been presented elsewhere (Arthern & Hindmarsh 2003), in which we describe how the linearized ice sheet model can be used to select the optimal weighting of observations. That paper describes how to obtain an optimal estimate of three important glaciological quantities: the rate of thickness change (*r*), the rate of snow accumulation (*a*), and the horizontal divergence of the depth-integrated volumetric ice flux, or drawdown rate (*D*).

In our algorithm, the continuity equation (4.1) is enforced as a constraint. The fields , are modelled as(7.1)(7.2)i.e. as zero-mean random fluctuations and , each superimposed upon a ‘background’ that is represented well by a linear sum of basis functions and . The background fields are composed of spatio-temporal patterns, multiplied by coefficient or , respectively. For the inversions described in this paper, we used a basis composed of biquartic (bi-directional, fourth-order) polynomials. Other possibilities are discussed in §8. The coefficients and have yet to be determined. To proceed further, we must introduce the assimilated observations.

Measurements () of the rate of thickness change are considered to be temporal averages of *r*, over some interval , corrupted by noise (). This noise is assumed to have zero-mean and known covariance. We write(7.3)Similar expressions denote measurements of accumulation rate, and of flow divergence, and their errors and . Grouping all such observations as vectors, , , and , the observation vector and noise vector are defined as(7.4)

An optimal estimate of the rate of thickness change , is derived as a linear combination of the observations:(7.5)with weights chosen such that is the best linear unbiased estimate for some temporal average () of the rate of thickness change(7.6)There is no guarantee that the measurement locations are coincident with the locations for which estimates are required, so there is a component of spatial interpolation to the estimation problem. Neither is there any guarantee that the measurement interval corresponds with the estimation interval , so there is also a component of temporal interpolation or extrapolation. A complete description of our algorithm for calculating the weights is given elsewhere (Arthern & Hindmarsh 2003). The optimization of the weights requires the covariance of the errors *ν* to be known. Furthermore, we need to know how the fluctuations and covary with themselves, and with each other, in both space and time. This spatio-temporal covariance can be estimated using the linearized model of the ice sheet flow, as described above, provided we have some idea of the length-scale *L*, the time-scale *T*, and the variance , characteristic of the fluctuations in accumulation rate.

Identical methods can be used to estimate temporal averages of accumulation rate (), or flow divergence (), allowing a full analysis of the mass balance of the ice sheet to be performed over any time interval .

Having formulated the estimation problem in this way, we are completely free to choose the time interval that our estimates , and represent. Arthern & Hindmarsh (2003) dealt mainly with estimation of the current state of the ice sheet. Simulations were performed showing that the algorithm could perform better than either of two more traditional methods: (i) direct estimation of volume change (*r*) from radar altimetry; or (ii) ice-budget calculations balancing separate observations of ice accumulation (*a*) and losses (*D*) integrated over each drainage basin.

Here, we aim to show the potential of the method for forecasting the future evolution of the ice sheet. Section 7*a* describes a simulation of such a forecast, and an investigation into how the different types of observation contribute information to the prediction.

### (a) Practicalities

Of all data assimilation methods applied to numerical weather forecasting, our approach corresponds most closely to optimal interpolation in *space* and *time*. The ice sheet model is used to evaluate the covariance matrix that will optimize this interpolation/extrapolation. We are using a vertically integrated ice flow model, so our method is two-dimensional-space, one-dimensional-time optimal interpolation. The vertical integration has important practical consequences. It reduces the number of variables that must be used to describe the state of the ice sheet, compared to vertically layered models.

Our model has approximately 10 000 state variables for a 50 km grid spanning grounded portions of the Antarctic ice sheet. This is small compared with the dimensionality of state-space for a numerical atmospheric model. Nevertheless, it is large enough to pose practical problems for the estimation problem. The size of matrices that we must manipulate to obtain an optimal forecast scales with the square of the number of state variables, i.e. 10^{4} by 10^{4}. We have found that this necessitates solution of (6.8) by specialized iterative methods that take account of the sparsity of ** F**, rather than through direct matrix manipulations. Krylov methods are particularly attractive because these allow approximate solutions to (6.8) to be obtained relatively cheaply, in a way that scales manageably with grid resolution (Li & White 2002). Once we have obtained a Krylov subspace approximation to the operator

**, this can also be used to approximate (e.g. Edwards**

*F**et al*. 1994) the exponentiated operator

**, which is needed to derive the temporal covariance of the state variables. In effect, we assimilate the observations into a reduced order model (Li & White 2002).**

*Φ*### (b) Analysis and forecasting

Some approaches to data assimilation make explicit distinction between an analysis step (i.e. estimation of the initial conditions for the forecast), and a forecast step (i.e. propagation of those initial conditions into the future via the dynamical model). Our method does not explicitly make this distinction. Nevertheless, as a necessary condition, it ought to estimate the present state, and indeed the past state, of the ice sheet reliably.

In this paper, we illustrate the application of the algorithm using simulations: no real observations are considered; no predictions of the actual behaviour of the ice sheet are made. In this section, we consider how accurately present, past and future states might be estimated from a plausible, but hypothetical, distribution of measurements, scattered across Antarctica. The hypothetical measurements are also distributed through time according to the historical exploration of the ice sheet. This includes the era of satellite altimetry and Synthetic Aperture Radar (approx. 1990–2000), field expeditions (approx. 1950–2000), and analysis of ice cores originally deposited thousands of years ago.

We performed a test of our algorithm using synthetically generated observations. A forward model of the ice flow was run for 20 000 years with random pattern of snowfall to produce a synthetic ice sheet. The details of this forward model calculation are described in Arthern & Hindmarsh (2003) for the smaller Berkner Island ice cap. Here, we have performed a model run for grounded portions of the Antarctic ice sheet, excluding the mountainous terrain of the Antarctic Peninsula. Bed topography is from BEDMAP (Lythe *et al*. 2001), long-term average accumulation rates were from Vaughan *et al*. (1999), with a random component added that has correlation length *L*=500 km, *T*=1 year, and , i.e. interannual variability at approximately 25% of the mean annual accumulation. The model grid spacing was 50 km.

From our synthetic ice sheet, we sample the rates of surface elevation change, snow accumulation, and ice flow divergence, just as the real ice sheet is observed by radar altimetry, ice-cores, and satellite measurements of ice flow. Adding noise to the sampled values simulates the effects of measurement errors. For the purposes of this simulation, the errors were assumed to be uncorrelated, with identical standard deviation (0.1 m a^{−1}) for all observations, , and . Table 1 summarizes the synthetic measurements. For convenience, we use calendar dates to refer to model output, although this is an arbitrary choice for such a simulation.

Armed with these synthetic datasets, we ask how well we can reconstruct the state of the synthetic ice sheet, and how the weighting applied to the different datasets depends upon whether we are trying to estimate the past state of the ice sheet, the present, or to forecast the future.

Figure 1 shows the reconstruction of the past (1900–2000), recent (1990–2000) and future (2000–2100) states of the ice sheet. The panels show the average rates of thickness change , accumulation , and drawdown , derived for our synthetic ice sheet. The lower panels show the fields , , and recovered by the assimilation algorithm. It is clear that, at least for these observations, the present-day rate of thickness change can be determined more accurately than the average rate that occurred over the past century, or will occur over the next.

Figure 2 shows estimates derived decade-by-decade using the same synthetic observations. The method described in this paper was used to forecast the evolution of the synthetic ice sheet 100 years into the future, beyond the end of the measurements. A prediction is shown for each decade, along with a 95% confidence interval. This confidence interval was derived using the posterior error covariance for best linear unbiased estimation, assuming a Gaussian distribution for errors and fluctuations (Kitanidis 1997). The confidence interval includes the effects of snowfall fluctuations. Inevitably, these contribute significant uncertainty to the estimation and prediction. Some of this uncertainty is attached to the estimation of the coefficients that multiply the basis functions in equation (7.2). These coefficients capture any underlying background trend in the thickness over time. Uncertainty in their estimation lessens the confidence in estimates towards the start and finish of the interval considered.

If the results shown in figure 2 were a real forecast, observing the unfolding evolution of the ice sheet from 2000 to 2100 would give little reason to doubt the hypothesis *H*. This demonstrates the potential of our algorithm for performing such a test. Of course, we cannot say yet whether the real ice sheet will remain within the forecast envelope, or depart from it.

It is interesting to examine the weighting applied to the various datasets. This is shown in figure 3 The weighting depends upon which target location we choose, and we cannot show weights for every observation and for every potential target location. Instead, we targeted every grid cell within the model domain in turn, multiplied the weights by the area of that cell, summed these products over all cells, and plotted this sum at the location appropriate to the observation. The result of this procedure is that the weights plotted are those needed to estimate the spatially integrated volume change, not the value at any particular location. Only the magnitude of the weights is plotted, not the sign. Observations that are weighted strongly, either positively or negatively, are shown in red. The figure shows how the hypothetical observations would contribute to assessing Antarctica's impact upon sea-level over various time intervals. Weights are shown for three different target intervals (top, 1900–2000; middle, 1990–2000; bottom, 2000–2100).

The optimal weight for any observation depends upon what other data are collected, where the data are collected, and when. The weights also depend upon what target interval is selected, e.g. past, recent, or future. For the estimation of the present day rate of thickness change (middle row), radar altimetry is weighted more strongly than other measurements. As the target interval extends to the past century, before the satellite era (top row), an ice-budget calculation is favoured (i.e. input from snow accumulation minus export from ice flux), so these observations receive stronger weighting. Extending the forecast a century into the future (lower row) also favours the ice-budget approach. This reflects the fact that thickness changes occurring over a single decade do not represent the long-term behaviour particularly well, because of snowfall fluctuations. Despite this, the weighting of the thickness change observations cannot be neglected, even for forecasts 100 years into the future. Under our hypothesis, the altimetry observations from 1990 to 2000 provide useful information for forecasting the behaviour of the ice sheet over the coming century.

Before carrying out these simulations we would have struggled to predict exactly how the various observations would be weighted by our algorithm. With hindsight, we can see that the weighting is qualitatively in accord with current glaciological thinking: that ice-budget calculations represent the longer-term behaviour of the ice sheet, while radar altimetry observations of thickness change are affected more by the year-to-year fluctuations in snowfall.

The exact weighting of the various observations is sensitive to their error covariance, and to the assumptions built into our hypothesis. It is also sensitive to the number and location of observations available. This means that the weighting shown here can only be seen as an example, indicative of how observations might be weighted in a more realistic situation. Nevertheless, our simulations show that all of the various types of observation that we have considered are potentially important in predicting changes in the mass of the Antarctic ice sheet over the coming century. This supports our contention that the best assessment of the mass balance of the ice sheet will come from an optimal combination of observations, and glaciological modelling, through data assimilation.

## 8. Challenges

The simulations described above demonstrate that a practical forecast of the behaviour of the ice sheet is possible, at least under our limited hypothesis. A number of challenges must be met though, if we are to progress towards operational forecasting of Antarctica's contribution to sea-level rise. Some of these challenges are identified below. The strong analogy between our glaciological problem and the problem of forecasting tomorrow's weather, for which data assimilation methods are further advanced, suggests approaches to many of these challenges that may prove profitable.

### (a) Uncertainty in parametrization

It is difficult to formulate mathematical models of complicated physical systems like the Earth's atmosphere, oceans or ice sheets. Often, the aspects of such a model that are most provisional are related to processes that operate on spatial scales finer than the resolution of the numerical model. If the resolution cannot be increased, the effects of these fine-scale processes must be quantified in terms of the model variables available at the coarser resolution. This approach has become known as parametrization. Examples include simple prescriptions of turbulence in atmospheric modelling, sub-grid-scale mixing in oceans, or the sliding of ice sheets over soft sediment.

Often there is a choice to be made between parametrization schemes, but insufficient experimental evidence to make this choice straightforward: the relevant scales, although small in terms of model resolution, are too large to perform full-scale experiments directly. This is certainly true of glacier sliding, where the stress transmitted to the bed varies on a wide variety of scales. Even if a scheme can be derived from small-scale experiments or from general physical principles, the precise values of the numerical parameters may be subject to great uncertainty. This model uncertainty makes any forecast provisional: a different prediction might have been made, had things been parametrized differently.

Some of the effects of model uncertainty can be assessed by varying model parameters, or swapping parametrization schemes entirely, then repeating the forecast. This can reveal how different the forecast might be for plausible alternatives; giving a Monte-Carlo simulation of how model uncertainty propagates into the prediction. The supporting observations discussed earlier can be varied in a similar way to bracket their uncertainty. After many such changes have been made, an ensemble of forecasts will be available. Their spread gives a measure of the uncertainty in the prediction.

Even if multiple forecasts are performed, there is a risk that predictions will be inaccurate, or uncertainties underestimated, simply because we have not employed a broad enough selection of parametrization schemes in our ensemble of forecasts. It is better to err on the side of caution, and include parametrizations that we have some reason to doubt, rather than leave out parametrizations that we are suspicious of, but that might perhaps be realistic. As we compare our predictions with the observed large-scale evolution of the ice sheet, some members of the ensemble of forecasts will no doubt become inconsistent with the observations. These members of the ensemble can be rejected at that stage, playing no further part in the forecast, or in assessing its likely error.

Although there is inevitably a degree of subjective judgement involved, it seems sensible to initially include as many parametrizations of ice rheology, basal sliding, and calving as is practical in our ensemble forecasts of the ice sheet. In this spirit, we plan to perform Monte-Carlo forecasts that include linear and nonlinear ice rheologies, with a wide range of stress dependence, and a wide variety of till rheologies. We will also need to consider a number of alternative boundary conditions representing grounding line migration and calving.

### (b) Limitations of the ice flow model

The isothermal, shallow-ice approximation used for these simulations provides one of the simplest ice flow models imaginable. All resistance to the motion of the ice is from sheer stress between vertically separated layers that are moving at different horizontal velocities. Recent studies (Pattyn 2003; Payne *et al*. 2004) suggest that it would be beneficial to include other stresses. Although it would add complication, a linearized model could incorporate stress from shearing at the side margins of glaciers and ice streams, as well as forces brought about by the longitudinal compression or extension of the ice. Identical methods to those described here could then be used to develop a forecast.

Variations in the flow of the ice sheet can also be produced by changes in the viscosity of the ice. These are most likely to happen because of warming or cooling, although anisotropy of the crystal structure might also be important. To include thermal effects, numerical models evolve the temperature field, coupled to the flow of the ice, via the viscosity. These models tend to exhibit nonlinear instabilities because of feedbacks between the surface-geometry, flow and temperature fields. It may prove difficult to derive a linearized thermomechanically coupled model that captures the most important processes. More likely, it will be necessary to adopt more general methods of data assimilation that are capable of dealing with nonlinear behaviour.

### (c) Nonlinear behaviour

The methods described in this paper are only applicable when the ice sheet responds linearly to climatic forcing. If changes in the ice sheet become large enough, then higher order effects could become important. More extreme types of nonlinear behaviour might be exhibited if the model equations allow multiple solutions for the ice flux. Simplified systems of equations representing ice flow can produce behaviour analogous to the surging of mountain glaciers if multiple solutions for the ice flux can exist (Fowler 1987; Schoof 2004). If ice flow models are to capture the variability of the ice sheets that is generated by internal processes then behaviour of this type must be incorporated. The linearized model, equation (6.3), would then be replaced with a more general nonlinear model(8.1)

Data assimilation techniques for general models of this form have been developed for numerical weather prediction (e.g. Kalnay 2002). One approach is to evolve an ensemble of forecasts using the nonlinear model, and estimate the covariance of the ice sheet state using the sampling statistics of that ensemble. Another method is to linearize the ice sheet model about the estimate of the ice sheet state. Then methods similar to those described in this paper can be employed to advance this estimate over some limited interval of time. This provides a new estimate of the ice sheet state, which the model can be linearized about, and so on. It is not yet clear which approach will prove most suitable for ice sheets.

### (d) Testing other hypotheses

The form of (7.2) is the most flexible description of the ice flow that can be solved by linear methods. Here, we consider some examples of the basis functions that might be useful in practice. A simple choice is *N*_{i}=*N*_{j}=1, with *ϕ*_{ri}(*x*, *t*)=1 and *ϕ*_{aj}(*x*, *t*)=1 everywhere. This corresponds to using the method of Ordinary Kriging (Kitanidis 1997) to interpolate the observations in space and in time, albeit under the extra constraint imposed by the continuity equation, and with the extra complication that a linear ice flow model is used to estimate the multi-variate covariance of the fluctuations *r*_{f} , *a*_{f} and *D*_{f} in space and time. This simple parametrization would correspond to a ‘background’ accumulation rate that is uniform in space and time, with the random component, specified according to the Markov model, superimposed upon it. The ice flow ‘model’ embodied by this parametrization is equally simple: that there is some uniform ‘background’ rate of thickness change, with the statistical fluctuations driven by the random snowfall pattern superimposed upon that.

It seems unlikely that the background accumulation rate, and rate of thickness change, would be uniform across the ice sheet. To allow spatial variations in these rates, we could include higher order polynomials to capture spatial trends, on whatever scale is thought reasonable. This approach was used in the simulations described in §7, where biquartic (bi-directional, fourth-order) polynomials were used. If temporal trends are also anticipated, terms that are linear, or higher-order, in time could also be included. These examples illustrate that before the observations are introduced, the underlying ice flow ‘model’ is really just a hypothesis of our choosing, tentatively proposed, and having whatever prior assumptions we feel are reasonable built in. By judging which of the coefficients and are significantly different from zero, we may obtain some idea which of the spatio-temporal patterns and are actually needed to explain the observations.

We might have independent reasons to anticipate thinning in certain parts of the ice sheet, perhaps from modelling the response of the ice to the warming that occurred at the transition from glacial to interglacial conditions. Then an effective strategy would include the predicted thinning rate as one of the , and the modelled accumulation that led to it as one of the . Performing the inversion, and judging whether the corresponding coefficients and differ significantly from zero, would provide a test of whether these patterns can actually be detected within the available observations. This would also allow the ongoing effects of deglaciation to be included in any forecast.

Any pattern of thinning or thickening could be included in a similar way. A further example would be to include the snow accumulation rates predicted by simulations of twentieth to twenty-first century climate from general circulation models as one of the , and their estimated contribution to the elevation rate as one of the . Any number of patterns could be included at once, to perform a ‘detection-attribution’ experiment for ice sheet mass balance. A similar approach has been used to investigate the various forcings of global temperature (IPCC, 2001). The close link between data assimilation and hypothesis testing means that even our simple forecasting scheme provides useful machinery for testing whether or not such patterns can be detected within the available observations.

### (e) Including further observations

It will be important to include other observations that have a bearing on the prediction problem. Measurements of the changing gravity field could be assimilated in a straightforward way. This may help to avoid confusing changes in density with changes in mass (Wahr *et al*. 2000). Information about past accumulation rates can be derived from layering within the ice sheet, as observed by airborne radar surveys. Assimilating the layering directly would constrain how ice has flowed in the past. Airborne imagery has been used to chart the evolution of the coastline over time (Cook *et al*. 2005). The grounding line can also be mapped sequentially (Rignot 1998). Assimilating these observations would mitigate some of the problems of parametrizing grounding line migration and calving. The simulations considered here provide a good indication of the relative importance of the various types of measurement in our hypothetical observation network. Similar computations using more realistic error assessments, and a more realistic spatio-temporal distribution of data, would allow targeted observations to be made, with locations and temporal sampling chosen specifically to reduce errors in the forecast.

## 9. Conclusions

The problem of forecasting the behaviour of the Antarctic ice sheet remains challenging. We have described the development and testing of a glaciological data assimilation scheme. This can take advantage of the many observations that are now being made in Antarctica. We hope this will eventually form a basis for useful forecasts of the ice sheet. We have demonstrated through numerical simulations that our algorithm can provide estimates of the past, present, and future state of the ice sheet, along with error estimates. At present though, this is only possible under a restricted hypothesis that requires certain assumptions, implicit in our model of the ice sheet flow, to be satisfied. Our approach is not yet capable of providing accurate forecasts of strongly nonlinear behaviour, or of unstable retreat.

Even if the main assumptions of our simple ice sheet model are satisfied, measurement errors and snowfall fluctuations can introduce considerable uncertainty. This applies to estimates of how the Antarctic ice sheet has contributed to sea-level over the past century, and to forecasts of its contribution over the next. Our data assimilation scheme provides a method for quantifying these uncertainties, and hence provides the machinery for testing a number of hypotheses concerned with the future behaviour of the ice sheet.

Further work will be needed to extend our approach, so that it can be applied under less restrictive assumptions. We have given some indications of the challenges that lie ahead. Many of these could be addressed by adopting methods from numerical weather forecasting.

In time, data assimilation schemes for ice sheets will be developed further, more observations will be made, and further research will improve our modelling of processes operating beneath Antarctic glaciers and at their coastal termini. No doubt these advances will spur each other, as they have for weather forecasting. Indeed, the morning weather forecast should provide routine encouragement as we seek to predict a highly nonlinear physical system, rife with instabilities, uncertain parameters, and mathematical intractabilities. We hope that a similarly determined approach to forecasting the ice sheets will provide the predictions of sea-level that are needed for planning flood defences.

## Acknowledgements

We would like to thank the reviewers and editor for their useful comments and suggestions.

## Footnotes

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

- © 2006 The Royal Society