Ongoing glacial isostatic contributions to observations of sea level change

Summary

Studies determining the contribution of water fluxes to sea level rise typically remove the ongoing effects of glacial isostatic adjustment (GIA). Unfortunately, use of inconsistent terminology between various disciplines has caused confusion as to how contributions from GIA should be removed from altimetry and GRACE measurements. In this paper, we review the physics of the GIA corrections applicable to these measurements and discuss the differing nomenclature between the GIA literature and other studies of sea level change. We then examine a range of estimates for the GIA contribution derived by varying the Earth and ice models employed in the prediction. We find, similar to early studies, that GIA produces a small (compared to the observed value) but systematic contribution to the altimetry estimates, with a maximum range of −0.15 to −0.5 mm yr⁻¹. Moreover, we also find that the GIA contribution to the mass change measured by GRACE over the ocean is significant. In this regard, we demonstrate that confusion in nomenclature between the terms ‘absolute sea level’ and ‘geoid’ has led to an overestimation of this contribution in some previous studies. A component of this overestimation is the incorrect inclusion of the direct effect of the contemporaneous perturbations of the rotation vector, which leads to a factor of ∼two larger value of the degree two, order one spherical harmonic component of the model results. Aside from this confusion, uncertainties in Earth model structure and ice sheet history yield a spread of up to 1.4 mm yr⁻¹ in the estimates of this contribution. However, even if the ice and Earth models were perfectly known, the processing techniques used in GRACE data analysis can introduce variations of up to 0.4 mm yr⁻¹. Thus, we conclude that a single-valued ‘GIA correction’ is not appropriate for sea level studies based on gravity data; each study must estimate a bound on the GIA correction consistent with the adopted data-analysis scheme.

Satellite geodesy, Sea level change, Time variable gravity

1 Introduction

A number of recent studies have investigated the sea level budget using three complementary measurements: altimetry (e.g. Jason-1), gravity (Gravity Recovery and Climate Experiment, GRACE) and thermosteric variations (Argo) (Willis et al. 2008; Cazenave et al. 2009; Leuliette & Miller 2009). The ongoing effects of glacial isostatic adjustment (GIA), that is, the continuing response of the viscoelastic Earth to the loading from the ice age, impact the first two of these measurements. In this regard, the contribution of GIA to altimetry is generally cited as −0.3 mm yr⁻¹ (thus, subtracting the GIA contribution implies adding 0.3 mm yr⁻¹ to the observed altimetry rate), following the value derived by Peltier (2001). Peltier (2009) found same value based upon ICE-5G(VM2) (Peltier 2004), but a slightly more negative value of −0.32 mm yr⁻¹ when averaging over a reduced latitude range of ± 66°. However, this rate is sensitive to the Earth model employed in the analysis, and the plausible range of values has not been established.

The impact of GIA on the estimates of ocean mass change derived from GRACE has caused much debate. Indeed, ocean mass balance studies have used estimates near −1 mm yr⁻¹ (Willis et al. 2008; Leuliette & Miller 2009) or −2 mm yr⁻¹ (Cazenave et al. 2009). These values are based on GIA model predictions developed by Paulson et al. (2007) and Peltier (2004), respectively. The significant discrepancy in the value adopted in previous studies is surprising given that the GIA predictions were both derived using the ICE-5G ice model with some form of VM2 Earth model (Peltier 2004; Paulson et al. 2007). Peltier (2009) derived a value of −1.8 mm yr⁻¹ from the ICE-5G(VM2) (Peltier 2004) and explored the sensitivity of this estimate to different smoothing values and exclusion of particular spherical harmonic components of the model prediction. It is interesting to note that several of the studies (Leuliette & Miller 2009; Cazenave et al. 2009; Peltier 2009) have claimed closure of the sea level budget using the three observation techniques despite using these significantly different estimates of the GIA contribution to the mass estimate derived from GRACE. In the effort to better constrain the mass flux into the oceans, it is important to understand the range of uncertainty in the GIA contribution to the GRACE observation.

