Inference of mantle viscosity from GRACE and relative sea level data

Summary

1 Introduction

Mantle viscosity is important not only for understanding the Earth's transient deformation (10−10⁵ yr) of the Earth, but also for the essential role it plays in studies of long-term mantle dynamics and the Earth's evolution. Two methods have been used to infer-mantle viscosity: studies of postglacial rebound (PGR) during the last 10⁵ yr (Haskell 1935; Cathles 1975; Peltier 1976, 1998; Wu & Peltier 1982; Mitrovica & Peltier 1995; Kaufmann & Lambeck 2002), and the interpretation of long-wavelength geoid anomalies in terms of mantle convection (Ricard et al. 1984; Richards & Hager 1984; Hager & Richards 1989). In PGR studies, mantle viscosity is inferred by fitting modelled signals to various types of observations: for example, relative sea level (RSL) histories in formerly glaciated regions, the secular change in the earth's oblateness determined from satellite laser ranging (SLR) measurements and true polar wander estimated from astrometric data over the last century (e.g. Peltier 1998), and the static gravity field in Canada and Scandinavia (e.g. Simons & Hager 1997). In long-wavelength geoid models, mantle density anomalies constructed from seismic tomography and the locations of subducted slabs, are used to load the Earth internally. The geoid anomalies caused by those loads depend on the mantle viscosity profile, and so the viscosity can be constrained by comparing with the observed geoid (e.g. Hager & Richards 1989).

The viscosity profiles inferred from these two approaches are not always consistent with one another. While the geoid studies suggest a significant (factor of 30 or more) increase in viscosity from the upper mantle to the lower mantle (e.g. Hager & Richards 1989), some PGR studies (e.g. Peltier 1998) favour a more uniform mantle viscosity, increasing by a factor of 3 or less. Other PGR studies (e.g. Han & Wahr 1995; Simons & Hager 1997) suggest that the PGR observations are also consistent with a profile where the lower mantle is significantly more viscous than the upper mantle. Studies that jointly invert both PGR and long-wavelength geoid observations suggest a more complicated mantle radial viscosity structure, but with an overall increase in viscosity towards the lower mantle (e.g. Forte & Mitrovica 1996; Mitrovica & Forte 2004). Steffen & Kaufmann (2005), who include vertical GPS data from Fennoscandia, have also shown that a viscosity model with a weak asthenosphere is the best fit to their data, but do not argue that this feature is necessarily resolvable.

Although recent studies have started to examine the effects of laterally varying (3-D) viscosity structure (Kaufmann & Wu 2002; Wu & van der Wal 2003; Zhong et al. 2003; Latychev et al. 2005; Paulson et al. 2005) and of non-Newtonian rheology (Wu 1992) on the PGR observations, an important question is to what extent, and how robustly, are the radially dependent (1-D) viscosity profiles constrained by PGR observations. This issue can be restated in more specific terms through the following questions: (1) knowing that the mantle viscosity is undoubtedly 3-D, what does 1-D radial viscosity inverted from PGR observations represent? (2) What PGR observations are most helpful in constraining 1-D mantle viscosity? (3) How much 1-D detail can the PGR observations truly resolve?Paulson et al. (2005) used synthetic studies to show that the 1-D radial viscosity inverted from PGR observables is a good representation of the horizontally averaged mantle viscosity structure below the ice loads. PGR observables in formerly glaciated areas provide the most powerful constraints on mantle viscosity. Those at the peripheral bulges are more difficult to interpret because they are sensitive to mantle viscosity both below the ice load and below the periphery. However, even with synthetic studies in which ice loading models are specified and known, Paulson et al. (2007) show that the PGR observables, including RSL data, time-variable gravity fields from Gravity Recovery And Climate Experiment (GRACE), , and true polar wander, are not able to confidently resolve more than two viscosity layers. The studies by Paulson et al. (2005, 2007) could partially explain the diverse mantle viscosity structures inferred in previous PGR studies. Kaufmann & Wu (2002) also found that 1-D model predictions are insensitive to trade-offs between lithospheric thickness and asthenospheric viscosity.

The purpose of this study is to continue to explore PGR constraints on mantle viscosity. This analysis extends the synthetic studies of Paulson et al. (2005, 2007) by using real PGR observations in and around Hudson Bay to invert for mantle viscosity. In addition to employing the RSL data that have been used extensively before, this study describes one of the first PGR applications of GRACE time-varying gravity data (see also Paulson 2006; Tamisiea et al. 2007). The secular trends in the GRACE gravity fields over North America show a feature centred around Hudson Bay, that is almost certainly caused by PGR. This feature includes two local maxima, one west and the other southeast of Hudson Bay. Tamisiea et al. (2007), published as this manuscript was in review, noted this two two-dome structure of the GRACE pattern. The GRACE trends are included here as additional constraints on mantle viscosity.

In this study, we seek the mantle viscosity structures that best fit the RSL and GRACE data in the Hudson Bay region. We find that the assumption of a homogeneous mantle unambiguously requires a viscosity between 1.4 × 10²¹ and 2.3 × 10²¹ Pa s. An inversion for two layers shows an ensemble of viscosity structures that are compatible with the data, with the lowest misfit occurring for upper- and lower-mantle viscosities of η_UM = 5.3 × 10²⁰ Pa s and η_LM = 2.3 × 10²¹ Pa s, when the boundary between the layers is at a depth of 670 km. However, a weaker upper mantle can be compensated by a stronger lower mantle, such that there exist other models which also provide a reasonable fit to the data. We also find that attempts to invert for more than two layers result in many viscosity models that fit the PGR data, thus preventing a meaningful constraint on the viscosity of any given layer. This conclusion is generally consistent with the synthetic study of Paulson et al. (2007).

