The sensitivity of glacial isostatic adjustment predictions to a low-viscosity layer at the base of the upper mantle

Earth and Planetary Science Letters 154 Ž1998. 265–278

The sensitivity of glacial isostatic adjustment predictions to a low-viscosity layer at the base of the upper mantle Glenn A. Milne a

a,)

, Jerry X. Mitrovica a , Alessandro M. Forte

b

Department of Physics, UniÕersity of Toronto, 60 St. George Street, Toronto, Ont. M5S 1A7, Canada b Institut de Physique du Globe de Paris, Departement de Sismologie, Paris, France ´ Received 10 March 1997; revised 10 October 1997; accepted 10 October 1997

Abstract Recent inferences of mantle viscosity based on convection-related observables Žlong-wavelength nonhydrostatic geoid and free-air gravity harmonics. indicate the presence of a relatively thin low-viscosity region at the base of the upper mantle. We explore the sensitivity of observables associated with glacial isostatic adjustment ŽGIA. to the viscosity and thickness of this region. These observables include post-glacial relative sea-level ŽRSL. variations within previously glaciated areas and in the ‘‘far field’’ of the Late Pleistocene ice sheets, as well as anomalies associated with the rotational state of the planet. The latter include present-day secular variations in the geopotential harmonic J2 Žor, equivalently, variations in the rotation rate. and polar wander speed. We find, in contrast to previous suggestions, that the GIA observables Žwith the exception of polar wander speed. are sensitive to fine-scale low-viscosity structure at the base of the upper mantle. For example, J˙2 predictions can be altered by as much as 2 = 10y11 yry1 by the inclusion of a low-viscosity zone. This signal is sufficient to significantly affect both inferences of mantle viscosity and constraints on the present-day melting of large polar ice sheets ŽAntarctic, Greenland. that are derived using this observable. Indeed, if a low-viscosity region exists, then models that do not include it can overestimate recent polar melting events by as much as 0.5 mmryr of equivalent eustatic sea-level rise. q 1998 Elsevier Science B.V. Keywords: mantle; isostasy; glacial rebound; viscosity; secular variations

1. Introduction Recent viscous flow modeling of the low-degree nonhydrostatic geoid indicates that these convection-related data are highly sensitive to radial viscosity variations at depths in the vicinity of the transition zone Že.g., w1–3x.. This sensitivity is a consequence of the high correlation between the seismically inferred 3-D density structure in this )

Corresponding author. E-mail: [email protected]. utoronto.ca

region and the low-degree nonhydrostatic geoid w1x. A number of authors have used this sensitivity to argue for a thin low-viscosity layer at the base of the upper mantle Že.g., w1,2,4x.. These studies have generally assumed a whole-mantle convection scenario and, in this context, there are several physical mechanisms that may explain the existence of a thin low-viscosity layer. These involve, for example, transformational superplasticity w3,5,6x or the latentheat release associated with upwelling material passing through the 670 km phase boundary w2x. If, on the other hand, mantle material convects as a two-

0012-821Xr98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 0 1 2 - 8 2 1 X Ž 9 7 . 0 0 1 9 1 - X

266

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

layer system with separate reservoirs in the upper and lower mantle, then a thermal boundary layer will develop at the 670 km discontinuity producing high temperature gradients and thus large variations in viscosity over short radial length scales Že.g., w2,7x.. It is notable, however, that the preference for a low-viscosity layer at 670 km depth is weaker in the case of geoid modeling which assumes a layered mantle system Že.g., w8x.. Observables used to infer viscosity structure within the Earth also include those associated with glacial isostatic adjustment ŽGIA. w9–17x. The response of the Earth to the Late Pleistocene glacial cycles is manifest in a number of geophysical phenomena, such as relative sea-level ŽRSL. change, secular variation in the Earth’s rotation vector, anomalies in the gravitational field and 3-D motions of the Earth’s crust. Attempts to infer viscosity on the basis of these observables have invariably been based on forward analyses characterized by a restricted class of viscosity models. Indeed, by far the most common choice is a three-layer model which includes a high-viscosity Želastic. lithosphere, and isoviscous upper and lower-mantle regions Žwith the boundary between the two regions taken at 670 km depth.. Almost no effort has been made to consider the sensitivity of the observables to fine-scale viscosity structure at depth. One of the few exceptions is the study of Peltier w7x, who used an observational constraint on the non-tidal acceleration of the Earth’s rotation rate Žor, effectively, the present-day secular variation in the degree-two zonal harmonic of the Earth’s geopotential, J˙2 . to infer the average viscosity of the lower mantle. In particular, Peltier w7x investigated the influence of a thin Ž50 km. high-viscosity layer located immediately beneath the 670 km boundary on J˙2 predictions; he argued that this ‘‘internal lithosphere’’ would be the net result of a thermal boundary layer at the 670 km discontinuity. ŽThe analyses w2,4x have, in contrast to w7x, argued that a thermal boundary layer at 670 km depth would also give rise to a pronounced low-viscosity region at this depth.. The results of Peltier’s sensitivity analysis show that J˙2 predictions are relatively insensitive to this type of fine-scale viscosity structure, with a maximum effect of ; 15% for Earth models in which the lower-mantle viscosity is ; 2–3 times greater than an assumed