This paper addresses three issues. First, GIA studies of sea level often use the term geoid interchangeably with absolute sea level or sea surface. We begin by reviewing these GIA calculations to rigorously describe the physical meaning of the predicted quantities. This discussion clarifies the GIA contribution to ongoing changes in sea level as measured by either altimetry and gravity missions. This review also demonstrates that a GIA correction to GRACE estimates of ocean mass balance based upon absolute sea level predictions (e.g. Peltier 2004) is inconsistent with the observation it is correcting and is thus in error. Secondly, for both altimetry and gravity observations we estimate a plausible range of values associated with uncertainties in Earth and ice sheet models. Understanding the uncertainty in these predicted contributions is vital to assessing the constraints imposed by the altimetry and gravity observations. Finally, we illustrate that the GIA correction to a GRACE estimate of ocean mass change will vary significantly depending upon the analysis techniques and averaging regions applied to the data. We will conclude that a universally applicable ‘GIA correction’ to mass change measurements over the ocean is neither possible nor appropriate.

2 Terminology

Studies of global sea level aim to quantify changes to both the total ocean volume and mass. If the Earth was rigid and the observations were made in a well-realized reference frame, then altimetry would measure changes in the sea surface (or ocean volume), while GRACE would measure changes in ocean mass. In practise, neither of these conditions is met. However, part of Earth′s non-rigid character is frequently incorporated into such studies. For example, GRACE observations of geopotential change are generally reported in terms of change in equivalent water height (EWH) (e.g. Wahr et al. 1998, eq. 14). This conversion provides estimates of the geopotential change in EWH as a result of changes in the surface load (from water and atmospheric pressure) at the Earth′s surface, after allowing for the Earth′s elastic response to the changing surface load.

Beyond this elastic response of the Earth to the changing distribution of the surface load around the world, GIA causes a response in both the solid Earth and the ocean. In particular, continuing deformation and flow of the crust and mantle beneath the ocean changes both the seafloor and the equipotential in the region. Note that because the mass change associated with GIA in oceanic regions is primarily due to the flow in the Earth, which has a much greater density than water, the conversion of geopotential change to EWH greatly amplifies the size of the change. To account for the GIA contribution to both of these measurements, predictions of both sea surface and geoid change are needed. However, as we will describe below, the GIA literature that uses the sea level equation frequently interchanges the terms ‘absolute sea level’ and ‘geoid’; a problem that is highlighted when comparing and discussing predicted GIA contributions to altimetry and GRACE observations. We review the relationship between ‘absolute sea level’ and ‘geoid’ in this context.

Classically, GIA models have focused on predictions of sea level because many of the time-series used as constraints in GIA modelling are from (broadly defined) paleoshoreline data. [A comprehensive discussion of the general concepts involved in the prediction of GIA-induced sea level changes, and the first modern theoretical treatment of these changes, may be found in the canonical work of Farrell & Clark (1976).] Given the relatively long time scale of GIA, sea level variations driven by this process are predicted under the assumption that the evolving ocean is in static equilibrium (surfaces of constant pressure and density are equipotentials). This static sea level theory treats the sea surface as an equipotential surface (i.e. no dynamic effects are taken into account). However, it is important to note that the value of the potential that defines the surface will be time dependent (e.g. Dahlen 1976; Farrell & Clark 1976), as becomes evident when one considers that sea level was over 120 m lower at the Last Glacial Maximum. It is the time dependence of the potential value that has lead to confusion of terminology in the past. Note that we also assume the density of water is constant, both spatially and temporally. These assumptions will hold throughout the paper.

The sea level predictions in GIA literature are typically a measure of ‘relative sea level’, S_R(θ, ϕ, t_j), which is a globally defined field at co-latitude θ, east-longitude ϕ and time t_j. As an example, if tide gauges could be deployed globally and were only affected by GIA, they would observe , where the dot indicates a time derivative and t_p is the present-day time (see Fig. 1a). The ocean is bounded by two surfaces, the sea surface and the crust, and relative sea level refers difference between these two boundaries (see Fig. 1b)

1

In this equation, S_A(θ, ϕ, t_j) is the absolute sea level or sea surface and R(θ, ϕ, t_j) is the height of the solid surface, assuming that these are measured relative to a common datum (e.g. the centre of mass of the Earth system, CM). Examples of and are shown in Figs 1(c) and (d). As indicated by the time derivative of eq. (1), in Fig. 1 the panels (a) = (c) − (d) The shoreline at t=t_j is a location where S_R(θ, ϕ, t_j) = 0, or, alternatively, where the height of the sea and solid surfaces are the same.