2 Ice Model and Observations

Our goal is to constrain mantle viscosity by comparing PGR model predictions with RSL and GRACE observations for the region around Hudson Bay. We do not include polar wander, the static gravity field over Canada, or SLR observations of , because those could be significantly impacted by processes unrelated to PGR. could have contributions from on-going changes in polar ice, polar wander could be affected by both polar ice and mantle convection, and the static gravity field over Canada is likely to have contributions from underlying density anomalies resulting from mantle convection. These contributions are not well known, but in each case they could be a significant fraction of the PGR signal.

2.1 Glaciation model

The shape, time evolution, and amplitude of the Late Pleistocene and Holocene ice sheets are critical in determining the PGR signal. In our PGR models we employ the most recent publicly available global ice history, ICE-5G from Peltier (2004). We use this ice history, unaltered, in all model runs, and vary the viscosity profile to best match the observations. Errors in the ice history could thus lead to errors in our viscosity estimates. The RSL and GRACE data scaling described in Sections 2.2 and 2.3 below, is an attempt to minimize the impact of ice history errors, especially of errors in the ice amplitudes.

Keeping the ice history fixed while varying the viscosity, is somewhat inconsistent with the methods used to generate the ice history to begin with. In a sense, ICE-5G was obtained through iteration with the viscosity profile to best match a set of observations. The final result was the combined ice history and viscosity profile that together provided a best fit to those observations. In principle, if the viscosity profile is changed (as we will do many times below), then ICE-5G would not be the optimal ice history for matching the original set of observations. The data scaling described below partially addresses this problem, but does not completely resolve it.

2.2 Relative sea levels

A change in RSL at a given location can occur when the land surface goes up or down, or the total amount of water in the ocean changes, or there is a change in shape of the local geoid. In the region surrounding Hudson Bay, changes in RSL since the last glacial maximum are dominated by vertical rebound of the crust. The PGR signal can be large, on the order of several cm yr⁻¹ within a few thousand years of deglaciation. This is large enough that any reliable RSL observation from this region can be confidently interpreted as a PGR signal, and so provides an unambiguous constraint on the rebound.

RSL indicators, though, are often hard to decipher. Here, we use the RSL histories provided by Art Dyke (personal communication) in a large database of shoreline elevations. These data include sea levels estimated from driftwood, ancient whale parts, plant materials, and a variety of shell samples. For consistency, we follow Walcott (1980) and use only tidal mollusk shell samples. This avoids complications caused by mixing sample types. Wood, for example, tends to float and so may be inconsistently displaced from tidal mollusk markers. Of sites throughout North America, we consider only those with substantial data at postglacial times. In accordance with Peltier (1998) (and as advocated by Art Dyke, personal communication), where the elevation errors are not provided in the data set we assume the error to be 5 per cent of the elevation, with minimum and maximum values of 0.5 and 2.5 m. The carbon-14 dates provided in the database are converted to sidereal time with the INTCAL04 calibration table, which is based on a sample set of dated tree rings, U-Th dated corals and varve-counted marine sediment (Stuiver et al. 1998). Finally, we remove from the data the bias introduced by glacial meltwater entering the oceans, as given by the ICE-5G glaciation model.

All the sites considered in this study (Fig. 1) are near the centre of the Laurentide ice-sheet. Although included in preliminary experimentation, sites at the periphery of the loaded region are omitted in our final analysis. These sites tend to have small amplitudes and more complicated RSL histories due to migration of the PGR forebulge and to non-local viscosity sensitivity (Paulson et al. 2005); they also tend to have much less constraining power on mantle viscosity than do the Hudson Bay sites.

Figure 1.

Sites of RSL measurements used in this study. Place names are: (a) James Bay, (b) Richmond Gulf, (c) Ottawa Island, (d) Churchill, Manitoba, (e) Boothia Peninsula and (f) Melville Peninsula.

Open in new tabDownload slide

Predictions from PGR models are sensitive to errors in the ice history. It is desirable to compare with data in a way that minimizes the impact of those errors. To this end, Mitrovica & Peltier (1995) propose two methods for analysing RSL data: (1) scaling away their amplitude, to give a time-series for what they termed ‘normalized relative sea level’, or normalized relative sea level (NRSL) and (2) fitting an exponential form to the RSL curve and extracting the decay time (see, also, Mitrovica 1996; Peltier 1996; Mitrovica & Forte 2004). The use of decay times is associated with several difficulties, demonstrated by the plethora of highly inconsistent estimates of decay times. For example, Mitrovica et al. (2000) in a meticulous study of the RSL data obtain decay times of 5.3 ± 1.3 and 2.4 ± 0.4 kyr at Richmond Gulf and James Bay, respectively, to replace earlier estimates of 3.4 kyr for all of southeastern Hudson Bay (Peltier 1998) or 7.6 kyr for Richmond Gulf (Peltier 1994). The reasons for the high variability in the best-fitting exponential form are plentiful: data points are often subjectively included or excluded from the RSL record, gaps in the RSL record occur frequently, methods of decay time estimation are themselves highly variable (compare, for example, Dyke & Peltier 2000; Mitrovica et al. 2000), RSL data at neighbouring sites have been grouped or binned inconsistently (Mitrovica et al. 2000), and the meltwater contribution to the RSL history imposes a dependence on the ice model used. For these reasons, we avoid the use of decay times for this study.