upper-mantle viscosity of 10 21 Pa s. In general, the influence of the high-viscosity layer does not exceed a few percent. Others w18x considered a similar class of models and obtained a consistent result. To our knowledge, a sensitivity analysis similar to that appearing in w7x has not been carried out for a low-viscosity layer at the base of the upper mantle. Recent results from long-wavelength geoid modeling described above suggest that such an analysis is warranted. There have been suggestions that GIA observables are relatively insensitive to such a layer Že.g., w19,20x. although no quantitative support for this assertion has appeared. Indeed, this assumption of insensitivity plays an important role in a recent analysis by Peltier w20x which considered an inversion of GIA data Žincluding RSL variations, J˙2 and polar wander speed.. Peltier w20x claimed that his inverted viscosity profile was compatible with an independent profile derived from geoid modeling despite the absence of a low-viscosity notch in the former. He argued that the addition of a low-viscosity ‘‘notch’’, which would be necessary to fit the geoid data, would have no effect on the GIA observables used in the inversion. Recent simultaneous joint inversions of GIA and convection-related observables w21,22x are characterized by a low-viscosity layer immediately above the 670 km boundary; however, it is not clear whether the GIA data played a role in constraining this feature. Our goal here is to explicitly examine the sensitivity of a suite of GIA observables to the viscosity and thickness of a lowviscosity layer in the transition zone region. We choose to consider three GIA observables: RSL variations, J˙2 and polar wander speed. RSL data are the most common GIA observables used to infer mantle viscosity Že.g., w13,16,17,23– 26x.. Predictions of these observables are, however, sensitive to details of the space–time geometry of the now disintegrated ice sheets. A number of authors have reduced this ambiguity by adopting carefully chosen parameterizations of RSL data from specific geographic areas. For example, Nakada and Lambeck w16x considered the height difference between Late Holocene sea-level highstands at nearby sites distant from the Late Pleistocene ice sheets. ŽIn the following we adopt the term ‘‘far field’’ for these sites. We refer the reader to Ref. w16x for a discussion of this term.. Furthermore, Mitrovica and

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

Peltier w26x estimated site-dependent decay times using RSL data obtained from regions near the center of previously ice-covered areas. In considering the sensitivity of RSL variations to a low-viscosity zone at 670 km depth we will focus on these parameterizations. Predictions of the GIA component of J˙2 and polar wander speed have been shown to be sensitive to lower-mantle viscosity and relatively insensitive to upper-mantle viscosity w27x. These data have, consequently, been used to infer deep mantle viscosity w7,15,28,29x with some models even including DY structure in the analysis w18x. The inference of viscosity is, however, complicated by the fact that these observables are also sensitive to, for example, present-day mass flux associated with ice sheets and mountain glaciers Že.g., w27,30,31x.. Clearly, the J˙2 and polar wander observables cannot be used to constrain mantle viscosity until these contributions are independently estimated. Alternatively, these observables might be applied to constrain this presentday mass balance once the viscosity is rigorously inferred from independent data sets Žsee, e.g., w22x.. Regardless of the application, it is important to determine whether GIA predictions of J˙2 and polar wander are sensitive to the presence of a thin lowviscosity region at the base of the upper mantle. Our goal is to quantify this sensitivity and compare it with the predicted range of signal expected from other geophysical processes. The latter is now understood to include not only present-day cryospheric mass changes, but also tectonic processes w32x, advection of mantle density heterogeneities w33x and pressure changes in the fluid core w34x.

2. Results and discussion As described in Section 1, a number of recent modeling studies of nonhydrostatic geoid data indicate the presence of a low-viscosity zone above the 670 km discontinuity. The inferred viscosity profile from one such study w2x, in which both geoid and free-air gravity data were considered, is shown as profile ‘‘thn’’ in Fig. 1 Žsolid line.. Note that these data are sensitive only to relative variations of viscosity. Hence, the right-hand axis of Fig. 1 Žwhich represents the logarithm of the viscosity variation.

267

Fig. 1. A subset of the viscosity profiles considered in this study. Profile ‘‘thn’’ Ž solid line . corresponds to the relative viscosity model proposed by w2x on the basis of nonhydrostatic geoid data, and it is thus arbitrarily scaled Žsee right-hand axis .. The left-hand axis shows the absolute scaling associated with this profile which w22x found provides a best fit to a suite of post-glacial RSL decay times. The thickness of the low-viscosity notch is 70 km. Profiles ‘‘thk’’ and ‘‘nlvl’’ are derived from profile ‘‘thn’’ by either increasing the thickness of the low-viscosity region Žto 220 km. or increasing the viscosity of the region to the value which characterizes the average upper-mantle viscosity. The model ‘‘thk8’’ is derived by varying the viscosity at the base of the upper mantle in model ‘‘thk’’ until a best fit to convection-related observables is obtained Žsee Fig. 2..

can be arbitrarily normalized. In this respect, we set the logarithm of the viscosity of the layer above the low-viscosity notch to be zero. The sensitivity of both the nonhydrostatic geoid data and the free-air gravity data to perturbations in the viscosity of the notch in profile ‘‘thn’’ is explored in Fig. 2. In particular, the figure shows the variance reduction associated with profile ‘‘thn’’ and with a suite of models in which the viscosity of the notch is increased over an order of magnitude from the value which characterizes this profile. Note that the optimum viscosity for this thickness of notch is Žnot surprisingly. the value obtained by Forte et al. w2x, which is ; 100 times lower than the viscosity of the adjacent upper-mantle material Žsee Fig. 1.. In this case, the variance reduction for both the geoid and free-air gravity harmonics is roughly 80%. It is clear that the predictions are sensitive to the viscosity of this thin layer, with the goodness of fit deterio-

268

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

Fig. 2. The variance reduction associated with model fits to the nonhydrostatic geoid Ž solid lines . and free-air gravity Ž dashed lines . harmonics for degrees 2 to 8 as a function of the perturbation to the logarithm of the viscosity of the notch in profile ‘‘thn’’ Ž circles . and profile ‘‘thk’’ Ž triangles .. The left-most abscissa value, ‘‘0.0’’, refers to either model ‘‘thn’’ Ž circles . or ‘‘thk’’ Ž triangles . in Fig. 1. The right-most value, ‘‘2.0’’, represents model ‘‘nlvl’’ in Fig. 1. Finally, the abscissa value 0.6, for the case of the curves denoted by triangles, is model ‘‘thk8’’.

rating rapidly as the notch viscosity is increased. We conclude that the existence of a significant low-amplitude notch is a robust requirement of the convection modeling when the viscosity structure in w2x is adopted. In calculations described below we adopt profile ‘‘thn’’ Žand variations on it. to investigate the sensitivity of a subset of GIA observables to the viscosity of the notch. Since GIA predictions are sensitive to the absolute value of viscosity within the Earth Žas well as depth variations of this parameter., we adopt the scaling determined by w22x Žsee left-hand axis in Fig. 1.. This scaling, which provides an optimal fit to a large set of decay times determined from RSL variations in Canada and Fennoscandia, yields a deep mantle viscosity of ; 10 22 Pa s. Although the geoid data seem to prefer a thin Ž; 100 km. low-viscosity layer Že.g., w4x., these data do not uniquely constrain the thickness or viscosity of this structure Že.g., w35x.. Accordingly, we also consider GIA predictions based on models in which the thickness of the low-viscosity layer is increased