Figure 1

Example model predictions of the ongoing GIA contributions to geodetic measurements, showing the (a) relative sea level change (tide gauges), (c) absolute sea level change (altimetry), (d) crustal motion (GNSS receivers), (e) equipotential change observed on Earth, (f) equipotential change observed by GRACE, (g) change in equivalent water height (GRACE). These examples use ICE-5G and VM2 (Peltier 2004) as input models. (b) Cartoon illustrating the crust and sea surface at an initial time (t_j−1, solid black line) and a later time (t_j, solid red line). The dotted red line illustrates the new position at time t_j of the potential surface that corresponded to the sea surface at time t_j−1. (h) Equivalent to (g), but for a Earth model with a lithospheric thickness of of 120 km, ν_UM= 0.8 × 10²¹ Pa s, and ν_LM= 3 × 10²¹ Pa s.

In GIA models, the solution to eq. (1) is obtained using the sea level equation introduced by Farrell & Clark (1976). The central principle to these solutions is that the volume of water in the Earth′s system must be conserved. Thus, the change in S_R between times t_j−1 and t_j integrated over the oceans must be equal to water that entered or left the oceans during that time period. The only way that this conservation of mass can be achieved is if the sea surface, the upper boundary of the ocean, moves from one equipotential to another. For GIA at present time, we assume that there is no water currently being exchanged between the oceans and the ice sheets, and Fig. 1(a) sums to zero.

Return to Fig. 1(b) for an illustration. We denote the equipotential associated with the sea surface at t=t_j−1 by ϕ(t_j−1), and the perturbation to the potential, from t_j−1 to t_j, on this surface as δϕ(θ, ϕ, t_j). Following equilibrium theory, the change in absolute sea level across the time interval is given by

2

where g is the surface gravitational acceleration. The first term on the right-hand-side of eq. (2) represents the shift in the radial position of the original equipotential ϕ(t_j−1) from t_j−1 to t_j, and the new position of ϕ(t_j−1) at t_j is represented by the dotted red line in Fig. 1(b). This term includes effects due to changes in mass distribution and changes in the centrifugal potential resulting from changes in Earth rotation. All changes in mass in a given region, such as the ocean, are reflected in this term. However, taking the average of this term over the Earth must equal zero because the Earth system did not lose or gain mass. For our present-day example, the time derivative of δϕ(θ, ϕ, t_p) is shown in Fig. 1(e).

The spatially invariant second term of eq. (2) is included because, as described earlier, the sea surface does not have to remain on the same equipotential as a function of time [i.e. Δϕ(t_j) =ϕ(t_j) −ϕ(t_j−1)](Dahlen 1976). The integrated change of δϕ(θ, ϕ, t_j) and δR(θ, ϕ, t_j) over the oceans may not equal the change in volume associated with water entering or leaving the ocean. Thus, the sea surface must shift to a different potential value, and this shift is represented by Δϕ(t_j). In other words, assuming a constant density of water, the change in mass of the water in the ocean must be balanced by the change in the difference between the sea surface height and crustal height over the ocean. The final position of the sea surface is illustrated in Fig. 1(b) by the solid red line. Rather than a change in potential due to changes in mass, which is completely represented by δϕ(θ, ϕ, t_j), Δϕ (t_j) is simply introduced to allow us to track the changes of the sea surface with time. The difference between Figs 1(c) and (e) is simply the spatially constant value Δϕ (t_p)/δt/g.

GIA nomenclature that utilizes the sea level equation has frequently equated the change in the sea surface height, the change in the absolute sea level, and the change in the geoid height. Indeed, this assumed equivalence is emphasized by the fact that, in many cases, S_A is often referred to as G. However, if the change in potential surface representing the geoid had a spatially constant term, this would represent a change in mass of the earth system. Therefore, it is important in geodetic analysis of GIA model predictions that a clear distinction is made between the change in a particular potential surface, δϕ(θ, ϕ, t_j)/g, and the change in sea surface, δS_A(θ, ϕ, t_j). The distinction becomes clear when interpreting the GIA prediction in the context of measurements made by Jason and GRACE.