Instead, we apply a method similar to Mitrovica & Peltier (1995)'s NRSL. NRSL involves scaling a modelled RSL curve so that it fits the amplitude of the RSL data, and then computing a misfit between the RSL data and the scaled model results. The overall amplitude of an RSL curve is strongly dependent on the thickness of deglaciated ice in that region, and so scaling the RSL curve removes that dependence. It is then primarily the curvature of the RSL curve that governs the misfit, and the curvature is much more dependent on the viscosity profile than on the ice history. This is the same rationale for using exponential decay times. However, NRSL is less sensitive to many of the difficulties that complicate the use of those decay times.

Our analysis method differs from Mitrovica & Peltier (1995)'s NRSL by imposing a limit on the amount of scaling allowed. Since typical RSL data often have a large time-dependent scatter, allowing an arbitrarily large or small scaling of the model output (that is, scaling of the RSL curves) can lead to reasonable fits to the data by absurdly exotic viscosity structures. For example, an extremely stiff mantle (greater than 10²⁴ Pa s) would lead to a very small amplitude and near-linear RSL history. For some sites with highly variable RSL data, a sufficiently high scaling of this model RSL curve can bring it into reasonable accordance with the data, even though it should be safe to reject models that would require several times thicker ice than is provided by the glaciation model. In this study, we use a normalized RSL approach that allows for scaling of the ice model only within some pre-defined limit, chosen to lie between 10 and 50 per cent. For example, suppose we believe the ice model's thickness is uncertain to 20 per cent. We then scale the model RSL curve by the optimal amount, but keeping the scaling within the range of 0.8–1.2, and then computing the misfit to the RSL data. ‘Optimal,’ in this case, means the scaling that provides the smallest misfit. Thus model data that are scaled to match the RSL amplitude end up fitting only the shape of the RSL curve, but with limits on the scaling so that we reject viscosity profiles that require an unrealistic ice thickness. Note that by scaling the model RSL curves we are accommodating an uncertainty in the ice model without actually changing the ice model, only the model output data. A similar technique is applied to GRACE data, as described below.

RSL data for the six sites are shown in Fig. 2. It is apparent that the scatter in the data is not consistent with the relatively small error bars. This appears to be a feature of most RSL data sets. The lines in Fig. 2 show two typical forward-model RSL curves for each site (these are the models labelled I and II in Section 5). The misfit of a given RSL curve to the data is computed in a manner similar to that used by Mitrovica et al. (2000). For each point in the RSL data, with time t_i and sea level s_i, the nearest point on the model curve is found (T_i, S_i), and that point is given an error of

1

where δs_i and δt_i are the errorbars on data point i, and N is the number of data points in the RSL data at the given site. A single number for the misfit for all RSL data, χ²_RSL, is computed as the average of the χ²_j for each of the six RSL sites around Hudson Bay.

Figure 2.

RSL data at the sites shown in Fig. 1 (points with error bars). The lines show the results of two forward models: model I (solid line) with η_UM = 5.3 × 10²⁰ Pa s and η_LM = 2.3 × 10²¹ Pa s for the upper and lower mantle viscosities, respectively; and model II (dashed line) with η_UM = 4.3 × 10¹⁹ Pa s and η_LM = 6.6 × 10²² Pa s.

Open in new tabDownload slide

2.3 Grace

The GRACE satellite mission, launched in March 2002, recovers global, monthly solutions for the Earth's gravity field down to scales of a few hundred kilometres (Tapley et al. 2004). These solutions can be used to estimate the secular change in the gravity field, which can then be compared with PGR predictions. This constraint is similar to the constraint provided by SLR. However, unlike for , the much shorter spatial resolution available from GRACE makes it possible to separate the secular signal over northern Canada from secular signals elsewhere.

We use fields based on Release 4 solutions from the University of Texas, for 53 months between April 2002 and December 2006. The solutions are in the form of spherical harmonic (Stokes) coefficients, C_lm and S_lm. These fields have been post-processed to reduce noise and artificial vertical stripes in the data, using the method described by Swenson & Wahr (2006). We simultaneously fit a constant, a linear trend, and an annually varying term to the time-series of each coefficient. We transform into the spatial domain to obtain secular changes in gravity on an equal-area grid. The results are smoothed to reduce noise using a Gaussian smoothing function of 400 km half-width (eqs 32–34 Wahr et al. 1998). Fig. 3(a) shows results in the vicinity of Hudson Bay. A large-amplitude anomaly, of about 2.5 μGal yr⁻¹, is spread out over ∼3000 km around Hudson Bay, with a broad maximum just west of Hudson Bay and a secondary maximum southeast of it. We interpret this anomaly as a PGR signal.

Figure 3.

Secular GRACE gravity signal over Hudson Bay (a) and its errors (b), and two forward models using ICE-5G (c) and ICE-3G (d). All plots have units of μGal yr⁻¹. The forward models are produced using the viscosity profiles optimized for ICE-5G (Peltier 2004) and ICE-3G (Tushingham & Peltier 1991), respectively. Hydrological contributions have been removed using output from the GLDAS/Noah land surface model (Rodell et al. 2004). The dashed line in Fig. 3(b) bounds the region used to compare with PGR model predictions; it follows the 1.0 μGal yr⁻¹ contour from Fig. 3(a).

Open in new tabDownload slide

GRACE has no vertical resolution, and so can not distinguish between a PGR signal in the solid Earth and a linear trend in water storage in this region. Since we are using only 4.5 yr of GRACE data, water storage variability at multiyear periods could affect our trend estimates. To reduce this problem we remove the water storage gravity trend for the same 4.5-yr period predicted from the GLDAS/Noah land surface model (Rodell et al. 2004), before constructing the smoothed results shown in Fig. 3(a). The predicted water storage trend in this region is small, only about 10 per cent of the total GRACE trend.