to ; 220 km Žor approximately the depth extent of the transition zone — see the ‘‘thk’’ class of models in Fig. 1.. A low-viscosity layer of this approximate thickness has been proposed in other inversion studies Že.g., w1,21,22x.. The sensitivity of nonhydrostatic geoid data and free-air gravity data to the viscosity of this thicker notch is shown in Fig. 2. Note that the variance reduction is a maximum when the change in the logarithm of the viscosity of this notch is ; 0.60 Žapproximately corresponding to a four-fold increase. relative to its value in profile ‘‘thn’’ or ‘‘thk8’’ Žwe denote this optimal model as ‘‘thk8’’ — see Fig. 1.. This result indicates that there is a trade-off between notch thickness and viscosity Žcompare models ‘‘thn’’ and ‘‘thk8’’.. For purposes of comparison, we will also consider a profile in which the notch viscosity is increased to the value associated with the main upper-mantle layer in Fig. 1. This profile, which we denote as model ‘‘nlvl’’ Žan abbreviation for ‘‘no low-viscosity layer’’., is thus characterized by an isoviscous upper mantle Ži.e., the low-viscosity feature has been removed. below the high-viscosity lithosphere — see Fig. 1. It is clear from Fig. 2 that model ‘‘nlvl’’ Žcorresponding to d log nnotch s 2.0. does an extremely poor job in fitting the convection-related observables. Variance reductions for this model are roughly y2.44 for the nonhydrostatic geoid and y1.44 for the free-air gravity harmonics. In the following, our predictions of GIA observables are generated using spherically symmetric, self-gravitating, compressible, Maxwell viscoelastic Earth models. The elastic structure of these models is defined to be that of the seismic model PREM w36x. The adopted ice load is based on the ICE-3G deglaciation model described in w37x. We include a glaciation phase by simply reversing the sign of the ICE-3G deglaciation increments and spacing these in time so as to produce a glaciation phase of duration ; 100 kyr. Since polar wander speed predictions are sensitive to the number of glacial cycles before present Že.g., w15x., our predictions of GIA induced anomalies in the Earth’s rotation vector were calculated using seven of these full glacial cycles. The complementary ocean load is modeled as having a eustatic Ži.e., geography independent. variation. The viscosity profiles shown in Fig. 1 do not include an elastic lithosphere Žalthough the viscosity of the top

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

80 km layer is more than an order of magnitude greater than the viscosity of the underlying upper mantle.. A thin elastic outer layer is strongly preferred on the basis of independent GIA constraints. Accordingly, our GIA predictions will adopt an elastic lithosphere of 80 km depth Žin practice this rheology is obtained by increasing the viscosity of the top layer in Fig. 1 to very high values.. 2.1. RelatiÕe sea leÕel: PreÕiously glaciated regions Fig. 3 investigates the sensitivity of predicted post-glacial RSL variations at four sites within previously glaciated regions to variations in viscosity within the transition zone. Oslo and Angerman River are situated in the region once covered by the Fennoscandian ice sheet, while Richmond Gulf and Ipik Bay were covered, respectively, by the Laurentide and the Arctic ice complexes. Predictions are shown for profiles ‘‘thn’’, ‘‘thk8’’ and ‘‘nlvl’’. As discussed in Section 1, previous studies have suggested that the decay times associated with the most recent portion of these RSL variations provide a robust constraint on viscosity w26x. The choice of time window is based on two requirements: Ž1. that the solid surface be in a state of free decay; and Ž2. that global eustatic sea-level changes do not contaminate the local decay time estimates. These requirements are met by considering only the last 6.5 kyr in Hudson Bay and Arctic Canada and the last 5.0 kyr in Fennoscandia w22x. Over this time window the site-dependent decay times, t i , are estimated via a Monte Carlo determination of the best-fitting exponential form: RSL i Ž t . s A i  exp Ž trt i . y 1 4

Ž 1.

through the predicted RSL time series RSL i Ž t .. The decay times determined from models ‘‘thn’’, ‘‘thk8’’ and ‘‘nlvl’’ are listed on each frame of Fig. 3. The 2 s observational uncertainty Ž Dtob . is also listed on each frame Žthese uncertainties are taken from w22x.. The two sites from each geographic region ŽHudson Bay and Fennoscandia. were chosen because they exhibit roughly the minimum and maximum effect of the low-viscosity notches on RSL predictions in these areas. The results shown in Fig. 3 indicate that the predicted RSL curves for three of the four chosen

269

sites are significantly influenced by the presence of the low-viscosity layer in profiles ‘‘thn’’ and ‘‘thk8’’. The maximum effect on RSL ranges from ; 12 m at the Ipik Bay site to ; 50 m at the Angerman River site over the time window considered in Fig. 3. This sensitivity is, in general, larger than the observational uncertainties at these sites Že.g., the error in RSL height for data obtained at Richmond Gulf is ; " 2.5 m for index points dating back ; 6 kyr w38x.. Therefore, inferences of viscosity based on raw RSL records from these regions may be biased if the class of viscosity models considered does not include fine-scale structure in the transition zone. Data from previously glaciated regions are commonly used to constrain Late Pleistocene ice thicknesses Že.g., w37x. and these analyses generally assume a simple threelayer viscosity profile Žlithosphere–upper mantle– lower mantle.. In the case of centrally located sites Že.g., Angerman River, Richmond Gulf. we conclude that the presence of a thin low-viscosity notch directly above 670 km depth would lead to significantly higher estimates of ice sheet thicknesses in these regions. That is, if a simple three-layer model fits the data, then the introduction of a low-viscosity notch would have to be compensated by a large increase in local ice thickness. Inferences of viscosity based on the RSL decay times t i ŽEq. Ž1.. are characterized by less sensitivity to the ice load history than the raw RSL time series. ŽThe amplitudes A i , in contrast, are predominantly governed by the load history.. It has been shown Žw22x, fig. 5., via calculations of Frechet kernels, that decay time data become progressively more sensitive to viscosity variations at greater depth as one considers sites closer to the edge of the ancient ice complex at its glacial maximum. This is consistent with the results in Fig. 3, which indicate that the decay times for central sites ŽAngerman River, Richmond Gulf. show more sensitivity to the existence of the low-viscosity notch than the decay times from Oslo and Ipik Bay, which are closer to the edge of their respective ice complexes. As one would expect, the predicted decay times are reduced by the presence of a low-viscosity layer relative to predictions based on profile ‘‘nlvl’’. As discussed by a number of authors Že.g., w39x., GIA scales naturally with the logarithm of the viscosity. Thus, a measure of the ‘‘strength’’ of the notch is