Altimeters measure the sea surface height, and thus sample δS_A. However, GRACE is only sensitive to changes in mass distribution and not Δϕ. Therefore, one must not use the absolute sea level predictions in place of geoid predictions when removing GIA estimates from GRACE observations. As an example of the difference between these two predictions, the time derivative of the term Δϕ(t_p)/g can range from −0.1 to −0.4 mm yr⁻¹ for the models considered below. Assume a value of −0.2 mm yr⁻¹. If one were to include this spatially constant (spherical harmonic degree 0) term in to the conversion of geopotential change to EWH (in this case, multiplying by the ratio of the average density of the Earth to water divided by three; e.g. eq. 14 of Wahr et al. 1998), this term would add −0.37 mm yr⁻¹ EWH to the predicted GIA contribution to the GRACE measurement. Failure to recognize this distinction by Peltier (2009) accounts for roughly half of the discrepancy between his estimate of the GIA contribution to GRACE observations and the estimate derived from Paulson et al. (2007) (Chambers et al. 2010).

Using δS_A predictions directly in GRACE analyses can lead to another error if the impact of contemporaneous perturbations in the rotation vector on sea level (i.e. rotational feedback) are included. Most recent GIA predictions are computed for an observer on the rotating Earth, and they therefore include the full contribution due to polar motion. However, the direct component of this contribution (i.e. the shift in the orientation of the centrifugal potential) would not be present in the GRACE data, because the observation is not made in the rotating frame. Therefore the direct component must be removed when calculating the GIA correction to a GRACE observation (e.g. one must remove the value of ‘1’ multiplying the δ-function in eq. (A11) of Mitrovica et al. (2001)). Removing this term causes the size of the degree two, order one coefficient to be reduced by a factor of slightly greater than two; see Appendix. The impact of removing this term is illustrated by the difference between Figs 1(e) and (f). The remainder of the discrepancy between estimates based on GIA predictions from Peltier (2004) and Paulson et al. (2007), after accounting for the Δϕ(t_p)/g term described above, is due to differences in the degree two, order one terms (Chambers et al. 2010). While this factor of slightly greater than two described in Appendix explains part of discrepancy, the reason for the remaining difference is unknown and thought to be an error in this component of δS_A of Peltier (2004; Chambers et al. 2010). An appropriately-modified prediction of the δϕ(θ, ϕ, t_j)/g will be represented by δS_G(θ, ϕ, t_j) in the discussion below.

3 Modelling

To derive estimates of the current changes in relative sea level and its two bounding surfaces from eq. (1), one can solve the ‘sea level equation’ which describes how these surfaces are related and change with time. In this regard, the equilibrium theory we adopt is revised from the Farrell & Clark (1976) approach to include the influence of Earth rotation changes as well as a coastline that evolves due to either relative sea level changes or variations in the extent of grounded, marine-based ice (e.g. Mitrovica et al. 2005; Kendall et al. 2005). This theory will generate predictions of the present-day (t=t_p) secular change of absolute sea level, the radial position of the solid surface, and relative sea level due to GIA [i.e. δS_A(θ, ϕ, t), δR (θ, ϕ, t) and δS_R(θ, ϕ, t) divided by δt]. We denote these fields by and , respectively (from now on, the θ, ϕ, t_p dependence of these quantities will be assumed).

The sea level theory requires two inputs: a global model of the Pleistocene history of ice sheet change, and a model for the Earth′s viscoelastic and density structure. For the former, we adopt a version of the ICE-5G deglaciation history (Peltier 2004). Because our purpose here is to present predictions of large-scale spatial averages of the absolute and relative sea level rates, our main concern in choosing an ice history is that it be roughly consistent with the total volume of grounded ice thought to have melted from the last glacial maximum to the end of the main deglaciation phase (4 kyr BP). We will briefly discuss comparative results based on ICE-3G (Peltier 1994).