The GRACE secular gravity estimates will be compared below with PGR predictions for many viscosity profiles, to help determine the best-fitting viscosity. To make the model predictions directly comparable to the GRACE estimates, each predicted gravity field is first destriped as described by Swenson & Wahr (2006) and then smoothed with a 400 km Gaussian, just as we do to obtain the GRACE results.

First, though, Figs 3(c) and (d) show the trend in gravity (after destriping and smoothing) predicted using de-glaciation models ICE-5G (Peltier 2004) and ICE-3G (Tushingham & Peltier 1991), together with the same viscosity profiles used in the construction of those models (viscosity model VM2 (Peltier 2004) is used here for the ICE-5G results). There is reasonable agreement between the ICE-5G and GRACE results, though the Gaussian smoothing has sharply reduced the amplitude of the ICE-5G maximum southeast of Hudson Bay. None of the GRACE data were used in the construction of ICE-5G, which makes this agreement particularly satisfying. The ICE-3G results, on the other hand, differ significantly from those of GRACE, particularly in the amplitude of the maximum just west of Hudson Bay. This illustrates one of the most promising future applications of GRACE for PGR studies: helping to constrain the ice model. That application is beyond the scope of this paper, however.

The GRACE-model comparisons described below are done by summing the difference in predicted and observed gravity rates over an equal-area grid of 200 points, distributed within the heavy dashed line shown in Fig. 3(b). That line marks the contour of 1.0 μGal yr⁻¹ in the GRACE trend (Fig. 3a). We use this restricted region to minimize contamination from external gravity trends unrelated to the Canadian PGR signal (e.g. variations in water/snow/ice in Canada and Alaska, and PGR and present-day ice mass signals in Greenland). The restriction to this relatively small region also makes our results less sensitive to errors in the spatial pattern of the deglaciation model. The disadvantage of using a small region, is that we are then also less able to distinguish between different spatial patterns predicted for different viscosity profiles. Instead, our GRACE constraint reduces mostly, though not entirely (see below), to a constraint on the amplitude of the gravity trend.

The GRACE-model comparisons require an estimate of the errors in the GRACE trends. First, we estimate the uncertainty in the monthly values of each individual Stokes coefficient, by removing a constant and an annually varying term from the time-series of each coefficient, and interpreting the residuals as a measure of the error (Wahr et al. 2006). This tends to be an overestimate of the errors, since some of the non-annual variability is certainly a real signal. The monthly uncertainties are then used to determine the uncertainty in the trend of each Stokes coefficient. By assuming the errors in different Stokes coefficients are uncorrelated, we are able to combine the Stokes coefficient uncertainties to obtain an uncertainty in the gravity field trend at each point in the latitude/longitude domain (see Wahr et al. 2006, eq. 4). The resulting spatial error field is shown in Fig. 3(b): a uniform increase in error with decreasing latitude. This pattern is an oversimplification of the true error pattern, due to our assumption that the errors in different Stokes coefficients are uncorrelated. The overall error magnitudes, though, are well represented by the results shown in Fig. 3(b).

Model misfits to the GRACE data are computed as described above for the RSL data (eq. 1): with

2

where σ_j is the error in the GRACE gravity trend at location j, and j = 1, …, M and M = 200. Since the model gravity fields are relatively smooth, as in Fig. 3(c), their average χ² is fairly independent of the number of points used in the summed field.

We have found that the GRACE data calculated from a forward model are quite sensitive to the ice load applied, and that the predominant sensitivity is to the thickness of that load. This issue is discussed further in Sections 4 and 5. To accommodate some level of uncertainty in the loading model, we have used an approach similar to the normalized RSL discussed above. First, a limit on the ice-amplitude uncertainty is chosen; for example, a 20 per cent uncertainty. The regional gravity signal from a given viscosity model is then computed, the magnitude of which is proportional to the ice-amplitude (for constant spatial distribution of ice load). The gravity field is then scaled up or down by the value within the range 0.8–1.2 that minimizes the misfit to the GRACE data. For scaling within this range, therefore, the ice-amplitude dependence of the gravity signal is removed, and the shape of the signal provides the misfit (e.g. the slope of signal reduction in moving away from the point of maximum amplitude). For a modelled gravity field with a maximum value outside the range of 0.8–1.2 times the maximum of the GRACE data, the misfit becomes dominated by this amplitude mismatch, which is appropriate if the ice model is known to be accurate to within 20 per cent.

3 Computational Methods and Models

In our inversion, we use a simple grid search for one- and two-layer inversions and a Monte Carlo method for multiple layer inversion. We calculate a large number of forward models for PGR with different viscosity structures. Calculation of PGR involves solution of the governing equations of mass and momentum conservation, along with gravitational perturbation via Poisson's equation and viscoelastic rheology (Wu & Peltier 1982). Our computational method is identical to that described in Paulson et al. (2005). We assume an incompressible Earth with self-gravitation, and a Maxwell solid mantle overlying an inviscid core. The mantle is assumed to have 1-D structure in its viscosity, density, and elastic properties. The mantle can have multiple layers with different viscosities and shear moduli. However, only three mantle density layers are allowed, with interfaces at 410- and 670-km depths. The model parameters are given in Table 1. The governing equations are solved via Love numbers in the Laplace transform domain using the collocation technique (Mitrovica & Peltier 1992; Wahr et al. 2001). We include the effects of centre-of-mass motion, dynamic ocean response through the Sea Level equation, and the new formulation of True Polar Wander described by Mitrovica et al. (2005).

Table 1.

Earth model parameters used in this study.

Open in new tabDownload slide

3.1 Compressibility test