270

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

Fig. 3. Predicted RSL curves at four different sites Žas labeled.. The dotted, solid and dashed lines refer to predictions calculated using the viscosity profiles ‘‘nlvl’’, ‘‘thn’’ and ‘‘thk8’’, respectively, in Fig. 1 Žthese line types correspond to those employed in Fig. 1 to specify the various viscosity models.. The parameters tnlvl , t thn and t thk8 are the decay times ŽEq. Ž1.. corresponding to each of these RSL curves. Dtob is the 2 s observational uncertainty in the decay time for each site.

the logarithm of the viscosity drop times the thickness of the region. Using this measure, the notch in profile ‘‘thk8’’ has approximately twice the ‘‘strength’’ of the notch which characterizes profile

‘‘thn’’. This is roughly in accord with the predictions in Fig. 3, where the difference Žlog t thn y log tnlvl . is on average ; 50% of the difference Žlog t thk8 y log tnlvl .. The relative effects vary from site to site

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

since the detailed depth-dependent sensitivity of decay time predictions to variations in viscosity is a function of the geographic location. To determine whether the variation of the predicted decay times in Fig. 3 is significant, it is necessary to compare the range in predicted values to the observational uncertainty at the four sites. ŽThe ‘‘edge’’ sites Oslo and Ipik Bay have relatively higher observational uncertainties for two reasons. First, the amplitude of the raw RSL variations at these sites is smaller than that of the central sites and the decay time estimates are thus less precise. Second, the observational uncertainties quoted in Fig. 3 are augmented to include a contribution associated with ice load uncertainties. This contribution generally increases as one consider sites closer to the edge of the ancient ice complexes.. The difference tnlvl y t thk8 for the sites Angerman River and Richmond Gulf is larger than the 2 s error. The predicted difference tnlvl y t thn exceeds the 2 s error only in the case of the Richmond Gulf site. Thus, the decay time predictions at central sites are sensitive to the presence of the thick and, to a lesser extent, the thin notch. 2.2. RelatiÕe sea leÕel: The far field Numerical predictions of GIA-induced RSL change in the far field of the Late Pleistocene ice sheets are characterized by a moderate Žsub-mmryr. sea-level fall subsequent to the end of the final deglaciation event. This fall leads to the occurrence of a sea-level highstand with an amplitude of a few meters around 4–5 kyr BP. In Table 1 we show the predicted highstands at a number of sites in the Australian region and its

Table 1 Sea-level highstands Žm. Viscosity model

Port Pirie, Australia Cape Spencer, Australia Halifax Bay, Australia Karumba, Australia Moruya, Australia Christchurch, New Zealand

thn

thk8

nlvl

4.9 2.8 3.3 4.3 3.3 2.9

4.5 2.9 3.1 3.9 2.9 2.1

5.4 3.0 3.9 5.1 3.6 3.0

271

vicinity. The predictions indicate that the calculated highstand is lowered at all sites by the inclusion of a low-viscosity notch. ŽOnce again, the effect of the thicker low-viscosity notch, relative to predictions using profile ‘‘nlvl’’, is roughly twice, on average, the effect associated with the thin notch.. The sensitivity varies from site to site, with the largest predicted variation of ; 1.2 m at Karumba, and the smallest of ; 0.1 m at Cape Spencer. There are three mechanisms that affect the value of the predicted highstands: ocean syphoning caused by the collapse of peripheral bulges w40x, continental levering associated with local ocean loading Že.g., w16,41x. and ongoing melting of global ice from the time of the highstand to present w16x. The first two of these mechanisms are dependent on viscosity structure and thus contribute to the sensitivity observed in Table 1. Since the peripheral bulge collapse produces a more geographically uniform sea-level fall than the ocean loading effect within the far-field regions, the latter likely accounts for most of the site-dependent variation evident in Table 1. Nakada and Lambeck w16x investigated the sensitivity of highstand amplitudes to variations in several model parameters Že.g., upper- and lower-mantle viscosity, lithospheric thickness.. The significance of the results in Table 1 is best illustrated by comparison to these earlier predictions. As an example, the analysis Žw16x, fig. 7. found that a two order of magnitude increase in lower-mantle viscosity produced less than a 2 m change in the predicted highstand at Halifax Bay. The introduction of a low-viscosity notch has an effect which reaches ; 30–40% of this value. Comparable Žsignificant. percentages are obtained for the other sites listed in Table 1. Nakada and Lambeck w16x used differential sealevel highstands Žrather than the observed individual highstand amplitudes. to infer mantle viscosity structure. Predicted differential sea-level highstands for various pairs of sites are listed in Table 2. These predictions are less sensitive to the presence of a low-viscosity notch than the predictions of the individual highstands. Indeed, comparing these values to the relevant observational uncertainties Že.g., w16x, fig. 11. leads us to conclude that previous viscosity inferences based on these differential highstands w16x would not be significantly affected by the existence

272

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

Table 2 Differential highstands Žm. Viscosity model

Karumba–Halifax Bay Port Pirie–Cape Spencer Halifax Bay–Moruya Moruya–Christchurch

thn

thk8

nlvl

1.0 2.1 0.0 0.4

0.8 1.6 0.1 0.8

1.2 2.4 0.3 0.6