The GIA predictions of present day rates of relative sea level and its bounding surfaces are sensitive to variations in the adopted Maxwell viscoelastic Earth model. In this regard, we have generated results for a suite of spherically symmetric, self-gravitating, compressible Earth models distinguished on the basis of the thickness of a purely elastic lithosphere (LT) and assumed constant viscosity within the sublithospheric upper (ν_UM) and the lower (ν_LM) mantle: 71 km < LT < 120 km, 10²⁰ < ν_UM < 10²¹ Pa s, and 2 × 10²¹ < ν_LM < 5 × 10²² Pa s. The boundary between these latter two regions is taken to be 670 km depth. The elastic and density structure of the Earth model is prescribed by the seismic model PREM (Dziewonski & Anderson 1981). Results are presented in the centre of Earth (CM) reference frame. GIA calculations are usually conducted in the centre of solid Earth (CE) reference frame (e.g. Farrell 1972). However, the difference between GIA predictions in the CE and CM is small because this difference is due only to the ongoing changes in the ocean water distribution due to GIA (Argus 2007; Klemann & Martinec 2009). (Note, the same statement could not be made for the large present-day changes in ice volumes.)

It should be noted that an ice sheet model used in the sea level equation is derived using a specified Earth model, VM2 in the case of ICE-5G (Peltier 2004). Varying the Earth model independently of the ice sheet model may introduce larger variations in the GIA estimates than would result if the Earth model and ice loading history were varied simultaneously. However, our suite of predictions will nevertheless provide an indication of the sensitivity of the GIA predictions to variations in the Earth model. In addition, the Earth is characterised by 3-D variations in mantle properties, and the uncertainties caused by the neglect of these variations cannot be adequately captured using a suite of 1-D profiles. Thus, while the range of model predictions we investigate is relatively large, narrowing the range of possible values without consideration of the other sources of uncertainty is premature.

4 Results

In this section, we start by examining the GIA contribution to sea level change as measured by satellite altimetry, . We also will plot the predictions for and to better illustrate the physical mechanisms responsible for the corrections. Then, we consider measurements of sea level variations over the ocean as measured by GRACE and show that GIA is responsible for a significant contribution due to the nature of the measurement.

Estimates of secular ‘global sea level change’ inferred from satellite-derived data sets (e.g. TOPEX/Poseidon, Jason-1 and Jason-2) are derived from the rate of change of the absolute sea level averaged over some large geographic area. In the case of these missions, this average is obtained across oceanic regions bounded to the north by 66°N and to the south by 66°S, excluding Hudson Bay and areas with ocean bathymetry less than 120 m (in the case of Leuliette et al. 2004). A correction of this estimate for GIA should be obtained by the removal of a similarly averaged map of absolute sea level rates (Peltier 2001).

Fig. 2 shows the average values of and for two cases: the global ocean and the ocean sampled by the altimetry missions as described above. For this plot, the lithospheric thickness is chosen to be 71 km. This choice of LT has a small impact on these results: generally less than a 15 per cent change, with a thicker lithosphere producing smaller predictions for an upper-mantle viscosity above 3 × 10²⁰ Pa s. We will discuss in more detail the sensitivity of the predictions to variations in lithospheric thickness when reviewing the GIA contribution to GRACE observations.

Figure 2

Comparison of GIA contribution to altimetry estimates of sea level rise. The contours, with units of mm yr⁻¹, show the variation caused by changing the upper- or lower-mantle viscosity. All plots assume a lithospheric thickness of 71 km. The first column shows the altimetry contribution from GIA, the second column shows the average change of crustal depth over the ocean, and the final column shows the difference between the first two. The results in the first row are generated by averaging over the global ocean, and the second restricts the average to latitudes between ± 66°.

As indicated by eq. (1), the third column in Fig. 2 represents the difference between the first two columns. Because the global average of over the ocean is equal to the change in its water volume, the present-day values for this plot are zero for GIA, and the plots for and are identical. However, if one restricts the average to a smaller region of the ocean, would be non-zero. Most of the GIA-induced vertical deformation of the crust in the ocean occurs in high latitudes near the loading centres (e.g. Mitrovica & Peltier 1991). The impact of ongoing GIA on oceanic volume is mainly due to the collapse of the forebulge near these regions. Thus, as an extreme example, one could limit averages to the mid-latitudes. In this case the magnitudes for the values of would be significantly smaller and the average of would be negative as water is being transferred to higher latitudes.