Since our models assume an incompressible Earth for computational efficiency, it is necessary to examine the potential effects of compressibility on our results. The PGR observables discussed above are computed with both our incompressible formulation and a compressible formulation for two different viscosity models, one uniform and the other with a viscosity contrast factor of 10 at a depth of 670 km. The compressible formulation does not use the collocation method, but explicitly finds a set of modes and their associated residues (Han & Wahr 1995). The rest of the calculation, summing the modes and convolving with the load, is performed as with the incompressible spectral code (Paulson et al. 2005). The compressible calculations use the PREM (Dziewonski & Anderson 1981) density and elastic parameter values. The incompressible calculations use the structure shown in Table 1, which is a three-layer approximation of PREM. Table 2 shows the χ² misfit between compressible and incompressible results, for the RSL and GRACE observables at the same sites and times as the data described in Section 2. The values in Table 2 show the misfit of the compressible results to the incompressible results (i.e. the incompressible results are treated as the ‘data’ and the compressible results are the ‘model’ for these calculations). The misfits reflect only the contribution from the incompressibility assumption made in our forward modelling. These misfits values are small compared to the model/data misfits shown in the following sections. This assures us that the use of incompressible code does not significantly affect the results given below.

Table 2.

χ² misfits between two versions of the spectral code, one compressible and the other incompressible. ‘Viscosity I’ compares models of a homogeneous mantle of 10²¹ Pa s. ‘Viscosity II’ compares models of η_UM = 10²¹ Pa s and η_LM = 10²² Pa s. The incompressible results are treated as the “data” and the compressible results are the “model” for these χ² calculations. All misfits are small compared to those considered in this study.

Open in new tabDownload slide

4 One-Layer Inversion

The first and most straightforward viscosity inversion is for a mantle with homogeneous viscosity. Using synthetic inversion studies, Paulson et al. (2005) found that GRACE and RSL observations over the load regions provide information on the log-average of mantle viscosity below the load. Thus, the results of a one-layer inversion presumably reflect a weighted average of the mantle viscosity profile beneath Hudson Bay. The forward models are run with the density and elastic parameters shown in Table 1, a lithospheric thickness of 120 km, and differing mantle viscosities.

When considering only one layer of viscosity, the GRACE data prefer either of two whole-mantle viscosities: 1.8 × 10²¹ or 2.3 × 10²² Pa s (when no ice-amplitude scaling is applied). This can be seen in Fig. 4(a), which shows two minima that do not change even when allowing the ice load an independence of up to 50 per cent in amplitude, though those minima are far less well defined in the 50 per cent case. Viscosity values between these two minima are rejected mostly due to a discrepancy in the gravity signal magnitude. This can be deduced by noting that the misfit in that range of viscosities decreases when the scaling factor is allowed to increase; that is, when the ice-amplitude is allowed to adjust itself to bring the modelled gravity signal magnitude more in line with the GRACE magnitude. Viscosity values outside of this range, both above it and below it, are ruled out no matter what scaling factor is allowed. They produce a much smoother gravity pattern than is evident in the GRACE results, and no amount of scaling can change that pattern.

Figure 4.

Dependence of χ² misfits for the GRACE data (a) and the RSL data (b), and the total misfit to both the GRACE and RSL data (i.e. average of the two) (c) on whole-mantle viscosity, and dependence of total misfit on lithospheric thickness (d) for homogeneous mantle viscosity models. In Figs 3(a) and (b), the different lines show how the results change when allowing uncertainty in the ice amplitude ranging from none (‘no scaling’) to 50 per cent. In Fig. 3(c), the dotted, dashed and solid lines are the misfits to GRACE, RSL, and both GRACE and RSL data, respectively for an ice uncertainty of 20 per cent. In Fig. 4(d), four different curves show misfits to four different types and combinations of data.

Open in new tabDownload slide

This two-minimum feature is similar to an inversion based on a single number, such as : since the value of is small for weak viscosity (where the relaxation completes shortly after deglaciation), is larger for intermediate viscosities, and is small again for very stiff viscosity (where the relaxation is very slow), the correct value of will usually be fit perfectly at two values in this roughly Gaussian-shaped (small-large-small) dependency. Similarly, we see the GRACE data being optimally fit for two viscosity values, that each give an overall gravity signal magnitude of the same magnitude as GRACE. The two values preferred by GRACE are well constrained, and, as will be shown, the RSL data can discriminate between the two.

Fig. 4(b) shows a similar plot for the RSL data. The RSL misfit in all cases tends to be much larger than the misfit to GRACE data, due predominantly to the quality (scatter) of the RSL data. There is only one misfit minimum for RSL data, occurring at a viscosity of 1.6 × 10²¹ Pa s when the ice-amplitude is taken without error (the curve labelled ‘no scaling’ in Fig. 4b). The bottom of this curve broadens and shifts to 1.0 × 10²¹ Pa s when the ice amplitude is allowed to scale by up to 50 per cent. The minima for these curves are wider than for GRACE data (note the different horizontal scale from Fig. 4), indicating a poorer constraint on mantle viscosity.

We combine these results for a single best measure of homogeneous mantle viscosity in Fig. 4(c), averaging the GRACE and RSL χ² with equal weight. The solid line shows the average of the GRACE and RSL misfit when allowing for a 20 per cent scaling in ice amplitude. It shows a minimum misfit at a viscosity of 1.7 × 10²¹ Pa s. The GRACE data have two nearly equal preferences for mantle viscosity, and the RSL data select one of them. Generally, for a mantle with homogeneous viscosity, the GRACE and RSL data require a viscosity between 1.4 × 10²¹ and 2.3 × 10²¹ Pa s.