certainly coincidental, that the minimum J˙2 amplitude predicted for the thick notch calculations Ž d log nnotch s 0.6. is obtained for the model which simultaneously maximizes the fit to the convectionrelated observables for this class of viscosity models Žthat is, model thk8; compare Fig. 2 — dashed lines with Fig. 4A — solid triangles.. Note that the two curves in Fig. 4A converge, as they must, when the viscosity of both notches reach that of the overlying upper mantle.

of the fine-scale viscosity structure in profiles ‘‘thn’’ and ‘‘thk8’’. The optimum viscosity model based on differential highstand data was used in w16x to predict the highstands at a number of geographic sites. The discrepancy between the predicted highstand amplitudes and the observed values was then used to estimate ongoing melting of the Antarctic ice sheet from the timing of the highstands to present. Since the highstand predictions are sensitive to a lowviscosity layer ŽTable 1., estimates of this melt signal will vary, depending on whether or not such a layer is included in the viscosity model parameterization. For example, on the basis of the results in Table 1, failure to account for the existence of a lowviscosity region would result in an over estimate of the Late Holocene melting from Antarctica by as much as 0.2 mmryr of equivalent eustatic sea-level rise Žor ; 1 m in 5 kyr.. 2.3. J˙2 and polar wander speed Next, we turn our attention to the sensitivity of predictions of present-day J˙2 and polar wander speed to the existence of a low-viscosity notch at the base of the upper mantle. In Fig. 4A we plot predictions of J˙2 as we increase the viscosity of the notch two orders of magnitude from the values which characterize profiles ‘‘thn’’ and ‘‘thk’’ to the case of the isoviscous upper-mantle profile ‘‘nlvl’’. The sensitivity of the predictions to the viscosity within either the thin or thick notch is dramatic. Relative to the results for profile ‘‘nlvl’’, the presence of a thin low-viscosity feature favored by modeling of convection-related observables Žprofile ‘‘thn’’ — see Fig. 2. reduces the J˙2 prediction by 40%. The analogous results for the thick notch show a maximum reduction of ; 30%. It is interesting, though

Fig. 4. ŽA. Predictions of the GIA-induced J˙2 signal as a function of the perturbation to the logarithm of the viscosity of the notch in profile ‘‘thn’’ Ž circles . and profile ‘‘thk’’ Ž triangles .. ŽB. As in frame ŽA., except for predictions of polar wander speed.

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

Fig. 4B shows analogous predictions of polar wander speed as a function of Žthin and thick. notch viscosity. In contrast to the J˙2 results, predictions of polar wander speed are relatively insensitive to the existence of a low-viscosity layer at the bottom of the upper mantle. Indeed, the prediction for model ‘‘thn’’ differs by only ; 9% from the result obtained using the no-notch profile ‘‘nlvl’’. The discrepancy between predictions based on the ‘‘thk8’’ and the ‘‘nlvl’’ viscosity profiles is even less at ; 3%. The sensitivity of J˙2 predictions to the existence of a low-viscosity notch at the base of the upper mantle has important implications for a variety of geophysical applications based upon this observable. One such application is the estimate of bounds on the present-day global cryospheric mass flux. This mass flux can be separated into two sources: small mountain glaciers and ice sheets, and the large polar ice caps ŽGreenland, Antarctica.. The retreat of 31 small mountain glaciers and ice sheets has been tabulated by w42x and the J˙2 signal associated with this forcing has been predicted to be ; 10y1 1 yry1 Že.g., w27,43,44x.. The difference between the J˙2 predictions for models ‘‘nlvl’’ and ‘‘thn’’ is 2 = 10y1 1 yry1 ŽFig. 4A., or twice the signal from Meier’s w42x sources. The J˙2 signal associated with mass flux in the Antarctic and Greenland ice sheets is ; 4 = 10y1 1 yry1 per mmryr of net eustatic sea-level rise w27x. ŽVariations in the mass of the Antarctic and Greenland ice sheets are equally efficient at exciting a J˙2 response, and one need not consider their independent contributions.. Therefore, the thin low-viscosity notch which characterizes model ‘‘thn’’ has an effect on J˙2 which is comparable to the signal associated with a polar melting event which yields a eustatic sea-level rise of 0.5 mmryr. We conclude that failure to account for the possible presence of such a low-viscosity layer could lead to a ; 0.5 mmryr overestimate of the ongoing melt signal associated with the large polar ice sheets. The observed J˙2 likely has contributions from a number of geophysical phenomena unrelated to cryospheric forcings. For example, pressure changes in the fluid outer core of the Earth may contribute as much as ; y1.3 = 10y1 1 yry1 w34x. ŽOther effects associated with the core–mantle boundary have also been described w45x which produce comparable sig-

273

nals.. This contribution is ; 60% of the signal associated with the presence of a thin low-viscosity zone immediately above the 670 km boundary. A second common application of the J˙2 observable involves the inference of lower-mantle viscosity Že.g., w15,46x.. This application has generally involved a standard three-layer parameterization of the radial viscosity structure which includes an elastic lithosphere and isoviscous upper- and lower-mantle regions. In Fig. 5 Ždotted line. we plot predictions of J˙2 for the case of an 80 km lithosphere, an uppermantle viscosity of 10 21 Pa s, and a lower-mantle viscosity which ranges from 10 21 to 10 23 Pa s Žas shown on the abscissa.. The ‘‘classic’’ inverted parabolic form of the predicted curve is now well understood; small amplitudes are associated both with low-viscosity models, which have relaxed close to equilibrium in the 5 kyr since the end of deglaciation, and high-viscosity models, which are character-

Fig. 5. Predictions of the GIA-induced J˙2 signal for a suite of viscosity models. In all cases a lithospheric thickness of 80 km is used. Dotted line — the upper-mantle viscosity is fixed to 10 21 Pa s and the lower-mantle viscosity is varied according to the abscissa. Solid line — same as the dotted line, with the exception that a 70 km thick zone of low- Ž10 19 Pa s. viscosity is introduced immediately above the 670 km discontinuity. Long-dashed-dotted line — same as the dotted line, with the exception that a 70 km thick zone of high- Ž10 23 Pa s. viscosity is introduced immediately below the 670 km discontinuity. Short-dashed-dotted line — same as the dotted line with the exception that both the high- and low-viscosity 70 km zones are introduced to generate a ‘‘dipole’’ viscosity structure at 670 km depth Žsee text..