The main result in Fig. 2 is the predicted average of the GIA-induced absolute sea level rate over regions sampled by the altimetry missions (bottom of first column). This average ranges from ∼− 0.15 to −0.5 mm yr⁻¹ for the Earth models we have considered. Thus, the mean (altimetry-averaged) absolute sea level rate obtained from altimetry data could be increased by 0.15 to 0.5 mm yr⁻¹ by a correction for the signal due to GIA. We note that this correction remains relatively constant as one considers averages over a larger subregion of the ocean. Thus our conclusion will hold for satellite-derived estimates with a sampling range that is different from the current Jason-type missions. It should be noted that there is only a small change in the results if ICE-3G is used instead of ICE-5G, with ICE-3G producing slightly larger values, peaking at ∼20 per cent larger for ν_LM= 2 × 10²¹ Pa s.

Next, we consider the GIA contribution to the GRACE observation of geoid height change. The GIA estimate that should be compared to the GRACE data differs from the in two significant ways. First, the geoid does not include the spatially invariant second term in eq. (2) or the direct effect of the perturbed centrifugal potential (e.g. Fig. 1f). Secondly, most GRACE results are quoted in terms of ‘equivalent water height’, the change in height of water that would cause the observed geopotential change on an elastic Earth. As discussed in Section 2, because the variations caused by GIA are due primarily to crustal motion a relatively small height variation will produce significantly larger apparent changes in EWH (e.g. Fig. 1g). Thus, the GIA contribution becomes a significant component in the global water balance estimates when reconciling total regional contributions to observations (Willis et al. 2008; Cazenave et al. 2009; Leuliette & Miller 2009; Peltier 2009).

Fig. 3 shows the variation of the GIA contribution to an ocean-averaged value of , with contour lines representing change in mm yr⁻¹ EWH, obtained in three ways. The first column assumes an average using the present-day ocean mask derived in the GIA calculations but truncated at spherical harmonic degree and order 60. (Truncation at degree 60 has a negligible impact on the averages; we refer to this mask as R= 0.) The second column applies a 500 km Gaussian smoothing to the mask (R= 500), a common processing step applied to GRACE data to reduce the impact on mass estimates of increasing error with degree (Wahr et al. 1998). However, with GRACE data, there is considerable leakage of the much larger gravity variations over the continents into the ocean averages if a border around the continents is not maintained (e.g. Chambers 2009). Thus, the third column uses a mask (B&R = 300) that maintains a 300 km border around the continents, limits the averaging to between ± 66°, and applies a 300 km Gaussian smoothing. This mask is a non-optimized version of that used in Leuliette & Miller (2009). The two rows in Fig. 3 show the results for the limiting values of lithospheric thickness.

Figure 3

Comparison of GIA contribution to GRACE estimates of sea level change over the ocean. The contours, with units of mm yr⁻¹, show how the water equivalent average varies with upper- and lower-mantle viscosity. The rows assume a lithospheric thickness of 71 and 120 km, respectively. The left column uses the present-day ocean mask, truncated to degree and order 60, from the GIA prediction; the middle column uses the same mask with an additional 500 km Gaussian smoothing; and the final column uses a non-optimised version of the mask employed by Leuliette & Miller (2009).

The variation between filters shown in Fig. 3 can be explained by considering the regions excluded. For example, B&R = 300 excludes the latitudes greater than 66° and regions within 300 km of the continent. Thus, it excludes many of the areas experiencing mass increase, such as Hudson Bay and the Antarctic ice shelves that would be included in the R= 0 mask, and this yields a much more negative value. The decrease in magnitude from the R= 0 to the R= 500 masks is due to the opposite effect. The smoothing effectively includes more of these regions of mass increase, thus lowering the average magnitude.

Generally, the estimated GIA contribution varies most strongly with lower mantle viscosity. For the mid-viscosity range, 8 × 10²¹ Pa s ≥ν_LM≤ 2 × 10²² Pa s for ν_UM= 10²¹ Pa s, the peripheral bulge surrounding Laurentia extends further into the ocean and the Earth has recovered less than for lower values of ν_LM. Thus, the average over the ocean is more negative. For the purpose of discussion, assume a lithospheric thickness of 71 km and an upper-mantle viscosity of 5 × 10²⁰ Pa s. Then, by varying the lower mantle viscosity, the resulting averages range from −0.82 to −1.72 mm yr⁻¹ using the ocean mask with no smoothing, and −0.86 to −2.10 mm yr⁻¹ for the B&R = 300 mask. As would be expected when the results are dominated by lower mantle viscosity, the lower-valued spherical harmonics contribute significantly. In particular, the degree 2 contribution is significantly enhanced in the B&R = 300 mask (between −0.3 and −0.7 mm yr⁻¹) when the higher latitudes are excluded from the average.

