This paper deals with present day gravity changes in response to the evolving Greenland ice sheet. We present a detailed computation from a 3 D thermomechanical ice sheet model that is interactively coupled with a self gravitating spherical viscoelastic bedrock model. The coupled model is run over the last two glacial cycles to yield the loading evolution over time. Based on both the ice sheet's long term history and its modern evolution averaged over the last 200 years, results are presented of the absolute gravity trend that would arise from a ground survey and of the corresponding geoid rate of change a satellite would see from space. The main results yield ground absolute gravity trends of the order of ±1 µgal yr−1 over the ice free areas and total geoid changes in the range between −0.1 and +0.3 mm yr−1. These estimates could help to design future measurement campaigns by revealing areas of strong signal and/or specific patterns, although there are uncertainties associated with the parameters adopted for the Earth's rheology and aspects of the ice sheet model. Given the instrumental accuracy of a particular surveying method, these theoretical trends could also be useful to assess the required duration of a measurement campaign. According to our results, the present day gravitational signal is dominated by the response to past loading changes rather than current mass changes of the Greenland ice sheet.
We finally discuss the potential of inferring the present day evolution of the Greenland ice sheet from the geoid rate of change measured by the future geodetic GRACE mission. We find that despite the anticipated high quality data from satellites, such a method is compromised by the uncertainties in the earth model, the dominance of isostatic recovery on the current bedrock signal, and other inaccuracies inherent to the method itself.
The bedrock adjustment caused by the changing load of evolving ice sheets has become a subject of great interest because of its now well established coupling with ice dynamics and its potential as a proxy of past and current ice sheet evolution (Oerlemans 1980; Le Meur & Huybrechts 1996, 1998; Tarasov & Peltier 1997). As a consequence, bedrock displacements are generally computed within ice sheet models in order to reproduce the specific ice/Earth dynamics using a broad range of methods ranging from simple parametrizations to elaborate coupled ice/bedrock models (Le Meur & Huybrechts 1996). Although often neglected among glaciologists, the gravitational perturbation associated with the process of bedrock adjustment is also of interest, as can be judged from the profuse literature following pioneering work some 60 years ago (e.g. Vening Meinesz 1937). Only recently were gravity changes given more consideration by glaciologists because of their role as a potential proxy for the current state of balance of the ice sheets and because they provide a wealth of information on the isostatic process itself (James & Ivins 1998; Bentley & Wahr 1998). A crucial problem in the interpretation of gravity signals is the ability to distinguish between the effects of current mass changes and the contamination caused by postglacial rebound as recorded in the isostatic memory of the bedrock.
Gravitational changes induced by an evolving ice sheet mainly originate from superficial mass exchanges between the ocean and the ice sheets, and internal mass displacements in the underlying Earth. Other geodynamic changes such as mantle convection also take part in these gravity changes, but they operate at a much larger timescale (characteristic times of the order of 106 yr) and therefore can be neglected over the memory period of the viscoelastic Earth response. The most straightforward way to quantify gravity changes is to feed an Earth response model with a plausible loading scenario (ice and water) in a forward computation scheme. This is basically the approach followed here except that the loading scenario is not prescribed beforehand but modelled through time along with the bedrock evolution. The big advantage is that the effects of bedrock adjustment on ice sheet dynamics can be properly taken into account as was already done in several previous studies (Le Meur & Huybrechts 1996, 1998; Tarasov & Peltier 1997). The properties of our bedrock model further allow one to split the gravitational signal from the deforming Earth into a long term component driven by the past history of the ice sheet, and an instantaneous one exclusively driven by the current evolution of the load distribution (Le Meur & Huybrechts 1998).
The main purpose of the calculations presented in this paper is to determine which orders of magnitude direct observations of gravity trends are likely to yield and to assess whether meaningful information regarding ice thickness changes can be retrieved. To this end, we present a model based example of the current total gravity trend as it could be obtained from a high precision gravity survey on the ice free ground. To obtain the gravitational changes exclusively driven by the evolving geometry of the ice sheet, an alternative calculation is presented based on the usual gravity corrections (Bouguer and terrain corrections). This provides more accurate results than those directly output by the bedrock model. For comparison with future satellite data all over Greenland, a different observable, namely the rate of change in geoid height, is more relevant. Although the evolving geopotential is also partly driven by ice mass changes, as is the case for time dependent gravity, it can be directly computed from the bedrock model with no loss in accuracy and the corresponding results are also shown.
The two models and how they interact
The bedrock model
The bedrock model used for this study (Le Meur 1996a,b; Le Meur & Hindmarsh 2000) belongs to the category of self gravitating spherical viscoelastic earth models (Peltier 1974; Wu & Peltier 1982; Lambeck et al. 1990; Spada et al. 1992). The spherical approach is based on a harmonic decomposition and considers the entire planet when solving the Earth's response to a surface load (Fig. 1). Unlike the half space approximation, the displacement field and the associated gravitational potential can be accurately computed from the three equations for elastic gravitational free oscillations of the Earth (Backus 1967). The mantle rheology is approximated by a Maxwell body according to the ‘correspondence principle’ (Biot 1954; Peltier 1974), yielding an ‘equivalent elastic problem’. The solution for the resulting boundary value system is then Laplace inverted (Peltier 1985) to give the time dependent bedrock response to the impulse point load under the form of dimensionless Love numbers (see for instance Peltier 1974). As a result of the inversion, these Love numbers split into an elastic term, hnE(r), and several (j=1, Nm) viscous modes, which are each characterized by a viscous amplitude hVn,j(r) and a decay time τn,j such that hn(a, t), the harmonic solution of degree n at the surface (r=a) and at time t, reads
In this expression, the time dependence consists of an instantaneous elastic term, weighted by the δ function, which is superimposed onto several time decaying viscous modes. This property is also conserved in the Green's function derived from these Love numbers. The convolution of the Green's function with a loading scenario then yields the splitting of the Earth's response into a viscous long term and an elastic instantaneous contribution.
Following the work of Spada et al. (1992), the bedrock model has been coded into a semi analytical form using the algebraic manipulator Mathematica. This is a more efficient and accurate approach as it circumvents the round off errors inherent in the stepwise numerical solution of most boundary value problems. We adopted an earth structure that comprises an inviscid core, a viscoelastic mantle and a purely elastic lithosphere of 100 km thickness, which is here assumed to be compressible (Fig. 1). Because of its compressibility, the approach followed for the lithosphere is slightly different from Spada et al. (1992). The required boundary value problem is solved here assuming uniform earth parameters for the upper 100 km, which implies a matrix as in Wu & Peltier (1982). The mantle consists of three layers, each characterized by its density, shear modulus and viscosity. Except for the lithospheric compressibility, the characteristic Earth properties used in the model are adopted from Spada et al. (1992) and are summarized in Table 1. The computation yields three different dimensionless Love numbers, hn, ln and kn (eq. 1), respectively related to the radial displacement, the tangential displacement and the gravity potential perturbation. The Earth's response for different observable quantities (displacement pattern, gravity anomalies, geoid height, etc.) can then be obtained by summing the appropriate Love numbers (or a combination of them) in a Legendre series to obtain the Green's function. The actual bedrock response to a specific load finally results from a double time–space convolution with the time–space distribution for this load (Le Meur & Hindmarsh 2000).
The ice sheet model
The ice sheet model consists of two main components which respectively describe the ice sheet system and the surface mass balance, the latter of which represents the main driving force of the system. The time dependent ice dynamics model solves the fully coupled thermomechanical equations for ice flow on a 3 D mesh and includes basal sliding as well as heat conduction in the underlying bedrock. This basically involves the simultaneous solution of conservation laws for momentum, mass and heat under appropriate simplifications, supplemented by Glen's flow law (Huybrechts & de Wolde 1999):
where τi is the vector in the i direction whose components are the stress tensor elements [τij, j=x, y, z], ρi is the ice density (910 kg m−3) and gi is the projection of the gravitational acceleration along the i direction. H is the ice thickness, q is the volume flux, v¯ is the average horizontal velocity, M is the mass balance and t is time. The thermodynamic eq. (4) accounts for vertical heat conduction, 3 D advection and heat generation by internal deformation. Here, T is temperature, V is the 3 D ice velocity, W is the internal heating and k and cp are the temperature dependent thermal conductivity and specific heat capacity of ice, respectively. The lower boundary condition is the geothermal heat flux, GH, of 42 mW m−2. The rate factor A(T) in Glen's flow law with exponent n=3 depends on the ice temperature according to the Arrhenius equation and furthermore allows for different mechanical characteristics of Holocene and ice age ice, the latter of which is made to deform three times faster for the same stress and temperature conditions. Such flow enhancement was empirically established in Greenland ice cores, and is related to a marked contrast in crystal size associated with varying concentrations of chloride and sulphate ions (Paterson 1991). ε˙ij are the strain rate components, τ′ij are the stress deviators and τ is the effective stress.
A schematic representation of the main components of the ice sheet model is shown in Fig. 2. A basic assumption is that the ice flows in direct response to pressure gradients set up by gravity. Longitudinal stresses are disregarded so that ice deformation results from shearing in horizontal planes. Sliding is of Weertman type (Weertman 1964) and restricted to areas where the basal temperature is within 1 °C of the pressure melting point. There is a free interaction between climatic input and ice thickness. Calving dynamics are not described explicitly. Instead, the contemporaneous coastline, which is a function of both eustatic sea level and local bed elevation, acts as a natural barrier to grounded ice, beyond which all ice is removed as calf ice. The treatment of the coastline allows for ice sheet expansion down to the continental shelf break during periods of maximal sea level depression in so far as the surface mass balance permits it. The horizontal resolution is 20 km and there are 31 layers in the vertical. Model input consists of bed topography, surface temperature, surface mass balance, thermal and rheological parameters and an initial state. The environmental forcing is made up of the global sea level stand and a prescribed change of background temperature, from which the mass balance components (snow accumulation and meltwater runoff) are calculated.
The mass balance model distinguishes between snow accumulation, rainfall, superimposed ice formation and runoff, the components of which are all parametrized in terms of temperature (Huybrechts & de Wolde 1999). Lacking a convincing alternative, the precipitation rate is based on its present distribution and perturbed in different climates according to sensitivities derived from ice core studies. The melt and runoff model is based on the degree day method. It takes into account ice and snow melt, the daily temperature cycle, random temperature fluctuations around the daily mean, liquid precipitation and refreezing of meltwater.
The ice dynamic model has been rigorously tested within the framework of the EISMINT intercomparison project and was extensively used to investigate the Greenland ice sheet on timescales ranging from ice sheet inception during the Tertiary to the behaviour during the glacial cycles to the response to future greenhouse warming (Letreguilly et al. 1991; Huybrechts 1996; Huybrechts et al. 1996).
Coupling of the two models
The coupling consists first of forcing the bedrock model with loading changes from the ice sheet model. These also include changes in the water loading over the ocean from both prescribed sea level forcing and ocean bottom changes (Le Meur & Huybrechts 1996, 1998). With these loading data, the bedrock model computes the corresponding new bedrock topography, which is then reinserted in the ice sheet model so that the effect of bedrock height changes on ice sheet dynamics can be fully accounted for. This is because bedrock elevation controls ice sheet surface elevation and consequently surface temperature and the surface mass balance (e.g. Weertman 1961; Oerlemans 1980; Tarasov & Peltier 1997). Additionally, bed elevation and sea level control the extent of the emerged continental platform over which the ice sheet can advance and retreat. The coupling is effectuated at a 100 yr time step. For a standard simulation over two glacial cycles, the coupled model needs about 50 hr CPU time on a CRAY C 90 computer. This computational burden precludes running a large number of numerical experiments, so that only the results from the standard experiment are discussed in this paper.
Characteristics of the viscoelastic bedrock response
At the heart of the bedrock model is the calculation of the viscoelastic response to a specified loading scenario. This response is obtained from the Love numbers computed by the bedrock model, which have to be convolved (in time and in space) with the space and time distribution of the ice/water load.
Convolution of the Green's function
The surface Green's function G(θ, t) represents the axisymmetric response of the Earth to a point impulse load at the pole (Fig. 1). It is obtained by summing a solution of the form shown by eq. (1) in a Legendre series according to
with GE(θ), the elastic term, written as
and GjV(θ, t), the Green's function time dependent expression for the jth viscous mode written as, according to eq. (1),
Here, Pn are the Legendre polynomials, θ is the colatitude between the central point load and the remote point and a/Me is a necessary scaling factor, a consequence of the dimensionless Love numbers, where a and Me are the Earth's radius and mass, respectively (see e.g. Wu & Peltier 1982). Practical computation of the resulting series implies the use of appropriate cut offs and a careful approach in the problematic computation of the elastic response at the origin (θ=0). The problem is fully addressed in Le Meur & Hindmarsh (2000), to which the reader is referred. Since the interest is in the surface response, only surface Love numbers will be considered, so that hnE, hn,jV will hereafter implicitly stand for hnE(a), hn,jV(a).
To obtain the response R(i, j, t) (the vertical displacement in metres) at any of the 83×411 nodes of our 20 km×20 km grid that covers Greenland, the Green's function G(θ, t) is convolved according to
where a discussion of the different terms and their significance can be found in Le Meur & Hindmarsh (2000). The radius of influence needed for determination of the subdomain Di,j is set so that all loading changes that occur within 1000 km of the point under consideration are taken into account. This requires us to extend the 83×141 numerical grid by 50 points in each direction. All these additional gridpoints are assumed to be at sea so that the loading changes are computed as the water depth evolution (sea level change minus bedrock change)×ρw, the water density. Despite the overall improvement, this is partly in error for the northwestern part of Greenland since the loading over the nearby north American continent was certainly different. However, the fast decreasing amplitude of the response with distance is believed to reduce this error to an acceptable level.
Time dependent properties of the bedrock response
which introduces the splitting between the elastic and viscous terms at time t. Note that since Δx=Δy, these two terms can be replaced by Δx2. The first part of the right hand side of eq. (8) can be further simplified according to the properties of delta functions as follows:
which shows how the elastic term is only driven by L(i, j, t), the current state of loading. Conversely, the viscous contribution accounts for all of the past loading contributions as expressed by the time integral.
Time integration of the viscous term
In practice, the second term of eq. (8) is only integrated over a memory period of Tmem=30 000 yr, which is sufficient to approach the exact viscous solution to within less than 2 per cent. This leads to the following numerical representation:
where k is the time index encompassing the entire memory period Tmem at a resolution of Δt=100 yr. Each of the resulting integrals is evaluated analytically, which according to the expression for GjV in eq. (6c) yields the following viscous response at point (i, j):
where Nh is the chosen harmonic cut off (Le Meur & Hindmarsh 2000). It should be noted that the time discretization used here does not account for the viscous relaxation driven by load changes occurring during the last 100 yr. An alternative computation was proposed by Ivins & James (1999) in which linear segments are considered between key epochs. Their formulation is still compatible with an analytical integration similar to eq. (11) and has the advantage of integrating the viscous contribution from very recent load changes (see their eq. 33). However, given the average viscous relaxation time of the order of several thousand years, the relaxation process is completed by less than 5 per cent after 100 yr and only sudden drastic changes during the last 100 yr would lead to significantly different results. Careful inspection of our ice sheet time series did not reveal such features, at least not over areas large enough to have a serious impact on the results.
The elastic term as a function of the current load
The elastic term has the same spatial properties but a much simpler time dependence. From eqs (9) and (6b), it can be expressed as
Obtaining the bedrock response rate of change
The expressions as given above refer to the bedrock response to loading changes with reference to an initial state where isostatic equilibrium is assumed. In order to obtain the current time trends, differentiation with respect to time is necessary. For the elastic term, as can be seen from eq. (12), the corresponding time derivative implies the same formula, where the current loading rate of change d[L(i1, j1, t)]/dt=L˙(i1, j1, t) replaces L(i1, j1, t). Because of the discrete character of the loading history function, a similar time derivation would not be meaningful to obtain the viscous response. Instead, as in Le Meur & Huybrechts (1998), the viscous trend R˙V(i, j, t) can only be calculated by replacing L(i1, j1, k) in eq. (11) by a finite difference equivalent, [L(i1, j1, k+1)−L(i1, j1, k)]/(Δt). Hereinafter, unless specified otherwise, all bedrock components (displacement, gravity changes) are to be understood as trends or time rates of change, and as such are denoted with an overdot. Moreover, previous examples made exclusively use of the hn bedrock Love number (related to vertical displacement). It is clear that these derivations also apply to any combinations of any other Love numbers (such as kn, for instance) as required for the computation of the different gravity trends.
The present day gravity trend
We first concentrate on the gravity rate of change, that is, the rate of change of the norm of the gravitational acceleration vector g that, for instance, a gravimeter would record. The geoid change is considered in Section 6.2. As demonstrated in Appendix A, possible changes in the horizontal component of the acceleration vector can be neglected, so that the gravity anomaly reduces to the radial gradient of the gravitational potential, −∂Φ/∂r (the same applies for the time trends). This simplification is not a priori obvious since in some gravity problems the horizontal component contributes to the ‘terrain correction’, for instance when the topography is steep (Turcotte & Schubert 1982). If our statement holds here, it is partly because the gravity perturbation is small as the present state is rather similar to the initial reference state. As previously mentioned, the gravitational changes originate from both the Earth deformation and the load distribution changes.
Surface gravity changes resulting from Earth deformation
The gravitational potential associated with the deformation of the planet is directly obtainable from the bedrock model. Its radial gradient can therefore be computed using the appropriate Love number combination (g0/a)(n+1)kn as in Wahr et al. (1995) or in James & Ivins (1998), but without the 1/2 term that results from the loading change contribution itself and that we compute separately. The corresponding rate of change is therefore obtained by substituting this combination for hn in eqs (11) and (12) and by applying the correct time differentiation as described in Section 3.2.3.
Surface gravity changes from the loading changes
The second part of the gravity trend follows from the direct acceleration provided by the mere presence of the current load. However this ‘pure loading contribution’ is ambiguous because the gravitational potential cannot be defined exactly at the Earth's surface where the point load is applied. Most authors use an artefact consisting of giving an artificial thickness to the point load and compute this potential at the exact middle (so that gravitational effects from the upper and lower parts of the load cancel each other, James, personal communication, 1999). This leads to (g0/a)×(1/2) as the corresponding Love number combination, which only intervenes in the elastic part since it expresses the current state of the load (Wahr et al. 1995; James & Ivins 1998). Unfortunately, owing to the bedrock model first order boundary conditions expressing the presence of the point load at the surface of the planet (Wu & Peltier 1982), this combination underestimates the gravitational contribution of the load, especially when loading changes occur at different heights from that of the measuring point. This point is fully addressed in Appendix B, to which the interested reader is referred (see also Agnew 1983). The gravity contribution from loading changes is therefore deliberately omitted in the Love number combination, but is computed afterwards independently of the bedrock model, following techniques traditionally used for topographic corrections in gravimetry.
Alternative computation of the load induced anomaly
The alternative computation has two steps, the first equivalent to a Bouguer gravity anomaly accounting for the gravitational effects due to the loading changes at the local grid cell and the second including the remote effects from more distant loading changes. The Bouguer term is computed by approximating the local cell by a disc of similar area (Rc=Δx/), and by integrating eq. 5 109 from Turcotte & Schubert (1982, p. 216), in which the thickness h is now the local load variation L˙(i, j) for 1 yr (in kg m−2) divided by the appropriate uniform ice density ρi, and for which the observation point is at the surface (b=0). The corresponding vertical acceleration rate reads
With a maximum load variation of the order of 1 m yr−1 and a density of 910 kg m−3 for the ice, the last two terms can be neglected, reducing eq. (13) to 2πGL˙(i, j), equivalent to the first term of eq. (5) in Wahr et al. (1995). Provided the grid size is large enough (of the order of a kilometre), this approximation holds and approximating the cell area with that of the corresponding disc does not alter the result. For the terrain term, we simply use the law of gravitation again by stating that a remote grid cell (i1, j1) of area Δx2 whose upper surface S(i1, j1) undergoes the loading change L˙(i1, j1) will determine a gravity change at point (i, j) as follows:
where D is the distance between the points, and u is a unit vector pointing towards (i, j) at elevation S(i, j). After projecting on the vertical by multiplying with cos(u, g0)=[S(i, j)−S(i1, j1)]/D and summing over all of the remote points in Di,j within a radius of influence R, we obtain the terrain correction in (i, j) as
Because of the proportionality in 1/D3, this term becomes insignificant for D larger than a few tens of kilometres. No correction for the motion of the measuring point is accounted for at this stage. This will be considered under a free air correction term in Section 5.3.5
Results of the simulation
This section describes the main characteristics of Greenland's evolution during the last two glacial cycles with a special emphasis on the derivation of the present day gravity trends.
Experimental set up
Because of the long response timescales of both the ice sheet and the underlying bed (of the order of 103–104 yr), it is essential to start the calculations at a time early enough for the coupled model to forget its initial start up conditions. To this end, the ice sheet and bedrock models are first run over the last two glacial cycles (Fig. 3). A steady state run under interglacial conditions served as initial conditions. To obtain this reference state, it was assumed that the observed present day bedrock was in isostatic equilibrium with the observed present day ice and water loading. The effects of this assumption, although only a first approximation, turned out not to be very crucial for the model results obtained at the end of the simulation.
The model was forced over the last 225 000 yr by prescribing a temperature change derived from the GRIP δ18O record and imposing a sea level history derived from the SPECMAP stack at a 100 yr resolution (Dansgaard et al. 1993; Imbrie et al. 1984). Although the GRIP δ18O record is known for its defects prior to about 100 kyr BP, these have a negligible effect on the present day ice sheet and bedrock evolution.
Time dependent results
Ice and bedrock height evolution
The simultaneous evolution of ice sheet volume and the corresponding mean bedrock elevation is shown in Fig. 3 (lower panel). The figure displays two complete cycles of growth and decay as documented in more detail for similar set ups in earlier papers (Letreguilly et al. 1991; Huybrechts 1994, 1996). Interestingly, the spatially averaged bedrock elevation over the Greenland continental platform appears to be well anti correlated with total ice volume. The exact correspondence breaks down however for the shorter timescales because of the time delayed viscous response, as is evident from the time evolution over the last 25 kyr displayed in Fig. 4. At present, the ice volume is almost stationary, but residual bedrock uplift still occurs in response to the strong deglaciation between 10 and 5 kyr BP. From the evolution of ice volume, it turns out that the Greenland ice sheet basically completed its retreat from the Last Glacial Maximum some 5000 yr ago and is now roughly in equilibrium with the present Holocene climate. The volume changes over the last 5 and 1 kyr are respectively −4.4×1013 and +3.3×1012 m3, corresponding to global sea level changes of +11 and −0.8 cm, respectively. Those contributions are small compared to the total Greenland volume and sea level change of about −6×1014 m3 and +1.4 m since the Last Glacial Maximum.
The associated ice sheet geometries have been compared against available geological and glaciological (palaeo ) field evidence where possible (Funder 1989). In particular, the retreat history of the ice sheet in central west Greenland is well constrained by field data. Van Tatenhove et al. (1995) have shown that the model is in reasonably good agreement with a succession of dated moraines along a transect parallel to 67°N extending from the present ice margin down to the continental shelf break. Not only did the maximum and minimum extents of the model ice sheet coincide well with the geology, but also the chronology of the modelled retreat agreed to within 500–1000 yr of the glacial–geological reconstructions, or about the uncertainty on the 14C ages of the dated moraines. We take this as an important validation of the model and the loading history it is able to provide, and believe that this lends more robustness to the viscous bedrock evolution the model produces for the present time.
Mean gravity change over the last 2 kyr
The simultaneous evolution of the average load and mean bedrock displacement over the continental platform enables us to calculate an approximate average gravity change in the infinite sheet approximation as a function of time,
where H¯(t) and b¯(t) are average ice equivalent thickness and mean bedrock elevation over the continental platform (above the −300 m contour), respectively, which are taken as the difference with respect to the initial state (−225 000 yr). ρe is a mean representative density for the outer Earth (3350 kg m−3) and ρi is the ice density (910 kg m−3). can be considered as a good indicator for the average degree of isostatic disequilibrium. It is zero by definition at −225 000 yr, when isostatic equilibrium was imposed as an initial condition.
Fig. 4 shows how this gravity anomaly stems from both ice volume and ensuing bedrock changes. The effect of the ice loading is dominant because the mean thickness change with respect to the initial value of 1272 m is generally pronounced. However, once the volume stabilizes over the last 5 kyr, the ongoing uplift makes up for almost all of the gravity trend, which still amounts to about 0.28 µgal yr−1 for the present day. This corresponds to a mean uplift of about 2 mm yr−1.
The remaining present day anomaly of −3.47 mgal is almost exclusively due to the average bedrock depression of 23 m (−3.23 mgal), whereas the mean ice thickness difference of about 6.5 m accounts for only −0.24 mgal. It should be stressed that these values correspond to a specific bedrock rheology (viscosity profile, lithospheric thickness) among many possible rheologies. Because of the calculation burden of our coupled approach, no additional experiments were performed to test the sensitivity for a plausible range of Earth parameters. Despite this restriction, we believe we can already provide a first estimate of the order of magnitude of the gravity anomaly. A more accurate computation of the gravity anomaly pattern as computed from the coupled ice sheet/bedrock model is presented below.
Present day evolution patterns
The definition of the time period over which present day patterns should be calculated deserves some comment as it remains ambiguous and somewhat arbitrary. Theoretically, and with respect to the elastic crustal response, it should be the real instantaneous change occurring at present time t, but for numerical and technical reasons (the time step in the calculations, discontinuous forcing), the ice sheet model cannot yield a meaningful instantaneous trend. Moreover, in reality, the relevance of such an instantaneous trend is questionable because a strong interannual to decadal variability in the surface mass balance generally overrides a more significant longer term ice sheet dynamic trend. In fact, the largest volume change of the Greenland ice sheet occurs between the beginning and the end of the summer season, when around 50 per cent of the total annual accumulation over the entire ice sheet is melted from the ablation area. In this study we follow the same approach as in previous analyses (Huybrechts 1994; Le Meur & Huybrechts 1998) and average the model outputs over the last 200 yr to obtain the present day evolution, seen as a fair compromise between the typically strong interannual to decadal variability, the time resolution of the external forcing (100 yr) and the relevant physical processes. The implication is that possible imbalances resulting from mass balance changes within the last century are discarded (or effectively cancel one another).
Ice thickness evolution
With the above definition in mind, the coupled ice sheet/bedrock model yields a Greenland ice sheet that as a whole is almost stationary. Over the last 200 yr, the corresponding mean ice thickness change is around −65 cm and the mean bed uplift about +35 cm only (Huybrechts & Le Meur 1999). This is equivalent to a global sea level rise of only +0.15 cm century−1 or an ice volume change of −6 km3 yr−1. Despite a near overall equilibrium of the entire ice sheet, the geographical pattern shows a clear distinction between a general thickening of the accumulation area and a mostly thinning ablation area (Fig. 5c). Current ice thickness changes are highest over southern Greenland, with rates in excess of 20 mm yr−1, whereas thinning rates locally go up to −100 mm yr−1, especially in the SW and NE parts of the ice sheet. The single most important explanation for this pattern is the recovery of the ice sheet from the Little Ice Age cooling, which ended about 200 yr ago, leading to both higher accumulation and higher runoff rates. Superimposed on this pattern are the effects of basal warming following the last glacial–interglacial transition, the downward propagation of the harder Holocene ice and heat conduction into the bedrock, as discussed in more detail in Huybrechts (1994).
Bedrock uplift is primarily driven by past and current loading changes over both the Greenland continent and the surrounding ocean. The longer term loading history is well represented by the loading difference between the Last Glacial Maximum around 18 000 yr BP and the present time (Fig. 5b). The pattern is dominated by marginal ice sheet retreat, particularly pronounced in the SW and NE parts of the ice sheet, and by a thickening of several hundred metres over the central accumulation area.
Figs 6(a) and (b) show the corresponding viscous and elastic bedrock uplift rates as discussed in Le Meur & Huybrechts (1998). They are basically smoothed imprints of the past and present loading changes shown above. In this model experiment, the viscous response clearly dominates over the elastic response, with a maximum viscous uplift rate of 6.25 mm yr−1, which is an order of magnitude larger than the maximum elastic uplift rate of 0.48 mm yr−1. Also clearly noticeable is the regional character of the viscous bedrock response, which takes place in the asthenosphere underneath the large scale bending of the lithospheric rigid plate. The slightly more local aspect in the elastic response probably results from the instantaneous response of the lithosphere, whose compressibility would lead to a more local imprint of the small scale variations from the overlying load.
Gravity anomalies resulting from Earth deformation
Figs 6(c) and (d) show the gravity trends as they can be derived directly from the bedrock model. These fields are calculated from eqs (11) and (12) using the appropriate Love number combination (g0/a)(n+1)knV and (g0/a)(n+1)knE for the viscous and elastic responses respectively (see Section 4.1), and in which a time differentiation replaces the loading function L. Here the gravitational effect of the load has not been included, which makes this gravity trend exclusively representative of the Earth deformation. The similarity in shape between these fields of gravity change and the respective uplift patterns is therefore striking. It illustrates the more general correspondence between the total displacement of mass within the Earth (for which surface displacements are a good approximation) and the resulting gravity changes. Due to the regional aspect of the response, the bedrock surface displacement b˙ at any gridpoint can roughly be considered as an infinite sheet of growing rate b˙, for which the corresponding gravity trend can be approximated with the ‘Bouguer formula’ as 2πGρ¯b˙, with ρ¯ a mean density of 3350 kg m−3 representative of the outflowing mantle. It therefore leads to an approximate linear relationship of 0.14 µgal mm−1 between the gravity trend (exclusively due to the Earth's deformation as produced by the model) and then surface uplift rate, a point already raised by Wahr et al. (1995) and James & Ivins (1998). The comparison between the viscous bedrock displacement and corresponding gravity anomaly pattern in Fig. 6 gives a ratio of about 0.17 µgal mm−1, which can be considered consistent with the crude approximation of the infinite sheet. These values are also in good agreement with the 0.16 µgal mm−1 obtained by James & Ivins (1998) (after they removed the −0.32 µgal mm−1 free air effect from their viscous −0.16 µgal mm−1 ratio), and with the 6.5 mm µgal−1 proposed by Wahr et al. (1995).
The elastic trends in Figs 6(b) and (d) also exhibit strikingly similar patterns but with a much lower factor of proportionality of around 0.06 µgal mm−1. This ratio again agrees well with the 0.05 µgal mm−1 found by James & Ivins (1998), but it seems difficult to invoke a similar explanation entirely based on the infinite sheet formula, because that would require far too small a density contrast at the crust/atmosphere interface.
Gravity changes induced by the load changes
This gravity trend is represented in Fig. 7 as computed according to Section 4.3. From the figures, one can see how the Bouguer term shows a very similar pattern to that for the current ice loading changes (Fig. 5c), a consequence of the proportionality between the two in the infinite sheet formula. One can also see that over the ice free ground, where no local loading changes occur, this anomaly does not exist.
The second mass correction term (‘terrain correction’) is generally about two orders of magnitude smaller than the first one, making it only interesting from a qualitative point of view. The exception is a fairly narrow band along the ice sheet margin in central west Greenland where it locally reaches 0.13 µgal yr−1. This is evident from the formulation of the mass correction term itself in eq. (15). Since it is inversely proportional to the distance D between the points cubed, the effect can only be important over short distances and is all the more pronounced as the topography is steep and the difference of altitude between the points is large (James & Ivins 1998; Dietrich et al. 1998). That is also the reason why the noticeable thickening of the ice sheet in southern Greenland does not significantly contribute to the mass correction, because the ice topography is very flat there.
Adding these two fields gives the total gravity trend due to the evolving load. Owing to one to two orders of magnitude difference between the two, the resulting sum is almost indiscernible from the Bouguer term, which is already a good representation as shown in Wahr et al. (1995) or James & Ivins (1998), so the total gravity trend is not displayed separately. A comparison can be made with the same gravity anomaly but now computed from the bedrock model with the (g0/a)×(1/2) Love number expression that was previously disregarded (Fig. 7c). Comparison of the two figures (Figs 7a and c) confirms that the latter approach (based on Love numbers) seriously underestimates the load contribution to the gravity rate of change. The reason comes from the load in the unit bedrock approach having no physical thickness (infinitely concentrated at r=a), which is unrealistic in terms of our gravity corrections (see Appendix B for justification). It confirms that a ‘notional load’ such as that embodied by the delta function used to force the unit bedrock model is not appropriate to compute the loading gravitational impact and therefore justifies the approach followed here with a separate computation.
The total gravity trend a gravimeter would measure at the Earth's surface
With today's instrumental accuracy of 1–2 µgal for the most recent FG5 absolute gravimeters (Sawasaga et al. 1995), gravity surveys usually have to be carried out for several years in order to detect trends of the order of a microgal per year. Such a requirement makes measurements over the ice or at sea technically unrealistic, and only absolute gravity surveys on the ice free ground at the periphery of the ice sheet seem capable of providing significant trends. Assuming one wants to reproduce what a gravimeter left on the ground would measure, the theoretical gravity trend (from both Earth deformation and the loading changes) has to include the effects of the vertical motion of the device under the form of a free air correction. Given a rate of vertical displacement b˙, this correction reads (−2g0/a)×b˙. This rate of displacement b˙ is directly obtainable from the bedrock model with the surface displacement Love number hn, which makes (−2g0/a)×hn the relevant combination for the free air term and gives (g0/a)×[(n+1)kn−2hn] as the final term to use in eqs (11) and (12). By multiplying this latter combination with the scaling factor a/Me, one obtains a formula similar to eq. (21) in James & Ivins (1998), except for the elastic term, where the pure loading contribution (the 1/2 term) has been deliberately omitted. The total theoretical gravity trend for any gravimeter on the ground is finally obtained from eq. (11) with (g0/a)×[(n+1)knV−2hnV] replacing hnV for the viscous term (with an appropriate finite differentiation instead of L) and eq. (12) with (g0/a)×[(n+1)knE−2hnE] replacing hnE (with L replaced by L˙) and to which the loading contribution (Bouger+terrain terms) is added.
The corresponding pattern for the free air correction rate of change is shown for the whole of Greenland in Fig. 8(a), where values at sea (computed from the geoid rate of change as if one had a gravimeter on board a ship) and over the ice (deduced from the ice thickness rate of change) are also displayed, but the latter are of little or no practical interest.
The final theoretical gravity trend is represented in Fig. 8(b). The insets on the right hand side summarize the different components over the main ice free area of west Greenland. Over this area, there is no Bouguer term and the terrain correction is only important in a narrow fringe along the ice margin. Therefore, the initial gravity trend arising from Earth deformation only needs to be corrected for the vertical uplift of the device (free air correction). This last correction is, however, significant, since it is about twice as large as the gravity trend from the model (inset c2), and moreover of opposite sign (compare insets c2 and c3). It is characterized by a −32 µgal mm−1 gradient (−2g0/a) that more than compensates the 16 µgal mm−1 gradient obtained from Section 5.3.3 for the ratio of the total gravity anomaly to the total uplift (deduced from weighted ratios representative of the respective viscous and elastic gravity trends). Based on our model output, ground absolute gravity trends in central west Greenland are typically of the order of −0.3 to −0.4 µgal yr−1, of which the largest part, apart from the free air correction, is made up of the viscous effect as discussed earlier. It confirms that absolute gravity measurements with an instrumental accuracy of 1–2 µgal require several years of continuous survey to detect signals of this amplitude.
Towards an inference of the secular evolution of large ice sheets?
The secular ice sheet evolution from the viscoelastic theory
The relevant quantity for sea level changes is the trend of ice mass change effective over at least several decades, rather than the actual evolution at exactly the present time, which is probably dominated by interannual variations in surface mass balance. In our modelling we have defined the current evolution as the ice mass trend averaged over the last 200 years, so the elastic bedrock time dependent term as considered in this study is a good reflection of the secular trend we are interested in. The caveat to make here is that our calculations only yield the century timescale background evolution resulting from changes in environmental forcing extending back into the last glacial period, but exclude the possible contribution associated with mass balance changes over the last 100 years. This effectively assumes that recent decadal mass balance perturbations are on average small compared to the ice sheet's residual response to past climate changes.
Assuming the existence of high quality observations for the bedrock response, it is therefore tempting to infer the corresponding elastic component by subtracting the viscous long term response computed by the bedrock model from the corresponding field data, and to deconvolve the result in order to retrieve the secular ice loading changes. However, owing to the regional character of the Earth's response, the deconvolution process at a given location requires integration of the bedrock elastic information over all of the area within the radius of influence around this point. As a consequence, local gravity surveys along the ice sheet margin, as for instance started in central west Greenland by Dietrich et al. (1998), are not sufficient for such a derivation.
The high resolution geopotential from future satellite missions
The lack of coverage can be overcome by satellite missions recording the time dependent geopotential. Several past missions such as Starlette (Cheng et al. 1989) and LAGEOS (Gegout & Cazenave 1993; Eanes & Bettadpur 1996) have already contributed to first estimates of the large scale geopotential rate of change by providing the first few harmonic terms (J˙2, J˙3, J˙4,…). New techniques such as those to be implemented for GRACE (Gravity Recovery and Climate Experiment, to be launched by NASA in 2001), a forthcoming low orbit satellite to satellite tracking mission, referred to as SST in USNRC (1997), are soon expected to investigate the geopotential rate of change at a much higher resolution (up to the 180° order spherical expansion term) such that the geoid height rate of change could theoretically be derived down to a resolution of the order of a few hundred kilometres. This latter field can also be computed from the bedrock model in our experiment, in which the appropriate Love number combination to apply in eqs (11) and (12) now reads (g0/a)(1+kn), together with the appropriate time differentiation for L. The corresponding results are depicted in Fig. 9. Like gravity changes, geoid changes are also controlled by both the Earth's deformation and current mass exchanges at its surface. In some areas such as in northwest and central west Greenland (Fig. 9), the effects of crustal uplift and ice sheet thinning partly compensate, and probably explain the relatively low values as compared to ice free postglacial rebound areas elsewhere. The large mass loss in central west Greenland is even responsible for an inversion of the sign of the geoid motion.
From its technical specifications (USNRC 1997), one can expect the future GRACE satellite mission to be sensitive to geoid rates of change to an accuracy of about 0.05 mm yr−1 at the scale of the major drainage basins (500 km side square) over the assumed 5 yr duration of the mission. According to the magnitude of this observable as computed here (Fig. 9), this is likely to yield discernible information.
The geopotential as a proxy for secular ice sheet evolution?
Successful application of the differencing procedure to infer the elastic part of the geoid rate of change requires a firm handle on error sources. Above all, such a method would suffer from the uncertainties in the computation of the viscous bedrock term. This uncertainty has not been rigorously addressed in the present study by performing a comprehensive sensitivity study for the full range of Earth parameters (mainly the viscosity profile and the lithospheric thickness) cited in the literature. However, a recent study on the Antarctic ice sheet by Kaufmann (2000) estimated the bedrock uncertainty to be around 0.5 mm yr−1, and the result for Greenland may be of similar magnitude. Uncertainties arising from the ice sheet model and the loading history it produces are believed to be smaller, but following Bentley & Wahr (1998) in another study on the Antarctic ice sheet, an additional error of about 0.5 mm yr−1 cannot be excluded. To this one should add the errors inherent in the observational data, as discussed more fully in Bentley & Wahr (1998). Therefore, given the inaccuracies in the viscous response as well as those from the data, the inferred elastic term obtained by subtracting the viscous component from measured data will hardly be distinguishable from zero, unless it is much larger than our simulations suggest.
Nevertheless, as pointed out by Bentley & Wahr (1998), an integrated approach combining satellite missions such as GRACE with future high precision altimetry and detailed GPS campaigns on rock exposures offers the prospect to constrain ice and bedrock models independently and reduce error bars on the data. One may thus expect that such combined studies will ultimately produce reliable trends for the evolution of the Greenland ice sheet.
This study discussed the results of a comprehensive computation of the present day gravity changes induced by the evolution of the Greenland ice sheet. These results are complementary to the bedrock surface displacement fields obtained in Le Meur & Huybrechts (1998) for the same experiment. The model run considered the ice sheet/bedrock system in the coupled mode, allowing for more reliable results as the mutual interactions between bedrock adjustments and ice sheet evolution can be fully accounted for. This study has highlighted the potential and limits of using measured gravity trends, whether on ice free ground or area wide from satellites, to infer information on the isostatic evolution in general and the problem of the current ice sheet evolution in particular. To demonstrate these issues and establish orders of magnitude, we have synthetically simulated the gravity anomaly trends one might expect as a result of the modelled past and current evolution of the Greenland ice sheet. We find that the gravitational bedrock signal is strongly dominated by its viscous component and thus by the past history of the ice sheet. Because of the uncertainties inherent in bedrock response models, we conclude that the problem of inferring the current state of balance of the Greenland ice sheet from measured gravity values is seriously underdetermined and requires, amongst other things, more constraints on the Earth's rheology. The prospect of good quality data from future space missions is likely to offer improved possibilities to address the imbalance problem of large ice sheets, but basic problems due to the interannual variability of the surface climate and the effects of atmospheric pressure changes on the measurements (Bentley & Wahr 1998) will remain. Awaiting a better idea of these sources of uncertainty, the principal merit of our results is to guide future measurement campaigns by demonstrating what orders of magnitude one can reasonably expect, to point to specific spatial and temporal patterns, and to help with the attribution of measured gravity signals to past loading changes or current mass changes.
Appendix A: Neglecting the horizontal component in the gravity anomaly trend
Let g0 be the zero order acceleration vector for the unperturbed state (our reference state) and assume that it defines the local vertical for an orthonormal reference frame (er, e, e). The perturbation in gravity that we call Δg is the opposite of the gradient in the gravitational perturbation potential Φ, which in our reference frame is expressed as
The total acceleration vector g as the sum of g0 and Δg can split into a vertical vector,
and a horizontal vector,
The square of the norm (g)2 is then the sum of those for its vertical and horizontal components (gver)2+(ghor)2. Differentiating with respect to time and dividing both sides by 2, we obtain
where g represents , the norm of the corresponding vector. Considering Δg as a low order term, we can approximate both g and gver by g0. After dividing both sides of eq. (A4) by g, we obtain
Considering that ġver and ġhor are of the same order and that ghor/g0≪1, as confirmed by the calculations (not shown), the right hand side of eq. (A5) can be approximated simply by ġver=−(∂Φ˙/∂r)er. The fact that Δg can be considered as a low order term follows from the fact that in our experiment the present day state is very similar to the initial reference state.
Appendix B: The gravitational contribution of the load as computed by the unit bedrock model
The Love number representation
The Love number expression for the direct gravitational contribution of the current load cannot be derived directly from the gravitational potential Φ given by the unit bedrock model. However, depending on where exactly we consider this contribution, it is rather straightforward to derive appropriate Love numbers. For this gravitational contribution expressed at the Earth's surface just below the load, one obtains
and for the same contribution just above the load, one obtains
The difference (g0/a)(2n+1) between these two expressions represents the harmonic expansion, which once summed in the usual harmonic series (6a), gives exactly 4πGδ(θ), where G is the universal gravitational constant. This is in fact equal to twice the ‘Bouguer’ correction for the unit point load δ(θ), the necessary correction to apply to the measurement when moving from just below to just above the load. This led several authors (Wahr et al. 1995; James & Ivins 1998) to adopt an intermediate position and use an average Love number expression for the gravity anomaly at the solid surface,
Such an approach is equivalent to giving an artificial thickness to the load and assuming that the measurement is performed exactly at the middle of the resulting layer so that load gravitational contributions from below and above exactly compensate (James, personal communication, 1999). It also means that the pure gravitational contribution of the load reduces to the Love number (g0/a)×(1/2), as (g0/a)×(n+1)kn represents the viscoelastic contribution from the Earth's deformation.
B Justification for a separate computation of the load gravitational contribution
By summing the preceding Love number expression [(g0/a×1/2)] in the harmonic series as in eq. (6a), one obtains
Replacing the resulting Legendre series by its trigonometric expression (Farrell 1972),
and noticing that g0/Me=G/a2, we eventually obtain the gravitational acceleration produced by the unit point load, gU, at the pole as a function of colatitude θ as follows:
It is interesting to note that the same expression can also be obtained directly from the law of gravitation. This is demonstrated in Fig. B1. The gravitational acceleration of a point unit load applied in P (θ=0) at a remote point A at colatitude θ can be represented as a vector in the direction of the pole with an amplitude G/D2. Here, D is the distance (A–P) between the point load P and the remote point A. The vertical projection of this vector yields the downward gravitational acceleration equal to [G sin(θ/2)]/D2. Using the sine rule to relate D to the Earth radius a,
and noting that sin θ=2 sin(θ/2) cos(θ/2) enables us to express 1/D2 as 1/[4a2 sin2(θ/2)] and finally to obtain the same expression as in eq. (B6) for the downward gravity component. Whilst this similarity gives justification for using (g0/a)×(1/2) for the load contribution, it also reveals the inaccuracy of such an approach. Indeed, when applying the law of gravitation, we considered that both points (A and P) were exactly at a distance r=a from the centre of the Earth. The result would have been totally different in the case of a difference of altitude between the two points, especially if the gravitational effects are pronounced when these points lie close to each other. The same result also follows from the Love number approach as a consequence of the boundary conditions for the unit bedrock model (Longman 1962; Farrell 1972). These boundary conditions are expressed to first order by positioning the load at r=a, without considering any surface deformation or existing topography. These two formulations only make sense for loads exactly at the same altitude, which considerably reduces this direct gravitational effect of surface masses (expressed in this way, the contribution from remote loads arises solely as a consequence of the Earth's curvature). A separate full computation is therefore necessary in order to account properly for the exact location where the mass changes occur.
We thank Professor Heinz Miller for stimulating us to undertake this research. PH would like to thank the Belgian Impulse Programme on Global Change (Prime Minister's Office—Federal Office for Scientific, Technical and Cultural Affairs, contract CG/DD/09D) and the EU Framework IV project ‘Climate Change and Sea Level’ (ENV4 CT095 0124) for their support. We also thank the two anonymous reviewers for their suggestions that markedly improved the presentation of the paper.