274

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

ized by slow rates of adjustment at all times. Note that the Žno notch. model ‘‘nlvl’’ ŽFig. 1. has an average lower-mantle viscosity of ; 5 = 10 21 Pa s. The J˙2 prediction for this profile, ; y5.3 = 10y1 1 yry1 ŽFig. 4A., is essentially the same as the prediction for the simple three-layer model in which the isoviscous lower-mantle viscosity is set to 5 = 10 21 Pa s Žsee Fig. 5, dotted line.. This suggests that when a pronounced low-viscosity layer is not present, the J˙2 predictions are largely determined by the average viscosity within the lower mantle, and so a simple three-layer parameterization is adequate. The solid line in Fig. 5 is analogous to the dotted line with the exception that a thin Ž70 km. lowviscosity notch has been added to the base of the upper mantle. To be consistent with the profile ‘‘thn’’, the notch viscosity is taken to be two orders of magnitude lower than the mean upper-mantle value Žin this case, the notch viscosity is 10 19 Pa s.. In the case of a lower-mantle viscosity of 5 = 10 21 Pa s, the introduction of the low-viscosity notch increases the value of the J˙2 prediction from y5.3 = 10y1 1 to y3.7 = 10y1 1 yry1 , which is consistent with the results in Fig. 4A. ŽThat is, it is consistent with the difference in J˙2 predictions based on models ‘‘thn’’ and ‘‘nlvl’’.. A comparison of the dotted and solid lines in Fig. 5 leads to two important conclusions. First, the existence of a thin low-viscosity layer may significantly alter inferences of lowermantle viscosity. In particular, the J˙2 signal associated with a mean lower-mantle viscosity of 5 = 10 21 Pa s and the presence of a thin low-viscosity notch Žsolid line. is the same as the signal derived for a mean lower-mantle viscosity of ; 2 = 10 21 Pa s in the absence of the notch Žsee Fig. 5.. We conclude that previous three-layer modeling of the type associated with the dotted line in Fig. 5 Žlithosphere, upper mantle, lower mantle. will have significantly underestimated the mean lower-mantle viscosity if finescale viscosity structure exists at the base of the upper mantle. Second, examination of Fig. 5 indicates that our conclusion that J˙2 predictions are significantly influenced by the presence of a lowviscosity notch at the base of the upper mantle is valid for a wide range of lower-mantle viscosities. Predictions of the J˙l Ž2 F l F 10. signal for a suite of three-layer viscosity profiles similar to those used to calculate the dotted line in Fig. 5 appear in

w27x; Fig. 3. The study w27x found that a reduction in the average viscosity of the upper mantle from 10 21 to 10 20 Pa s reduced the maximum amplitude of the J˙2 prediction by ; 0.8 = 10y1 1 yry1. This is approximately half the effect of introducing a thin Ž70 km. layer of viscosity 10 19 Pa s within an upper mantle of viscosity 10 21 Pa s. As we have discussed, previous studies w7,18x have considered the influence on J˙2 predictions of a thin highly viscous layer located at 670 km depth. The dashed-dotted lines in Fig. 5 revisit this issue. The long-dashed-dotted line is analogous to the solid with the exception that a thin Ž70 km. high-viscosity Ž10 23 Pa s. layer is introduced just below 670 km depth Žin this case the abscissa refers to the viscosity of the remaining portion of the lower mantle.. The results are consistent with those appearing in w7,18x and they suggest that the influence of a high-viscosity layer is significantly less than the effect of a low-viscosity layer at the upper mantle–lower mantle boundary. Indeed, in the case of a mean lowermantle viscosity of 5 = 10 21 Pa s, the influence of the latter is roughly four times the former; this factor rises to an order of magnitude for a lower-mantle viscosity of 10 22 Pa s. The influence of a ‘‘dipole’’ viscosity structure at 670 km depth on convection-related observables has recently been considered w4x. This structure is characterized by a thin low-viscosity layer above the 670 km discontinuity and a high-viscosity zone of equal thickness below the discontinuity. The short-dasheddotted line in Fig. 5 represents such a model. It is the same as the model used to generate the solid line with the exception that a 70 km thick layer of viscosity 10 23 Pa s is added below the 670 km boundary. Comparison of the long-dashed-dotted and short-dashed-dotted lines on the figure indicates that the high-viscosity layer within the dipole structure has a relatively minor effect on the J˙2 predictions. Our demonstration that a low-viscosity zone at the base of the upper mantle has a significant effect on ŽGIA. predictions of J˙2 suggests that it would be of interest to consider the effect on higher-degree zonal harmonics. In Fig. 6 we examine this issue for the case of J˙l Ž3 F l F 6.. The dotted line in each frame corresponds to the three-layer viscosity parameterization adopted to calculate the dotted curve in Fig. 5 Ž80 km lithosphere, 10 21 Pa s upper-mantle viscosity

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

275

Fig. 6. Predictions of the GIA-induced J˙l Ž3 F l F 6. coefficients for a suite of viscosity models. All predictions adopt an 80 km thick lithosphere. The dotted line is computed using an upper-mantle viscosity fixed to 10 21 Pa s and a lower-mantle viscosity varied according to the abscissa. The solid line uses the same viscosity structure with the exception that a 70 km thick zone of low viscosity is introduced immediately above the 670 km boundary.

and a lower-mantle viscosity given by the abscissa.. The solid line uses the same viscosity structure with the exception that a 70 km thick low-viscosity zone is introduced immediately above 670 km Žthis suite of models is identical to the one adopted to generate the solid line in Fig. 5.. Clearly, the inclusion of a low-viscosity zone has a significant effect on the secular variations of the higher-degree zonal harmon-

ics. As an example, the maximum amplitude of the J˙4 prediction drops by ; 30% when the thin lowviscosity structure is introduced. We also note that the low-viscosity zone acts to reduce the slope of the high- Žlower mantle. viscosity limb on each of the predicted curves Žas in Fig. 5 for the case of J˙2 .. We have also calculated the J˙l Ž3 F l F 6. signal for a sequence of models in which a thin high-viscosity