We use a 20 per cent scaling limit when generating our preferred viscosity estimates, not only in this one-layer case, but in the two-layer case discussed below. This choice is somewhat arbitrary, but (as can be inferred from Fig. 4 in the one-layer case), it gives about the same results as smaller scaling, and it is more-or-less equivalent to the assumption that the maximum ice thickness in ICE-5G is accurate to ±1 km.

Given the results of studies such as Mitrovica & Forte (1997) and Zhong et al. (2003), we do not expect lithospheric thickness to affect PGR observables near the loading centre. Fig. 4(d) demonstrates this, showing misfits to the data for a lithosphere ranging from 60 to 250 km thickness, overlying a mantle with a uniform viscosity of 1 × 10²¹ Pa s. The figure shows the misfit for GRACE, RSL and ‘total’ (the average of the GRACE and RSL misfits), as well as the lithospheric dependence of ‘coastal RSL’ sites. The latter are RSL data for sites along the North American east coast and arctic coast, on the periphery of the loaded region, and are shown here to reinforce the conclusions of Zhong et al. (2003), that these periphery sites are sensitive to variations in lithospheric thickness. The data used for inversion in this study, GRACE and RSL over the load region, are insensitive to these variations.

5 Two-Layer Inversion

We turn next to the search for a two-layer viscosity model to fit the GRACE and RSL data. The forward models are again run with the parameters given in Table 1 and a lithospheric thickness of 120 km. A test of the inversion is performed with synthetically generated observables equivalent to those used in this study and, as in a previous study (Paulson et al. 2007), is found to unambiguously recover the correct two-layer viscosity structure.

We first consider two-layer inversions with a boundary at a depth of 670 km. We consider misfit measures that accommodate a 20 per cent uncertainty in the ice amplitude, as discussed above. The χ² misfits for GRACE data, RSL data, and a ‘total’ misfit (an average of the two) are shown for a wide range of two-layer viscosity models in Figs 5(a), (b) and (c), respectively.

Figure 5.

Dependence of χ² misfits to the GRACE data (a) and RSL data (b) and total misfit to both the GRACE and RSL data (c) on upper mantle and lower mantle viscosity for two-layer models with layer boundary at 670 km depth, and the total misfit for two-layer models with layer boundary at 1170 km (d). Figures (e) and (f) are equivalent to Figure (c), except that different amount of ice amplitude scaling is allowed: 10 per cent (e) and 40 per cent (f). The color scales give the χ² misfit for each type of data set considered. White regions indicate models with χ² above the values on the scale. In Fig. 5(c), the dotted lines show the locus of models where η_LM = η_UM, η_LM = 10η_UM and η_LM = 100η_UM, as indicated on each. The circles indicates the model of lowest misfit. In Fig. 5(d), the square indicates model VM2-2, with which the ICE-5G model was developed.

Open in new tabDownload slide

The most obvious feature of the GRACE misfits is their ring-shape (Fig. 5a). This shape is the 2-D extension of the two-minimum feature seen in Fig. 4(a)—the latter being the profile along a line from the lower left to the upper right corner of the two-layer misfit plot. The low-misfit ring can be understood as a contour of constant gravity-signal amplitude in the space of two-layer models shown. Models in the white centre of the ring fail to fit the GRACE data primarily because they have too large an amplitude, and models in the white space outside the ring have too small an amplitude. This low-misfit ring becomes thinner, but remains a ring, if we allow smaller ice amplitude scaling than 20 per cent, and becomes thicker for larger allowed scaling.

The RSL misfit plot (Fig. 5b) shows a broad region of low misfit near the uniform viscosity models, and a ‘tail’ of low misfit along models of weaker upper mantle and stronger lower mantle. As in the one-layer inversion, the RSL misfits are all larger than for the GRACE data, and the minimum is broader.

Combining the GRACE and RSL misfits by averaging them together (Fig. 5c), we obtain a measure of ‘total’ misfit to our PGR data. The ring shape from GRACE remains, but the best-fitting models are now restricted to the overlapping regions of small misfit: those models of near-uniform viscosity and in a ‘tail’ of models of weaker upper mantle and stronger lower mantle. The best misfit is at near-uniform mantle viscosity. In the vicinity of the lowest misfit, other models of low misfit extend in a roughly linear region of negative slope on the plot. This indicates a trade-off effect in the two layers: decreased viscosity in one layer can be compensated by an increased viscosity in the other. The best-fitting model appears at η_UM≈ 5.3 × 10²⁰ Pa s and η_LM≈ 2.3 × 10²¹ Pa s.

To examine the effect of a layer boundary at a greater depth, the ‘total’ misfit plot is repeated for a two-layer model with a boundary at a depth of 1170 km. This depth is chosen to coincide with the best two-layer approximation to Peltier's viscosity model VM2 (Peltier 1996), the model that was developed jointly with ICE-5G. This two-layer approximation shall hereafter be referred to as model VM2-2. Fig. 5(d) shows the two-layer inversion with a boundary at that depth. The essential difference between the misfit plots with boundaries at 670 and 1170 km is the roughly vertical and horizontal orientations (respectively) of the low misfit regions: in the vicinity of the best-fitting model, the slope of the low-misfit region is steeper for the shallower boundary (670 km) since a larger change would be required in the thin upper layer to compensate for an equivalent change in the larger lower layer. For a boundary at this depth, the best-fitting model has η_UM≈ 6.7 × 10²⁰ Pa s and η_LM≈ 6.4 × 10²¹ Pa s. This minimum is reasonably close to model VM2-2 (at 0.9 and 3.6 × 10²¹ Pa s, respectively, shown by the small white box in Fig. 5d), which was used during the construction of ICE-5G (Peltier 2004). In fact the misfit for VM2-2 is not notably larger than the minimum misfit. We find that the misfit plots for boundary depths other than the two shown in Figs 5(c) and (d) produce low-misfit regions that vary smoothly between (and beyond) these two examples.