One can draw two conclusions from these results. First, uncertainty in Earth model could have a significant impact on the estimated contribution that would be observed by GRACE. As a visual example, Figs 1(g) and (h) show that even switching from VM2 to an Earth model that provided good results in Tamisiea et al. (2007) can produce significant changes to the average, −0.77 to −0.97 mm yr⁻¹. Secondly even if the GIA prediction were known with absolute certainty, the different averaging regions and processing techniques applied to the GRACE data could change the estimate significantly. Thus, it is impossible to have a single-valued GIA ‘correction’ for GRACE observations, as each analysis could sample the GIA contribution differently.

Fig. 4 shows how the averages change when ICE-3G (Tushingham & Peltier 1991) is used instead of ICE-5G. These two ice models are significantly different (ICE-3G has more ice in Antarctica and less in the Northern Hemisphere), but were derived using somewhat similar Earth models. The most significant differences between the results occur at higher values of ν_UM and lower values of ν_LM. In absolute terms, the differences are generally less than 0.3 mm yr⁻¹, or 25 per cent of the total.

Figure 4

Comparison of GIA contribution to GRACE estimates of sea level change over the ocean. The contours, with units of mm yr⁻¹, show how the water equivalent average varies with upper- and lower-mantle viscosity. The results in each row were generated using either ICE-5G (Peltier 2004) or ICE-3G (Tushingham & Peltier 1991). The columns are described in Fig. 3. The lithospheric thickness for all of the results is 71 km.

For all model predictions, the estimated contribution is negative. So, for example, assume that GRACE observations observed no increase in ocean mass. The results from the first column of Fig. 3 would imply anywhere from 0.5 to 1.9 mm yr⁻¹ must be contributed from other sources, say melting glaciers or ice sheets, to offset the GIA contribution to this estimate. If we were to limit the set of Earth models to those preferred in Tamisiea et al. (2007) using ICE-5G, then the range estimated would be between −0.9 and −1.5 mm yr⁻¹ for the full ocean mask with no smoothing or between −1.1 and −1.8 mm yr⁻¹ for B&R = 300. However, these estimates are strongly biased by the mantle structure under Hudson Bay, and much of the sea level change is driven by the collapse of the peripheral bulge in the ocean near the loading centres. Thus, a different average mantle viscosity may alter the estimated range of values. Therefore, it may be premature to limit the range further without explicit application of other data sets.

Finally, if one looks at only the Earth model (VM2) used to derive ICE-5G, we find a values of −0.96, − 0.77 and −1.09 mm yr⁻¹ for the three filters shown in Fig. 3. The difference between these estimates and the results of Chambers et al. (2010) using the Paulson et al. (2007) model is partially due to the filters employed in the estimation procedure and partially due to the difference in the model predictions. For ICE-3G and an Earth model with 120 km, ν_UM= 1 × 10²¹ Pa s, and ν_UM= 2 × 10²¹ Pa s, we find −1.11, − 0.84 and −1.21 mm yr⁻¹.

5 Conclusions

Efforts to estimate global mean sea level variations that might be associated with recent climate require a correction for the signal due to ongoing GIA. However, the nature of this correction depends on the sea level signal being estimated. Although the GIA contribution to both altimetry and GRACE estimates of sea level rise are related to changes in geopotential height, satellite altimetry (e.g. Jason) measures changes in the position of the absolute sea level (sea surface) while GRACE measures the total mass change over the ocean. The sea surface change caused by the crustal motion due to GIA is relatively small. However, when the geoid change is expressed in terms of equivalent water height, the larger density associated with crustal material causes the GIA contribution to become of equal magnitude to the expected variation due to water fluxes into the ocean.