276

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

layer was introduced immediately below 670 km depth Žas in the dashed-dotted lines in Fig. 5.. Our results Žnot shown in Fig. 6. indicate that the inclusion of the high-viscosity feature has a maximum effect of ; 10% for models with a moderate increase in viscosity Ž2–3 times. across the 670 km boundary. There is only a few percent effect, at most, for models with higher lower-mantle viscosities. 3. Concluding remarks We have shown that a thin low-viscosity notch at the base of the upper mantle, recently inferred on the basis of convection-related long-wavelength nonhydrostatic geoid and free-air gravity harmonics, can have a significant influence on a suite of observables associated with GIA. Perhaps the most important effect we have isolated involves the sensitivity of predicted J˙2 Žor rotation rate. variations to the viscosity of the notch region. This sensitivity can bias constraints on either the mean lower-mantle viscosity or the present-day mass balance of large polar ice sheets which have been derived from simple models that do not include the presence of fine-scale lowviscosity structure at 670 km depth. These results counter recent suggestions that GIA observables are insensitive to such structure Že.g., w19,20x.. Our predictions of RSL variations for previously glaciated regions indicate that decay time estimates can be significantly reduced by the presence of a low-viscosity zone, with sites nearer the center of the rebound area showing greater sensitivity. Also, the presence of a thin low-viscosity layer above the 670 km discontinuity reduces the predicted amplitude of Late Holocene highstands in the far field of the Late Pleistocene ice sheets. While this sensitivity will have bearing on estimates of Late Holocene melting events, inferences of lower-mantle viscosity based on differential highstands from these far-field regions w16x will not be altered by the presence or absence of the low-viscosity notch. RSL decay time data and long-wavelength free-air gravity harmonics have been simultaneously inverted w21,22x for radial mantle viscosity structure. The solution profile obtained in these studies is characterized by a low-viscosity region immediately above the 670 km boundary which extends through the

transition zone. The convection-related harmonics played an important role in constraining this lowviscosity region ŽFig. 2.. The present results indicate that the site-dependent RSL decay times also provided constraints on this feature of the inversion. J˙2 , polar wander speed and RSL decay time data associated with GIA have been used by Peltier w20x to invert for the radial profile of mantle viscosity. The inverted profile was then compared to an independent viscosity model w4x inferred on the basis of low-degree nonhydrostatic geoid data alone. The latter model was scaled to provide a qualitative fit between the two predictions and the general agreement led Peltier w20x to conclude that the rebound-inferred viscosity profile simultaneously reconciled both GIA and convection-related observables. The inference w4x was characterized by a thin low-viscosity zone at the base of the upper mantle. The absence of such a feature in the GIA inversion w20x was not considered to be of significance. Indeed, Peltier w20x assumed that the fine-scale low-viscosity feature Žw20x, p.1364.: does not represent a significant departure from the GIA profile, because the GIA data are not particularly sensitive to its presence. We conclude that this assumption is invalid. We have demonstrated that a subset of RSL decay times, and in particular the J˙2 datum, are significantly sensitive to the presence of the low-viscosity region. The determination of a viscosity model which reconciles both GIA and convection observables requires a simultaneous consideration of both data sets w21,22x. References w1x S.D. King, G. Masters, An inversion for radial viscosity structure using seismic tomography, Geophys. Res. Lett. 19 Ž1992. 1551–1554. w2x A.M. Forte, A.M. Dziewonski, R.L. Woodward, Aspherical structure of the mantle, tectonic plate motions, nonhydrostatic geoid, and topography of the core–mantle boundary, in: J.-L. LeMouel, ¨ D.E. Smylie, T. Herring ŽEds.., Dynamics of the Earth’s Deep Interior and Earth Rotation, Am. Geophys. Union, Washington, DC, Geophys. Monogr. Ser. 72 Ž1993. 135–166. w3x S.V. Panasyuk, B.H. Hager, A.M. Forte, Understanding the effects of mantle compressibility on geoid kernels, Geophys. J. Int. 124 Ž1996. 121–133. w4x G. Pari, W.R. Peltier, The heat flow constraint on mantle

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

w5x

w6x

w7x

w8x

w9x w10x

w11x

w12x w13x

w14x

w15x

w16x

w17x w18x

w19x

w20x w21x

w22x

tomography-based convection models: Towards a geodynamically self-consistent inference of mantle viscosity, J. Geophys. Res. 100 Ž1995. 12731–12751. C.G. Sammis, J.L. Dein, On the possibility of transformational superplasticity in the Earth’s mantle, J. Geophys. Res. 79 Ž1974. 2961–2965. S.V. Panasyuk, B.H. Hager, A model of transformational superplasticity of the upper mantle, Geophys. J. Int. Ž1996, submitted.. W.R. Peltier, The LAGEOS constraint on deep mantle viscosity: Results from a new normal mode method for the inversion of viscoelastic relaxation spectra, J. Geophys. Res. 90 Ž1985. 9411–9421. V. Corrieu, Y. Ricard, C. Froidevaux, Converting mantle tomography into mass anomalies to predict the Earth’s radial viscosity, Phys. Earth. Planet. Inter. 84 Ž1994. 3–13. N.A. Haskell, The motion of a fluid under a surface load, 1, Physics 6 Ž1935. 265–269. R.K. McConnell, Viscosity of the mantle from relaxation time spectra of isostatic adjustment, J. Geophys. Res. 73 Ž1968. 7089–7105. R.J. O’Connell, Pleistocene glaciation and the viscosity of the lower mantle, Geophys. J.R. Astron. Soc. 23 Ž1971. 299–327. L.M. Cathles, The Viscosity of the Earth’s Mantle, Princeton University Press, Princeton, NJ, 1975. W.R. Peltier, J.T. Andrews, Glacial isostatic adjustment, I. The forward problem, Geophys. J.R. Astron. Soc. 46 Ž1976. 605–646. R. Sabadini, D.A. Yuen, E.V. Boschi, A comparison of the complete and truncated versions of the polar wander equations, J. Geophys. Res. 87 Ž1982. 2885–2903. P. Wu, W.R. Peltier, Pleistocene deglaciation and the Earth’s rotation: A new analysis, Geophys. J.R. Astron. Soc. 76 Ž1984. 753–792. M. Nakada, K. Lambeck, Late Pleistocene and Holocene sea-level change in the Australian region and mantle rheology, Geophys. J. Int. 96 Ž1989. 497–517. J.X. Mitrovica, Haskell w1935x revisited, J. Geophys. Res. 101 Ž1996. 555–569. E.R. Ivins, C.G. Sammis, C.F. Yoder, Deep mantle viscous structure with prior estimate and satellite constraint, J. Geophys. Res. 98 Ž1993. 4579–4609. W.R. Peltier, A.M. Forte, J.X. Mitrovica, Mantle viscosity and global dynamics: The compatibility of inferences from glacial isostasy and convection related constraints, Eos ŽTrans. Am. Geophys. Union. 73 Ž43. Ž1992. 62. W.R. Peltier, Mantle viscosity and ice-age ice sheet topography, Science 273 Ž1996. 1359–1364. A.M. Forte, J.X. Mitrovica, New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and post-glacial rebound data, Geophys. Res. Lett. 23 Ž1996. 1147–1150. J.X. Mitrovica, A.M. Forte, The radial profile of mantle viscosity: Results from the joint inversion of convection and post-glacial rebound observables, J. Geophys. Res. 102 Ž1997. 2751–2769.