Figs 5(e) and (f) are included to demonstrate the effects of scaling. They are equivalent to Fig. 5(c) except that 10 and 40 per cent uncertainty in the ice amplitude is allowed (respectively). While the best-fitting models are the same as before, greater scaling of the ice amplitude can accommodate many of the models rejected in the case of 20 per cent scaling.

Although for any choice of a boundary depth there is a model with an absolute minimum misfit, there also exists an ensemble of alternative viscosity models that fit the data reasonably well. These models generally have weaker upper mantle and stronger lower-mantle viscosities. To demonstrate this, we show how the data are fit for two models with the layer boundary at 670 km depth: model I produces the overall best-fit, with η_UM = 5.3 × 10²⁰ Pa s and η_LM = 2.3 × 10²¹ Pa s; model II has η_UM = 4.3 × 10¹⁹ Pa s and η_LM = 6.6 × 10²² Pa s. Model II appears in the ‘tail’ of low misfit discussed earlier (Fig. 5c). Both models I and II reproduce the GRACE data reasonably well (Fig. 6a for model I, and Fig. 6b for model II). Figs 6(c) and (d) show the misfits to the GRACE data for models I and II, respectively. Model I has a misfit to the GRACE data of χ² = 0.8, and model II has a misfit of χ² = 1.7—note these are the average values of χ² over the region. Neither model requires an ice-amplitude scaling beyond 20 per cent. The same two models fit the RSL data with misfits of χ² = 5.3 and 6.1, respectively. Fig. 2 shows the RSL fit for these two models.

Figure 6.

Model GRACE results for two-layer model I (a) and model II (b) in μGal yr⁻¹, and the misfit to GRACE data of model I (c) and II (d). Model I has η_UM = 5.3 × 10²⁰ Pa s and η_LM = 2.3 × 10²¹ Pa s. Model II has η_UM = 4.3 × 10¹⁹ Pa s and η_LM = 6.6 × 10²² Pa s. Both models have a layer boundary at 670 km depth. Model RSL for these two models at Melville Peninsula are shown in Fig. 1. The dotted lines in parts c and d show the region over which χ² is computed (as in Fig. 3b). Contours above χ² = 7 have been omitted for clarity, as they lie outside the region of interest.

Open in new tabDownload slide

The two-layer inversion results can be reported in a form amenable to use with an internal loading gravity study. In such a study, long-wavelength geoid anomalies are predicted using a convection model and density estimates inferred from seismic tomography, and provide an inference of the viscosity contrast across a layer boundary in the mantle (Hager & Richards 1989). Since such a study provides the ratio of layer viscosities, it could resolve the ambiguity remaining in the two-layer inversion. For example, if one adopts a factor of 10 viscosity contrast across the 670-km boundary, then the allowed viscosity models would lie along the line log₁₀η_LM = log₁₀η_UM+ 1 in Fig. 5(c) (labelled ‘10’). Fig. 7 shows the combined GRACE and RSL misfit along this profile, as well as along profiles for other viscosity contrasts. Comparing these curves, the minimum misfit occurs for a viscosity contrast of 3, but minimum misfits for other contrasts are not much greater. The figure legend shows the upper- and lower-mantle viscosities that best fit each contrast factor. For example, the factor of 30 viscosity contrast suggested by Hager & Richards (1989) requires η_UM = 1.4 × 10²⁰ Pa s and η_LM = 4.3 × 10²¹ Pa s to fit the GRACE and RSL data.

Figure 7.

The total misfits versus lower-mantle viscosity for a given viscosity contrast between the upper and lower mantles for two-layer mantle with a layer boundary at 670 km depth. The curves represent lines of unit slope through Fig. 5(c). The legend gives the viscosity contrast, and the upper- and lower-mantle viscosities of minimum misfit for that contrast in units of 10²¹ Pa s. For example, if we assume η_LM = 10η_UM (the × symbol), the best-fitting two-layer model is at η_UM = 0.28 × 10²¹ Pa s and η_LM = 2.85 × 10²¹ Pa s.

Open in new tabDownload slide

As a final note to the two-layer inversion, we produced (but do not reproduce here) misfit plots equivalent to those shown in Fig. 5(c), but with the two viscosity layers residing entirely above 1800 km depth. If the viscosity below 1800 km depth is given the constant value of 10²² Pa s or greater, and the other two layers are allowed to vary as before, we find that the models of low misfit are not appreciably different from the results presented here. That is, the two-layer inversion is providing information only about the viscosity structure above 1800 km. A similar conclusion for these types of PGR observables appears in the inversions of Mitrovica & Peltier (1992) and Mitrovica & Forte (2004).

6 Three-Layer Inversion

We attempt a final inversion for three viscosity layers. Due to the reduced sensitivity to very deep viscosity, we choose the three layers to all lie above 1800 km depth and to have equal thickness. This provides layer boundaries close to both 670 and 1200 km depth. The viscosity below 1800 km is set to 10²² Pa s. As expected from the synthetic studies of Paulson et al. (2007), the viscosities in these three layers are able to vary greatly under a ‘trade-off’ effect, where viscosity in one layer may be increased provided that the viscosity of a neighbouring layer is simultaneously reduced, all while preserving the fit to the PGR data.

The results are shown in Fig. 8(a). Every point represents a forward run in the inversion, with its χ² misfit to the PGR data on the horizontal axis, and with the vertical axis showing a measure of the average viscosity difference between that model and the model of lowest χ² misfit. The average viscosity difference between two models η₁(r) and η₂(r), labelled Δ, is computed by