In the past, confusion has resulted from the use of the term ‘geoid’ employed in the GIA literature. Because the GIA theory is mainly concerned with the sea surface and assumes equilibrium adjustment of the ocean, the equipotential surface corresponding to the sea surface is often referred to as the geoid height, though the associated equipotential value will change with time. Thus, altimetry measurements are appropriately compared to these estimated values. In the case of the current GIA contribution, the non-zero global average of sea surface change is due to change in the ocean volume due to crustal motion. However, for the estimated GRACE contribution, this non-zero global contribution must not be included because it does not represent a change of mass in the ocean; doing so would greatly increase the estimated mass change that would be observed by GRACE. Improperly including this term has caused some estimates (Cazenave et al. 2009; Peltier 2009) to be too large (Chambers et al. 2010). In addition, we have demonstrated that the degree two, order one coefficients of the sea surface, , would be a factor of two larger than those of , which can also lead to an overestimate.

Some altimetry data sets are limited to oceanic regions between ±66° latitude. Our predictions of the average absolute sea level rate induced by GIA range from −0.15 mm yr⁻¹ to −0.5 mm yr⁻¹, where the range reflects the suite of Earth and ice models we considered. We furthermore find that this range is not significantly altered if we consider averaging zones that are more or less confined in latitude. Our results for the GIA contribution to altimetry observations are in good agreement with predictions from earlier studies (Peltier 2001, 2009) generated using the same Earth and ice sheet models. The spread in our results demonstrate that the uncertainty in ice history and Earth structure will not lead to a significant increase in the uncertainty of the global mean rate of . Current uncertainties in the global mean are dominated by errors in the tide gauge calibration, and are at the level of 0.4 mm yr⁻¹ (Leuliette et al. 2004).

Unfortunately, the spread in the model predictions for the GIA contribution to the GRACE observations is quite large. Assuming a straight average over the oceans (R= 0), a range of estimated produced by varying the Earth model was −0.5 to −1.9 mm yr⁻¹, while the difference between estimates calculated using the two ice sheet models was up to 25 per cent. This range has significant implications for estimates of the present-day mass contribution to sea level rise. If GRACE observed close-to-zero mass change over the ocean, then this contribution would have to be compensated for by mass loss from ice sheets, glaciers and water stored on land. Given the importance of GRACE to help constrain the mass-flux contribution to sea level rise, it becomes vitally important to narrow this range as much as possible. Restricting the Earth models to those preferred in Tamisiea et al. (2007) reduces this range to −0.9 to −1.5 mm yr⁻¹. However, these predicted values are driven by forebulge collapse in oceanic regions, and it is not clear that a restriction based on constraints derived from Laurentian adjustment are appropriate.

Even if the ice and Earth models are known perfectly, however, our results indicate that a single-valued GIA ‘correction’ is not possible. The average value over the ocean can vary significantly depending upon the averaging area and processing techniques applied to prediction. This variation can be as much as 0.4 mm yr⁻¹. Thus, it is important that each study process the GIA prediction in the same manner as the GRACE data to appreciate how GIA will contribute to their GRACE estimate.

Acknowledgments

I would like to thank Chris Hughes, Jerry Mitrovica and Philip Woodworth for comments on the manuscript and Eric Leuliette for comments and providing his averaging mask used in GRACE processing. Two anonymous reviewers provided suggestions that greatly improved the manuscript. This work was funded by NERC as part of Oceans 2025.

Appendix

Appendix: relationship between the degree two, order one coefficients of AND

The degree two, order one coefficients of the and differ due to the fact that they are observed differently. (Note that we have shifted the subscripts A and G to superscripts in the appendix for clearer presentation when subscripts are added for the degree and order.) is observed on the Earth, and thus the direct effect of the shift in the orientation of the centrifugal potential must be included. Because is observed from a non-rotating frame, this direct effect would not be present. It is possible to derive a relationship between the degree two, order one coefficients for these two observations.

As is illustrated in Chambers et al. (2010),

3

where is change in the angular pole position, are the products of inertia, and C and A are the polar and equatorial moments of inertia. Using the normalization for the spherical harmonics presented in Matsuyama et al. (2006), can be expressed as

4

where M is the Earth′s mass and a is the mean radius. Thus, can be expressed as

5

Now, if we add the direct affect of the perturbation to the rotational potential, , to to obtain , while simultaneously subtracting the term to balance the equation, we have the expression

6

From Mitrovica et al. (2001), we can relate to

7

Thus, rearranging terms gives

8

where k_f is the fluid Love number, given by (3 G/a⁵/Ω²)(C−A). Finally, degree two, order one coefficients of and can be related by

9

The factor in the parentheses above evaluates to 2.06.