277

w23x P. Wu, W.R. Peltier, Glacial isostatic adjustment and the free air gravity anomaly as a constraint upon deep mantle viscosity, Geophys. J.R. Astron. Soc. 74 Ž1983. 377–449. w24x K. Lambeck, Glacial rebound of the British Isles, II. A high-resolution, high-precision model, Geophys. J. Int. 115 Ž1993. 960–990. w25x D. Breur, D. Wolf, Deglacial land emergence and lateral upper-mantle heterogeneity in the Svalbard Archipelago, I. First results for simple load models, Geophys. J. Int. 121 Ž1995. 775–788. w26x J.X. Mitrovica, W.R. Peltier, A new formalism for inferring mantle viscosity based on estimates of post glacial decay times: Application to RSL variations in N.E. Hudson Bay, Geophys. Res. Lett. 20 Ž1993. 2183–2186. w27x J.X. Mitrovica, W.R. Peltier, Present-day secular variations in the zonal harmonics of the Earth’s geopotential, J. Geophys. Res. 98 Ž1993. 4509–4526. w28x R. Sabadini, W.R. Peltier, Pleistocene deglaciation and the Earth’s rotation: Implications for mantle viscosity, Geophys. J.R. Astron. Soc. 66 Ž1981. 553–578. w29x D.A. Yuen, R. Sabadini, E.V. Boschi, Viscosity of the lower mantle as inferred from rotational data, J. Geophys. Res. 87 Ž1982. 10745–10762. w30x C.F. Yoder, E.R. Ivins, Changes in the Earth’s gravity field from Pleistocene deglaciation and present-day glacial melting, Eos ŽTrans. Am. Geophys. Union. 66 Ž18. Ž1985. 245. w31x A.S. Trupin, M.F. Meier, J.M. Wahr, Effect of melting glaciers on the Earth’s rotation and gravitational field: 1965– 1984, Geophys. J. Int. 108 Ž1992. 1–15. w32x L.L.A. Vermeersen, R. Sabadini, G. Spada, N.J. Vlaar, Mountain building and Earth rotation, Geophys. J. Int. 117 Ž1994. 610–624. w33x B.M. Steinberger, R.J. O’Connell, Change of the Earth’s rotation axis inferred from advection of mantle density heterogeneities, Nature ŽLondon. 387 Ž1997. 169–173. w34x M. Fang, B.H. Hager, T.A. Herring, Surface deformation caused by pressure changes in the fluid core, Geophys. Res. Lett. 23 Ž1996. 1493–1493. w35x S.D. King, Radial models of mantle viscosity: Results from a genetic algorithm, Geophys. J. Int. 122 Ž1995. 725–734. w36x A.M. Dziewonski, D.L. Anderson, Preliminary Reference Earth Model ŽPREM., Phys. Earth Planet. Inter. 25 Ž1981. 297–356. w37x A.M. Tushingham, W.R. Peltier, ICE-3G: A new global model of late Pleistocene deglaciation based on geophysical predictions of post-glacial relative sea level change, J. Geophys. Res. 96 Ž1991. 4497–4523. w38x C. Hillaire-Marcel, Multiple component postglacial emerŽEd.., gence, eastern Hudson Bay, Canada, in: N. Morner ¨ Earth Rheology, Isostasy, and Eustasy, Wiley, New York, NY, 1980, pp. 215–230. w 39 x W .R. Peltier, Glacial isostatic adjustment, II. The inverse problem, Geophys. J.R. Astron. Soc. 46 Ž1976. 669–706. w40x J.X. Mitrovica, W.R. Peltier, On postglacial geoid subsidence over the equatorial oceans, J. Geophys. Res. 96 Ž1991. 20053–20071.

278

G.A. Milne et al.r Earth and Planetary Science Letters 154 (1998) 265–278

w41x J.A. Clark, W.E. Farrell, W.R. Peltier, Global changes in postglacial sea level: A numerical calculation, Quat. Res. 9 Ž1978. 265–287. w42x M.F. Meier, Contribution of small glaciers to global sea level, Science 226 Ž1984. 1418–1421. w43x R. Sabadini, D.A. Yuen, P. Gasperini, Mantle rheology and satellite signatures from present-day glacial forcings, J. Geophys. Res. 93 Ž1988. 437–447.

w44x W.R. Peltier, Global sea level and earth rotation, Science 240 Ž1988. 895–901. w45x M. Lefftz, H. Legros, Influence of viscoelastic coupling on the axial rotation of the Earth and its fluid core, Geophys. J. Int. 108 Ž1992. 725–739. w46x W.R. Peltier, X. Jiang, The procession constant of the Earth: Variations throughout the ice age, Geophys. Res. Lett. 21 Ž1994. 2299–2302.