3

where the integral is over the whole mantle. Thus Δ = 1 implies that the two viscosity models differ, on average, by an order of magnitude. We see that the best-fitting models, those on the leftmost edge, vary among themselves by nearly an order of magnitude in their averaged viscosity. We choose three models from the lowest misfit (left edge) of the plot to display in Fig. 8(b). These three models have very different but correlated viscosity structures, exhibiting the trade-off of viscosity in neighbouring layers.

Figure 8.

The viscosity difference Δ versus misfit for the best-fitting three-layer models (a) and viscosity versus depth for three models with similarly small misfits (b). In Fig. 8(a), calculation of Δ is described in the text, and only the viscosity structure above 1800 km depth is included. Among the models of the lowest misfit (the circles on the left edge of the scatter plot) are a wide range of viscosity structures, implying very poorly constrained viscosity for any given layer. Three models, labelled a, b and c, are shown in Fig. 8(b). In Fig. 8(b), although there are significant differences in their viscosity structure, the alternating viscosity values in neighbouring layers indicates that it is the ‘trade-off effect’ that prevents a well-constrained three-layer inversion.

Open in new tabDownload slide

We conclude, as in Paulson et al. (2007), that these PGR data cannot constrain more than two layers of viscosity. Including multiple layers in the inversion allows large variability in viscosity structure, as any layer may compensate for the high viscosity of its neighbour with a smaller value of its own.

7 Conclusions and Discussions

The main results of this paper can be summarized as follows:

• 1

  The secular component of the GRACE data reveals a positive change in gravity with a maximum of 2.5 μGal yr⁻¹ over a broad region (>3000 km in diameter) near Hudson Bay. The pattern and amplitude of the gravity change are in good agreement with PGR predictions from the ICE-5G ice model and viscosity profile, suggesting a PGR origin for this gravity change. The GRACE—ICE-5G agreement is significantly better than that between GRACE and predictions from ICE-3G, indicating that ICE-5G is a significant improvement over ICE-3G in representing the overall deglaciation history of the Laurentide ice sheet. The good agreement with ICE-5G is especially meaningful, because no GRACE data were used in the construction of that ice model.

• 2

  For a mantle with homogeneous viscosity, the GRACE and RSL data in Hudson Bay region require a viscosity between 1.4 × 10²¹ and 2.3 × 10²¹ Pa s. An inversion for a two-layer mantle viscosity structure shows that with the upper/lower mantle boundary at 670 km depth, η_UM = 5.3 × 10²⁰ Pa s and η_LM = 2.3 × 10²¹ Pa s yields the best fit to the GRACE and RSL data. However, a weaker upper mantle together with a stronger lower mantle may also provide a reasonable fit to the data.

• 3

  The GRACE and RSL data are insensitive to mantle viscosity below 1800 km depth, similar to that from previous studies with only RSL data. However, the current data used in this study cannot provide meaningful constraints on the top 1800 km of the mantle with three viscosity layers due to trade-off effects.

While the GRACE and RSL data are sensitive to the loading model used, we can remove some of this dependence by scaling the data to allow for a possible error in the ice amplitude. This approach is analogous to the ‘NRSLs’, shown by Mitrovica & Peltier (1995) to be more sensitive to viscosity structure than to the ice amplitude. We apply the same methods in the use of RSL data around Hudson Bay. When seeking the viscosity model that best fits the GRACE and RSL data, we allowed for a 20 per cent uncertainty in the magnitude of the ice load.

The limited resolving power of the GRACE and RSL data leaves some ambiguity in the inverted two-layer mantle viscosity structure. Two ways to remove the ambiguity are: (1) adding additional PGR-related observations, and (2) including observations of non-PGR processes in the inversion. Observations of the long-wavelength geoid, for example, help constrain radial viscosity variations (Hager & Richards 1989) and have already been included with PGR observables in joint inversion studies (e.g. Forte & Mitrovica 1996; Mitrovica & Forte 2004). However, a Monte Carlo study such as that employed by Paulson et al. (2007) for PGR is useful to better understand the radial resolving power of the geoid.

In the past, polar wander, RSL data at far-field sites, and SLR observations of have been used to infer 1-D mantle viscosity (e.g. Peltier 1998; Lambeck et al. 1990). Both and polar wander are sensitive to long-wavelength PGR signals, and so can help constrain lower mantle viscosity. However, their use can also introduce additional ambiguities, since they can be affected by processes that are unrelated to PGR and not well understood.

There is also ambiguity in using far-field RSL data, including those from sites close to the peripheral bulge or further away from the ice load region. This ambiguity results from the significant sensitivity of far-field RSL signals to mantle viscosity structures below both the observation sites and the ice load regions Paulson et al. (2005, 2007). The viscosity structures in those locations could be significantly different from one another for a realistic 3-D Earth (e.g. oceanic mantle versus continental mantle). This difficulty is not limited to far-field RSL data, but is relevant to any long-wavelength data, including and polar wander. This is because these signals sample the mantle globally, not just those regions below the ice-loaded region that are sampled by the local GRACE and RSL data (Paulson et al. 2007).

Ongoing vertical and horizontal motion in the ice load regions determined using GPS or absolute gravity measurements, can provide important constraints on mantle viscosity (e.g. Larson & van Dam 2000). More studies are needed to examine the sensitivity and resolving power of such PGR observations for mantle viscosity.

Acknowledgments

This research is supported by NSF grants EAR-0087567 and EAR-0134939, and by a grant from the David and Lucile Packard Foundation. We thank Art Dyke for providing an extensive relative sea level database. We also thank reviewers Georg Kaufmann and Konstantin Latychev for their constructive comments.

References

Author notes